Heliyon 10 (2024) e35038

Contents lists available at ScienceDirect

Heliyon
journal homepage: www.cell.com/heliyon

Research article

In-host density-dependent model of high-risk HPV virions, basal

cells, lymphocytes t-cells incorporating functional responses
Elosy Makena a, * , Cyrus Gitonga Ngari b , Patrick Mwangi Kimani a ,
Jeremiah Savali Kilonzi c
a
Department of Mathematics and Statistics, University of Embu, 6-60100 Embu, Kenya
b
Department of Pure and Applied Sciences, Kirinyaga University, 143-10300 Kerugoya, Kenya
c
Department of Mathematics, Meru University of Science and Technology, 972-60200 Meru, Kenya

A R T I C L E I N F O A B S T R A C T

Keywords: Cervical cancer is one of the most common types of cancer and it is caused mostly by high-risk
In-host model Human Papillomavirus (HPV) and continues to spread at an alarming rate. While HPV impacts
Functional responses have been investigated before, there are currently only a scanty number of mathematical models
Stability analysis
that account for HPV′s dynamic role in cervical cancer. The objectives were to develop an in-host
Simulation and reproduction number
density-dependent deterministic model for the dynamics implications of basal cells, virions, and
lymphocytes incorporating immunity and functional responses. Analyze the model using tech­
niques of epidemiological models such as basic reproduction number and simulate the model
using Matlab ODE solver. Six compartments are considered in the model that is; Susceptible cells
(S), Infected cells (I), Precancerous cells (P), Cancerous cells (C), Virions (V), and Lymphocytes
(L). Next generation matrix (NGM), survival function, and characteristic polynomial method were
used to determine the basic reproduction number denoted as R0 . R0 was obtained using three
methods because NGM has some weaknesses hence the need for the other two methods. The
findings from this research indicated that Disease-Free Equilibrium point is locally asymptotically
stable whenever R*0 < 1 and globally asymptotically stable if R*0 ≤ 1 and the Endemic Equilibrium
is globally asymptotically stable if R*0 > 1. The results obtained shows that the progression rate of
precancerous cells to cancerous cells (θ) has the most direct impact on the model. The model was
able to estimate the longevity of a patient as 10 days when (θ) increases by 0.08. The findings of
this research will help healthcare providers, public health authorities, and non-governmental
health groups in creating effective prevention strategies to slow the development of cervical
cancer. More research should be done to determine the exact number of cancerous cells that can
lead to the death of a cervical cancer patient since this paper estimated a proportion of 75%.

1. Introduction

According to Ref. [1], cervical cancer arises from the uncontrolled, invasive proliferation of epithelial cells in the cervix, the area of
the uterus that attaches to the vagina. Human Papillomavirus (HPV) is the main cause of cervical cancer, and in over 90 % of instances,

* Corresponding author.
E-mail addresses: elseakesh@gmail.com (E. Makena), ngaricyrus15@gmail.com (C.G. Ngari), patkim32@gmail.com (P.M. Kimani),
jeremiahkilonzi66@gmail.com (J.S. Kilonzi).

https://doi.org/10.1016/j.heliyon.2024.e35038
Received 16 February 2024; Received in revised form 6 July 2024; Accepted 22 July 2024
Available online 25 July 2024
2405-8440/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4.0/).
E. Makena et al. Heliyon 10 (2024) e35038

HPV infection contributes to the development of cancer [1]. Small double-stranded DNA viruses with a diameter of 52–55 nm are
known as human papillomaviruses (HPV) [2]. There are more than 100 varieties of Human Papillomavirus, but the types that cause
cervical cancer are high-risk HPV 16 and 18 [3,4]. Sexual contact between people can spread the sexually transmitted virus HPV [5].
This is often through a cut, abrasion, or a small tear in the skin and sexual activity without protection [6]. Some of the risk factors of
cervical cancer include tobacco usage, harmful alcohol consumption, overweight, obesity, age, the individual’s sex, and genetic or
inherited factors, exposure to carcinogens in the environment, such as chemicals, radiation, and infectious agents [7] (see Figs. 7–15).
The dynamics of HPV to cervical cancer are clearly shown in Fig. 1 [8].
One of the main causes of morbidity and death worldwide is cervical cancer. According to a 2018 World Health Organization
(WHO) fact sheet, approximately 311,000 women worldwide lost their lives to cervical cancer in 2018, making it the fourth most
common malignancy among women globally and a global health concern [9]. In 2020 the estimated number of deaths worldwide was
342,000. Focusing on women’s health is necessary to achieve Sustainable Development Goals (SDG) targets 3.4 and 5, which are to
lower premature mortality from non-communicable diseases (NCDs) by a third by 2030 compared to 2015 levels, to enhance mental
health and well-being, and to achieve gender equality and empower all women and girls, respectively [10]. The estimated number of
secondary infections triggered by an index case in a community that is fully susceptible is known as the basic reproduction number, or
R0 [11]. Methods used to calculate R0 include the survival function, the next generation method, the Jacobian matrix’s eigenvalues,
the presence of the endemic equilibrium, and the characteristic polynomial’s constant term [12]. A description of a system using
mathematical terminology and instruments is called a mathematical model. In the natural sciences, such as biology and epidemiology,
mathematical models play a crucial role. They support our efforts to discover new information about a system, arrange and interpret
biological data, and ascertain the system’s reaction behavior. Daniel Bernoulli’s smallpox model from the 1760s marked the beginning
of mathematical modeling of infectious illness. Since then, numerous infectious diseases have been simulated through the development
of mathematical models, including influenza, HIV, TB, malaria, and tuberculosis to mention a few. These mathematical models were
designed to address an array of unresolved questions [13].
To provide some insights into this study, this section brings to light relevant literature. In particular, literature that is based on
mathematical modeling of Human Papillomavirus dynamics. Finally, the focus is directed on current research that will attempt to fill
the gap in the literature.
[1]formulated an HPV-induced cervical cancer model with six compartments including, Susceptible cells (S), Infected cells (I),
Cancerous cells (C), Dendritic cell population (D), CTL population (T), and HPV (V) to explain how a combination of drugs can treat
cervical cancer [14].developed an SIPVC model consisting of five compartments which are susceptible, infected, precancerous,
cancerous cells, and HPV to the features of HPV infectious and cancer disease.
Other compartments which were used to model the dynamics of HPV to cervical cancer include; SIVPC, SIVPC, SIVPCM, SITR,
SIHPV IUC ICC , S(t)V(t)E(t)Iu (t)C(t)R(t), SIVPC where S-susceptible, I-infected, V-virus, P-precancerous, C- cancerous, T-treated and R-
recovered, E− exposed [6,15–20], respectively.
Those studies missed critical aspects such as; a density-dependent population, functional responses, and the maximum proportion
of cancerous cells that can cause an individual to die and also the trio-interaction of basal cells, virions, and lymphocytes. To achieve
this, a new mathematical model is formulated; an in-host density-dependent deterministic model of basal cells, virions, and
lymphocyte dynamics and their relation to cervical cancer. The study analyzes the model using epidemiological techniques and
simulates it using an inbuilt Matlab ODE solver. The paper is organized as follows: Section 2 is on Model Formulation, Section 3 is on
Model analysis and Section 4 is on Conclusion and recommendations. The findings of this research will help healthcare providers,
public health authorities, and non-governmental health groups in creating effective prevention strategies to slow the development of
cervical cancer.

Fig. 1. Progression of cervical cancer by HPV.

2
E. Makena et al. Heliyon 10 (2024) e35038

2. Model formulation

In this study, first-order nonlinear ordinary differential equations are used to create a logistic deterministic model for the effects of
the Human Papillomavirus. The model of HPV infection and cancer development has considered 6 compartments namely: S-Sus­
ceptible (normal) cells, I-Infected cells, V-Free virus, P-Precancerous cells, C-Cancerous cells, and L- Lymphocytes cells. The model is
analyzed in terms of positivity, equilibrium points, and their stabilities, and the basic reproduction number is determined using the
next-generation matrix method, survival function, and the constant term of a characteristic polynomial.
This study assumes that: the number of cervical epithelial cells remains roughly constant, the epithelium is replaced every 4–5 days
[16], all the epithelium cells are susceptible, the basal cell population grows logistically, with an intrinsic growth rate and a carrying
capacity. There is significance to the recovery from natural immunity and Human Papillomavirus infection is the strongest risk factor
for cervical cancer. Besides the common assumptions of any epidemiological model, the following are additional assumptions.

τ2 S
i. S+BV is Arditi-Ginzburg function predator (virions) density can also influence individual (Prey/Basal Cells) consumption rate, an
effect termed predator dependence. The ratio of prey density to predator density will determine how virions attack Basal cells.
CP2
ii. D2 +P2
is a saturation function in which precancerous cells progress to cancerous cells based on Holling type III response.
EM
iii. F+M is based on Holling type II response on the assumption that the rate of infected basal cell recovery increases with the
immunity level, but saturates at some maximum level M.

The parameters and state variables are summarized in the table below.
Some parameter values in Table 1 were obtained from different literature while others were estimated using the maximum like­
lihood estimation method (see Table 2).
The equations are derived from the model flow chart as follows;
( )
dS N τ2 SV EM
= r1 S 1 − − + I − δ1 S, (1)
dt K1 S + BV F + M
( )
dI N τ2 SV EM
= r1 I 1 − + − βI − I − η1 I − δ1 I, (2)
dt K1 S + BV F+M
( )
dP N CP2
= r1 P 1 − + βI − ϴ 2 − η2 P − δ1 P, (3)
dt K1 D + P2
{ ( ) }( )
dC N CP2 C
= r1 C 1 − +ϴ 2 − η 3 C − δ 1 C 1 − , (4)
dt K1 D + P2 K4

Table 1
Parameter and state variables description.
Parameter Description Value Reference

r1 Mitosis division rate of basal cells 0.8506 [21]

r2 Division rate virions 0.011 Estimated
r3 Division rate lymphocytes 0.902 Estimated
K1 Carrying capacity of basal cells 10, 000, 000 Estimated
K2 Carrying capacity of virions 100, 000 Estimated
K3 Carrying capacity of lymphocytes 50, 000 Estimated
K4 The population of cancerous cells that causes a person to die 6, 500, 000 Estimated
τ2 The rate of infection of susceptible cells by the virions 0.000001 [15]
β Progression rate from infected cells to precancerous cells 0.0082 [17]
ϴ Progression rate from precancerous cells to cancerous cells 0.01 [21]
η1 = η2 = η3 Induced death rate of I, P and C 0.005 [21]
δ1 The apoptosis rate of basal cells 0.0048 [22]
δ2 The apoptosis rate of virions 0.0021 [22]
δ3 The apoptosis rate of lymphocytes 0.0005 Estimated
D Half-saturation concentration for the progression from P to C 105 [23]
ω1 = ω2 = ω3 Number of virions that are produced by I, P and C 1000 [24]
M The autoimmune rate of I to S 0.00003333 [22]
α1 Contact rate 0.0002 Estimated
ϵ Rate at which the virions are killed by the lymphocytes 0.90005 Estimated
τ1 Rate at which virions are eliminated as they attack normal cells 0.0005 Estimated
S(t) Susceptible basal cells with time
I(t) Infected basal cells with time
P(t) Precancerous basal cells with time
C(t) Cancerous basal cells with time
V(t) Virions with time
L(t) Lymphocytes with time

3
E. Makena et al. Heliyon 10 (2024) e35038

Table 2
Table of sensitivity indices.
Parameter Next Generation Matrix. Survival Function. Constant term of the characteristic polynomial.

F 0.993351 0.0000346813344295621 0.993351

M − 0.99463 − 0.000034681334429562096 − 0.99463
K1 − 0.0000176205 0.026674154887992983 − 0.000158879
r1 1 0.9879380508173314 3
r2 0 5.883416249394072 × 10− 7 1
B − 1.05469 × 10− 10 − 0.000009212399239203786 − 1.05469 × 10− 10

δ1 6.17395 × 10− 6 − 0.22125138109693612 − 0.979586

δ2 0 − 0.12494859329913026 0.00799353
η1 6.43119 × 10− 6 − 0.047374582684402054 6.43119×10− 10
η2 0 − 0.29268294620149876 − 0.510204
η3 0 − 0.1764099056748374 − 0.510204
S − 176.208 − 0.011490780279154997 − 529.625
E 1 0 − 1
L 0 0 0.00135933
K4 0 − 0.0006712604976524636 − 6.36864 × 10− 10

ϵ 0 − 2.188287319016677 × 10− 11 6.36864 × 10− 10

τ1 0 − 0.011499992678394202 − 1.00799
I 0 0.005348560226721884 0
ω3 0 0.001337140056680471 0
ω2 0 0.004813704204049696 0
ω1 0 0.005348560226721884 0
K2 0 0.00029487161100197416 0
D 0 0.04629077249233771 0
P 0 − 0.01364474132726996 0
C 0 − 0.025823926384812002 0
θ 0 − 0.027832326961018047 0
E 0 − 0.00003491252020486956 0
β 0 − 0.08646595437424326 0
V 0 − 0.011241531993419006 0
N 0 − 0.026674154887992983 0

( )
dV V L
= {r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C} 1 − − τ1 SV − ϵ V − δ2 V, (5)
dt K2 L + K4
( )
dL L
= r3 L 1 − − α1 LV − δ3 L. (6)
dt K3

with initial conditions S(0) ≥ 0, I(0) ≥ 0, P(0) ≥ 0, C(0) ≥ 0, V(0) ≥ 0, L(0) ≥ 0.;

3. Results and discussion

This section is organized as follows: 3.1. Model analysis and 3.2 Model simulation.

3.1. Model analysis

This sub-section is organized as follows: 3.1.1 Existence and positive in variance, 3.1.2. Boundedness of the solution, 3.1.3
Disease-free equilibrium point, 3.1.4 Basic reproduction number, 3.1.5 Endemic equilibrium, 3.1.6 Stability analysis, and 3.1.7
Bifurcation analysis.

3.1.1. Existence and positive invariance

Theorem 1. Solutions of the model equations (1)–(6) together with the initial conditions S(0) ≥ 0, I(0) ≥ 0, P(0) ≥ 0, C(0) ≥ 0,
V(0) ≥ 0, L(0) ≥ 0 are always positive or the model variables S (t), I(t), P(t), C(t), V(t) and L(t) all positive for all t and will remain in R+
6
[6].
Proof. .
For the sake of analysis
dS τ2 SV
= r1 S(1 − k1 N) − + Ω1 I − δ1 S, (7)
dt S + BV

4
E. Makena et al. Heliyon 10 (2024) e35038

dI τ2 S
= r1 I(1 − k1 N) + − Ω2 I, (8)
dt S + BV

dP CP2
= r1 P(1 − k1 N) + βI − ϴ 2 − Ω3 P, (9)
dt D + P2
{ }
dC CP2
= r1 C(1 − k1 N) + ϴ 2 − Ω4 C (1 − k4 C), (10)
dt D +P 2

dV L
= {r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C}(1 − k2 V) − τ1 SV − ϵ V − δ2 V, (11)
dt L + K4

dL
= r3 L(1 − k3 L) − α1 LV − δ3 L, (12)
dt
1
Where K1 = k1 ; K12 = k2 ; K13 = k3 ; F+M
EM
= Ω1 ; Ω2 = β + Ω1 + η1 + δ1 ; Ω3 = η2 + δ1 ;
/Ω4 = η3 + δ1 , to reduce the number of parameters, therefore;
(
For t > 0, let W = (s(t), i(t), p(t), n(t), v(t), l(t))T and F(W) = F1(W), F2(W), F3(W), F4(W), F5(W), F6(W))T ,where F1(W) =
{
τ2 SV τ2 SV 2
r1 S(1 − k1 N) − S+BV + Ω1 I − δ1 S, F2(W) = r1 I(1 − k1 N) + S+BV − Ω2 I,F3(W) = r1 P(1 − k1 N) + βI − ϴ DCP
2 +P2 − Ω3 P,F4(W) = r1 C(1 −
}
2
k1 N) + ϴ DCP2 +P2 −
L
Ω4 C (1 − k4 C),F5(W) = {r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C}(1 − k2 V) − τ1 SV − ϵ L+K 4
V − δ2 V,F6(W) = r3 L(1 − k3 L) −
6
α1 LV − δ3 L. Then, system (1)–(6) can be written as dQ 6
dt = F(W) where F : C+ →(R)+ with W(0) = W0 ϵ R+ . Thus, the function W is locally
Lipschitzian and completely stable on R6+ .
Therefore, the solution of the system with non-negative initial conditions exists and is unique. It is also evident that these solutions
exist for all t > 0 and are non-negative, hence the region R6+ is an invariant domain of the system [25].

3.1.2. Boundedness of the solution

For the system to be mathematically meaningful, it is necessary to show that its state variables are positive and bounded for all t.
That is, the solution of the system with a positive initial value will remain positive for all t ≥ 0.
Theorem 2. The positive solutions of the system of model equations (1)–(6) are bounded. That is, the model variables S(t), I(t), P(t), C
(t), V(t) and L(t) are bounded for all t [6].
Proof.
N(t) = S(t) + I(t) + P(t) + C(t) (13)

By differentiating equation (13) gives,

dN dS dI dP dC
= + + + (14)
dt dt dt dt dt
( ) ( ) ( ) ( ) ( )
dN N N N N C
= r1 S 1 − − δ1 S + r 1 I 1 − − η1 I − δ1 I + r1 P 1 − − η2 P − δ1 P + r1 C 1 − − η3 C − δ1 C 1 − (15)
dt K1 K1 K1 K1 K4
( ) ( ( ))
dN N C
= r1 N 1 − − δ1 S + η1 I + δ1 I + η2 P + δ1 P + η3 C + δ1 C 1− (16)
dt K1 K4
( )
dN N
≤ r1 N 1 −
dt K1
By using the separation of variables of inequality, we have
dN
( ) ≤ r1 dt (17)
N 1− N
K1

On integrating both sides (17) and applying the initial conditions to get the value of A and finally substituting the value of A, we
have
K1
N(t) ≤ ( ) . (18)
K1
1+ N(0)
− 1 e− rt

5
E. Makena et al. Heliyon 10 (2024) e35038

Introducing limits, lim N(t) ≤ K1 .Implying that 0 ≤ N(t) ≤ K1 .

t→∞
Also for, V(t).
( )
dV V L
= {r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C} 1 − − τ1 SV − ϵ V − δ2 V. (19)
dt K2 L + K4

Let r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C = r4 V (20)
( ) ( )
dV V L dV V
= r4 V 1 − − τ1 SV − ϵ V − δ2 V, then, ≤ r4 V 1 − (21)
dt K2 L + K4 dt K2
By separating variables (21), integrating and applying initial conditions and limits we obtained lim V(t) ≤ K2 .Therefore; V(t) ≤ K2 .
t→∞
Implying that, 0 ≤ V(t) ≤ K2 .
( )
K3
Similarly, L(t), dL
dt ≤ r3 L 1 −
L
K3 . Therefore; L(t) ≤ ( ) . On applying the limits as t→∞, it follows that; 0 ≤ L(t) ≤ K3 .
K3
1+ L(0)
− 1 e− rt

3.1.3. Disease-free equilibrium point (DFE)

( ) ( )
Theorem 3. The system of equations (1)–(6) has disease-free equilibrium point E0 obtained as; E0 = S0 , I0 , P0 , C0 , V 0 , L0 =
( )
K1 (r1 − δ1 )
r1 , 0, 0, 0, 0, Kr33 (r3 − δ3 ) .

Proof. .
DFE of system (1)–(6) is obtained by setting the right-hand side to zero and equating the infectious classes to zero I = 0, P = 0, C =
0, A = 0 and V = 0 [26].
Solving, we obtained. S0 = 0 and S0 = K1 (rr11− δ1 )
L0 = 0 and L0 = K3 (r3 − δ3 )
r3 . The set, S0 = L0 = 0 was not biologically meaningful
because it is not feasible to have a cervix with no basal cells and Lymphocytes. Hence S0 = K1 (rr11− δ1 )
and L0 = K3 (rr33− δ3 )
. Hence the disease
free equilibrium point of the system of equations (1)–(8) was obtained as:
( )
( ) K1 (r1 − δ1 ) K3
E0 = S0 , I0 , P0 , C0 , V 0 , L0 = , 0, 0, 0, 0, (r3 − δ3 ) .
r1 r3

3.1.4. Basic reproduction number (R0 )

There are numerous controversies surrounding the method of calculating basic reproduction number (R0 ) because it has been
proven that each method produces a unique estimate of (R0 ) hence posing a challenge to stakeholders on how best to control the
dynamics of a disease [12]. Although most studies have evaluated R0 using the next-generation method, it estimates R0 as an average
regardless of whether the population is of human or host cells. It also lacks some uniqueness [12]. To address the challenges of those
methods, this study determined R0 using three methods for comparative purpose: Sub-Sub section 3.1.4.1 next-generation matrix,
3.1.4.2 Survival Function and 3.1.4.3 constant term of characteristic polynomial.

(M+F)(K1 − S0 )γ 1
3.1.4.1. Using next generation matrix. Theorem 4. The basic reproduction number , R0 = K1 (Mβ+Fβ+ME+Mδ1 +Fδ1 +M η by Next-
1 +F η1 )
generation method.
Proof. .
To obtain basic reproduction number, we use next generation matrix method proposed by Kermack and Diekmann [27]. Let
non-negative matrix f and non-singular matrix v represent new infections terms and transfer of infections terms respectively. Thus
⎡ ( ) ⎤ ⎡ ⎤
N τ2 SV EM
r1 I 1 − + βI + I + η I + δ I
⎢
⎢ K1 S + BV ⎥
⎥
⎢
⎢ F+M 2 1 ⎥
⎥
⎢ ( ) ⎥ ⎢ ⎥
⎢ N ⎥ ⎢ θCP 2 ⎥
⎢ r 1 P 1 − ⎥ ⎢ − βI + + η P+δ 1 P ⎥
⎢
⎢ K 1
⎥
⎥
⎢
⎢ D2 + P2 2 ⎥
⎥
f =⎢ ( )( ) ⎥; v = ⎢ ( )( ) ⎥
⎢ N C ⎥ ⎢ θCP 2
C ⎥
⎢ r 1 C 1 − 1 − ⎥ ⎢ − + η 3 C + δ 1 C 1 − ⎥
⎢
⎢ K1 K4 ⎥
⎥
⎢
⎢ D2 + P2 K4 ⎥⎥
⎢ ( )⎥ ⎢ ⎥
⎣ V ⎦ ⎣ ϵLV ⎦
r1 V + ω1 η1 I + ω2 η2 P + ω3 η3 C 1− τ1 SV + + δ2 V
K2 L + K 4

At E0 point, Jacobian matrices of f and v was evaluated to find out matrices F and V respectively,

6
E. Makena et al. Heliyon 10 (2024) e35038

⎡( ) ⎤
S0 − βS0 V τ2 ⎡ ⎤
1 − γ 0 0 EM
⎢ K1 1 S2 ⎥ β + + δ 1 + η 1 0 0 0
⎢
⎢ ( )
⎥
⎥
⎢
⎢ F+M ⎥
⎥
⎢ S0 ⎥ ⎢
δ1 + η2 0 0
⎥
⎢
⎢ 0 1− γ 0 0 ⎥
⎥
⎢
⎢
− β ⎥
⎥
F=⎢ K1 1 ⎥;V = ⎢ ⎥
⎢
⎢ ( ) ⎥
⎥
⎢
⎢ 0 0 δ1 + η3 0 ⎥
⎥
⎢ S0 ⎥ ⎢ ⎥
⎢ 0 0 1− γ1 0 ⎥ ⎣ L0
ϵ ⎦
⎣ K1 ⎦ 0 0 0 + δ 2 + τ 1 S0
V
L0 + K 5
η1 ω 1 η2 ω2 η3 ω 3 γ2
⎡ ( ) ⎤
S 0
1− r1
⎢ K1 ⎥
⎢
⎢ MQ 0 0 0 ⎥
⎥
⎢ + β + δ1 + η1 ⎥
⎢ ⎥
⎢F + M ⎥
⎢ ( ) ⎥
⎢
⎢ S 0 ⎥
⎥
⎢ 1− r1 ⎥
⎢ K1 ⎥
⎢ 0 0 0 ⎥
At E point, FV , was obtained as, ⎢
0 − 1 ⎢ δ1 + η2 ⎥;
⎥
⎢ ( ) ⎥
⎢ S0 ⎥
⎢ 1− r1 ⎥
⎢ K 1
⎥
⎢
⎢ 0 0 0 ⎥
⎥
⎢ δ1 + η3 ⎥
⎢ ⎥
⎢ r2 ⎥
⎢ 0 0 0 ⎥
⎣ L0 ϵ ⎦
+ δ2 + τ 1 S V
0
L0 + K5
(K1 − S0 )γ (K1 − S0 )γ (M+F)(K1 − S0 )γ1 (L0 +K4 )γ2
The four eigenvalues of FV − 1 at E0 were: K1 (δ1 +η )1 , K1 (δ1 +η )1 , K1 (Mβ+Fβ+ME+Mδ1 +Fδ1 +M η +Fη ) and ϵL0 +δ2 L0 +K4 δ2 .
3 2 1 1

By inspection method (which was be verified by numerical method too) the dominant eigenvalue represents the R0 [27]. Hence
( )
(M + F) K1 − S0 γ 1
R0 =
K1 (Mβ + Fβ + ME + Mδ1 + Fδ1 + Mη1 + Fη1 )

3.1.4.2. Basic reproduction number by survival function. Theorem 5. The basic reproduction number
( ) ( ) ( ) ( )
τ2 SV
r1 1 − KN1 + S+BV r1 1 − KN1 r1 1 − KN1 (r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C) 1 − V
K2
R0 = + +( ( )) +
EM
β + F+M + η1 + δ1 D2θCP + η2 + δ2 ϵL
+ δ2 + τ1 SV
+P2 (η3 + δ1 ) 1 − KC4 L+K4

by the method proposed by Ref. [28]. Where N(0), S(0), V(0), P(0), C(0), I(0) and L(0) are assumed to be constant terms at the
beginning of the epidemic. For numerical computation, in this study, we assumed those constant values to be the initial conditions of
the model.
Proof. .
Evaluation of R0 using the survival function method is considered to be a more accurate method. One benefit of the survival
function approach is that it consistently yields the average number of secondary basal cells infected by one infected basal cell in the
same class [12]. It also takes into account the different attributes of the population, therefore, the basic reproduction number using the
∫∞
survival function is given by; R0 = 0 (k ×b ×p)dt where:k = rate at which an individual in that class causes an infection.,b =
probability at which an infected individual remains in the same class to cause an infection. p = probability that an infected case will
enter that class [28].
Considering the infectious classes (I, P, C, and V), and N, S, V, P, C, I and L are assumed to be constant terms at the beginning of the
epidemic, evaluating R0 of equations (2)–(5), the basic reproduction number was obtained by summing them [12] that is,
( ) ( )
{ ( ) } ∫ ∞ − β+ EM + η + δ T { ( )} ∫ ∞ − θCP
+ η2 + δ2 TA
N τ2 SV F+M 1 1 A N D2 +P2
r1 1 − + e dTA + r1 1 − e dTA
K1 S + BV 0 K1 0
( ( )) ( )
{ (
N
)} ∫ ∞ − (η +δ1 ) 1− C
3 K4
TA
{ (
V
)} ∫ ∞ − ϵL
L+K4
+ δ2 + τ 1 SV TA
+ r1 1 − e dTA + (r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C) 1 − e dTA
K1 0 K2 0

⎡ ( ) ⎤∞ ⎡ ( ) ⎤∞
{ ( ) } − β+
EM
F+M
+ η1 + δ1 TA { ( )} ⎢ −
θCP
D2 +P2
+ η2 + δ2 TA ⎥ { ( )}
⎢
τ2 SV ⎢
⎥ ⎢ ⎥
On integration, r1 1 − N
K1
+S+BV ⎣−
e
EM
β+F+M+η1 +δ1
⎥ + r1 1 −
⎦
N
K1 ⎢− e
θCP +η +δ ⎥ + r1 1 − N
K1
⎣ 2 2 ⎦
D2 +P2

0 0

7
E. Makena et al. Heliyon 10 (2024) e35038

⎡ ( ( )) ⎤∞ ⎡ ( ) ⎤∞
⎢ − (η3 +δ1 ) 1−
C
K4 ⎥ TA { ( ϵL
)}⎢ − L+K 4
+ δ 2 + τ 1 SV TA ⎥
⎢ )) ⎥ ⎢e ⎥
⎢− e ( ( ⎥ + (r2 V +ω1 η1 I +ω2 η2 P +ω3 η3 C) 1 − V
K2 ⎢ ϵL +δ +τ SV
2 1
⎥ . Putting the limits the basic repro­
⎣ C
⎦ ⎣ L+K4 ⎦
(η3 +δ1 ) 1− K4
0 0
duction number is obtained as:
( ) ( ) ( ) ( )
τ2 SV
r1 1 − KN1 + S+BV r1 1 − KN1 r1 1 − KN1 (r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C) 1 − V
K2
R0 = + +( ( )) +
EM
β + F+M + η1 + δ1 D2θCP + η2 + δ2 ϵL
+ δ2 + τ1 SV
+P2 (η3 + δ1 ) 1 − KC4 L+K4

3.1.4.3. Basic reproduction number by evaluating the constant term of the characteristic polynomial. Theorem 6. The basic reproduction
number
3
(F+M)(S0 − K1 ) (L0 +K4 )r3 r2
R0 = K3 (Fβ+M(Q+β)+(F+M)(δ +η ))(δ +η )(δ +η ) 1L0 ϵ+ L0 +K (δ +Sτ ) by method evaluating the constant term of the characteristic polynomial
1 1 1 1 2 1 3 ( ( 4) 2 1 )

[12].
Proof. .
For comparison purposes, there is a need to determine the basic reproduction number using the constant term of a characteristic
polynomial. When λ max = 0, the constant term of the characteristic polynomial will be zero [12]. However, the reverse is not true, as
the polynomial could have both zero and positive roots. The characteristic polynomial by Next generation matrix (FV − 1 ) is of the form
b4 λ4 + b3 λ3 + b2 λ2 + b1 λ + b0 = 0 , where the expressions of b4 , b3 , b2 , b1 and b0 are found in the appendix in section 4. The study
[12] proposes the following conditions.
Let b0 represent R0 , then b0 = 0 is a threshold if bj ≥ 0 for all j. Some sufficient conditions include; non-constant coefficients all
being positive and bj ≥ 0 under the constraint b0 = 0 (so that the largest eigenvalue at b0 = 0 is 0). This method is significantly easier to
use than finding the largest eigenvalue, although verifying that bj = 0 necessarily corresponds to the largest eigenvalue can become
complicated for some models [12].
( )3 ( )
(F + M) S0 − K1 L0 + K4 r13 r2
R0 = ( ( ) ) .
K31 (Fβ + M(Q + β) + (F + M)(δ1 + η1 ))(δ1 + η2 )(δ1 + η3 ) L0 ϵ + L0 + K4 (δ2 + Sτ1 )

3.1.5. Endemic equilibrium

By setting the system of equations to zero and evaluating the state variables, the endemic equilibrium points would be in the form:

(EEP) = (S* , I* , P* , C* , V * , L* )
{ ( ) }( )
2
C* P * C*
From; r1 C* 1 − N
K1 +ϴ 2 − η3 C* − δ1 C* 1− K4 = 0;
D2 +P*
*
It follows that; C = K4 ;
The endemic equilibrium point exists at C* = K4 and since the equations are highly non-linear, it was not tractable to solve them
explicitly.
( )
* *2
From; r1 P* 1 − KN1 + βI * − ϴ C2 P * 2 − η2 P* − δ1 P* = 0;
D +P
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
− CθK1 ± C2 θ2 K21 − 4(βK1 − Nr1 +K1 r1 − K1 δ1 − K1 η2 )(D2 βK1 − D2 Nr1 +D2 K1 r1 − D2 K1 δ1 − D2 K1 η2 )
It implies that; P* = 2(− βK +Nr − K r +K δ +K η ) ;
( 1 )1 1 1 1 1 1 2
* *
From; {r2 V * + ω1 η1 I* + ω2 η2 P* + ω3 η3 C* } 1 − KV2 − τ1 SV − ϵ L* L+K4 V * − δ2 V * = 0
⎛ √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ⎞
( )2
K2 ⎝KLϵ4 − r2 + JηK12ω1 + PηK22ω2 + TηK32ω3 + 4r2 (Jη1 ω1 +Pη2 ω2 +Tη3 ω3 )
K2
+ − K4 + r2 − K2 − K2 − K2
Lϵ Jη1 ω1 Pη2 ω2 Tη3 ω3 ⎠

V* =
2r2
Substituting the value of P , C and V to; * * *
( )
τ2 S0 V M
r1 I* 1 − KN1 + S+BV * EM * * * * *
* − βI − F+MI − η1 I − δ1 I > 0 and solving for I , we get the value of I = −
( ( τ2( ) );
EM X
S+BV) − F+M− β+ 1− K1 r1 − δ1 − η1
( )
L*
From r3 L* 1 − K3 − α 1 L* V * − δ 3 ;

K3 (r3 − Vδ1 − δ3 )
L* =
r3

8
E. Makena et al. Heliyon 10 (2024) e35038

Theorem 7. The necessary and sufficient conditions for existence is R0 > 1, dI dP dC

dt > 0, dt > 0, dt > 0 and
dV
dt > 0 [29].

3.1.6. Stability analysis

3.1.6.1. Local stability of disease-free equilibrium point. Theorem 2. The Disease Free Equilibrium of the system (1)-(6) is locally
asymptomatic stable whenever R0 < 1.
Proof. .
The relative stability of the system can be determined by the Routh-Hurwitz criterion of stability without having to solve each
equation.
Determining Jacobian Matrix of system (1)–(6) at Disease Free Equilibrium is obtained;

[ R1 R2 R3 R4 R5 R6 ]T
[ ]
( ) ( ) ( ) ( ) τ2 B
R1 = − δ1 − k1 r1 S0 1 − k1 S0 Ω1 − k1 r1 S0 1 − k1 S0 − k1 r1 S0 1 − k1 S0 − k1 r1 S0 1 − k1 S0 0
S0
[ ]
( ) τ2 B
R2 = 0 1 − k1 S0 r1 − Ω2 0 0 − 0
S0
[ ( ) ]
R3 = 0 β 1 − k1 S0 r1 − Ω3 0 0 0
[ ( ) ]
R4 = 0 0 0 1 − k1 S0 r1 − Ω4 0 0
[ ]
ϵL0
R5 = 0 η1 ω1 η2 ω2 η3 ω3 − + r2 − δ2 − τ1 S0 0
L0 + K4
[ ( )]
R6 = 0 0 0 0 − α 1 L0 − δ3 − k3 r3 L0 1 − k3 L0

The characteristic polynomial is obtained as a0 λ6 + a1 λ5 + a2 λ4 + a3 λ3 + a4 λ2 + a5 λ + a6 = 0, where the expression ai , i = 1 for

system (1)–(6) are in appendix section 2. The Routh table for the coefficients was also derived in the appendix in section 1.
From the characteristic polynomial, the values a2 ,a4 ,a6 ,b1 ,c1 , d1 and e1 are determined using Mathematica Software and expressed
in terms of Ro. By Routh-Hurwitz Criteria for stability, system (1)–(6) is locally asymptotically stable at DFE whenever Ro < 1 if and
only if a2 > 0, a4 > 0, a6 > 0, b1 > 0, c1 > 0, d1 > 0, e1 > 0 are satisfied and otherwise unstable [26]. The values of a2 , a4 , a6 , b1 , c1 ,
d1 and e1 have been expressed in the appendix.

3.1.6.2. Global stability of disease-free equilibrium point. The global stability of disease-free equilibrium is investigated using the
Castillo-Chavez Metzler Matrix method. dX dZ
dt = F(X, Z); dt = G(X, Z), G(X, 0) = 0 Where; X = (S, L) εR2 denote non-infectious cervical
+

cancer compartments and Z = (I, P, C, V) ε R4 denote the infectious cervical cancer compartments Eo = (X*, 0) represents the disease-
+

free equilibrium of the system if this point satisfies following conditions.

dX
i. dt = (X, 0), Where X* is globally asymptotically stable.
dZ
ii. dt= DzG(X, 0)Z − G(X, Z) ≥ 0 for all X,Z ∈ Ω, then we can conclude that Eo is locally asymptotically stable if the following the­
orems hold.

Theorem; The equilibrium point Eo = (X*, 0) of the system [1–6] is globally asymptotically stable if R*0 ≤ 1 and the conditions (i)
and (ii) are satisfied, otherwise unstable. From equation (1) two vectors function G (X, Z) and F (X,Z), we consider systems dX
dt = (X, 0) =
0 Letting A = Dz G(X*, 0), which is the Jacobian of Ĝ(X,Z) taken in (I, P, C, V) and evaluated at (X*,0) such that the matrix A is given by;
⎡( ) ⎤ ⎡ ⎤
β1 S0 − k1 E + η1 β1 S0 I + η2 β1 S0 T + η3 β1 S0 C (λ1 + λ2 )S − k1 E
⎢ ρE − k 2 I ⎥ ⎢ ρE − k2 I ⎥
AZ = ⎢ ⎥ G(X, Z) = ⎢ ⎥
⎣ ωI − k3 T ⎦ ⎣ ωI − k3 T ⎦
γI + αT − k4 C γI + αT − k4 C
⎤ ⎡
G
⎢ 1
ˆ
(X, Y) ⎥
⎡( ( ( )
β1 E + η1 I + η2 T + η3 C) S0 − S − λ2 S
⎤
⎢ Gˆ (X, Y) ⎥ ⎢ 0 ⎥
But Gˆ(X, Z) = AZ − G(X, Z), This reduces to Gˆ(X, Z) = ⎢ ˆ2
⎢ ⎥ ⎥;
⎥ =⎢
⎢ G3 (X, Y) ⎥ ⎣ 0 ⎦
⎣ ⎦
G4 (X, Y)
ˆ 0

Thus, if then the disease-free equilibrium (E0 ) is globally stable and unstable otherwise. The susceptible is bounded as, S ≤ S0 .Thus,

9
E. Makena et al. Heliyon 10 (2024) e35038

( ( )
DFE, E0 is globally asymptomatically stable if and only if β1 (E +η1 I +η2 T +η3 C) S0 − S ≥ λ2 S [26,29,30].

3.1.6.3. Global stability analysis of endemic equilibrium point. Lyapunov functions are mathematical tools used to study the stability of
dynamical systems. Different Lyapunov functions are employed for global stability analysis such as; quadratic, non-quadratic, radial
basis functions, composite and piecewise Lyapunov functions. This study adopted composite Lyapunov function.
∑ ( )
L= bi xi − x*i ln xi

The function L is said to be positive definite if it satisfies the following conditions.

i) Strict positivity: L > 0 for all x ∕

= 0.
ii) Zero at origin: L(0) = 0.

Where bi the constant is selected such that bi > 0, xi is the population of the ith compartment and x*i is the endemic equilibrium
point.

L = b1 (S − S* ln S) + b2 (I − I* ln I) + b3 (P − P* ln P) + b4 (C − C* ln C) + b5 (V − V * ln V) + b3 (L − L* ln L)
( ) ( ) ( ) ( ) ( ) ( )
dL S* dS I* dI C* dC P* dP V * dV L* dL
= b1 1 − + b2 1 − + b3 1 − + b4 1 − + b5 1 − + b6 1 −
dt S dt I dt C dt P dt V dt L dt
( ){ } ( ){
dL S* τ2 S I* τ2 S
= b1 1 − r1 S(1 − k1 N) − + Ω1 I − δ1 S + b2 1 − r1 I(1 − k1 N) +
dt S S + BV I S + BV
} ( ){{ } } ( ){
C* CP2 P* CP2
− Ω2 I + b3 1 − r1 C(1 − k1 N) + ϴ 2 − Ω4 C (1 − k4 C) + b4 1 − r1 P(1 − k1 N) + βI − ϴ 2
C D +P 2 P D + P2
} ( ){ } ( )
V *
L L *
− Ω3 P + b5 1 − {r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C}(1 − k2 V) − τ1 SV − ϵ V − δ2 V + b6 1 − {r3 L(1 − k3 L)
V L + K4 L
− α1 LV − δ3 L}

After solving we obtain the value of X and Y as follows.

bS* τ2 θCP2
X = b1 r1 S + b1 Ω1 I + bS* r1 K1 N + + bS* δ1 + b2 r1 I + b2 τ2 S + b2 I* K1 N + b2 I* Ω2 + b3 r1 C + b3 2
S + BV D + P2
θC* CP2 K4
+ b3 r1 C2 K1 K4 N + b3 Ω4 C2 K4 + b3 K1 NC* + b3 Ω4 C* + b3 K4 CC* r1 + b3 2 + b4 r1 P + b4 β I + b4 r1 PK1 NK4 C
D + P2
( (
θCP 2
P *
θCP 2
+ b4 2 K4 C + b4 Ω4 C2 K4 + b4 r1 PK1 N + 2 + Ω4 C + r1 PK4 C + βICK4 + b5 r2 V + ω1 η1 I + ω2 η2 P
D +P 2 P D +P 2
( )
V* ϵLV L* ( )
+ ω3 η3 C + ω2 η2 PK2 V + b5 r2 K2 V + ω1 η1 IK2 V + ω3 η3 CK2 V + τ1 SV +
2
+ δ2 V + b6 r3 L + b6 r3 L2 K3 + α1 LV + δ3 L
V L + K4 L

bSτ2 bS* Ω1 I I* Sτ2

Y = − br1 SK1 N − − bSδ1 − bS* r1 − − b2 r1 IK1 N − b2 Ω2 I − b2 I* r1 − b2 − b3 r1 CK1 N − b3 Ω4 C − b3 r1 C2 K4
S + BV S I S + BV
( )
θC2 P2 K4 θC* P2 θCP2
− b3 2 − b3 C* − b3 2 − b3 C* Ω4 CK4 − b3 C* K1 CK4 r1 − b4 r1 PK1 N + 2 + Ω4 C + r1 PK4 C + β I K4 C
D +P 2 D +P 2 D +P 2
( ) ( )
P* θCP2 ϵLV
− b4 r1 P + β I + r1 PK1 NK4 C + 2 K4 C − b5 r2 K2 V 2 + ω1 η1 IK2 V + ω3 η3 CK2 V + τ1 SV + + δ2 V
P D +P 2 L + K4
V* ( )
2
− b5 (r2 V + ω1 η1 I + ω2 η2 P + ω3 η3 C + ω2 η2 PK2 V) − b6 r3 L K3 + α1 LV + δ3 L − b6 L r3 *
V

By inspection method X > Y, therefore this result shows that cervical cancer would persist whenever X > Y irrespective of the
initial conditions, and if Y > X the disease will die out irrespective of initial conditions. The global stability for EEP exists for people
without underlying health conditions, hence implying that the global stability for EEP for system (1–6) [29].

3.1.7. Bifurcation analysis

Most mathematical models often undergo bifurcation which makes the control of most diseases difficult. Utilizing the center
manifold theory, the likelihood of population hopf bifurcation was investigated. The renaming of variables is done simply by letting,
S = x1 , I = x2, P = x3 , C = x4, V = x5, L = x6 ;
Using vector notation;
x = x1 , x2, x3, x4, x5, x6

System (1)–(6) is written as,

10
E. Makena et al. Heliyon 10 (2024) e35038

( )T
dx
dt = F(x) where F = f1 , f2, f3 , f4 , f5 , f6 , It follows that:
( )
dx1 N τ2 f1 f5 EM
= p1 = r1 f1 1 − − + f2 − δ1 f1
dt K1 f1 + Bf5 F + M
( )
dx2 N τ2 f1 f5 EM
= p2 = r1 f2 1 − + − βf2 − f2 − η1 f2 − δ1 f2
dt K1 f1 + Bf5 F+M
( )
dx3 N f4 f3 2
= p3 = r1 f3 1 − + βf2 − ϴ − η2 f3 − δ1 f3 .
dt K1 D + f3 2
2

{ ( ) }( )
dx4 N f4 f3 2 f4
= p4 = r1 f4 1 − +ϴ − η3 f4 − δ1 f4 1− .
dt K1 D2 + f3 2
K4

( )
dx5 { } f5 f6
= p5 = r2 f5 + ω1 η1 f2 + ω2 η2 f3 + ω3 η3 f4 1 − − τ1 f1 f5 − ϵ f5 − δ2 f5 .
dt K2 f6 + K4
( )
dx6 f6
= p6 = r3 f6 1 − − α1 f6 f5 − δ3 f6 .
dt K3
Jacobian solved at DFE,
( )
E0 = (S0, I0, P0, C0, V0, L0) = K1 (r1 − δ1 )
r1 ,0,0,0,0, Kr33 (r3 − δ3 ) , in a case, where R*0 = 1 and further, suppose that r1 = r1* is a bifurcation

parameter, then solving for r1* from R*0 = 1 we get,

K1 (MQ + Fβ + Mβ + Fδ1 + Mδ1 + Fη1 + Mη1 )

r1* = ( )
(F + M) K1 − S0

Gives
[ ]
A1 A2
A3 A4

Where;
⎡( ) ⎤
S0 * r1* S0 ME r* S0 r1* S0
1 − r − − δ1 − 1 −
⎢
⎢ K1 1 K1 F+M K1 K1 ⎥
⎥
⎢ ( ) ⎥
⎢ ⎥
⎢ ME S0 * ⎥
A1 = ⎢ 0 − − β+ 1− r − δ1 − η1 0 ⎥
⎢
⎢ F+M K1 1 ⎥
⎥
⎢ ( ) ⎥
⎣ S0 * ⎦
0 β 1− r − δ1 − η2
K1 1
It can easily be shown that, Jacobian of the system
⎡ ⎤
r* S0 τ2 B
⎢− 1 0 ⎥ ⎡ ⎤
⎢ K1 S0 ⎥ 0 0 0
⎢ ⎥
A2 = ⎢
⎢ 0 τ2 B ⎥ A3 = ⎣ 0 η1 ω1 η2 ω2 ⎦
⎢ − 0 0⎥ ⎥ 0 0 0
⎣ S ⎦
0 0 0
⎡( ) ⎤
S0 *
⎢ 1− r − δ1 − η3 0 0 ⎥
⎢ K1 1 ⎥
⎢ ⎥
⎢ ⎥
⎢ ϵL0 ⎥
A4 = ⎢ η3 ω3 + r2 − δ2 − τ1 S0 0 ⎥
⎢
⎢ L0 + K4 ⎥
⎥
⎢ ( ) ⎥
⎣ L0 r 3 L0 ⎦
0 − α 1 L0 1− r3 − − δ3
K3 K3

It has been proved that at least one of the eigenvalues of the matrix is a simple zero eigenvalue [26]. Therefore, the bifurcation of
the system can be evaluated using the Castillo-Chavez theorem.
( )T
Let u = u1 , u2, u3 , u4 , u5 , u6 be the right eigenvector and v = (v1 , v3 , v4 , v5 , v6 )T be the left eigenvector linked with zero eigen­
values of the Jacobian matrix near r1 = r1* , of the system.

11
E. Makena et al. Heliyon 10 (2024) e35038

Solving the system of the equations, we obtain;

u3 β − u5 (η1 ω1 ) − u5 (η3 ω3 ) − u5 (η3 ω3 )
u1 = 0, u2 = ( ) , u3 = ( ) , u4 = ( )
Ω1 + β + S0 k1 − 1 − δ1 − η1 1 − S0 k1 r1 − η3 1 − S0 k1 r1 − η4

τ2 B τ2 B
u2 − + u6 α1 L0 u2 −
u5 = S0
= S0
, u6 = 0
ϵL0
L0 +k4
+ r2 − δ2 − τ1 S0 ϵL0
L0 +k4
+ r2 − δ2 − τ1 S0

( )
− η1 − r1 S0 k1 v2 + v3 r1 S0 k1 − r1 S0 k1 v4 + v5 τS20B v3 r1 S0 k1 + v5 τS20B
v1 = ( ) =( ) , v2 = 0
1 − S k1 r1 − r1 S k1 − δ1
0 0 1 − S0 k1 r1 − r1 S0 k1 − δ1
(( ) )
S0 k1 − 1 r1 + η1+δ1 + Ω1 + β v2 + v5 τS20B v5 τS20B − (η ω2 )v3 − (η3 ω3 )v4
v3 = ( ) =( ) , v4 = 0, v5 = εL0 2
1 − S0 k1 r1 − η3 1 − S0 k1 r1 − η3 K4 +L0
+ r2 − δ2 − τ1 S0
− (η2 ω2 )v3 v5 α 1 L
= , v6 = ( ) .
εL0
+ r2 − δ2 − τ1 S0 ( ) r3 L0
K4 +L0 1− L0 k1 r3 − k3
− δ3

Let pk be the kth component of p and

n
∑ ∂2 pk ∑ n
∂2 pk
a= vk ui uj (0, 0) and b = vk ui (0, 0)
k,ij=1
∂fi ∂fj k,ij=1
∂fi ∂r1*

then the local dynamics of the system around the equilibrium point (0,0) is totally determined by the signs of a and b [31].
On evaluating the values of a and b which are found in the appendix in section 3, we conclude that since, a < 0 and b < 0, when r1* <
⃒ ⃒
0 with ⃒r1* ⃒ ≪ 1, (0, 0) is unstable; when 0 < r1* ≪ 1, (0, 0) is then asymptotically stable and there exists a positive unstable equilibrium.

3.2. Numerical simulation

MATLAB2019a is utilized in numerical simulation to illustrate the non-linear ODE’s dynamic behavior in system (1)–(6). The
simulations are run with the initial conditions and parameter values (taken from the literature review and graphically depicted) in
Table 1.
The graphs obtained from simulations are interpreted as follows.
The immunity level for the infected cells varied from M = 0.0000333 to M = 0.000000000000001333 and M = 0.9333 while the
other parameters were constant. From Fig. 2 it shows that increasing the immunity level it reduces the number of infected cells. This
indicates that eating a balanced diet can boost the immunity level which fights the virus reducing the number of infected cells.
The infection rate for the infected cells was varied from γ 2= 0.000001 to γ2= 0.5000001 while the other parameters were held
constant. From Fig. 3 it indicates that when the rate of infection was increased the number of infected cells increased.
The Progression rate from Precancerous cells to Cancerous cells varied from θ = 0.000001 to θ = 0.101 while the other parameters
were held constant. Fig. 4 shows that when the progression rate was decreased it increased the pre-cancerous cells while when it was
increased it significantly reduced the number of pre-cancerous cells.

Fig. 2. A graph of variation of Immunity level against time.

12
E. Makena et al. Heliyon 10 (2024) e35038

Fig. 3. A graph of variation of the rate of infection against time

Fig. 4. A graph of variation of progression rate from pre-cancerous to cancerous against time.

Fig. 5. A graph of variation of progression rate from pre-cancerous to cancerous against time

13
E. Makena et al. Heliyon 10 (2024) e35038

Fig. 6. A graph of variation of ∈ against time for Virions.

Fig. 7. A graph of variation of Ω1 against time for Virions

Fig. 8. A graph of Variation of η3 against time for Virons

14
E. Makena et al. Heliyon 10 (2024) e35038

Fig. 9. A graph of Variation of η2 against time for Virons.

Fig. 10. A graph of Variation of τ1 against time

Fig. 11. A graph of variation of η1 against time for Virons

15
E. Makena et al. Heliyon 10 (2024) e35038

Fig. 12. A graph of Variation of θ for Virons against Lympocytes.

Fig. 13. A graph of Susceptible cells against Virons with variation of contact rate (τ1)

Fig. 14. A plot of variation of θ for a graph of pre-cancerous cells against cancerous cells.

The progression rate from Precancerous cells to Cancerous cells was varied for Cancerous cells while the other parameters were
held constant. The results as represented by Fig. 5indicate that when the rate was reduced it decreased the number of cancerous cells
while when it was increased, the population of cancerous cells increased (see Fig. 6).

3.3. Sensitivity analysis of the model

The relationship below provides the sensitivity index of the model parameter;

16
E. Makena et al. Heliyon 10 (2024) e35038

Fig. 15. A plot of variation of α1 for a graph of Virons against Lympocytes

∂Ro X
SRXo = *
∂X Ro
Sensitivity analysis for the whole model is as follows.
From the sensitivity indices table above, r1 is the most positive parameter implying that it is directly related to the dynamics of
high-risk HPV to cervical cancer. To reduce the risk of cervical cancer r1 should be decreased while, S being the most negative means
that it is inversely related to the dynamics of HPV to cervical cancer. Increasing the value of S reduces the risk of the disease.

4. Discussion

Based on the analysis done, we discovered that the model was bounded and fell within the positive region; the DFE persisted even in
the absence of the disease; the endemic equilibrium point was present when R0 > 1; and the local stability of the DFE is unstable when
R*0 > 1 and locally asymptotically stable when R*0 < 1. When R0 > 1, the global stability of EEP is asymptotically stable. Additional
examination was conducted utilizing the simplified system of the model (1)–(6), demonstrating that the Disease Free Equilibrium’s
global stability is asymptotically stable if R*0 < 1 and the lack of backward bifurcation indicates that is feasible to completely eradicate
cervical cancer. The reproduction number was found to have a numerical value of − 7.2485210 × 10− 6 , which was used to simulate
the model. It was observed that increasing parameter γ2 increased the number of infected cells hence increasing the risk to cervical
cancer while decreasing M reduced the number of infected cells hence the need to focus on a balanced diet so as to boost the immunity
levels. It was shown that θ has the most direct impact and the model was able to estimate that a patient can die within 10 days when
θ = 0.01, however, the patient can live up to 20 days when θ = 0.09.

5. Conclusion

In the presence of immunity and functional responses, this work aimed to construct an in-host density-dependent deterministic
model for the dynamics of basal cells, virions, and lymphocytes and their consequences for cervical cancer. The general (SIVPC) model
by Ref. [1] was adjusted for our investigation to include the maximum number of malignant cells that can result in a patient’s death. A
system of six first-order non-linear ordinary differential equations that explain the dynamics of HPV to cervical cancer helped to
achieve this goal in section 2. In this work, the survival function, method of characteristic polynomial, next-generation matrix, pos­
itivity and boundedness of the solution, equilibrium points, and basic reproduction number were used to examine the behavior of the
deterministic model. The Routh-Hurwitz criteria for stability were used to determine the local stability of the DFE, while the
Castillo-Chavez approach was used to determine the global stability of the EEP. Bifurcation analyses were also produced. Additionally,
sensitivity indices for the model’s system were calculated using the next-generation matrix, survival function, and characteristic
polynomial, depending on the reproduction number. Also, sensitivity indices for the system of the model was performed based on the
reproduction number using three methods which are next-generation matrix, survival function, and characteristic polynomial.
The model suggests a proportion of 75% of cancerous cells that can lead to the death of a cervical cancer patient however, future
studies should focus to obtaining the real data for the proportion of cancerous cells that can lead to the death of a patient.

Funding statement

We recognize the University of Embu for granting me a full scholarship to pursue my master’s degree.

17
E. Makena et al. Heliyon 10 (2024) e35038

Data availability statement

Data used was a secondary data which was obtained from literature review.

CRediT authorship contribution statement

Elosy Makena: Writing – review & editing, Writing – original draft, Methodology. Cyrus Gitonga Ngari: Supervision. Patrick
Mwangi Kimani: Supervision. Jeremiah Savali Kilonzi: Software.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Acknowledgement

The authors acknowledge the University of Embu for their kind support.

Appendix

Section 1.

The characteristic polynomial of the above matrix is given as;

a0 λ6 + a1 λ5 + a2 λ4 + a3 λ3 + a4 λ2 + a5 λ + a6 = 0,

Where constant a0 , a1, a2, a3, a4, a5 , a6 are determined using Mathematica software as;

Routh Table
Label

λ6 1 a2 a4 a6
λ5 a1 a3 a5 0
λ4 b1 b2 b3 0
λ3 c1 c2 0 0
λ2 d1 d2 0 0
λ1 e1 0 0 0
λ0 f1 0 0 0

Where,
⃒ ⃒ ⃒ ⃒ ⃒ ⃒
⃒1 a2 ⃒⃒ ⃒1 a4 ⃒⃒ ⃒1 a6 ⃒⃒
− ⃒⃒ − ⃒⃒ − ⃒⃒
a1 a3 ⃒ a1 a2 − a3 a1 a5 ⃒ a1 a4 − a5 a1 0⃒
b1 = = b2 = = b3 = = a6
a1 a1 a1 a1 a1
⃒ ⃒ ⃒ ⃒ ⃒ ⃒
⃒a a3 ⃒⃒ ⃒a a5 ⃒⃒ ⃒b b2 ⃒⃒
− ⃒⃒ 1 − ⃒⃒ 1 − ⃒⃒ 1
b1 b2 ⃒ b1 a3 − a1 b2 b1 b3 ⃒ b1 a5 − a1 b3 c1 c2 ⃒ c1 b2 − b1 c2
c1 = = c2 = = d1 = =
b1 b1 b1 b1 c1 c1
⃒ ⃒ ⃒ ⃒
⃒b b3 ⃒⃒ ⃒c c2 ⃒⃒ ⃒d d ⃒
− ⃒⃒ 1 − ⃒⃒ 1 ⃒ 1 2⃒
c1 0⃒ d1 d2 ⃒ d1 c2 − c1 d2 ⃒
⃒e 0 ⃒ d2 e1
⃒
d2 = = b3 e1 = = f1 = − ⃒ 1 ⃒=
c1 d1 d1 ⃒ d1 ⃒ d1
⃒ ⃒

18
E. Makena et al. Heliyon 10 (2024) e35038

Section 2.

1 ( ( ( )( ( ) )) ( )( (
a1 = ( ) S0 ϵL0 − L0 r2 − K4 r2 + L0 δ2 + K4 δ2 + − L0 − K4 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3 − S0 L0 + K4 1
S0 L0 + K4
) ) ( )( ( ) ʹ ʹ ʹ))
− S0 k1 r1 − Ω4 + S0 L0 + K4 δ1 + δ3 + k3 1− L0 k3 r3 L0 + k1 r1 S0 − S0 k21 r1 S0

( (
1 ( ) ( )( )( ( ) ) (
a2 = ( 0 ) − AB − L0 − K4 η1 ω1 + S0 − L0 − K4 − r1 + S0 k1 r1 + Ω2 1 − S0 k1 r1 − Ω3 + − L0 ϵ + L0 r2 + K4 r2
S L + K4
) (
)( ( ) ) ( )( (
− L0 δ2 − K4 δ2 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3 − S0 L0 ϵ − L0 r2 − K4 r2 + L0 δ2 + K4 δ2 + − L0 − K4 r1 − S0 k1 r1 + 1
) ( (
) ) (( ) ) ( )( ( )
− S0 k1 r1 − Ω2 − Ω3 1 − S0 k1 r1 − Ω4 + S0 L0 ϵ − L0 r2 − K4 r2 + L0 δ2 + K4 δ2 + − L0 − K4 r1 − S0 k1 r1 + 1 − S0 k1 r1
) )
) ( )( ( ) ) ( ( ) ʹ ʹ ʹ) ( )(
− Ω2 − Ω3 − S0 L0 + K4 1 − S0 k1 r1 − Ω4 δ1 + δ3 + k3 1− L0 k3 r3 L0 + k1 r1 S0 − S0 k21 r1 S0 + S0 L0 + K4 − δ3
)
( ) ʹ)( ʹ ʹ)
− k3 1− L0 k3 r3 L0 − δ1 − k1 r1 S0 + S0 k21 r1 S0

( ( ( ) ( ) ) )) ( )
1
a3 = ( ) τ2 B β L0 + K4 η2 ω2 + − L0 − K4 η1 ω1 ((1− S0 k1 r1 − Ω3 + S0
− L 0
ϵ + L 0
r 2 + K4 r 2 − L 0
δ 2 − K4 δ 2
S0 L0 + K4
( )(( ) ) ( ( ) )( )
− r1 + S0 k1 r1 + Ω2 1 − S0 k1 r1 − Ω3 + τ2 B − L0 − K4 η1 ω1 − S0 ((− L0 − K4 − r1 + S0 k1 r1 + Ω2
(( ) ) ( )( ( ) )))
1 − S0 k1 r1 − Ω3 + − L0 ϵ + L0 r2 + K4 r2 − L0 δ2 − K4 δ2 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3
(( ) ) ( ( ) )( )(( ) )
1 − S0 k1 r1 − Ω4 + − AB − L0 − K4 η1 ω1 + S0 ((− L0 − K4 − r1 + S0 k1 r1 + Ω2 1 − S0 k1 r1 − Ω3
( )( ( ) ))
+ − L0 ϵ + L0 r2 + K4 r2 − L0 δ2 − K4 δ2 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3
( ( )( ( ) ))(( ) ))
− S0 L0 ϵ − L0 r2 − K4 r2 + L0 δ2 + K4 δ2 + − L0 − K4 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3 1 − S0 k1 r1 − Ω4
( ( ) ) ( ( ( )( (
ʹ ʹ ʹ
δ1 + δ3 + k3 1− L0 k3 r3 L0 + k1 r1 S0 − S0 k21 r1 S0 + S0 L0 ϵ − L0 r2 − K4 r2 + L0 δ2 + K4 δ2 + − L0 − K4 r1 − S0 k1 r1 + 1
) )) ( )(( ) ))( ( ) )( ))
ʹ ʹ ʹ
− S0 k1 r1 − Ω2 − Ω3 − S0 L0 + K4 1 − S0 k1 r1 − Ω4 − δ3 − k3 1− L0 k3 r3 L0 − δ1 − k1 r1 S0 + S0 k21 r1 S0

19
E. Makena et al. Heliyon 10 (2024) e35038

( ) ( ) ) )) ( )
1
a4 = ( ) (( − AB(β L0 + K4 η2 ω2 + − L0 − K4 η1 ω1 ((1− S0 k1 r1 − Ω3 − S 0
− L 0
ϵ + L 0
r 2 + K4 r2 − L 0
δ 2 − K 4 δ2
S0 L0 + K4( )(( ) ))(( ) ) ( ( ( )
− r1 + S0 k1 r1 + Ω2 1 − S0 k1 r1 − Ω3 1 − S0 k1 r1 − Ω4 + AB β L0 + K4 η2 ω2 +
( ) ) )) ( )( )(( )
− L0 − K4 η1 ω1 ((1− S0 k1 r1 − Ω3 + S0 − L0 ϵ + L0 r2 + K4 r2 − L0 δ2 − K4 δ2 − r1 + S0 k1 r1 + Ω2 1 − S0 k1 r1
) ( ( ) )( )(( ) )
− Ω3 + AB − L0 − K4 η1 ω1 − S0 ((− L0 − K4 − r1 + S0 k1 r1 + Ω2 1 − S0 k1 r1 − Ω3
( )( ( ) )))(( ) ))
+ − L0 ϵ + L0 r2 + K4 r2 − L0 δ2 − K4 δ2 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3 1 − S0 k1 r1 − Ω4
( ( ) ) ( ( ) )( )((
ʹ ʹ
δ1 + δ3 + k3 1− L0 k3 r3 L0 + k1 r1 S0 − S0 k21 r1 S0 + − AB − L0 − K4 η1 ω1 + S0 ((− L0 − K4 − r1 + S0 k1 r1 + Ω2 1
) ) ( )( ( ) ))
− S0 k1 r1 − Ω3 + − L0 ϵ + L0 r2 + K4 r2 − L0 δ2 − K4 δ2 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3
( ( )( ( ) ))(( ) ))
− S0 L0 ϵ − L0 r2 − K4 r2 + L0 δ2 + K4 δ2 + − L0 − K4 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3 1 − S0 k1 r1 − Ω4
( ( ) )( ))
ʹ ʹ ʹ
− δ3 − k3 1− L0 k3 r3 L0 − δ1 − k1 r1 S0 + S0 k21 r1 S0

( ( ) ( ) ) )) ( )
1
a5 = ( )((− τ2 AB β L0 +K4 η2 ω2 + − L0 − K4 η1 ω1 ((1− S0 k1 r1 − Ω3 − S 0
− L 0
ϵ+L 0
r 2 +K4 r 2 − L 0
δ 2 − K 4 δ2
S0 L0 + K4
( )(( ) ))(( ) )( ( ) )
ʹ ʹ ʹ
− r1 +S0 k1 r1 +Ω2 1 − S0 k1 r1 − Ω3 1 − S0 k1 r1 − Ω4 δ1 + δ3 + k3 1− L0 k3 r3 L0 +k1 r1 S0 − S0 k21 r1 S0
( ( ( ) ( ) ) )) ( )( )
+ τ2 B β L0 +K4 η2 ω2 + − L0 − K4 η1 ω1 ((1− S0 k1 r1 − Ω3 + S0 − L0 ϵ +L0 r2 +K4 r2 − L0 δ2 − K4 δ2 − r1 +S0 k1 r1 + Ω2
(( ) ) ( ( ) )( )(( ) )
1 − S0 k1 r1 − Ω3 + τ2 B − L0 − K4 η1 ω1 − S0 ((− L0 − K4 − r1 + S0 k1 r1 +Ω2 1 − S0 k1 r1 − Ω3
( )( ( ) )))(( ) ))
+ − L0 ϵ+L0 r2 +K4 r2 − L0 δ2 − K4 δ2 r1 − S0 k1 r1 + 1 − S0 k1 r1 − Ω2 − Ω3 1 − S0 k1 r1 − Ω4
( ( ) )( ))
ʹ ʹ ʹ
− δ3 − k3 1− L0 k3 r3 LL0 − δ1 − k1 r1 S0 S0 + S0 k21 r1 S0

( ( ( ) ( ) ) )) (
1
a6 = ( ) − τ2 B β L0 + K4 η2 ω2 + − L0 − K4 η1 ω1 ((1− S0 k1 r1 − Ω3 − S0
− L0 ϵ + L0 r2 + K4 r2 − L0 δ2
S0 L0 + K4
)( )(( ) ))(( ) )( ( ) )(
ʹ
− K4 δ2 − r1 + S0 k1 r1 + Ω2 1 − S0 k1 r1 − Ω3 1 − S0 k1 r1 − Ω4 − δ3 − k3 1− L0 k3 r3 L0 − δ1
)
ʹ ʹ
− k1 r1 S0 + S0 k21 r1 S0

( )
By Routh-Hurwitz criteria for stability, system (1)–(6) is locally asymptomatically stable at disease-free equilibrium E0 if and only if
a1 > 0, a2 > 0, a4 > 0, a6 > 0, b1 > 0, c1 > 0, d1 > 0 and e1 > 0 are satisfied and otherwise unstable.

Section 3.

The bifurcation results of a and b values are as follows;

20
E. Makena et al. Heliyon 10 (2024) e35038

( ( ) ( ) )
2(− AB + u2 )v3 η2 ω2 − ϵf6
f6 +K4
+ 1− f5
K2
r2 − δ2 + 1 − Kf52 η1 ω1 − f5 r2 +f2 η1 ω1 +fK32 η2 ω2 +f4 η3 ω3 (βu3 − u5 [η1 ω1 ])
a= − ( )( )
ϵ
1+k4
+ r 2 − δ 2
ε
1+K4
+ r 2 − δ2 ( − 1 + β + k1 − δ1 − η1 + Ω1 )
( ( ) ( ) )
2(− AB + u2 )v3 η2 ω2 − ϵf6
f6 +K4
+ 1 − Kf52 r2 − δ2 + 1 − Kf52 η2 ω2 − f5 r2 +f2 η1 ω1 +fK32 η2 ω2 +f4 η3 ω3 u5 [η3 ω3 ]
+ ( )( )
ϵ
1+k4
+ r 2 − δ 2
ε
1+K4
+ r 2 − δ 2 ((1 − k1 )r1 − η3 )
( ( ) ( ) )
f5 f5 f5 r2 +f2 η1 ω1 +f3 η2 ω2 +f4 η3 ω3
2(− AB + u2 )v3 η2 ω2 − ϵf6
f6 +K4
+ 1 − K2
r2 − δ 2 + 1 − K2
η 3 ω3 − K2
u5 [η3 ω3 ]
+ ( )( )
ϵ
1+k4
+ r2 − δ2 ε
1+K4
+ r2 − δ2 ((1 − k1 )r1 − η4 )
( [ ] )
2ϴf33 f4 2ϴf3 f4
2ABv5 β + 2 − D2 +f32
+1 1− N
K1
r1 − δ1 − η2 (βu3 − u5 [η1 ω1 ])u5 [η3 ω3 ]
(D2 +f32 )
−
((1 − k1 )r1 − η3 )2 ( − 1 + β + k1 − δ1 − η1 + Ω1 )
( )
(( ) ( ) )
ϴf32
2AB β − D2 +f 2 v5 (βu3 − u5 [η1 ω1 ])u5 [η3 ω3 ] 2v3 η ω2 1 − f5
η ω 1 + 1 − f5
η ω 2 (βu3 − u5 [η1 ω1 ])u5 [η3 ω3 ]
3 2 K2 1 K2 2
− + ( )
((1 − k1 )r1 − η3 ) ε
1+K4
+ r2 − δ2 ((1 − k1 )r1 − η3 )( − 1 + β + k1 − δ1 − η1 + Ω1 )
(( ) ( ) )
2v3 η2 ω2 1 − Kf52 η1 ω1 + 1 − Kf52 η3 ω3 (βu3 − u5 [η1 ω1 ])u5 [η3 ω3 ]
+ ( )
ε
1+K4
+ r2 − δ2 ((1 − k1 )r1 − η4 )( − 1 + β + k1 − δ1 − η1 + Ω1 )
( [ ] )
ϴf32 2ϴf33 f4
(( ) ( ) )
2ABv5 − D2 +f 2 + 2 2 2 − D2 +f 2 + 1 1 − K1 r1 − δ1 − η2 u5 [η3 ω3 ]2 2v η ω
2ϴf3 f4 N
1 − Kf52 η2 ω2 + 1 − Kf52 η3 ω3 u5 [η3 ω3 ]2
3 (D +f3 ) 3 3 2 2
+ − ( )
((1 − k1 )r1 − η3 )2 ((1 − k1 )r1 − η4 ) ε
+ r2 − δ2 ((1 − k1 )r1 − η3 )((1 − k1 )r1 − η4 )
1+K4
⎛ [ ]⎞
f +f +f +f
r1 f1ʹ 1− 1 2K 3 4
⎫
⎜ eM τ2 Bf1
1 ⎟⎪
2(− τ2 B + u2 )( − k1 r1 v4 + τ2 Bv5 + v2 (k1 r1 − η1 ))u5 [η1 ω1 ]⎜ +
⎝F+M (f1 +Bf5 ) 2 − K1
⎟⎪
⎪
⎠⎪
⎪
⎪
⎪
⎬
−
((1 − k1 )r1 − k1 r1 − δ1 ) ⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭

21
E. Makena et al. Heliyon 10 (2024) e35038

( [ ])
ABf1 f1 +f2 +f3 +f4
(− AB + u2 )(k1 r1 v3 + ABv5 ) 2 + f1 1 − K1 [
2ϴf33 f4 2ϴf3 f4 (f1 +Bf5 )
b=( )2 − D2 + f 2 − δ1 − η2 + + f3 1
D2 + f32 3 ((1 − k1 )r1 − k1 r1 − δ1 )
[ ] ( ( ) )
] AB(− AB + u2 )v5 f3 1 − f1 +f2K+f3 +f4 AB β + f 3 1 − f1 +f2 +f3 +f4
− f3 r1
v5 (βu3 − u5 [η1 ω1 ])
f1 + f2 + f3 + f4 1 K1 K1
− +( ) +
K1 ((1 − k1 )r1 − η3 )
ϵ
1+k4
+ r2 − δ2 ((1 − k1 )r1 − η3 )
⎛ [ ]⎞
[ ] r1 f1ʹ 1− f1 +f2K+f3 +f4
⎜ f1 +f2 +f3 +f4
1 ⎟
(k1 r1 v3 + ABv5 )u5 [η3 ω3 ]⎜
⎝f1 1 − K1
− K1
⎟
⎠

−
((1 − k1 )r1 − k1 r1 − δ1 )
⎛ [ ]⎞
[ ] f1 +f2 +f3 +f4
r1 f1ʹ 1− K1
⎜ f1 +f2 +f3 +f4
⎟
(k1 r1 v3 + ABv5 )u5 [η3 ω3 ]⎜
⎝f1 1 − K1
− K1
⎟
⎠

−
((1 − k1 )r1 − k1 r1 − δ1 )
⎛ [ ]⎞
[ ] f1 +f2 +f3 +f4
r1 f1 ʹ 1− K1
⎜ eM f1 +f2 +f3 +f4
⎟
(k1 r1 v3 + ABv5 )(βu3 − u5 [η1 ω1 ])⎜
⎝F+M
+ f1 1− K1
− K1
⎟
⎠

+
((1 − k1 )r1 − k1 r1 − δ1 )
⎛ [ ]⎞ ⎛ [ ]⎞
[ ] f3ʹ 1− f1 +f2K+f3 +f4 [ ] r1 f3ʹ 1−
f1 +f2 +f3 +f4
K1
⎜ f1 +f2 +f3 +f4
1 ⎟ ⎜ ϴf32 f1 +f2 +f3 +f4
⎟
ABr1 v5 u5 [η3 ω3 ]⎜
⎝1 1 − K1
− K1
⎟
⎠ ABv5 u5 [η3 ω3 ]⎜
⎝− D2 +f32
+ f3 1 − K1
− K1
⎟
⎠

− −
((1 − k1 )r1 − η3 )2 ((1 − k1 )r1 − η3 )

Section 4.

b4 = 1;
⎛ ⎞
⎛ ⎞
⎜ 1 ⎟
S⎝ MQ +β+δ
1
+ 1
δ1 +η2
+ δ1 +η3 ⎠⎟
⎜ F+M 1 +η1
⎜
⎜ 1 1 1 ⎟
⎟ r2
b3 = r1 ⎜ − − − + ⎟−
⎜ MQ
F+M
+ β + δ1 + η1 δ1 + η2 δ1 + η3 K1 ⎟ Lϵ
L+K4
+ δ2 + Sτ1
⎝ ⎠

⎛ ⎞
⎛ ⎞
⎛ ⎜ 1 1 1
⎟ ⎞
MQ S⎝− − δ1 +η2 − δ1 +η3 ⎠
⎜ F+M
+β+δ1 +η1 ⎟
⎜ ⎜ ⎟⎟
⎜ ⎟⎟
⎜
⎜ r2 ⎜ MQ +β+δ +η + δ1 +η + δ1 +η +
1 1 1
K1 ⎟⎟
⎜F+M 1 1 2 3 ⎟⎟
⎜ ⎝ ⎠⎟
⎜
⎜(S − K )2 r (Fβ + M(Q + β) + (F + M)(3δ + η + η + η )) ⎟
⎜ 1 1 1 1 2 3 ⎟
b2 = r1 ⎜ 2 + ⎟
⎜ K1 (Fβ + M(Q + β) + (F + M)(δ1 + η1 ))(δ1 + η2 )(δ1 + η3 ) Lϵ
L+K4
+ δ2 + S τ 1 ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎝ ⎠

22
E. Makena et al. Heliyon 10 (2024) e35038

r13 S3 r13
b1 = − ( ) + ( )
MQ
F+M
+ β + δ1 + η1 (δ1 + η2 )(δ1 + η3 ) K31 MQ
F+M
+ β + δ1 + η1 (δ1 + η2 )(δ1 + η3 )

3S2 r13 3Sr13

− ( ) + ( )
K21 MQ
F+M
+ β + δ1 + η1 (δ1 + η2 )(δ1 + η3 ) K1 F+M
MQ
+ β + δ1 + η1 (δ1 + η2 )(δ1 + η3 )

r12 r2 S2 r2 r
− ( ) ( )− ( ) 1 2 ( )
MQ
F+M
+ β + δ1 + η1 Lϵ
(δ1 + η2 ) L+K 4
+ δ 2 + Sτ 1 K2 MQ
1 F+M + β + δ1 + η 1 (δ 1 + η 2 ) Lϵ
L+K4
+ δ 2 + S τ 1

2Sr2 r r2 r
+ ( ) 1 2 ( )− ( ) 1 2 ( )
K1 MQ
F+M
Lϵ
+ β + δ1 + η1 (δ1 + η2 ) L+K 4
+ δ2 + S τ 1 MQ
F+M
Lϵ
+ β + δ1 + η1 (δ1 + η3 ) L+K 4
+ δ2 + Sτ1

S2 r 2 r 2Sr2 r
− ( ) 1 2 ( )+ ( ) 1 2 ( )
K21 MQ
F+M
Lϵ
+ β + δ1 + η1 (δ1 + η3 ) L+K 4
+ δ2 + Sτ 1 K1
MQ
F+M
+ β + δ 1 + η1 (δ 1 + η3 ) Lϵ
L+K 4
+ δ2 + Sτ 1

r12 r2 S2 r12 r2 2Sr12 r2

− ( )− ( )+ ( ).
(δ1 + η2 )(δ1 + η3 ) L+KLϵ
4
+ δ 2 + Sτ 1 K 2
1 (δ1 + η 2 )(δ 1 + η3 ) Lϵ
L+K4
+ δ2 + Sτ 1 K 1 (δ 1 + η 2 )(δ 1 + η3 ) Lϵ
L+K4
+ δ2 + Sτ 1

( )3 ( )
(F + M) S0 − K1 L0 + K4 r13 r2
b0 = ( ( ) )
K31 (Fβ + M(Q + β) + (F + M)(δ1 + η1 ))(δ1 + η2 )(δ1 + η3 ) L0 ϵ + L0 + K4 (δ2 + Sτ1 )

References

[1] S. Chakraborty, X.Z. Li, P.K. Roy, How can HPV-induced cervical cancer be controlled by a combination of drug therapy? A mathematical study, Int. J. Biomath.
(IJB) 12 (6) (2019), https://doi.org/10.1142/S1793524519500700.
[2] J.C. Sierra-Rojas, R. Reyes-Carreto, C. Vargas-De-León, J.F. Camacho, Modeling and mathematical analysis of the dynamics of HPV in cervical epithelial cells:
transient, acute, latency, and chronic infections, Comput. Math. Methods Med. 2022 (2022), https://doi.org/10.1155/2022/8650071.
[3] R.P. Insinga, E.J. Dasbach, E.H. Elbasha, Epidemiologic natural history and clinical management of Human Papillomavirus (HPV) Disease: a critical and
systematic review of the literature in the development of an HPV dynamic transmission model, BMC Infect. Dis. 9 (2009) 1–26, https://doi.org/10.1186/1471-
2334-9-119.
[4] A.A. Alsaleh, A.B. Gumel, Dynamics analysis of a vaccination model for HPV transmission 22 (4) (2014), https://doi.org/10.1142/S0218339014500211.
[5] L. Ribassin, A SIS Model for Human Papillomavirus Transmission to Cite This Version :, 2012.
[6] E.D. Gurmu, P.R. Koya, Impact of chemotherapy treatment on SITR compartmentalization and modeling of human papilloma virus 15 (3) (2019) 17–29, https://
doi.org/10.9790/5728-1503011729.
[7] G. F., M. F., A. M., and B. Z., Comprehensive knowledge about cervical cancer is low among women in Northwest Ethiopia, BMC Cancer 13 (2013), https://doi.
org/10.1186/1471-2407-13-2%5Cn [Online]. Available: http://www.embase.com/search/results?subaction=viewrecord&from=export&id=L52379507%5Cn
http://www.biomedcentral.com/1471-2407/13/2%5Cn http://wt3cf4et2l.search.serialssolutions.com?sid=EMBASE&issn=14712407&id=doi.
[8] H. Stark, A. Zivković, HPV vaccination: prevention of cervical cancer in Serbia and in europe, Acta Fac. Med. Naissensis 35 (1) (2018) 5–16, https://doi.org/
10.2478/afmnai-2018-0001.
[9] T.A. Kessler, Cervical cancer: prevention and early detection, Semin. Oncol. Nurs. 33 (2) (2017) 172–183, https://doi.org/10.1016/j.soncn.2017.02.005.
[10] P. Online, P.M. Ezzati, Health Policy NCD Countdown 2030 : Pathways to Achieving Sustainable Development Goal Target 3 . 4, vol. 396, 2020, https://doi.org/
10.1016/S0140-6736(20)31761-X.
[11] R.W. Menzel, Further notes on the albino catfish, J. Hered. 49 (6) (1958) 159–178, https://doi.org/10.1093/oxfordjournals.jhered.a106828.
[12] R.J. Smith, J. Li, D. Blakeley, The failure of R0, Comput. Math. Methods Med. 2011 (2011), https://doi.org/10.1155/2011/527610.
[13] E. Dadi Gurmu, P. Rao Koya, Sensitivity analysis and modeling the impact of screening on the transmission dynamics of human papilloma virus (HPV), Am. J.
Appl. Math. 7 (3) (2019) 70, https://doi.org/10.11648/j.ajam.20190703.11.
[14] V.V. Akimenko, F. Adi-Kusumo, Stability analysis of an age-structured model of cervical cancer cells and HPV dynamics, Math. Biosci. Eng. 18 (5) (2021)
6155–6177, https://doi.org/10.3934/mbe.2021308.
[15] S. Chakraborty, X. Cao, S. Bhattyacharya, P.K. Roy, The role of HPV on cervical cancer with several functional response: a control based comparative study,
Comput. Math. Model. 30 (4) (2019) 439–453, https://doi.org/10.1007/s10598-019-09469-4.
[16] T.S.N. Asih, et al., The dynamics of HPV infection and cervical cancer cells, Bull. Math. Biol. 78 (1) (2016) 4–20, https://doi.org/10.1007/s11538-015-0124-2.
[17] M.Z. Ndii, Mathematical Model of Cervical Cancer Treatment Using Chemotherapy Drug, February, 2020, pp. 10–16, https://doi.org/10.14421/
biomedich.2019.81.11-15.
[18] P. Pongsumpun, Mathematical Model of Cervical Cancer Due to Human Papillomavirus Infection, 2012, pp. 157–161, no. January.
[19] U. Sulaiman, I.I. Adamu, A. Tahir, Mathematical Model on the Impact of Natural Immunity, Vaccination, Screening and Treatment on the Dynamics of the
Human Papillomavirus Infection and Cervical Cancer, 2016. February 2018.
[20] K. Allali, Stability analysis and optimal control of HPV infection model with early-stage cervical cancer, Biosystems 199 (December 2020) (2021) 104321,
https://doi.org/10.1016/j.biosystems.2020.104321.
[21] B. Ma, J.J. Thibodeaux, Well-Posedness and a Finite Difference Approximation for a Mathematical Model of HPV-Induced Cervical Cancer, 2023, pp. 151–172,
https://doi.org/10.4236/am.2023.143009.
[22] R.P. Mondaini, C. Lectures, R. De Janeiro, Trends in Biomathematics: Chaos and Control in Epidemics, Ecosystems, and Cells, 2021, https://doi.org/10.1007/
978-3-030-73241-7.
[23] S.H. Erwin, Mathematical models of immune responses to infectious diseases [Online]. Available: http://search.ebscohost.com/login.aspx?
direct=true&db=ddu&AN=49D017A07C5B90D3&site=ehost-live, 2017.

23
E. Makena et al. Heliyon 10 (2024) e35038

[24] M. Verma, et al., Modeling the mechanisms by which HIVAssociated immunosuppression influences HPV persistence at the Oral Mucosa, PLoS One 12 (1)
(2017) 1–20, https://doi.org/10.1371/journal.pone.0168133.
[25] B. Belew, D. Melese, Modeling and analysis of predator-prey model with fear effect in prey and hunting cooperation among predators and harvesting, J. Appl.
Math. 2022 (2022), https://doi.org/10.1155/2022/2776698.
[26] G.K. Mutua, C.G. Ngari, G.G. Muthuri, D.M. Kitavi, Mathematical modeling and simulating of Helicobacter pylori treatment and transmission implications on
stomach cancer dynamics, Commun. Math. Biol. Neurosci. 2022 (2022) 1–29, https://doi.org/10.28919/cmbn/7542.
[27] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in
heterogeneous populations, J. Math. Biol. 28 (4) (1990) 365–382, https://doi.org/10.1007/BF00178324.
[28] C.L. Shaw, D.A. Kennedy, What the reproductive number R0 can and cannot tell us about COVID-19 dynamics, Theor. Popul. Biol. 137 (2021) 2–9, https://doi.
org/10.1016/j.tpb.2020.12.003.
[29] J.S. Kilonzi, C.G. Ngari, S. Karanja, Modelling the Impact of Screening , Treatment and Underlying Health Conditions on Dynamics of Covid-19, January 2023,
2024.
[30] J. Ochwach, M.O. Okongo, C. Ngari, Mathematical Modelling of Cholera Incorporating the Dynamics of the Induced Achlorhydria Condition and Treatment,
November, 2022.
[31] N.C. Gitonga, Modelling Childhood Pneumonia and its Implications Regarding its Control Using Kenyan Data, 2017, pp. 1–120.

24