Analysis of a mathematical model of human papillomavirus transmission dynamics with optimal control for cervical cancer prevention

This study presents a mathematical model that explains the transmission dynamics of high-risk Human Papillomavirus (HPV) and how it progresses to cervical cancer. The model particularly includes a population group for Cervical Intraepithelial Neoplasia (CIN), which is a critical precancerous stage that is usually neglected in previous models, and brings in a corresponding treatment compartment to fit realistic clinical interventions. Control measures such as vaccination, public campaign awareness, and regular screening are merged into the framework to analyze their unified impact on disease spread and cancer prevention. Using a system of non-linear differential equations, the model is assessed for positivity, boundedness, and stability. The basic reproduction number $R_o$ is computed using the next-generation matrix approach, and both local and global sensitivity analyses are carried out to determine the parameters that are most influential. An optimal control structure that is based on Pontryagin's Maximum Principle is applied to evaluate suitable intervention strategies. Simulation outcome reveals that increasing awareness and screening coverage can remarkably decrease both HPV infection and cervical cancer cases, especially in settings with inadequate vaccine coverage. This model provides practical insights for public health decision-making and strengthens the importance of integrated intervention in HPV-induced cervical cancer control.

1 Introduction

Human papillomavirus (HPV) is a sexually transmitted infection that is primarily spread through intimate skin-to-skin contact. Unlike many other sexually transmitted infections, the transmission of HPV does not require the exchange of bodily fluids [1]. This implies that HPV cannot be transmitted through blood or semen. However, this in no way lessens the virus's ability to spread, as skin-to-skin contact—whether penetrative or non-penetrative, including genital-genital, anal-genital, oral-genital, and oral-anal contact—can easily result in transmission from an infected individual [2]. HPV has more than 150 strains, and it is possible for an individual to be infected with multiple HPV strains simultaneously [3, 4]. According to the International Agency for Research on Cancer, thirteen HPV strains are capable of causing cervical cancer, and HPV-16 and HPV-18—known as high-risk types—account for approximately 70% of cervical cancer cases worldwide [5, 6].

In 2019 alone, HPV gave rise to an estimated 620,000 cancer cases in women and 70,000 cancer cases in men [7]. The World Health Organization stated that prophylactic vaccination against HPV can prevent these cancers. They also added that HPV screening and the treatment of precancerous lesions are effective ways to prevent cervical cancer [7]. In light of this, the model formulated in this study incorporates a screening parameter that facilitates the treatment of precancerous lesions to assess its impact on cervical cancer cases. The formulated model also integrates an awareness campaign parameter that promotes vaccination uptake.

Over 90% of sexually active men and over 80% of sexually active women will be infected with the virus in their lifetime [8]. Even though vaccines are available to prevent this infection, it should be noted that the vaccines are not effective for those who already have the virus. Hence, it is important to examine suitable control strategies to manage existing infections and prevent their progression to cervical cancer, which is the most common HPV-associated cancer [9].

Before the development of cervical cancer in women, it is usually preceded by cervical dysplasia, which is characterized by the presence of precancerous lesions in the cervix. If left untreated, this precancerous condition may develop into cervical cancer after ten years or more [10]. Additionally, it takes about five to ten years for cervical cells infected with HPV to develop into precancerous lesions [11].

Mathematical models have been instrumental in describing the transmission dynamics of various viral diseases. Since studies have shown that 80% of sexually active women will be infected with the virus in their lifetime, there is a need to develop a mathematical model that helps to understand the transmission dynamics of high-risk HPV and its impact on cervical cancer cases. Currently, Africa has the highest incidence of cervical cancer [12]. Records show that about 660,000 cervical cancer cases and 350,000 related deaths were reported globally in 2022 [13], and approximately 125,700 African women were diagnosed with cervical cancer in the same year [14]. Estimates indicate that more than 266,000 women die from cervical cancer annually, with 87% of these deaths occurring in developing countries, particularly within the WHO African region. This further underscores the need to optimize control strategies to identify effective solutions that can minimize cervical cancer incidence in Africa and globally.

In reviewing models on HPV-related cervical cancer, it was observed that, over time, these models have improved, with each new model building upon the previous one to align more effectively with biological and clinical evidence.

“Early models like the one developed by Lee and Tameru [15] came up with an idea of incorporating the progression of HPV transmission to Cervical cancer. They also showed the potential effect of vaccination. Their result suggested that HPV can be curbed if prevention and treatment are intensified [15]. Zhang et al. took an important step forward by developing an ODE-based model that combined both screening and vaccination in its formulation [16]. The result of their study shows that increasing the vaccination rate of the system is the most effective way to reduce the basic reproduction number. Despite this, however, a critical limitation persisted in their model- a separate class for treating individuals with precancerous conditions (CIN) that screening exercises detect was not presented in the model. The absence of this compartment happens to be frequent in subsequent models that incorporate screening in their formulation. Other Authors introduced a more realistic characterization in their formulation. For example, Ziyadi et al. utilized a two-sex framework [17], while Ren et al. developed an advanced age-structured model that is in agreement with WHO's elimination approach [18], and Omame et al. studied the complexities of HPV co-infection with another disease-causing organism, specifically Chlamydia trachomatis [19]. Lankoande et al. proposed a two-sex model that explains HPV transmission and its progression to cervical cancer [20]. The result of their research stresses the importance of addressing the transmission pathways of same-sex and vertical transmission in the public health sector in order to effectively reduce HPV spread and mortalities driven by Cervical Cancer [20].”

An important development in this research area has been the implementation of optimal control theory to identify effective public health intervention strategies. Studies conducted by authors such as Arunkumar et al. and Saldana et al. applied this technique to analyze the efficacy of intervention strategies, focusing specifically on awareness and vaccination programs, respectively [21, 22]. This present study uniquely incorporates an optimal control approach that simultaneously achieves an effective balance between awareness campaigns and screening programs, which promote vaccination uptake and treatment, respectively, since existing studies on optimal control have not simultaneously optimized the combined impact of awareness programs that promote vaccination and screening efforts, which drive the treatment of CIN.

Regardless of the developments so far in existing models, none have incorporated an explicit compartment dedicated to treating individuals with cervical dysplasia (T_Q), whose increase is driven by a screening parameter. This omission overlooks a crucial aspect of cervical cancer intervention recommended by WHO [7].

The Table 1 provides a summary of the focus of existing studies and the corresponding gap relative to this current work.

Table 1

Table 1. Comparison of existing HPV models with current study.

In order to address these shortcomings, this study formulates an innovative compartmental model that distinctively incorporates a compartment dedicated to treating women with precancerous conditions. The entry rate at which individuals enter this treatment compartment is driven by a specific parameter that represents screening programs. It also integrates an awareness parameter that drives vaccination uptake. Additionally, an optimal control problem is developed and solved to determine the approach that is most effective between awareness campaigns and screening programs in terms of resource allocation. In light of this, it is safe to say that this research is dedicated to providing a tool for public health that is more clinically realistic and practically useful to assist their decision-making aimed at combating HPV infection and HPV-driven Cervical Cancer in the long run.

This study is organized as follows: Section 2 presents the model assumptions and formulation, displaying the compartmental structure and flow of the HPV-Cervical Cancer model. In Section 3, the model is analyzed mathematically by establishing the positivity, boundedness, and invariant region of the system. The basic reproduction number is derived using the next-generation matrix approach, and the local stability of the disease-free and endemic equilibrium points is analyzed. Section 4 details the estimation of initial conditions of the dependent variable and model parameters based on the epidemiological data that are relevant to the Nigerian context. Section 5 performs both local and global sensitivity analyses to identify the most influential parameters of the disease transmission. In Section 6, optimal control strategies are formulated and analyzed using Pontryagin's Maximum Principle to assess the impact of awareness, vaccination, and screening interventions. Section 7 presents the numerical simulation and results, and finally, the study concludes with key findings and policy recommendations in Section 8.

2 Model assumptions and formulation

This model is formulated based on the recent control strategies available to curb the spread of HPV and reduce cervical cancer incidence.

2.1 Assumptions

The formulation of this model is based on the following assumptions:

• Total population considered is not constant;

• All compartments experience a similar natural death rate μ;

• Susceptible women and men are 9 years and older, and any increase in population is due to aging;

• Gardasil 9 vaccine is considered, and it is assumed that all who move to V compartment have taken the required three-dose series, and re-infection is not possible;

• Infectious classes are made up of individuals who have only HPV 16 who are 15 years and above;

• Women in the T_Q compartment still undergo a routine screening called the pap smear test to detect any relapse;

• Q compartment houses women with CIN 1, CIN 2, and CIN3;

• Infectious men have persistent HPV-16 infection.

• The compartment C_c comprises women with early-stage cervical cancer, and thus, there is no mortality rate yet.

• Infections with HPV occur through interactions between susceptible and infectious populations.

Justifications of some assumptions: Studies reveal that the HPV strain most likely to persist compared to other types is HPV-16 [23]. The male population is assumed to have persistent HPV infection in this study because studies regarding human papillomavirus in men show that when they are infected with one strain of HPV, their risk of reinfection with the same type of HPV increases substantially. In fact, it has been discovered that men who have been infected with high-risk HPV are twenty times more likely to be reinfected within 1 year [24]. Moreover, studies reveal that men have a lower seroconversion rate than females of the same age, which indicates that males have a greater tendency to be reinfected throughout life [25].

A heterosexual infection transmission is assumed in this study since, globally, 80% of the population identifies as heterosexual [26]. Particularly, in a concentrated study at a Nigerian University, 96% of students identify as heterosexual [27]. Also, over 90% identified as heterosexual in another research conducted amongst Nigerian students in another university [28]. Furthermore, epidemiological studies have indicated that heterosexual transmission is the predominant mode of HPV spread in a lot of populations where same-sex relationships are under-reported and not prevalent [29, 30].

All parameters are positive in the model. The force of infection utilized is mass action, and the total population N is given as:

$\begin{array}{rr}\begin{array}{r}N(t)=S_F(t)+I_F(t)+R_F(t)+Q(t)+C_C(t)+T_Q(t)+S_M(t)\\ +I_M(t)\end{array}& (1)\end{array}$

The Table 2 explains the different compartments used in this model.

Table 2

Table 2. Dependent variables and their physical meanings.

2.2 Model formulation

This model is an SIR based extended compartmental model that individually accounts for the male and female population while still incorporating control measures to reduce the occurrence of cervical cancer in women.

2.2.1 The schematic diagram

The flow diagram of the model is provided in Figure 1.

Figure 1

Figure 1. Schematic diagram of HPV-cervical cancer model.

2.2.2 The system of equations

$\begin{array}{lll}\frac{d{S}_F}{dt}& =& {\text{Λ}}_1-\mu S_F-(1-\epsilon ){\beta }_1{S}_F{I}_M-\epsilon \tau S_F\\ \frac{d{I}_F}{dt}& =& (1-\epsilon ){\beta }_1{S}_F{I}_M-(\mu +\sigma +c)I_F\\ \frac{dQ}{dt}& =& c{I}_F-(\alpha \rho +\varphi +\mu )Q-(1-\alpha )\rho Q\\ \frac{d{T}_Q}{dt}& =& \alpha \rho Q-\mu T_Q\\ \frac{d{C}_C}{dt}& =& (1-\alpha )\rho Q-\mu C_C\\ \frac{d{R}_F}{dt}& =& \sigma I_F-\mu R_F+\varphi Q\\ \frac{d{S}_M}{dt}& =& {\text{Λ}}_2-{\beta }_2{S}_M{I}_F-\mu S_M-d{S}_M\\ \frac{d{I}_M}{dt}& =& {\beta }_2{S}_M{I}_F-\mu I_M\\ \frac{dV}{dt}& =& -\mu V+\epsilon \tau S_F+d{S}_M& (2)\end{array}$

Table 3 explains the physical meanings of the parameters used in the model formulation.

Table 3

Table 3. Physical interpretation of the model parameters.

3 Model analysis

The analysis of the model comprises the following in this section: The positivity and the boundedness of solutions; The disease-free and the endemic steady states; Stability of steady states and the basic reproduction number of the model.

3.1 Positivity of model solutions and boundedness

3.1.1 Positivity of model solutions

Here, we show that the system of Equation 2 is well defined and biologically meaningful in the realistic sense since it represents a population system. It is important that we do not obtain negative values of the population classes, and thus, we must prove the positivity of the solutions.

Let the initial conditions of the sub-populations be:

S_F(0)≥0, I_F(0)≥0, Q(0)≥0, T_Q(0)≥0, C_C(0)≥0, V(0)≥0, S_M(0)≥0, I_M(0)≥0 and R_F(0)≥0

Positivity of the solutions of the model implies that the state variables remain non-negative for all values of t, provided the initial conditions are non-negative [16].

Mathematical Proof: From the model (Equation 2), for any compartment, say X, it will be observed that whenever X = 0, the derivative $\frac{dX}{dt}\ge 0$

This hereby ensures that all the different state variables remain non-negative as they change with time.

3.1.2 Boundedness of solutions

It is important that we show that the total population N(t) does not grow unrealistically. This is because in the realistic sense, natural deaths and other factors will not allow the total population to grow infinitely.

Recall the total population N(t) given as:

$\begin{array}{r}N(t)=S_F(t)+I_F(t)+R_F(t)+Q(t)+C_C(t)+T_Q(t)+S_M(t)\\ +I_M(t)+V(t)\end{array}$

The total population with respect to time is:

$\begin{array}{r}\frac{dN}{dt}=\frac{d{S}_F}{dt}+\frac{d{I}_F}{dt}+\frac{d{R}_F}{dt}+\frac{dQ}{dt}+\frac{d{C}_C}{dt}+\frac{d{T}_Q}{dt}+\frac{d{S}_M}{dt}+\\ \frac{d{I}_M}{dt}+\frac{dV}{dt}\end{array}$

The sum of all the equations in the model (Equation 2) gives:

$\begin{array}{ll}\frac{dN}{dt}=\text{Λ}-\mu N& (3)\end{array}$

Where Λ = Λ₁+Λ₂

Solving (Equation 3) we have:

$\begin{array}{ll}N(t)=\frac{\text{Λ}}{\mu }+C{e}^{-\mu t}& (4)\end{array}$

As time t tends to ∞, the second term on the RHS, Ce^−μt, approaches zero.

This implies that:

$\begin{array}{ll}{\displaystyle \underset{t\to \infty }{lim sup}}N(t)\le \frac{\text{Λ}}{\mu }& (5)\end{array}$

The physical interpretation is that as time t tends to ∞, the population will not grow beyond the fixed limit, $\frac{\text{Λ}}{\mu }$ and thus N(t) will not explode to infinity since $\frac{\text{Λ}}{\mu }$ is the largest value it can take as time t approaches ∞.

3.1.3 Invariant region

This concept ensures that the model solutions continue to maintain a realistic set of values for all time t≥0.

If we define a region K:

$\begin{array}{r}K=\{(S_F(t),I_F(t),R_F(t),Q(t),C_C(t),T_Q(t),S_M(t),I_M(t),V(t))\\ \in {ℝ}_{+}^9|S_F(t)+I_F(t)+R_F(t)+Q(t)+C_C(t)+T_Q(t)+\\ S_M(t)+I_M(t)+V(t)\le \frac{\text{Λ}}{\mu }\}\end{array}$

Where Λ = Λ₁+Λ₂.

With Equation 5, it can be concluded that N(t) = S_F(t)+I_F(t)+R_F(t)+Q(t)+C_C(t)+T_Q(t)+S_M(t)+I_M(t)+V(t)remains bounded which implies that model (Equation 2) is positively invariant in the defined region K.

This implies that if the system starts in the region K, it remains in the region K. The physical interpretation of this is that, at any time t, the total population will not exceed $\frac{\text{Λ}}{\mu }$, making sure that all the solutions of the model stay realistic or biologically meaningful.

3.2 Basic reproduction number

Basic reproduction number $R_O$ is the average number of secondary infections that will occur when one infectious individual is introduced into an environment filled with only susceptible individuals [31]. Jones [31] described a method for calculating R_O which involves the formation of the next generation matrix, FV⁻¹. The spectral radius, ρ(FV⁻¹) represents the basic reproduction number.

$\begin{array}{ll}F={\left[\frac{\partial f_i}{\partial x_j}\right]}_{x=x_0},V={\left[\frac{\partial v_i}{\partial x_j}\right]}_{x=x_0}& (6)\end{array}$

Physically,

• F measures how new infections arise owing to contact with infected individuals.

• V accounts for disease progression, recovery, and natural or disease-induced deaths within infected compartments.

and:

• f_i = rate of new infections in the i^th infected compartment,

• v_i = net rate of transfer into and out of the i^th infected compartment (with new infections excluded),

• x_j = state variable of the j^th infected compartment (I_F, Q, C_C, T_Q, and I_M ),

• x₀ = Disease-Free Equilibrium (DFE) point,

• $\frac{\partial f_i}{\partial x_j}$ = partial derivative of f_i with respect to x_j, evaluated at the DFE,

• $\frac{\partial v_i}{\partial x_j}$ = partial derivative of v_i with respect to x_j, evaluated at the DFE.

$R_O$ is usually computed at the disease free equilibrium point- DFE (also known as the infected free equilibrium points [32]). The disease-free equilibrium point x₀ of the model (Equation 2) refers to that point when there is no disease in the system. Thus, at DFE $S_F^{(0)}\ne 0,V^{(0)}\ne 0,S_M^{(0)}\ne 0$ and $I_F^{(0)}=Q^{(0)}=T_Q^{(0)}=C_C^{(0)}=I_M^{(0)}=R_F^{(0)}=0$.

$\begin{array}{rr}\begin{array}{r}{x}_0=(S_F{(t)}^{(0)},I_F{(t)}^{(0)},R_F{(t)}^{(0)},Q{(t)}^{(0)},C_C{(t)}^{(0)}{T}_Q{(t)}^{(0)},S_M{(t)}^{(0)},\\ I_M{(t)}^{(0)},V{(t)}^{(0)})\end{array}& (7)\end{array}$

The equilibrium points are obtained when $\frac{d{S}_F}{dt}=\frac{d{I}_F}{dt}=\frac{dQ}{dt}=\frac{d{T}_Q}{dt}=\frac{d{C}_C}{dt}=\frac{dV}{dt}=\frac{d{S}_M}{dt}=\frac{d{I}_M}{dt}=\frac{d{R}_F}{dt}=0$. Thus, from the system (Equation 2), we have:

$\begin{array}{ll}\begin{array}{l}{x}_0=(\frac{{\text{Λ}}_1}{\mu +\epsilon \tau },0,0,0,0,0,\frac{{\text{Λ}}_2}{\mu +d},0,\\ \frac{1}{\mu }[\frac{\epsilon \tau {\text{Λ}}_1(\mu +d)+d{\text{Λ}}_2(\mu +\epsilon \tau )}{(\mu +\epsilon \tau )(\mu +d)}])\end{array}& (8)\end{array}$

The DFE x₀ is non-negative since:

$S_F(0)=\frac{{\text{Λ}}_1}{\mu +\epsilon \tau }>0,$

$S_M(0)=\frac{{\text{Λ}}_2}{\mu +d}>0,$

$V(0)=\frac{\epsilon \tau {\text{Λ}}_1(\mu +d)+d{\text{Λ}}_2(\mu +\epsilon \tau )}{\mu (\mu +\epsilon \tau )(\mu +d)}>0,$

and all other compartments are zero. This condition stands because all the parameters involved: Λ₁, Λ₂, μ, ε, τ, , and d are positive, to ensure that the DFE is meaningful, biologically.

To compute $R_O$, the compartments I_F, Q, T_Q, C_C and I_M from the model (Equation 2) are made use of.

$f=\left[\begin{array}{c}(1-\epsilon ){\beta }_1{S}_F{I}_M\\ 0\\ 0\\ 0\\ {\beta }_2{S}_M{I}_F\end{array}\right]$

$\begin{array}{l}-\text{v}=\\ \left[\begin{array}{c}-(\mu +\sigma +c)I_F\\ ϵI_F-(\alpha \rho +\varphi +\mu )-(1-\alpha )\rho Q\\ \alpha \rho Q-\mu T_Q\\ (1-\alpha )\rho Q-\mu C_C\\ -\mu I_M\end{array}\right]\end{array}$

$\begin{array}{rr}\begin{array}{r}F=\left[\begin{array}{rrrrr}0& 0& 0& 0& (1-\epsilon ){\beta }_1{S}_F\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ {\beta }_2{S}_M& 0& 0& 0& 0\end{array}\right]\\ -V=\\ \left[\begin{array}{ccccc}-(\mu +\sigma +c)& 0& 0& 0& 0\\ \epsilon & -\left[(\alpha \rho +\varphi +\mu )+(1-\alpha )\rho \right]& 0& 0& 0\\ 0& \alpha \rho & -\mu & 0& 0\\ 0& (1-\alpha )\rho & 0& -\mu & 0\\ 0& 0& 0& 0& -\mu \end{array}\right]\end{array}& (9)\end{array}$

$\begin{array}{ll}F{V}^{-1}=\left[\begin{array}{ccccc}0& 0& 0& 0& \frac{(1-\epsilon ){\beta }_1{\text{Λ}}_1}{\mu (\mu +\epsilon \tau )}\\ \epsilon & 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ \frac{{\beta }_2{\text{Λ}}_2}{(\mu +d)(\mu +m\sigma +c)}& 0& 0& 0& 0\end{array}\right]& (10)\end{array}$

The basic reproduction number is obtained as:

$\begin{array}{ll}\rho (F{V}^{-1})=R_O=\sqrt{\frac{{\beta }_1{\beta }_2{\text{Λ}}_1{\text{Λ}}_2(1-\epsilon )}{\mu (\mu +\sigma +c)(\mu +\epsilon \tau )(\mu +d)}}& (11)\end{array}$

Also, we can write that:

$\begin{array}{ll}{R}_{O1}=\sqrt{\frac{{\beta }_1{\text{Λ}}_1(1-\epsilon )}{(\mu +\sigma +c)(\mu +\epsilon \tau )}}& (12)\end{array}$

Where $R_{O1}$ represents the basic reproduction number of the female population.

$\begin{array}{ll}{R}_{O2}=\sqrt{\frac{{\beta }_2{\text{Λ}}_2}{\mu (\mu +d)}}& (13)\end{array}$

Where $R_{O2}$ stands for the basic reproduction number of the male population.

When the basic reproduction number is less than 1, it implies that the disease will fizzle out with time, but will spread when it is greater than 1 [33].

3.3 Endemic equilibrium point EEP

The endemic equilibrium points x_* is obtained when the infected population is not zero [34, 35]. Thus,

$\begin{array}{r}\begin{array}{r}{S}_F^{(*)}\ne 0,I_F^{(*)}\ne 0,Q^{(*)}\ne 0,T_Q^{(*)}\ne 0,C_C^{(*)}\ne 0,V^{(*)}\ne 0,S_M^{(*)}\ne 0,\\ I_M^{(*)}\ne 0,R_F^{(*)}\ne 0.\end{array}\end{array}$

As previously stated, the equilibrium points are obtained when $\frac{d{S}_F}{dt}=\frac{d{I}_F}{dt}=\frac{dQ}{dt}=\frac{d{T}_Q}{dt}=\frac{d{C}_C}{dt}=\frac{dV}{dt}=\frac{d{S}_M}{dt}=\frac{d{I}_M}{dt}=\frac{d{R}_F}{dt}=0$.

Working rigorously with the S_F, S_M, I_F, and I_M equations from the system of Equation 2, I_M is presented in terms of parameters only as written in Equation 14 below. The other endemic equilibrium points are also presented in terms of I_M as well.

$\begin{array}{lll}{I}_M^{(*)}& =\frac{{\beta }_2{\text{Λ}}_2(1-\epsilon ){\beta }_1{\text{Λ}}_1-\mu (\mu +d)(\mu +\sigma +c)(\mu +\epsilon \tau )}{\mu (\mu +d)(\mu +\sigma +c)(1-\epsilon ){\beta }_1+\mu {\beta }_2(1-\epsilon ){\beta }_1{\text{Λ}}_1}& (14)\end{array}$

$\begin{array}{lll}{S}_F^{(*)}& =\frac{{\text{Λ}}_1}{\mu +(1-\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau }& (15)\end{array}$

$\begin{array}{lll}{I}_F^{(*)}& =\frac{(1-\epsilon ){\text{Λ}}_1{\beta }_1{I}_M^{(*)}}{(\mu +\sigma +c)(\mu +(1-\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau )}& (16)\end{array}$

$\begin{array}{lll}{Q}^{(*)}& =\frac{(1-\epsilon ){\text{Λ}}_1{\beta }_1{I}_M^{(*)}\epsilon }{\begin{array}{c}\left[(\alpha \rho +\mu )+(1-\alpha )\rho \right][(\mu +\sigma +c)(\mu +(1-\epsilon ){\beta }_1{I}_M^{(*)}\\ +\epsilon \tau )]\end{array}}& (17)\end{array}$

$\begin{array}{lll}{T}_Q^{(*)}& =\frac{(1-\epsilon ){\text{Λ}}_1{\beta }_1{I}_M^{(*)}c\alpha \rho }{\begin{array}{c}\mu \left[(\alpha \rho +\mu )+(1-\alpha )\rho \right][(\mu +\sigma +c)(\mu +(1\\ -\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau )]\end{array}}& (18)\end{array}$

$\begin{array}{lll}{C}_C^{(*)}& =\frac{(1-\alpha )(1-\epsilon ){\text{Λ}}_1{\beta }_1{I}_M^{(*)}c\rho }{\begin{array}{c}\mu \left[(\alpha \rho +\mu )+(1-\alpha )\rho \right][(\mu +\sigma +c)(\mu \\ +(1-\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau )]\end{array}}& (19)\end{array}$

$\begin{array}{lll}{S}_M^{(*)}& =\frac{\begin{array}{c}{\beta }_2(1-\epsilon ){\text{Λ}}_1{\beta }_1{I}_M^{(*)}+(\mu +d)(\mu +\sigma +c)(\mu +(1\\ -\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau )\end{array}}{{\text{Λ}}_2(\mu +\sigma +c)(\mu +(1-\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau )}& (20)\end{array}$

$\begin{array}{ll}{V}^{(*)}=\frac{\begin{array}{c}(\mu +(1-\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau )[\epsilon \tau {\text{Λ}}_1{\text{Λ}}_2(\mu +\sigma +c)\\ +d({\beta }_2(1-\epsilon ){\text{Λ}}_1{\beta }_1{I}_M^{(*)}+(\mu +d)(\mu +\sigma +c))]\end{array}}{\begin{array}{c}\mu (\mu +(1-\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau ){\text{Λ}}_2(\mu +\sigma +c)(\mu \\ +(1-\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau )\end{array}}& (21)\end{array}$

$\begin{array}{l}{R}_F^{(*)}=\frac{\begin{array}{c}[(\alpha \rho +\mu )+(1-\alpha )\rho ]·(\mu +\sigma +c)(\mu +(1-\epsilon ){\beta }_1{I}_M^{(*)}\\ +\epsilon \tau )\end{array}}{(\mu +\sigma +c)(\mu +(1-\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau )}\\ \times \frac{\left[\sigma (1-\epsilon ){\text{Λ}}_1{\beta }_1{I}_M^{(*)}+\varphi (1-\epsilon ){\text{Λ}}_1{\beta }_1\epsilon I_M^{(*)}\right]}{(\alpha \rho +\varphi +\mu )+(1-\alpha )\rho [(\mu +\sigma +\epsilon ){\beta }_1{I}_M^{(*)}+\epsilon \tau ]}& (22)\end{array}$

Rewriting Equation 14, we have:

$\begin{array}{l}{I}_M^{(*)}=\frac{\mu (\mu +d)(\mu +\sigma +c)(\mu +\epsilon \tau )[\frac{{\beta }_2{\beta }_1{\text{Λ}}_1{\text{Λ}}_2(1-\epsilon )}{\mu (\mu +d)(\mu +\sigma +c)(\mu +\epsilon \tau )}-1]}{{\beta }_1(1-\epsilon )\left[\mu (\mu +d)(\mu +\sigma +c)+\mu {\beta }_2{\text{Λ}}_1\right]}\end{array}$

I_M(*) can then be written in terms of $R_O$ obtained in Equation 11.

$\begin{array}{l}{I}_M(*)=\mu (\mu +d)(\mu +\sigma +c)(\mu +\epsilon \tau )\\ \frac{\left[R_0^2-1\right]}{{\beta }_1(1-\epsilon )\left[\mu (\mu +d)(\mu +\sigma +c)+\mu {\beta }_2{\text{Λ}}_1\right]}.\end{array}$

I_F(*) can also be written in terms of $R_0$ as below:

$\begin{array}{l}{I}_F(*)=\frac{(1-\epsilon ){\text{Λ}}_1{\beta }_1{c}_3(R_0^2-1)}{(\mu +\epsilon \tau )(\mu +\sigma +c)+(1-\epsilon )(\mu +\sigma +c){\beta }_1{c}_3(R_0^2-1)},\end{array}$

Where $c_3=\frac{\mu (\mu +d)(\mu +\sigma +c)(\mu +\epsilon \tau )}{{\beta }_1(1-\epsilon )\left[\mu (\mu +d)(\mu +\sigma +c)+\mu {\beta }_2{\text{Λ}}_1\right]}$.

The two equations above show I_F and I_M written in terms of the basic reproduction number (11); It can be observed that if $R_0>1$, then $I_M^{*}>0$ and $I_F^{*}>0$. This implies that an endemic state exists (that is, the disease persists) when the basic reproduction number is greater than 1, the DFE becomes unstable, and infection grows [34, 36].

However, if $R_0<1$, then $I_M^{(*)}<0$ and $I_F^{(*)}<0$. This is non-physical, and it indicates the absence of any biologically meaningful positive endemic solution.

3.4 Local stability of disease free equilibrium DFE

The Jacobian of the model (Equation 2) is obtained as:

$J=\left[\begin{array}{ccccccccc}-\mu -(1-\epsilon ){\beta }_1{I}_M-\epsilon \tau & 0& 0& 0& 0& 0& 0& -(1-\epsilon ){\beta }_1{S}_F& 0\\ (1-\epsilon ){\beta }_1{I}_M& -(\mu +\sigma +c)& 0& 0& 0& 0& 0& (1-\epsilon ){\beta }_1{S}_F& 0\\ 0& \epsilon & -(\alpha \rho +\varphi +\mu )-(1-\alpha )\rho & 0& 0& 0& 0& 0& 0\\ 0& 0& \alpha \rho & -\mu & 0& 0& 0& 0& 0\\ 0& 0& (1-\alpha )\rho & 0& -\mu & 0& 0& 0& 0\\ \epsilon \tau & 0& 0& 0& 0& -\mu & d& 0& 0\\ 0& -{\beta }_2{S}_M& 0& 0& 0& 0& -{\beta }_2{I}_F-\mu -d& 0& 0\\ 0& {\beta }_2{S}_M& 0& 0& 0& 0& {\beta }_2{I}_F& -\mu & 0\\ 0& \sigma & \varphi & 0& 0& 0& 0& 0& -\mu \end{array}\right]$

The Jacobian at DFE J(x₀) is written as:

$J(x_0)=\left[\begin{array}{ccccccccc}-\mu -\epsilon \tau & 0& 0& 0& 0& 0& 0& \frac{-(1-\epsilon ){\beta }_1{\text{Λ}}_1}{\mu +\epsilon \tau }& 0\\ 0& -(\mu +\sigma +c)& 0& 0& 0& 0& 0& \frac{(1-\epsilon ){\beta }_1{\text{Λ}}_1}{\mu +\epsilon \tau }& 0\\ 0& \epsilon & -(\alpha \rho +\varphi +\mu )-(1-\alpha )\rho & 0& 0& 0& 0& 0& 0\\ 0& 0& \alpha \rho & -\mu & 0& 0& 0& 0& 0\\ 0& 0& (1-\alpha )\rho & 0& -\mu & 0& 0& 0& 0\\ \epsilon \tau & 0& 0& 0& 0& -\mu & d& 0& 0\\ 0& \frac{-{\beta }_2{\text{Λ}}_2}{\mu +d}& 0& 0& 0& 0& -\mu -d& 0& 0\\ 0& \frac{{\beta }_2{\text{Λ}}_2}{\mu +d}& 0& 0& 0& 0& {\beta }_2{I}_F& -\mu & 0\\ 0& \sigma & \varphi & 0& 0& 0& 0& 0& -\mu \end{array}\right]$

The determinant of |J(x₀)−λI| is:

$|J(x_0)-\lambda I|=-{(\mu -\lambda )}^4(-\mu -d-\lambda )(-\mu -\epsilon \tau -\lambda )\left|\begin{array}{ccc}-(\mu +\sigma +c)-\lambda & 0& \frac{(1-\epsilon ){\beta }_1{\text{Λ}}_1}{\mu +\epsilon \tau }\\ \epsilon & -(\alpha \rho +\varphi +\mu )-(1-\alpha \rho )-\lambda & 0\\ \frac{{\beta }_2{\text{Λ}}_2}{\mu +d}& 0& -\mu -\lambda \end{array}\right|=0$

The first part of the right-hand side is equated to zero,

$\begin{array}{ll}-{(\mu -\lambda )}^4(-\mu -d-\lambda )(-\mu -\epsilon \tau -\lambda )=0& (23)\end{array}$

λ_{1, 2, 3, 4} = −μ, λ = −(μ+d) and λ = −(μ+ετ) and the second part gives the following quadratic equation:

$\begin{array}{ll}{\lambda }^2+(2\mu +\sigma +c)\lambda +\mu (\mu +\sigma +c)-\frac{(1-\epsilon ){\beta }_1{\text{Λ}}_1{\beta }_2{\text{Λ}}_2}{(\mu +\epsilon \tau )(\mu +d)}=0& (24)\end{array}$

which can also be written as:

$\begin{array}{r}{\lambda }^2+(2\mu +\sigma +c)\lambda +\mu (\mu +\sigma +c)\\ \left[1-\frac{(1-\epsilon ){\beta }_1{\text{Λ}}_1{\beta }_2{\text{Λ}}_2}{\mu (\mu +\sigma +c)(\mu +\epsilon \tau )(\mu +d)}\right]=0& (25)\end{array}$

But from the basic reproduction number (Equation 11),

$\begin{array}{ll}{R}_O^2=\frac{{\beta }_1{\beta }_2{\text{Λ}}_1{\text{Λ}}_2(1-\epsilon )}{\mu (\mu +\sigma +c)(\mu +\epsilon \tau )(\mu +d)}& (26)\end{array}$

$J(x_{*})=\left[\begin{array}{ccccccccc}-\mu -(1-\epsilon ){\beta }_1{I}_M^{(*)}-\epsilon \tau & 0& 0& 0& 0& 0& 0& -(1-\epsilon ){\beta }_1{S}_F^{(*)}& 0\\ (1-\epsilon ){\beta }_1{I}_M^{(*)}& -(\mu +\sigma +c)& 0& 0& 0& 0& 0& (1-\epsilon ){\beta }_1{S}_F^{(*)}& 0\\ 0& \epsilon & -(\alpha \rho +\varphi +\mu )-(1-\alpha )\rho & 0& 0& 0& 0& 0& 0\\ 0& 0& \alpha \rho & -\mu & 0& 0& 0& 0& 0\\ 0& 0& (1-\alpha )\rho & 0& -\mu & 0& 0& 0& 0\\ \epsilon \tau & 0& 0& 0& 0& -\mu & d& 0& 0\\ 0& -{\beta }_2{S}_M^{(*)}& 0& 0& 0& 0& -{\beta }_2{I}_F^{(*)}-\mu -d& 0& 0\\ 0& {\beta }_2{S}_M^{(*)}& 0& 0& 0& 0& {\beta }_2{I}_F^{(*)}& -\mu & 0\\ 0& \sigma & \varphi & 0& 0& 0& 0& 0& -\mu \end{array}\right]$

Which means that (25) can be written as:

$\begin{array}{ll}{\lambda }^2+(2\mu +\sigma +c)\lambda +\mu (\mu +\sigma +c)(1-R_O^2)=0& (27)\end{array}$

If $R_O<1$, then the polynomial (Equation 27) will have negative roots. From this, coupled with the fact that in Equation 23, all the eigenvalues are negative (since all the parameters are positive), we can therefore conclude that, provided $R_O<1$, then all the eigenvalues of |J(x₀)−λI| are negative. This implies that the system is locally asymptotically stable at DFE [36].

3.5 Local stability of endemic equilibrium point EEP

The Jacobian (at EEP) of the model (Equation 2) is obtained as:

Where $S_F^{(*)},S_M^{(*)},I_F^{(*)}$ and $I_M^{(*)}$ are exactly as obtained in Section 3.3. The determinant |J(x_*)−λI| obtained as:

$\begin{array}{l}|J(x_{*})-\lambda I|=-{(\mu -\lambda )}^4·(-(\alpha \rho +\varphi +\mu )-(1-\alpha )\rho -\lambda )·\left|\begin{array}{cccc}-\mu -(1-\epsilon ){\beta }_1{I}_M^{(*)}-\epsilon \tau -\lambda & 0& 0& -(1-\epsilon ){\beta }_1{S}_F^{(*)}\\ (1-\epsilon ){\beta }_1{I}_M^{(*)}& -(\varphi +\mu +\epsilon )-\lambda & 0& (1-\epsilon ){\beta }_1{S}_F^{(*)}\\ 0& -{\beta }_2{S}_M^{(*)}& -{\beta }_2{I}_F^{(*)}-\mu -d-\lambda & 0\\ 0& {\beta }_2{S}_M^{(*)}& {\beta }_2{I}_F^{(*)}& -\mu -\lambda \end{array}\right|=0\\ & (28)\end{array}$

This implies that −(μ−λ)⁴(−(αρ+ϕ−μ)−(1−α)ρ−λ) = 0 which in turn boils down to:

λ_{1, 2, 3, 4} = −μ and λ₅ = −((αρ+ϕ+μ) + (1−α)ρ), which are all negative λs since all parameters are positive and α lies between 0 and 1.

The matrix term on the right side of Equation 28 becomes:

$\begin{array}{ll}{a}_0{\lambda }^4+a_1{\lambda }^3+a_3{\lambda }^2+a_4\lambda +a_5=0& (29)\end{array}$

Where

$\begin{array}{lll}{a}_0& =1& (30)\end{array}$

$\begin{array}{lll}{a}_1& =A+\mu +B+G& (31)\end{array}$

$\begin{array}{lll}{a}_2& =A(\mu +B+G)+(B\mu +G\mu +BG-EC)& (32)\end{array}$

$\begin{array}{lll}{a}_3& =A(B\mu +G\mu +BG-EC)+JHC-(B\mu G+ECF-ECG)& (33)\end{array}$

$\begin{array}{lll}{a}_4& =A(B\mu G+F-ECG)+JHCF+JHCG& (34)\end{array}$

and A = μ+(1−ε)β₁I_M+ετ, B = μ+σ+c, C = β₂S_M, G = β₂I_Fε+μ+d, F = β₂I_F, E = (1−ε)β₁S_F, H = E = (1−ε)β₁S_F, J = (1−ε)β₁I_M and F = β₂I_F. Using the Routh-Hurwitz Criterion, we construct an array for Equation 29;

$\begin{array}{cccc}{\lambda }^4& a_0& a_2& a_4\\ {\lambda }^3& a_1& a_3& 0\\ {\lambda }^2& \frac{a_4{a}_2-a_3}{a_4}& a_4& 0\\ {\lambda }^1& a_3& \frac{a_4{a}_1}{b_1}& 0\\ {\lambda }^0& a_4& 0& 0\end{array}$

For the system to be stable at EEP, the Routh Criterion demands that all the signs of each term in the first column of the array have to be the same. If this condition holds for the system, then it implies that HPV will invade the population when R₀>1 [37].

4 Estimation of state variable values

In this work, initial population and parameter estimates are made using information that as much as possible reflects the actual state of HPV-Cervical Cancer incidence in Nigeria, focusing more on HPV-16 infections because, according to Healthline (2023), HPV-16 is the most prevalent high-risk type of HPV [38]. Thus, this research focuses on the number of men and women infected with HPV 16.

• According to the 2023 report summary, the population pyramid of Nigeria found in ICO/IARC HPV information center document [39], if we assume an even distribution of males aged 5–9 years old and using the result of this assumption to obtain the number of males who are 9 years old, adding the result obtained to the rest of the population with higher ages, we will have the total number of male population aged 9+ years to be approximately 77,889,575 and the number of males aged 15 years and above to be 60,867,936. Doing the same for the female population, the total number of females aged 9+ will be approximately 76,081,014 and for those aged 15+, we will have 59,621,056. It is worth noting here that individuals from the age of 9+ years are considered in this work as the total population in order to maintain consistency with the age of girls who were targeted to get the HPV vaccine during the Nigerian vaccination program [40]

• The number of males and females infected with HPV 16 will be calculated using the number of individuals in the population aged 15 years and above because of the recorded increase in the proportion of those aged 15 and 19 years engaging in sexual acts in Nigeria [41].

Although Ohihoin et al. conducted research amongst five 570 participants in Nigeria to determine that only 3% of that sample had HPV-16 [42], the small sample size however is not enough to represent the broader population dynamics of HPV-16 prevalence in Nigeria as a result of its vast population and regional health differences. Moreover, global epidemiological standards suggests that bigger and more diverse samples are better used in order to obtain better estimates [43]. Consequently, we will be using the global prevalence instead to estimate the number of men infected with this strain. This approach will reduce any potential biases that are usually associated with limited local sampling.

According to Bruni et al. [39], the global prevalence of HPV-16 in men is 5%. With the absence of Nigeria-specific data for HPV 16 in men, this percentage is used alongside our estimated total number of males aged 15 years and above to obtain the approximate value of men with HPV 16 to be 3,043,397 (I_M). To obtain the approximate number of HPV-16 infected females I_F, we will also make use of a specific estimate of the global prevalence of HPV 16 which is 1.17%. Using this on the total number of females aged 15 and above, earlier estimated, we have approximately 697,566 females infected with only HPV-16.

• From Bruni et al. [39], the number of women with cervical cancer C_C is 12 075. This was the estimation of the number of women who are diagnosed annually with cervical cancer in Nigeria [39].

• According to an article from Johns Hopkins Medicine (n.d), it says that without intervention, approximately 90% of HPV-16 infections clear out on their own [44]. This implies that the body's immune system is capable of getting rid of this virus 6 to 18 months after infection, 90% of the time. Mathematically, 0.90 of infected males and females who recover naturally is approximately 2,739,057 and 627 809 respectively, in number. Only the number of females who recover R_F is, however, required in this research article.

• If 90% of females recover from HPV, it implies that the rest of the 10% do not recover but will go ahead to develop precancerous cells in their cervix. This number (Q) is approximately 69 757.

• The number of individuals in the V compartment is obtained from [45]. It says that over 12 million girls have received the HPV vaccine.

• There is, however, currently no record of any vaccination program for males in Nigeria, but it is possible that some of them must have gotten it through private health facilities, international schools, or during travel abroad. To reflect this rare but possible access, we conservatively assume that about 0.01% of the total male population (7,789) from 9 years above got the vaccine. This estimation is not based on national data but merely acknowledges the few individuals who had access to the vaccines through these expensive means highlighted above. Thus to avoid over stating we would take the total number of individuals in V compartment as 12,007,789 in number.

• A study was conducted in Nigeria across government facilities located in Kaduna, Rivers State, and Lagos in 2021. It was recorded that 96% of women who were diagnosed with Cervical Intraepithelial Neoplasia (Precancerous lesions) were treated [46]. Also, around June 2022, over 200,000 women were screened for cervical cancer, and those who were diagnosed with cervical lesions were treated on site [46]. Thus, we would be estimating the number of people in the T_Q compartment by assuming that 96% of those who have precancerous lesions were treated; 0.96 of Q is approximately 66,967.

All this information is summarized in Table 4:

Table 4

Table 4. Estimated initial values for state variables.

The values used for the state variables in this study were estimated/calculated with a focus on information from around the year 2022. Where data from 2022 were not available, values from nearby years were utilized, and global estimations were used where local information was unavailable. Reasonable assumptions were made to make sure that the estimates still reflected required conditions.

The numbers in Table 4 have been converted to manageable versions by rescaling to a population of 100,000 people, as shown in Table 5.

Table 5

Table 5. Scaled values of the state variables to a population of 100,000.

4.1 Estimation of the values of model parameters

To be able to analyze and simulate the dynamics of this model, it is important that we assign values that are realistic to the parameters of the model since these parameters describe the key biological processes involved in the transmission dynamics of HPV as it progresses to early cervical cancer. The estimation process was done in such a way that it reflects a realistic context of HPV transmission in Nigeria, particularly in 2022. Where national data were unavailable, global estimates and assumptions were made. Information from nearby years was utilized where data from 2022 was not available.

• In October 2023, Nigeria introduced the vaccine into its regular immunization program. The program's target was to make sure that 7.7 million girls aged 9 to 14 years old get vaccinated. The program, however, exceeded its target by vaccinating over 12 million girls [40, 47]. This initiative was introduced to significantly reduce the risk of these girls developing cervical cancer in the future. The vaccination rate τ is obtained by dividing the total vaccinated by the product of the target population and 2. This will result in 0.779 per year. Since the vaccination of girls is still ongoing [48], the total number vaccinated was annualized over 2 years to obtain a yearly vaccination probability (τ = 0.779 per year).

• Even though approximately 90% of HPV infections clear naturally within 1–2 years after infection [49], it has been reported to persist longer [50]. Therefore, an average duration of infection of approximately 2.2 years is assumed, and this corresponds to a recovery rate of σ = 0.45 per year, to reflect the slower clearance of high-risk HPV. Mathematically, $\sigma =\frac{1}{2.2yrs}.$

• To get an estimate of the probability that a susceptible female is enlightened about the HPV vaccine ε, a conservative estimate of 20% will be assumed due to results from recent studies carried out in different regions of the country. For instance, among university students in Niger state, it was recorded that only 25% were aware of the HPV vaccine [51]. Also, in Southwest Nigeria, it was revealed in a study that only about 27.1% of women were aware of the HPV vaccine [52]. Additionally, 26% of students from a tertiary institution in Plateau state, Nigeria, are aware of the HPV vaccine [51]. It is also important to add here that ε cannot be 0 because it would imply that no female in Nigeria knows about the vaccine, which is untrue. Similarly, ε cannot be equal to 1 because it would mean that every female in Nigeria is aware of the vaccine. This interesting scenario is further presented in Figures 16, 17.

• The rate at which individuals die naturally μ in Nigeria is 0.01243 per year. This means that about 12.43 deaths occur in 1,000 people yearly [53].

• The rate at which susceptible males move to the vaccination compartment d is assumed to be 0.0001 as previously stated in Section 4.0. This is to avoid overstating. Also, a conservative value has been assumed because there is currently no vaccination program for males in Nigeria.

• To obtain the probability that someone with CIN is screened successfully α, an estimate is done with respect to the information provided in Lawson et al. [46] and Clinton Health Access Initiative [54] as first mentioned in Section 4 (last bullet point). A 0.96 per year value for α will be used in order to avoid over estimating. Moreover, α lies between 0 and 1 to reflect realistic conditions, since α = 0 would mean that no one has been screened for Cervical Cancer or the pre-cancerous lesions, and α = 1 would imply that every woman has been screened. Both boundaries are unrealistic, and this study intends to, as much as possible, reflect the real world. The impact of these scenarios is shown in Figures 18, 19.

• To obtain ρ, the total rate at which CIN progresses out of the CIN compartment (Q), the reciprocal of the average duration an individual remains in the CIN stage (Q compartment) is computed. Since Q houses all the stages of CIN, the value of ρ has to mirror the average progression time across the three stages combined. Generally, CIN 1 takes over 10 years to progress. Studies have shown that it takes approximately 23.5 years for CIN2/3 to progress [6, 55]. In this research article, slow and fast progression rate is accounted for, and thus an average progression duration that an individual remains in the Q compartment is assumed to be 12 years. The reciprocal of this is 0.0833 per year.

• Moscicki et al. reported that 60% of CIN 1 cases clears out within 1 year [56]. A 2022 study based on vaccines also reported that 76.4% of vaccinated women experienced regression of CIN 1 lesions and 68.1% of CIN 1 cases among unvaccinated women cleared out [57]. Maintaining a middle ground in this study, an estimate of 0.7 per year is assumed. This is the estimate of the rate at which CIN 1 recovers phi in this work naturally.

• The cumulative progression rate from HPV-16 infection to CIN 1 has been reported to be around 20.7% [55]. This means that c = 0.207 per year, which will be approximated to 0.2 per year in this study.

• The basic reproduction number $R_O$ estimate for HPV 16 ranges from approximately 1.03 to 4.1, based on studies conducted by Ribassin-Majed et al. [58]; Li et al. [59] and Gao et al. [60]. This model will, however, adopt a conservative but realistic estimate of 1.5, considering the low vaccination coverage and enlightenment level in Nigeria. This value will be used in Equations 12, 13 to compute β₁ and β₂. The resulting estimate will be approximately 0.00017 and 0.0000016, respectively.

• A stable population usually has its age distribution spread evenly across all age groups, especially in a large population [61]. Since Nigeria has a very large population, it will be assumed in this model that the age group is spread evenly. In 2023, Statista reported that the life expectancy at birth in Nigeria is approximately 61.79 years, where 64 years is for the females and 60 years is for the males [62]. This model assumes that individuals enter the susceptible compartments through aging (when they turn 9). The recruitment rate of females into the susceptible class S_F is 1,345,154 girls per year if we assume an even age distribution, while the recruitment rate of males into S_M is 1,456,448 boys per year. Rescaling this to per 100,000, we obtain 1,818.2 girls/year and 1,960.8 boys/year. Mathematically, $\text{Λ}=\frac{\mathrm{\text{Number of individuals aged from 9 up to life expectancy age}}}{\mathrm{\text{Life expectancy age}}-9}$, and the result was rescaled to per 100,000 as $\frac{\text{Λ}\times 100,000}{\mathrm{\text{Number of individuals aged from 9 up to life expectancy age}}}$.

The Table 6 summarizes the parameter values.

Table 6

Table 6. Estimated parameter values used in the model.

It is important to re-emphasize that the parameter estimations were done in such a way that they try as much as possible to reflect the context of HPV transmission in Nigeria, particularly in 2022. Where national or local information was not available, global estimates and assumptions were made. Information from nearby years was also utilized, where data from 2022 were unavailable.

5 Sensitivity analysis

Sensitivity analysis was carried out in model (2) in order to promote understanding of how the parameters in the model affect the spread of HPV-16. A common approach in epidemiological modeling, known as the normalized forward sensitivity index method, was used as described in Chitnis et al. [63] and Rodrigues et al. [64]. This method will help us to determine which parameter influences the basic reproduction number the most. The formula is given as:

Sensitivity Index $P_k=\frac{k}{R_O}\frac{d{R}_O}{dk}$

Where $R_O$ is the basic reproduction number obtained in Equation 11 and k could be any of the parameters in the $R_O$ expression (Equation 11).

The sensitivity expressions are given below:

$\begin{array}{r}{P}_{{\beta }_1}=\frac{1}{2},P_{{\beta }_2}=\frac{1}{2},P_{{\text{Λ}}_1}=\frac{1}{2},P_{{\text{Λ}}_2}=\frac{1}{2}\end{array}$

$P_{\mu }=-\frac{1}{2}\left[1+\frac{\mu }{\mu +\sigma +c}+\frac{\mu }{\mu +\epsilon \tau }+\frac{\mu }{\mu +d}\right]$

$P_{\sigma }=\frac{-\sigma ({\mu }^3+{\mu }^{2}d+\epsilon {\mu }^2\tau +\epsilon \tau d\mu )}{2\mu (\mu +\sigma +c)(\mu +\epsilon \tau )(\mu +d)}$

$P_{\tau }=\frac{-\epsilon \tau }{\mu +\epsilon \tau }$

$P_d=\frac{-d}{2(\mu +d)}$

$P_c=\frac{-c}{2(\mu +\sigma +c)}$

$P_{\epsilon }=\frac{-\epsilon \tau }{\mu +\epsilon \tau }$

The Table 7 gives us more information about the sensitivity index of various parameters in the model. Values were obtained based on the estimated values provided in Table 6.

Table 7

Table 7. Sensitivity indices of model parameters with respect to $R_O$.

If P_k < 0, it physically means that the disease will die out because increasing the parameter k reduces $R_O$, but if P_k>0, then the disease will spread out because parameter k increase will amplify the basic reproduction number [64].

To further understand how $R_O$ responds to variations in key parameters, 3D mesh plots were generated to illustrate the relationship between the basic reproductive number and some parameters.

The first is a mesh plot to show the relationship between $R_O$, awareness probability ε, and β₁, which is the contact rate for women.

In Figure 2, the x-axis represents the contact rate β₁ for women, the y-axis represents awareness probability (ε), and the z-axis is $R_0$. Red-brown zones represent higher values of $R_0$. It can be observed that a higher basic reproduction number implies a zero or low awareness level and a higher contact rate (β₁). For instance, when the awareness level is zero, $R_0$ is 24.29 and 31.17 when β₁ is equal to 0.0001327 and 0.0002184, respectively. Indicating that higher contact rates result in higher basic reproduction numbers. Observing values close to the yellow and blue colors, on the other hand, it can be seen that higher values of awareness increase the basic reproduction number. For instance, when awareness (ε) has the values 0.8571 and 0.9592, $R_0$ has the values 1.562>1 and 0.835 < 1 respectively.

Figure 2

Figure 2. $R_O$ sensitivity to ε and β₁.

In Figure 3, the x-axis is τ, the y-axis stands for Λ₁, and the z-axis represents $R_0$. Regions of high $R_0$ (toward the yellow area) are observed to have a low vaccination rate τ. For instance, 0 vaccination rate has $R_0$ value of 25.26 and τ = 0.6735 gives a lower value of $R_0=7.284$. When τ increases further to 1.469, there is a much lower $R_0=3.73$ (regions of darker blue shades). This further stresses the impact of vaccination programs in curbing HPV infections.

Figure 3

Figure 3. $R_O$ sensitivity to Λ₁ and τ.

5.1 Global sensitivity analysis

In order to fully grasp the impact of the model parameters on HPV transmission and progression to cervical cancer, it is pertinent to transcend the local sensitivity analysis, owing to the fact that it only explains how small changes in parameters influence the model around a specific point. Global Sensitivity Analysis GSA), however, addresses this constraint because it analyses the impact of parameters across their entire possible ranges. It helps to understand the behavior of the model when subjected to different conditions, making it possible to identify factors that are most influential in disease spread [65].

The Partial Rank Correlation Coefficient (PRCC) method is used to carry out GSA in this study primarily because it accurately measures how input parameters and corresponding model outcomes relate, even in cases of non-linear relationship [66].

Figure 4 below displays the tornado plot of PRCCs of the parameters of the model (Equation 2). In the figure, varepsilon is the name given to the symbol ε, which stands for the awareness probability. It is observed that when parameters with negative PRCCs are elevated, the basic reproduction number R₀ decreases, while an increase in parameters with positive PRCCs amplifies the basic reproduction number $R_O$.

Figure 4

Figure 4. PRCC analysis of model parameters on the basic reproductive number.

6 Optimal control

6.1 Optimal control problem formulation

Promoting awareness about the existence of HPV and its vaccines, and the availability of Screening programs, all play a vital role in reducing the spread of HPV and HPV induced cervical cancer cases. In this study, reducing the number of infected females and cervical cancer cases is our target. This will be accomplished by making use of: (i) screening as a control measure v₁, (ii) awareness, control parameter v₂. Therefore, the optimal control problem below is considered.

$\begin{array}{ll}minR(v_1^{*}{v}_2^{*})={\displaystyle {\int }_0^t}(f_1{C}_C+f_2{I}_F+f_3{v_1}^2+f_4{v_2}^2))dt& (35)\end{array}$

with f₁, f₂, f₃ and f₄ as positive weights which align with the advantages of reducing the spread of HPV and cervical cancer cases over a specific period of time T. The aim is to reduce the number of infected female individuals and cases of Cervical Cancer. Consequently, we seek an optimal control $(v_1^{*},v_2^{*})$ in such a way that

$\begin{array}{ll}R(v_1^{*}{v}_2^{*})={\displaystyle \underset{v_iϵV}{min}}R(v_1,v_2),fori=1,2& (36)\end{array}$

Where v₁ and v₂ are the screening (α) and awareness (ε) parameters, respectively, as provided in our model Equation 2, which is now rewritten below as system (Equation 37):

$\begin{array}{l}\begin{array}{ll}\frac{d{S}_F}{dt}& ={\text{Λ}}_1-\mu S_F-(1-v_2){\beta }_1{S}_F{I}_M-v_2\tau S_F\\ \frac{d{I}_F}{dt}& =(1-v_2){\beta }_1{S}_F{I}_M-(\mu +\sigma +c)I_F\\ \frac{dQ}{dt}& =c{I}_F-(v_1\rho +\varphi +\mu )Q-(1-v_1)\rho Q\\ \frac{d{T}_Q}{dt}& =v_1\rho Q-\mu T_Q\\ \frac{d{C}_C}{dt}& =(1-v_1)\rho Q-\mu C_C\end{array}\end{array}$

$\begin{array}{ll}\begin{array}{ll}\frac{d{R}_F}{dt}& =\sigma I_F-\mu R_F+\varphi Q\\ \frac{d{S}_M}{dt}& ={\text{Λ}}_2-{\beta }_2{S}_M{I}_F-\mu S_M-d{S}_M\\ \frac{d{I}_M}{dt}& ={\beta }_2{S}_M{I}_F-\mu I_M\\ \frac{dV}{dt}& =-\mu V+v_2\tau S_F+d{S}_M\end{array}& (37)\end{array}$

Where,

$V=(v_1(t),v_2(t)ϵ{[L^2(0,T)]}^2):0\le v_i(t)\le 1,tϵ\left[0,T\right],$ for i = 1, 2

We will be using Pontryagin's Maximal Principle to get the most effective control solution [67]. This principle requires showing that an optimal solution exists and also that optimal conditions are defined as described in Aguegboh et al. [68] and Mohsen and Naji [69].

We shall first proceed to show the existence of optimal control.

6.1.1 Optimal control existence

The existence of optimal controls $v_1^{*}$ and $v_2^{*}$ follows from standard optimal control theory as shown in Kirk [70]. The conditions for existence are satisfied for the following reasons:

• The first being that the control set

$V=\left\{(v_1,v_2)\in {\left[L^2(0,T)\right]}^2:0\le v_i\le 1\right\}$

• is closed and convex.

• The control variables appear in the model equations only as linear terms. Also, all population sizes stay within realistic bounds.

• The objective functional we intend to minimize is convex in the controls as a result of the quadratic cost terms: ${v_1}^2$ and ${v_2}^2$.

6.1.2 Optimality control system

The Optimality system will be established in this section; With v₁(t) and v₂(t) possessing a value in [0, 1], consider that v₁(t) and v₂(t), be control functions which are measurable on [0, T] [68]. The objective function $R$ of 35 will then be minimized using optimal controls v₁(t) and v₂(t) where

$\begin{array}{r}{v_1}^{*}=max\left\{min\left\{\widetilde{v_1},1\right\},0\right\},\\ {v_2}^{*}=max\left\{min\left\{\widetilde{v_2},1\right\},0\right\}\\ \widetilde{v_1}=\frac{({\zeta }_5-{\zeta }_4)\rho Q}{2f_3},\\ \widetilde{v_2}=\frac{({\zeta }_2-{\zeta }_1){\beta }_1{S}_F{I}_M+({\zeta }_1-{\zeta }_9)\tau S_F}{2f_4}\end{array}$

such that S_F, I_F, R_F, V, S_M, I_M, C_C, T_Q, Q now corresponds to the solution of the system of Equation 37.

To prove this, there is a need to write out the Hamiltonian $H$ of the resulting optimal control problem:

$H=f_1{C}_C+f_2{I}_F+f_3{v_1}^2+f_4{v_2}^2+{\zeta }_1({\text{Λ}}_1-\mu S_F-(1-v_2){\beta }_1{S}_F{I}_M-v_2\tau S_F)$

+ζ₂((1−v₂)β₁S_FI_M−(μ+σ+c)I_F)

+ζ₃(cI_F−(v₁ρ+ϕ+μ)Q−(1−v₁)ρ)

+ζ₄(v₁ρQ−μT_Q)

+ζ₅((1−v₁)ρQ−μC_C)

+ζ₆(σI_F−μR_F+ϕQ)

+ζ₇(Λ₂−β₂S_MI_F−μS_M−dS_M)

+ζ₈(β₂S_MI_F−μI_M)

+ζ₉(μV+v₂τS_F+dS_M)

such that

ζ_i(T), where i = 1, 2, 3, 4, 5, 6, 7, 8, 9 refers to the adjoint variables that satisfies ζ_i(T) = 0, provided by the canonical equations below:

$\begin{array}{l}\frac{d{\zeta }_1}{dt}=-\frac{\partial H}{\partial S_F}=\\ -\left[{\zeta }_1(-\mu -(1-v_2){\beta }_1{I}_M-v_2\tau )+{\zeta }_2(1-v_2){\beta }_1{I}_M+{\zeta }_9{v}_2\tau \right]\\ \frac{d{\zeta }_2}{dt}=-\frac{\partial H}{\partial I_F}=\\ -\left[f_2-{\zeta }_2(\mu +\sigma +c)+{\zeta }_{3}c-{\zeta }_7{\beta }_2{S}_M+{\zeta }_8{\beta }_2{S}_M+{\zeta }_6\sigma \right]\\ \frac{d{\zeta }_3}{dt}=-\frac{\partial H}{\partial Q}=\\ -\left[-{\zeta }_3(v_1\rho +\varphi +\mu )+{\zeta }_4{v}_1\rho +{\zeta }_5(1-v_1)\rho +{\zeta }_6\varphi \right]\\ \frac{d{\zeta }_4}{dt}=-\frac{\partial H}{\partial T_Q}=-\left[-{\zeta }_4\mu \right]\\ \frac{d{\zeta }_5}{dt}=-\frac{\partial H}{\partial C_C}=-\left[f_1-{\zeta }_5\mu \right]\\ \frac{d{\zeta }_6}{dt}=-\frac{\partial H}{\partial R_F}=-\left[-{\zeta }_6\mu \right]\\ \frac{d{\zeta }_7}{dt}=-\frac{\partial H}{\partial S_M}=-\left[-{\zeta }_7{\beta }_2{I}_F-{\zeta }_7\mu -{\zeta }_{7}d+{\zeta }_8{\beta }_2{I}_F+{\zeta }_{9}d\right]\\ \frac{d{\zeta }_8}{dt}=-\frac{\partial H}{\partial I_M}=-\left[-{\zeta }_1(1-v_2){\beta }_1{S}_F+{\zeta }_2(1-v_2){\beta }_1{S}_F-{\zeta }_8\mu \right]\\ \frac{d{\zeta }_9}{dt}=-\frac{\partial H}{\partial V}=-\left[{\zeta }_9\mu \right]& (38)\end{array}$

Using the Pontryagin Maximal Principle, the following optimum conditions are established:

$\begin{array}{ll}\frac{\partial H}{\partial v_1}=2f_3{v}_1+{\zeta }_4\rho Q-{\zeta }_5\rho Q=0& (39)\end{array}$

Solving Equation 39 using the state and the adjoint variables, we have:

$\begin{array}{ll}\widetilde{v_1}=\frac{({\zeta }_5-{\zeta }_4)\rho Q}{2f_3}& (40)\end{array}$

In a similar manner, for $\frac{\partial H}{\partial v_2}=0$, we have:

$\begin{array}{ll}\widetilde{v_2}=\frac{({\zeta }_2-{\zeta }_1){\beta }_1{S}_F{I}_M+({\zeta }_1-{\zeta }_9)\tau S_F}{2f_4}& (41)\end{array}$

For the most effective control measures to implement; v₁ and v₂ the signs and control constraints of $\frac{\partial H}{\partial v_1}$ and $\frac{\partial H}{\partial v_2}$, are considered.

Demonstrating this:

$\begin{array}{l}{v_1}^{*}=\{\begin{array}{cc}0& if\frac{\partial H}{\partial v_1}<0\\ \\ \widetilde{v_1}& if\frac{\partial H}{\partial v_1}=0\\ \\ 1,& if\frac{\partial H}{\partial v_1}>0\end{array}\end{array}$

and ${v_1}^{*}=max\left\{min\left\{\widetilde{v_1},1\right\},0\right\}$ where $\widetilde{v_1}=\frac{({\zeta }_5-{\zeta }_4)\rho Q}{2f_3}$

$\begin{array}{l}{v_2}^{*}=\{\begin{array}{cc}0& if\frac{\partial H}{\partial v_2}<0\\ \\ \widetilde{v_2}& if\frac{\partial H}{\partial V_2}=0\\ \\ 1,& if\frac{\partial H}{\partial v_2}>0\end{array}\end{array}$

and $\widetilde{v_2}=\frac{({\zeta }_2-{\zeta }_1){\beta }_1{S}_F{I}_M+({\zeta }_1-{\zeta }_9)\tau S_F}{2f_4}$

The optimal solution will be obtained when ${v_1}^{*}$ and ${v_2}^{*}$ are substituted into the system (Equation 37),

7 Numerical simulation

In both simple and complex disease models, numerical methods help us to study the behavior of the disease under different or varying conditions. Numerical Simulation provides researchers with an effective way of exploring the dynamics of disease transmission, treatment impacts, and the effects of intervention techniques like screening and vaccination [71].

In this study, the mathematical model 2 is solved using the Runge-Kutta 4th order method (RK4) implemented in MATLAB. This method happens to be an effective method of solving ordinary differential equations (ODEs), providing a balance between a high accuracy and stability [72].

The Runge-Kutta 4th order method creates a precise forecast of the dynamics of the population in each compartment by approximating the solutions of the ordinary differential equations at small time steps. In biological models where accurate tracking of Susceptible, Infected, and other compartments of interest is important for understanding the transmission and control of diseases, this approach happens to be particularly effective [73].

In this research, the simulation was carried out using recent data available for the year 2022 as initial values for each compartment. Every estimate made was based on available information, which, as much as possible, reflected the actual state of HPV-Cervical Cancer incidence in Nigeria during that year. Some parameters were also assumed, and reasons were provided for some of the assumptions. The model is then solved over a 10-year period to study the way the disease progresses and the effect of various intervention strategies, such as awareness, vaccination, and screening, on the population classes of interest.

The following graphical results were obtained using information contained in Table 7.

7.1 Graphical results and discussion

Figure 5 the susceptible class S_F (blue- solid line) exhibits a decrease during a 10-year period due to their movement into the infectious female I_F compartment as a result of contact with the Infectious male population I_M. The decrease is also due to their natural death rate μ = 0.01243 and the vaccination rate τ = 0.779. Also, with $R_0>1$ as estimated in this study, it is only reasonable that the susceptible compartment experiences a decrease because of the disease persistence.

Figure 5

Figure 5. Population dynamics of all the compartments over time when $R_0$ (0 to 10 years).

The orange curve (dashed line), which represents the infectious female compartment I_F, displays a slight increase due to the transmission rate β₁ = 0.00017 with $R_0>1$. It, however, starts to decrease after 2.5 years at a slow rate.

The compartment that houses women with cervical dysplasia Q (the yellow-orange-dash dot line) experiences a slight increase at first due to the entry of I_F individuals into the compartment, then it starts reducing gradually as they move to the T_Q compartment.

The T_Q compartment (purple–dotted line) experiences a marginal growth over time due to a high screening success rate α = 0.96, which facilitates a transition from Q to T_Q. This correctly captures this study's emphasis on successful screening for Q individuals, which is consistent with realistic expectations considering the high value of α = 0.96 that was used in this study. It can be observed that after an initial increase, it stays at a constant level after six years.

The C_C compartment (green- solid line) stays flat in this plot. This is as a result of the high screening rate α = 0.96 chosen in this study. Individuals would tend to move more to the treatment compartment due to the high screening efficacy in preventing cervical cancer advancement.

Recovered females R_F (Blue– dashed lines curve) go through a slight growth. This is primarily due to the value of the natural recovery rate from I_F, which is (σ = 0.45), coupled with the high natural recovery rate of 0.7 for the Q individuals.

The Susceptible male individuals S_M (dark red–dash dot line curve) experience a gradual increase in the space of 10 years because of the assumed low male vaccination rate d = 0.0001. It implies that very few male individuals get to move into the V compartment to the point that no significant change or decrease (in this particular case) is seen in the S_M compartment, instead they just keep growing as a result of the recruitment rate Λ₂ = 1960.8

The infectious male individuals I_M, represented by the dark blue (dotted line) curve, grow slowly, reflecting the small value of the transmission rate β₂. They continue to grow at a very slow rate over the ten-year period. For instance, at year 8, I_M = 5871 and at the 9th year, infected male individuals rise to 6,074 in number.

The vaccination individuals V (Red Orange/brownish—solid line) displays a steady increase due to male and female vaccination. This emphasizes the advantage of vaccination programs for combating infections.

In Figure 6, the I_F compartment experiences an initial increase with high β₁ values. However, a reduction is seen when β₁ = 0.000006 (blue curve), a considerably low transmission rate. This indicates that a decrease in male-to-female transmission rate reduces the rate at which susceptible females become infected with HPV-16. A sharp increase is displayed when β₁ = 0.00017, our initial estimate. This suggests a rather aggressive spread of the HPV-16 in the I_F compartment. This result is consistent with realistic expectations because, in reality, a higher contact rate between susceptible females and infectious males will definitely give rise to more HPV-16 infections among females. A gradual decrease can also be observed after some time in this graph for all values of β₁ as I_F individuals proceed to the Q compartments while some recover naturally and move to R.

Figure 6

Figure 6. Population dynamics of infectious female population over time with variations in β₁ (0 to 10 years).

In Figure 7, consistent with realistic expectations are the resulting curves that display an increase in cases of cervical cancer as β₁ increases. It also shows that the initial estimated value β₁ = 0.00017, which caused an increase in the infectious women population, also led to an increase in cervical cancer cases. The curve shows that the smaller the contact rate, the smaller the Cancer cases.

Figure 7

Figure 7. Population dynamics of Cervical Cancer female population over time with variations in β₁ (0 to 10 years).

Figure 8 matches the realistic expectation that a larger β₂ (that is, higher sexual contact rate between infected females and susceptible males) gives rise to higher HPV-16 infections in the male population. At the initial estimate (0.0000016), however, which happens to be the lowest value of β₂ in this graph, the number of infected males grows gradually, but at a lower rate compared to the curves with higher contact rates, that is, for β₂ = 0.00007 − 0.0009. Higher values of β₂ (0.00007–0.0009) give rise to a quick increase but become steady at a high level. It can also be observed that at values of β₂ = 0.0007 and 0.0009, the curves appear to overlap, indicating that increasing β₂ any further at some point does not affect the number of infected males. This means that at certain higher values of β₂, all individuals who can be infected have been infected, and the others are probably protected through vaccination. This result shows that implementing early intervention to decrease the $R_0$ is very crucial because control measures might become less effective when the infection has already spread widely.

Figure 8

Figure 8. Population dynamics of Infectious Male population over time with variations in β₂ (0 to 10 years).

Figure 9 shows that an increase in the contact rate β₂ also gives rise to an increase in the number of infectious females. This is most likely due to the fact that there will be a lot of infected males around to infect the susceptible females. It is logical to look at it this way because a higher β₂ means that a lot of males are getting infected (also evident in the previous Figure 8). It is observed that the initial estimated value of β₂ = 0.0000016 gradually increases I_F while the curves with higher values of β₂= (0.00007–0.0009) initially lead to a sharper increase of I_F individuals before it reduces and stays at a steady level. For all the different values of β₂ implemented here, all the curves reduce and stays at an almost constant level of about 2,760 individuals at approximately 8.49 years due to their movement to Q, R_F, and V compartments.

Figure 9

Figure 9. Population dynamics of the Infectious female population over time with variations in β₂ (0 to 10 years).

Figure 10 above displays the effect of screening (α) variations on the T_Q compartment. Since Clinton Health Access Initiative [54] and Lawson et al. [46] has actually shown that a very large proportion of those who were diagnosed of precancerous lesions during different screening programs are usually given treatment often in the same day, this graphical result relates the effect of screening on the treatment compartment in the sense that, a large value of the screening probability α gives rise to more individuals receiving treatment. This is consistent with realistic expectations.

Figure 10

Figure 10. Population dynamics of Treated CIN cases T_Q over time with variations in α (0 to 10 years).

In Figure 11, the effect of screening on Cervical Cancer cases can be observed. When α is as low as 0.2 (green curve), there is an increase in cervical cancer cases. A high screening probability value α = 0.96 (blue curve), however, gives rise to a welcome decrease in the number of cervical cancer cases. This result is consistent with realistic expectations because high screening implies that more people get treated for Precancerous lesions and move into the T_Q compartment, depopulating the C_c Class.

Figure 11

Figure 11. Population dynamics of Cervical cancer cases C_C over time with variations in α (0 to 10 years).

Next, we investigate the impact of varying the enlightenment value ε on Cervical cancer cases, HPV infections, and other compartments.

This graphical display (Figure 12) is biologically realistic because a high enlightenment probability value ε = 0.9 (green curve), gave rise to a remarkable reduction of the susceptible female population due to an increased movement into the vaccinated compartment. On the other hand, a lower value of ε = 0.2, denoted by the blue curve, causes a decrease in S_F population over time into the vaccinated compartment but at a slower rate.

Figure 12

Figure 12. Population dynamics of susceptible females S_F over time with variations in ε (0 to 10 years).

In Figure 13, the vaccinated compartment experiences a huge increase (green curve) when ε = 0.9. The initial estimate ε = 0.2 (blue curve), however, gives rise to an increase as well but at a slower rate- this outcome is biologically realistic.

Figure 13

Figure 13. Population dynamics of Vaccinated compartment V over time with variations in ε (0 to 10 years).

Very little increase is seen in Figure 14 (green curve) when the enlightenment value is very high (ε = 0.9). This is to be expected because enlightening more females about the existence of HPV and its vaccine will increase their tendency to partake in vaccination programs. Thus, the susceptible females move into the vaccinated compartment, and a smaller number of people become infected. A lower enlightenment value (initial estimate of 0.2) gives rise to high infection of females because fewer become enlightened and move to the vaccinated class.

Figure 14

Figure 14. Population dynamics of Infected female compartment I_F over time with variations in ε (0 to 10 years).

Figure 15 records a realistic behavior of the number of cervical cases with respect to variations in ε. From the curves, it can be seen that higher enlightenment probability values do not really give rise to much increase in cervical cancer cases compared to the curves produced by lower ε values.

Figure 15

Figure 15. Population dynamics of Cervical Cancer cases C_C over time with variations in ε (0 to 10 years).

To mathematically study what happens to the system during unrealistic conditions, that is, when α = 0, 1 and when ε = 0, 1, plots of I_F and V when ε = 0, 1 are presented, and what happens in the T_Q and C_c compartments when α = 0, 1 are displayed graphically as well.

Figure 16 displays what happens in the Infected I_F and Vaccinated V compartment when awareness ε = 0. The Vaccinated compartment (blue curve) reduces very slowly from its initial number of 7,799 individuals, primarily due to the natural death of Individuals in that compartment, since there is no one coming in from the susceptible female class because of awareness. The Infected compartment I_F increases drastically and gets to a peak of 14,000 before coming to a gradual decrease due to natural recovery, natural mortality, and their movement to the Q compartment.

Figure 16

Figure 16. Population dynamics of Infectious Female I_F and Vaccinated class V over time with awareness ε = 0 (0 to 10 years).

Figure 17 shows that when awareness or enlightenment = 1, the vaccinated compartment (blue curve) keeps rising, and the Infected I_F compartment stays at an almost zero level.

Figure 17

Figure 17. Population dynamics of infectious female I_F and vaccinated class V over time with awareness ε = 1 (0 to 10 years).

In Figure 18, the Early Cervical Cancer (blue curve) keeps rising (from the initial value of 8) as a result of the absence of screening programs, while the T_Q compartment stays at a constant initial value of 43, which it started with, without experiencing any change.

Figure 18

Figure 18. Population dynamics of Treated T_Q and Early Cervical Cancer C_C over time with screening α = 0 (0 to 10 years).

Figure 19 shows that T_Q (red curve) compartment keeps rising until it gets to the highest value (if the number of years is extended from 10 years) when screening is 1. It also shows the Early Cervical Cancer (blue curve) appearing to constantly be at the same level as the initial value it first started with.

Figure 19

Figure 19. Population dynamics of treated T_Q and early cervical cancer C_C over time with screening α = 1 (0 to 10 years).

In order to decide the most impactful approaches for reducing both Cervical cancer cases and HPV infections, an optimal control problem was developed and solved. The optimal controls for screening (v₁) and awareness (v₂), which are time-dependent, were numerically computed using a method called the Forward-Backward Sweep Method (FBSM) that is premised on Pontryagin's Maximum Principle (PMP). The following Figures 20, 21 demonstrate the projected dynamics of Infected females (I_F) and Early cervical cancer (C_C) in the presence and absence of these optimal control profiles over a period of ten years.

Figure 20

Figure 20. Population dynamics of treated I_F and Early Cervical Cancer C_C with controls (0 to 10 years).

Figure 21

Figure 21. Population dynamics of treated I_F and Early Cervical Cancer C_C without controls (0 to 10 years).

In Figure 20, the Cervical Cancer curve (blue) starts from its initial value of 8 and decreases very slowly while the C_C compartment reduces from 453 until it is almost touching zero. What this means is that the combined effect of the optimal controls most likely caused an increase in screening and awareness, leading to increased treatment and vaccination, respectively, resulting in the decrease in C_C and I_F individuals.

Figure 21 shows an initial increase of I_F population without controls. It signifies the negative effect the absence of controls has on the Infected and Cancer population. I_F (red line) increases and gets to the peak of 14,000 before it begins to decrease, probably due to their movement to the Q compartment, natural recovery, and mortality. Cervical cancer cases (blue curve) also appear to continuously rise in the absence of the combined effect of the controls.

This shows that the combined effect of awareness and screening plays a vital role in reducing HPV infections in females and Cervical Cancer.

8 Conclusion

In this research study, nine ordinary differential equations were formulated to explain the transmission dynamics of high-risk human papillomavirus, showing clearly how it leads to cervical cancer in women. Analysis showed that the formulated model is biologically meaningful and epidemiologically well posed, evident from the boundedness proof as well as proof of the existence of an invariant region, k. The basic reproduction number, which was obtained using the next generation matrix approach, provided an insight regarding the transmission or spread of the virus; A close study of this threshold parameter emphasizes the importance of interventions such as enlightenment, vaccination, and screening. It is observed that when the parameters that match these interventions increase, the basic reproduction number approaches a value less than 1. This is consistent with epidemiological studies that reveal that if $R_O<1$, the disease fizzles out with time [55]. Furthermore, it was also shown that the number of infected female individuals and cervical cancer cases increases as the basic reproduction number of the female population ($R_{O1}$) increases, this is also consistent with saying that if $R_O>1$, then the diseases continues to spread with time [55]. The sensitivity analysis shed more light on the fact that when parameters that represent enlightenment, vaccination rates, and natural recovery of females are increased, HPV-16 infections fizzled out with time because an increase in these parameters will reduce the basic reproduction number. Using Pontryagin's Maximum principle to carry out the optimal control of the model provided a great insight into effective strategies to curb the spread of HPV and reduce cervical cancer cases. The incorporation of the cervical dysplasia compartment in the model, which is often ignored in existing models, revealed a very important stage for intervention before cervical cancer kicks in. Including a treatment class (facilitated by screening) for the CIN cases (which is not found in existing models) was also a deliberate attempt to reduce the number of cervical cancer cases. This novel approach highlights the need and effectiveness of screening programs.

8.1 Contributions to public health and future research

The implication of our findings suggests that an increase in awareness campaigns would lead to a tremendous increase in the number of vaccinated female individuals. This will also go a long way to reinforce public health policies around vaccination against HPV. Another implication of our findings highlights that screening programs could substantially mitigate disease impact, specifically, reducing the number of cervical cancer cases in Nigeria and possibly in the world at large. Thus, integrating routine cervical cancer screening and awareness campaigns into national HPV prevention programs is important. This is very much in line with WHO's recommendation as referenced earlier in the introductory section [7]. Our research findings recommend that Health authorities should focus on educating communities in order to inform and encourage early detection and reduce the stigma that surrounds screening. Increasing access to affordable health screenings and prompt care for cervical dysplasia is crucial, particularly in underprivileged regions. Although vaccination remains important, especially where controlling the disease on a long-term basis is concerned, intensified awareness and screening efforts promote HPV and Cervical Cancer reduction. This model provides a practical tool for decision makers and professionals in the public health sector to analyze the probable results of different control measures and also help them to plan data-driven decisions to decrease the strain of HPV-induced Cervical Cancer.

In summary, the result of this study shows that a reduction in contact rates and an increase in awareness campaigns reduce HPV infection, while intensified screening efforts reduce Cervical cancer cases.

Future researchers could extend this model by investigating the possible impact of male vaccination on the model. This would further provide insight into the role of the male population in HPV transmission. They could also extend the work by integrating age-based modeling and stochastic effects to address real-world changes in infection rates and access to intervention strategies.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

OR: Conceptualization, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. JA: Conceptualization, Supervision, Validation, Writing – review & editing. PK: Conceptualization, Supervision, Validation, Writing – review & editing. OC: Conceptualization, Supervision, Validation, Visualization, Writing – review & editing. NA: Conceptualization, Supervision, Validation, Writing – review & editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work was supported by the Pan African University Institute for Basic Sciences Technology and Innovation (PAUSTI).

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. National Library of Medicine. Cervical Cancer: Learn more-Human Papillomavirus (2024). Available online at: https://www.ncbi.nlm.nih.gov/books/NBK279260/ (Accessed June 15, 2025).

Google Scholar

2. Planned Parenthood. Human Papillomavirus (HPV) (2024). Available online at: https://www.plannedparenthood.org/learn/stds-hiv-safer-sex/hpv (Accessed June 15, 2025).

Google Scholar

3. NYU Langone Health. Types of Human Papillomavirus (n.d). Available online at: https://nyulangone.org/conditions/human-papillomavirus/types (Accessed June 15, 2025).

Google Scholar

4. World Health Organization. Human Papillomavirus and Cancer (2024). Available online at: https://www.who.int/news-room/fact-sheets/detail/human-papilloma-virus-and-cancer (Accessed June 15, 2025).

Google Scholar

5. CDC. Basic Information about HPV and Cancer (n.d.). Available online at: https://www.cdc.gov/cancer/hpv/basic-information.html (Accessed June 15, 2025).

Google Scholar

6. National Cancer Institute. Cervical Cancer Causes, Risk Factors, Prevention (2024). Available online at: https://www.cancer.gov/types/cervical/causes-risk-prevention (Accessed June 17, 2025).

Google Scholar

7. World Health Organization. Human papillomavirus (HPV) and Cancer (n.d.). Available online at: https://www.who.int/news-room/fact-sheets/detail/human-papilloma-virus-and-cancer (Accessed July 2, 2025).

Google Scholar

8. M Department of Health. Quick Facts - HPV associated cancer (2024). Available online at: https://www.health.state.mn.us/data/mcrs/data/qfhpv.html (Accessed July 2, 2025).

Google Scholar

9. CDC. HPV Vaccine Recommendations. (n.d.). Available online at: https://www.cdc.gov/hpv/hcp/vaccination-considerations/index.html (Accessed July 2, 2025).

Google Scholar

10. Canadian Cancer Society. Canadian Cancer Society Website. (2024). Available online at: https://cancer.ca/en/cancer-information/cancer-types/cervical/what-is-cervical-cancer/precancerous-conditions (Accessed July 2, 2025).

Google Scholar

11. National Cancer Institute. HPV and Cancer. (2023). Available online at: https://www.cancer.gov/about-cancer/causes-prevention/risk/infectious-agents/hpv-and-cancer (Accessed July 2, 2025).

Google Scholar

12. Fennie M, Toefy Y, Sewram V. Barriers to cervical cancer screening in Africa: a systematic review. BMC Public Health. (2024) 24:525. doi: 10.1186/s12889-024-17842-1

PubMed Abstract | Crossref Full Text | Google Scholar

13. World Health Organization. Cervical cancer. (2024). Available online at: https://www.who.int/news-room/fact-sheets/detail/cervical-cancer (Accessed October 10, 2025).

Google Scholar

14. International Agency for Research on Cancer (IARC). GLOBOCAN 2022: Africa – Fact Sheet (2022). Available online at: https://gco.iarc.who.int/media/globocan/factsheets/populations/903-africa-fact-sheet.pdf (Accessed July 2, 2025).

Google Scholar

15. Lee SL, Tameru AM. A mathematical model of human papillomavirus (HPV) in the United States and its impact on cervical cancer. J Cancer. (2012) 3:62–268. doi: 10.7150/jca.4161

PubMed Abstract | Crossref Full Text | Google Scholar

16. Zhang K, Wang X, Liu H, Ji Y, Pan Q, Wei Y, et al. Mathematical analysis of a human papillomavirus transmission model with vaccination and screening. Mathem Biosci Eng. (2020) 17:5449–76. doi: 10.3934/mbe.2020294

PubMed Abstract | Crossref Full Text | Google Scholar

17. Ziyadi N, A. male-female mathematical model of human papillomavirus (HPV) in African American population. Mathem Biosci Eng. (2017) 14:339–58. doi: 10.3934/mbe.2017022

Crossref Full Text | Google Scholar

18. Ren H, Xu R, Zhang J. HPV transmission and optimal control of cervical cancer in China. Sci Rep. (2025) 15:21354. doi: 10.1038/s41598-025-05514-y

PubMed Abstract | Crossref Full Text | Google Scholar

19. Omame A, Nnanna CU, Inyama SC. Optimal control and cost-effectiveness analysis of an HPV-chlamydia trachomatis co-infection model. Acta Biotheoret. (2021) 69:185–223. doi: 10.1007/s10441-020-09401-z

PubMed Abstract | Crossref Full Text | Google Scholar

20. Lankoande J, Sawadogo WO, Kiemtoré A. Mathematical modeling and numerical simulation of the dynamics of human papillomavirus (hpv) and cervical cancer in burkina faso. Asia Pac J Math. (2025)12:1–32. doi: 10.28924/APJM/12-8

Crossref Full Text | Google Scholar

21. Arunkumar M, Rajan PK, Murugesan K. Mathematical analysis of an optimal control problem for mitigating HPV transmission and cervical cancer progression through educational campaigns. Iranian J Sci. (2025) 49:1305–1325. doi: 10.1007/s40995-025-01804-2

Crossref Full Text | Google Scholar

22. Saldaña F. Vaccination strategies in a pair formation model for human papillomavirus infection: an optimal control approach. J Theor Biol. (2025) 597:111994. doi: 10.1016/j.jtbi.2024.111994

PubMed Abstract | Crossref Full Text | Google Scholar

23. Miranda PM, Silva NNT, Pitol BCV, Silva IDCG, Lima-Filho JL, Carvalho RF, et al. Persistence and clearance of human papillomavirus infection and the role of cofactors. BioMed Res Int. (2013) 2013:782563. doi: 10.1155/2013/578276

Crossref Full Text | Google Scholar

24. Rajan PK, Kuppusamy M, Egbelowo OF. A mathematical model for Human Papillomavirus and its impact on cervical cancer in India. J Appl Mathem Comput. (2023) 69:853–770. doi: 10.1007/s12190-022-01767-2

Crossref Full Text | Google Scholar

25. Sexually Transmitted Infections Education Foundation. HPV: Overview, Epidemiology and Transmission (2024). Available online at: https://guidelines.stief.org.nz/hpv-overview-epidemiology-transmission (Accessed August 9, 2025).

Google Scholar

26. World Population Review. Which Country Has the Largest LGBTQI+ Population? (2025) Available online at: https://worldpopulationreview.com/country-rankings/lgbtq-population-by-country (Accessed August 9, 2025).

Google Scholar

27. World Health Organization. Public Health Special Edition 2020. Public Health. (2020) Special Edition.

Google Scholar

28. Office for National Statistics. Sexual orientation, UK: 2023. UK: Office for National Statistics (2023). Available online at: https://www.ons.gov.uk/peoplepopulationandcommunity/culturalidentity/sexuality/bulletins/sexualidentityuk/2023 (Accessed August 9, 2025).

Google Scholar

29. Centers for Disease Control and Prevention (CDC). HPV and Cancer (2021). Available online at: https://www.cdc.gov/hpv/parents/cancer.html (Accessed May 14, 2025).

Google Scholar

30. World Health Organization (WHO). Human papillomavirus (HPV) and cervical cancer (2017). Available online at: https://www.who.int/news-room/fact-sheets/detail/human-papillomavirus-(hpv)-and-cervical-cancer (Accessed May 14, 2025).

Google Scholar

31. Jones JH. Notes on R₀. Department of Anthropological Sciences, Stanford University. (2007).

Google Scholar

32. Yaseen RM, Ali NF, Mohsen AA, Khan A, Abdeljawad T. The modeling and mathematical analysis of the fractional-order of Cholera disease: dynamical and simulation. Part Differ Equat Appl Mathem. (2024) 12:100978. doi: 10.1016/j.padiff.2024.100978

Crossref Full Text | Google Scholar

33. Aronson JK, Brassey J, Mahtani KR. “When will it be over?”: an introduction to viral reproduction numbers, R₀ and R_e (2020). Available online at: https://www.cebm.net/covid-19/when-will-it-be-over-an-introduction-to-viral-reproduction-numbers-r0-and-re/ (Accessed May 15, 2025).

Google Scholar

34. in Biology Group M. Dynamical systems analysis of epidemiological models (2020). Available online at: https://modelinginbiology.github.io/Dynamical-System-Analysis-of-Epidemiological-Models (Accessed May 15, 2025).

Google Scholar

35. Amanso O, Aguegboh NS, Achimugwu PU, Okeke CA, Chukwuemeka BE, Nnamaga KC, et al. Mathematical model of the early incidence and spread of COVID-19 in Nigeria combined with control measure. Int J Sci Eng Res. (2020) 11:1110.

Google Scholar

36. Naim M, Zeb A, Mohsen AA, Sabbar Y, Yildiz M. Local and global stability of a fractional viral infection model with two routes of propagation, cure rate and non-lytic humoral immunity. Mathem Model Numer Simul Applic. (2024) 4:94–115. doi: 10.53391/mmnsa.1517325

Crossref Full Text | Google Scholar

37. Peet MM. Lecture 10: Routh-Hurwitz Stability Criterion (n.d.). Available online at: https://control.asu.edu/Classes/MAE318/318Lecture10.pdf (Accessed May 15, 2025).

Google Scholar

38. Healthline. HPV Types: What to Know About Diagnosis, Outlook, and Prevention (n.d.). Available online at: https://www.healthline.com/health/sexually-transmitted-diseases/hpv-types (Accessed May 15, 2025).

Google Scholar

39. Bruni L, Albero G, Serrano B, Mena M, Collado JJ, Gómez D, et al. Human Papillomavirus and Related Diseases in Nigeria: Summary Report 10 March 2023 (2023). Available online at: http://www.hpvcentre.net (Accessed May 15, 2025).

Google Scholar

40. Nigeria Health Watch. Protecting girls: Nigeria's successful HPV vaccine rollout (2023). Available online at: https://nigeriahealthwatch.medium.com/protecting-girls-nigerias-successful-hpv-vaccine-rollout-d8a4a2795c9a (Accessed May 15, 2025).

Google Scholar

41. Akarolo-Anthony SN, Famooto AO, Dareng EO, Olaniyan OB, Offiong R, Wheeler CM, et al. Age-specific prevalence of human papilloma virus infection among Nigerian women. BMC Public Health. (2014) 14:656. doi: 10.1186/1471-2458-14-656

PubMed Abstract | Crossref Full Text | Google Scholar

42. Ohihoin AG, Okwuraiwe PA, Musa AZ, Olorunfemi G, Onwuamah CK, Ige F, et al. Prevalence and predictors of high-risk HPV in Nigeria. Adv Inf Dis. (2022) 12:745–57. doi: 10.4236/aid.2022.124052

Crossref Full Text | Google Scholar

43. Kasiulevicius V, Šapoka V, Filipaviciute R. Sample size calculation in epidemiological studies. Gerontologija. (2006) 7:225–31.

Google Scholar

44. Johns Hopkins Medicine. HPV: 5 things all women should know (n.d.). Available online at: https://www.hopkinsmedicine.org/health/conditions-and-diseases/human-papillomavirus-hpv/hpv-5-things-all-women-should-know (Accessed May 15, 2025).

Google Scholar

45. Gavi the Vaccine Alliance. A leader's promise: Protecting Nigeria's girls from cervical cancer. (2023). Available online at: https://www.gavi.org/vaccineswork/leaders-promise-protecting-nigerias-girls-cervical-cancer (Accessed May 15, 2025).

Google Scholar

46. Lawson O, Ameyan L, Tukur Z, Dunu S, Kerry M, Okuyemi OO, et al. Cervical cancer screening outcomes in public health facilities in three states in Nigeria. BMC Public Health. (2023) 23:1688. doi: 10.1186/s12889-023-16539-1

PubMed Abstract | Crossref Full Text | Google Scholar

47. United Nations. Nigeria to vaccinate 7.7 million girls against leading cause of cervical cancer (2023). Available online at: https://www.un.org/africarenewal/magazine/october-2023/nigeria-vaccinate-77-million-girls-against-leading-cause-cervical-cancer (Accessed May 15, 2025).

Google Scholar

48. Bolaji A. 70% and counting: how Lagos is fending off misinformation to get more girls protected from cancer (2025). Available online at: https://www.gavi.org/vaccineswork/70-and-counting-how-lagos-fending-misinformation-get-more-girls-protected-cancer (Accessed October 10, 2025).

Google Scholar

49. Huber J, Mueller A, Siler M, Regidor PA. Human Papillomavirus persistence or clearance after infection in reproductive age. What is the status? Review of the literature and new data of a vaginal gel containing silicate dioxide, citric acid, and selenite. Women's Health. (2021) 17:17455065211020702. doi: 10.1177/17455065211020702

PubMed Abstract | Crossref Full Text | Google Scholar

50. Castle PE, Rodríguez AC, Burk RD, Herrero R, Wacholder S, Hildesheim A, et al. Long-term persistence of prevalently detected human papillomavirus infections in the absence of detectable cervical precancer and cancer. J Infect Dis. (2011) 203:814–22. doi: 10.1093/infdis/jiq116

PubMed Abstract | Crossref Full Text | Google Scholar

51. Ogbolu MO, Kozlovszky M. Assessment of HPV knowledge and awareness among students and staff at IBB University, Niger State, Nigeria: implications for health education and prevention. Healthcare. (2024) 12:665. doi: 10.3390/healthcare12060665

PubMed Abstract | Crossref Full Text | Google Scholar

52. Akanbi OA, Iyanda A, Osundare FA, Opaleye OO. Perceptions of Nigerian women about human papilloma virus, cervical cancer, and HPV vaccine. J Cancer Educ. (2015) 30:678–83. doi: 10.1155/2015/285702

PubMed Abstract | Crossref Full Text | Google Scholar

53. Statista. Nigeria: Country dossier. (2025). Cites World Bank data; p. 14, “Death rate in Nigeria... 12.43 per 1,000 inhabitants”. Statista (PDF country dossier). Available online at: https://eksportas.inovacijuagentura.lt/upload/files/2025/04/29/study-id25738-nigeria-statista-dossier.pdf (Accessed June 15, 2025).

Google Scholar

54. Clinton Health Access Initiative. How CHAI and Unitaid Are Helping Countries Like Nigeria Build a Future Free of Cervical Cancer (2022). Available online at: https://www.clintonhealthaccess.org/blog/how-chai-and-unitaid-are-helping-countries-like-nigeria-build-a-future-free-of-cervical-cancer (Accessed May 15, 2025).

Google Scholar

55. Insinga RP, Dasbach EJ, Elbasha EH. Epidemiologic natural history and clinical management of human papillomavirus (HPV) disease: a critical and systematic review of the literature in 2003–2005. Infect Agents Cancer. (2007) 9:119. doi: 10.1186/1471-2334-9-119

Crossref Full Text | Google Scholar

56. Moscicki AB, Ma Y, Wibbelsman C, Darragh TM, Powers A, Farhat S, et al. Rate of and risks for regression of cervical intraepithelial neoplasia 2 in adolescents and young women. Obstetr Gynecol. (2004) 104:243–50. doi: 10.1097/AOG.0b013e3181fe777f

PubMed Abstract | Crossref Full Text | Google Scholar

57. De Vincenzo R, Conte C, Ricci C, Caruso S, Scambia G, Capelli G. Low-grade cervical intraepithelial neoplasia (CIN1) evolution: analysis of opportunistic HPV vaccination impact. Vaccines. (2023) 11:284. doi: 10.3390/vaccines11020284

Crossref Full Text | Google Scholar

58. Baussano I, Lazzarato F, Ronco G, Lehtinen M, Dillner J, Franceschi S. Different challenges in eliminating hpv16 compared to other types: a modeling study. J Infect Dis. (2017) 216:336–44. doi: 10.1093/infdis/jix299

PubMed Abstract | Crossref Full Text | Google Scholar

59. Li Y, Wang H, Zhang J. A two-sex model of human papillomavirus infection. medRxiv. (2021).

Google Scholar

60. Abuduwaili S, Wang L, Teng Z, Ailawaer A, Rifhat R. Estimating real-time reproduction number for hpv infection in xinjiang, china. Adv Contin Discrete Models. (2025) 2025:33. doi: 10.1186/s13662-025-03896-x

Crossref Full Text | Google Scholar

61. United United Nations Department of Economic and Social Affairs Population Division. Manual III: Methods for Estimating Populations: Chapter III – Stable Population Model (2003). Available online at: https://www.un.org/en/development/desa/population/publications/pdf/manuals/model/stablepopulation/chap3.pdf (Accessed May 15, 2025).

Google Scholar

62. Statista Research Department. Nigeria: Life expectancy at birth by gender 2022 (2024). Available online at: https://www.statista.com/statistics/1122851/life-expectancy-in-nigeria-by-gender/ (Accessed May 15, 2025).

Google Scholar

63. Chitnis N, Hyman JM, Cushing JM. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull Math Biol. (2008) 70:1272–96. doi: 10.1007/s11538-008-9299-0

PubMed Abstract | Crossref Full Text | Google Scholar

64. Rodrigues HS, Monteiro MTT, Torres DFM. Sensitivity analysis in a dengue epidemiological model. In: Conference Papers in Mathematics. (2013). p. 721406. doi: 10.1155/2013/721406

Crossref Full Text | Google Scholar

65. Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global Sensitivity Analysis: The Primer. Chichester, UK: John Wiley & Sons. (2008). doi: 10.1002/9780470725184

Crossref Full Text | Google Scholar

66. Blower SM, Dowlatabadi H. Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int Stat Rev. (1994) 62:229–43. doi: 10.2307/1403510

Crossref Full Text | Google Scholar

67. Okongo W, Okelo Abonyo J, Kioi D, Moore SE, Aguegboh SN. Mathematical modeling and optimal control analysis of monkeypox virus in contaminated environment. Model Earth Syst Environ. (2024) 10:3969–94. doi: 10.1007/s40808-024-01987-4

Crossref Full Text | Google Scholar

68. Aguegboh NS, Okongo W, Boubacar D, Munkaila D, Nnamaga KC, Nnaji DU, et al. A novel approach to modeling malaria with treatment and vaccination as control strategies in Africa using the Atangana-Baleanu derivative. Model Earth Syst Environ. (2025) 11:110. doi: 10.1007/s40808-024-02273-z

Crossref Full Text | Google Scholar

69. Mohsen AA, Naji RK. Dynamical analysis within-host and between-host for HIV/AIDS with the application of optimal control strategy. Iraqi J Sci. (2020) 61:1173–89. doi: 10.24996/ijs.2020.61.5.25

Crossref Full Text | Google Scholar

70. Kirk DE. Optimal Control Theory: An Introduction. North Chelmsford: Courier Corporation. (2004).

Google Scholar

71. Press WH Teukolsky SA Vetterling WT Flannery BP. Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge: Cambridge University Press. (2007).

Google Scholar

72. Butcher JC. Numerical Methods for Ordinary Differential Equations. New York: 2nd ed John Wiley & Sons. (2016). doi: 10.1002/9781119121534

Crossref Full Text | Google Scholar

73. Iskandar D, Ong CT. The application of Runge-Kutta fourth order method in SIR model for simulation of COVID-19 cases. In: Proceedings of Science and Mathematics. vol. 10. UTM Science (2022). p. 61–70.

Google Scholar

74. Desta HD, Eshetu D, Damtie M. Mathematical model of human papillomavirus (HPV) dynamics with double-dose vaccination and its impact on cervical cancer. Discr Dyn Nat Soc. (2024) 2024:9971859. doi: 10.1155/ddns/9971859

Crossref Full Text | Google Scholar

75. Okware T, Apima J, Wanjara G, Mutwiwa J. Mathematical modelling of human papillomavirus dynamics with vaccination incorporating optimal control analysis. Asian Res J Mathem. (2023) 19:36–51. doi: 10.9734/arjom/2023/v19i11751

Crossref Full Text | Google Scholar

76. National Cancer Institute. Cervical Cancer Prevention (PDQ^®)–Health Professional Version (n.d.). Available online at: https://www.cancer.gov/types/cervical/hp/cervical-prevention-pdq (Accessed May 15, 2025).

Google Scholar