East African Journal of

Biophysical and Computational

Sciences

( E A J B C S )

ISSN (Online): 2789-3618 and ISSN (Print): 2789-360X

COLLEGE OF NATURAL AND COMPUTATIONAL SCIENCES,

HAWASSA UNIVERSITY

June, 2026

Volume 7 Issue 1

East African Journal of Biophysical and

Computational Sciences

ISSN (Online): 2789-3618 and ISSN (Print): 2789-360X

East African Journal of

Biophysical and Computational

Sciences (EAJBCS)

ISSN (Online): 2789-3618 and ISSN (Print): 2789-360X

Volume 7 Issue 1

College of Natural and Computational Sciences

Hawassa University

June, 2026

Table of Contents

East African Journal of Biophysical

and Computational Sciences

(EAJBCS)

East African Journal of Biophysical and Computational Sciences

(EAJBCS) whose ISSN (Online): 2789-3618 and ISSN (Print): 2789-

360X is a double-blind peer-reviewed open-access journal published by

Hawassa University, College of Natural & Computational Sciences. is

Journal is a multi and interdisciplinary journal that is devoted to

attracting high-quality, latest, and valuable advancements in the ﬁelds of

natural sciences. e Journal invites publications from diﬀerent

geographical contexts and disciplines to advance the depths of

knowledge related to physics, chemistry, geology, biology, & veterinary

medicine. e manuscript originated from other sciences such as

biotechnology, sport science, statistics, and mathematics can also be

accepted based on their adjunct nature. e Journal encourages

publications of both scholarly and industrial papers on various themes

with the aim of giving innovative solutions to natural sciences. It

encourages the publishing of open access academic journals on a

regular basis (presumably biannual). e Journal publishes original

research articles, critical reviews, mini-reviews, short communications,

case reports related to the speciﬁc theme & a variety of special issues in

English. e Journal, published under the Creative Commons open

access license (CC BY-NC-ND), doesn't charge fees for publishing an

article and hence oﬀers an opportunity to all social classes regardless of

their economic statuses. is helps to promote academic research

published by resource-poor researchers as a mechanism to give back to

society.

e journal is already indexed on known databases like AJOL, DOAJ

and CABI ABSTRACTS

Editorial Team

Andualem G/Silase, MSc.

Editor-in-Chief

Assistant

Professor

of

Teaching

Physical

Education, HU

Admasu Tadesse, PhD

Assistant Professor of Mathematics , HU

admasut@hu.edu.et, +251-91-3267054

andualem.g@hu.edu.et, +251-911-082383

Wogene Talelign, MSc.

Editorial Manager

Lecturer of Engineering Geology, HU

wegenetalelign@hu.edu.et, +251-913-939058

Abnet Woldesenbet, PhD

Assistant Professor of Aquatic Ecology, HU

abnetm@hu.edu.et, +251-911811819

Dereje Danbe, PhD

Assistant Professor of Applied Statistics, HU

derejedanbe@hu.edu.et, +251-913-927596

Associate Editors

Abebe Getachew, PhD

Assistant Professor of Solid State Physics, HU

abebeg@hu.edu.et, +251-911-362198

Yonnas Shuke, PhD

Assistant Professor of Applied Statistics, HU

yonasshuke@hu.edu.et, +251-910-191357

Prof. Desie Sheferaw

Professor of Veterinary Epidemiology, HU

desies@hu.edu.et, +251-916-832419

Firew Kebede, PhD

Associate Professor of Botanical Sciences, HU

ﬁrew@hu.edu.et, +251-911-342084

Prof. Sisay Tadesse, PhD

Professor of Physical Chemistry, HU

sisaytad@hu.edu.et, +251-922-598889

Abrham Mikru, PhD

Assistant Professor of Applied Microbiology,

HU

abrahammikru@hu.edu.et, +251-916-867353

Tegene Tesfaye, PhD

Associate Professor of Organic Chemistry, HU

tegenetesfaye@hu.edu.et, +251-942-495546

Kiros Gebreargawi, PhD

Girma Tilahun, PhD

Assistant Professor of Applied Mathematics,

HU

kirosg@hu.edu.et, +251-926-528484

Associate Professor of Limnology, HU

girma@hu.edu.et, +251-932-206985

Sintayehu Tesfa, PhD

Advisory Board

Professor Zinabu G/mariam

Professor Freshwater Ecology, HU, Ethiopia

luzinabu@gmail.com

Assistant Professor of Quantum Optics, Jazan

University, Saudi Arabia

sint_tesfa@yahoo.com,

Prof. Natarajan Pavanasam

Zeytu Gashaw, PhD

Professor of aquatic science and aquaculture,

IVRI, India

drpnatarajan123@gmail.com,

Associate Professor of Applied Statistics, AAU,

Ethiopia

zeytu.gashaw@aau.edu.et

Prof. Legesse Kassa

Dr. Yifat Denbarga

Professor of Statistics, University of South

Africa, SA

debuslk@unisa.ac.za,

Associate Professor of Tropical Veterinary

Medicine, HU, Ethiopia

yifatd@hu.edu.et,

Prof. Endrias Zewdu

Prof. Abebe Geletu

Professor of Veterinary Public Health, Ambo U,

Ethiopia

endrias.zewdu@gmail.com,

Professor of Process Optimization, Technical

University of Ilmenau, Germany

abebe.geletu@tu-ilmenau.de,

Edessa Negera Gobena, PhD

Prof. Bekele Megersa

Senior Lecturer at School of Health, Sport and

Bioscience, University of East London, UK

E.N.Gobena@uel.ac.uk,

Professor of Veterinary Epidemiology, AAU,

Ethiopia

bekelebati@gmail.com,

Layout Designer

Geda Hoka Homecho

Prof. Abiy Yenesew

Professor of Natural Product Chemistry, Nairobi

University, Kenya

ayenesew@uonbi.ac.ke,

Lecturer,

Department of Geography and

Environmental Studies

PhD candidate, Wondogenet College of Forestry

and Natural Resources

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

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

ARTICLE

Computational Study of MHD Blood Flow

through Bifurcated Artery Using

Caputo-Fabrizio Fractional Derivative,

Thermal Radiation, and Magnetic Field for

Tumor Therapies

ARTICLE INFO

Volume 7(1), 2026

https:

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

Isah Abdullahi¹, Dauda Gulibur Yakubu^1,∗, Muhammad

ARTICLE HISTORY

Shamsuddeen Dauda², Mahmood Abdulhameed³,Saidu Abubakar

Kadas⁴,Mohammed Abdulhameed⁵, and Garba Tahiru Adamu⁶

Received: March 10, 2026

Accepted: 23 May, 2026

Published Online: 10 June, 2026

CITATION

¹Department of Mathematical Sciences, Abubakar Tafawa Balewa University, Bauchi, Nigeria

²Department of Biological Sciences, Abubakar Tafawa Balewa University, Bauchi, Nigeria

³Department of Electrical Electronic Engineering, Abubakar Tafawa Balewa University, Bauchi, Nigeria

⁴Department of Obstetrics Gynaecology, ATBU, Teaching Hospital, Bauchi, Nigeria

⁵School of Science and Technology, The Federal Polytechnic Bauchi, Nigeria

Abdullahi et.al (2026). Computational

Study of MHD Blood Flow through

Bifurcated Artery Using Caputo-Fabrizio

Fractional Derivative, Thermal

Radiation, and Magnetic Field for Tumor

Therapies. East African Journal of

Biophysical and Computational

⁶Department of Mathematical Sciences, Bauchi State University, Gadau, Bauchi, Nigeria

^∗Corresponding author: dgyakubu@atbu.edu.ng

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

doi.org/10.4314/eajbcs.v7i1.2S.18-26

Abstract

OPEN ACCESS

This study investigates the impact of heat sources, thermal radiation, and chemical reactions on the

magnetohydrodynamic blood flow through a bifurcated artery in the presence of a slanted magnetic

field. Using Laplace transform and the method of undetermined coefficients, the constitutive equations

for the mathematical model of Caputo-Fabrizio fractional derivative order were solved. Blood flow

velocity, temperature distribution, and concentration were found to have analytical expressions. The

effects of certain physical parameters on blood velocity, temperature and concentration are graphically

represented, and these representations accurately depict the flow disturbances. We discovered that

the bifurcation apex of the artery with a symmetrical divider has steady blood flow. This may lead

to significant shear stresses on either side of the bifurcation. Near the apex, when the flow is

substantially different, obstruction may result from the formation of boundary layers on the inner walls

of the bifurcation. Sluggish flow also occurred along the outer walls of the bifurcation. It has also

been discovered that the temperature distribution, concentration, and arterial blood flow velocity are

significantly influenced by the fractional order parameter, the slanted magnetic fields, the heat source,

and the chemical reaction parameter. This study offers significant benefits for medical applications in

biomechanical engineering, biomedical engineering, and medicine.

This work is licensed under the Creative

Commons open access license (CC

BY-NC 4.0).

East African Journal of Biophysical and

Computational Sciences (EAJBCS) is

already indexed on known databases

like AJOL, DOAJ, CABI ABSTRACTS and

FAO AGRIS.

Keywords: Chemical reaction; Heat source; MHD Blood low; Slanted magnetic ield; Thermal

radiation

(Shit & Majee, 2015). Understanding many facets (aspects) of the

medical sciences, such as homeostasis, treating cancerous tumors, and

administering medication using magnetic particles, depends on the

study of biomagnetic ﬂuid dynamics (Shaw & Murthy, 2010). Blood’s

hemoglobin molecules are regarded as biomagnetic ﬂuids with magnetic

properties. Blood can also be considered as a Newtonian ﬂuid if it ﬂows

through the bigger arteries at a high shear rate. These arteries are thought

to be homogenous whose ﬂow behavior can be described by a Newtonian

model (see, Caro et al., 2011; MacDonald, 1979). Numerous researchers

have looked into various possibilities for studying physiological ﬂuids

1 Introduction

Bio-magnetic ﬂuid dynamics (BFD), which is the study of bio-ﬂuid ﬂow in

the presence of magnetic ﬁeld, is a rapidly developing subject of study in

ﬂuid mechanics (Tzirtzilakis, 2005). This ﬁeld of study is of tremendous

importance to the ﬁeld of medical science and has the possibility

to be utilized in a diversity of domains, including the delivery of

medications via the utilization of magnetized particles, the management

of severe bleeding, and the assistance in dealing with malignant cancers

Abdullahi et.al (2026)

18

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

using porous media (see, Bhatti & Lu, 2019; Dash et al., 1996; Ramesh &

Devakar, 2015; Shit & Roy, 2015) developed blood ﬂow model via porous

medium. Based on Darcy’s law, Bhargava et al. (2007) and Ghasemi et al.

(2015) investigated the pulsatile ﬂow and mass transfer of an electrically

conducting Newtonian bio-ﬂuid via a channel comprising porous media

using blood as the porous medium ﬂuid. Bhatti et al. (2018) developed a

mathematical model to investigate heat transfer, mass transfer, and blood

ﬂow in a porous medium channel while accounting for the integrated

Darcy-Brinkman-Forchheimer model. Blood behaves non-Newtonian

even in larger arteries at low shear rates, as demonstrated by Liepsch

(1986). When blood ﬂows through arteries at a low shear rate, it can be

treated as Cassons ﬂuid (Srivastava & Srivastava, 1984). Many researches

have supported the Casson ﬂuid model for blood ﬂow via tiny arteries at

low shear rates (see,Hayat et al., 2016; Nagarani et al., 2006; Venkatesan

et al., 2013 ). Many authors (see, Abdulhameed et al., 2017; Misra & Shit,

2009; Mondal & Shit, 2017; Yakubu et al., 2020; Zeeshan et al., 2017 )

have regarded blood as a non-Newtonian ﬂuid, because of its electrical

conductivity, displays magneto hydrodynamic behavior.

Therefore, it wasn’t until the last few decades that a signiﬁcant number

of scholars started to highlight the fact that diﬀerential equations and

fractional derivatives have numerous applications in a variety of domains

(see, Abdulhameed et al., 2023; Imoro et al., 2024). These days, fractional

derivative order diﬀerential equation problems are the most eﬀective

and successful ways to model the nonlinear processes that emerge in

many domains of applied study, including biology, chemistry, ecology,

engineering, and many other application areas. Several mathematical

models have shown that they oﬀer a more realistic depiction of the

phenomenon under research. Examples of these models include those

employed in biomedical engineering, viscoelastic mechanics, boundary

layers, electromagnetic, and porous media. Bansi et al. (2018) investigated

a fractional blood ﬂow model in the oscillatory artery with magnetic ﬁeld

and heat radiation eﬀects. With the aid of fractional time derivative,

(Yakubu et al., 2021) examined the eﬀects of pressure gradient, body

acceleration, and magnetic ﬁeld on blood ﬂow through artery. The

eﬀects of blood ﬂow parameters, Caputo’s time fractional derivatives, and

the external magnetic ﬁeld on the cylindrical domain were studied by

(Shah et al., 2016). Ali et al. (2017) solved a fractional order model for

Cassons ﬂuid ﬂow using the Hankel transform and Laplace transform

techniques to determine the exact solutions. He and collaborators (2019)

used the fractional order Caputo derivative to investigate the complexity

of blood in arteries under various forces. In the ﬁeld of medicine, magneto

hydrodynamic ﬂow plays a crucial role. It is considered for the reduction

of bleeding from wounds and for the treatment of malignant tumors.

Kumar et al. (2021) employed a chemical reaction, heat source, and

inclined magnetic ﬁeld to cure malignancies.

Many authors considered the examination of the heat and mass

transfer occurrences generated from these processes to be a highly

relevant element with respect to modeling physiological processes

(Prasad et al., 2025) and industrial processes (Sademaki et al., 2026).

Electromagnetohydrodynamics is the study of ﬂuids whose motion is

constantly aﬀected by externally applied magnetic ﬁeld and electric

ﬁeld. In order to comprehend the impact of magnetohydrodynamic

(MHD) and electrohydrodynamic (EHD) forces on the ﬂow of normal

ﬂuids, including blood, several studies have mostly concentrated on

the theoretical, computational, and experimental aspects of these forces.

Cell-based therapies, medication delivery, and biological processes are

just a few of the ﬁelds where the application of (EHD) has shown notable

advancement. Additional force components, primarily the Lorentz and

Coulomb forces, are incorporated into momentum equations and have a

direct eﬀect on ﬂuid velocity. Magnetohydrodynamics or MHD has been

used in a wide variety of biomedical applications (Vardanyan, 1973).

The fractional order time derivative of MHD blood ﬂow via a bifurcated

artery in the presence of a slanted magnetic ﬁeld, as well as the

coupling impact of heat transfer and blood ﬂow concentration, are

described here using Newtonian ﬂuid. The goal is to investigate

how magneto-hydrodynamic blood ﬂow through a bifurcated artery is

aﬀected by thermal radiation and a slanted magnetic ﬁeld during tumor

treatments. The Laplace transform and the indeterminate coeﬃcients

approach were used to ﬁnd the exact solutions, which were then

simulated to produce graphical outputs and the implications of several

important parameters on the outcomes were explored. The study was

motivated by the fact that there is currently very little information

available on the ﬂow in arterial bifurcation since the phenomena is

currently not stringent to mathematical analysis or precise experimental

measurement. The present investigation shows that the vast number of

variables involved are the main challenge in both situations.

The heat transmission and magnetohydrodynamic (MHD) blood ﬂow in

a restricted artery were studied by Majee and Shit (2017). Akbar and Butt

(2017) considered ferromagnetic blood ﬂow in a restricted, smaller artery

with a porous wall. The radiant heat transfer that takes place in the blood

vessels must also be considered while treating hyperthermia. Oncology

professionals are familiar with the medical practice of using heat therapy

to cancer patients. Chinyoka and Makinde (2014) investigated the eﬀects

of magnetic ﬁelds and heat radiation on arterial blood ﬂow. Sinha and

Shit (2015) investigated the magnetic hydrodynamic blood ﬂow in the

presence of thermal radiation. Tabi et al. (2017) studied the combined

eﬀects of magnetic ﬁelds and external radiation on blood ﬂow in the major

blood arteries. Yakubu et al. (2022) examined blood ﬂow of Oldroyd-B

ﬂuids in order to investigate the erratic ﬂow, with magnetic ﬁeld applied

perpendicular to the ﬂow direction. Heat transfer processes were studied

in the peristaltic ﬂow of blood with variable viscosity particle-liquid

suspensions by Bhatti et al. (2016). Blood ﬂow is greatly aﬀected when

the human body is exposed to a vibratory environment, as occurs when

operating machines or traveling in spacecraft. When the human body

undergoes body acceleration, a number of health problems might arise,

such as an elevated heart rate and vision loss. In the study of the impact of

body acceleration, a number of researches have produced mathematical

simulations of oscillatory blood ﬂow (see, Chaturani & Palanisamy, 1990;

Ghasemi et al., 2016; Sud & Sekhon, 1984 ). Bhatti and Lu (2019)

investigated the propagation of a hydro elastic single wave in a channel

with uniform ﬂow. Blood ﬂow characteristics have been discovered

to promote blood velocity in a vibratory environment using fractional

order derivative diﬀerential equation problems. Fractional diﬀerential

equations are the most used method for modeling natural phenomena.

This is due to the fact that equations oﬀer the possibility for a system to

either retain memory or to be hereditary with the properties of its history,

similar to how dynamic systems work (Syed et al., 2026).

2 Methods

2.1 Physical Structure and Mathematical

formulation

Blood considered in this study, is Newtonian, incompressible,

homogeneous, sticky ﬂuid that ﬂows from the trunk to the branches.

A mass stream’s rate at any cross-section that is perpendicular to its

direction is equal to m = 2bv, where b is the stream’s radius and v is

its mean speed. The mass stream’s speed at any cross section of the

extended channel is equal to m/2, and the bifurcating divider (internal

apical curve) has no eﬀect on this (see Figure 1). The magnetic ﬁeld is

applied to the ﬂow at an angle (φ) since the evaluated magnetic Reynolds

number is low. Therefore, it is believed that the magnetic and electric

ﬁelds produced by blood ﬂow are insigniﬁcant, the angle of bifurcation

is set to zero (Θ = 0), that is, the blood ﬂow region is bifurcated into two

streams that ﬂow parallel to the principal artery (the trunk). Figure 2

demonstrates the smooth muscle ﬁbers of the three concentric layers that

make up the walls of a typical elastic artery. These ﬁbers are controlled

by the sympathetic nervous system to contract or relax.

Fractional order derivatives have been applied in many ﬁelds of study,

including the complicated dynamics and rheological properties of

diﬀerent kinds of ﬂuids. The behavior of ﬂuid ﬂow is well depicted by

substituting fractional-order derivatives for the ordinary time derivative

in the constitutive equations (see, Atangana and Baleanu, 2016; Caputo

and Fabrizio, 2015; Samko et al., 1993). The concept of fractional calculus

was initially proposed by L’Hôpital in 1695, more than four centuries ago.

Abdullahi et.al (2026)

19

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

where D is the diﬀusion coeﬃcient and G = −k₁(C − C_∞) represents

chemical reaction rate in the ﬂuid ﬂow. It is important to mention that the

eﬀect of an electric ﬁeld in the concentration equation was also ignored.

θ = `^−z˙

,

u = `^−z˙

,

v = `^−z˙

,

C = `^−z˙at y = −1,

1

(5)

and θ → 0, u → 0, v → 0, C → 0 at y = 1.

However, by using the proper normalizing factors, the governing

equations (1)–(4) can be converted to dimensionless form. We present

the non-dimensional parameters as follows:

Figure 1: Physical ﬂow diagram of the bifurcated artery with zero angle of

bifurcation

x¯

b

y¯

b

u¯

v¯

dp¯/dx¯

u_HSηm/2b³ρ

x = , y = , u =

,

v =

, h(x, t) =

mu_HS/2bρ

ꢂ

ꢃ

ꢂ

ꢃ

3

2

3

2

¯

θ

C

2b

mηu

ρ

2b

mηu

ρ

¯

η

ρ

¯

t

k

ρ/η

ꢀ

ꢁ

k= ^ꢀ

and

^ꢁ, C =

, θ =

, τ =

,

t =

,

2

b

ρ/η

HS

−16δT₀³

¯

∂T

∂y¯

q¯ =

(6)

3k⁰

Applying (5) and (6) to eqns. (1)- (4) and removing the bars we obtain:

Figure 2: The structure of an artery walls (Transverse section through an artery)

∂u

∂t

∂²u

∂y²

ꢄ

u

k

2.2 Fundamental flow equations and their

solutions

+ h =

+ gβθ + gβ⁰C − M²sin²φ −

(7)

(8)

ꢅ

ꢄ

ꢅ

∂θ

=

∂t

1

∂²θ

∂y²

S

τp_r

+ R

+

θ

τp_r

Blood ﬂows through a porous media as two-limit layers when it is

subjected to magnetic ﬁelds and heat, with the assumptions made in the

numerical deﬁnition guiding its movement. In the stream ﬁeld headings

of x and y at time t , let u and v be the speed parts, η and ρ denote

blood density and thickness, respectively. Blood pressure is represented

by p , warm conductivity (K_T), C_pis the speciﬁc heat capacity at steady

strain, hotness is represented by Q , temperature is represented by T ,

the volumetric development boundary is represented by β , the angle of

the slanted (inclined) magnetic ﬁeld is represented by φ, and porosity

parameter is represented by K . With these, we have the equations

provided by (see, Ali et al., 2017; He et al., 2019; Kumar et al., 2021), with

some additional terms as follows:

∂u

∂x

∂v

∂y

+

= 0

(9)

∂C

∂t

1 ∂²C

S_C∂y²

=

− ωC

(10)

where

σB²

⁰, p_r=

16δT₀³

3k⁰τ

Qb²

k_T

τ

k₁b²

τ

ρC

k_T

M²

=

^ρ, R =

, S =

, S_C

=

, ω =

ρ

bD

2

∂u

∂t

1 ∂p

ρ ∂x

η ∂²u

cB_αsin φ

v

k

+

=

+ gβ(T − T_∞) + gβ⁰(C − C_∞) −

v −

2

ρ

∂y

are the magnetic ﬁeld parameter, Prandtl number, thermal radiation

parameter, heat source parameter, Schmidt number, chemical reaction

parameter and θ is the temperature conveyance. Now, using the

Caputo-Fabrizio fractional derivative as stated in (Caputo and Fabrizio

2015), we consider the time fractional momentum equations as:

(1)

∂²T

Q

ρC_o

∂q

∂y

¯

k_T

∂T

∂t

=

+

(T − T_o) −

(2)

(3)

∂y²

ρC_o

ꢄ

ꢅ

1

∂u(y, τ)

∂τ

α(t − τ)

t

^CFD_t^αu(y, t) =

exp

−

dτ, 0 < α < 1

0

1 − α

∂u

∂x

∂v

∂y

= 0

n

o

su(y, s) − u(y, 0)

L

^CFD_t^αu(y, t)

=

(1 − α)s + α

∂

∂t

(11)

where,

is the material time derivative. On the other hand, the

dimensionless concentration equation for medication (concentration)

delivery in magneto hydrodynamic blood ﬂow through permeable

bifurcated artery is provided by,

u(y, 0) =

(12)

∂²C

∂y²

The Caputo-Fabrizio derivative corresponding to equations (7), (8) and

(10) are as follows:

∂C

∂ f

= D

+ G

(4)

Abdullahi et.al (2026)

20

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

(

)

A₇+ A cosh A₂y + B sinh A₂y + A₈cos A₆y + A₉sin A₆y +

F =

.

∂²u

∂y²

ꢄ

u

k

+ A₁₀cosh Λy + A₁₁sinh Λy

^CFD_t^αu(y, t) + h =

+ gβθ + gβ⁰C − M²sin²φ −

(13)

(14)

(15)

(30)

ꢅ

1

τp_r

∂²θ

∂y²

S

τp_r

^CFD_t^αθ(y, t) =

+ R

+ (

)θ

We now have blood velocity in the axial direction, using equation (30) and

equation (19) as

∂C

∂t

1 ∂²C

S_C∂y²

^CFD_t^αC (y, t) =

=

− ωC.

!

A₇+ A cosh A₂y¯ + B sinh A₂y¯ + A₈cos A₆y¯ +

1

u¯(y, s) =

s + λ²

A₉sin A₆y¯ + A₁₀cosh Λy¯ + A₁₁sinh Λy¯

Applying Laplace transform to equations (13)-to-(15), and using the

boundary condition in equation (12) we have;

(31)

Equation (26) and equation (21) together provide the usual direction of

blood velocity in the bifurcated artery, which is given by

su(y, s)

(1 − α)s + α

∂²u¯

∂y¯²

u¯

k

0

2

¯

+ h =

+ gβθ + gβ C − M sin φ −

(16)

(17)

(18)

1

v¯(y, s) = A₁

(32)

s + λ²

ꢄ

ꢅ

2

¯

su(y, s)

(1 − α)s + α

1

τp_r

∂ θ

S

τp_r

¯

)θ

=

+ R

+ (

∂y¯²

According to equations (28) and (20), the temperature distribution in the

bifurcated artery is given as follows:

2

¯

su(y, s)

(1 − α)s + α

1 ∂ C

¯

=

− ωC.

S_C∂y¯²

ꢄꢄ

ꢅ

cos A₆y¯

2 cos A₆

sin A₆y¯

1

¯

θ(y, s) =

−

(33)

2 sin A₆s + λ²

2.3 Exact solutions

Equations (22) and (29) provide the drug’s concentration in the ﬂowing

blood in the carotid artery as follows:

Here we assume the following as the arbitrary solutions of equations (9),

(16), (17) and (18),

1

¯

u¯ = F(y)

ꢄ

ꢅ

,

(19)

(20)

(21)

(22)

s + λ²

cosh A₉y¯

2 cosh A₉

sinh A₉y¯

1

¯

C(y, s) =

−

.

2 sinh A₉s + λ²

1

¯

θ = H(y)

,

s + λ²

1

¯

v¯ = G(y)

,

(34)

s + λ²

1

Equations (31) through (34) yield the inverse Laplace transform.

Using Mathcad software, we simulated the given solutions using

Gaver-Stehfest’s algorithm, and the results are shown graphically in the

next section.

¯

C = I(y)

,

s + λ²

then the boundary conditions in eqns. (5) reduce to;

H = 1, I = 1, F = 1, at y = −1,

H → 0, I → 0, F → 0, at y = 1.

3 Results and Discussion

(23)

To get the ﬂow information, we simulated the solutions of equations

(31), (33), and (34) using Mathcad software. The inﬂuence of the

fractional-parameter (α) on velocity, temperature and blood concentration

are displayed graphically and discussed. Axial ﬂuid velocity, temperature

distribution, and concentration are explored as functions of several

dimensionless factors, including: slanted (inclined) magnetic ﬁeld

parameter (M), radiation parameter (R), fractional parameter (α), heat

source parameter (S), and Schmidt number (S_C). In all the dimensionless

parameter calculations, we vary the value of the fractional parameter (α),

but we maintain other values constant, such as, t = 1, S_C= 0.5, ω = 0.5,

S = 1, P_r= 2, K = 2, R = 0.5, h = 0.5, β = 0.5, φ = 30^◦.

As a result, the following are the simpliﬁed governing equations of

motions with arbitrary solutions:

2

¯

d F

0

¯

− A₂F = A₃− gβH − gβ I,

(24)

dy²

2

¯

d H

2

¯

+ A₆H = 0,

(25)

(26)

(27)

dy²

G = A₁(constant),

2

¯

d I

dy²

2

¯

− A₉I = 0.

3.1 Velocity Profile

Equation (23)’s boundary conditions are used to solve equations (24)

through (27) and the following solutions are obtained:

Consequently, the magnetic ﬁeld always has a greater inﬂuence on the

blood velocity proﬁle. The application of the magnetic ﬁeld to the system,

as shown in Figure 3, increases the Lorentz force, a resistive force that

primarily restricts the ﬂow of ﬂuid Bunonyo and Ebiwareme (2023) and

Vardanyan (1973). For fractional order (α = 0.4), as the magnetic ﬁeld

parameter’s strength increases, as seen in Figure 3(a), the blood velocity

reduces sharply, whereas it declines gradually for (α = 1) as shown in

Figure 3b.

ꢄ

ꢅ

cos A₆y¯

2 cos A₆

sin A₆y¯

2 sin A₆

¯

H =

−

(28)

(29)

cosh A₉y¯

2 cosh A₉

sinh A₉y¯

2 sinh A₉

¯

I =

−

Abdullahi et.al (2026)

21

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

However, for the fractional order parameter in Figure 6a and 6b, the ﬂow

velocity vanishes between angles of (80^◦to 85^◦).

Because several of the related parameters’ values have changed, the

graphs in Figure 6c and 6d behave very diﬀerently from one another.

Using Figure 6c as an example, y = 0.004, p = 4, and the fractional

parameter (α

=

0.2), whereas Figure 6d also includes the fractional

parameter (α = 0.4) and y = 0.004, p = 3.

(a)

(b)

Figure 3: Axial velocity proﬁle for diﬀerent values of magnetic parameter: (a)

α = 0.4 (b) α = 1, ω = 0.5, S = 1, P = 2, K = 2, R = 0.5, t = 1, h = 0.5,

r

β = 0.5, φ = 30^◦, Sc = 0.5.

For therapeutic purposes and treatment procedures related to

atherosclerosis, bone fractures, controlled tissue damage, and malignant

tumors, to mention a few, a regulated magnetic ﬁeld can therefore be a

useful tool (Imoro et al., 2024). For both the fractional parameter (α = 0.4)

and the integer order model blood ﬂow (α = 1), Figure 4 shows the

variation in blood ﬂow at diﬀerent heat source parameter (S) values. It is

clear that an increase in the heat source has an impact on blood velocity

and the fractional ﬂuid parameter (α = 0.4) (see Figure 4a). As seen in

Figure 4b, the axial velocity does, however, drop symmetrically as the

heat source parameter increases.

(a)

(b)

(c)

(d)

Figure 6: Axial velocity proﬁle for various angles of inclination of the magnetic ﬁeld:

(a) α = 0.4 (b) α = 1 (c) y = 0.004, p = 4 and α = 0.2 (d) y = 0.004, p = 3, α = 0.4,

ω = 0.5, S = 1, P = 2, K = 2, M = 0.5, t = 1, h = 0.5, β = 0.5, φ = 30^◦, Sc = 0.5.

r

(a)

(b)

Figure 4: Proﬁle of axial velocity for various heat source parameter: (a) α = 0.4 (b)

The blood ﬂow velocity proﬁle at two independent times, t

=

0.01

α = 1, ω = 0.5, M = 0.5, P = 2, K = 2, R = 0.5, t = 1, h = 0.5, β = 0.5, φ = 30^◦,

r

and 0.5, is shown in Figure 7 with ﬁve diﬀerent values of the fractional

parameter (α = 0.2, 0.4, 0.6, 0.8, and 1). It has been observed that the

fractional parameter (α) plays a critical role in regulating blood velocity.

The fractional derivative ﬂuid velocity initially moves faster than the

integer order ﬂuid model when the time is relatively small (t = 0.01).

However, for a longer period of time (t = 0.5), the reverse behavior

is seen, that is, ﬂuids with integer order have a faster velocity than

those with fractional order parameters. Naturally, this results from the

system’s stability, which can improve over longer timescales. For both

fractional order derivative ﬂuid models and integer order derivative ﬂuid

models, it is often observed that blood velocity increases with increasing

time t. Figure 7a shows the evolution of the primary velocity proﬁle,

showing how the ﬂow develops into fully formed Poiseuille ﬂow, which

is distinguished by the typical parabolic proﬁle. Figure 7b clearly depicts

the velocity proﬁle at the fork section for various values of the fractional

parameter. It is signiﬁcant to note that when (α = 1), a zone of sluggish

ﬂow appears along the outside wall and gets worse as time goes on, as

was previously observed by Gade et al. (2026).

Sc = 0.5.

Figure 5 displays the velocity distribution based on diﬀerent thermal

radiation parameters (R). The blood velocity increases as the radiation

parameter (R) increases, as indicated by both the fractional parameter

(α

=

0.4) and the classical order parameter (α

=

1). Remarkably,

comparable results for a related ﬂuid model were discussed in Tabi et al.

(2017). According to Yakubu et al. (2022), heat radiation possesses the

capability to modify the eﬀective viscosity of ﬂuids, hence potentially

causing an indirect inﬂuence on the velocity proﬁle.

(a)

(b)

Figure 5: Proﬁle of axial velocity for various heat source parameter: (a) α = 0.4 (b)

α = 1, ω = 0.5, M = 0.5, P = 2, K = 2, R = 0.5, t = 1, h = 0.5, β = 0.5, φ = 30^◦,

r

Sc = 0.5.

(a)

(b)

The applied magnetic ﬁeld parameter for various tilted values is

displayed in Figure 6. Blood ﬂow is reduced over the aﬀected area when

the applied magnetic ﬁeld’s angle of inclination is increased for both the

fractional order parameter and the classical order parameter (α = 1).

Figure 7: Axial velocity proﬁle for diﬀerent values of α at: (a) t = 0.01 (b) t = 0.5

ω = 0.5, S = 1, P = 2, K = 2, M = 0.5, t = 1, h = 0.5, β = 0.5, φ = 30^◦, Sc = 0.5

r

Abdullahi et.al (2026)

22

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

3.2 Temperature profile

It’s noteworthy to see in Figure 9b that the temperature at the middle line

of the channel decreases as the heat source’s values rise. The temperature

exhibits oscillating behavior for diﬀerent amounts of the heat source in

Figures 9c and 9d. As the values of the heat source rise, the temperature

is maximum at the center, decreases, and ﬁnally approaches zero at the

artery walls. The diﬀerent values of the fractional parameter really cause

a shift in the temperature distribution, as shown in Figure 10. It illustrates

how temperature increases with increasing fractional parameter . This

implies that the fractional order ﬂuid model’s temperature distribution is

more faster and higher over a longer period of time, which is what causes

the variation shown in Figures 10a and 10b, as mentioned earlier. The

temperature gradually decreases toward the artery’s axis in Figures 10c

and 10d, eventually tending to align with the axis.

Temperature proﬁles for diﬀerent radiation parameters (R), fractional

parameters , and heat source parameters (S) are shown in Figures 8–??.

The temperature change for diﬀerent values of the radiation parameter

R, as shown in Figure 8. It is evident that when the thermal radiation

increases, temperature increases for both fractional and integer order

derivatives. The temperature varies near the center line for both the

integer order and the fractional order derivative, as shown in Figure 8a

and 8b. Consequently, it is more visible in the graphs of Figure 8b.

(a)

(b)

Figure 8: Temperature proﬁle for diﬀerent values of radiation parameter: (a) α = 0.4

(b) α = 1, ω = 0.5, S = 1, P = 2, K = 2, M = 0.5, t = 1, h = 0.5, β = 0.5, φ = 30^◦,

(a)

(b)

r

Sc = 0.5

During hyperthermia, the temperature distribution is very important. It

is commonly recognized that hyperthermia results from a breakdown

in thermoregulation, which takes place when the body absorbs heat

from outside sources like radiation or a body temperature that is being

generated or absorbed. When a person has hyperthermia, the blood’s

internal temperature increases without damaging the tissues around

the blood vessel. We have not taken into account the temperature

exchange at the artery wall to account for this, meaning that the wall’s

temperature is zero. In light of this, the blood temperature in the current

model is low at the artery wall and high at the midline for classical

ﬂuid. Numerous theoretical and experimental studies for Newtonian

and non-Newtonian ﬂuids of integer order reported similar phenomena,

for example in (Ramesh & Devakar, 2015). Similar to radiation, the heat

source (S) another crucial factor, has a large impact on the bloodstream’s

temperature distribution. More mitochondria per cell increase the

thermal activity involved with the heat production process, as seen

in Figure 9, which raises the system’s temperature. The heat source

improves the temperature distribution and supplies more heat to the

blood ﬂow system even though the wall temperature must remain zero in

order to meet the boundary conditions. The temperature distribution at

the channel walls, which decreases and becomes more ﬂattened toward

the channel’s center line when the heat source is increased as shown in

Figure 9a, which is ampliﬁed to maintain a constant ﬂow rate, as seen in

Figure 9a.

(c)

(d)

Figure 10: Temperature proﬁle for diﬀerent values of α at: (a) t = 0.05, y = 0.004,

p = 1 and α = 1, t = 0.1 (b) t = 0.25, ω = 0.5, S = 1, P = 2, K = 2, M = 0.5, t = 1,

r

h = 0.5, β = 0.5, φ = 30^◦, Sc = 0.5, y = 0.004, p = 3, α = 0.4. (c) y = 0.004, p = 2

and α = 1, t = 0.1 (d) y = 0.004, p = 10, α = 0.01, t = 0.1.

3.3 Concentration profile

The concentration proﬁle for diﬀerent values of fractional order

parameter (α) , Schmidt number (Cs) , and chemical process (ω) is

shown in Figures 11 to 13. There is a relationship between the blood

concentration and the quantity of blood cells ﬂoating in the plasma. Red

blood cells (RBCs) are important blood cells because of their size and

density in the bloodstream. RBCs assembled at the center of the vessel,

where there is a greater concentration of solutes, due to their revolving

nature. However, because the oﬀ-axis zone is an area predominantly

represented by cells that carry plasma, the solute concentration there

decreases to a minimum. This observation is displayed in all of the

concentration graphs in this section. The fractional model ﬂuid in Figure

11 reaches a greater concentration more quickly than the integer order

ﬂuid (Imoro et al., 2024).

(a)

(b)

(c)

(d)

(a)

(b)

Figure 9: Temperature proﬁle for diﬀerent values of heat source parameter: (a)

α = 0.4 (b) α = 1, ω = 0.5, R = 0.5, P = 2, K = 2, M = 0.5, t = 1, h = 0.5,

Figure 11: Concentration proﬁle for diﬀerent values of α at: (a) t = 0.1 (b) t = 0.5

r

β = 0.5, φ = 30^◦, Sc = 0.5

Sc = 0.5, S = 1, R = 0.5, P = 2, K = 2, M = 0.5, t = 1, h = 0.5, β = 0.5, φ = 30^◦

r

Abdullahi et.al (2026)

23

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

This is because a fractional order derivative that restricts ﬂuid ﬂow

is included in the model. The Schmidt number exhibits the opposite

pattern. The blood cells show an additional force of the temperature

gradient in the presence of the Schmidt number, as seen in Figure 12,

which further increases the concentration. Therefore, lower Schmidt

number values, for example in industrial applications, physically

represent hydrogen gas as the species diﬀusing (Sademaki et al., 2026).

4 Conclusion

Currently, a fractional-order model of the magneto hydrodynamic blood

ﬂow via a bifurcated artery under the inﬂuence of thermal radiation,

a slanted magnetic ﬁeld, and a heat source during tumor treatment is

being developed. The Laplace transform and the combined methods

of indeterminate coeﬃcients were used to solve the mathematical

models. The fractional order parameter has a major eﬀect on the

blood velocity proﬁles, concentration, and temperature distribution.

It has been noted that ﬂuids with fractional order can occasionally

move faster than those with integer order. Fractional model ﬂuid

ﬂow is slower than integer-order ﬂuid ﬂow over longer dimensionless

durations. The impact of ﬂuid velocity is demonstrated by the fact

that the rate of increase in ﬂuid velocity is slower at larger levels of

the magnetic ﬁeld parameter. As the chemical reaction parameter rises,

the blood ﬂow concentration falls. The blood concentration rises as the

Schmidt number rises. As the fractional parameter and the heat source

increase, the blood ﬂow’s dimensionless temperature rises, which also

aﬀects the radiation parameter. We noted that the outcomes will be

intriguing to comprehend and evaluate throughout cancer therapy using

hyperthermia. Additionally, it will be useful in understanding the drug

particle concentration phenomena for applications and administrations

involving drug delivery. Our research’s ﬁndings should serve as a

foundation for the study of increasingly sophisticated blood ﬂow models

and also serve as a basis for in vitro and in vivo testing, particularly

in the application areas like medicine, biomedical engineering, biology,

pathology, and other related domains.

(a)

(b)

Figure 12: Concentration proﬁle for diﬀerent values of Schmidt number: (a) α = 0.4

(b) α = 1, ω = 0.5, S = 1, R = 0.5, P = 2, K = 2, M = 0.5, t = 1, h = 0.5, β = 0.5,

r

φ = 30^◦

Figure 12 illustrates how the species’ chemical molecular diﬀusivity

decreases dramatically with increasing Schmidt number (Sc) , making

it easier for the species to enter the ﬂow ﬁeld and raising the mass

transfer function. Higher Schmidt number compounds can enhance

mass transfer and dispersion properties in the bloodstream, especially

for pharmaceutical diﬀusion in pulse blood ﬂow. The amplitude of the

blood ﬂow concentration is larger for the integer order derivative. As

demonstrated in Figure 13a and 13b, this phenomena is clearly seen along

the ﬂow axis (0 ≤ y ≤ 0.5) and slowly declines in the region (0 ≤ y ≤ 1)

for both fractional and integer order derivatives, respectively. As can be

observed from all of the graphs in Figure 13, the blood ﬂow decreases

along the distensible tube’s length where the graphs begin to ﬂuctuate

because of the size of the chemical reaction parameter’s peak value

(pressure gradient) (see, Abdul-Wahab & Al-Saif, 2024). Furthermore,

it has been demonstrated that the variation is more pronounced in the

larger section of the artery wall, permitting the ﬂow to pass without

producing a perceptible pressure gradient. Nevertheless, the substantial

pressure gradient is usually required to maintain a consistent ﬂow rate as

it passes through the constrictions in the artery.

Acknowledgments

This work was supported by Tertiary Education Trust Fund (TETFund)

Ref. No. TETF/ DR&D/CE /UNI /BAUCHI /IBR /2025/ VOL.1.

Therefore, the authors gratefully acknowledged the ﬁnancial support of

the TETFUND. The authors would like to express their sincere gratitude

to the handily editor and the reviewers for their helpful and informative

comments, which have enhanced the manuscript.

Declaration of Generative AI

The authors declare that they do not used generative AI in the scientiﬁc

writing.

References

Abdulhameed, M., Babagana, B., Markus, S., Yakubu, D. G., & Adamu,

G. T. (2023). The eﬀects of fractional relaxation time and

magnetic ﬁeld on blood ﬂow through arteries along with

nanoparticles. Defect and Diﬀusion Forum, 424, 59–76.

Abdulhameed, M., Vieru, D.,

&

Roslan, R. (2017). Modeling

electro-magnetohydrodynamic thermo-ﬂuidic transport of

bioﬂuids with new trend of fractional derivative without

singular kernel. Physica A, 484, 233–252.

(a)

(b)

Abdul-Wahab, M. S., & Al-Saif, A. S. J. A. (2024). A new method for

studying blood ﬂow through stenotic artery in the presence of

a magnetic ﬁeld. Intern. J. Appl. Comput. Math., 10(49). http s :

//doi.org/10.1007/s40819-024-01684-x

Akbar, N. S., & Butt, A. W. (2017). Entropy generation analysis in

convective ferromagnetic nano blood ﬂow through a composite

stenosed arteries with permeable wall. Commun. Theor. Phys.,

67, 554–560.

Ali, F., Sheikh, N. A., Khan, I., & Saqib, M. (2017). Magnetic ﬁeld eﬀect on

blood ﬂow of casson ﬂuid in axisymmetric cylindrical tube: A

fractional model. J. Magn. Magn. Mater., 423, 327–336.

Atangana, A., & Baleanu, D. (2016). New fractional derivative with

non-local and non-singular kernel: Theory and application to

heat transfer model. Thermal Science, 21(2), 761–766. http s ://

doi.org/10.2298/TSCI160517180A

(c)

(d)

Figure 13: Concentration proﬁle for diﬀerent values of the chemical reaction

parameter: Sc = 0.5, S = 1, P = 2, K = 2, M = 0.5, t = 1, h = 0.5, β = 0.5,

φ = 30^◦, R = 0.5 (a) α = 0.4 (b) ^rα = 1.

Abdullahi et.al (2026)

24

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

Bansi, C. D. K., Tabi, C. B., Motsumi, T. G., & Mohamadoud, A. (2018).

Fractional blood ﬂow in oscillatory arteries with thermal

radiation and magnetic ﬁeld eﬀects. J. Magn. Magn. Mater., 456,

38–45.

Bhargava, R., Rawat, S., Takhar, H. S., & Bég, O. A. (2007). Pulsatile

magnetobioﬂuid ﬂow and mass transfer in a non-darcian

porous medium channels. Meccanica, 42, 247–262.

Bhatti, M. M., & Lu, D. Q. (2019). Analytical study of the head-on collision

process between hydroelastic solitary waves in the presence of

a uniform current. Symmetry, 11, 333.

of magnetic ﬁeld for drug delivery application. J. Magn. Magn.

Mater., 442, 319–328.

Nagarani, P., Sarojamma, G., & Jayaraman, G. (2006). Exact analysis

of unsteady convective diﬀusion in casson ﬂuid ﬂow in an

annulus-Application to catheterized artery. Acta Mechanica,

187, 189–202.

Prasad, K. V., Vaidya, H., Choudhari, R., Tripathi, D., Karanth, S., &

Hanumantha. (2025). Advancing blood ﬂow in stenotic arteries

through magnetohydrodynamic peristaltic motion of hybrid

nanoparticles. Chinese Journal of Physics, 96, 1144–1163.

Ramesh, K., & Devakar, M. (2015). Magneto hydrodynamic peristaltic

transport of couple stress ﬂuid through porous medium in an

inclined asymmetric channel with heat transfer. J. Magn. Magn.

Mater., 394, 335–348.

Sademaki, L. J., Reddy, B. P., & Matao, P. M. (2026). Dissipative and

radiative consequences on diﬀusional reactive MHD nanoﬂuid

ﬂow over an inclined vertical cone in a porous medium with

reactive species: FEM study. Partial Diﬀ. Equat. Appl. Math., 18,

101365.

Bhatti, M. M., Zeeshan, A., & Ellahi, R. (2016). Heat transfer analysis

on peristaltically induced motion of particle-ﬂuid suspension

with variable viscosity: Clot blood model. Comput. Math. Prog.

Biomed., 137, 115–124.

Bhatti, M. M., Zeeshan, A., Ellahi, R.,

&

Shit, G. C. (2018).

Mathematical modeling of heat and mass transfer eﬀects on

MHD peristaltic propulsion of two-phase ﬂow through a

Darcy-Brinkman-Forchheimer porous medium. Adv. Powder

Tech., 29, 1189–1197.

Bunonyo, K. W., & Ebiwareme, L. (2023). Mathematical analysis of a

magnetic and conducting ﬂuid ﬂow through blood vessel

along with an inclination and chemical radiation. European J.

Theoretical and Applied Sciences, 1(6), 3–15.

Caputo, M., & Fabrizio, M. (2015). A new deﬁnition of fractional

derivative without singular kernel. Progress in Fractional

Diﬀerentiation and Applications, 1(2), 73–85.

Caro, C. G., Pedley, T. J., Schroter, R. C., & Seed, W. A. (2011). The mechanics

of the circulation. Cambridge University Press.

Samko, S. G., Kilbas, A. A., & Marichev, O. I. (1993). Fractional integrals

and derivatives: Theory and applications. Gordon; Breach Science

Publishers.

Shah, N. A., Vieru, D., & Fetecau, C. (2016). Eﬀects of the fractional order

and magnetic ﬁeld on the blood ﬂow in cylindrical domains. J.

Magn. Magn. Mater., 409, 10–19.

Shaw, S., & Murthy, P. V. S. N. (2010). Magnetic drug targeting in

the permeable blood vessel - The eﬀect of blood rheology. J.

Nanotechnol. Eng. Med., 1(2), 021001–11.

Chaturani, P., & Palanisamy, V. (1990). Casson ﬂuid model for pulsatile

ﬂow of blood under periodic body acceleration. Biorheology,

27(5), 619–630.

Chinyoka, T., & Makinde, O. D. (2014). Computational dynamics of

arterial blood ﬂow in the presence of magnetic ﬁeld and

thermal radiation therapy. Adv. Math. Phys., 2014, 915640.

Dash, R. K., Mehta, K. N., & Jayaraman, G. (1996). Casson ﬂuid ﬂow in

a pipe ﬁlled with homogeneous porous medium. Int. J. Engg.

Sci., 34, 1146–1156.

Gade, M. R., Kalakuntla, S. R., Adigoppula, R., & Itikela, S. (2026).

Prediction of micropolar ﬂuid ﬂow characteristics in a stenosed

bifurcated artery using feed-forward neural networks trained

by the Levenberg Marquardt Algorithm. Partial Diﬀ. Equat.

Appl. Math., 18, 101366.

Ghasemi, S. E., Hatami, M., Hatami, J., Sahebi, S. A. R., & Ganji, D. D.

(2016). An eﬃcient approach to study the pulsatile blood ﬂow

in femoral and coronary arteries by diﬀerential quadrature

method. Physica A, 443, 406–414.

Ghasemi, S. E., Hatami, M., Sarokolaie, A. K., & Ganji, D. D. (2015).

Study on blood ﬂow containing nanoparticles through porous

arteries in presence of magnetic ﬁeld using analytical methods.

Physica E, 70, 146–156.

Hayat, T., Asad, S., & Alsaedi, A. (2016). Flow of casson ﬂuid with

nanoparticles. Appl. Math. Mech., 37(4), 479–470.

He, S., Fataf, N. A. A., Banerjee, S., & Sun, K. (2019). Complexity

in the muscular blood vessel model with variable fractional

derivative and external disturbances. Physica A, 526, 120904.

Imoro, I., Etwire, C. J., & Musah, R. (2024). MHD ﬂow of blood-based

hybrid nanoﬂuid through a stenosed artery with thermal

radiation eﬀect. Case Studies in Thermal Engin., 59, 104418.

Kumar, D., Satyanarayana, B., Rajesh, K., Narendra, D., & Sanjeev, K.

(2021). Application of heat source and chemical reaction

in magnetohydrodynamic blood ﬂow through permeable

bifurcated arteries with inclined magnetic ﬁeld in tumor

treatments [1-13]. Results in Applied Mathematics, 10, 100151.

Liepsch, D. (1986). Flow in tubes and arteries - A comparison. Biorheology,

23, 395–433.

Shit, G. C., & Majee, S. (2015). Pulsatile ﬂow of blood and heat

transfer with variable viscosity under magnetic and vibration

environment. J. Magn. Magn. Mater., 388, 106–115.

Shit, G. C., & Roy, M. (2015). Eﬀect of slip velocity on peristaltic transport

of a magneto-micropolar ﬂuid through a porous non-uniform

channel. Int. J. App. Compt. Math., 1, 121–141.

Sinha, A., & Shit, G. C. (2015). Electromagnetohydrodynamic ﬂow of

blood and heat transfer in a capillary with thermal radiation. J.

Magn. Magn. Mater., 378, 143–151.

Srivastava, L., & Srivastava, V. (1984). Peristaltic transport of blood:

Casson model-11. J. Biomech., 17(11), 821–829.

Sud, V. K., & Sekhon, G. S. (1984). Blood ﬂow subject to a single cycle of

body acceleration. Bull. Math. Biol, 46, 937–949.

Syed, M. H., Mustansar, S. H. S., Hi az, A., Nazar, T., Wasim, J.,

Mohamed, R. E., et al. (2026). Thermal characteristics of

magnetic blood-based hexa-hybrid nanoﬂuids in stenotic

arteries with heat source/sink by applying Caputo-Fabrizio

fractional derivatives [In Press]. Results in Surfaces and Interfaces.

Tabi, C. B., Motsumi, T. G., Kamdem, C. D. B., & Mohamadou, A. (2017).

Nonlinear excitations of blood ﬂow in large vessels under

thermal radiations and uniform magnetic ﬁeld. Commun. Nonl.

Sci. Numer. Simul., 49, 1–8.

Tzirtzilakis, E. E. (2005). A mathematical model for blood ﬂow in

magnetic ﬁeld. Phys. Fluids, 17(7), 077103.

Vardanyan, V. A. (1973). Eﬀect of magnetic ﬁeld on blood ﬂow. Bioﬁzika,

18, 491–496.

Venkatesan, J., Sankar, D., Hemalatha, K.,

&

Yatim, Y. (2013).

Mathematical analysis of casson ﬂuid model for blood

rheology in stenosed narrow arteries. J. Appl. Math., 2013, 1–11.

Yakubu, D. G., Abdulhameed, M., Adamu, G. T., & Kwami, A. M. (2020).

A study of fractional relaxation time on blood ﬂow in arteries

with magnetic radiation eﬀects. Diﬀ. Found., 26, 126–144.

Yakubu, D. G., Abdulhameed, M., Adamu, G. T., Roslan, R., Issakhov, A.,

Rahimi-Gorji, M., & Bakouri, M. (2021). Towards the exact

solution of Burger’s ﬂuid ﬂow through arteries with fractional

time derivative magnetic ﬁeld and thermal radiation eﬀects. J.

Proce. Mech. Eng., 235, 1618–1627.

MacDonald, D. A. (1979). On steady ﬂow through modeled vascular

stenosis. J. Biomech., 12(1), 13–20.

Majee, S., & Shit, G. C. (2017). Numerical investigation of MHD ﬂow of

blood and heat transfer in a stenosed arterial segment. J. Magn.

Magn. Mater., 424, 137–147.

Misra, J. C., & Shit, G. C. (2009). Flow of a biomagnatic visco-elastic ﬂuid

in a channel with stretching walls. J. Appl. Mech., 76(6), 061006.

Mondal, A., & Shit, G. C. (2017). Transport of magneto-nanoparticles

drugging electro-osmotic ﬂow in a micro-tube in the presence

Yakubu, D. G., Mohammed, A., Garba, T. A., Usman, H., & Muhammad,

L. K. (2022). Construction of the exact solution of blood ﬂow

of Oldroyd-B ﬂuids through arteries with eﬀects of fractional

derivative magnetic ﬁeld and heat transfer. J. Mech. Med. Biol.,

22(10), 2250068.

Zeeshan, A., Bhatti, M. M., Akbar, N. S.,

&

Sajjad, Y. (2017).

Hydromagnetic blood ﬂow of Sisko-ﬂuid in a non-uniform

channel induced by peristaltic wave. Commun. Theor. Phys., 68,

103–110.

Abdullahi et.al (2026)

25

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

Appendix

1

k

s

τp_r

Rτp_r+ 1

A₁= 1, A₂= M²sin²φ −

−

,

A₃= h(s + λ²), A₄=

,

(1 − α)s + α

ꢄ

ꢅ

p

τp_r

Rτp_r+ 1

S

τp_r

S

A₃

A

A₃= h(s + λ²), A₄=

,

A₅=

−

(s + λ²), A₆= A₄A₅, A₇=

,

(1 − α)s + α

q

−gβ

gβ

A₉= Sc(ω + s((1 − α)s + α)), A₁₀

=

,

A₁₁

=

2(A²₉− A₂) cosh A₉

2(A²₉− A₂) sinh A₉

Abdullahi et.al (2026)

26

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

ART ICLE

Prevalence and Determinant Factors of

Malaria Infection among Patients Attending

Gimbichu Primary Hospital, Soro District,

Central Ethiopia, Ethiopia

ARTICLE INFO

Volume 7(1), 2026

https:

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

Melese Birmeka ^1,∗�, Gebremedhin Gebrezgabiher², Tekleweyni

Asayehegn³, and Mohammed Kasso⁴

ARTICLE HISTORY

Received: 12 June, 2025

¹Department of Biology, Hawassa University, Hawassa, P. O. Box 05,

Accepted: 21 January, 2026

Published Online: 10 June, 2026

²Department of Veterinary Medicine, College of Veterinary Medicine and Animal Sciences, Samara

University, P. O. Box 132, Samara, Afar, Ethiopia

³Department of Aquatic Sciences, Fisheries and Aquaculture, Hawassa University, Hawassa, P. O. Box 05,

⁴Department of Biology, Hawassa University, Hawassa, P. O. Box 05.

CITATION

Birmeka M. et.al (2026). Prevalence and

Determinant Factors of Malaria Infection

among Patients Attending Gimbichu

Primary Hospital, Soro District, Central

Ethiopia, Ethiopia. East African Journal

of Biophysical and Computational

^∗Corresponding author: melesebirmeka@yahoo.com

Abstract

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

doi.org/10.4314/eajbcs.v7i1.3S.27-33

In the world, particularly in Ethiopia, malaria has a great influence on human health and economy.

This study intended to determine the prevalence, trends and associated risk factors of malaria patients

visiting Gimbichu Primary Hospital, Ethiopia. To assess the trend and parasitological examination,

a hospital-based cross-sectional study was carried out. To determine factors that significantly

associated with infection, a bivariate and multivariable logistic regression analyses were performed

with statistical significance set at p < 0.05. The findings of the study revealed the overall malaria

prevalence of 72.4% among suspected patients. The study also revealed that males (AOR = 3.5, 95%

CI; 1.5 - 3.8, p<0.001), individuals under five years (AOR = 2.8, 95% CI: 1.13 - 2.2), 5-20 years (AOR

= 1.75, 95% CI: 1.1 - 1.91) and 21-45 years (AOR = 1.65, 95% CI: 1.01 - 1.49) were at higher risk.

Additionally, study participants living close to mosquito breeding sites (AOR = 2.54, 95% CI: 2.53 -

4.14), rural (AOR = 2.13, 95% CI: 1.01 - 2.6), houses with thatch roof (AOR = 1.43, 95% CI: 1.01-2.30),

not using bed nets (AOR = 1.51, 95% CI: 2.01 - 4.1), homes with wall openings (AOR = 1.6, 95%

CI: 1.13 - 2.57), monthly income of less than 1,000 Ethiopian Birr (AOR = 2.93, 95% CI: 1.3 - 4.6),

and pregnant women (AOR = 1.6, 95% CI: 1.13 - 2.57) had maximum risk for malaria infection. The

analysis from the retrospective data showed the overall decreasing trend in malaria infection rates,

despite the fluctuations recorded between 2015 - 2021. The study indicates that malaria is persistent

and a significant public health challenge which is driven by a complex interrelationship of demographic,

social and environmental factors. Plasmodium vivax infection is the most prevalent species known to

cause malaria in the study area. These findings necessitate targeted interventions focusing on housing

improvements, economic support, and vector control measures.

OPEN ACCESS

This work is licensed under the Creative

Commons open access license (CC

BY-NC 4.0).

East African Journal of Biophysical and

Computational Sciences (EAJBCS) is

already indexed on known databases

like AJOL, DOAJ, CABI ABSTRACTS and

FAO AGRIS.

Keywords: Determinants; Malaria; Prevalence; Soro District, Trend

both mother and the child. Similarly, children under ﬁve are at higher

risk because their immune systems are not fully developed, with a child

dying of malaria every 45 seconds worldwide.

1 Introduction

Malaria is a contagious disease that is caused by parasitic proatozoa

(Ferede et al., 2013). It has great impact on world population health and

economy (Baird, 2013). It is dominantly caused by Plasmodium falciparum

and Plasmodium vivax. Out of the two species P. falciparum is the stronger

pathogenic species that causes most deaths by malaria diseases at world

scale, accounting for more than 90% of the world malaria mortality

(Baird, 2013; Ferede et al., 2013). The malaria disease is a severe disease

particularly in children and pregnant women. Pregnant women are more

susceptible due to decreased immunity during pregnancy, endangering

As an infectious vector-borne disease, malaria continues to be a key public

health challenge in the country, with transmission patterns varying across

regions depending on climatic conditions, rainfall, and altitude. In

Ethiopia, malaria is known to be dominantly caused by P. falciparum (60%)

and P. vivax (40%) (FMOH (Federal Ministry of Health), 2018). In the

country, about 75% of areas located below 2,000 meters above sea level are

susceptible to malaria epidemics and the persistent risk of transmission

(Girum et al., 2019). Annually, approximately 4,782,000 reported cases

Birmeka M., et.al (2026)

27

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

and related deaths, with morbidity and mortality increasing markedly

during epidemic periods were recorded in the country (Alemayehu et al.,

2014). Of this, the large-scale epidemics tend to occur every ﬁve to eight

years although smaller, localized outbreaks are reported annually (Tsige

et al., 2011). An estimated 68% of Ethiopians, a country with a population

of over 100 million people, are at risk of contracting malaria (WHO (World

Health Organization), 2016).

Exclusion criteria: Malaria suspected patients who were not to give

consent for participation in this study.

Sampling and Sample Size Determination All patients suspected of

having malaria were consecutively selected during their visits to the

outpatient department of Gimbichu Primary Hospital till the compulsory

sample size was achieved. The sample size was estimated using Daniel’s

formula (Daniel, 2004).

z²(1 − p)

In Ethiopia, malaria transmission shows a considerable variation across

seasons, years, and geographic settings. The high impact of malaria is

particularly pronounced in rural areas (Donnelly et al., 2005), largely

due to proximity to mosquito breeding sites, limited coverage of control

interventions, widespread poverty, low literacy levels, land-use practices

and poor housing conditions (Stratton et al., 2008). Exceptionally, a

yearly based transmission observation is evident in the southwestern

lowland regions bordering neighboring countries (Zhou et al., 2016).

The communities with lower socioeconomic status are known to be

disproportionately aﬀected (WHO (World Health Organization), 2012),

because it limits access to medical care and preventative measures like

indoor spraying, bed nets treatment and eﬃcient antimalarial therapy

(Yamamoto et al., 2010). The drug-resistant strains of P. falciparum and

P. vivax have also emerged and spread, posing a signiﬁcant challenge to

the control of malaria which contributed to the recent increase in malaria

cases in the nation (Yarcho, 2010).

N =

(1)

d²

Where - p = 50%, because of the absence of previous malaria prevalence

studies in the area, - d = margin of error at 5% and - z = 1.96 at 95% CI

Consequently, the sample size was determined to be 384.

2.5 Data Collection

The structured pretested questionnaires were used to collect information

on socio-demographic and economic status of study participants. Blood

sample collection was done by ﬁnger prick by healthcare professional,

and on the same slide both thick and thin blood smears were prepared.

Throughout the data collection process, continuous monitoring and

supervision were maintained. The activities performed by laboratory

technicians, interviewers and nurses were closely overseen. Additionally,

retrospective data spanning for seven years (2015 -2021) was retrieved

from hospital registration records.

The Government of Ethiopian has made signiﬁcant progresses since

2005 in malaria control interventions such as diagnostic testing, rapid

case treatment, and prevention strategies for pregnant women through

intermittent preventive therapy. High eﬀorts also implemented on the

distribution of IRS and ITNs. However, the widespread emergence of

drug resistance in parasites and insecticide resistance in vectors have

obstructed eﬀorts of malaria eradication (Abeku et al., 2015; Tafese

et al., 2018), particularly in the Hadiya Zone of central Ethiopia. This

situation underscores the need for continuous evaluation and monitoring

of malaria control interventions to address existing gaps.

2.6 Data Analysis

Following a completeness check, the data was analyzed by using SPSS

version 24. Logistic regression studies were performed to identify the

relationship between a few possible risk variables and malaria infection.

To determine the existence and strength of a connection, AOR at 95% CI

were calculated; if p < 0.05, statistical signiﬁcance was proclaimed.

2 Materials and Methods

2.1 Study Area

2.7 Ethical Consideration

The Institutional Research Ethics Review Committee of CNCS of Hawassa

University examined and approved the study proposal and ethical

clearance was received (Ref.no. IRB/279/13). Additional, permission

was also granted by the Hadiya Zone Health Department and the Soro

District Primary Hospital. Conﬁdentiality and privacy were strictly

upheld, and participation in the study was entirely voluntary. After

awareness made on the objectives of the study, participants gave their

consent participation. Conﬁdentiality was also maintained.

The study was conducted at Gimbichu Primary Hospital which provides

care for the Soro District in Hadiya Zone, central Ethiopia region. The

district is about 264 kilometers south of the nation’s capital, Addis Ababa.

Soro District is home to a substantial population of 233,015 people, nearly

evenly split between genders (115,825 men and 117,190 women), giving

the hospital a wide and diverse community to serve.

2.2 Study Design and Period

3 Results

An institution-based cross-sectional study was carried out between

October 2022 and January 2023.

3.1 Characteristics of Study Participants

in Retrospective Study of Malaria in

Soro District, 2015 - 2021

2.3 Study Population

All individuals who presented to Gimbichu Primary Hospital with

suspected malaria during the data collection period and satisﬁed the

eligibility requirements were involved in the study population

A total of 65,211 clients were registered in the laboratory logbooks of

Gimbichu Primary Hospital. Of these, 36,132(55.4%) were males and

29,089(44.6%) were females. Between 2015 and 2021, 65,211 blood ﬁlms

were microscopically examined. The majority of the cases were males

accounting 17,813(49.3%). Although malaria prevalence ﬂuctuated from

2016 to 2021, there was an overall decreasing trend. Over the seven year

period, a considerable malaria cases were recorded in the age 15 - 24 years

old (18,718 cases, 28.7%), followed by 5 - 14 years old (15,131 cases, 23.2%).

The lowest number of cases was over 54 years (8,740 cases, 13.4%) (Table

1).

2.4 Eligiblity Criteria

Inclusion criteria: Malaria suspected patients who were consented to

participate in the study.

Birmeka M., et.al (2026)

28

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

Table 1: The social and demographic characteristics of microscopically examined suspected patients in Soro District, 2015 - 2021.

Socio-demographic

variables