
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)11069-9
10.1016/j.heliyon.2024.e35038
e35038
Research Article
In-host density-dependent model of high-risk HPV virions, basal cells, lymphocytes t-cells incorporating functional responses
Makena Elosy elseakesh@gmail.com
a⁎
Ngari Cyrus Gitonga ngaricyrus15@gmail.com
b
Kimani Patrick Mwangi patkim32@gmail.com
a
Kilonzi Jeremiah Savali jeremiahkilonzi66@gmail.com
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
⁎ Corresponding author. elseakesh@gmail.com
25 7 2024
30 8 2024
25 7 2024
10 16 e3503816 2 2024
6 7 2024
22 7 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Cervical cancer is one of the most common types of cancer and it is caused mostly by high-risk Human Papillomavirus (HPV) and continues to spread at an alarming rate. While HPV impacts have been investigated before, there are currently only a scanty number of mathematical models that account for HPVˈs dynamic role in cervical cancer. The objectives were to develop an in-host density-dependent deterministic model for the dynamics implications of basal cells, virions, and lymphocytes incorporating immunity and functional responses. Analyze the model using techniques 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 R0*<1 and globally asymptotically stable if R0*≤1 and the Endemic Equilibrium is globally asymptotically stable if R0*>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%.

Keywords

In-host model
Functional responses
Stability analysis
Simulation and reproduction number
==== Body
pmc1 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, 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 Fig. 7, Fig. 8, Fig. 9, Fig. 10, Fig. 11, Fig. 12, Fig. 13, Fig. 14, Fig. 15).

The dynamics of HPV to cervical cancer are clearly shown in Fig. 1 [8].Fig. 1 Progression of cervical cancer by HPV.

Fig. 1

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,SIHPVIUCICC,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], [16], [17], [18], [19], [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.

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-Susceptible (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.i. τ2SS+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.

ii. CP2D2+P2 is a saturation function in which precancerous cells progress to cancerous cells based on Holling type III response.

iii. EMF+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 likelihood estimation method (see Table 2).Table 1 Parameter and state variables description.

Table 1Parameter	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			

Table 2 Table of sensitivity indices.

Table 2Parameter	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	

The equations are derived from the model flow chart as follows;(1) dSdt=r1S(1−NK1)−τ2SVS+BV+EMF+MI−δ1S,

(2) dIdt=r1I(1−NK1)+τ2SVS+BV−βI−EMF+MI−η1I−δ1I,

(3) dPdt=r1P(1−NK1)+βI−ϴCP2D2+P2−η2P−δ1P,

(4) dCdt={r1C(1−NK1)+ϴCP2D2+P2−η3C−δ1C}(1−CK4),

(5) dVdt={r2V+ω1η1I+ω2η2P+ω3η3C}(1−VK2)−τ1SV−ϵLL+K4V−δ2V,

(6) dLdt=r3L(1−LK3)−α1LV−δ3L.

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), (2), (3), (4), (5), (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 R6+ [6].

Proof .

For the sake of analysis(7) dSdt=r1S(1−k1N)−τ2SVS+BV+Ω1I−δ1S,

(8) dIdt=r1I(1−k1N)+τ2SS+BV−Ω2I,

(9) dPdt=r1P(1−k1N)+βI−ϴCP2D2+P2−Ω3P,

(10) dCdt={r1C(1−k1N)+ϴCP2D2+P2−Ω4C}(1−k4C),

(11) dVdt={r2V+ω1η1I+ω2η2P+ω3η3C}(1−k2V)−τ1SV−ϵLL+K4V−δ2V,

(12) dLdt=r3L(1−k3L)−α1LV−δ3L,

Where 1K1=k1;1K2=k2;1K3=k3;EMF+M=Ω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)=r1S(1−k1N)−τ2SVS+BV+Ω1I−δ1S, F2(W)=r1I(1−k1N)+τ2SVS+BV−Ω2I,F3(W)=r1P(1−k1N)+βI−ϴCP2D2+P2−Ω3P,F4(W)={r1C(1−k1N)+ϴCP2D2+P2−Ω4C}(1−k4C),F5(W)={r2V+ω1η1I+ω2η2P+ω3η3C}(1−k2V)−τ1SV−ϵLL+K4V−δ2V,F6(W)=r3L(1−k3L)−α1LV−δ3L. Then, system (1)–(6) can be written as dQdt=F(W) where F:C+→(R)+6 with W(0)=W0ϵR+6. Thus, the function W is locally Lipschitzian and completely stable on R+6.

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 R+6 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), (2), (3), (4), (5), (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 (13) N(t)=S(t)+I(t)+P(t)+C(t)

By differentiating equation (13) gives,(14) dNdt=dSdt+dIdt+dPdt+dCdt

(15) dNdt=r1S(1−NK1)−δ1S+r1I(1−NK1)−η1I−δ1I+r1P(1−NK1)−η2P−δ1P+r1C(1−NK1)−η3C−δ1C(1−CK4)

(16) dNdt=r1N(1−NK1)−(δ1S+η1I+δ1I+η2P+δ1P+η3C+δ1C(1−CK4))

dNdt≤r1N(1−NK1)

By using the separation of variables of inequality, we have(17) dNN(1−NK1)≤r1dt

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(18) N(t)≤K11+(K1N(0)−1)e−rt.

Introducing limits, limt→∞N(t)≤K1.Implying that 0≤N(t)≤K1.

Also for, V(t).(19) dVdt={r2V+ω1η1I+ω2η2P+ω3η3C}(1−VK2)−τ1SV−ϵLL+K4V−δ2V.

(20) Letr2V+ω1η1I+ω2η2P+ω3η3C=r4V

(21) dVdt=r4V(1−VK2)−τ1SV−ϵLL+K4V−δ2V,then,dVdt≤r4V(1−VK2)

By separating variables (21), integrating and applying initial conditions and limits we obtained limt→∞V(t)≤ K2.Therefore; V(t)≤K2 .Implying that, 0≤V(t)≤K2.

Similarly, L(t), dLdt≤r3L(1−LK3). Therefore; L(t)≤K31+(K3L(0)−1)e−rt. On applying the limits as t→∞, it follows that; 0≤L(t)≤K3.

3.1.3 Disease-free equilibrium point (DFE)

Theorem 3 The system of equations (1), (2), (3), (4), (5), (6) has disease-free equilibrium point (E0) obtained as; E0=(S0,I0,P0,C0,V0,L0)=(K1(r1−δ1)r1,0,0,0,0,K3r3(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=0andV=0 [26].

Solving, we obtained. S0=0 and S0=K1(r1−δ1)r1 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(r1−δ1)r1 and L0=K3(r3−δ3)r3. Hence the disease free equilibrium point of the system of equations (1), (2), (3), (4), (5), (6), (7), (8) was obtained as:E0=(S0,I0,P0,C0,V0,L0)=(K1(r1−δ1)r1,0,0,0,0,K3r3(r3−δ3)).

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.

3.1.4.1 Using next generation matrix

Theorem 4 The basic reproduction number ,R0=(M+F)(K1−S0)γ1K1(Mβ+Fβ+ME+Mδ1+Fδ1+Mη1+Fη1) by Next-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. Thusf=[r1I(1−NK1)+τ2SVS+BVr1P(1−NK1)r1C(1−NK1)(1−CK4)r1V+ω1η1I+ω2η2P+ω3η3C(1−VK2)];v=[βI+EMF+MI+η2I+δ1I−βI+θCP2D2+P2+η2P+δ1P(−θCP2D2+P2+η3C+δ1C)(1−CK4)τ1SV+ϵLVL+K4+δ2V]

At E0 point, Jacobian matrices of f and v was evaluated to find out matrices F and V respectively,F=[(1−S0K1)γ100−βS0Vτ2S20(1−S0K1)γ10000(1−S0K1)γ10η1ω1η2ω2η3ω3γ2];V=[β+EMF+M+δ1+η1000−βδ1+η20000δ1+η30000L0ϵL0+K5+δ2+τ1S0V]

At E0 point, FV−1, was obtained as, [(1−S0K1)r1MQF+M+β+δ1+η10000(1−S0K1)r1δ1+η20000(1−S0K1)r1δ1+η30000r2L0ϵL0+K5+δ2+τ1S0V];

The four eigenvalues of FV−1 at E0 were: (K1−S0)γ1K1(δ1+η3), (K1−S0)γ1K1(δ1+η2),(M+F)(K1−S0)γ1K1(Mβ+Fβ+ME+Mδ1+Fδ1+Mη1+Fη1)and(L0+K4)γ2ϵL0+δ2L0+K4δ2.

By inspection method (which was be verified by numerical method too) the dominant eigenvalue represents the R0 [27]. HenceR0=(M+F)(K1−S0)γ1K1(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 numberR0=r1(1−NK1)+τ2SVS+BVβ+EMF+M+η1+δ1+r1(1−NK1)θCPD2+P2+η2+δ2+r1(1−NK1)((η3+δ1)(1−CK4))+(r2V+ω1η1I+ω2η2P+ω3η3C)(1−VK2)ϵLL+K4+δ2+τ1SV

by the method proposed by Ref. [28]. Where N(0),S(0),V(0),P(0),C(0),I(0)andL(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,IandL are assumed to be constant terms at the beginning of the epidemic, evaluating R0 of equations (2), (3), (4), (5), the basic reproduction number was obtained by summing them [12] that is,{r1(1−NK1)+τ2SVS+BV}∫0∞e−(β+EMF+M+η1+δ1)TAdTA+{r1(1−NK1)}∫0∞e−(θCPD2+P2+η2+δ2)TAdTA+{r1(1−NK1)}∫0∞e−((η3+δ1)(1−CK4))TAdTA+{(r2V+ω1η1I+ω2η2P+ω3η3C)(1−VK2)}∫0∞e−(ϵLL+K4+δ2+τ1SV)TAdTA

On integration, {r1(1−NK1)+τ2SVS+BV}[−e−(β+EMF+M+η1+δ1)TAβ+EMF+M+η1+δ1]0∞+{r1(1−NK1)}[−e−(θCPD2+P2+η2+δ2)TAθCPD2+P2+η2+δ2]0∞+{r1(1−NK1)}[−e−((η3+δ1)(1−CK4))TA((η3+δ1)(1−CK4))]0∞+{(r2V+ω1η1I+ω2η2P+ω3η3C)(1−VK2)}[e−(ϵLL+K4+δ2+τ1SV)TAϵLL+K4+δ2+τ1SV]0∞. Putting the limits the basic reproduction number is obtained as:R0=r1(1−NK1)+τ2SVS+BVβ+EMF+M+η1+δ1+r1(1−NK1)θCPD2+P2+η2+δ2+r1(1−NK1)((η3+δ1)(1−CK4))+(r2V+ω1η1I+ω2η2P+ω3η3C)(1−VK2)ϵLL+K4+δ2+τ1SV

3.1.4.3 Basic reproduction number by evaluating the constant term of the characteristic polynomial

Theorem 6 The basic reproduction number

R0=(F+M)(S0−K1)3(L0+K4)r13r2K13(Fβ+M(Q+β)+(F+M)(δ1+η1))(δ1+η2)(δ1+η3)(L0ϵ+(L0+K4)(δ2+Sτ1)) by method evaluating the constant term of the characteristic polynomial [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,b1andb0 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].R0=(F+M)(S0−K1)3(L0+K4)r13r2K13(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*)

From; {r1C*(1−NK1)+ϴC*P*2D2+P*2−η3C*−δ1C*}(1−C*K4)=0;

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.

From; r1P*(1−NK1)+βI*−ϴC*P*2D2+P*2−η2P*−δ1P*=0;

It implies that; P*=−CθK1±C2θ2K12−4(βK1−Nr1+K1r1−K1δ1−K1η2)(D2βK1−D2Nr1+D2K1r1−D2K1δ1−D2K1η2)2(−βK1+Nr1−K1r1+K1δ1+K1η2);

From; {r2V*+ω1η1I*+ω2η2P*+ω3η3C*}(1−V*K2)−τ1SV−ϵL*L*+K4V*−δ2V*=0V*=K2(LϵK4−r2+Jη1ω1K2+Pη2ω2K2+Tη3ω3K2+4r2(Jη1ω1+Pη2ω2+Tη3ω3)K2+(−LϵK4+r2−Jη1ω1K2−Pη2ω2K2−Tη3ω3K2)2)2r2

Substituting the value of P*,C*andV* to;

r1I*(1−NK1)+τ2S0VS+BV*−βI*−EMF+MI*−η1I*−δ1I*>0 and solving for I*, we get the value of I*= −τ2M(S+BV)(−EMF+M−β+(1−XK1)r1−δ1−η1);From r3L*(1−L*K3)−α1L*V*−δ3;L*=K3(r3−Vδ1−δ3)r3

Theorem 7 The necessary and sufficient conditions for existence is R0>1, dIdt>0,dPdt>0,dCdt>0anddVdt>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;[R1R2R3R4R5R6]T

R1=[−δ1−k1r1S0(1−k1S0)Ω1−k1r1S0(1−k1S0)−k1r1S0(1−k1S0)−k1r1S0(1−k1S0)τ2BS00]

R2=[0(1−k1S0)r1−Ω200−τ2BS00]

R3=[0β(1−k1S0)r1−Ω3000]

R4=[000(1−k1S0)r1−Ω400]

R5=[0η1ω1η2ω2η3ω3−ϵL0L0+K4+r2−δ2−τ1S00]

R6=[0000−α1L0−δ3−k3r3L0(1−k3L0)]

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, d1ande1 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, d1ande1 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. dXdt=F(X,Z); dZdt=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.i. dXdt=(X,0), Where X∗ is globally asymptotically stable.

ii. dZdt=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 theorems hold.

Theorem; The equilibrium point Eo=(X*,0) of the system [[1], [2], [3], [4], [5], [6]] is globally asymptotically stable if R0*≤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 dXdt=(X,0)=0 Letting A=DzG(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;AZ=[(β1S0−k1)E+η1β1S0I+η2β1S0T+η3β1S0CρE−k2IωI−k3TγI+αT−k4C]G(X,Z)=[(λ1+λ2)S−k1EρE−k2IωI−k3TγI+αT−k4C]

But Gˆ(X,Z)=AZ−G(X,Z), This reduces to Gˆ(X,Z)=[G1ˆ(X,Y)G2ˆ(X,Y)G3ˆ(X,Y)G4ˆ(X,Y)]=[(β1(E+η1I+η2T+η3C)(S0−S)−λ2S000];

Thus, if then the disease-free equilibrium (E0) is globally stable and unstable otherwise. The susceptible is bounded as, S≤S0.Thus, DFE, E0 is globally asymptomatically stable if and only if (β1(E+η1I+η2T+η3C)(S0−S)≥λ2S [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−xi*lnxi)

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 xi* is the endemic equilibrium point.L=b1(S−S*lnS)+b2(I−I*lnI)+b3(P−P*lnP)+b4(C−C*lnC)+b5(V−V*lnV)+b3(L−L*lnL)

dLdt=b1(1−S*S)dSdt+b2(1−I*I)dIdt+b3(1−C*C)dCdt+b4(1−P*P)dPdt+b5(1−V*V)dVdt+b6(1−L*L)dLdt

dLdt=b1(1−S*S){r1S(1−k1N)−τ2SS+BV+Ω1I−δ1S}+b2(1−I*I){r1I(1−k1N)+τ2SS+BV−Ω2I}+b3(1−C*C){{r1C(1−k1N)+ϴCP2D2+P2−Ω4C}(1−k4C)}+b4(1−P*P){r1P(1−k1N)+βI−ϴCP2D2+P2−Ω3P}+b5(1−V*V){{r2V+ω1η1I+ω2η2P+ω3η3C}(1−k2V)−τ1SV−ϵLL+K4V−δ2V}+b6(1−L*L){r3L(1−k3L)−α1LV−δ3L}

After solving we obtain the value of X and Y as follows.X=b1r1S+b1Ω1I+bS*r1K1N+bS*τ2S+BV+bS*δ1+b2r1I+b2τ2S+b2I*K1N+b2I*Ω2+b3r1C+b3θCP2D2+P2+b3r1C2K1K4N+b3Ω4C2K4+b3K1NC*+b3Ω4C*+b3K4CC*r1+b3θC*CP2K4D2+P2+b4r1P+b4βI+b4r1PK1NK4C+b4θCP2D2+P2K4C+b4Ω4C2K4+b4P*P(r1PK1N+θCP2D2+P2+Ω4C+r1PK4C+βICK4+b5(r2V+ω1η1I+ω2η2P+ω3η3C+ω2η2PK2V+b5V*V(r2K2V2+ω1η1IK2V+ω3η3CK2V+τ1SV+ϵLVL+K4+δ2V)+b6r3L+b6L*L(r3L2K3+α1LV+δ3L)

Y=−br1SK1N−bSτ2S+BV−bSδ1−bS*r1−bS*Ω1IS−b2r1IK1N−b2Ω2I−b2I*r1−b2I*ISτ2S+BV−b3r1CK1N−b3Ω4C−b3r1C2K4−b3θC2P2K4D2+P2−b3C*−b3θC*P2D2+P2−b3C*Ω4CK4−b3C*K1CK4r1−b4(r1PK1N+θCP2D2+P2+Ω4C+r1PK4C+βIK4C)−b4P*P(r1P+βI+r1PK1NK4C+θCP2D2+P2K4C)−b5(r2K2V2+ω1η1IK2V+ω3η3CK2V+τ1SV+ϵLVL+K4+δ2V)−b5V*V(r2V+ω1η1I+ω2η2P+ω3η3C+ω2η2PK2V)−b6(r3L2K3+α1LV+δ3L)−b6L*r3

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,

dxdt=F(x) where F=(f1,f2,f3,f4,f5,f6)T, It follows that:dx1dt=p1=r1f1(1−NK1)−τ2f1f5f1+Bf5+EMF+Mf2−δ1f1

dx2dt=p2=r1f2(1−NK1)+τ2f1f5f1+Bf5−βf2−EMF+Mf2−η1f2−δ1f2

dx3dt=p3=r1f3(1−NK1)+βf2−ϴf4f32D2+f32−η2f3−δ1f3.

dx4dt=p4={r1f4(1−NK1)+ϴf4f32D2+f32−η3f4−δ1f4}(1−f4K4).

dx5dt=p5={r2f5+ω1η1f2+ω2η2f3+ω3η3f4}(1−f5K2)−τ1f1f5−ϵf6f6+K4f5−δ2f5.

dx6dt=p6=r3f6(1−f6K3)−α1f6f5−δ3f6.

Jacobian solved at DFE,

E0 = (S0, I0, P0, C0, V0, L0) = (K1(r1−δ1)r1,0,0,0,0,K3r3(r3−δ3)), in a case, where R0*=1 and further, suppose that r1=r1* is a bifurcation parameter, then solving for r1* from R0*=1 we get,r1*=K1(MQ+Fβ+Mβ+Fδ1+Mδ1+Fη1+Mη1)(F+M)(K1−S0)

Gives[A1A2A3A4]

Where;A1=[(1−S0K1)r1*−r1*S0K1−δ1MEF+M−r1*S0K1−r1*S0K10−MEF+M−β+(1−S0K1)r1*−δ1−η100β(1−S0K1)r1*−δ1−η2]

It can easily be shown that, Jacobian of the systemA2=[−r1*S0K1τ2BS000−τ2BS00000]A3=[0000η1ω1η2ω2000]

A4=[(1−S0K1)r1*−δ1−η300η3ω3ϵL0L0+K4+r2−δ2−τ1S000−α1L0(1−L0K3)r3−r3L0K3−δ3]

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.

Let u=(u1,u2,u3,u4,u5,u6)T be the right eigenvector and v=(v1,v3,v4,v5,v6)T be the left eigenvector linked with zero eigenvalues of the Jacobian matrix near r1=r1*, of the system.

Solving the system of the equations, we obtain;u1=0,u2=u3β−u5(η1ω1)Ω1+β+(S0k1−1)−δ1−η1,u3=−u5(η3ω3)(1−S0k1)r1−η3,u4=−u5(η3ω3)(1−S0k1)r1−η4

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

v1=−(η1−r1S0k1)v2+v3r1S0k1−r1S0k1v4+v5τ2BS0(1−S0k1)r1−r1S0k1−δ1=v3r1S0k1+v5τ2BS0(1−S0k1)r1−r1S0k1−δ1,v2=0

v3=((S0k1−1)r1+η1+δ1+Ω1+β)v2+v5τ2BS0(1−S0k1)r1−η3=v5τ2BS0(1−S0k1)r1−η3,v4=0,v5=−(η2ω2)v3−(η3ω3)v4εL0K4+L0+r2−δ2−τ1S0=−(η2ω2)v3εL0K4+L0+r2−δ2−τ1S0,v6=v5α1L((1−L0k1)r3−r3L0k3−δ3).

Let pk be the kth component of p anda=∑k,ij=1nvkuiuj∂2pk∂fi∂fj(0,0)andb=∑k,ij=1nvkui∂2pk∂fi∂r1*(0,0)

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.Fig. 2 A graph of variation of Immunity level against time.

Fig. 2

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.Fig. 3 A graph of variation of the rate of infection against time

Fig. 3

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. 4 A graph of variation of progression rate from pre-cancerous to cancerous against time.

Fig. 4

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).Fig. 5 A graph of variation of progression rate from pre-cancerous to cancerous against time

Fig. 5

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

Fig. 6

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

Fig. 7

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

Fig. 8

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

Fig. 9

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

Fig. 10

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

Fig. 11

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

Fig. 12

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

Fig. 13

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

Fig. 14

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

Fig. 15

3.3 Sensitivity analysis of the model

The relationship below provides the sensitivity index of the model parameter;SXRo=∂Ro∂X*XRo

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 R0*>1 and locally asymptotically stable when R0*<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 R0*<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, positivity 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 10.13039/100018051 University of Embu for granting me a full scholarship to pursue my master's degree.

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.

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 Routh TableLabel	
λ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,b1=−|1a2a1a3|a1=a1a2−a3a1b2=−|1a4a1a5|a1=a1a4−a5a1b3=−|1a6a10|a1=a6

c1=−|a1a3b1b2|b1=b1a3−a1b2b1c2=−|a1a5b1b3|b1=b1a5−a1b3b1d1=−|b1b2c1c2|c1=c1b2−b1c2c1

d2=−|b1b3c10|c1=b3e1=−|c1c2d1d2|d1=d1c2−c1d2d1f1=−|d1d2e10d1|=d2e1d1

Section 2 a1=1S0(L0+K4)(S0(ϵL0−L0r2−K4r2+L0δ2+K4δ2+(−L0−K4)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))−S0(L0+K4)((1−S0k1)r1−Ω4)+S0(L0+K4)(δ1+δ3+k3(1−L0k3)r3L0′+k1r1S0′−S0k12r1S0′))

a2=1S(L0+K4)(−AB(−L0−K4)η1ω1+S0((−L0−K4)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))−S0(L0ϵ−L0r2−K4r2+L0δ2+K4δ2+(−L0−K4)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))((1−S0k1)r1−Ω4)+(S0(L0ϵ−L0r2−K4r2+L0δ2+K4δ2+(−L0−K4)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))−S0(L0+K4)((1−S0k1)r1−Ω4))(δ1+δ3+k3(1−L0k3)r3L0′+k1r1S0′−S0k12r1S0′)+S0(L0+K4)(−δ3−k3(1−L0k3)r3L0′)(−δ1−k1r1S0′+S0k12r1S0′))

a3=1S0(L0+K4)(τ2B(β(L0+K4)η2ω2+(−L0−K4)η1ω1((1−S0k1)r1−Ω3))+S0(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(τ2B(−L0−K4)η1ω1−S0((−L0−K4)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3)))((1−S0k1)r1−Ω4)+(−AB(−L0−K4)η1ω1+S0((−L0−K4)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))−S0(L0ϵ−L0r2−K4r2+L0δ2+K4δ2+(−L0−K4)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))((1−S0k1)r1−Ω4))(δ1+δ3+k3(1−L0k3)r3L0′+k1r1S0′−S0k12r1S0′)+(S0(L0ϵ−L0r2−K4r2+L0δ2+K4δ2+(−L0−K4)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))−S0(L0+K4)((1−S0k1)r1−Ω4))(−δ3−k3(1−L0k3)r3L0′)(−δ1−k1r1S0′+S0k12r1S0′))

a4=1S0(L0+K4)((−AB(β(L0+K4)η2ω2+(−L0−K4)η1ω1((1−S0k1)r1−Ω3))−S0(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3))((1−S0k1)r1−Ω4)+(AB(β(L0+K4)η2ω2+(−L0−K4)η1ω1((1−S0k1)r1−Ω3))+S0(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(AB(−L0−K4)η1ω1−S0((−L0−K4)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3)))((1−S0k1)r1−Ω4))(δ1+δ3+k3(1−L0k3)r3L0′+k1r1S0′−S0k12r1S0)+(−AB(−L0−K4)η1ω1+S0((−L0−K4)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))−S0(L0ϵ−L0r2−K4r2+L0δ2+K4δ2+(−L0−K4)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3))((1−S0k1)r1−Ω4))(−δ3−k3(1−L0k3)r3L0′)(−δ1−k1r1S0′+S0k12r1S0′))

a5=1S0(L0+K4)((−τ2AB(β(L0+K4)η2ω2+(−L0−K4)η1ω1((1−S0k1)r1−Ω3))−S0(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3))((1−S0k1)r1−Ω4)(δ1+δ3+k3(1−L0k3)r3L0′+k1r1S0′−S0k12r1S0′)+(τ2B(β(L0+K4)η2ω2+(−L0−K4)η1ω1((1−S0k1)r1−Ω3))+S0(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(τ2B(−L0−K4)η1ω1−S0((−L0−K4)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3)+(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(r1−S0k1r1+(1−S0k1)r1−Ω2−Ω3)))((1−S0k1)r1−Ω4))(−δ3−k3(1−L0k3)r3LL0′)(−δ1−k1r1S0S0′+S0k12r1S0′))

a6=1S0(L0+K4)(−τ2B(β(L0+K4)η2ω2+(−L0−K4)η1ω1((1−S0k1)r1−Ω3))−S0(−L0ϵ+L0r2+K4r2−L0δ2−K4δ2)(−r1+S0k1r1+Ω2)((1−S0k1)r1−Ω3))((1−S0k1)r1−Ω4)(−δ3−k3(1−L0k3)r3L0′)(−δ1−k1r1S0′+S0k12r1S0′)

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>0ande1>0 are satisfied and otherwise unstable.

Section 3 The bifurcation results of a and b values are as follows;a=−2(−AB+u2)v3η2ω2(−ϵf6f6+K4+(1−f5K2)r2−δ2+(1−f5K2)η1ω1−f5r2+f2η1ω1+f3η2ω2+f4η3ω3K2)(βu3−u5[η1ω1])(ϵ1+k4+r2−δ2)(ε1+K4+r2−δ2)(−1+β+k1−δ1−η1+Ω1)+2(−AB+u2)v3η2ω2(−ϵf6f6+K4+(1−f5K2)r2−δ2+(1−f5K2)η2ω2−f5r2+f2η1ω1+f3η2ω2+f4η3ω3K2)u5[η3ω3](ϵ1+k4+r2−δ2)(ε1+K4+r2−δ2)((1−k1)r1−η3)+2(−AB+u2)v3η2ω2(−ϵf6f6+K4+(1−f5K2)r2−δ2+(1−f5K2)η3ω3−f5r2+f2η1ω1+f3η2ω2+f4η3ω3K2)u5[η3ω3](ϵ1+k4+r2−δ2)(ε1+K4+r2−δ2)((1−k1)r1−η4)−2ABv5(β+2ϴf33f4(D2+f32)2−2ϴf3f4D2+f32+1[1−NK1]r1−δ1−η2)(βu3−u5[η1ω1])u5[η3ω3]((1−k1)r1−η3)2(−1+β+k1−δ1−η1+Ω1)−2AB(β−ϴf32D2+f32)v5(βu3−u5[η1ω1])u5[η3ω3]((1−k1)r1−η3)+2v3η2ω2((1−f5K2)η1ω1+(1−f5K2)η2ω2)(βu3−u5[η1ω1])u5[η3ω3](ε1+K4+r2−δ2)((1−k1)r1−η3)(−1+β+k1−δ1−η1+Ω1)+2v3η2ω2((1−f5K2)η1ω1+(1−f5K2)η3ω3)(βu3−u5[η1ω1])u5[η3ω3](ε1+K4+r2−δ2)((1−k1)r1−η4)(−1+β+k1−δ1−η1+Ω1)+2ABv5(−ϴf32D2+f32+2ϴf33f4(D2+f32)2−2ϴf3f4D2+f32+1[1−NK1]r1−δ1−η2)u5[η3ω3]2((1−k1)r1−η3)2((1−k1)r1−η4)−2v3η2ω2((1−f5K2)η2ω2+(1−f5K2)η3ω3)u5[η3ω3]2(ε1+K4+r2−δ2)((1−k1)r1−η3)((1−k1)r1−η4)−2(−τ2B+u2)(−k1r1v4+τ2Bv5+v2(k1r1−η1))u5[η1ω1](ⅇMF+M+τ2Bf1(f1+Bf5)2−r1f1′[1−f1+f2+f3+f4K1]K1)((1−k1)r1−k1r1−δ1)}

b=2ϴf33f4(D2+f32)2−2ϴf3f4D2+f32−δ1−η2+(−AB+u2)(k1r1v3+ABv5)(ABf1(f1+Bf5)2+f1[1−f1+f2+f3+f4K1])((1−k1)r1−k1r1−δ1)+f3[1−f1+f2+f3+f4K1]+AB(−AB+u2)v5f3[1−f1+f2+f3+f4K1](ϵ1+k4+r2−δ2)((1−k1)r1−η3)+AB(β+f3(1−f1+f2+f3+f4K1)−f3r1K1)v5(βu3−u5[η1ω1])((1−k1)r1−η3)−(k1r1v3+ABv5)u5[η3ω3](f1[1−f1+f2+f3+f4K1]−r1f1′[1−f1+f2+f3+f4K1]K1)((1−k1)r1−k1r1−δ1)−(k1r1v3+ABv5)u5[η3ω3](f1[1−f1+f2+f3+f4K1]−r1f1′[1−f1+f2+f3+f4K1]K1)((1−k1)r1−k1r1−δ1)+(k1r1v3+ABv5)(βu3−u5[η1ω1])(ⅇMF+M+f1[1−f1+f2+f3+f4K1]−r1f1′[1−f1+f2+f3+f4K1]K1)((1−k1)r1−k1r1−δ1)−ABr1v5u5[η3ω3](1[1−f1+f2+f3+f4K1]−f3′[1−f1+f2+f3+f4K1]K1)((1−k1)r1−η3)2−ABv5u5[η3ω3](−ϴf32D2+f32+f3[1−f1+f2+f3+f4K1]−r1f3′[1−f1+f2+f3+f4K1]K1)((1−k1)r1−η3)

Section 4 b4=1;

b3=r1(−1MQF+M+β+δ1+η1−1δ1+η2−1δ1+η3+S(1MQF+M+β+δ1+η1+1δ1+η2+1δ1+η3)K1)−r2LϵL+K4+δ2+Sτ1

b2=r1((S−K1)2r1(Fβ+M(Q+β)+(F+M)(3δ1+η1+η2+η3))K12(Fβ+M(Q+β)+(F+M)(δ1+η1))(δ1+η2)(δ1+η3)+r2(1MQF+M+β+δ1+η1+1δ1+η2+1δ1+η3+S(−1MQF+M+β+δ1+η1−1δ1+η2−1δ1+η3)K1)LϵL+K4+δ2+Sτ1)

b1=−r13(MQF+M+β+δ1+η1)(δ1+η2)(δ1+η3)+S3r13K13(MQF+M+β+δ1+η1)(δ1+η2)(δ1+η3)−3S2r13K12(MQF+M+β+δ1+η1)(δ1+η2)(δ1+η3)+3Sr13K1(MQF+M+β+δ1+η1)(δ1+η2)(δ1+η3)−r12r2(MQF+M+β+δ1+η1)(δ1+η2)(LϵL+K4+δ2+Sτ1)−S2r12r2K12(MQF+M+β+δ1+η1)(δ1+η2)(LϵL+K4+δ2+Sτ1)+2Sr12r2K1(MQF+M+β+δ1+η1)(δ1+η2)(LϵL+K4+δ2+Sτ1)−r12r2(MQF+M+β+δ1+η1)(δ1+η3)(LϵL+K4+δ2+Sτ1)−S2r12r2K12(MQF+M+β+δ1+η1)(δ1+η3)(LϵL+K4+δ2+Sτ1)+2Sr12r2K1(MQF+M+β+δ1+η1)(δ1+η3)(LϵL+K4+δ2+Sτ1)−r12r2(δ1+η2)(δ1+η3)(LϵL+K4+δ2+Sτ1)−S2r12r2K12(δ1+η2)(δ1+η3)(LϵL+K4+δ2+Sτ1)+2Sr12r2K1(δ1+η2)(δ1+η3)(LϵL+K4+δ2+Sτ1).

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

Acknowledgement

The authors acknowledge the 10.13039/100018051 University of Embu for their kind support.
==== Refs
References

1 Chakraborty S. Li X.Z. Roy P.K. How can HPV-induced cervical cancer be controlled by a combination of drug therapy? A mathematical study Int. J. Biomath. (IJB) 12 6 2019 10.1142/S1793524519500700
2 Sierra-Rojas J.C. Reyes-Carreto R. Vargas-De-León C. Camacho J.F. 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 10.1155/2022/8650071
3 Insinga R.P. Dasbach E.J. Elbasha E.H. 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 10.1186/1471-2334-9-119 19144106
4 Alsaleh A.A. Gumel A.B. Dynamics analysis of a vaccination model for HPV transmission 22 4 2014 10.1142/S0218339014500211
5 Ribassin L. A SIS Model for Human Papillomavirus Transmission to Cite This Version : 2012
6 Gurmu E.D. Koya P.R. Impact of chemotherapy treatment on SITR compartmentalization and modeling of human papilloma virus 15 3 2019 17 29 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 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 Stark H. Zivković A. HPV vaccination: prevention of cervical cancer in Serbia and in europe Acta Fac. Med. Naissensis 35 1 2018 5 16 10.2478/afmnai-2018-0001
9 Kessler T.A. Cervical cancer: prevention and early detection Semin. Oncol. Nurs. 33 2 2017 172 183 10.1016/j.soncn.2017.02.005 28343836
10 Online P. Ezzati P.M. Health Policy NCD Countdown 2030 : Pathways to Achieving Sustainable Development Goal Target 3 . 4 vol. 396 2020 10.1016/S0140-6736(20)31761-X
11 Menzel R.W. Further notes on the albino catfish J. Hered. 49 6 1958 159 178 10.1093/oxfordjournals.jhered.a106828
12 Smith R.J. Li J. Blakeley D. The failure of R0 Comput. Math. Methods Med. 2011 2011 10.1155/2011/527610
13 Dadi Gurmu E. Rao Koya P. 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 10.11648/j.ajam.20190703.11
14 Akimenko V.V. Adi-Kusumo F. Stability analysis of an age-structured model of cervical cancer cells and HPV dynamics Math. Biosci. Eng. 18 5 2021 6155 6177 10.3934/mbe.2021308 34517528
15 Chakraborty S. Cao X. Bhattyacharya S. Roy P.K. The role of HPV on cervical cancer with several functional response: a control based comparative study Comput. Math. Model. 30 4 2019 439 453 10.1007/s10598-019-09469-4
16 Asih T.S.N. The dynamics of HPV infection and cervical cancer cells Bull. Math. Biol. 78 1 2016 4 20 10.1007/s11538-015-0124-2 26676766
17 Ndii M.Z. Mathematical Model of Cervical Cancer Treatment Using Chemotherapy Drug February 2020 10 16 10.14421/biomedich.2019.81.11-15
18 Pongsumpun P. Mathematical Model of Cervical Cancer Due to Human Papillomavirus Infection 2012 157 161 no. January
19 Sulaiman U. Adamu I.I. Tahir A. 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 Allali K. Stability analysis and optimal control of HPV infection model with early-stage cervical cancer Biosystems 199 December 2020 2021 104321 10.1016/j.biosystems.2020.104321
21 Ma B. Thibodeaux J.J. Well-Posedness and a Finite Difference Approximation for a Mathematical Model of HPV-Induced Cervical Cancer 2023 151 172 10.4236/am.2023.143009
22 Mondaini R.P. Lectures C. De Janeiro R. Trends in Biomathematics: Chaos and Control in Epidemics, Ecosystems, and Cells 2021 10.1007/978-3-030-73241-7
23 Erwin S.H. 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
24 Verma M. Modeling the mechanisms by which HIVAssociated immunosuppression influences HPV persistence at the Oral Mucosa PLoS One 12 1 2017 1 20 10.1371/journal.pone.0168133
25 Belew B. Melese D. Modeling and analysis of predator-prey model with fear effect in prey and hunting cooperation among predators and harvesting J. Appl. Math. 2022 2022 10.1155/2022/2776698
26 Mutua G.K. Ngari C.G. Muthuri G.G. Kitavi D.M. Mathematical modeling and simulating of Helicobacter pylori treatment and transmission implications on stomach cancer dynamics Commun. Math. Biol. Neurosci. 2022 2022 1 29 10.28919/cmbn/7542
27 Diekmann O. Heesterbeek J.A.P. Metz J.A.J. 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 10.1007/BF00178324 2117040
28 Shaw C.L. Kennedy D.A. What the reproductive number R0 can and cannot tell us about COVID-19 dynamics Theor. Popul. Biol. 137 2021 2 9 10.1016/j.tpb.2020.12.003 33417839
29 Kilonzi J.S. Ngari C.G. Karanja S. Modelling the Impact of Screening , Treatment and Underlying Health Conditions on Dynamics of Covid-19 January 2023, 2024
30 Ochwach J. Okongo M.O. Ngari C. Mathematical Modelling of Cholera Incorporating the Dynamics of the Induced Achlorhydria Condition and Treatment November, 2022
31 Gitonga N.C. Modelling Childhood Pneumonia and its Implications Regarding its Control Using Kenyan Data 2017 1 120
