Further evidence: Mathematical Assessment of the Roles of Vaccination and Pap Screening on the Burden of HPV and Related Cancers in Korea

Daniel Rossi; Sigrid Bergström

doi:10.58301/y93nc380

Home Journals Meridian Journal of Mathematical SciencesVolume 11, Issue 410.58301/y93nc380

Meridian Journal of Mathematical Sciences

Article7 November 2025

Further evidence: Mathematical Assessment of the Roles of Vaccination and Pap Screening on the Burden of HPV and Related Cancers in Korea

Daniel Rossi^1* (unauthenticated), Sigrid Bergström²

¹Department of Mathematics, Sapienza Mediterranea, Italy

²Faculty of Education, Anatolian Technical University, Turkey

* Correspondence: daniel.rossi@mediterranea.example

Meridian J. Math. Sci. 2025, 11(4), 319-327; https://doi.org/10.58301/y93nc380

Abstract

This study is based on using a novel sex-structured mathematical model to assess the effectiveness of vaccination and Pap screening against HPV and related cancers in South Korea. In addition to its disease-free equilibrium (DFE) being locally-asymptotically stable when the associated control reproduction number is less than one, the model could have one or three endemic equilibria, for a special case with negligible disease-induced mortality, if the reproduction number exceeds one. It’s shown, using a Krasnoselskii sublinearity argument, that this special case has a unique and locally-asymptotically stable endemic equilibrium, when the reproduction number is larger than one, if, additionally, the HPV vaccine is assumed to be perfect. The DFE of a simplified version of the model, which is calibrated using HPV-related cancer data in Korea, is globally-asymptotically stable when its reproduction number is less than one. Simulations of the full model showed that, although vaccine-derived herd immunity (needed for HPV elimination) cannot be achieved in Korea under the current vaccination coverage of females (of 88%), it can be achieved if, additionally, at least 65% of males are vaccinated at steady-state. While the current combined vaccination-screening strategy (termed Strategy A) will fail to eliminate HPV, extended strategies that include increased coverage of female vaccination (termed Strategy B) or additionally vaccinating boys (termed Strategy C) could lead to such elimination in Korea. The implementation of boys-only vaccination strategy induces a significant spillover benefit in reducing cervical cancer burden, which exceeds the corresponding spillover benefit achieved by implementing a girls-only vaccination strategy.

Keywords:Keywords HPVCervical cancerVaccinationKoreaReproduction numberAsymptotic stability

Introduction

While cervical cancer is the fourth most frequent female cancer globally (there were 661,021 new cervical cancer cases and 348,189 fatalities globally in the year 2022 alone (Bray et al. 2024)), it is the third most frequent in the 15-44-year female age group in South Korea (there were 3,218 new cervical cancer cases and 1,018 fatalities in South Korea in the year 2020 (Bruni et al. 2023)). Cervical cancer is one of the top three cancers in females under the age of 45 globally, and the average ages for cervical cancer diagnosis and death are 53 years and 59 years, respectively (Arbyn et al. 2020). The burden of the cervical cancer disproportionately affects resource-limited countries that have lack of availability and access to screening, vaccination, and treatment measures (Bray et al. 2024).

Virtually all cervical cancer cases are caused by HPV(Human Papillomavirus) (Walboomers et al. 1999), which is transmitted primarily through skin-to-skin or mucosa-to-mucosa sexual contacts (Schiffman et al. 2007). Among the more than 200 serotypes of HPV, the most prevalent serotypes are HPV-16 and HPV-18, which are high-risk oncogenic types that are found in 70% of cervical cancer cases (Smith et al. 2007). The oncogenic types HPV-31/33/35/45/52/58 are the next most common serotypes that account for additional 20% of cervical cancers globally (Clifford et al. 2006). These high-risk oncogenic types cause the most pre-cancerous lesions and cancers in humans. On the other hand, HPV-6 and HPV-11 are low-risk non-oncogenic types that are responsible for 90% of genital warts (Walboomers et al. 1999). Although 70-90% of sexually active adults contract HPV at least once during their lifetime (Bosch and de Sanjose 2003; Chesson et al. 2014), 90% of HPV infections clear naturally within 2 years (Schiffman et al. 2007). HPV persists in the remaining 10%, increasing the risk of developing pre-cancerous lesions and cancer in these individuals (Bulkmans et al. 2007; Schiffman et al. 2007). Empirical studies in the US and Netherlands show that males recover faster from HPV infection than females (Giuliano et al. 2008a; Moscicki et al. 2004, 1998), and persistent infection is more common in females than in males (Van Doornum et al. 1994). This explains, in part, the higher HPV prevalence and number of HPV-related cancer cases in females than in males in South Korea (National Cancer Center Korea Central Cancer Registry 2022; Shin et al. 2004). The pre-cancer stage of the natural history of HPV is referred to as intraepithelial neoplasia. The cervical intraepithelial neoplasia (CIN) stage is divided into three grades, namely CIN1, CIN2, and CIN3, depending on the extent of the abnormalities (IARC 1995), and the risk of developing invasive cervical cancer is higher for higher grades of CIN (IARC 1995; Schiffman et al. 2007). Similarly, intraepithelial neoplasia in males (INM) is divided into three stages (INM1, INM2, and INM3), and leads to various cancers in males (such as anal and penile cancers) (IARC 1995).

Although the risk posed by the cervical cancer has been the major driving force of the studies on the HPV dynamics, HPV is also responsible for 88% of anal cancer, 50% of penile cancer, and 30.8% of oropharyngeal cancer cases globally, exerting burden on both females and males (de Martel et al. 2017). In South Korea, HPV-related male cancer cases tripled in the past 20 years, reaching 573 new cases in 2020 (National Cancer Center Korea Central Cancer Registry 2022). An empirical study in the US, Mexico, and Brazil reported that 65.2% of asymptomatic males in the age group of 18-70 have HPV DNA in the genital area (Giuliano et al. 2008b). Major risk factors for HPV infection include anal intercourse and contraction of other sexually-transmitted infections (STIs) (Chin-Hong Peter et al. 2004; Hernandez et al. 2008; Kyo et al. 1994; Lee et al. 2016; Marra et al. 2018), hence men who have sex with men (MSM) is one of the main risk groups of HPV infection among sexually active individuals (Chin-Hong Peter et al. 2004; Lee et al. 2016; Marra et al. 2018).

Three HPV vaccines have been approved for use in females of age 9-45 and males of age 9-26 in South Korea; the bivalent (Cervarix^®, preventing HPV 16/18), quadrivalent (Gardasil^®, preventing HPV 16/18/6/11), and 9-valent (Gardasil9^®, preventing HPV 16/18/6/11/31/33/45/52/58) vaccines, which were approved in 2008, 2007, and 2016, respectively. The 9-valent vaccine (Gardasil9^®) provides the broadest coverage of HPV types, covering 9 HPV types, and costs about $287.54 per dose in the United States in 2024 (CDC 2024) (it costs about $150 per dose in South Korea (Choi et al. 2022)). All three vaccines protect against HPV infection with over 90% efficacy if used prior to HPV exposure, and offer protection against precancerous lesions (Cheng et al. 2020). People of age 9-14 are recommended to get 2 doses, while those of age 15-26 are recommended to get 3 doses of the vaccine (World Health Organization 2022). In South Korea, the HPV National Immunization Program (NIP) began in June 2016, offering 2 doses of bivalent(Cervarix^®) or quadrivalent(Gardasil^®) vaccines to girls of age 12. The routine HPV immunization program in Korea was expanded to cover girls of ages 12-17 in 2022. In 2016, before the onset of the NIP, 23.1% of Korean females aged 9-59 years received the first dose of HPV vaccine (Choi et al. 2022). In 2018 (two years after the implementation of the NIP), however, 82.3% of 12-year-old girls received the first dose of HPV vaccine in Korea (Kim et al. 2018). This is still relatively low compared to the 91-98% coverage of infant vaccination series against childhood diseases covered by the NIP in South Korea (KCDC Center for Infectious Disease Control 2023).

Papanicolaou test (Pap test) is widely used as the primary method of screening for cervical cancer in countries with organized screening programs, allowing for early detection and treatment of pre-cancerous lesions, thereby reducing the incidence and mortality of cervical cancer (Devesa et al. 1989; Lim and Yoo 2019; Min et al. 2015). Some empirical studies suggest the utility of using HPV DNA test (such as, Hybrid Capture 2 and PCR-based methods) as an alternative or adjunctive screening method to the Pap test (Reid et al. 1991; Denny and Wright 2005). However, South Korea does not recommend it as a standalone test, due to its low specificity, high cost, and high variability of its performance (Lim and Yoo 2019). In contrast, Pap test (conventional Pap smear or liquid-based cytology) is implemented in South Korea, with the Korean Society for Cytopathology entrusted with managing the quality control of its nationwide implementation (Lim and Yoo 2019; Lee et al. 2008). South Korea started offering free Pap screening, as a part of the National Cancer Screening Program (NCSP), to females over age 30 every 2 years in 2002, and then extended to females over age 20 since 2016 (Shin et al. 2022). The national cervical cancer screening guideline in South Korea (2015) suggests that females over age 20 receive Pap test every 3 years, with the option of combining Pap test with HPV test in consideration of clinical decision for individual risk and preference (Min et al. 2015). The guideline also states that the screening can terminate at the age of 74 if more than 3 consecutive cytology results are negative within 10 years (Min et al. 2015). The Korean National Cancer Screening Survey reported that the Pap screening coverage in Korea during the period 2005-2020 was 61% on average (Shin et al. 2022).

Numerous mathematical models, typically in the form of a deterministic system of nonlinear differential equations, have been developed and used to study the natural history of HPV (Elbasha and Galvani 2005; Gao et al. 2021) and its progression to associated cancers (Alsaleh and Gumel 2014; Choi and Shim 2024; Elbasha et al. 2007; Goldie et al. 2007, 2008; Javame and Gumel 2016; Malik et al. 2013; Milwid et al. 2018). For instance, Elbasha and Galvani (Elbasha and Galvani 2005) examined the impact of vaccination on the dynamics of two types of HPV in a community. (Gao et al. 2021) developed and analyzed a model for HPV dynamics that includes vaccination in both the heterosexual and the MSM populations. (Alsaleh and Gumel 2014) developed and analyzed a 29-dimensional risk-structured model for HPV dynamics, computing optimal vaccination coverage level. Some mathematical models were also developed and used to investigate the cost-effectiveness of different HPV vaccination strategies (Choi and Shim 2024; Cody et al. 2021; Elbasha et al. 2007; Goldie et al. 2007, 2008). Furthermore, while (Javame and Gumel 2016) investigated the effects of Pap screening on the prevalence of HPV and related cancers, other studies (such as those in (Malik et al. 2013; Milwid et al. 2018)) examined the effects of both the vaccination and Pap screening programs against HPV and related cancers in a community.

The purpose of the current study is to extend some of the previous modeling studies to quantify the combined population-level impacts of vaccination and Pap screening against HPV and HPV-related cancers in South Korea. Specifically, the main objective of this study is to assess whether or not the existing control resources (i.e., vaccination and Pap screening, under current levels of coverage and efficacy) will be sufficient to significantly curtail the incidence and burden of HPV and HPV-related cancers in the main high-risk groups in South Korea. The effectiveness of various control strategies will also be theoretically assessed. The paper is organized as follows. The vaccination and Pap screening model for the transmission of HPV is formulated in Section 2. Detailed qualitative analysis (with respect to the existence and asymptotic stability of its disease-free and endemic equilibria) of the model are carried out in this section. In order to calibrate the model and estimate its unknown parameters, a simplified version of the model without Pap screening and with the three CIN and INM stages lumped into one compartment, is considered in Section 3. The calibration and qualitative analysis of this simplified version of the model are also described and presented in this section. Detailed global sensitivity analysis for the full model with vaccination and Pap screening (based on the values of the unknown parameters estimated from the fitting of the reduced model) is carried out in Section 4 to determine the most influential parameters of the full vaccination and Pap screening model. Numerical simulations of this model are carried out in Section 5. Concluding remarks and discussion are presented in Section 6.

Formulation of Vaccination and Pap Screening Model

In this section, a mechanistic model for the heterosexual transmission of HPV will be developed and used to assess the combined population-level impact of routine HPV vaccination and Pap cytology screening on the incidence and burden of HPV and HPV-related cancers in South Korea. In addition to incorporating the life history features of HPV and associated dysplasia, the model to be developed stratifies the total sexually active population (which is assumed to be of age 17 or older (Nakazawa et al. 2023; Shin et al. 2004)) based on gender (male or female). Specifically, the total sexually active population of South Korea at time t, denoted by N(t), is split into the total number of females at time t (denoted by

N_{f} (t)

), and the total number of males at time t (denoted by

N_{m} (t)

), so that

N (t) = N_{f} (t) + N_{m} (t)

. Furthermore, the population of sexually active females is stratified into the 15 mutually exclusive compartments of unvaccinated susceptible (

S_{f} (t)

), vaccinated susceptible (

V_{f} (t)

), exposed/pre-symptomatic infectious (

E_{f} (t)

), symptomatic infectious (

I_{f} (t)

), persistently infected (

P_{f} (t)

) females, females having undetected CIN (

Q_{iu}

, with

i = 1, 2, 3

accounting for the three CIN stages or grades), detected CIN (

Q_{id}, i = 1, 2, 3

), undetected cervical cancer (

C_{u} (t)

), detected cervical cancer (

C_{d} (t)

), females who recovered from cervical cancer (

R_{f}^{C} (t)

), and females who recovered from HPV infection before developing cervical cancer (

R_{f} (t)

). Hence,

\begin{matrix} N_{f} (t) = & S_{f} (t) + V_{f} (t) + E_{f} (t) + I_{f} (t) + P_{f} (t) \\ + \sum_{i = 1}^{3} (Q_{iu} (t) + Q_{id} (t)) + C_{u} (t) + C_{d} (t) + R_{f}^{C} (t) + R_{f} (t) . \end{matrix}

Similarly, the total population of sexually active males is split into the 11 mutually exclusive compartments of unvaccinated susceptible (

S_{m} (t)

), vaccinated susceptible (

V_{m} (t)

), exposed/pre-symptomatic infectious (

E_{m} (t)

), symptomatic infectious (

I_{m} (t)

), persistently infected (

P_{m} (t)

) males, males having INM (

J_{im},

with

i = 1, 2, 3

accounting for the three INM stages), males having HPV-related cancers (

C_{m}

), males who recovered from HPV-related cancers (

R_{m}^{C}

), and males who recovered from HPV infection before developing HPV-related cancers (

R_{m}

), so that

\begin{matrix} N_{m} (t) = S_{m} (t) + V_{m} (t) + E_{m} (t) + I_{m} (t) + P_{m} (t) + \sum_{i = 1}^{3} J_{im} (t) + C_{m} (t) + R_{m}^{C} (t) + R_{m} (t) . \end{matrix}

The population of sexually active unvaccinated susceptible females (

S_{f}

) and males (

S_{m}

) acquire HPV infection at the rate

λ_{mf}

and

λ_{fm}

, respectively, given by:

\begin{matrix} λ_{mf} & = (β_{mf}, c_{mf}) [\frac{I_{m} + η_{m} E_{m} + θ_{m} [P_{m} + J_{1 m} + κ_{m} (J_{2 m} + J_{3 m})]}{N_{m}}], \\ λ_{fm} & = (β_{fm}, c_{fm}) [\frac{I_{f} + η_{f} E_{f} + θ_{f} [P_{f} + Q_{1 u} + κ_{f} (Q_{2 u} + Q_{3 u}) + ν {Q_{1 d} + κ_{fd} (Q_{2 d} + Q_{3 d})}]}{N_{f}}] . \end{matrix}

2.1

In (2.1),

β_{mf}

is the probability of HPV transmission from symptomatic infectious males to unvaccinated susceptible females per partnership, and

c_{mf}

is the average number of male sexual partners a sexually active female has per unit time (for this study, time is measured in years). The modification parameter

0 \leq η_{m} < 1

accounts for the assumed reduced transmission probability of HPV by exposed/pre-symptomatic infectious males, in comparison to the transmission probability of symptomatic infectious males. Furthermore,

0 \leq θ_{m} < 1

measures the assumed reduction of infectiousness of males who are at, or past, the persistent HPV infection stage, in comparison to the infectiousness of symptomatic infectious males, and

κ_{m} \geq 1

is the modification parameter accounting for the assumed increased probability of HPV transmission from males with INM2 and INM3, in comparison to males with INM1 and persistent HPV infection. Similarly,

β_{fm}

is the transmission probability from symptomatic infectious females to unvaccinated susceptible males per partnership, and

c_{fm}

represents the average number of female sexual partners a male has per unit time. The modification parameter

0 \leq η_{f} < 1

accounts for the assumed reduced probability of HPV transmission by exposed/pre-symptomatic infectious females, in comparison to the HPV transmission by symptomatic infectious females. Furthermore, the parameter

0 \leq θ_{f} < 1

represents the assumed reduction of the infectiousness of females that are at, or past, the persistent HPV infection stage, in comparison to the infectiousness of symptomatic infectious females, and

κ_{f} \geq 1

accounts for the assumed increased probability of HPV transmission by females with undetected CIN2 and CIN3, in comparison to females with undetected CIN1 and persistent HPV infection. The parameter

0 \leq ν < 1

accounts for the assumed reduction in the number of sexual contacts by females detected with pre-cancerous lesions, in comparison to females with undetected pre-cancerous lesions and persistent infection. Finally, the parameter

κ_{fd} \geq 1

represents the assumed increase in the probability of HPV transmission by females with detected CIN2 and CIN3, in comparison to females with detected CIN1.

It should be stated that, since the model to be developed in this section only considers heterosexual transmission of HPV (i.e., transmission only occurs via sexual contacts between males and females), the following conservation law of total number of sexual partnerships per unit time must hold (i.e., the total number of sexual partnerships sexually active females have (with males) in the community must balance those by sexually active males (with females)):

\begin{matrix} c_{mf} N_{f} (t) = c_{fm} (t) N_{m} (t) . \end{matrix}

2.2

It is assumed that females can always find enough number of male partners in the community to have sexual contacts with, so

c_{mf}

is set to be a constant (this assumption is justified by the fact that males generally have a higher desired number of female sexual partners than females (Fenigstein and Preston 2007), and tend to overreport their lifelong number of female sexual partners (Brewer et al. 2000), suggesting a higher demand among males for multiple female sexual partners. In addition, an empirical study using the famous online dating app, Tinder, shows that women attain a larger number of men matches in a short time (Tyson et al. 2016), which supports the claim that it is easier for women to find sexual partners than men). Hence, based on all this, the parameter

c_{fm}

(for the average number of female sexual partners a male has) can be computed from (2.2) as (i.e., the rate

c_{fm}

depends on the ratio of the total sub-populations of males and females at each time)

\begin{matrix} c_{fm} (t) = \frac{c_{mf} N_{f} (t)}{N_{m} (t)} . \end{matrix}

2.3

Using (2.2) (with (2.3)), it follows that the expressions for the forces of infection,

λ_{mf}

and

λ_{fm}

, given in (2.1), can now be re-written, respectively, as:

\begin{matrix} λ_{mf} & = (β_{mf}, c_{mf}) [\frac{I_{m} + η_{m} E_{m} + θ_{m} [P_{m} + J_{1 m} + κ_{m} (J_{2 m} + J_{3 m})]}{N_{m}}], and \\ λ_{fm} & = (β_{fm}, c_{mf}) [\frac{I_{f} + η_{f} E_{f} + θ_{f} [P_{f} + Q_{1 u} + κ_{f} (Q_{2 u} + Q_{3 u}) + ν {Q_{1 d} + κ_{fd} (Q_{2 d} + Q_{3 d})}]}{N_{m}}] . \end{matrix}

2.4

Based on the above formulation, derivations (for the conservation law (2.3) and forces of infection (2.4)), and assumptions, the sex-structured model for the heterosexual transmission of HPV and related cancers, in the presence of HPV vaccination and Pap screening, in South Korea is given by the following deterministic system of nonlinear differential equations:

\begin{matrix} Females \{\begin{matrix} {\dot{S}}_{f} & = & (1 - q_{f}) π_{f} - (λ_{mf} + ξ_{f} + μ) S_{f}, \\ {\dot{V}}_{f} & = & q_{f} π_{f} + ξ_{f} S_{f} - (1 - ε_{v}) λ_{mf} V_{f} - μ V_{f}, \\ {\dot{E}}_{f} & = & λ_{mf} S_{f} + (1 - ε_{v}) λ_{mf} V_{f} - (σ_{f} + μ) E_{f}, \\ {\dot{I}}_{f} & = & σ_{f} E_{f} - (ψ_{f} + μ) I_{f}, \\ {\dot{P}}_{f} & = & (1 - b_{f}) ψ_{f} I_{f} - (α_{f} + μ) P_{f}, \\ {\dot{Q}}_{1 u} & = & (1 - k_{f}) α_{f} P_{f} + h_{u 2} y_{f} Q_{2 u} - (g_{f} + μ) Q_{1 u}, \\ {\dot{Q}}_{2 u} & = & (1 - d_{u 1} - d_{u 3}) g_{f} Q_{1 u} + q_{u 2} z_{f} Q_{3 u} - (y_{f} + μ) Q_{2 u}, \\ {\dot{Q}}_{3 u} & = & (1 - h_{u 1} - h_{u 2} - h_{u 3}) y_{f} Q_{2 u} - (z_{f} + μ) Q_{3 u}, \\ {\dot{Q}}_{1 d} & = & d_{u 3} g_{f} Q_{1 u} + h_{d 2} ζ_{2} Q_{2 d} - (ζ_{1} + μ) Q_{1 d}, \\ {\dot{Q}}_{2 d} & = & h_{u 3} y_{f} Q_{2 u} + (1 - d_{d 1}) ζ_{1} Q_{1 d} + q_{d 2} ζ_{3} Q_{3 d} - (ζ_{2} + μ) Q_{2 d}, \\ {\dot{Q}}_{3 d} & = & q_{u 3} z_{f} Q_{3 u} + (1 - h_{d 1} - h_{d 2}) ζ_{2} Q_{2 d} - (ζ_{3} + μ) Q_{3 d}, \\ {\dot{C}}_{u} & = & (1 - q_{u 1} - q_{u 2} - q_{u 3}) z_{f} Q_{3 u} - (w_{f} + r_{u} + δ + μ) C_{u}, \\ {\dot{C}}_{d} & = & (1 - q_{d 1} - q_{d 2}) ζ_{3} Q_{3 d} + w_{f} C_{u} - (r_{d} + δ_{d} + μ) C_{d}, \\ {\dot{R}}_{f}^{C} & = & r_{u} C_{u} + r_{d} C_{d} - μ R_{f}^{C}, \\ {\dot{R}}_{f} & = & b_{f} ψ_{f} I_{f} + k_{f} α_{f} P_{f} + d_{u 1} g_{f} Q_{1 u} + h_{u 1} y_{f} Q_{2 u} + q_{u 1} z_{f} Q_{3 u} + d_{d 1} ζ_{1} Q_{1 d} \\ + h_{d 1} ζ_{2} Q_{2 d} + q_{d 1} ζ_{3} Q_{3 d} - μ R_{f}, \end{matrix}) \end{matrix}

\begin{matrix} \begin{matrix} Males \{\begin{matrix} {\dot{S}}_{m} & = & (1 - q_{m}) π_{m} - (λ_{fm} + ξ_{m} + μ) S_{m}, \\ {\dot{V}}_{m} & = & q_{m} π_{m} + ξ_{m} S_{m} - (1 - ε_{v}) λ_{fm} V_{m} - μ V_{m}, \\ {\dot{E}}_{m} & = & λ_{fm} S_{m} + (1 - ε_{v}) λ_{fm} V_{m} - (σ_{m} + μ) E_{m}, \\ {\dot{I}}_{m} & = & σ_{m} E_{m} - (ψ_{m} + μ) I_{m}, \\ {\dot{P}}_{m} & = & (1 - b_{m}) ψ_{m} I_{m} - (α_{m} + μ) P_{m}, \\ {\dot{J}}_{1 m} & = & (1 - k_{m}) α_{m} P_{m} + h_{m 2} y_{m} J_{2 m} - (g_{m} + μ) J_{1 m}, \\ {\dot{J}}_{2 m} & = & (1 - d_{m}) g_{m} J_{1 m} + q_{m 2} z_{m} J_{3 m} - (y_{m} + μ) J_{2 m}, \\ {\dot{J}}_{3 m} & = & (1 - h_{m 1} - h_{m 2}) y_{m} J_{2 m} - (z_{m} + μ) J_{3 m}, \\ {\dot{C}}_{m} & = & (1 - q_{m 1} - q_{m 2}) z_{m} J_{3 m} - (r_{m} + μ) C_{m}, \\ {\dot{R}}_{m}^{C} & = & r_{m} C_{m} - μ R_{m}^{C}, \\ {\dot{R}}_{m} & = & b_{m} ψ_{m} I_{m} + k_{m} α_{m} P_{m} + d_{m} g_{m} J_{1 m} + h_{m 1} y_{m} J_{2 m} + q_{m 1} z_{m} J_{3 m} - μ R_{m} . \end{matrix}) \end{matrix} \end{matrix}

2.5

In the model (2.5), the population of sexually active unvaccinated susceptible females (

S_{f}

) is generated by the recruitment of new sexually active individuals into the community at a rate

(1 - q_{f}) π_{f}

(where the proportion

q_{f}

is vaccinated against HPV), and is reduced by the acquisition of HPV infection at the rate

λ_{mf}

, administration of the catch-up HPV vaccination (at a rate

ξ_{f}

), and natural death (at a rate

μ

; individuals in all epidemiological class are assumed to die naturally at this rate). The population of sexually active vaccinated susceptible females (

V_{f}

) is increased by the cohort (at the rate

q_{f} π_{f}

) and catch-up (at the rate

ξ_{f}

) vaccination. Breakthrough infection occurs at the rate

(1 - ε_{v}) λ_{mf}

, where

0 \leq ε_{v} \leq 1

is the efficacy of the vaccine to protect against the acquisition of HPV infection. Exposed/pre-symptomatic infectious females (i.e., those in the

E_{f}

class) develop clinical symptoms of HPV at a rate

σ_{f}

. Symptomatic infectious females (i.e., those in the

I_{f}

class) leave this class at a rate

ψ_{f}

, a proportion,

b_{f}

, of these individuals recover (at a rate

b_{f} ψ_{f}

), and the remaining proportion,

1 - b_{f}

, progress to persistent HPV infection (at a rate

(1 - b_{f}) ψ_{f}

). Females with persistent HPV infection (i.e., those in the

P_{f}

class) leave this class at a rate

α_{f}

, a proportion,

k_{f}

, of these individuals recover (at a rate

k_{f} α_{f}

), and the remaining proportion,

1 - k_{f}

, move to the undetected CIN1 class (at a rate

(1 - k_{f}) α_{f}

).

Females with undetected CIN1 (i.e., those in the $Q_{1 u}$ class) transition out of this class at a rate $g_{f}$ , and a proportion, $d_{u 1}$ , of these individuals recover (at a rate $d_{u 1} g_{f}$ ), another proportion, $d_{u 3}$ , are detected via Pap screening and move to the detected CIN1 class (at a rate $d_{u 3} g_{f}$ ), and, finally, the remaining proportion, $1 - (d_{u 1} + d_{u 3})$ , progress to the undetected CIN2 class (at a rate $(1 - d_{u 1} - d_{u 3}) g_{f}$ ). Females with undetected CIN2 (i.e., those in the $Q_{2 u}$ class) leave this class at a rate $y_{f}$ . A proportion, $h_{u 1}$ , of these individuals recover (at a rate $h_{u 1} y_{f}$ ), another proportion, $h_{u 2}$ , regress to the undetected CIN1 class (at a rate $h_{u 2} y_{f}$ ), a proportion, $h_{u 3}$ , are detected via Pap screening and move to the detected CIN2 class (at a rate $h_{u 3} y_{f}$ ), and the remaining proportion, $1 - (h_{u 1} + h_{u 2} + h_{u 3})$ , progress to the undetected CIN3 class (at a rate $(1 - h_{u 1} - h_{u 2} - h_{u 3}) y_{f}$ ). Females with undetected CIN3 (i.e., those in the $Q_{3 u}$ class) leave this class at a rate $z_{f}$ . A proportion, $q_{u 1}$ , of these individuals recover (at a rate $q_{u 1} z_{f}$ ), a proportion, $q_{u 2}$ , regress to the undetected CIN2 class (at a rate $q_{u 2} z_{f}$ ), a proportion, $q_{u 3}$ , are detected via Pap screening and move to the detected CIN3 class (at a rate $q_{u 3} z_{f}$ ), and the remaining proportion, $1 - (q_{u 1} + q_{u 2} + q_{u 3})$ , progress to the undetected cervical cancer class (at a rate $(1 - q_{u 1} - q_{u 2} - q_{u 3}) z_{f}$ ). Females with detected CIN1 (i.e., those in the $Q_{1 d}$ class) transition out of this class at a rate $η_{1}$ . A proportion, $d_{d 1}$ , of these individuals recover (at a rate $d_{d 1} η_{1}$ ), and the remaining proportion, $1 - d_{d 1}$ , progress to the detected CIN2 class (at a rate $(1 - d_{d 1}) η_{1}$ ). Females with detected CIN2 (i.e., those in the $Q_{2 d}$ class) leave this class at a rate $η_{2}$ . A proportion, $h_{d 1}$ , of these individuals recover (at a rate $h_{d 1} η_{2}$ ), a proportion, $h_{d 2}$ , regress to the detected CIN1 class (at a rate $h_{d 2} η_{2}$ ), and the remaining proportion, $1 - (h_{d 1} + h_{d 2})$ , progress to the detected CIN3 class (at a rate $(1 - h_{d 1} - h_{d 2}) η_{2}$ ). Females with detected CIN3 (i.e., those in the $Q_{3 d}$ class) leave this class at a rate $η_{3}$ . A proportion, $q_{d 1}$ , of these individuals recover (at a rate $q_{d 1} η_{3}$ ), a proportion, $q_{d 2}$ , regress to the detected CIN2 class (at a rate $q_{d 2} η_{3}$ ), and the remaining proportion, $1 - (q_{d 1} + q_{d 2})$ , progress to the detected cervical cancer class (at a rate $(1 - q_{d 1} - q_{d 2}) η_{3}$ ). Females with undetected cervical cancer (i.e., those in the $C_{u}$ class) recover from cancer at a rate $r_{u}$ , are detected via Pap screening and move to the detected cervical cancer class at a rate $w_{f}$ , or die of cancer at a rate $δ$ . Females with detected cervical cancer (i.e., those in the $C_{d}$ class) recover from cervical cancer at a rate $r_{d}$ , or die of cancer at a rate $δ_{d}$ .

The population of sexually active unvaccinated susceptible males ( $S_{m}$ ) is increased by recruitment at a rate $(1 - q_{m}) π_{m}$ (where the proportion $q_{m}$ is vaccinated against HPV), and is decreased by the acquisition of HPV infection (at the rate $λ_{fm}$ ), catch-up vaccination (at a rate $ξ_{m}$ ), and natural death (at the rate $μ$ ). The population of vaccinated susceptible males ( $V_{m}$ ) is increased by the cohort (at the rate $q_{m} π_{m}$ ) and catch-up (at the rate $ξ_{m}$ ) vaccination. Breakthrough infection occurs at the rate $(1 - ϵ_{v}) λ_{mf}$ . Exposed/pre-symptomatic infectious males (i.e., those in the $E_{m}$ class) develop clinical symptoms of HPV at a rate $σ_{m}$ . Symptomatic infectious males (i.e., those in the $I_{m}$ class) leave this class at a rate $ψ_{m}$ . A proportion, $b_{m}$ , of these individuals recover (at a rate $b_{m} ψ_{m}$ ), and the remaining proportion, $1 - b_{m}$ , progress to persistent infection (at a rate $(1 - b_{m}) ψ_{m}$ ). Males with persistent HPV infection (i.e., those in the $P_{m}$ class) leave this class at a rate $α_{m}$ . A proportion, $k_{m}$ , of these individuals recover (at a rate $k_{m} α_{m}$ ), and the remaining proportion, $1 - k_{m}$ , move to the undetected CIN1 class (at a rate $(1 - k_{m}) α_{m}$ ). Males with INM1 (i.e., those in the $J_{1 m}$ class) transition out of this class at a rate $g_{m}$ . A proportion, $d_{a}$ , of these individuals recover (at a rate $d_{m} g_{m}$ ), and the remaining proportion, $1 - d_{m}$ , progress to the INM2 class (at a rate $(1 - d_{m}) g_{m}$ ). Males with INM2 (i.e., those in the $J_{2 m}$ class) leave this class at a rate $y_{m}$ . A proportion, $h_{m 1}$ , of these individuals recover (at a rate $h_{m 1} y_{m}$ ), a proportion, $h_{m 2}$ , regress to the INM1 class (at a rate $h_{m 2} y_{m}$ ), and the remaining proportion, $1 - (h_{m 1} + h_{m 2})$ , progress to the INM3 class (at a rate $(1 - h_{m 1} - h_{m 2}) y_{m}$ ). Males with INM3 (i.e., those in the $J_{3 m}$ class) leave this class at a rate $z_{m}$ . A proportion, $q_{m 1}$ , of these individuals recover (at a rate $q_{m 1} z_{m}$ ), a proportion, $q_{m 2}$ , regress to the INM2 class (at a rate $q_{m 2} z_{m}$ ), and the remaining proportion, $1 - (q_{m 1} + q_{m 2})$ , progress to the HPV-related cancer class (at a rate $(1 - q_{m 1} - q_{m 2}) z_{m}$ ). Males with HPV-related cancer (i.e., those in the $C_{m}$ class) recover from cancer at a rate $r_{m}$ .

The main assumptions made in the formulation of the model (2.5) include (inter alia):

(i) Homogeneous mixing (i.e., every sexually active heterosexual individual is equally likely to mix with, and have sexual contact, with every other sexually active heterosexual individual)
(ii) Large population size (so that random or stochastic fluctuations can be ignored)
(iii) Exponentially-distributed waiting times in each epidemiological compartment
(iv) No prior immunity against HPV infection (i.e., no maternal immunity) and recovery from previous HPV infection induces permanent immunity against the acquisition of future infections (Elbasha and Galvani 2005; Malik et al. 2013)
(v) Immunity derived from HPV vaccine does not wane during the period of consideration for the current study (i.e., during the simulation period considered in this study) (Elbasha 2007)
(vi) Individuals in the detected CIN stages ( $Q_{id} ; i = 1, 2, 3$ ) are females who got detected for their pre-cancerous lesions during the regular Pap screening performed in the population (South Korea) every 2 years (Shin et al. 2022)
(vii) Males do not die of HPV-related cancers (although males die of anal, penile, or oropharyngeal cancers, the rates of occurrence of these deaths are quite negligible compared to deaths due to HPV-related cancers in females (Bray et al. 2024; de Martel et al. 2017)).

The model (2.5) is an extension of many HPV models in the literature, such as those in (Alsaleh and Gumel 2014; Elbasha 2007; Elbasha and Galvani 2005; Javame and Gumel 2016; Malik et al. 2013; Milwid et al. 2018) by (inter alia):

including both cohort and catch-up vaccination for sexually active unvaccinated susceptible females and males (only cohort vaccination is considered in (Alsaleh and Gumel 2014; Elbasha 2007; Elbasha and Galvani 2005; Javame and Gumel 2016; Malik et al. 2013))
incorporating multiple grades of pre-cancerous lesions (CIN for females, INM for heterosexual males) in the HPV dynamics (pre-cancerous stages are not included in (Elbasha 2007; Elbasha and Galvani 2005), and only CIN stages for females are included in (Malik et al. 2013))
allowing for regression from grade i to grade $i - 1$ ( $i = 2$ and 3) of the pre-cancerous CIN and INM stages for females and males, respectively (regression of detected pre-cancerous lesions are not included in (Javame and Gumel 2016), regression from grade 3 to grade 2 is not included in (Malik et al. 2013), and regression of undetected pre-cancerous lesions from grade 3 to grade 2 are not included in (Milwid et al. 2018))
allowing for HPV transmission by individuals in the pre-cancerous stages (this is not included in (Alsaleh and Gumel 2014; Malik et al. 2013), and HPV transmission by individuals with CIN3 is not included in (Milwid et al. 2018))
including a class of exposed (pre-symptomatic infectious) individuals (females and males), and allowing for HPV transmission by these individuals (no exposed class is included in (Elbasha 2007; Elbasha and Galvani 2005; Malik et al. 2013; Milwid et al. 2018))

A flow diagram of the model (2.5) is depicted in Figure 1 (the state variables and parameters of the model are described in Tables 1-3).

**Fig. 1** Schematic diagram of the model (2.5).

Basic qualitative properties

The basic qualitative properties of the model (2.5) will now be explored to assess its well-posedness. Specifically, the non-negativity, boundedness, and invariance of its solutions will be explored. We claim the following result (the proof of Lemma 2.1 is given in Appendix A).

Lemma 2.1

Consider the model (2.5) with non-negative initial conditions and strict positivity condition for the initial values of $Q_{2 u}, Q_{2 d}$ , and $J_{2 m}$ (i.e., $Q_{2 u} (0) > 0, Q_{2 d} (0) > 0$ , and $J_{2 m} (0) > 0$ ). Then, all solutions of the model remain non-negative for all time $t > 0$ .

Consider, next, the following biologically feasible region for the model: (2.5):

\begin{matrix} D_{S} = D_{S, f} \times D_{S, m} \subset R_{+}^{15} \times R_{+}^{11}, \end{matrix}

where,

D_{S, f} = {(S_{f}, V_{f}, E_{f}, I_{f}, P_{f}, Q_{1 u}, Q_{2 u}, Q_{3 u}, Q_{1 d}, Q_{2 d}, Q_{3 d}, C_{u}, C_{d}, R_{f}^{C}, R_{f}) \in R_{+}^{15} : 0 \leq N_{f} \leq \frac{π_{f}}{μ}}

and

D_{S, m} = {(S_{m}, V_{m}, E_{m}, I_{m}, P_{m}, J_{1 m}, J_{2 m}, J_{3 m}, C_{m}, R_{m}^{C}, R_{m}) \in R_{+}^{11} : 0 \leq N_{m} \leq \frac{π_{m}}{μ}} .

We claim the following result.

Lemma 2.2

The region $D_{S}$ is positively invariant with respect to the flow generated by the model (2.5), and attracts all solutions of the model. Furthermore, all solutions of the model are bounded for all time $t \geq 0$ .

Proof

Adding the equations for the rates of change of each of the epidemiological compartment of the model (2.5) shows that the equations for the rate of change of the total female and male populations are given, respectively, by

\begin{matrix} {\dot{N}}_{f} (t) = π_{f} - μ N_{f} - δ C_{u} - δ_{d} C_{d} \leq π_{f} - μ N_{f} and {\dot{N}}_{m} (t) = π_{m} - μ N_{m}, \end{matrix}

2.6

from which it, respectively, follows, using a standard comparison theorem (Lakshmikantham et al. 2016), that

\begin{matrix} N_{f} (t) \leq N_{f} (0) e^{- μ t} + \frac{π_{f}}{μ} [1 - e^{- μ t}] and N_{m} (t) = N_{m} (0) e^{- μ t} + \frac{π_{m}}{μ} [1 - e^{- μ t}] . \end{matrix}

2.7

Consequently,

N_{f} (t) \leq \frac{π_{f}}{μ}

and

N_{m} (t) \leq \frac{π_{m}}{μ}

for all

t \geq 0

provided that

N_{f} (0) \leq \frac{π_{f}}{μ}

and

N_{m} (0) \leq \frac{π_{m}}{μ}

. In addition, Lemma 2.1 shows that solutions with non-negative initial conditions stays non-negative for all

t \geq 0

. Thus, these results (for the upper bounds of

N_{f} (t)

and

N_{m} (t)

and the non-negativity of the solutions of the model) show that any solution of the model with initial condition inside

D_{S}

stays inside

D_{S}

for

t \geq 0

, and the region

D_{S}

is positively invariant with respect to the model (2.5). Moreover, it follows from (2.6) that

{\dot{N}}_{f} (t) < 0

whenever

N_{f} (t) > \frac{π_{f}}{μ}

, and

{\dot{N}}_{m} (t) < 0

whenever

N_{m} (t) > \frac{π_{a}}{μ}

. Thus, any solution of the model (2.5) either enters

D_{S}

in a finite time, or

N_{f} (t)

and

N_{m} (t)

approach

\frac{π_{f}}{μ}

and

\frac{π_{m}}{μ}

, respectively. Hence, all solutions of the model (2.5) in

R_{+}^{15} \times R_{+}^{11}

are bounded for all time

t \geq 0

, and

D_{S}

is attracting (i.e., attracts all solutions of the model (2.5) in

R_{+}^{15} \times R_{+}^{11}

).

□

Since the region $D_{S}$ is positively invariant and attracting, it can be concluded that the model (2.5) is well-posed mathematically and epidemiologically in the region $D_{S}$ . Hence, it is sufficient to study the dynamics of the model in the invariant and attracting region $D_{S}$ (Hethcote 2000).

Existence and local asymptotic stability of equilibria of the model (2.5)

Disease-free equilibrium (DFE)

The model (2.5) has a unique disease-free equilibrium (DFE) given by:

\begin{matrix} E_{0 s} = (S_{f}^{*}, V_{f}^{*}, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, S_{m}^{*}, V_{m}^{*}, 0, 0, 0, 0, 0, 0, 0, 0, 0), \end{matrix}

where,

S_{f}^{*} = \frac{(1 - q_{f}) π_{f}}{ξ_{f} + μ}, V_{f}^{*} = \frac{π_{f} (ξ_{f} + μ q_{f})}{μ (ξ_{f} + μ)}, S_{m}^{*} = \frac{(1 - q_{m}) π_{m}}{ξ_{m} + μ}, V_{m}^{*} = \frac{π_{m} (ξ_{m} + μ q_{m})}{μ (ξ_{m} + μ)}

(note that

q_{f}

and

q_{m}

are fractions, hence

S_{f}^{*} > 0

and

S_{m}^{*} > 0

). The stability of the DFE is explored using the next generation operator method (Diekmann et al. 1990; Van den Driessche and Watmough 2002). In particular, using the notation in (Van den Driessche and Watmough 2002), the associated matrices

F_{s}

(for the new infection terms near the DFE) and

V_{s}

(for the linear transition terms in the infected compartments) are given, respectively, by (note that 0

_{n \times n}

is the zero matrix of order n):

F_{s} = [\begin{matrix} 0_{11 \times 11} & e_{2}^{T} f_{2} \\ e_{1}^{T} f_{1} & 0_{7 \times 7} \end{matrix}] and V_{s} = [\begin{matrix} V_{1} & 0_{7 \times 7} \\ 0_{11 \times 11} & V_{2} \end{matrix}],

where,

\begin{matrix} e_{1} & = (1, 0, 0, 0, 0, 0, 0), \\ f_{1} & = \frac{β_{fm} c_{mf} p_{m}^{*}}{N_{m}^{*}} (η_{f}, 1, θ_{f}, θ_{f}, θ_{f} κ_{f}, θ_{f} κ_{f}, θ_{f} ν, θ_{f} ν, θ_{f} ν κ_{fd}, 0, 0), \\ e_{2} & = (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), f_{2} = \frac{β_{mf} c_{mf} p_{f}^{*}}{N_{m}^{*}} (η_{m}, 1, θ_{m}, θ_{m}, θ_{m} κ_{m}, θ_{m} κ_{m}, 0), \\ V_{1} & = [\begin{matrix} k_{1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ - σ_{f} & k_{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & - l_{1} & k_{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - l_{2} & k_{4} & - h_{u 2} y_{f} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - l_{3} & k_{5} & q_{u 2} z_{f} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - l_{4} & k_{6} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - d_{u 3} g_{f} & 0 & 0 & k_{7} & - h_{d 2} ζ_{2} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - h_{u 3} y_{f} & 0 & - l_{6} & k_{8} & - q_{d 2} ζ_{3} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - q_{u 3} z_{f} & 0 & - l_{7} & k_{9} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - l_{5} & 0 & 0 & 0 & k_{10} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - l_{8} & - w_{f} & k_{11} \end{matrix}], \end{matrix}

and,

\begin{matrix} V_{2} & = [\begin{matrix} e_{1} & 0 & 0 & 0 & 0 & 0 & 0 \\ - σ_{m} & e_{2} & 0 & 0 & 0 & 0 & 0 \\ 0 & - (1 - b_{m}) ψ_{m} & e_{3} & 0 & 0 & 0 & 0 \\ 0 & 0 & - (1 - k_{m}) α_{m} & e_{4} & 0 & 0 & 0 \\ 0 & 0 & 0 & - (1 - d_{m}) g_{m} & e_{5} & 0 & 0 \\ 0 & 0 & 0 & 0 & - f_{1} & e_{6} & 0 \\ 0 & 0 & 0 & 0 & 0 & - f_{2} & e_{7} \end{matrix}], \end{matrix}

with,

p_{f}^{*} = S_{f}^{*} + (1 - ε_{v}) V_{f}^{*} = \frac{π_{f}}{μ} (1 - ε_{v} \frac{ξ_{f} + μ q_{f}}{ξ_{f} + μ}), p_{m}^{*} = S_{m}^{*} + (1 - ε_{v}) V_{m}^{*} = \frac{π_{m}}{μ} (1 - ε_{v} \frac{ξ_{m} + μ q_{m}}{ξ_{m} + μ}), k_{1} = σ_{f} + μ, k_{2} = ψ_{f} + μ, k_{3} = α_{f} + μ, k_{4} = g_{f} + μ, k_{5} = y_{f} + μ, k_{6} = z_{f} + μ, k_{7} = ζ_{1} + μ, k_{8} = ζ_{2} + μ, k_{9} = ζ_{3} + μ, k_{10} = w_{f} + r_{u} + δ + μ, k_{11} = r_{d} + δ_{d} + μ, l_{1} = (1 - b_{f}) ψ_{f}, l_{2} = (1 - k_{f}) α_{f}, l_{3} = (1 - d_{u 1} - d_{u 3}) g_{f}, l_{4} = (1 - h_{u 1} - h_{u 2} - h_{u 3}) y_{f}, l_{5} = (1 - q_{u 1} - q_{u 2} - q_{u 3}) z_{f}, l_{6} = (1 - d_{d 1}) ζ_{1}, l_{7} = (1 - h_{d 1} - h_{d 2}) ζ_{2}, l_{8} = (1 - q_{d 1} - q_{d 2}) ζ_{3}, e_{1} = σ_{m} + μ, e_{2} = ψ_{m} + μ, e_{3} = α_{m} + μ, e_{4} = g_{m} + μ, e_{5} = y_{m} + μ, e_{6} = z_{m} + μ, e_{7} = r_{m} + μ, f_{1} = (1 - h_{m 1} - h_{m 2}) y_{m}

, and

f_{2} = (1 - q_{m 1} - q_{m 2}) z_{m}

. It is convenient to define the following quantity:

\begin{matrix} R_{vs} = ρ (F_{s} V_{s}^{- 1}) = \sqrt{R_{f} R_{m}}, \end{matrix}

2.8

where,

\begin{matrix} R_{f} = & \frac{β_{fm} c_{mf} μ p_{m}^{*}}{π_{m}} [\frac{η_{f}}{k_{1}} + \frac{σ_{f}}{k_{1} k_{2}} + \frac{θ_{f} σ_{f} l_{1}}{k_{1} k_{2} k_{3}} +) \\ + (\frac{θ_{f} σ_{f} l_{1} l_{2} (ι_{f 2} + κ_{f} l_{3} k_{6} + κ_{f} l_{3} l_{4})}{k_{1} k_{2} k_{3} ι_{f 1}} + \frac{θ_{f} ν σ_{f} l_{1} l_{2} (j_{f 1} + κ_{fd} j_{f 2} + κ_{fd} j_{f 3})}{k_{1} k_{2} k_{3} ι_{f 1} ι_{f 3}}], \\ R_{m} = & \frac{β_{mf} c_{mf} μ p_{f}^{*}}{π_{m}} [\frac{η_{m}}{e_{1}} + \frac{σ_{m}}{e_{1} e_{2}} + \frac{θ_{m} σ_{m} ψ_{m} (1 - b_{m})}{e_{1} e_{2} e_{3}}) \\ (+ \frac{θ_{m} σ_{m} ψ_{m} (1 - b_{m}) α_{m} (1 - k_{m}) {ι_{m 2} + κ_{m} g_{m} (1 - d_{m}) e_{6} + κ_{m} g_{m} (1 - d_{m}) f_{1}}}{e_{1} e_{2} e_{3} ι_{m 1}}], \end{matrix}

with,

\begin{matrix} ι_{f 1} & = k_{4} k_{5} k_{6} - l_{3} h_{u 2} y_{f} k_{6} - k_{4} l_{4} q_{u 2} z_{f}, ι_{f 2} = k_{5} k_{6} - l_{4} q_{u 2} z_{f}, \\ ι_{f 3} & = k_{7} k_{8} k_{9} - l_{6} h_{d 2} ζ_{2} k_{9} - k_{7} l_{7} q_{d 2} ζ_{3}, \\ j_{f 1} & = d_{u 3} g_{f} (k_{5} k_{6} - l_{4} q_{u 2} z_{f}) (k_{8} k_{9} - l_{7} q_{d 2} ζ_{3}) + l_{3} h_{u 3} y_{f} h_{d 2} ζ_{2} k_{6} k_{9} \\ + l_{3} l_{4} q_{u 3} z_{f} q_{d 2} ζ_{3} h_{d 2} ζ_{2}, \\ j_{f 2} & = d_{u 3} g_{f} l_{6} k_{9} (k_{5} k_{6} - l_{4} q_{u 2} z_{f}) + l_{3} h_{u 3} y_{f} k_{6} k_{7} k_{9} + l_{3} l_{4} q_{u 3} z_{f} q_{d 2} ζ_{3} k_{7}, \\ j_{f 3} & = d_{u 3} g_{f} l_{6} l_{7} (k_{5} k_{6} - l_{4} q_{u 2} z_{f}) + l_{3} h_{u 3} y_{f} l_{7} k_{6} k_{7} + l_{3} l_{4} q_{u 3} z_{f} (k_{7} k_{8} - l_{6} h_{d 2} ζ_{2}), \\ ι_{m 1} & = e_{4} e_{5} e_{6} - (1 - d_{m}) g_{m} h_{m 2} y_{m} e_{6} - e_{4} f_{1} q_{m 2} z_{m}, ι_{m 2} = e_{5} e_{6} - f_{1} q_{m 2} z_{m} . \end{matrix}

It can be shown that the expression for

R_{vs}

, given in (2.8), is non-negative (see Appendix B for the proof). The result below follows from Theorem 2 of (Van den Driessche and Watmough 2002).

Theorem 2.3

The DFE, $E_{0 s}$ , of the model (2.5) is locally-asymptotically stable if $R_{vs} < 1$ , and unstable if $R_{vs} > 1 .$

The epidemiological implication of Theorem 2.3 is that a small influx of HPV-infected individuals into the community (i.e., initial number of infected individuals lie within the basin of attraction of the DFE) will not generate a large outbreak of HPV in the community if the interventions implemented in the community (vaccination and Pap screening in this case) can bring, and maintain,

R_{vs}

to a value less than unity (and the disease ultimately dies out). The threshold quantity

R_{vs}

is the control reproduction number of the model (2.5). It measures the average number of new HPV cases generated by a typical HPV-infected individual (female or male) if introduced into a community where a certain proportion is vaccinated and/or receiving Pap screening. It is convenient to let

\begin{matrix} v_{f}^{*} : = \frac{V_{f}^{*}}{N_{f}^{*}} = \frac{ξ_{f} + μ q_{f}}{ξ_{f} + μ}, v_{m}^{*} : = \frac{V_{m}^{*}}{N_{m}^{*}} = \frac{ξ_{m} + μ q_{m}}{ξ_{m} + μ} \end{matrix}

2.9

be the fraction of females and males vaccinated at steady-state, respectively. It follows, by using (2.9) in (2.8), that the control reproduction number

R_{vs}

can be re-written in terms of these steady-state vaccination fractions as:

\begin{matrix} R_{vs} (v_{f}^{*}, v_{m}^{*}) = R_{cs} \sqrt{(1 - ε_{v} v_{m}^{*}) (1 - ε_{v} v_{f}^{*})}, \end{matrix}

2.10

where,

\begin{matrix} R_{cs} = R_{vs} |_{q_{f} = ξ_{f} = q_{m} = ξ_{m} = 0} \end{matrix}

2.11

is the control reproduction number of the model in the absence of vaccination (i.e., the control reproduction number without vaccination but with screening). It is convenient to define

R_{0} = R_{vs} |_{q_{f} = ξ_{f} = q_{m} = ξ_{m} = d_{u 3} = h_{u 3} = q_{u 3} = ν = 0} = \sqrt{R_{0 f} R_{0 m}},

where,

\begin{matrix} R_{0 f} = & β_{fm} c_{mf} [\frac{η_{f}}{k_{1}} + \frac{σ_{f}}{k_{1} k_{2}} + \frac{θ_{f} σ_{f} l_{1}}{k_{1} k_{2} k_{3}} + \frac{θ_{f} σ_{f} l_{1} l_{2} ({\tilde{ι}}_{f 2} + κ_{f} \tilde{l_{3}} k_{6} + κ_{f} \tilde{l_{3}} \tilde{l_{4}})}{k_{1} k_{2} k_{3} {\tilde{ι}}_{f 1}}], \\ R_{0 m} = & \frac{β_{mf} c_{mf} π_{f}}{π_{m}} [\frac{η_{m}}{e_{1}} + \frac{σ_{m}}{e_{1} e_{2}} + \frac{θ_{m} σ_{m} ψ_{m} (1 - b_{m})}{e_{1} e_{2} e_{3}}) \\ (+ \frac{θ_{m} σ_{m} ψ_{m} (1 - b_{m}) α_{m} (1 - k_{m}) {ι_{a 2} + κ_{m} g_{m} (1 - d_{m}) e_{6} + κ_{m} g_{m} (1 - d_{m}) f_{1}}}{e_{1} e_{2} e_{3} ι_{a 1}}], \end{matrix}

with,

{\tilde{ι}}_{f 1} = ι_{f 1} |_{d_{u 3} = h_{u 3} = 0} = k_{4} k_{5} k_{6} - \tilde{l_{3}} h_{u 2} y_{f} k_{6} - k_{4} \tilde{l_{4}} q_{u 2} z_{f}, {\tilde{ι}}_{f 2} = ι_{f 2} |_{h_{u 3} = 0} = k_{5} k_{6} - \tilde{l_{4}} q_{u 2} z_{f}, \tilde{l_{3}} = l_{3} {|_{d_{u 3} = 0} = (1 - d_{u 1}) g_{f}, \tilde{l_{4}} = l_{4} |}_{h_{u 3} = 0} = (1 - h_{u 1} - h_{u 2}) y_{f}

. The quantity

R_{0}

is the basic reproduction number of the model (which represents the average number of new HPV cases generated by a typical HPV-infected individual introduced into a completely susceptible population). Then, the quantity

R_{cs}

, given in (2.11), can be expressed as

R_{cs} = \sqrt{R_{csf} R_{0 m}},

where

\begin{matrix} R_{csf} = & R_{0 f} - (1 - \frac{{\tilde{ι}}_{f 1}}{ι_{f 1}}) β_{fm} c_{mf} [\frac{θ_{f} σ_{f} l_{1} l_{2} ({\tilde{ι}}_{f 2} + κ_{f} \tilde{l_{3}} k_{6} + κ_{f} \tilde{l_{3}} \tilde{l_{4}})}{k_{1} k_{2} k_{3} {\tilde{ι}}_{f 1}}] \\ - \frac{θ_{f} σ_{f} l_{1} l_{2} κ_{f} (\tilde{l_{3}} h_{u 3} y_{f} + d_{u 3} g_{f} \tilde{l_{4}} - d_{u 3} g_{f} h_{u 3} y_{f})}{k_{1} k_{2} k_{3} ι_{f 1}} \\ + β_{fm} c_{mf} [\frac{θ_{f} ν σ_{f} l_{1} l_{2} (j_{f 1} + κ_{fd} j_{f 2} + κ_{fd} j_{f 3})}{k_{1} k_{2} k_{3} ι_{f 1} ι_{f 3}}] . \end{matrix}

Furthermore, it can be seen, using (2.10), that the gradient of

R_{vs}

is:

\begin{matrix} (\frac{\partial R_{vs}}{\partial v_{f}^{*}}, \frac{\partial R_{vs}}{\partial v_{m}^{*}}) = \frac{- ε_{v} R_{cs}}{2 \sqrt{(1 - ε_{v} v_{m}^{*}) (1 - ε_{v} v_{f}^{*})}} (1 - ε_{v} v_{m}^{*}, 1 - ε_{v} v_{f}^{*}), \end{matrix}

from which it follows that both

\frac{\partial R_{vs}}{\partial v_{f}^{*}}

and

\frac{\partial R_{vs}}{\partial v_{m}^{*}}

are strictly negative for

0 < ε_{v} < 1

and

0 \leq v_{f}^{*}, v_{m}^{*} \leq 1

. Thus,

R_{vs}

is a decreasing function of

v_{f}^{*}

and

v_{m}^{*}

for any fixed values of the other parameters (this is depicted in a 3-dimensional plot in Figure 2 (a), for

R_{vs}

as a function of the two vaccination coverages). Since a reduction in the control reproduction number implies a reduction in the disease burden (as measured in terms of the number of HPV infections, HPV-related cancer cases and mortality etc.), the above analysis shows that the HPV vaccination program implemented in Korea will have a positive impact in reducing HPV burden in the country. In other words, a vaccination program (using an HPV vaccine with protective efficacy

0 < ε_{v} < 1

) combined with Pap screening will result in a decrease in HPV burden (and related cancers) regardless of the values of the proportions of females and males vaccinated at steady-state. This result is also illustrated in a contour plot of the control reproduction number,

R_{vs}

, as a function of

v_{f}^{*}

and

v_{m}^{*}

, depicted in Figure 2 (b). This figure shows that

R_{vs}

decreases (i.e., disease burden decreases), with increasing values of the proportion of females and males vaccinated at steady-state. Furthermore, this figure shows that, using the model (2.5) with the baseline values in Table 9, the use of an HPV vaccine (with protective efficacy of 95%) and baseline Pap screening coverage can lead to the effective control of HPV in Korea for a certain combination of the vaccination coverage for females and males (e.g., 88% of females vaccinated combined with 65% of males vaccinated at steady-state; since this combination of efficacy and coverages will bring

R_{vs}

to a value less than one, which is needed for the local asymptotic stability of the disease-free equilibrium, in line with Theorem 2.3). It can also be seen from this contour plot that

R_{vs}

can be brought to a value less than one if nearly 100% of females or males are vaccinated (but this is generally not realistically attainable).

Table 9 Baseline values of the parameters for the full model (2.5)

Parameter	Baseline Value (per year)	Source
$π_{f}$ , $π_{m}$	358000, 372000	(Statistics Korea 2000)
$μ$	1/65	(Malik et al. 2013)
$q_{f}$ , $ξ_{f}$	0.05 (dimensionless), 0.01	(Choi et al. 2022)
$q_{m}$ , $ξ_{m}$	0 (dimensionless), 0	Assumed
$c_{mf}$	2	(Podder and Gumel 2009)
$c_{fm}$	$\frac{c_{mf} N_{f} (t)}{N_{m} (t)} = \frac{2 N_{f} (t)}{N_{m} (t)}$	From value of $c_{mf}$ above
$β_{mf}$	0.72 (dimensionless)	Derived from (Podder and Gumel 2009)
$β_{fm}$	$\frac{β_{fm} c_{fm}}{c_{fm} (t)} = \frac{2.5}{c_{fm} (t)}$ (dimensionless)	Derived from (Podder and Gumel 2009)
$η_{f}, η_{m}$	0.5 (dimensionless)	(Javame and Gumel 2016)
$θ_{f}, θ_{m}$	0.9 (dimensionless)	(Javame and Gumel 2016)
$κ_{f}, κ_{m}$	1.5 (dimensionless)	(Javame and Gumel 2016)
$ν$	0.05 (dimensionless)	Assumed
$κ_{fd}$	1.5 (dimensionless)	(Javame and Gumel 2016)
$ε_{v}$	0.95 (dimensionless)	(Malik et al. 2013)
$σ_{f}, σ_{m}$	3	(Javame and Gumel 2016)
$(1 - b_{f}) ψ_{f}$	0.017	(Javame and Gumel 2016)
$(1 - b_{m}) ψ_{m} = r_{1} (1 - b_{f}) ψ_{f}$	0.0076	(Javame and Gumel 2016) and fitted value of $r_{1}$
$(1 - k_{f}) α_{f}$	0.07	(Javame and Gumel 2016)
$(1 - k_{m}) α_{m} = r_{1} (1 - k_{f}) α_{f}$	0.03	(Javame and Gumel 2016) and fitted value of $r_{1}$
$(1 - d_{u 1} - d_{u 3}) g_{f}$	0.33	Derived from model (3.1)
$(1 - d_{d 1}) ζ_{1}$	0.14	(Elbasha et al. 2007)
$(1 - d_{m}) g_{m} = r_{1} (1 - d_{u 1} - d_{u 3}) g_{f}$	0.14	Derived from model (3.1) and fitted value of $r_{1}$
$(1 - h_{u 1} - h_{u 2} - h_{u 3}) y_{f}$	0.33	Derived from model (3.1)
$(1 - h_{d 1} - h_{d 2}) ζ_{2}$	0.14	(Elbasha et al. 2007)
$(1 - h_{m 1} - h_{m 2}) y_{m} = r_{1} (1 - h_{u 1} - h_{u 2} - h_{u 3}) y_{f}$	0.14	Derived from model (3.1) and fitted value of $r_{1}$
$(1 - q_{u 1} - q_{u 2} - q_{u 3}) z_{f}$	0.98	Derived from model (3.1)
$(1 - q_{d 1} - q_{d 2}) ζ_{3}$	0.42	(Elbasha et al. 2007)
$(1 - q_{m 1} - q_{m 2}) z_{m} = r_{1} (1 - q_{u 1} - q_{u 2} - q_{u 3}) z_{f}$	0.43	Derived from model (3.1) and fitted value of $r_{1}$
$h_{u 2} y_{f}, h_{d 2} ζ_{2}$	0.133	(Elbasha et al. 2007)
$h_{m 2} y_{m}$	0.133	(Elbasha et al. 2007)
$q_{u 2} z_{f}, q_{d 2} ζ_{3}$	0.03	(Elbasha et al. 2007)
$q_{m 2} z_{m}$	0.03	(Elbasha et al. 2007)
$d_{u 3} g_{f}$ , $h_{u 3} y_{f}$ , $q_{u 3} z_{f}$	0.18, 0.24, 0.24	Derived from (Shin et al. 2022; Malik et al. 2013)
$w_{f}$	0.24	Derived from (Shin et al. 2022; Malik et al. 2013)
$b_{f} ψ_{f}$	0.495	(Malik et al. 2013)
$b_{m} ψ_{m} = r_{2} b_{f} ψ_{f}$	0.625	(Malik et al. 2013) and fitted value of $r_{2}$
$k_{f} α_{f}$	0.05	(Javame and Gumel 2016)
$k_{m} α_{m} = r_{2} k_{f} α_{f}$	0.063	(Javame and Gumel 2016) and fitted value of $r_{2}$
$d_{u 1} g_{f}, h_{u 1} y_{f}, q_{u 1} z_{f}$	0.01	(Javame and Gumel 2016)
$d_{d 1} ζ_{1}, h_{d 1} ζ_{2}, q_{d 1} ζ_{3}$	0.13	(Malik et al. 2013)
$d_{m} g_{m} = r_{2} d_{u 1} g_{f}, h_{m 1} y_{m} = r_{2} h_{u 1} y_{f}, q_{m 1} z_{m} = r_{2} q_{u 1} z_{f}$	0.013	(Javame and Gumel 2016) and fitted value of $r_{2}$
$r_{u}, r_{d}$	0.5, 0.63	(Javame and Gumel 2016)
$r_{m} = r_{2} r_{u}$	0.63	(Javame and Gumel 2016) and fitted value of $r_{2}$
$δ, δ_{d}$	0.23, 0.023	(Ministry of Health and Welfare 2021)

Using the baseline values of the parameters of the model (2.5), given in Table 9, it can be seen that the the control ( $R_{vs}$ ) and basic ( $R_{0}$ ) reproduction numbers for South Korea before the implementation of the National Immunization Program for HPV in 2016 are $R_{vs} = 3.0764$ and $R_{0} = 3.9858$ (this value of $R_{0}$ is within the range of $R_{0}$ for HPV reported in (Baussano et al. 2018)), respectively. In other words, based on the HPV vaccination coverage for females in Korea before the implementation of the National Immunization Program (estimated to be about 42% (Choi et al. 2022)) and the baseline Pap screening coverage (estimated to be about 61%), HPV (and associated cancers) will continue to persist in Korea (since $R_{0}$ is about 4, and the pre-NIP vaccination and Pap screening only reduce $R_{0}$ by about 23%). However, for post-NIP vaccination (where the vaccination coverage for females in Korea is estimated to be 88% (Choi et al. 2022; Kim et al. 2018)), the baseline value of $R_{vs}$ is 1.6181, which represents a 59% reduction of the basic reproduction number, $R_{0}$ . Here, too, HPV will persist (since $R_{vs} > 1$ ), albeit HPV burden will be significantly reduced, in comparison to the pre-NIP era.

Computation of vaccine-derived herd immunity threshold: A precise expression for the vaccine-derived herd immunity threshold (i.e., the minimum proportion of susceptible individuals in the community that needs to be vaccinated in order to protect the ones that cannot be vaccinated) can be derived by setting the control reproduction number (

R_{vs}

) to one and solving for the proportions of females and males vaccinated at steady-state (

v_{f}^{*}

and

v_{m}^{*}

) (Elbasha and Gumel 2021; Gumel et al. 2006; Malik et al. 2013; McLean and Blower 1993; Pant and Gumel 2024). Let

R_{cs} > 1

. Then, setting the control reproduction number

R_{vs}

, given by Equation (2.10), to one (and simplifying) shows that (noting that

0 < ε_{v} < 1

,

0 \leq v_{f}^{*}, v_{m}^{*} \leq 1

, and

R_{cs} > 1

):

\begin{matrix} (1 - ε_{v} v_{m}^{*}) (1 - ε_{v} v_{f}^{*}) = \frac{1}{{(R_{cs})}^{2}}, 0 \leq v_{f}^{*}, v_{m}^{*} \leq 1, and R_{cs} > 1 . \end{matrix}

2.12

It follows from the level curve

R_{vs} (v_{f}^{*}, v_{m}^{*}) = 1

, represented by the white curve in the contour plot depicted in Figure 2 (b), that increasing the steady-state vaccination coverage of one gender group (females or males) will reduce the coverage level of the other group needed to bring (and maintain)

R_{vs}

to a value below one. Some representative

(v_{f}^{*}, v_{m}^{*})

pairs are illustrated in Figure 3 (a), showing that if no susceptible males are vaccinated at steady-state, then at least 99% of the sexually active females need to be vaccinated to bring the reproduction threshold

R_{vs}

below one (so that the disease can be effectively controlled and/or eliminated). However, this figure further shows that, if half of the males in the community are vaccinated, then 93% of the females need to be vaccinated in order to bring the reproduction threshold below one. Furthermore, Figure 3 (a) shows that if the vaccination coverage of males is increased (from 50%) to 70%, then the vaccination coverage needed for the female population decreases to 85%. Finally, if the post-NIP vaccination coverage of females in Korea is maintained at the

88 %

level (i.e.,

v_{f}^{*} = 0.88

), then it follows from Figure 3 (a) that vaccinating

65 %

of males would be sufficient to bring

R_{vs}

below one in Korea. In other words, this study shows that HPV can be effectively controlled and/or eliminated in Korea if the current HPV vaccination coverage in females is complemented with HPV vaccination of males with at least 65% coverage at steady-state. Solving for

v_{m}^{*}

from (2.12) gives:

\begin{matrix} v_{m}^{*} = \frac{1}{ε_{v}} [1 - \frac{1}{{(R_{cs})}^{2} (1 - ε_{v} v_{f}^{*})}] = : {(v_{m}^{*})}_{c}, with 0 < ε_{v} < 1, and 0 \leq v_{f}^{*} \leq 1 . \end{matrix}

2.13

It follows from (2.10) and (2.13) that the control reproduction number,

R_{vs}

, is less than one if and only if

v_{m}^{*} > {(v_{m}^{*})}_{c} = f (v_{f}^{*} ; ε_{v}, R_{0})

, and it equals one if

v_{m}^{*} = {(v_{m}^{*})}_{c}

. Hence, the result of Theorem 2.3 can be re-written in terms of the vaccine-derived herd immunity threshold (

{(v_{f}^{*})}_{c}

) as follows:

**Fig. 2** (a) Three-dimensional and (b) contour plot of the control reproduction number, $R_{vs}$ , of the model (2.5), as a function of vaccine coverage of females at steady-state ( $v_{f}^{*}$ ) and vaccine coverage of males at steady-state ( $v_{m}^{*}$ ), for vaccine efficacy $ε_{v} = 0.95$ . The white curve represents the 1-level set $R_{vs} (v_{f}^{*}, v_{m}^{*}) = 1$ . The other parameter values used to generate the plots are as given in Table 9.

Theorem 2.4

The DFE, $E_{0 s}$ , of the model (2.5) is locally-asymptotically stable if $v_{m}^{*} > {(v_{m}^{*})}_{c}$ (i.e., $R_{vs} < 1$ ), and unstable if $v_{m}^{*} < {(v_{m}^{*})}_{c}$ (i.e., $R_{vs} > 1$ ).

Theorem 2.4 implies that HPV will be eliminated in Korea if the proportion of males vaccinated at steady-state exceeds the herd immunity threshold value ${(v_{m}^{*})}_{c}$ . In other words, Theorem 2.4 shows that, using the post-NIP vaccination coverage for females (i.e., $v_{f}^{*} = 0.88$ ) and the other parameters of the model (2.5) maintained at their baseline values given in Table 9, vaccine-derived herd immunity can be achieved in Korea by vaccinating 65% of males (i.e., ${(v_{m}^{*})}_{c} = 0.65$ for Korea). However, if males are not vaccinated (i.e., $v_{m}^{*} = 0$ ), then, $v_{f}^{*} = 0.99$ (i.e., in this case, at least 99% of females must be vaccinated to achieve herd immunity, and, consequently, eliminate the disease).

Optimization of total vaccination coverage at steady-state: The objective here is to determine the optimal values of the steady-state vaccination coverage pair,

(v_{f}^{*}, v_{m}^{*})

, needed to minimize the total vaccination coverage (for both females and males in Korea) that results in

R_{vs} < 1

in the country. In other words, the aim here is to minimize the following objective functional, for the total vaccination coverage in Korea:

\begin{matrix} \frac{V_{f}^{*} + V_{m}^{*}}{N_{f}^{*} + N_{m}^{*}} = (\frac{N_{f}^{*}}{N_{f}^{*} + N_{m}^{*}}) \frac{V_{f}^{*}}{N_{f}^{*}} + (\frac{N_{m}^{*}}{N_{f}^{*} + N_{m}^{*}}) \frac{V_{m}^{*}}{N_{m}^{*}} = : n_{f}^{*} v_{f}^{*} + n_{m}^{*} v_{m}^{*}, \end{matrix}

subject to reducing the value of

R_{vs}

below one. More formally, the aim is to

\begin{matrix} minimize n_{f}^{*} v_{f}^{*} + n_{m}^{*} v_{m}^{*} subject to 0 \leq v_{f}^{*}, v_{m}^{*} \leq 1, and R_{vs} \leq 1, \end{matrix}

2.14

where

R_{vs}

is given by (2.10) with fixed vaccine efficacy

ε_{v} = 0.95

. The solution to this nonlinear optimization problem is obtained by using the geometric approach introduced by Elbasha and Gumel (Elbasha and Gumel 2021). The approach entails defining the feasible set (that satisfies the constraints of the optimization problem) and finding the points on the lines

n_{f}^{*} v_{f}^{*} + n_{m}^{*} v_{m}^{*} = k

(with

k \in R

, and with slope of the line,

d v_{m}^{*} / d v_{f}^{*}

, equal to

- n_{f}^{*} / n_{m}^{*}

) which intersect the feasible set. The optimal points will then lie on the line with the lowest

y -

intercept on the

(v_{f}^{*}, v_{m}^{*})

-plane. The feasible set for the optimization problem (2.14) is

{(v_{f}^{*}, v_{m}^{*}) \in {[0, 1]}^{2} | v_{m}^{*} \geq {(v_{m}^{*})}_{c} = f (v_{f}^{*})},

and is depicted by the shaded region in Figure 3 (b)). Since the line

n_{f}^{*} v_{f}^{*} + n_{m}^{*} v_{m}^{*} = k

has slope

- n_{f}^{*} / n_{m^{*}} = - π_{f} / π_{m} > - 1

, it follows (from the fact that the feasible set is symmetric with respect to the line

v_{m}^{*} = v_{f}^{*}

(i.e., the inequality

v_{m}^{*} \geq f (v_{f}^{*})

remains the same if the two variables are interchanged) and Figure 3 (b)) that k is minimized when the line intersects the feasible set uniquely on the x-axis. The intersection occurs at the unique point on the

v_{f}^{*}

-intercept of the line

v_{m}^{*} = f (v_{f}^{*})

, given by (where

\bar{v_{f}^{*}}

denotes the optimal value for

v_{f}^{*}

):

\begin{matrix} (v_{f}^{*}, v_{m}^{*}) = (\bar{v_{f}^{*}}, 0) = (\frac{1}{ε_{v}} (1 - \frac{1}{R_{cs}^{2}}), 0) . \end{matrix}

2.15

Since the value of

\bar{v_{f}^{*}}

, computed using the baseline values of the parameters in the expression for

R_{0}

given in Table 9 with

ε_{v} = 0.95

, is 0.99, it follows from (2.15) that the total vaccination coverage is optimized if 99% of females and 0% of males are vaccinated at steady-state. Thus, this study shows that the current (post-NIP) HPV vaccination coverage (of 88% vaccination of females) is insufficient to lead to vaccine-derived herd immunity. More (11%) unvaccinated females need to be vaccinated to achieve such immunity (with minimum number of vaccinated individuals). If 88% (i.e., the current vaccination coverage of females) is the maximum achievable vaccine uptake in Korea in each group (females and males) due to factors including vaccine hesitancy, limited availability and access to public health resources, and individual medical conditions, then the optimal steady-state vaccination coverage is 88% of females and 65% of males.

**Fig. 3** (a) The level set $R_{vs} (v_{f}^{*}, v_{m}^{*}) = 1$ (shown in solid blue curve) for the model (2.5), as a function of steady-state vaccination coverage for females ( $v_{f}^{*}$ ) and males ( $v_{m}^{*}$ ), showing some of the, among the infinitely many, $(v_{f}^{*}, v_{m}^{*})$ pairs that fall on the level set. (b) Geometric interpretation of the constrained optimization problem (2.14) associated with the model. The optimal solution in the feasible (shaded gray) set that minimizes $n_{f}^{*} v_{f}^{*} + n_{m}^{*} v_{m}^{*}$ is the point $(\bar{v_{f}^{*}}, 0)$ (i.e., the point in the shaded region that intersects the blue line, which is the line closest to the origin among the infinitely many lines with slope $- n_{f}^{*} / n_{m}^{*}$ , which also intersect the feasible set). The parameter values used to generate this figure are as given in Table 9.

Existence and local asymptotic stability of endemic equilibria: special case

In this section, conditions for the existence of endemic equilibria (EE) of the model (2.5) will be explored for a special case with negligible mortality for females due to cervical cancer (i.e., the model with

δ = δ_{d} = 0

will now be considered). Although this assumption (for the negligible disease-induced mortality) is made for mathematical tractability, it can also be epidemiologically justified for the dynamics of HPV and its related cancers in South Korea, where the average number of deaths caused by cervical cancer each year ranges from 1, 000 to 1, 500 (which represents about

0.006 %

to

0.009 %

of the sexually active females in South Korea). Setting

δ = δ_{d} = 0

in the model (2.5) shows that

N_{f} (t) \to π_{f} / μ

and

N_{m} (t) \to π_{m} / μ

as

t \to \infty

. For the purpose of the analysis in this subsection, the total populations

N_{f} (t)

and

N_{m} (t)

will now be replaced by their limiting values

N_{f}^{*} = π_{f} / μ

and

π_{m} / μ

, respectively. It should be mentioned that, since females with cervical cancer (i.e., individuals in the

C_{u}

and

C_{d}

compartments) do not transmit HPV, the control reproduction number

R_{vs}

of the special case of the model (2.5) is the same as that for the full model (2.5). Let

\begin{matrix} E_{1 s} = ( & S_{f}^{* *}, V_{f}^{* *}, E_{f}^{* *}, I_{f}^{* *}, P_{f}^{* *}, Q_{1 u}^{* *}, Q_{2 u}^{* *}, Q_{3 u}^{* *}, Q_{1 d}^{* *}, Q_{2 d}^{* *}, Q_{3 d}^{* *}, C_{u}^{* *}, C_{d}^{* *}, R_{f}^{C * *}, R_{f}^{* *}, \\ S_{m}^{* *}, V_{m}^{* *}, E_{m}^{* *}, I_{m}^{* *}, P_{m}^{* *}, J_{1 m}^{* *}, J_{2 m}^{* *}, J_{3 m}^{* *}, C_{m}^{* *}, R_{m}^{C * *}, R_{m}^{* *}) \end{matrix}

represent any arbitrary endemic equilibrium (EE) of the aforementioned special case of the heterosexual model (2.5). In addition, substituting

N_{f} (t) = π_{f} / μ

and

N_{m} (t) = π_{m} / μ

into the expressions for the forces of infections, given in (2.4), shows that:

\begin{matrix} λ_{mf}^{* *} = & (β_{mf}, c_{mf}, μ) [\frac{I_{m}^{* *} + η_{m} E_{m}^{* *} + θ_{m} [P_{m}^{* *} + J_{1 m}^{* *} + κ_{m} (J_{2 m}^{* *} + J_{3 m}^{* *})]}{π_{m}}], \\ λ_{fm}^{* *} = & (β_{fm}, c_{mf}, μ) [\frac{I_{f}^{* *} + η_{f} E_{f}^{* *} + θ_{f} [P_{f}^{* *} + Q_{1 u}^{* *} + κ_{f} (Q_{2 u}^{* *} + Q_{3 u}^{* *}) + ν {Q_{1 d}^{* *} + κ_{fd} (Q_{2 d}^{* *} + Q_{3 d}^{* *})}]}{π_{m}}] . \end{matrix}

2.16

It follows then that the special case of the model (2.5) at the endemic equilibrium satisfies the following steady-state relations:

\begin{matrix} S_{f}^{* *} & = \frac{(1 - q_{f}) π_{f}}{λ_{mf}^{* *} + ξ_{f} + μ}, V_{f}^{* *} = \frac{π_{f}}{(1 - ε_{v}) λ_{mf}^{* *} + μ} [q_{f} + \frac{(1 - q_{f}) ξ_{f}}{λ_{mf}^{* *} + ξ_{f} + μ}], \\ E_{f}^{* *} & = \frac{λ_{mf}^{* *}}{k_{1}} \cdot \frac{π_{f} (1 - ε_{v}) λ_{mf}^{* *} + π_{f} μ (1 - q_{f} ε_{v}) + π_{f} (1 - ε_{v}) ξ_{f}}{(λ_{mf}^{* *} + ξ_{f} + μ) [(1 - ε_{v}) λ_{mf}^{* *} + μ]}, \\ I_{f}^{* *} & = \frac{σ_{f} E_{f}^{* *}}{k_{2}}, P_{f}^{* *} = \frac{l_{1} I_{f}^{* *}}{k_{3}}, Q_{1 u}^{* *} = \frac{l_{2} P_{f}^{* *} + h_{u 2} y_{f} Q_{2 u}^{* *}}{k_{4}}, Q_{2 u}^{* *} = \frac{l_{3} Q_{1 u}^{* *} + q_{u 2} z_{f} Q_{3 u}^{* *}}{k_{5}}, \\ Q_{3 u}^{* *} & = \frac{l_{4} Q_{2 u}^{* *}}{k_{6}}, Q_{1 d}^{* *} = \frac{d_{u 3} g_{f} Q_{1 u}^{* *} + h_{d 2} ζ_{2} Q_{2 d}^{* *}}{k_{7}}, Q_{2 d}^{* *} = \frac{h_{u 3} y_{f} Q_{2 u}^{* *} + l_{6} Q_{1 d}^{* *} + q_{d 2} ζ_{3} Q_{3 d}^{* *}}{k_{8}}, \\ Q_{3 d}^{* *} & = \frac{q_{u 3} z_{f} Q_{3 u}^{* *} + l_{7} Q_{2 d}^{* *}}{k_{9}}, C_{u}^{* *} = \frac{l_{5} Q_{3 u}^{* *}}{k_{10}}, C_{d}^{* *} = \frac{l_{8} Q_{3 d}^{* *} + w_{f} C_{u}^{* *}}{k_{11}}, R_{f}^{C * *} = \frac{r_{u} C_{u}^{* *} + r_{d} C_{d}^{* *}}{μ}, \\ R_{f}^{* *} & = \frac{1}{μ} (b_{f} ψ_{f} I_{f}^{* *} + k_{f} α_{f} P_{f}^{* *} + d_{u 1} g_{f} Q_{1 u}^{* *} + h_{u 1} y_{f} Q_{2 u}^{* *} + q_{u 1} z_{f} Q_{3 u}^{* *} + d_{d 1} ζ_{1} Q_{1 d}^{* *} \\ + h_{d 1} ζ_{2} Q_{2 d}^{* *} + q_{d 1} ζ_{3} Q_{3 d}^{* *}), \\ S_{m}^{* *} & = \frac{(1 - q_{m}) π_{m}}{λ_{fm}^{* *} + ξ_{m} + μ}, V_{m}^{* *} = \frac{π_{m}}{(1 - ε_{v}) λ_{fm}^{* *} + μ} [q_{m} + \frac{(1 - q_{m}) ξ_{m}}{λ_{fm}^{* *} + ξ_{m} + μ}], \\ E_{m}^{* *} & = \frac{λ_{fm}^{* *}}{e_{1}} \cdot \frac{π_{m} (1 - ε_{v}) λ_{fm}^{* *} + π_{m} μ (1 - q_{m} ε_{v}) + π_{m} (1 - ε_{v}) ξ_{m}}{(λ_{fm}^{* *} + ξ_{m} + μ) [(1 - ε_{v}) λ_{fm}^{* *} + μ]}, \\ I_{m}^{* *} & = \frac{σ_{m} E_{m}^{* *}}{e_{2}}, P_{m}^{* *} = \frac{(1 - b_{m}) ψ_{m} I_{m}^{* *}}{e_{3}}, \\ J_{1 m}^{* *} & = \frac{(1 - k_{m}) α_{m} P_{m}^{* *} + h_{m 2} y_{m} J_{2 m}^{* *}}{e_{4}}, J_{2 m}^{* *} = \frac{(1 - d_{m}) g_{m} J_{1 m}^{* *} + q_{m 2} z_{m} J_{3 m}^{* *}}{e_{5}}, J_{3 m}^{* *} = \frac{f_{1} J_{2 m}^{* *}}{e_{6}}, \\ C_{m}^{* *} & = \frac{f_{2} J_{3 m}^{* *}}{e_{7}}, R_{m}^{C * *} = \frac{r_{m} C_{m}^{* *}}{μ}, \\ R_{m}^{* *} & = \frac{1}{μ} (b_{m} ψ_{m} I_{m}^{* *} + k_{m} α_{m} P_{m}^{* *} + d_{m} g_{m} J_{1 m}^{* *} + h_{m 1} y_{m} J_{2 m}^{* *} + q_{m 1} z_{m} J_{3 m}^{* *}) . \end{matrix}

2.17

Substituting (2.17) into (2.16), and simplifying, gives

\begin{matrix} λ_{mf}^{* *} = \frac{λ_{fm}^{* *} R_{m} π_{m} μ (ξ_{f} + μ) {(1 - ε_{v}) λ_{fm}^{* *} + μ (1 - q_{m} ε_{v}) + (1 - ε_{v}) ξ_{m}}}{π_{f} [μ (1 - q_{f} ε_{v}) + (1 - ε_{v}) ξ_{f}] (λ_{fm}^{* *} + ξ_{m} + μ) [(1 - ε_{v}) λ_{fm}^{* *} + μ]}, \end{matrix}

2.18

and,

\begin{matrix} λ_{fm}^{* *} = \frac{λ_{mf}^{* *} R_{f} π_{f} μ (ξ_{m} + μ) {(1 - ε_{v}) λ_{mf}^{* *} + μ (1 - q_{f} ε_{v}) + (1 - ε_{v}) ξ_{f}}}{π_{m} [μ (1 - q_{m} ε_{v}) + (1 - ε_{v}) ξ_{m}] (λ_{mf}^{* *} + ξ_{f} + μ) [(1 - ε_{v}) λ_{mf}^{* *} + μ]} . \end{matrix}

2.19

Substituting (2.19) into (2.18) leads to the following quartic in

λ_{mf}^{* *}

:

\begin{matrix} \sum_{i = 0}^{4} a_{i} {(λ_{mf}^{* *})}^{i} = 0, \end{matrix}

2.20

where the coefficients

a_{i}

(with

i = 0, 1, \dots, 4

) are given in Appendix C. It follows from the expressions for

a_{i}

(in Appendix C) that the coefficient

a_{4}

is always non-negative (given that the parameters of the model (2.5) are all non-negative), and the coefficients

a_{1}, a_{2}

, and

a_{3}

are all non-negative if

R_{vs} \leq 1

. Furthermore,

a_{0} > 0

if and only if

R_{vs} < 1,

and

a_{0} = 0

if

R_{vs} = 1 .

Thus, it follows, using the Descartes’ Rule of Signs, that the special case of the model (2.5) can have one or three positive (endemic) equilibria whenever the control reproduction number

R_{vs}

exceeds unity, and no endemic equilibria otherwise. This result is summarized below.

Theorem 2.5

The special case of the model (2.5), with $δ = δ_{d} = 0$ , has one or three positive endemic equilibria whenever $R_{vs} > 1,$ and no endemic equilibrium whenever $R_{vs} \leq 1 .$

The various possibilities for the number of equilibria for the special case of the model are tabulated in Table 4, from which it can be deduced that the special case of the model can have one or three endemic equilibria when the associated control reproduction exceeds one. Table 4 also shows that the special case of the model does not have multiple endemic equilibria when the control reproduction number

R_{vs}

is less than one, thus ruling out the existence of backward bifurcation (a dynamic behavior characterized by the co-existence of multiple stable attractors (one disease-free and one endemic) when the reproduction number of the model is less than unity) in the special case of the model. The presence of backward bifurcation in the transmission dynamics of an infectious disease makes its effective control or elimination more difficult (Alsaleh and Gumel 2014; Brauer 2004; Dushoff et al. 1998; Gumel 2012; Malik et al. 2013; Sharomi et al. 2007). This result (for the non-existence of backward bifurcation in the model (2.5) when the disease-induced mortality is negligible) is consistent with the result in (Gumel 2012) (Theorem 6, see also (Sharomi et al. 2007)).

Table 4 Number of possible positive real roots of the polynomial (2.20) for values of the control reproduction number, $R_{vs}$ , greater, equal to, or less than one

Case		$a_{4}$	$a_{3}$	$a_{2}$	$a_{1}$	$a_{0}$	Number of sign changes	Number of possible positive real roots
$R_{vs} > 1$	1	$+$	$+$	$+$	$+$	−	1	1
	2	$+$	$+$	$+$	−	−	1	1
	3	$+$	$+$	−	$+$	−	3	1 or 3
	4	$+$	$+$	−	−	−	1	1
	5	$+$	−	$+$	$+$	−	3	1 or 3
	6	$+$	−	$+$	−	−	3	1 or 3
	7	$+$	−	−	$+$	−	3	1 or 3
	8	$+$	−	−	−	−	1	1
$R_{vs} = 1$	9	$+$	$+$	$+$	$+$	0	0	0
$R_{vs} < 1$	10	$+$	$+$	$+$	$+$	$+$	0	0

In order to explore the local asymptotic stability of an endemic equilibrium of the aforementioned special case of the model (2.5), an additional assumption of perfect vaccine efficacy ( $ε_{v} = 1$ ) is made. Although also made for mathematical tractability, this additional assumption can be justified considering the fact that the HPV vaccines currently being used in Korea (Cervarix^® and Gardasil^®) are highly efficacious, with efficacy above 90% (Cheng et al. 2020).

It is convenient to define ${\tilde{R}}_{vs} : = R_{vs} |_{ε_{v} = 1}$ . We claim the following result.

Theorem 2.6

The special case of the model (2.5) with no cancer-induced mortality in females (i.e., $δ = δ_{d} = 0$ ) and with perfect vaccine protective efficacy (i.e., $ε_{v} = 1$ ) has a unique endemic equilibrium that is locally-asymptotically stable whenever ${\tilde{R}}_{vs} > 1$ .

The proof of Theorem 2.6, based on using a Krasnasolskii sub-linearity argument (Esteva and Vargas 2000; Hethcote and Thieme 1985; Thieme 1985), is given in Appendix D. The epidemiological implication of Theorem 2.6 is that, for the special case of the model (2.5) with

δ = δ_{d} = 0

and

ε_{v} = 1

, HPV will persist in the community whenever

{\tilde{R}}_{vs} > 1

. The result of this theorem is illustrated in Figure 4, showing convergence of initial solutions to the unique endemic equilibrium when

δ = δ_{d} = 0

,

ε_{v} = 1

, and

{\tilde{R}}_{vs} > 1

.

**Fig. 4** Simulations of the special case of the model (2.5) with $δ = δ_{d} = 0$ and $ε_{v} = 1$ , showing convergence of initial solutions to the unique endemic equilibrium when ${\tilde{R}}_{vs} > 1$ . Parameter values used in the simulations are as given in Table 9, with $δ = δ_{d} = 0$ and $ε_{v} = 1$ (so that, ${\tilde{R}}_{vs} = 3.0212$ ).

Calibration and Analysis of Simplified Model: No Screening

In order to fit the model (2.5) to available data and estimate its unknown parameters, a simplified version of the model with no Pap screening and with the three CIN and INM stages lumped into one CIN and INM compartment, respectively, is considered. The resulting model, tagged as the simplified model, is given by the following deterministic system of nonlinear differential equations:

\begin{matrix} Females \{\begin{matrix} {\dot{S}}_{f} & = (1 - q_{f}) π_{f} - (λ_{mf} + ξ_{f} + μ) S_{f}, \\ {\dot{V}}_{f} & = q_{f} π_{f} + ξ_{f} S_{f} - (1 - ε_{v}) λ_{mf} V_{f} - μ V_{f}, \\ {\dot{E}}_{f} & = λ_{mf} S_{f} + (1 - ε_{v}) λ_{mf} V_{f} - (σ_{f} + μ) E_{f}, \\ {\dot{I}}_{f} & = σ_{f} E_{f} - (ψ_{f} + μ) I_{f}, \\ {\dot{P}}_{f} & = (1 - b_{f}) ψ_{f} I_{f} - (α_{f} + μ) P_{f}, \\ {\dot{Q}}_{f} & = (1 - k_{f}) α_{f} P_{f} - ({\tilde{g}}_{f} + μ) Q_{f}, \\ {\dot{C}}_{f} & = (1 - d_{f}) {\tilde{g}}_{f} Q_{f} - (r_{f} + \tilde{δ} + μ) C_{f}, \\ {\dot{R}}_{f}^{C} & = r_{f} C_{f} - μ R_{f}^{C}, \\ {\dot{R}}_{f} & = b_{f} ψ_{f} I_{f} + k_{f} α_{f} P_{f} + d_{f} {\tilde{g}}_{f} Q_{f} - μ R_{f}, \end{matrix}) \\ Males \{\begin{matrix} {\dot{S}}_{m} & = (1 - q_{m}) π_{m} - (λ_{fm} + ξ_{m} + μ) S_{m}, \\ {\dot{V}}_{m} & = q_{m} π_{m} + ξ_{m} S_{m} - (1 - ε_{v}) λ_{fm} V_{m} - μ V_{m}, \\ {\dot{E}}_{m} & = λ_{fm} S_{m} + (1 - ε_{v}) λ_{fm} V_{m} - (σ_{m} + μ) E_{m}, \\ {\dot{I}}_{m} & = σ_{m} E_{m} - (ψ_{m} + μ) I_{m}, \\ {\dot{P}}_{m} & = (1 - b_{m}) ψ_{m} I_{m} - (α_{m} + μ) P_{m}, \\ {\dot{J}}_{m} & = (1 - k_{m}) α_{m} P_{m} - ({\tilde{g}}_{m} + μ) J_{m}, \\ {\dot{C}}_{m} & = (1 - {\tilde{d}}_{m}) {\tilde{g}}_{m} J_{m} - (r_{m} + μ) C_{m}, \\ {\dot{R}}_{m}^{C} & = r_{m} C_{m} - μ R_{m}^{C}, \\ {\dot{R}}_{m} & = b_{m} ψ_{m} I_{m} + k_{m} α_{m} P_{m} + {\tilde{d}}_{m} {\tilde{g}}_{m} J_{m} - μ R_{m}, \end{matrix}) \end{matrix}

3.1

where the forces of infection,

λ_{mf}

and

λ_{fm}

, defined in (2.4), are now given by (note that the conservation law (2.3) is also applied in these expressions):

\begin{matrix} λ_{mf} = (β_{mf}, c_{mf}) [\frac{I_{m} + η_{m} E_{m} + {\tilde{θ}}_{m} (P_{m} + {\tilde{κ}}_{m} J_{m})}{N_{m}}], \\ λ_{fm} = (β_{fm}, c_{mf}) [\frac{I_{f} + η_{f} E_{f} + {\tilde{θ}}_{f} (P_{f} + {\tilde{κ}}_{f} Q_{f})}{N_{m}}] . \end{matrix}

The state variables of the simplified model (3.1) are as described in Table 1, except now that the compartments for females in the undetected and detected CIN stages (

Q_{iu} (t)

and

Q_{id} (t)

, with

i = 1, 2, 3

) are lumped into the

Q_{f} (t)

compartment of females having CIN. Similarly, the populations

C_{u} (t)

and

C_{d} (t)

, of females having undetected and detected cervical cancer, are now lumped into the compartment

C_{f} (t)

of females with cervical cancer. The compartments for males in the three INM stages,

J_{im} (t), i = 1, 2, 3

, are also lumped into the compartment

J_{m} (t)

of males having INM. The parameters of the simplified model (3.1) are as described in Tables 2 and 3, with a few others described in Table 5.

Table 5 Description of parameters of the simplified model (3.1) with modified description.

Parameter	Description
${\tilde{θ}}_{f} ({\tilde{θ}}_{m})$	Modification parameter for infectiousness of $P_{f}, Q_{f} (P_{m}, J_{m})$ , in comparison to the infectiousness of $I_{f} (I_{m})$
${\tilde{κ}}_{f} ({\tilde{κ}}_{m})$	Modification parameter for infectiousness of $Q_{f} (J_{m})$ in comparison to $P_{f} (P_{m})$
${\tilde{g}}_{f} ({\tilde{g}}_{m})$	Rate at which females (males) leave the CIN (INM) class
$r_{f}$	Recovery rate of cervical cancer for females
$\tilde{δ}$	Mortality rate of females due to cervical cancer
$d_{f}$	Proportion of females with CIN that recover
$1 - d_{f}$	Proportion of females with CIN that progress to cervical cancer
${\tilde{d}}_{m}$	Proportion of males with INM that recover
$1 - {\tilde{d}}_{m}$	Proportion of males with INM that progress to HPV-related cancer

Table 1 Description of the state variables of the model (2.5).

State variable ( $i = 1, 2, 3$ )	Description ( $i = 1, 2, 3$ )
$S_{f} (S_{m})$	Number of unvaccinated susceptible females (males)
$V_{f} (V_{m})$	Number of vaccinated susceptible females (males)
$E_{f} (E_{m})$	Number of exposed/pre-symptomatic infectious females (males)
$I_{f} (I_{m})$	Number of symptomatic infectious females (males)
$P_{f} (P_{m})$	Number of females (males) with persistent infection
$Q_{iu}$	Number of females having undetected CIN of grade i
$Q_{id}$	Number of females having detected CIN of grade i
$J_{im}$	Number of males having INM of grade i
$C_{u}$	Number of females having undetected cervical cancer
$C_{d}$	Number of females having detected cervical cancer
$C_{m}$	Number of males having HPV-related cancers
$R_{f}^{C} (R_{m}^{C})$	Number of females (males) recovered from cancer
$R_{f} (R_{m})$	Number of females (males) recovered from HPV infection before developing cancer

Table 2 Description of the parameters of the model (2.5).

Parameter	Description ( $i = 1, 2, 3$ )
$π_{f} (π_{m})$	Recruitment rate of sexually active females (males)
$q_{f} (q_{m})$	Proportion of new sexually active females (males) vaccinated
$ξ_{f} (ξ_{m})$	Catch-up vaccination rate of females (males)
$1 / μ$	Average duration of sexual activity
$c_{mf} (c_{fm})$	Average number of sexual partners of females (males) per year
$β_{mf} (β_{fm})$	Transmission probability from infected males (females) to susceptible females (males) per partnership
$η_{f} (η_{m})$	Modification parameter for infectiousness of $E_{f} (E_{m})$ in comparison to $I_{f} (I_{m})$
$θ_{f} (θ_{m})$	Modification parameter for infectiousness of $P_{f}, Q_{iu}, Q_{id} (P_{m}, J_{im})$ in comparison to $I_{f} (I_{m})$
$κ_{f} (κ_{m})$	Modification parameter for infectiousness of $Q_{2 u}, Q_{3 u} (J_{2 m}, J_{3 m})$ in comparison to $P_{f}, Q_{1 u} (P_{m}, J_{1 m})$
$ν$	Modification parameter for infectiousness of $Q_{1 d}, Q_{2 d}, Q_{3 d}$ in comparison to $P_{f}, Q_{iu}$
$κ_{fd}$	Modification parameter for infectiousness of $Q_{2 d}, Q_{3 d}$ in comparison to $Q_{1 d}$
$ε_{v}$	Efficacy of bivalent Cervarix^® (quadrivalent Gardasil^®) vaccine
$σ_{f} (σ_{m})$	Rate of symptom development of females (males)
$ψ_{f} (ψ_{m})$	Rate at which females (males) leave the symptomatic infectious class
$α_{f} (α_{m})$	Rate at which females (males) leave the persistent infection class
$g_{f}, y_{f}, z_{f}$	Rate at which females leave the undetected CIN1, CIN2, CIN3 class
$ζ_{1}, ζ_{2}, ζ_{3}$	Rate at which females leave the detected CIN1, CIN2, CIN3 class
$g_{m}, y_{m}, z_{m}$	Rate at which males leave the INM1, INM2, INM3 class
$w_{f}$	Detection rate of cancer in females
$r_{u} (r_{d})$	Recovery rate of undetected (detected) cancer for females
$r_{m}$	Recovery rate of cancer for males
$δ (δ_{d})$	Mortality rate for females with undetected (detected) cervical cancer

Table 3 Epidemiological proportions associated with some of the parameters of the model (2.5) described in Table 2.

Parameter	Description
$b_{f} (b_{m})$	Proportion of symptomatic infectious females (males) that recover
$1 - b_{f}$	Proportion of symptomatic infectious females that progress to persistent infection
$1 - b_{m}$	Proportion of symptomatic infectious males that progress to persistent infection
$k_{f} (k_{m})$	Proportion of persistently infected females (males) that recover
$1 - k_{f}$	Proportion of persistently infected females that progress to undetected CIN1
$1 - k_{m}$	Proportion of persistently infected males that progress to INM1
$d_{u 1}$	Proportion of females with undetected CIN1 that recover naturally
$d_{u 3}$	Proportion of females with undetected CIN1 that get detected for CIN
$1 - d_{u 1} - d_{u 3}$	Proportion of females with undetected CIN1 that progress to undetected CIN2
$h_{u 1}$	Proportion of females with undetected CIN2 that recover naturally
$h_{u 2}$	Proportion of females with undetected CIN2 that regress to undetected CIN1
$h_{u 3}$	Proportion of females with undetected CIN2 that get detected for CIN
$1 - h_{u 1} - h_{u 2} - h_{u 3}$	Proportion of females with undetected CIN2 that progress to undetected CIN3
$q_{u 1}$	Proportion of females with undetected CIN3 that recover naturally
$q_{u 2}$	Proportion of females with undetected CIN3 that regress to undetected CIN2
$q_{u 3}$	Proportion of females with undetected CIN3 that get detected for CIN
$1 - q_{u 1} - q_{u 2} - q_{u 3}$	Proportion of females with undetected CIN3 that progress to undetected cancer
$d_{d 1}$	Proportion of females with detected CIN1 that recover
$1 - d_{d 1}$	Proportion of females with detected CIN1 that progress to detected CIN2
$h_{d 1}$	Proportion of females with detected CIN2 that recover
$h_{d 2}$	Proportion of females with detected CIN2 that regress to detected CIN1
$1 - h_{d 1} - h_{d 2}$	Proportion of females with detected CIN2 that progress to detected CIN3
$q_{d 1}$	Proportion of females with detected CIN3 that recover
$q_{d 2}$	Proportion of females with detected CIN3 that regress to detected CIN2
$1 - q_{d 1} - q_{d 2}$	Proportion of females with detected CIN3 that progress to detected cancer
$d_{a}$	Proportion of males with INM1 that recover
$1 - d_{m}$	Proportion of males with INM1 that progress to INM2
$h_{m 1}$	Proportion of males with INM2 that recover
$h_{m 2}$	Proportion of males with INM2 that regress to INM1
$1 - h_{m 1} - h_{m 2}$	Proportion of males with INM2 that progress to INM3
$q_{m 1}$	Proportion of males with INM3 that recover
$q_{m 2}$	Proportion of males with INM3 that regress to INM2
$1 - q_{m 1} - q_{m 2}$	Proportion of males with INM3 that progress to HPV-related cancer

Here, too, it can be shown (using the approach in Section 2.1) that all non-negative solutions of the simplified model (3.1) stay non-negative for all time t. Furthermore, using the approach in Section 2.1, the following region

\begin{matrix} D = D_{f} \times D_{m} \subset R_{+}^{9} \times R_{+}^{9}, \end{matrix}

with,

D_{f} = {(S_{f}, V_{f}, E_{f}, I_{f}, P_{f}, Q_{f}, C_{f}, R_{f}^{C}, R_{f}) \in R_{+}^{9} : 0 \leq N_{f} \leq \frac{π_{f}}{μ}}

and

D_{m} = {(S_{m}, V_{m}, E_{m}, I_{m}, P_{m}, J_{m}, C_{m}, R_{m}^{C}, R_{m}) \in R_{+}^{9} : 0 \leq N_{m} \leq \frac{π_{m}}{μ}}

, can be shown to be positively invariant with respect to the simplified model (3.1). Hence, the simplified model is well-posed in the region

D

. The disease-free equilibrium of the simplified model (3.1) is given by:

\begin{matrix} E_{0 n} = (S_{f}^{*}, V_{f}^{*}, 0, 0, 0, 0, 0, 0, 0, S_{m}^{*}, V_{m}^{*}, 0, 0, 0, 0, 0, 0, 0), \end{matrix}

where, now,

\begin{matrix} S_{f}^{*} = \frac{(1 - q_{f}) π_{f}}{ξ_{f} + μ}, V_{f}^{*} = \frac{π_{f} (ξ_{f} + μ q_{f})}{μ (ξ_{f} + μ)}, S_{m}^{*} = \frac{(1 - q_{m}) π_{m}}{ξ_{m} + μ}, \\ V_{m}^{*} = \frac{π_{m} (ξ_{m} + μ q_{m})}{μ (ξ_{m} + μ)} . \end{matrix}

3.2

It is convenient to define the following quantity (the control reproduction number for the simplified model):

\begin{matrix} R_{vn} = \sqrt{R_{fn} R_{mn}}, \end{matrix}

3.3

where,

\begin{matrix} R_{fn} & = \frac{β_{fm} c_{mf} μ p_{m}^{*}}{π_{m} (σ_{f} + μ)} [η_{f} + \frac{σ_{f}}{ψ_{f} + μ} + \frac{{\tilde{θ}}_{f} σ_{f} (1 - b_{f}) ψ_{f}}{(ψ_{f} + μ) (α_{f} + μ)}) \\ (+ \frac{{\tilde{θ}}_{f} {\tilde{κ}}_{f} σ_{f} (1 - b_{f}) ψ_{f} (1 - k_{f}) α_{f}}{(ψ_{f} + μ) (α_{f} + μ) ({\tilde{g}}_{f} + μ)}], \\ R_{mn} & = \frac{β_{mf} c_{mf} μ p_{f}^{*}}{π_{m} (σ_{m} + μ)} [η_{m} + \frac{σ_{m}}{ψ_{m} + μ} + \frac{{\tilde{θ}}_{m} σ_{m} (1 - b_{m}) ψ_{m}}{(ψ_{m} + μ) (α_{m} + μ)}) \\ (+ \frac{{\tilde{θ}}_{m} {\tilde{κ}}_{m} σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m}}{(ψ_{m} + μ) (α_{m} + μ) ({\tilde{g}}_{m} + μ)}], \end{matrix}

with,

p_{f}^{*} = \frac{π_{f}}{μ} (1 - ε_{v} \frac{ξ_{f} + μ q_{f}}{ξ_{f} + μ})

and

p_{m}^{*} = \frac{π_{m}}{μ} (1 - ε_{v} \frac{ξ_{m} + μ q_{m}}{ξ_{m} + μ})

. The basic reproduction number of the simplified model is given by:

\begin{matrix} R_{0 n} = R_{vn} |_{q_{f} = ξ_{f} = q_{m} = ξ_{m} = 0} . \end{matrix}

3.4

The result below follows from Theorem 2 of (Van den Driessche and Watmough 2002) (see Appendix E.1 for detailed proof).

Theorem 3.1

The DFE, $E_{0 n}$ , of the simplified model (3.1) is locally-asymptotically stable if $R_{vn} < 1$ and unstable if $R_{vn} > 1 .$

This local asymptotic stability result (Theorem 3.1) is extended to the global asymptotic stability of the DFE of the special case of the simplified model with negligible disease-induced mortality (i.e., the model (3.1) with $\tilde{δ} = 0$ ). It is convenient to define the DFE of the special case of the simplified model (3.1) by ${\tilde{E}}_{0 n} = E_{0 n} |_{\tilde{δ} = 0}$ . The following result will be needed in proving the global asymptotic stability of the DFE, ${\tilde{E}}_{0 n}$ .

Lemma 3.2

The region

\begin{matrix} D^{*} = & {(S_{f}, V_{f}, E_{f}, I_{f}, P_{f}, Q_{f}, C_{f}, R_{f}^{C}, R_{f}, S_{m}, V_{m}, E_{m}, I_{m}, P_{m}, J_{m}, C_{m}, R_{m}^{C}, R_{m}) \in D : \\ S_{f} \leq S_{f}^{*}, V_{f} \leq V_{f}^{*}, S_{m} \leq S_{m}^{*}, V_{m} \leq V_{m}^{*}} \end{matrix}

is positively-invariant, and attracts all solutions of the simplified model (3.1) with

\tilde{δ} = 0

.

Proof

Consider the special case of the simplified model (3.1) with

\tilde{δ} = 0

. Substituting

\tilde{δ} = 0

into the equations for the rate of change of the total populations of females and males of the simplified model shows that

N_{i} (t) \to \frac{π_{i}}{μ}

as

t \to \infty

, for

i = f, m

, respectively. Thus, from now on,

N_{i} (t)

will be replaced by their respective disease-free equilibrium values (or the respective maximum total populations in the absence of disease),

\frac{π_{i}}{μ}

(for

i = f, m

). Then, it follows from the first and 10th equations of the simplified model that for

i, j \in {f, m}

and

i \neq j

,

\begin{matrix} {\dot{S}}_{i} & = (1 - q_{i}) π_{i} - (λ_{ji} + ξ_{i} + μ) S_{i} \\ \leq (1 - q_{i}) π_{i} - (ξ_{i} + μ) S_{i}, since λ_{ji} \geq 0 \\ = (ξ_{i} + μ) (S_{i}^{*} - S_{i}), \end{matrix}

so that, by a standard comparison theorem (Lakshmikantham et al. 2016),

S_{i} (t) \leq S_{i} (0) e^{- (ξ_{i} + μ) t} + S_{i}^{*} [1 - e^{- (ξ_{i} + μ) t}]

. Thus,

\begin{matrix} S_{i} (t) \leq S_{i}^{*} for all t \geq 0 provided that S_{i} (0) \leq S_{i}^{*} . \end{matrix}

3.5

Similarly, from the second and 11th equations of the simplified model (noting that

N_{i} (t) \leq \frac{π_{i}}{μ}

in

D

):

\begin{matrix} {\dot{V}}_{i} & = q_{i} π_{i} + ξ_{i} S_{i} - (1 - ε_{v}) λ_{ji} V_{i} - μ V_{i} \\ \leq q_{i} π_{i} + ξ_{i} (\frac{π_{i}}{μ} - V_{i} - E_{i} - I_{i} - P_{i} - Q_{i} - C_{i} - R_{i}^{C} - R_{i}) - μ V_{i} \\ \leq \frac{(μ q_{i} + ξ_{i}) π_{i}}{μ} - (ξ_{i} + μ) V_{i} \\ = (ξ_{i} + μ) (V_{i}^{*} - V_{i}) . \end{matrix}

Hence, by the comparison theorem,

V_{i} (t) \leq V_{i} (0) e^{- (ξ_{i} + μ) t} + V_{i}^{*} [1 - e^{- (ξ_{i} + μ) t}],

from which it follows that

\begin{matrix} V_{i} (t) \leq V_{i}^{*} for all t \geq 0 provided that V_{i} (0) \leq V_{i}^{*} . \end{matrix}

3.6

Hence, based on the bounds in (3.5) and (3.6), the region

D^{*}

is positively-invariant with respect to the flow generated by the special case of the simplified model (3.1) with

\tilde{δ} = 0

. In fact, it can also be shown that each solution in

D

either enters

D^{*}

in finite time or limits to

{\tilde{E}}_{0 n}

(Gumel et al. 2006).

□

We claim the following result (the proof, which is based on using a standard comparison theorem (Lakshmikantham et al. 2016), is given in Appendix E.2):

Theorem 3.3

The DFE, ${\tilde{E}}_{0 n}$ , of the special case of the simplified model (3.1) with negligible disease-induced mortality (i.e., $\tilde{δ} = 0$ ) is globally-asymptotically stable in $D^{*}$ whenever $R_{vn} < 1$ .

Unlike in the case of Theorem 3.1 (where the effective control of the disease when $R_{vn} < 1$ depends on the initial number of infected individuals being within the basin of attraction of the DFE of the simplified model), the epidemiological implication of Theorem 3.3 is that bringing the control reproduction number of the aforementioned special case of the simplified model below one is necessary and sufficient for the effective control of the disease, regardless of the size of the initial number of infected individuals.

Data fitting and parameter estimation for the simplified model

Although the values of most of the parameters of the full model (2.5) are reliably known from the literature (as tabulated in Table 9) and demographic data from Korea, the values of some of the parameters are not known. These unknown parameters will be estimated by fitting the simplified version of the model (3.1) to the observed cumulative HPV-related cancer cases and mortality in the Republic of Korea, for both females and males, for the period from 1999 to 2020 (National Cancer Center Korea Central Cancer Registry 2022). For fitting purposes, it is assumed that males are not vaccinated (i.e., $q_{m} = ξ_{m} = 0$ ) and the vaccination coverage of females remains low ( $q_{f} = 0.05$ per year, $ξ_{f} = 0.01$ per year (Choi et al. 2022)) before South Korea started the National Immunization Program (NIP) on HPV in 2016. Furthermore, since HPV is not the only cause of some cancer types (de Martel et al. 2017), an estimate for the average number of cancer cases attributable to HPV (needed for the data fitting) is obtained by taking the weighted average of cancer cases from different cancer sites with their HPV-attributable fraction given in Table 6.

Table 6 HPV-attributable fraction by cancer sites (de Martel et al. 2017)

Cancer sites (ICD-10 code)	HPV-attributable fraction (%)	Cancer sites (ICD-10 code)	HPV-attributable fraction (%)
Cervix uteri (C53)	100	Oropharynx (C01, C09-10)	30.8
Anus (C21)	88	Vulva (C51)	24.9
Vagina (C52)	78	Larynx (C32)	2.4
Penis (C60)	50	Oral cavity (C02-06)	2.2

It should be noted that, although the parameters related to disease progression and recovery in females are well-known in the literature, the corresponding parameters for males are unknown. In other words, although the values of the female progression and recovery parameters (i.e.,

(1 - b_{f}) ψ_{f}, (1 - k_{f}) α_{f}

,

(1 - d_{f}) {\tilde{g}}_{f}

,

b_{f} ψ_{f}, k_{f} α_{f}

, and

d_{f} {\tilde{g}}_{f}

) are known (and their baseline values are tabulated in Table 8), the values of the male-related progression and recovery parameters (i.e.,

(1 - b_{m}) ψ_{m}, (1 - k_{m}) α_{m}, (1 - {\tilde{d}}_{m}) {\tilde{g}}_{m}

,

b_{m} ψ_{m}, k_{m} α_{m}

, and

{\tilde{d}}_{m} {\tilde{g}}_{m}

) are unknown. These unknown male-related parameters will be estimated by fitting for the ratio of these rates for males, in relation to those for females (this is to take advantage of the fact that males recover faster from HPV infection and are less likely to have persistent infection than females (Giuliano et al. 2008a; National Cancer Center Korea Central Cancer Registry 2022; Moscicki et al. 2004, 1998; Shin et al. 2004; Van Doornum et al. 1994)). To do this, we assume that the ratio of the rate at which males progress from the

I_{m}

class to the

P_{m}

class (i.e.,

(1 - b_{m}) ψ_{m}

) to the rate at which females progress from the

I_{f}

class to the

P_{f}

class (i.e.,

(1 - b_{f}) ψ_{f}

), defined as

r_{1} \in (0, 1]

, equals the ratio for the progression rate of males from the

P_{m}

class to the

J_{m}

compartment (i.e.,

(1 - k_{m}) α_{m}

) to the progression rate of females from the

P_{f}

compartment to the

Q_{f}

compartment (i.e.,

(1 - k_{f}) α_{f}

). The justification for

r_{1} \in (0, 1]

stems from the fact that HPV prevalence in males is lower than in females (National Cancer Center Korea Central Cancer Registry 2022; Shin et al. 2004; Van Doornum et al. 1994). Similarly, we assume that this ratio is the same as that for the progression rate of males from the

J_{m}

class to the

C_{m}

class (i.e.,

(1 - {\tilde{d}}_{m}) {\tilde{g}}_{m}

) to that for the progression rate of females from the

Q_{f}

class to the

C_{f}

class (i.e.,

(1 - d_{f}) {\tilde{g}}_{f}

). That is, we set

\begin{matrix} \frac{(1 - b_{m}) ψ_{m}}{(1 - b_{f}) ψ_{f}} = \frac{(1 - k_{m}) α_{m}}{(1 - k_{f}) α_{f}} = \frac{(1 - {\tilde{d}}_{m}) {\tilde{g}}_{m}}{(1 - d_{f}) {\tilde{g}}_{f}} = : r_{1} . \end{matrix}

3.7

The objective then is to use the observed data for cumulative cancer cases and cancer mortality in South Korea (National Cancer Center Korea Central Cancer Registry 2022) to fit the simplified model for

r_{1}

(this entails replacing every occurrence of the terms in the numerators of (3.7) in the simplified model (3.1) with their respective product of their denominators with

r_{1}

). In other words, once the fitted value of

r_{1}

is obtained, equation (3.7) is then be used, together with the baseline values of the female-related parameters in Table 8, to obtain the estimated values of the three unknown male-related progression parameters (namely,

(1 - b_{m}) ψ_{m}, (1 - k_{m}) α_{m}

, and

(1 - {\tilde{d}}_{m}) {\tilde{g}}_{m}

). Similarly, it is assumed that the ratio of the recovery rate of males from the

I_{m}

class (i.e.,

b_{m} ψ_{m}

) to that of females from the

I_{f}

class (i.e.,

b_{f} ψ_{f}

), defined as

r_{2} > 1

, equals that for the recovery rate of males from the

P_{m}

class (i.e.,

k_{m} α_{m}

) to that of females from the

P_{f}

class (i.e.,

k_{f} α_{f}

). This ratio is also assumed to be the same as that for the recovery rate of males from the

J_{m}

class (i.e.,

{\tilde{d}}_{m} {\tilde{g}}_{m}

) to that of females from the

Q_{f}

class (i.e.,

d_{f} {\tilde{g}}_{f}

). In other words, we set

\begin{matrix} \frac{b_{m} ψ_{m}}{b_{f} ψ_{f}} = \frac{k_{m} α_{m}}{k_{f} α_{f}} = \frac{{\tilde{d}}_{m} {\tilde{g}}_{m}}{d_{f} {\tilde{g}}_{f}} = : r_{2} . \end{matrix}

3.8

Here, too, once the fitted value of

r_{2}

is obtained, equation (3.8) can then be used, along with the baseline values of the female-related parameters in Table 8, to obtain estimated values of the three unknown male-related recovery parameters (namely,

b_{m} ψ_{m}, k_{m} α_{m}

, and

{\tilde{d}}_{m} {\tilde{g}}_{m}

).

Table 8 Baseline values of the parameters for the simplified model (3.1)

Parameter	Baseline Value (per year)	Source / Comment
$π_{f}$ , $π_{m}$	358000, 372000	(Statistics Korea 2000)
$μ$	1/65	(Malik et al. 2013)
$q_{f}$ , $ξ_{f}$	0.05 (dimensionless), 0.01	(Choi et al. 2022)
$q_{m}$ , $ξ_{m}$	0 (dimensionless), 0	Assumed
$c_{mf}$	2	(Podder and Gumel 2009)
$c_{fm}$	$\frac{c_{mf} N_{f} (t)}{N_{m} (t)} = \frac{2 N_{f} (t)}{N_{m} (t)}$	From the value of $c_{mf}$ above
$β_{mf}$	0.72 (dimensionless)	Derived from (Podder and Gumel 2009)
$β_{fm}$	$\frac{β_{fm} c_{fm}}{c_{fm} (t)} = \frac{2.5}{c_{fm} (t)}$ (dimensionless)	Derived from (Podder and Gumel 2009)
$η_{f}, η_{m}$	0.5 (dimensionless)	(Javame and Gumel 2016)
${\tilde{θ}}_{f}, {\tilde{θ}}_{m}$	0.9 (dimensionless)	(Javame and Gumel 2016)
${\tilde{κ}}_{f}, {\tilde{κ}}_{m}$	1.5 (dimensionless)	(Javame and Gumel 2016)
$ε_{v}$	0.95 (dimensionless)	(Malik et al. 2013)
$σ_{f}, σ_{m}$	3	(Javame and Gumel 2016)
$(1 - b_{f}) ψ_{f}$	0.017	(Javame and Gumel 2016)
$(1 - b_{m}) ψ_{m} = r_{1} (1 - b_{f}) ψ_{f}$	0.0076	(Javame and Gumel 2016) and fitted value of $r_{1}$
$(1 - k_{f}) α_{f}$	0.07	(Javame and Gumel 2016)
$(1 - k_{m}) α_{m} = r_{1} (1 - k_{f}) α_{f}$	0.03	(Javame and Gumel 2016) and fitted value of $r_{1}$
$(1 - d_{f}) {\tilde{g}}_{f}$	0.14	(Malik et al. 2013)
$(1 - {\tilde{d}}_{m}) {\tilde{g}}_{m} = r_{1} (1 - d_{f}) {\tilde{g}}_{f}$	0.06	(Malik et al. 2013) and fitted value of $r_{1}$
$b_{f} ψ_{f}$	0.495	(Malik et al. 2013)
$b_{m} ψ_{m} = r_{2} b_{f} ψ_{f}$	0.625	(Malik et al. 2013) and fitted value of $r_{2}$
$k_{f} α_{f}$	0.05	(Javame and Gumel 2016)
$k_{m} α_{m} = r_{2} k_{f} α_{f}$	0.063	(Javame and Gumel 2016) and fitted value of $r_{2}$
$d_{f} {\tilde{g}}_{f}$	0.01	(Javame and Gumel 2016)
${\tilde{d}}_{m} {\tilde{g}}_{m} = r_{2} d_{f} {\tilde{g}}_{f}$	0.013	(Javame and Gumel 2016) and fitted value of $r_{2}$
$r_{f}$	0.5	(Javame and Gumel 2016)
$r_{m} = r_{2} r_{f}$	0.63	(Javame and Gumel 2016) and fitted value of $r_{2}$
$\tilde{δ}$	0.23	(Ministry of Health and Welfare 2021)

Thus, the unknown parameters (or combination of parameters) of the simplified model (3.1) (and, by extension, the unknown parameters or combination of parameters for the full model (2.5)) to be estimated from the data fitting are:

the effective transmission rate from infected males to susceptible females ( $c_{mf} β_{mf})$ ;
the effective transmission rate from infected females to susceptible males ( $c_{fm} β_{fm}$ );
the ratio ( $r_{1} \in (0, 1]$ ) of the progression rates of males (from $I_{m}$ to $P_{m}$ , $P_{m}$ to $J_{m}$ , and from $J_{m}$ to $C_{m}$ ) in comparison to progression rates of females (from $I_{f}$ to $P_{f}$ , $P_{f}$ to $Q_{f}$ , and from $Q_{f}$ to $C_{f}$ ). This allows us to estimate $(1 - b_{m}) ψ_{m}, (1 - k_{m}) α_{m}$ , and $(1 - {\tilde{d}}_{m}) {\tilde{g}}_{m}$ ;
the ratio ( $r_{2} > 1$ ) of the recovery rates of males (from $I_{m}$ , $P_{m}$ , and $J_{m}$ to $R_{m}$ ) in comparison to recovery rates of females (from $I_{f}$ , $P_{f}$ , and $Q_{f}$ to $R_{f}$ ). This allows us to estimate $b_{m} ψ_{m}, k_{m} α_{m}$ , and ${\tilde{d}}_{m} {\tilde{g}}_{m}$ .

The simplified model (3.1) is fitted to the data using a standard nonlinear least squares approach. Specifically, the "lsqcurvefit" function in-built in MATLAB is used to minimize the sum of the squared differences between each observed data point (for the cumulative cancer cases and mortality) and the corresponding cumulative cancer cases and mortality projection obtained from the simulation of the simplified model (3.1). The baseline values of the aforementioned eight unknown parameters obtained from the data fitting, together with their corresponding 95% confidence intervals, are tabulated in Table 7.

Table 7 Baseline values of estimated (fitted) parameters for the simplified model (3.1)

Parameter	Estimated Value	$95 %$ CI
$β_{mf} c_{mf}$	1.45 per year	(0.91, 1.48)
$β_{fm} c_{fm}$	2.5 per year	(2.49, 5.57)
$r_{1}$	0.42 (dimensionless)	(0.38, 0.44)
$r_{2}$	1.24 (dimensionless)	(1.21, 1.58)

The results of the data fitting of the simplified model is depicted in Figure 5, showing a good fit ( $R^{2} = 0.99997$ ). The goodness of the fitting is further corroborated by simulating the simplified model (3.1) using the estimated and fixed parameters (in Tables 7 and 8) to predict the yearly cancer cases and mortality in both females and males, and compare these with the corresponding observed data. The results obtained, depicted in Figure 6, also show a very good prediction. It should be stated that, using the baseline values of the fixed and estimated parameters given in Tables 7 and 8, Equations (3.3) and (3.4) show that the vaccination reproduction number ( $R_{vn}$ ) and the basic reproduction number ( $R_{0 n}$ ) of the simplified model (3.1) take the values $R_{vn} = 3.1206$ , and $R_{0 n} = 4.0389$ , respectively. Thus, this study shows that the pre-NIP vaccination coverage level (which is what we used to fit the simplified model) will not lead to the effective control of HPV in Korea (since it fails to bring $R_{vn}$ to a value less than one).

**Fig. 5** Data fitting of the simplified model (3.1), using cumulative data on the incidence and mortality of HPV-related cancers in females and males in Korea from 1999 to 2020 (National Cancer Center Korea Central Cancer Registry 2022). (a) Cumulative cervical cancer cases in females, (b) Cumulative cervical cancer mortality in females, (c) Cumulative HPV-related cancer cases in males. The first 17 years of data are used for data fitting, the next 5 years of data are used for cross-validation, and the calibrated model (3.1) predicts the next 5 years.

**Fig. 6** Simulation of the calibrated model (3.1) using parameters estimated from the cumulative data on the incidence and mortality of HPV-related cancers in females and males in Korea (National Cancer Center Korea Central Cancer Registry 2022). The fixed and fitted parameter values used are as given in Table 8. The shaded regions represent 95% confidence bands propagated from the bootstrap samples of fitted parameters in Table 7. (a) Annual cervical cancer cases in females, (b) Annual cervical cancer mortality in females, (c) Annual HPV-related cancer cases in males.

The baseline values for the fixed and estimated parameters of the simplified model (given in Tables 7 and 8) are used to acquire the baseline values of the remaining parameters of the full model (2.5), as tabulated in Table 9. For instance, the detection rates of females with CIN1, CIN2, and CIN3 (denoted by $d_{u 3} g_{f}$ , $h_{u 3} y_{f}$ , and $q_{u 3} z_{f}$ , respectively) are obtained by multiplying the Pap screening coverage (assumed to be 61% (Shin et al. 2022)) with sensitivity of the Pap test (assumed to be 60% for CIN1 and 78% for CIN2/CIN3 (Van de Velde et al. 2007)), and dividing by the frequency of the Pap test (which, in Korea, is recommended to be every two years for females of age 20 to 74 (Shin et al. 2022)).

Sensitivity Analysis for the Model (2.5)

The model (2.5) contains 62 parameters. Although the baseline values of these parameters are known (either available directly from the literature, or are obtained by fitting the simplified model (3.1) to the available data, as tabulated in Table 9), uncertainties are expected to arise in the estimate of their values. Therefore, it is crucial to evaluate the impact of these uncertainties on the overall outcome of the numerical simulations of the model, and determine which parameters have the most influence on a chosen response function (Blower and Dowlatabadi 1994; Marino et al. 2008; Safdar et al. 2022). The sensitivity analysis entails defining distributions of each parameter and finding the partial rank correlation coefficients (PRCCs) for each parameter in relation to the response function (Blower and Dowlatabadi 1994; Marino et al. 2008; Safdar et al. 2022). PRCC values range from $- 1$ to 1, and a higher magnitude of the PRCC value (typically greater or equal to 0.5 in magnitude) implies a stronger correlation with the response function (Safdar et al. 2022). In addition, parameters with positive PRCC values are positively correlated, while those with negative values are negatively correlated with the response function. To implement the method, a response function is chosen first of all (specifically, the control reproduction number $R_{vs}$ is chosen as the response function), and distribution of each parameter is chosen as uniform distribution with range defined by taking 20% to the left and 20% to the right of the baseline values given in Table 9 (Safdar et al. 2022). It should be mentioned that, since the response function (i.e., the control reproduction number, $R_{vs}$ ) contains 56 of the 62 parameters of the model (2.5), the remaining six parameters of the model (namely, $w_{f}, r_{u}, r_{d}, r_{m}, δ, δ_{d}$ ) do not feature in the sensitivity analysis with respect to $R_{vs}$ . Moreover, the factor fixing approach (Saltelli 2008) was used to determine the parameters that are not influential to the value of $R_{vs}$ . This analysis revealed that the value of $R_{vs}$ is not sensitive to 31 of the 56 parameters in $R_{vs}$ . Consequently, the sensitivity analysis was conducted for the remaining 25 parameters, after the monotonicity of $R_{vs}$ was confirmed in each of them by generating respective scatter plots of $R_{vs}$ as a function of each of the parameters (as depicted in Figure 12 in Appendix F). Specifically, the parameter ranges are divided into 1, 000 sub-intervals of equal length, and parameter sets are chosen from the sub-intervals without replacement. This leads to a $1, 000 \times 25$ matrix, and each row of the matrix is used to calculate the response function $R_{vs}$ . Then, the Pearson correlation coefficient is calculated for each parameter (which implies the degree of influence of each parameter on the response function in the absence of the influence of the other parameters) (Blower and Dowlatabadi 1994; Marino et al. 2008).

**Fig. 12** Scatter plots of the control reproduction number, $R_{vs}$ , of the model (2.5), with respect to each of the parameters in $R_{vs}$ included in the sensitivity analysis. For each plot, the other parameter values used are as given in Table 9.

Figure 7 depicts the PRCC values of the most influential (i.e., PRCC value greater or equal to 0.5 in magnitude) parameters of the response function obtained from the sensitivity analysis. This figure shows that the top-six parameters that have the highest impact on the response function (

R_{vs}

) are:

the average number of sexual partners for females per year ( $c_{mf}$ ; with PRCC = +0.997). It should be stated that the parameter for the average number of sexual partners for males per year ( $c_{fm}$ ) is also highly correlated with the response function, since it is proportional to $c_{mf}$ , by the conservation law (2.3).
the transmission probability of HPV from infected females to susceptible males ( $β_{fm}$ ; with PRCC = +0.976).
the recovery rate of symptomatic infectious females from HPV infection ( $b_{f} ψ_{f}$ ; with PRCC $= - 0.959$ ).
the catch-up vaccination rate of females ( $ξ_{f}$ ; with PRCC $= - 0.919$ ).
the recovery rate of symptomatic infectious males from HPV infection ( $b_{m} ψ_{m}$ ; with PRCC $= - 0.918$ ).
the transmission probability of HPV from infected males to susceptible females ( $β_{mf}$ ; with PRCC $= + 0.917$ ).

These sensitivity analysis results indicate that, with respect to the reproduction number (

R_{vs}

) as the response function, the parameters related to HPV transmission (i.e.,

c_{mf}, c_{fm}, β_{fm}, β_{mf}

), recovery of symptomatic infectious individuals (i.e.,

b_{f} ψ_{f}, b_{m} ψ_{m}

), and vaccination (i.e.,

ξ_{f}

) have the most influence on the control reproduction number

R_{vs}

(hence, on the burden of HPV and HPV-related cancers in Korea). Thus, public health control and mitigation efforts against HPV infection and associated cancers should be targeted towards these highly influential parameters (i.e., increasing (decreasing) their value if their PRCC is negative (positive)). Specifically:

(i) the values of the transmission-related parameters ( $c_{mf}, c_{fm}, β_{fm},$ and $β_{mf}$ ) can be decreased by implementing public health strategies that encourage safer sexual practices, such as condom use (which reduces $β_{fm}$ and $β_{mf}$ ), and reduction of average number of sexual partners (which reduces $c_{mf}$ and $c_{fm}$ ).
(ii) the values of the recovery-related parameters ( $b_{f} ψ_{f}, b_{m} ψ_{m}$ , and also $k_{f} α_{f}$ ) can be increased by increasing the testing, detection and treatment of individuals with clinical symptoms of HPV.
(iii) the values of the vaccine-related parameters ( $ξ_{f}$ , and also $q_{f}, ε_{v}$ ) can be increased by increasing the coverage of the HPV vaccine among the target population (the efficacy of both the Cervarix^® and Gardasil^® vaccine are already very high (Cheng et al. 2020)).

It should be mentioned, based on the PRCC values in Figure 7, that other parameters, such as the rate of symptom development of males (

σ_{m}

, with PRCC = -0.773), the progression rate of infectious symptomatic females to persistent infection (

(1 - b_{f}) ψ_{f}

, with PRCC = +0.671), the rate of symptom development of females (

σ_{f}

, with PRCC = -0.644), and the modification parameter for the infectiousness of pre-symptomatic males compared to symptomatic males (

η_{m}

, with PRCC = +0.540), are also highly correlated to the response function,

R_{vs}

.

**Fig. 7** Partial rank correlation coefficients (PRCCs) of the most significant coefficients (with $| PRCC | \geq 0.5$ ) in the chosen response function ( $R_{vs}$ ) of the model (2.5). The baseline values of the parameters used are as given in Table 9, and their ranges are taken to be 20% to the left and 20% to the right of their baseline values, assuming a uniform distribution.

Numerical Simulation and Assessment of Control Strategies

The full model (2.5) will now be simulated, using the estimated and fixed parameters (given in Tables 7 and 9), to assess the population-level impact of the combined vaccination and Pap screening strategies currently being implemented in an effort to effectively control the spread of HPV and related cancers in Korea. The impact of other extended versions of these strategies will also be assessed. Specifically, the following strategies will be assessed:

Strategy A: This strategy is based on the post-NIP vaccination strategy (which entails the routine annual vaccination of 80% of girls in the 12-17 year age group), the annual catch-up vaccination of 30,000 females in the 18-26 year age group (which corresponds to the baseline catch-up vaccination coverage of females, calculated using the baseline value of the catch-up vaccination rate parameter for females, $ξ_{f}$ , given in Table 9), and the routine Pap screening strategy (which involves the cervical cancer screening of 61% of adult females in the age group 20-74 (Min et al. 2015)). This strategy does not include the vaccination of boys or males.
Strategy B: This is the same as Strategy A (and boys and males are not vaccinated), but with the vaccination coverage for girls increased to 95% (as against the 80% in Strategy A) and the vaccination coverage for the 18-26 year old females increased to 175,000 (as against 30,000 in Strategy A). This increase allows the vaccination coverage of females at the disease-free equilibrium $(v_{f}^{*} = \frac{ξ_{f} + μ q_{f}}{ξ_{f} + μ})$ , obtained by solving for $v_{f}^{*}$ from equation (2.12) for the case when $v_{m}^{*} = 0$ , to reach the herd immunity threshold ( $v_{f}^{*} = 0.99$ , as described in Section 2.2.1).
Strategy C: This is the same as Strategy A, but with the annual routine vaccination of 80% of boys in the 12-17 year age group (the catch-up vaccination rate for males, $ξ_{m}$ , is set to zero in the numerical simulations).

Simulations will also be carried out for the aforementioned three strategies with the Pap screening coverage increased to 90%.

Disease burden corresponding to each strategy

The control reproduction number, $R_{vs}$ , corresponding to each of the three strategies are tabulated in Table 10, for the cases where Pap screening is at 61% baseline and when it is increased to 90%. It follows from this table that, while the reproduction numbers for Strategies B and C are below one, the reproduction number for Strategy A is above one. In other words, based on the parameter values (tabulated in Table 9 and the caption of Table 10) used in these simulations, HPV can be eliminated under Strategies B and C, but will persist under Strategy A. Thus, the increase in the annual vaccination coverage for girls (from the current post-NIP baseline of 80% to 95%) combined with the increase in the catch-up vaccination coverage for females (from 30,000 to 175,000), corresponding to Strategy B, will significantly alter the trajectory and burden of the disease in Korea (from persistence, as shown for Strategy A, to possible elimination). Hence, extending the current post-NIP routine vaccination in Korea to also include catch-up vaccination for females of ages 18-26 will strongly enhance HPV elimination prospect in Korea. Similarly, Table 10 shows that the additional vaccination of 80% of boys (corresponding to Strategy C) can change the trajectory and burden from persistence (for Strategy A) to possible elimination. Here, too, extending the post-NIP vaccination strategy to include the vaccination of boys will significantly enhance the aforementioned elimination prospect. These increases in vaccination coverage (for girls, adult females, and even boys) are generally realistically-attainable in Korea (considering the high coverages for infant vaccination series covered by the NIP in Korea (KCDC Center for Infectious Disease Control 2023)). Table 10 also shows that increasing the screening coverage from the 61% baseline to 90% only has marginal impact in the values of the reproduction numbers for the three strategies (this may be due to the assumption in the model that screened or detected HPV-infected individuals have reduced transmissibility, compared to undetected individuals).

Although the National Immunization Program for HPV in Korea started in 2016 (providing routine vaccination to 12-year old girls) (Kwak and Hwang 2024), its real effect did not begin to be felt in the population until the vaccinated girls reach the age of commencement of sexual activity (which is assumed to be 17 years in Korea). Consequently, the simulations to be carried out in the study will start in the year 2021 (when the vaccinated girls have started to become sexually-active).

Impact of vaccination and pap screening on cancer cases and mortality

The first set of simulations carried out are for assessing the impact of Pap screening and vaccination (of girls, females, and boys; corresponding to Strategies A, B, and C, as described above) on the yearly incidence and mortality caused by HPV-related cancers in Korea for the period between 2021 and 2130. These simulation results, obtained using the parameter values in Table 9 and the caption of Figure 8 to simulate the model (2.5), show a pronounced decline in the number of cervical cancer cases and mortality in females (for any of the three strategies) under both the baseline (61%) and improved (90%) Pap screening coverages, as depicted in Figure 8 (subplots (a)-(d)). Specifically, Figure 8 shows that while Strategy A induces the highest cervical cancer burden in females (red curve in subplots (a)-(d)), Strategies B and C incur much lower (and nearly the same) burden over the same period. These figures also show that vaccinating boys at 80% coverage (Strategy C) induces very strong reduction in cervical cancer burden in females, at a level almost the same as the reduction recorded under Strategy B, which entails allocating all the vaccine resources to girls and females only (compare the blue and gold curves in subplots (a)-(d)). These results suggest a strong spillover benefit (i.e., the vaccination of boys reduces the incidence and mortality due to cervical cancer in females). Furthermore, this figure shows that although the increase in Pap screening coverage (from 61% to 90%) has only marginal effect on the incidence of cervical cancer (compare subplots (a) and (b)), such an increase has a significant impact in reducing cervical cancer mortality (compare subplots (c) and (d)). It should be noted from these plots that while the burden of cervical cancer persists under Strategy A, such burden dies out with time under Strategies B and C (this is because the control reproduction number for Strategy A exceeds unity, while those for Strategies B and C are below unity, as tabulated in Table 10; the time to elimination of the HPV burden, for Strategies B and C, are given in Table 11).

**Fig. 8** Simulation of the model (2.5) to assess the combined impact of Pap screening and vaccination (of girls, females, and boys) on the burden of cervical cancer in females and HPV-related cancers in males, for the aforementioned Strategies A, B, and C and varying Pap screening coverages. (a) Yearly cervical cancer cases in females as a function of time from 2021 to 2130 for Pap screening coverage at 61%, (b) Yearly cervical cancer cases in females from 2021 to 2130 for Pap screening coverage at 90%, (c) Yearly cervical cancer-induced deaths in females from 2021 to 2130 for Pap screening coverage at 61%, (d) Yearly cervical cancer-induced deaths in females from 2021 to 2130 for Pap screening coverage at 90%, (e) Yearly HPV-related cancer cases in males from 2021 to 2130 for Pap screening coverage at 61%, (f) Yearly HPV-related cancer cases in males from 2021 to 2130 for Pap screening coverage at 90%. Parameter values used in these simulations are as given in Table 9, but with the following changes: (i) Strategy A: $q_{f} = 0.8$ , (ii) Strategy B: $q_{f} = 0.95$ , $ξ_{f} = 0.06$ per year, (iii) Strategy C: $q_{f} = q_{m} = 0.8$ . In all simulations for the 90% Pap screening coverage, the following values were used: $d_{u 3} g_{f} = 0.27, h_{u 3} y_{f} = q_{u 3} z_{f} = w_{f} = 0.35$ per year.

Table 10 Values of the control reproduction number, $R_{vs}$ , of the model (2.5) based on (a) Strategy A ( $q_{f} = 0.8$ ), (b) Strategy B ( $q_{f} = 0.95, ξ_{f} = 0.06$ ), and (c) Strategy C ( $q_{f} = q_{m} = 0.8$ ), under (i) baseline Pap screening coverage of 61% and (ii) increased Pap screening coverage of 90% ( $d_{u 3} g_{f} = 0.27, h_{u 3} y_{f} = q_{u 3} z_{f} = w_{f} = 0.35$ ). The rest of the parameter values used are as given Table 9.

$R_{vs}$	(a) Strategy A	(b) Strategy B	(c) Strategy C
(i) Baseline Pap screening coverage (61%)	1.6181	0.9728	0.7927
(ii) Increased Pap screening coverage (90%)	1.6134	0.97	0.7904

Table 11 Simulation of the model (2.5), showing time-to-elimination of cervical cancer in females in Korea for Strategies B and C and varying Pap screening coverages (using parameter values such that the control reproduction number, $R_{vs}$ , is less than one, as given in Table 10, so that the disease-free equilibrium of the model is locally-asymptotically stable, by Theorem 2.3). Parameter values used in these simulations are as given in Table 9, but with the following changes: (i) Strategy B: $q_{f} = 0.95$ , $ξ_{f} = 0.06$ per year, (ii) Strategy C: $q_{f} = q_{m} = 0.8$ . In all simulations for the 90% Pap screening coverage, the following values were used: $d_{u 3} g_{f} = 0.27, h_{u 3} y_{f} = q_{u 3} z_{f} = w_{f} = 0.35$ per year

Control strategy	Strategy B		Strategy C
Pap screening coverage	$61 %$ (baseline)	$90 %$	$61 %$ (baseline)	$90 %$
Time of cervical cancer elimination (year)	2088	2085	2097	2094

The simulations depicted in Figure 8 also show that while the number of HPV-related cancer cases in males initially increases during the first 10-15 years (for any of the three strategies), it sharply declines reaching elimination levels, except for Strategy A, which rebounds after about 70 years of implementation (see subplots (e) and (f) in Figure 8). This figure further shows that vaccinating girls and females has a significant spillover effect in reducing the burden of HPV-related cancers in males (see red and gold curves in subplots (e) and (f)). In these subplots, Strategy C consistently does better than Strategy B for most of the simulation duration, suggesting that if the objective is to reduce cancer cases in males, then allocating more vaccination resources to boys (Strategy C) would be more effective than allocating these resources to girls and females (Strategy B). Furthermore, this figure shows that increasing Pap screening coverage of females (from 61% to 90%) induces very marginal spillover effect on the incidence of HPV-related cancers in males (compare Figure 8 subplots (e) and (f)).

Overall, the simulations for the impact of Pap screening and vaccination strategies implemented in Korea (depicted in Figure 8) show that increasing the vaccination coverage for girls and females (Strategy B) or the vaccination of boys at 80% coverage (Strategy C) will have a significant effect in reducing (and eliminating) cervical cancer burden in Korea (with Strategy B mostly marginally better than Strategy C, particularly during the first 60 years of implementation), and that vaccinating boys induces strong spillover effect in reducing the cervical cancer burden in females (subplots (a)-(d) of Figure 8). Furthermore, increasing the vaccination coverage for girls and females (Strategy B) or the vaccination of boys at 80% coverage (Strategy C) will have a significant effect in reducing (and eliminating) HPV-related cancer burden in males in Korea (with Strategy C always marginally better than Strategy B), and that vaccinating girls and females induces a strong spillover effect in reducing the HPV-related cancer burden in males (subplots (e) and (f) of Figure 8). The results in Figure 8 are also displayed using bar charts in Figure 9, showing similar trends and profiles, but also allowing for clear comparison of the impact between baseline (61%) and improved (90%) Pap screening coverage under each strategy. It should also be mentioned that similar trends (of the aforementioned simulation results in Figure 8 or Figure 9, for the yearly incidence and mortality of HPV-related cancers in Korea) are observed with respect to the reductions in the cumulative incidence and mortality for the HPV-related cancers, when Strategies B and C are compared with the prevailing Strategy A in Korea (as depicted in Figure 10). The reductions in the cumulative incidence and mortality for HPV-related cancers are also depicted for various vaccination coverages for boys (set at 25%, 50%, and 75%) in addition to Strategy A, and the results (depicted in Figure 13 of Appendix G.1) showed lower reductions in disease burden (compare Figure 13 with Figure 10). Similarly, simulations for Strategy A with increased vaccination of girls (set at 85%, 90%, 95%) show greater reductions in cancer incidence and mortality compared to Strategy A (compare Figure 14 in Appendix G.2 with Figure 10).

**Fig. 9** Simulation of the model (2.5) to assess the combined impact of Pap screening and vaccination (of girls, females, and boys) on the burden of HPV-related cancers in females and males during the year 2080, for Strategies A, B, and C and varying Pap screening coverages. The panels depict the number of (a) cervical cancer cases in females, (b) cervical cancer-induced deaths in females, and (c) HPV-related cancer cases in males, during 2080. Parameter values used in these simulations are as given in Table 9, but with the following changes: (i) Strategy A: $q_{f} = 0.8$ , (ii) Strategy B: $q_{f} = 0.95$ , $ξ_{f} = 0.06$ per year, (iii) Strategy C: $q_{f} = q_{m} = 0.8$ . In all simulations for the 90% Pap screening coverage, the following values were used: $d_{u 3} g_{f} = 0.27, h_{u 3} y_{f} = q_{u 3} z_{f} = w_{f} = 0.35$ per year.

**Fig. 10** Simulation of the model (2.5) to assess the combined impact of Pap screening and vaccination (of girls, females, and boys) on the burden of HPV-related cancers in females and males in Korea by the year 2080, for Strategies B and C in comparison to Strategy A under baseline Pap screening coverage. The panels show reduction in the cumulative number of (a) cervical cancer cases in females, (b) cervical cancer-induced deaths in females, and (c) HPV-related cancer cases in males, by 2080. Parameter values used in these simulations are as given in Table 9, but with the following changes: (i) Strategy A: $q_{f} = 0.8$ , (ii) Strategy B: $q_{f} = 0.95$ , $ξ_{f} = 0.06$ per year, (iii) Strategy C: $q_{f} = q_{m} = 0.8$ .

**Fig. 13** Simulation of the model (2.5) to assess the combined impact of Pap screening and vaccination (of girls, females, and boys) on the burden of HPV-related cancers in females and males in Korea by the year 2080, for three vaccination coverages for boys in addition to Strategy A, in comparison to Strategy A under baseline Pap screening coverage. The panels depict reduction in the cumulative (a) cervical cancer cases in females, (b) cervical cancer-induced deaths in females, and (c) HPV-related cancer cases in males, by 2080. Parameter values used in these simulations are as given in Table 9, but with the following changes: (i) Strategy A: $q_{f} = 0.8$ , (ii) 25% boys vaccination: $q_{f} = 0.8$ , $q_{m} = 0.25$ , (iii) 50% boys vaccination: $q_{f} = 0.8$ , $q_{m} = 0.5$ , (iv) 75% boys vaccination: $q_{f} = 0.8$ , $q_{m} = 0.75$ .

**Fig. 14** Simulation of the model (2.5) to assess the combined impact of Pap screening and vaccination (of girls and females) on the burden of HPV-related cancers in females and males in Korea by the year 2080, for three vaccination coverages for girls in addition to Strategy A, in comparison to Strategy A under baseline Pap screening coverage. The panels depict reduction in the cumulative (a) cervical cancer cases in females, (b) cervical cancer-induced deaths in females, and (c) HPV-related cancer cases in males, by 2080. Parameter values used in these simulations are as given in Table 9, but with the following changes: (i) Strategy A: $q_{f} = 0.8$ , (ii) 85% girls vaccination: $q_{f} = 0.85$ , (iii) 90% girls vaccination: $q_{f} = 0.9$ , (iv) 95% girls vaccination: $q_{f} = 0.95$ .

Since the simulations in Figure 8 show convergence to the disease-free equilibrium under strategies B and C (which is also in line with the corresponding results in Table 10, showing that the reproduction number, $R_{vs}$ , for each of the two strategies is always less than one, so that, by Theorem 2.3, the disease will die out if the initial number of HPV-infected individuals in the community is low enough), it is instructive to determine the expected time-to-elimination associated with the implementation of these strategies (B and C) under the long-term scenario. For simulation purposes, we consider elimination to occur when the total number of cervical cancer cases in Korea is less than 10. This figure is much lower than the figure corresponding to the World Health Organization’s elimination standard of four cervical cancer cases per 100,000 females (World Health Organization 2020) (which translates to 1,020 cervical cancer cases in Korea). Table 11 depicts the simulation results of the model (2.5), under the long-term scenario, showing the first time the total number of females with cervical cancer in Korea fall below 10. It follows from this table that cervical cancer elimination could be achieved faster under Strategy B with improved Pap screening coverage (with such elimination achieved by 2085) than under Strategy B with baseline Pap screening coverage (this achieves elimination by 2088) or Strategy C (which achieves elimination by 2097 or 2094 for baseline and improved Pap screening coverages, respectively). In other words, this study shows that, based on the parameter values used in our simulations, the combined anti-HPV vaccination and Pap screening strategy can lead to the elimination of cervical cancer in Korea in about 60 years if the routine vaccination coverage for girls (12-17 years of age) and the catch-up vaccination coverage for females (18-26 years of age), are increased from their current levels of 80% and 30,000 per year, respectively, to 95% and 175, 000 per year, respectively.

Spillover effect of vaccination-only strategy

In this section, the model (2.5) will be simulated for the case where vaccination is the only anti-HPV intervention implemented in Korea (i.e., no Pap screening), with the aim of quantifying and evaluating the potential spillover effect of vaccinating individuals of one gender on reducing the burden of HPV and related cancers in the other gender. Specifically, the following hypothetical scenarios of the vaccination-only strategy are considered:

(i) Girls-only vaccination strategy: This entails the annual vaccination of 99% of girls in the 12-17 year age group for the period between 2021 and 2060, and no vaccination of females (of age 18-26), boys, or males.
(ii) Boys-only vaccination strategy: This involves the annual vaccination of 99% of boys in the 12-17 year age group, and no vaccination of males (of age 18-26), girls, or females.

The model is simulated, using the baseline values of the parameters in Table 9 and in the caption of Figure 11, for the period from 2021 to 2060. The simulation results obtained are then compared with those generated under the worst-case scenario (when neither vaccination nor Pap screening is implemented).

**Fig. 11** Simulation of the model (2.5) to assess the spillover benefit of vaccinating adolescents of one gender on the burden of HPV-related cancers in the other gender, in comparison to the worst-case scenario (with no vaccination or screening). The panels depict reduced proportion of cumulative (a) cervical cancer cases in females, (b) cervical cancer-induced deaths in females, and (c) HPV-related cancer cases in males, by 2060. Parameter values used in these simulations are as given in Table 9, but with the following changes: (i) Girls-only vaccination: $q_{f} = 0.99$ , $q_{m} = 0$ , (ii) Boys-only vaccination: $q_{m} = 0.99$ , $q_{f} = 0$ . In all simulations, the following values are used: $ξ_{f} = ξ_{m} = 0$ per year (no catch-up vaccination), $d_{u 3} g_{f} = h_{u 3} y_{f} = q_{u 3} z_{f} = w_{f} = 0$ per year (no routine Pap screening).

The simulation results obtained, depicted in Figure 11, show that vaccinating boys only or girls only reduces cumulative cervical cancer cases in females, in comparison to the worst-case scenario without vaccination or screening. In particular, boys-only vaccination causes 24% reduction in cervical cancer cases by the year 2060 compared to the worst-case scenario, clearly showing a significant spillover benefit in vaccinating boys in reducing cumulative number of cervical cancer cases in females in Korea. Similarly, girls-only vaccination causes about 26% reduction in cumulative number of cervical cancer cases by the year 2060 in comparison to the worst-case scenario (see blue curve of Figure 11 (a)), which, as expected, is slightly higher than the spillover reduction generated by the boys-only vaccination (red curve in Figure 11 (a)). A similar trend is observed with respect to reductions in cumulative mortality due to cervical cancer in females (Figure 11(b)), where boys-only vaccination reduces 26% of cumulative cervical cancer-induced deaths in females by 2060, in comparison to the worst-case scenario (which is quite close to the 29% cumulative reduction obtained by vaccinating girls only at 99% coverage, in comparison to the worst-case). Finally, boys-only vaccination is, as expected, more effective in reducing the cumulative number of HPV-related cancer cases in males, than the spillover related to the implementation of girls-only vaccination (Figure 11(c)). In particular, while boys-only vaccination causes an 18% reduction of cumulative HPV-related cancers in males by the year 2060, in comparison to the baseline (where there is no vaccination or screening), girls-only vaccination only causes 11% reduction of these male cancers by the same time period, in comparison to the baseline (compare blue and red curves in Figure 11(c)). It is worth noting that the spillover effect of vaccinating boys in reducing cervical cancer cases and mortality in females is larger than that of vaccinating girls in reducing HPV-related cancer cases in males (compare subplots (a) and (b) with subplot (c) in Figure 11). This may be due to the fact that HPV-infected females are more likely to develop cancer, compared to HPV-infected males (Giuliano et al. 2008a; National Cancer Center Korea Central Cancer Registry 2022; Moscicki et al. 2004, 1998; Van Doornum et al. 1994).

Discussion and Conclusion

HPV (Human Papillomavirus) accounts for virtually all cervical cancer cases and significantly contributes to anal, penile, vaginal, and head and neck cancers (de Martel et al. 2017). Although cervical cancer is one of the few human cancers that is vaccine-preventable and curable if detected early and treated, it continues to pose major health burden globally (causing an average of over 300,000 deaths annually (Bray et al. 2024; Schiffman et al. 2007)). It is the third most common female cancer in the Republic of Korea in the 15-44 year age group (causing an average of over 1,000 deaths per year (Bruni et al. 2023)). In order to curtail or mitigate the burden of cervical cancer (and other HPV-related cancers in females and males), Korea started the National Immunization Program (NIP) for HPV in 2016. This program provides routine HPV immunization to girls in the 12-17 year age group (Kwak and Hwang 2024). Korea also has the National Cancer Screening Program (NCSP), which offers free Pap screening to females over age 30 every 2 years since 2002, and to females over age 20 since 2016 (Shin et al. 2022). The purpose of the current study was to use mathematical modeling approaches, together with data analytics and numerical simulations, to assess the population-level impact of the aforementioned vaccination and Pap screening programs on the transmission dynamics of HPV and its related cancers in South Korea. A sex-structured mathematical model for the heterosexual transmission of HPV in Korea, in the presence of HPV vaccination and Pap screening, was formulated. Although only females are vaccinated in Korea under the current NIP, the model developed in this study allows for the vaccination of males (particularly boys of age group 12-17 years). This allows for the hypothetical assessment of the potential impact of male vaccination on reducing the burden of HPV and related cancers in Korea.

The disease-free equilibrium of this model was shown to be locally-asymptotically stable when a certain epidemiological threshold, known as the control reproduction number, is less than one. The epidemiological implication of this result is that the combined vaccination and Pap screening program implemented in South Korea could lead to the elimination of HPV and related cancers if the control reproduction number could be brought to (and maintained at) a value less than one, provided that the initial number of HPV-infected individuals is small enough. Numerical simulations of the model showed that vaccine-derived herd immunity cannot be achieved in Korea under the current NIP coverage (which essentially amounts to 88% of females vaccinated at steady-state). The herd immunity can be achieved if, additionally, at least 65% of males are vaccinated at steady-state (if males are not vaccinated, which is the current situation under the NIP, then at least 99% of females need to be vaccinated to achieve herd immunity). Further rigorous analyses showed that the special case of the model with negligible disease-induced mortality could have one or three endemic equilibria when the associated reproduction number exceeds one. Using a Krasnoselskii sub-linearity argument, it was shown that, for this special case and with vaccine efficacy assumed to be perfect, the unique endemic equilibrium of the model is locally-asymptotically stable whenever its associated reproduction number exceeds one. The implication of this result is that vaccination and Pap screening, at these coverage and efficacy levels, are inadequate to lead to the elimination of HPV and related cancers in Korea (since, in this case, HPV and related cancers will persist in the population).

In order to estimate some of the unknown parameters of the model, a simplified version of the model without Pap screening and with the three cervical intraepithelial neoplasia (CIN) and intraepithelial neoplasia in males (INM) stages lumped into one epidemiological compartment, was calibrated using observed data for cumulative cases and mortality for HPV-related cancers in Korea for the period from 1999 to 2020 (National Cancer Center Korea Central Cancer Registry 2022). The disease-free equilibrium of the model is shown to be locally-asymptotically stable when the reproduction number of the model is less than unity. The epidemiological implication of this result is that the interventions implemented (i.e., vaccination and Pap screening) can lead to the effective control and elimination of HPV when the reproduction number of the model is less than one provided the initial number of infected individuals is small enough (i.e., fall within the basin of attraction of the disease-free equilibrium). The local asymptotic stability was extended to the global asymptotic stability of the disease-free equilibrium for a special case when the disease-induced mortality and the coverage for catch-up vaccination (for both females and males) in Korea are negligible, when its associated reproduction number is less than one. Epidemiologically-speaking, this result show that HPV (and, hence, HPV-related cancers) will be eliminated in Korea if the associated control reproduction number can be brought down to (and maintained at) a value less than one, regardless of the initial size of the number of HPV-infected individuals in the community.

Using the estimated values of the unknown parameters obtained from this calibration, the full model was then used to conduct global parameter sensitivity analysis to determine the parameters of the model that have the highest impact on a chosen response function, namely the control reproduction number $R_{vs}$ . The results of the sensitivity analysis revealed that parameters related to the transmission of HPV, recovery from symptomatic HPV infection, and HPV vaccine coverage and efficacy were the most influential in determining the value of the control reproduction number (hence, the burden of HPV and related cancers in Korea). The implication of these results are that effective control and mitigation strategies could be formulated by targeting these identified parameters.

The full model was simulated to assess the population-level impact of the combined vaccination and Pap screening strategies currently being implemented in Korea (defined as Strategy A), and to evaluate the potential impact of two hypothetical extended versions of these strategies (namely Strategies B and C). These strategies are formally defined below:

Strategy A: This is the post-NIP vaccination strategy that entails the annual vaccination of 80% of girls (12-17 years of age) and annual catch-up vaccination of 30,000 females (18-26 years of age) with baseline Pap screening (of 61% of adult females in the age group 20-74 years of age).
Strategy B: This contains Strategy A together with the increase in the annual vaccination coverage for girls from 80% to 95% and for females from 30,000 to 175,000.
Strategy C: This contains Strategy A together with the annual vaccination of 80% of boys (12-17 years of age).

The simulation results obtained showed that while all three strategies will lead to a significant reduction in the disease burden (measured in terms of yearly cases and mortality of HPV-related cancers), only Strategies B and C will lead to the elimination of HPV and related cancers in Korea. This is because only the latter two strategies will result in bringing the associated control reproduction number of the model (denoted by

R_{vs}

) below one. In particular, it was shown that, based on the parameter values used in our simulations, that cervical cancer can be eliminated in Korea using Strategy B in about 60 years (Strategy C could achieve elimination in about 70 years). These simulations suggest that increasing the current levels of the vaccination and Pap screening coverages in Korea, as well as extending the vaccination coverage to adult females and males (to specifically attain the levels in Strategies B and C) could lead to a dramatic reduction (and even elimination) of the burden of HPV and related cancers in Korea.

The model was further simulated to assess the potential impact of vaccinating one gender (only) on reducing the burden of HPV and related cancers in the other gender for the hypothetical case where Pap screening is not implemented. In other words, the model was used to explore and quantify the potential spillover benefit of vaccination of one gender on the other (when Pap screening is not implemented). These simulation results obtained, which are based on vaccinating 99% of one gender only, showed that vaccinating 99% of boys only induced a significant spillover benefit in reducing the cumulative number of cervical cancer cases and mortality in females in Korea, in comparison to the worst-case scenario where no vaccination or Pap screening is implemented. This reduction was larger than the spillover benefit achieved (in reducing the disease burden in males) by vaccinating 99% of girls only. Thus, this study showed that extending the current post-NIP strategy in Korea to also include the vaccination of boys (and catch-up vaccination of adult females) will significantly reduce the burden of HPV and related cancers in both females and males in Korea.

Some of the limitation of the model developed in this study include not accounting for homosexual transmission of HPV, reinfection of HPV and waning of infection-acquired and vaccine-acquired immunity against HPV. Adding the dynamics of MSM (men who have sex with men) into the model might affect the level or size of some of the epidemiological thresholds presented in this study. The authors hope to address these limitations in a future study. It should, however, be mentioned that the modeling framework presented in this study, for the dynamics of HPV and related cancers in Korea, can be adapted and applied to other countries or regions or jurisdictions interested in assessing the population-level impact of vaccination and Pap screening on reducing the burden of HPV and related cancers, by adjusting demographic parameters, the vaccination and Pap screening parameters, and the parameter for number of sexual contacts pertinent to the area. Overall, this study showed that the prospect of the effective control or elimination of HPV and related cancers in Korea, using the combined vaccination and Pap screening program, is promising, provided that the coverage levels are high enough and the program is extended to include the vaccination of boys and adult females (with high enough coverage).

Appendices

Proof of Lemma 2.1

Proof

Consider the model (2.5) with

Q_{2 u} (0)

,

Q_{2 d} (0)

,

J_{2 m} (0) > 0

, and all other state variables being non-negative at time

t = 0

. The equation for the rate of change of the compartment for unvaccinated susceptible females (

S_{f} (t)

) can be rewritten as:

\begin{matrix} \frac{d}{dt} [S_{f} (t) exp ((ξ_{f} + μ) t + \int_{0}^{t} λ_{m_{a} f} (s) d s)] \\ = (1 - q_{f}) π_{f} exp ((ξ_{f} + μ) t + \int_{0}^{t} λ_{m_{a} f} (s) d s), \end{matrix}

and integrating both sides with respect to t gives,

\begin{matrix} S_{f} (t) exp ((ξ_{f} + μ) t + \int_{0}^{t} λ_{m_{a} f} (s) d s) \\ = S_{f} (0) + (1 - q_{f}) π_{f} \int_{0}^{t} exp ((ξ_{f} + μ) t + \int_{0}^{t} λ_{m_{a} f} (s) d s) d t, \end{matrix}

from which it follows that (noting that

0 \leq q_{f} \leq 1

),

\begin{matrix} S_{f} (t) & = exp (- (ξ_{f} + μ) t - \int_{0}^{t} λ_{m_{a} f} (s) d s) \\ [S_{f} (0) + (1 - q_{f}) π_{f} \int_{0}^{t} (exp (\int_{0}^{t} (ξ_{f} + μ) t + λ_{m_{a} f} (s) d s) d t] \\ \geq 0, for all t > 0 . \end{matrix}

Next, given the non-negativity of

S_{f} (t)

, the equation for the rate of change of the vaccinated female population (

V_{f} (t)

) satisfies the following:

\begin{matrix} \frac{d V_{f}}{dt} = q_{f} π_{f} + ξ_{f} S_{f} - (1 - ε_{v}) λ_{m_{a} f} V_{f} - μ V_{f} \geq q_{f} π_{f} - (1 - ε_{v}) λ_{m_{a} f} V_{f} - μ V_{f}, \end{matrix}

which can be rewritten as,

\begin{matrix} \frac{d}{dt} [V_{f} (t) exp (μ t + \int_{0}^{t} (1 - ε_{v}) λ_{m_{a} f} (s) d s)] \\ \geq q_{f} π_{f} exp (μ t + \int_{0}^{t} (1 - ε_{v}) λ_{m_{a} f} (s) d s), \end{matrix}

from which it follows (by integrating both sides) that,

\begin{matrix} V_{f} (t) exp (μ t + \int_{0}^{t} (1 - ε_{v}) λ_{m_{a} f} (s) d s) \\ \geq V_{f} (0) + q_{f} π_{f} \int_{0}^{t} exp (μ t + \int_{0}^{t} (1 - ε_{v}) λ_{m_{a} f} (s) d s) d t . \end{matrix}

Hence,

V_{f} (t) \geq 0

for all

t > 0

. Similarly, using this approach, it can be shown that the state variables

E_{f} (t), I_{f} (t), P_{f} (t), S_{m} (t), V_{m} (t), E_{m} (t), I_{m} (t)

and

P_{m} (t)

of the model (2.5) are also non-negative for all

t > 0

. The next step is to prove the non-negativity of the state variables

Q_{iu}, Q_{id}, J_{im} (i = 1, 2, 3)

of the model (2.5). This is done by contradiction as follows. First of all, the equations for the rate of change of

Q_{iu} (i = 1, 2, 3)

can be rewritten as:

\begin{matrix} \frac{d}{dt} [Q_{1 u} exp ((g_{f} + μ) t)] & = [(1 - k_{f}) α_{f} P_{f} + h_{u 2} y_{f} Q_{2 u}] exp ((g_{f} + μ) t), \\ \frac{d}{dt} [Q_{2 u} exp ((y_{f} + μ) t)] & = [(1 - d_{u 1} - d_{u 3}) g_{f} Q_{1 u} + q_{u 2} z_{f} Q_{3 u}] exp ((y_{f} + μ) t), \\ \frac{d}{dt} [Q_{3 u} exp ((z_{f} + μ) t)] & = [(1 - h_{u 1} - h_{u 2} - h_{u 3}) y_{f} Q_{2 u}] exp ((z_{f} + μ) t), \end{matrix}

A.1

which, by integrating both sides with respect to t, gives

\begin{matrix} Q_{1 u} (t) & = Q_{1 u} (0) + [\int_{0}^{t}, {(1 - k_{f}) α_{f} P_{f} (s) + h_{u 2} y_{f} Q_{2 u} (s)}) \\ (exp ((g_{f} + μ) s) d s] exp ((- g_{f} - μ) t), \\ Q_{2 u} (t) & = Q_{2 u} (0) + [\int_{0}^{t}, {, (1 - d_{u 1} - d_{u 3}), g_{f}, Q_{1 u}, (s)) \\ (+ q_{u 2} z_{f} Q_{3 u} (s)} exp ((y_{f} + μ) s) d s] exp ((- y_{f} - μ) t), \\ Q_{3 u} (t) & = Q_{3 u} (0) + [(1 - h_{u 1} - h_{u 2} - h_{u 3}), y_{f}, \int_{0}^{t}, Q_{2 u}, (s)) \\ (exp ((z_{f} + μ) s) d s] exp ((- z_{f} - μ) t) . \end{matrix}

A.2

Suppose, now, that

Q_{2 u} (t) \leq 0

for some t in

(0, \infty) .

Since

Q_{2 u} (0) > 0

(by assumption), there exists a

t_{0} > 0

such that

t_{0} = inf {t > 0 : Q_{2 u} (t) = 0} .

Note that

P_{f} (t) \geq 0

for all

t \geq 0

from the proof above. Furthermore, all the model parameters, as well as the proportions

1 - k_{f}, 1 - d_{u 1} - d_{u 3}, 1 - h_{u 1} - h_{u 2} - h_{u 3}

, are non-negative (as defined in Section 2). Hence, it follows from (A.2) and the definition of

t_{0},

that both

Q_{1 u} (t)

and

Q_{3 u} (t)

are non-negative for any

t \in [0, t_{0}] .

Furthermore, it can be seen by evaluating the second equation of (A.2) at

t = t_{0}

, that

\begin{matrix} Q_{2 u} (t_{0}) = & Q_{2 u} (0) + [\int_{0}^{0} {(1 - d_{u 1} - d_{u 3}) g_{f} Q_{1 u} (s) + q_{u 2} z_{f} Q_{3 u} (s)} exp ((y_{f} + μ) s) d s] \\ exp ((- y_{f} - μ) t_{0}) . \end{matrix}

Using the non-negativity of

Q_{1 u} (t)

and

Q_{3 u} (t)

on

[0, t_{0}],

and the strict positivity of

Q_{2 u} (0)

, the above expression for

Q_{2 u} (t_{0})

is strictly positive. This contradicts the definition of

t_{0}

(which requires

Q_{2 u} (t_{0}) = 0

). Thus,

Q_{2 u} (t) > 0

for all

t \geq 0 .

Using the fact that

Q_{2 u} (t)

is strictly positive for all

t \geq 0

, it can be seen from the first and third equations of (A.2), that

Q_{1 u} (t)

and

Q_{3 u} (t)

are also strictly positive for all

t > 0 .

Similarly, the non-negativity of

Q_{id} (t)

, with

i = 1, 2, 3

, can be proved, using the fact that

Q_{1 u} (t), Q_{2 u} (t)

and

Q_{3 u} (t)

are all strictly positive for

t > 0

. Furthermore, the proof by contradiction can also be used to show that

J_{im} (t)

, with

i = 1, 2, 3

, are strictly positive for all

t > 0

. Finally, the non-negativity of the remaining state variables,

C_{u} (t), C_{d} (t), R_{f}^{C} (t), R_{f} (t), C_{m} (t), R_{m}^{C} (t)

, and

R_{m} (t)

, follows from the non-negativity of the state variables established above. Hence, all non-negative initial conditions of the model (2.5), with

Q_{2 u} (0), Q_{2 d} (0), J_{2 m} (0) > 0

, remain non-negative for all

t > 0

.

□

Non-negativity of the Control Reproduction Number ( $R_{vs}$ )

Recall from Section 2.2.1 that the control reproduction number of the model (2.5) is given by

\begin{matrix} R_{vs} = ρ (F_{s} V_{s}^{- 1}) = \sqrt{R_{f} R_{m}}, \end{matrix}

where (with all the associated parameters and quantities as defined in Section 2.2.1),

\begin{matrix} R_{f} = & \frac{β_{fm} c_{mf} μ p_{m}^{*}}{π_{m}} [\frac{η_{f}}{k_{1}} + \frac{σ_{f}}{k_{1} k_{2}} + \frac{θ_{f} σ_{f} l_{1}}{k_{1} k_{2} k_{3}} +) \\ + (\frac{θ_{f} σ_{f} l_{1} l_{2} {ι_{f 2} + κ_{f} l_{3} k_{6} + κ_{f} l_{3} l_{4}}}{k_{1} k_{2} k_{3} ι_{f 1}} + \frac{θ_{f} ν σ_{f} l_{1} l_{2} (j_{f 1} + κ_{fd} j_{f 2} + κ_{fd} j_{f 3})}{k_{1} k_{2} k_{3} ι_{f 1} ι_{f 3}}], \\ R_{m} = & \frac{β_{mf} c_{mf} μ p_{f}^{*}}{π_{m}} [\frac{η_{m}}{e_{1}} + \frac{σ_{m}}{e_{1} e_{2}} + \frac{θ_{m} σ_{m} ψ_{m} (1 - b_{m})}{e_{1} e_{2} e_{3}}) \\ (+ \frac{θ_{m} σ_{m} ψ_{m} (1 - b_{m}) α_{m} (1 - k_{m}) {ι_{m 2} + κ_{m} g_{m} (1 - d_{m}) e_{6} + κ_{m} g_{m} (1 - d_{m}) f_{1}}}{e_{1} e_{2} e_{3} ι_{m 1}}] . \end{matrix}

Since the parameters of the model are all non-negative and the proportions related to some of the parameters are between zero and one, the non-negativity of

R_{f}

and

R_{m}

(hence that of

R_{vs}

) depend on the non-negativity of the expressions for

ι_{f 1}, ι_{f 2}, ι_{f 3}, j_{f 1}, j_{f 2}, j_{f 3}, ι_{m 1}

, and

ι_{m 2}

. In order to prove the non-negativity of

j_{f 1}, j_{f 2},

and

j_{f 3}

, it is convenient to rewrite their expressions by factoring out

g_{f}

(which is automatically non-negative) and organizing the inner expressions in terms of

d_{u 3}

(for

j_{f 1}

) and

μ

(for

j_{f 2}

and

j_{f 3}

) in descending order (note that

d_{u 3}

and

μ

are non-negative), as follows:

\begin{matrix} j_{f 1} & = g_{f} (ι_{f 4} d_{u 3} + ι_{f 5}), \\ j_{f 2} & = g_{f} (ι_{f 6} μ^{3} + ι_{f 7} μ^{2} + ι_{f 8} μ + ι_{f 9}), \\ j_{f 3} & = g_{f} (ι_{f 10} μ^{2} + ι_{f 11} μ + ι_{f 12}), \end{matrix}

where,

\begin{matrix} ι_{f 4} = & {y_{f} z_{f} (1 - q_{u_{2}} (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}})) + μ (y_{f} + z_{f}) + μ^{2}} \\ {ζ_{2} ζ_{3} (1 - q_{d_{2}} (1 - h_{d_{1}} - h_{d_{2}})) \\ + μ (ζ_{2} + ζ_{3}) + μ^{2}}, \\ ι_{f 5} = & y_{f} ζ_{2} h_{d_{2}} h_{u_{3}} (1 - d_{u_{1}} - d_{u_{3}}) (z_{f} + μ) (ζ_{3} + μ) \\ + y_{f} z_{f} ζ_{2} ζ_{3} q_{u_{3}} q_{d_{2}} h_{d_{2}} (1 - d_{u_{1}} - d_{u_{3}}) (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}}), \\ ι_{f 6} = & d_{u_{3}} ζ_{1} (1 - d_{d_{1}}) + h_{u_{3}} y_{f} (1 - d_{u_{1}} - d_{u_{3}}), \\ ι_{f 7} = & d_{u_{3}} ζ_{1} (1 - d_{d_{1}}) (y_{f} + ζ_{3} + z_{f}) \\ + (1 - d_{u_{1}} - d_{u_{3}}) h_{u_{3}} y_{f} (ζ_{1} + ζ_{3} + z_{f}), \\ ι_{f 8} = & d_{u_{3}} ζ_{1} {y_{f} ζ_{3} + y_{f} z_{f} (1 - (1 - h_{u_{1}} - h_{u_{2}} \\ - h_{u_{3}}) q_{u_{2}}) + ζ_{3} z_{f}} (1 - d_{d_{1}}) + (1 - d_{u_{1}} - d_{u_{3}}) \\ h_{u_{3}} y_{f} (ζ_{1} ζ_{3} + ζ_{1} z_{f} + ζ_{3} z_{f}) + (1 - d_{u_{1}} \\ - d_{u_{3}}) (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}}) q_{d_{2}} q_{u_{3}} y_{f} ζ_{3} z_{f}, \\ ι_{f 9} = & d_{u_{3}} y_{f} ζ_{1} ζ_{3} z_{f} (1 - d_{d_{1}}) (1 - q_{u_{2}}) + h_{u_{3}} y_{f} ζ_{1} ζ_{3} z_{f} (1 - d_{u_{1}} - d_{u_{3}}) \\ + (1 - d_{u_{1}} - d_{u_{3}}) \\ (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}}) q_{d_{2}} q_{u_{3}} y_{f} ζ_{1} ζ_{3} z_{f} + (1 - d_{d_{1}}) \\ (h_{u_{1}} + h_{u_{2}} + h_{u_{3}}) d_{u_{3}} q_{u_{2}} y_{f} ζ_{1} ζ_{3} z_{f}, \\ ι_{f 10} = & d_{u_{3}} ζ_{1} ζ_{2} (1 - h_{d_{1}} - h_{d_{2}}) (1 - d_{d_{1}}) \\ + h_{u_{3}} y_{f} ζ_{2} (1 - d_{u_{1}} - d_{u_{3}}) (1 - h_{d_{1}} - h_{d_{2}}) \\ + y_{f} z_{f} q_{u_{3}} (1 - d_{u_{1}} - d_{u_{3}}) (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}}), \\ ι_{f 11} = & d_{u_{3}} ζ_{1} ζ_{2} (1 - h_{d_{1}} - h_{d_{2}}) (1 - d_{d_{1}}) (y_{f} + z_{f}) + h_{u_{3}} y_{f} ζ_{2} (ζ_{1} + z_{f}) \\ (1 - d_{u_{1}} - d_{u_{3}}) (1 - h_{d_{1}} - h_{d_{2}}) \\ + q_{u_{3}} y_{f} z_{f} (ζ_{1} + ζ_{2}) (1 - d_{u_{1}} - d_{u_{3}}) (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}}), \\ ι_{f 12} = & d_{u_{3}} ζ_{1} ζ_{2} y_{f} z_{f} (1 - d_{d_{1}}) (1 - h_{d_{1}} - h_{d_{2}}) \\ (1 - q_{u_{2}} (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}})) \\ + h_{u_{3}} y_{f} ζ_{1} ζ_{2} z_{f} (1 - d_{u_{1}} - d_{u_{3}}) (1 - h_{d_{1}} - h_{d_{2}}) \\ + q_{u_{3}} y_{f} z_{f} ζ_{1} ζ_{2} (1 - d_{u_{1}} - d_{u_{3}}) (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}}) (1 - h_{d_{2}} (1 - d_{d_{1}})) . \end{matrix}

B.1

Then, it immediately follows from the non-negativity of the parameters of the model (and the proportions in equations (B.1) being less than one) that

ι_{f 4}, ι_{f 5}, \dots, ι_{f 12}

are all non-negative, hence

j_{f 1}, j_{f 2},

and

j_{f 3}

are non-negative. The non-negativity of the remaining terms,

ι_{f 1}, ι_{f 2}, ι_{f 3}, ι_{m 1}

, and

ι_{m 2}

, can be shown by rewriting their expressions as follows:

\begin{matrix} ι_{f 1} = & μ^{3} + μ^{2} (g_{f} + y_{f} + z_{f}) + μ {g_{f} y_{f} (1 - h_{u 2} (1 - d_{u 1} - d_{u 3})) \\ + y_{f} z_{f} (1 - q_{u 2} (1 - h_{u 1} - h_{u 2} - h_{u 3})) \\ + g_{f} z_{f}} + g_{f} y_{f} z_{f} {(1 - q_{u 2}) (1 - h_{u 2}) \\ + h_{u 2} (d_{u 1} + d_{u 3}) + q_{u 2} (h_{u 1} + h_{u 3})}, \\ ι_{f 2} = & (1 - q_{u 2} (1 - h_{u_{1}} - h_{u_{2}} - h_{u_{3}})) y_{f} z_{f} + μ (y_{f} + z_{f} + μ), \\ ι_{f 3} = & μ^{3} + μ^{2} (ζ_{1} + ζ_{2} + ζ_{3}) + μ {ζ_{1} ζ_{2} (1 - h_{d_{2}} + d_{d_{1}} h_{d_{2}}) \\ + ζ_{2} ζ_{3} (1 - q_{d_{2}} (1 - h_{d_{1}} - h_{d_{2}})) + ζ_{1} ζ_{3}} \\ + ζ_{1} ζ_{2} ζ_{3} {(1 - h_{d_{2}}) (1 - q_{d_{2}}) + q_{d_{2}} h_{d_{1}} + d_{d_{1}} h_{d_{2}}}, \\ ι_{m 1} = & μ^{3} + μ^{2} (g_{m} + y_{m} + z_{m}) + μ {g_{m} y_{m} (1 - h_{m 2} (1 - d_{m})) \\ + y_{m} z_{m} (1 - q_{m 2} (1 - h_{m 1} - h_{m 2})) \\ + g_{m} z_{m}} + g_{m} y_{m} z_{m} {(1 - q_{m 2}) (1 - h_{m 2}) \\ + h_{m 2} d_{m} + q_{m 2} h_{m 1}}, \\ ι_{m 2} = & μ^{2} + μ (y_{m} + z_{m}) + y_{m} z_{m} (1 - q_{m 2} (1 - h_{a_{1}} - h_{a_{2}})) . \end{matrix}

B.2

Thus, since all the parameters of the model are non-negative and the proportions in equations (B.2) are less than or equal to one, it follows from (B.2) that

ι_{f 1}, ι_{f 2}, ι_{f 3}, ι_{m 1},

and

ι_{m 2}

are all non-negative. Hence, the control reproduction number,

R_{vs}

, of the model (2.5) is non-negative.

Coefficients of the Polynomial (2.20)

The expressions for the coefficients of the quartic (2.20), given in in Section 2.2.2, are given by:

\begin{matrix} a_{4} = & (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) π_{f} {(1 - ε_{v})}^{2} {R_{f} π_{f} (ξ_{m} + μ) μ + π_{m} (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))} \\ {(1 - ε_{v}) R_{f} π_{f} (ξ_{m} + μ) μ + π_{m} (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) μ} \geq 0, \\ a_{3} = & (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) π_{f} (1 - ε_{v}) (ξ_{m} + μ) [2 (1 - ε_{v}) \\ R_{f}^{2} π_{f}^{2} (ξ_{m} + μ) μ^{2} (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) \\ + R_{f} π_{f} π_{m} μ (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{m} + μ) \\ + μ) + 2 π_{m}^{2} μ {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} ((1 - ε_{v}) (ξ_{f} + μ) + μ) + R_{f} π_{f} π_{m} μ \\ (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{f} + μ) + μ) ((1 - ε_{v}) (ξ_{m} + μ) + μ)] - R_{f} R_{m} π_{f} π_{m} \\ {(1 - ε_{v})}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} {R_{f} π_{f} (1 - ε_{v}) (ξ_{m} + μ) μ \\ + π_{m} {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2}} \\ \geq & (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) π_{f} (1 - ε_{v}) (ξ_{m} + μ) [π_{m}^{2} μ (μ (1 - ε_{v} q_{m}) \\ + ξ_{m} (1 - ε_{v}))^{2} ((1 - ε_{v}) (ξ_{f} + μ) + μ) \\ + R_{f} π_{f} π_{m} μ (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{f} + μ) + μ) ((1 - ε_{v}) (ξ_{m} + μ) + μ)] \\ - R_{vs} π_{f} π_{m} {(1 - ε_{v})}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} {R_{f} π_{f} (1 - ε_{v}) (ξ_{m} + μ) μ \\ + π_{m} {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2}} \\ = & {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} π_{f} π_{m}^{2} (1 - ε_{v}) (ξ_{m} + μ) μ {(μ (1 - ε_{v} q_{f}) \\ + ξ_{f} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{f} + μ) + μ) \\ - R_{vs} (1 - ε_{v}) (ξ_{f} + μ) μ} + R_{f} π_{f}^{2} π_{m} (1 - ε_{v}) (ξ_{m} + μ) μ {(μ (1 - ε_{v} q_{f}) \\ + ξ_{f} (1 - ε_{v})) (μ (1 - ε_{v} q_{m}) \\ + ξ_{m} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{f} + μ) + μ) ((1 - ε_{v}) (ξ_{m} + μ) + μ) \\ - R_{vs} {(1 - ε_{v})}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2}} \\ \geq & (1 - R_{vs}) {{(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} π_{f} π_{m}^{2} (1 - ε_{v}) (ξ_{m} + μ) μ {(1 - ε_{v})}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} \\ + R_{f} π_{f}^{2} π_{m} (1 - ε_{v}) (ξ_{m} + μ) μ {(1 - ε_{v})}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2}} \geq 0 (if R_{vs} \leq 1), \\ a_{2} = & (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) π_{f} (ξ_{m} + μ) [(1 - ε_{v}) R_{f}^{2} π_{f}^{2} (ξ_{m} + μ) μ^{2} {(μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v}))}^{2} \\ + R_{f} π_{f} π_{m} μ (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{m} + μ) + μ) \\ ((1 - ε_{v}) (ξ_{f} + μ) + μ) + π_{m}^{2} μ {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} {((1 - ε_{v}) (ξ_{f} + μ) + μ)}^{2} \\ + π_{m} (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) (ξ_{f} + μ) μ^{2} {R_{f} π_{f} (1 - ε_{v}) \\ ((1 - ε_{v}) (ξ_{m} + μ) + μ) + 2 π_{m} (1 - ε_{v}) \\ (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}] - R_{f} R_{m} π_{f} π_{m} (1 - ε_{v}) (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} \\ {2 R_{f} π_{f} (1 - ε_{v}) (ξ_{m} + μ) μ \\ (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) + π_{m} (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} \\ + π_{m} {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} ((1 - ε_{v}) (ξ_{f} + μ) + μ)} \\ \geq & R_{f} (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) π_{f}^{2} π_{m} μ {(μ (1 - ε_{v} q_{m}) \\ + ξ_{m} (1 - ε_{v})) (ξ_{f} + μ) (1 - ε_{v}) μ ((1 - ε_{v}) \\ (ξ_{m} + μ) + μ) + (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{f} + μ) + μ) \\ ((1 - ε_{v}) (ξ_{m} + μ) + μ) - 2 R_{vs} {(1 - ε_{v})}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2}} + (μ (1 - ε_{v} q_{m}) \\ + ξ_{m} (1 - ε_{v}))^{2} π_{f} π_{m}^{2} \\ (ξ_{m} + μ) {(μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) μ {((1 - ε_{v}) (ξ_{f} + μ) + μ)}^{2} - R_{vs} (1 - ε_{v}) (ξ_{f} + μ) \\ ((1 - ε_{v}) (ξ_{f} + μ) + μ) μ^{2} + (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) (1 - ε_{v}) (ξ_{f} + μ) μ^{2} (2 - R_{vs})} \\ \geq & (1 - R_{vs}) [R_{f} (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) π_{f}^{2} π_{m} μ {(1 - ε_{v})}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} \\ + (μ (1 - ε_{v} q_{m}) \\ + ξ_{m} (1 - ε_{v}))^{2} π_{f} π_{m}^{2} (ξ_{m} + μ) {(1 - ε_{v}) (ξ_{f} + μ) ((1 - ε_{v}) (ξ_{f} + μ) + μ) μ^{2} + (μ (1 - ε_{v} q_{f}) \\ + ξ_{f} (1 - ε_{v})) (1 - ε_{v}) (ξ_{f} + μ) μ^{2}}] \geq 0 (if R_{vs} \leq 1), \\ a_{1} = & (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) π_{f} π_{m} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} [R_{f} π_{f} \\ (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{m} + μ) + μ) + 2 π_{m} (μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) ((1 - ε_{v}) \\ (ξ_{f} + μ) + μ)] - R_{f} R_{m} π_{f} π_{m} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} \\ {π_{m} (1 - ε_{v}) (ξ_{f} + μ) μ {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} \\ + R_{f} π_{f} (1 - ε_{v}) (ξ_{m} + μ) μ {(μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v}))}^{2} + π_{m} (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) \\ {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} ((1 - ε_{v}) (ξ_{f} + μ) + μ)} \\ \geq & {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} π_{f} π_{m}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} {2 (μ (1 - ε_{v} q_{f}) \\ + ξ_{f} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{f} + μ) \\ + μ) - R_{vs} (1 - ε_{v}) (ξ_{f} + μ) μ - R_{vs} (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{f} + μ) + μ))} \\ + R_{f} {(μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v}))}^{2} π_{f}^{2} π_{m} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} {(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v})) \\ ((1 - ε_{v}) (ξ_{m} + μ) + μ) - R_{vs} (1 - ε_{v}) (ξ_{m} + μ) μ} \\ \geq & (1 - R_{vs}) [{(μ (1 - ε_{v} q_{m}) + ξ_{m} (1 - ε_{v}))}^{2} π_{f} π_{m}^{2} (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} {(1 - ε_{v}) (ξ_{f} + μ) μ \\ + (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) ((1 - ε_{v}) (ξ_{f} + μ) + μ)} + R_{f} (μ (1 - ε_{v} q_{f}) \\ + ξ_{f} (1 - ε_{v}))^{2} π_{f}^{2} π_{m} \\ (ξ_{f} + μ) (ξ_{m} + μ) μ^{2} (1 - ε_{v}) (ξ_{m} + μ) μ] \geq 0 (if R_{vs} \leq 1), \\ a_{0} = & (μ (1 - ε_{v} q_{f}) + ξ_{f} (1 - ε_{v})) (μ (1 - ε_{v} q_{m}) \\ + ξ_{m} (1 - ε_{v}))^{2} π_{f} π_{m}^{2} {(ξ_{f} + μ)}^{2} (ξ_{m} + μ) μ^{3} [1 - R_{f} R_{m}] \\ \geq & 0 (if and only if R_{vs} \leq 1) . \end{matrix}

It should be mentioned that the inequalities above use the fact that (for

i = f, m

):

μ (1 - ε_{v} q_{i}) + ξ_{i} (1 - ε_{v}) \geq (μ + ξ_{i}) (1 - ε_{v}),

since

0 \leq q_{i} \leq 1

.

Proof of Theorem 2.6

Proof

Consider the special case of the model (2.5) with negligible disease-induced mortality (

δ = δ_{d} = 0

) and a perfect vaccine (

ε_{v} = 1

). Then, the coefficients of the polynomial (2.20) satisfy

a_{4} = a_{3} = 0

and

a_{2} \geq 0

, hence the quartic equation becomes quadratic in

λ_{mf}^{* *}

. In this case, the special case of the model (with

δ = δ_{d} = 0

and

ε_{v} = 1

) has a unique endemic equilibrium if

a_{0} < 0

(i.e.,

{\tilde{R}}_{vs} > 1

), by the Descartes’ Rule of Signs. It is convenient to define this unique endemic equilibrium by

{\tilde{E}}_{1 s} = E_{1 s} |_{ε_{v} = 1}

. The local asymptotic stability of the associated unique endemic equilibrium (EE) of the special case of the model can be established, when

{\tilde{R}}_{vs} > 1

, using the Krasnoselskii sub-linearity approach (as in (Esteva and Vargas 2000; Hethcote and Thieme 1985; Pant and Gumel 2024; Thieme 1985)). The method is based on the fact that, for a given

n -

dimensional system of nonlinear differential equations of the form

\dot{x} = f (x)

, its equilibrium

E

is locally-asymptotically stable if its corresponding linearized system,

\dot{Z} = D f (E) Z

does not have a solution of the form

Z (t) = Z e^{ω t}

, where

Z \in C^{n} \ {0}

,

ω \in C

, and Re

(ω) \geq 0

. Since

δ = δ_{d} = 0

, it follows that the total female (

N_{f} (t)

) and male (

N_{m} (t)

) populations at time t are asymptotically constant. That is,

(N_{f} (t), N_{m} (t)) \to (N_{f}^{*}, N_{m}^{*})

, as

t \to \infty

. Substituting

N_{m}^{*}

and

N_{f}^{*}

for

N_{m} (t)

and

N_{f} (t)

in the model shows that,

\begin{matrix} S_{f} = N_{f}^{*} - V_{f} - E_{f} - I_{f} - P_{f} - Q_{1 u} - Q_{2 u} - Q_{3 u} - Q_{1 d} - Q_{2 d} \\ - Q_{3 d} - C_{u} - C_{d} - R_{f}^{C} - R_{f}, \\ S_{m} = N_{m}^{*} - V_{m} - E_{m} - I_{m} - P_{m} - J_{1 m} - J_{2 m} - J_{3 m} - C_{m} - R_{m}^{C} - R_{m} . \end{matrix}

D.1

Using (D.1), the special case of the model (2.5), with

δ = δ_{d} = 0

and

ε_{v} = 1

, can be rewritten as:

\begin{matrix} \begin{matrix} Females \{\begin{matrix} {\dot{V}}_{f} & = & q_{f} π_{f} + ξ_{f} (N_{f}^{*} - V_{f} - E_{f} - I_{f} - P_{f} - Q_{1 u} - Q_{2 u} - Q_{3 u} - Q_{1 d} \\ - Q_{2 d} - Q_{3 d} - C_{u} - C_{d} - R_{f}^{C} - R_{f}) - μ V_{f}, \\ {\dot{E}}_{f} & = & λ_{mf} (N_{f}^{*} - V_{f} - E_{f} - I_{f} - P_{f} - Q_{1 u} - Q_{2 u} - Q_{3 u} - Q_{1 d} \\ - Q_{2 d} - Q_{3 d} - C_{u} - C_{d} - R_{f}^{C} - R_{f}) - k_{1} E_{f}, \\ {\dot{I}}_{f} & = & σ_{f} E_{f} - k_{2} I_{f}, \\ {\dot{P}}_{f} & = & l_{1} I_{f} - k_{3} P_{f}, \\ {\dot{Q}}_{1 u} & = & l_{2} P_{f} + h_{u 2} y_{f} Q_{2 u} - k_{4} Q_{1 u}, \\ {\dot{Q}}_{2 u} & = & l_{3} Q_{1 u} + q_{u 2} z_{f} Q_{3 u} - k_{5} Q_{2 u}, \\ {\dot{Q}}_{3 u} & = & l_{4} Q_{2 u} - k_{6} Q_{3 u}, \\ {\dot{Q}}_{1 d} & = & d_{u 3} g_{f} Q_{1 u} + h_{d 2} ζ_{2} Q_{2 d} - k_{7} Q_{1 d}, \\ {\dot{Q}}_{2 d} & = & h_{u 3} y_{f} Q_{2 u} + l_{6} Q_{1 d} + q_{d 2} ζ_{3} Q_{3 d} - k_{8} Q_{2 d}, \\ {\dot{Q}}_{3 d} & = & q_{u 3} z_{f} Q_{3 u} + l_{7} Q_{2 d} - k_{9} Q_{3 d}, \\ {\dot{C}}_{u} & = & l_{5} Q_{3 u} - k_{10} C_{u}, \\ {\dot{C}}_{d} & = & l_{8} Q_{3 d} + w_{f} C_{u} - k_{11} C_{d}, \\ {\dot{R}}_{f}^{C} & = & r_{u} C_{u} + r_{d} C_{d} - μ R_{f}^{C}, \\ {\dot{R}}_{f} & = & b_{f} ψ_{f} I_{f} + k_{f} α_{f} P_{f} + d_{u 1} g_{f} Q_{1 u} + h_{u 1} y_{f} Q_{2 u} + q_{u 1} z_{f} Q_{3 u} + d_{d 1} ζ_{1} Q_{1 d} \\ + h_{d 1} ζ_{2} Q_{2 d} + q_{d 1} ζ_{3} Q_{3 d} - μ R_{f}, \end{matrix}) \\ Males \{\begin{matrix} {\dot{V}}_{m} & = & q_{m} π_{m} + ξ_{m} (N_{m}^{*} - V_{m} - E_{m} - I_{m} - P_{m} - J_{1 m} - J_{2 m} - J_{3 m} \\ - C_{m} - R_{m}^{C} - R_{m}) - μ V_{m}, \\ {\dot{E}}_{m} & = & λ_{fm} (N_{m}^{*} - V_{m} - E_{m} - I_{m} - P_{m} - J_{1 m} - J_{2 m} - J_{3 m} \\ - C_{m} - R_{m}^{C} - R_{m}) - e_{1} E_{m}, \\ {\dot{I}}_{m} & = & σ_{m} E_{m} - e_{2} I_{m}, \\ {\dot{P}}_{m} & = & (1 - b_{m}) ψ_{m} I_{m} - e_{3} P_{m}, \\ {\dot{J}}_{1 m} & = & (1 - k_{m}) α_{m} P_{m} + h_{m 2} y_{m} J_{2 m} - e_{4} J_{1 m}, \\ {\dot{J}}_{2 m} & = & (1 - d_{m}) g_{m} J_{1 m} + q_{m 2} z_{m} J_{3 m} - e_{5} J_{2 m}, \\ {\dot{J}}_{3 m} & = & f_{1} J_{2 m} - e_{6} J_{3 m}, \\ {\dot{C}}_{m} & = & f_{2} J_{3 m} - e_{7} C_{m}, \\ {\dot{R}}_{m}^{C} & = & r_{m} C_{m} - μ R_{m}^{C}, \\ {\dot{R}}_{m} & = & b_{m} ψ_{m} I_{m} + k_{m} α_{m} P_{m} + d_{m} g_{m} J_{1 m} + h_{m 1} y_{m} J_{2 m} + q_{m 1} z_{m} J_{3 m} - μ R_{m}, \end{matrix}) \end{matrix} \end{matrix}

D.2

where the associated forces of infection,

λ_{mf}

and

λ_{fm}

, obtained by substituting

N_{m} = N_{m}^{*}

in the special case of the model, are given by:

\begin{matrix} λ_{mf} = (β_{mf}, c_{mf}) [\frac{I_{m} + η_{m} E_{m} + θ_{m} [P_{m} + J_{1 m} + κ_{m} (J_{2 m} + J_{3 m})]}{N_{m}^{*}}], \\ λ_{fm} = (β_{fm}, c_{mf}) \\ [\frac{I_{f} + η_{f} E_{f} + θ_{f} [P_{f} + Q_{1 u} + κ_{f} (Q_{2 u} + Q_{3 u}) + ν {Q_{1 d} + κ_{fd} (Q_{2 d} + Q_{3 d})}]}{N_{m}^{*}}], \end{matrix}

and

k_{1} = σ_{f} + μ, k_{2} = ψ_{f} + μ, k_{3} = α_{f} + μ, k_{4} = g_{f} + μ, k_{5} = y_{f} + μ, k_{6} = z_{f} + μ, k_{7} = ζ_{1} + μ, k_{8} = ζ_{2} + μ, k_{9} = ζ_{3} + μ, k_{10} = w_{f} + r_{u} + δ + μ, k_{11} = r_{d} + δ_{d} + μ, l_{1} = (1 - b_{f}) ψ_{f}, l_{2} = (1 - k_{f}) α_{f}, l_{3} = (1 - d_{u 1} - d_{u 3}) g_{f}, l_{4} = (1 - h_{u 1} - h_{u 2} - h_{u 3}) y_{f}, l_{5} = (1 - q_{u 1} - q_{u 2} - q_{u 3}) z_{f}, l_{6} = (1 - d_{d 1}) ζ_{1}, l_{7} = (1 - h_{d 1} - h_{d 2}) ζ_{2}, l_{8} = (1 - q_{d 1} - q_{d 2}) ζ_{3}, e_{1} = σ_{m} + μ, e_{2} = ψ_{m} + μ, e_{3} = α_{m} + μ, e_{4} = g_{m} + μ, e_{5} = y_{m} + μ, e_{6} = z_{m} + μ, e_{7} = r_{m} + μ, f_{1} = (1 - h_{m 1} - h_{m 2}) y_{m}, f_{2} = (1 - q_{m 1} - q_{m 2}) z_{m}

are as given in Section 2.2.1. Suppose that the linearized system

\begin{matrix} \dot{Z} = J ({\tilde{E}}_{1 s}) Z \end{matrix}

D.3

around the unique endemic equilibrium

{\tilde{E}}_{1 s}

has a solution of the form:

\begin{matrix} Z (t) = Z e^{ω t}, \end{matrix}

D.4

with

Z = (Z_{1}, Z_{2}, \dots, Z_{24}) \in C^{24} \ {0}, ω \in C

(i.e., assume that Z(t) is a nontrivial solution). Substituting the solution of the form (D.4) into the linearized system (D.3) gives

\begin{matrix} ω Z_{1} = & (- ξ_{f} - μ) Z_{1} - ξ_{f} \sum_{i = 2}^{14} Z_{i}, \\ ω Z_{2} = & - λ_{mf}^{* *} Z_{1} + (- λ_{mf}^{* *} - k_{1}) Z_{2} \\ - λ_{mf}^{* *} \sum_{i = 3}^{14} Z_{i} + a_{1} Z_{16} + a_{2} Z_{17} + a_{3} Z_{18} + a_{3} Z_{19} \\ + a_{4} Z_{20} + a_{4} Z_{21}, \\ ω Z_{3} = & σ_{f} Z_{2} - k_{2} Z_{3}, \\ ω Z_{4} = & l_{1} Z_{3} - k_{3} Z_{4}, \\ ω Z_{5} = & l_{2} Z_{4} - k_{4} Z_{5} + h_{u 2} y_{f} Z_{6}, \\ ω Z_{6} = & l_{3} Z_{5} - k_{5} Z_{6} + q_{u 2} z_{f} Z_{7}, \\ ω Z_{7} = & l_{4} Z_{6} - k_{6} Z_{7}, \\ ω Z_{8} = & d_{u 3} g_{f} Z_{5} - k_{7} Z_{8} + h_{d 2} ζ_{2} Z_{9}, \\ ω Z_{9} = & h_{u 3} y_{f} Z_{6} + l_{6} Z_{8} - k_{8} Z_{9} + q_{d 2} ζ_{3} Z_{10}, \\ ω Z_{10} = & q_{u 3} z_{f} Z_{7} + l_{7} Z_{9} - k_{9} Z_{10}, \\ ω Z_{11} = & l_{5} Z_{7} - k_{10} Z_{11}, \\ ω Z_{12} = & l_{8} Z_{10} + w_{f} Z_{11} - k_{11} Z_{12}, \\ ω Z_{13} = & r_{u} Z_{11} + r_{d} Z_{12} - μ Z_{13}, \\ ω Z_{14} = & b_{f} ψ_{f} Z_{3} + k_{f} α_{f} Z_{4} + d_{u 1} g_{f} Z_{5} + h_{u 1} y_{f} Z_{6} + q_{u 1} z_{f} Z_{7} \\ + d_{d 1} ζ_{1} Z_{8} + h_{d 1} ζ_{2} Z_{9} \\ + q_{d 1} ζ_{3} Z_{10} - μ Z_{14}, \\ ω Z_{15} = & (- ξ_{m} - μ) Z_{15} - ξ_{m} \sum_{i = 16}^{24} Z_{i}, \\ ω Z_{16} = & b_{1} Z_{2} + b_{2} Z_{3} + b_{3} Z_{4} + b_{3} Z_{5} + b_{4} Z_{6} + b_{4} Z_{7} + b_{5} Z_{8} + b_{6} Z_{9} + b_{6} Z_{10} \\ - λ_{fm}^{* *} Z_{15} + (- λ_{fm}^{* *} - e_{1}) Z_{16} - λ_{fm}^{* *} \sum_{i = 17}^{24} Z_{i}, \\ ω Z_{17} = & σ_{m} Z_{16} - e_{2} Z_{17}, \\ ω Z_{18} = & (1 - b_{m}) ψ_{m} Z_{17} - e_{3} Z_{18}, \\ ω Z_{19} = & (1 - k_{m}) α_{m} Z_{18} - e_{4} Z_{19} + h_{m 2} y_{m} Z_{20}, \\ ω Z_{20} = & (1 - d_{m}) g_{m} Z_{19} - e_{5} Z_{20} + q_{m 2} z_{m} Z_{21}, \\ ω Z_{21} = & f_{1} Z_{20} - e_{6} Z_{21}, \\ ω Z_{22} = & f_{2} Z_{21} - e_{7} Z_{22}, \\ ω Z_{23} = & r_{m} Z_{22} - μ Z_{23}, \\ ω Z_{24} = & b_{m} ψ_{m} Z_{17} + k_{m} α_{m} Z_{18} + d_{m} g_{m} Z_{19} + h_{m 1} y_{m} Z_{20} + q_{m 1} z_{m} Z_{21} - μ Z_{24}, \end{matrix}

D.5

where

a_{1} = \frac{β_{mf} c_{mf} η_{m} S_{f}^{* *}}{N_{m}^{*}}

,

a_{2} = \frac{β_{mf} c_{mf} S_{f}^{* *}}{N_{m}^{*}}

,

a_{3} = \frac{β_{mf} c_{mf} θ_{m} S_{f}^{* *}}{N_{m}^{*}}

,

a_{4} = \frac{β_{mf} c_{mf} θ_{m} κ_{m} S_{f}^{* *}}{N_{m}^{*}}

,

b_{1} = \frac{β_{fm} c_{mf} η_{f} S_{m}^{* *}}{N_{m}^{*}}

,

b_{2} = \frac{β_{fm} c_{mf} S_{m}^{* *}}{N_{m}^{*}}

,

b_{3} = \frac{β_{fm} c_{mf} θ_{f} S_{m}^{* *}}{N_{m}^{*}}

,

b_{4} = \frac{β_{fm} c_{mf} θ_{f} κ_{f} S_{m}^{* *}}{N_{m}^{*}}, b_{5} = \frac{β_{fm} c_{mf} θ_{f} ν S_{m}^{* *}}{N_{m}^{*}}, b_{6} = \frac{β_{fm} c_{mf} θ_{f} ν κ_{fd} S_{m}^{* *}}{N_{m}^{*}}

. The third to the 14th equations of the system (D.5) can be rewritten, by moving their negative terms to the left-hand side and dividing by their respective coefficients, as:

\begin{matrix} (1 + \frac{ω}{k_{2}}) Z_{3} = & \frac{σ_{f}}{k_{2}} Z_{2}, (1 + \frac{ω}{k_{3}}) Z_{4} = \frac{l_{1}}{k_{3}} Z_{3}, \\ (1 + \frac{ω}{k_{4}}) Z_{5} = & \frac{l_{2}}{k_{4}} Z_{4} + \frac{h_{u 2} y_{f}}{k_{4}} Z_{6}, (1 + \frac{ω}{k_{5}}) Z_{6} = \frac{l_{3}}{k_{5}} Z_{5} + \frac{q_{u 2} z_{f}}{k_{5}} Z_{7}, \\ (1 + \frac{ω}{k_{6}}) Z_{7} = & \frac{l_{4}}{k_{6}} Z_{6}, (1 + \frac{ω}{k_{7}}) Z_{8} = \frac{d_{u 3} g_{f}}{k_{7}} Z_{5} + \frac{h_{d 2} ζ_{2}}{k_{7}} Z_{9}, \\ (1 + \frac{ω}{k_{8}}) Z_{9} = & \frac{h_{u 3} y_{f}}{k_{8}} Z_{6} + \frac{l_{6}}{k_{8}} Z_{8} + \frac{q_{d 2} ζ_{3}}{k_{8}} Z_{10}, (1 + \frac{ω}{k_{9}}) \\ Z_{10} = \frac{q_{u 3} z_{f}}{k_{9}} Z_{7} + \frac{l_{7}}{k_{9}} Z_{9}, \\ (1 + \frac{ω}{k_{10}}) Z_{11} = & \frac{l_{5}}{k_{10}} Z_{7}, (1 + \frac{ω}{k_{11}}) Z_{12} = \frac{l_{8}}{k_{11}} Z_{10} + \frac{w_{f}}{k_{11}} Z_{11}, \\ (1 + \frac{ω}{μ}) Z_{13} = & \frac{r_{u}}{μ} Z_{11} + \frac{r_{d}}{μ} Z_{12}, \\ (1 + \frac{ω}{μ}) Z_{14} = & \frac{b_{f} ψ_{f}}{μ} Z_{3} + \frac{k_{f} α_{f}}{μ} Z_{4} + \frac{d_{u 1} g_{f}}{μ} Z_{5} + \frac{h_{u 1} y_{f}}{μ} \\ Z_{6} + \frac{q_{u 1} z_{f}}{μ} Z_{7} + \frac{d_{d 1} ζ_{1}}{μ} Z_{8} \\ + \frac{h_{d 1} ζ_{2}}{μ} Z_{9} + \frac{q_{d 1} ζ_{3}}{μ} Z_{10} . \end{matrix}

D.6

Similarly, the 17th through the 24th equations of the system (D.5) can be rewritten as:

\begin{matrix} \begin{matrix} (1 + \frac{ω}{e_{2}}) Z_{17} & = \frac{σ_{m}}{e_{2}} Z_{16}, (1 + \frac{ω}{e_{3}}) Z_{18} = \frac{(1 - b_{m}) ψ_{m}}{e_{3}} Z_{17}, \\ (1 + \frac{ω}{e_{4}}) Z_{19} & = \frac{(1 - k_{m}) α_{m}}{e_{4}} Z_{18} + \frac{h_{m 2} y_{m}}{e_{4}} Z_{20}, \\ (1 + \frac{ω}{e_{5}}) Z_{20} & = \frac{(1 - d_{m}) g_{m}}{e_{5}} Z_{19} + \frac{q_{m 2} z_{m}}{e_{5}} Z_{21}, \\ (1 + \frac{ω}{e_{6}}) Z_{21} & = \frac{f_{1}}{e_{6}} Z_{20}, (1 + \frac{ω}{e_{7}}) Z_{22} = \frac{f_{2}}{e_{7}} Z_{21}, \\ (1 + \frac{ω}{μ}) Z_{23} & = \frac{r_{m}}{μ} Z_{22}, \\ (1 + \frac{ω}{μ}) Z_{24} & = \frac{b_{m} ψ_{m}}{μ} Z_{17} + \frac{k_{m} α_{m}}{μ} Z_{18} + \frac{d_{m} g_{m}}{μ} Z_{19} + \frac{h_{m 1} y_{m}}{μ} Z_{20} + \frac{q_{m 1} z_{m}}{μ} Z_{21} . \end{matrix} \end{matrix}

D.7

To get a similar expression from the first and second equations of the system (D.5), all equations in (D.6) are rewritten in terms of

Z_{2}

as:

\begin{matrix} Z_{3} & = \frac{σ_{f}}{ω + k_{2}} Z_{2} & = : A_{1} Z_{2}, \\ Z_{4} & = \frac{σ_{f} l_{1}}{(ω + k_{2}) (ω + k_{3})} Z_{2} & = : A_{2} Z_{2}, \\ Z_{5} & = \frac{σ_{f} l_{1} l_{2} (ω + k_{5}) [ω^{2} + (k_{5} + k_{6}) ω + ι_{f 2}]}{(ω + k_{2}) (ω + k_{3}) B} Z_{2} & = : A_{3} Z_{2}, \\ Z_{6} & = \frac{σ_{f} l_{1} l_{2} l_{3} (ω + k_{6})}{(ω + k_{2}) (ω + k_{3}) B} Z_{2} & = : A_{4} Z_{2}, \\ Z_{7} & = \frac{σ_{f} l_{1} l_{2} l_{3} l_{4}}{(ω + k_{2}) (ω + k_{3}) B} Z_{2} & = : A_{5} Z_{2}, \\ Z_{8} & = (\frac{d_{u 3} g_{f} [ω^{2} + (k_{8} + k_{9}) ω + (ι_{f 3} + l_{6} h_{d 2} ζ_{2} k_{9}) / k_{7}]}{C}, A_{3}) \\ (+ \frac{(ω + k_{9}) h_{u 3} y_{f} h_{d 2} ζ_{2}}{C} A_{4} + \frac{q_{u 3} z_{f} q_{d 2} ζ_{3} h_{d 2} ζ_{2}}{C} A_{5}) Z_{2} & = : A_{6} Z_{2}, \\ Z_{9} & = \frac{(ω + k_{9}) h_{u 3} y_{f} A_{4} + q_{d 2} ζ_{3} q_{u 3} z_{f} A_{5} + (ω + k_{9}) l_{6} A_{6}}{ω^{2} + (k_{8} + k_{9}) ω + (ι_{f 3} + l_{6} h_{d 2} ζ_{2} k_{9}) / k_{7}} Z_{2} & = : A_{7} Z_{2}, \\ Z_{10} & = (\frac{q_{u 3} z_{f}}{ω + k_{9}} A_{5} + \frac{l_{7}}{ω + k_{9}} A_{7}) Z_{2} & = : A_{8} Z_{2}, \\ Z_{11} & = \frac{σ_{f} l_{1} l_{2} l_{3} l_{4} l_{5}}{(ω + k_{2}) (ω + k_{3}) B (ω + k_{10})} Z_{2} & = : A_{9} Z_{2}, \\ Z_{12} & = (\frac{l_{8}}{ω + k_{11}} A_{8} + \frac{w_{f}}{ω + k_{11}} A_{9}) Z_{2} & = : A_{10} Z_{2}, \\ Z_{13} & = (\frac{r_{u}}{ω + μ} A_{9} + \frac{r_{d}}{ω + μ} A_{10}) Z_{2} & = : A_{11} Z_{2}, \\ Z_{14} & = (\frac{b_{f} ψ_{f}}{ω + μ} A_{1} + \frac{k_{f} α_{f}}{ω + μ} A_{2} + \frac{d_{u 1} g_{f}}{ω + μ} A_{3} + \frac{h_{u 1} y_{f}}{ω + μ} A_{4} + \frac{q_{u 1} z_{f}}{ω + μ} A_{5}) \\ (+ \frac{d_{d 1} ζ_{1}}{ω + μ} A_{6} + \frac{h_{d 1} ζ_{2}}{ω + μ} A_{7} + \frac{q_{d 1} ζ_{3}}{ω + μ} A_{8}) Z_{2} & = : A_{12} Z_{2}, \end{matrix}

D.8

where,

\begin{matrix} B & = (ω + k_{4}) (ω + k_{5}) (ω + k_{6}) - (ω + k_{4}) l_{4} q_{u 2} z_{f} - l_{3} h_{u 2} y_{f} (ω + k_{6}) \\ = ω^{3} + (k_{4} + k_{5} + k_{6}) ω^{2} + (ι_{f 1} + k_{4} l_{4} q_{u 2} z_{f}) ω / k_{6} + (ι_{f 2} + k_{6} k_{4}) ω + ι_{f 1}, \\ C & = (ω + k_{7}) (ω + k_{8}) (ω + k_{9}) - (ω + k_{7}) l_{7} q_{d 2} ζ_{3} - l_{6} h_{d 2} ζ_{2} (ω + k_{9}) \\ = ω^{3} + (k_{7} + k_{8} + k_{9}) ω^{2} + (ι_{f 3} + k_{7} l_{7} q_{d 2} ζ_{3}) ω / k_{9} \\ + (ι_{f 3} + l_{6} h_{d 2} ζ_{2} k_{9}) ω / k_{7} + k_{9} k_{7} ω + ι_{f 3}, \end{matrix}

and

ι_{f 1} = k_{4} k_{5} k_{6} - l_{3} h_{u 2} y_{f} k_{6} - k_{4} l_{4} q_{u 2} z_{f}, ι_{f 2} = k_{5} k_{6} - l_{4} q_{u 2} z_{f}, ι_{f 3} = k_{7} k_{8} k_{9} - l_{6} h_{d 2} ζ_{2} k_{9} - k_{7} l_{7} q_{d 2} ζ_{3}

are as given in Section 2.2.1. It should be noted that

ι_{f 1}

,

ι_{f 2}

, and

ι_{f 3}

are non-negative (as shown in Appendix B). Hence (note also the aforementioned definitions of

k_{4}, k_{5}, k_{6}, k_{7}, k_{8}, k_{9}, l_{3}, l_{4}, l_{6}, l_{7}

), B and C have non-negative coefficients. Then, solving for

Z_{1}

from the first equation of the system (D.5) gives:

\begin{matrix} Z_{1} = - (\frac{ξ_{f}}{ω + ξ_{f} + μ}) \sum_{i = 2}^{14} Z_{i} . \end{matrix}

D.9

The second equation of (D.5) can be rewritten (by moving its negative terms to the left-hand side and dividing with

k_{1}

) as:

\begin{matrix} \frac{ω}{k_{1}} Z_{2} + \frac{λ_{mf}^{* *}}{k_{1}} Z_{1} + (\frac{λ_{mf}^{* *}}{k_{1}} + 1) Z_{2} \\ + \frac{λ_{mf}^{* *}}{k_{1}} \sum_{i = 3}^{14} Z_{i} = \frac{a_{1}}{k_{1}} Z_{16} + \frac{a_{2}}{k_{1}} Z_{17} + \frac{a_{3}}{k_{1}} (Z_{18} + Z_{19}) + \frac{a_{4}}{k_{1}} (Z_{20} + Z_{21}) . \end{matrix}

D.10

Then, substituting (D.8) and (D.9) into the left-hand side of (D.10) gives:

\begin{matrix} [1 + \frac{ω}{k_{1}} + \frac{λ_{mf}^{* *} (ω + μ)}{(ω + ξ_{f} + μ) k_{1}} (1 + \sum_{i = 1}^{12} A_{i})] Z_{2} \\ = \frac{a_{1}}{k_{1}} Z_{16} + \frac{a_{2}}{k_{1}} Z_{17} + \frac{a_{3}}{k_{1}} (Z_{18} + Z_{19}) + \frac{a_{4}}{k_{1}} (Z_{20} + Z_{21}) . \end{matrix}

D.11

To get a similar expression from the 15th and 16th equations of the system (D.5), the equations in (D.7) are rewritten in terms of

Z_{16}

as follows:

\begin{matrix} Z_{17} & = \frac{σ_{m}}{ω + e_{2}} Z_{16} & = : A_{13} Z_{16}, \\ Z_{18} & = \frac{σ_{m} (1 - b_{m}) ψ_{m}}{(ω + e_{2}) (ω + e_{3})} Z_{16} & = : A_{14} Z_{16}, \\ Z_{19} & = \frac{σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} [ω^{2} + (e_{5} + e_{6}) ω + ι_{m 2}]}{(ω + e_{2}) (ω + e_{3}) D} Z_{16} & = : A_{15} Z_{16}, \\ Z_{20} & = \frac{σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} (1 - d_{m}) g_{m} (ω + e_{6})}{(ω + e_{2}) (ω + e_{3}) D} Z_{16} & = : A_{16} Z_{16}, \\ Z_{21} & = \frac{σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} (1 - d_{m}) g_{m} f_{1}}{(ω + e_{2}) (ω + e_{3}) D} Z_{16} & = : A_{17} Z_{16}, \\ Z_{22} & = \frac{A_{17} f_{2}}{(ω + e_{7})} Z_{16} & = : A_{18} Z_{16}, \\ Z_{23} & = \frac{A_{17} f_{2} r_{m}}{(ω + e_{7}) (ω + μ)} Z_{16} & = : A_{19} Z_{16}, \\ Z_{24} & = (\frac{b_{m} ψ_{m}}{ω + μ} A_{13} + \frac{k_{m} α_{m}}{ω + μ} A_{14} + \frac{d_{m} g_{m}}{ω + μ}) \\ (A_{15} + \frac{h_{m 1} y_{m}}{ω + μ} A_{16} + \frac{q_{m 1} z_{m}}{ω + μ} A_{17}) Z_{16} & = : A_{20} Z_{16}, \end{matrix}

D.12

where,

\begin{matrix} D & = (ω + e_{4}) (ω + e_{5}) (ω + e_{6}) - (ω + e_{4}) f_{1} q_{m 2} z_{m} - (1 - d_{m}) g_{m} h_{m 2} y_{m} (ω + e_{6}) \\ = ω^{3} + (e_{4} + e_{5} + e_{6}) ω^{2} + (ι_{m 1} + e_{4} f_{1} q_{m 2} z_{m}) ω / e_{6} + (ι_{m 2} + e_{6} e_{4}) ω + ι_{m 1}, \end{matrix}

and

ι_{m 1} = e_{4} e_{5} e_{6} - (1 - d_{m}) g_{m} h_{m 2} y_{m} e_{6} - e_{4} f_{1} q_{m 2} z_{m}, ι_{m 2} = e_{5} e_{6} - f_{1} q_{m 2} z_{m}

are as given in Section 2.2.1. It should be noted that

ι_{m 1}

and

ι_{m 2}

are non-negative (as shown in Appendix B). Hence (note also the aforementioned definitions of

e_{4}, e_{5}, e_{6}, f_{1}

), all the coefficients of D are non-negative. Then, solving for

Z_{15}

from the 15th equation of the system (D.5) gives:

\begin{matrix} Z_{15} = - (\frac{ξ_{m}}{ω + ξ_{m} + μ}) \sum_{i = 16}^{24} Z_{i} . \end{matrix}

D.13

The 16th equation of (D.5) can also be rewritten (by moving its negative terms to the left-hand side and dividing with

e_{1}

) as:

\begin{matrix} \frac{ω}{e_{1}} Z_{16} + \frac{λ_{fm}^{* *}}{e_{1}} Z_{15} + (\frac{λ_{fm}^{* *}}{e_{1}} + 1) Z_{16} + \frac{λ_{fm}^{* *}}{e_{1}} \sum_{i = 17}^{24} Z_{i} = \frac{b_{1}}{e_{1}} Z_{2} + \frac{b_{2}}{e_{1}} Z_{3} + \frac{b_{3}}{e_{1}} (Z_{4} + Z_{5}) \\ + \frac{b_{4}}{e_{1}} (Z_{6} + Z_{7}) + \frac{b_{5}}{e_{1}} Z_{8} + \frac{b_{6}}{e_{1}} (Z_{9} + Z_{10}) . \end{matrix}

D.14

Then, substituting (D.12) and (D.13) into the left-hand side of (D.14) gives:

\begin{matrix} [1 + \frac{ω}{e_{1}} + \frac{λ_{fm}^{* *} (ω + μ)}{(ω + ξ_{m} + μ) e_{1}} (1 + \sum_{i = 13}^{20} A_{i})] Z_{16} = \frac{b_{1}}{e_{1}} Z_{2} + \frac{b_{2}}{e_{1}} Z_{3} + \frac{b_{3}}{e_{1}} (Z_{4} + Z_{5}) \\ + \frac{b_{4}}{e_{1}} (Z_{6} + Z_{7}) + \frac{b_{5}}{e_{1}} Z_{8} + \frac{b_{6}}{e_{1}} (Z_{9} + Z_{10}) . \end{matrix}

D.15

Hence, equations (D.6), (D.7), (D.11), and (D.15) can be rewritten in the following matrix-vector form (where

Y = (Z_{2}, Z_{3}, \dots, Z_{12}, Z_{16}, Z_{17}, \dots, Z_{22})

is the projection of the vector

Z

from

C^{24} \ {0}

to

C^{18} \ {0}

):

\begin{matrix} [1 + F_{i} (ω)] Y_{i} = {(H Y)}_{i}, i = 1, 2, \dots, 18, \end{matrix}

D.16

where,

\begin{matrix} \begin{matrix} F_{1} (ω) & = \frac{ω}{k_{1}} + \frac{λ_{mf}^{* *} (ω + μ)}{(ω + ξ_{f} + μ) k_{1}} (1 + \sum_{i = 1}^{12} A_{i}), F_{j} (ω) = \frac{ω}{k_{j}} (j = 2, 3, \dots, 11), \\ F_{12} (ω) & = \frac{ω}{e_{1}} + \frac{λ_{fm}^{* *} (ω + μ)}{(ω + ξ_{m} + μ) e_{1}} (1 + \sum_{i = 13}^{20} A_{i}), \\ F_{j} (ω) = \frac{ω}{e_{j - 11}} (j = 13, 14, \dots, 18), \end{matrix} \end{matrix}

and H is the non-negative

18 \times 18

matrix given below:

\begin{matrix} [\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{a_{1}}{k_{1}} & \frac{a_{2}}{k_{1}} & \frac{a_{3}}{k_{1}} & \frac{a_{3}}{k_{1}} & \frac{a_{4}}{k_{1}} & \frac{a_{4}}{k_{1}} & 0 \\ \frac{σ_{f}}{k_{2}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{l_{1}}{k_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{l_{2}}{k_{4}} & 0 & \frac{h_{u 2} y_{f}}{k_{4}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{l_{3}}{k_{5}} & 0 & \frac{q_{u 2} z_{f}}{k_{5}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{l_{4}}{k_{6}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{d_{u 3} g_{f}}{k_{7}} & 0 & 0 & 0 & \frac{h_{d 2} ζ_{2}}{k_{7}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{h_{u 3} y_{f}}{k_{8}} & 0 & \frac{l_{6}}{k_{8}} & 0 & \frac{q_{d 2} ζ_{3}}{k_{8}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{q_{u 3} z_{f}}{k_{9}} & 0 & \frac{l_{7}}{k_{9}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{l_{5}}{k_{10}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{l_{8}}{k_{11}} & \frac{w_{f}}{k_{11}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{b_{1}}{e_{1}} & \frac{b_{2}}{e_{1}} & \frac{b_{3}}{e_{1}} & \frac{b_{3}}{e_{1}} & \frac{b_{4}}{e_{1}} & \frac{b_{4}}{e_{1}} & \frac{b_{5}}{e_{1}} & \frac{b_{6}}{e_{1}} & \frac{b_{6}}{e_{1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{σ_{m}}{e_{2}} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{(1 - b_{m}) ψ_{m}}{e_{3}} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{(1 - k_{m}) α_{m}}{e_{4}} & 0 & \frac{h_{m 2} y_{m}}{e_{4}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{(1 - d_{m}) g_{m}}{e_{5}} & 0 & \frac{q_{m 2} z_{m}}{e_{5}} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{f_{1}}{e_{6}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{f_{2}}{e_{7}} & 0 \end{matrix}] . \end{matrix}

Furthermore, it can be shown that the projection of the endemic equilibrium

{\tilde{E}}_{1 s}

, defined by

\begin{matrix} {\tilde{E}}_{1 s}^{p} : = (E_{f}^{* *}, I_{f}^{* *}, P_{f}^{* *}, Q_{1 u}^{* *}, Q_{2 u}^{* *}, Q_{3 u}^{* *}, Q_{1 d}^{* *}, Q_{2 d}^{* *}, Q_{3 d}^{* *}, C_{u}^{* *}, C_{d}^{* *}, \\ E_{m}^{* *}, I_{m}^{* *}, P_{m}^{* *}, J_{1 m}^{* *}, J_{2 m}^{* *}, J_{3 m}^{* *}, C_{m}^{* *}), \end{matrix}

satisfies

\begin{matrix} {\tilde{E}}_{1 s}^{p} = H {\tilde{E}}_{1 s}^{p} . \end{matrix}

D.17

Since all components of

{\tilde{E}}_{1 s}^{p}

are strictly positive and

Y

is a nonzero vector (otherwise, (D.5) shows that

Z

is also a zero vector, which contradicts the assumption), there exists a minimal positive number s, depending on

Y

, such that

\begin{matrix} | Y | \leq s {\tilde{E}}_{1 s}^{p}, \end{matrix}

D.18

where

| Y | = (| Z_{2} |, | Z_{3} |, \dots, | Z_{12} |, | Z_{16} |, | Z_{17} |, \dots, | Z_{22} |)

and

| \cdot |

is the norm in

C

. Next, to show that Re

(ω) < 0

, we suppose the contrary. That is, we suppose that Re

(ω) \geq 0

, and consider the following two cases (Esteva and Vargas 2000; Pant and Gumel 2024).

Case 1: $ω = 0$ .

Since plugging

ω = 0

into (D.4) gives

Z (t) = Z

and hence

\dot{Z} (t) = 0

, it follows that (D.3) is a homogeneous linear system

J ({\tilde{E}}_{1 s}) Z = 0

, and its determinant is given by

\begin{matrix} Δ & = k_{1} k_{2} k_{3} ι_{f 1} ι_{f 3} e_{1} e_{2} e_{3} ι_{m 1} (1 - R_{vs}^{2} \frac{S_{m}^{* *} S_{f}^{* *}}{S_{m}^{*} S_{f}^{*}}) (μ + ξ_{f}) (μ + ξ_{m}) k_{10} k_{11} e_{7} μ^{4} \\ + [λ_{fm}^{* *} G_{m} k_{1} k_{2} k_{3} ι_{f 1} ι_{f 3} (μ + ξ_{f}) + λ_{mf}^{* *} G_{f} e_{1} e_{2} e_{3} ι_{m 1} (μ + ξ_{m})] k_{10} k_{11} e_{7} μ^{4} \\ + λ_{fm}^{* *} λ_{mf}^{* *} G_{m} G_{f}, \end{matrix}

D.19

where,

\begin{matrix} G_{m} & = e_{2} e_{3} ι_{m 1} μ + σ_{m} e_{3} ι_{m 1} (μ + b_{m} ψ_{m}) \\ + σ_{m} (1 - b_{m}) ψ_{m} ι_{m 1} (μ + k_{m} α_{m}) + σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} ι_{m 2} (μ + d_{m} g_{m}) \\ + σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} (1 - d_{m}) g_{m} {(e_{6} + f_{1}) μ \\ + f_{1} q_{m 1} z_{m} + h_{m 1} y_{m} e_{6} + f_{1} f_{2}}, \\ G_{f} & = k_{2} k_{3} ι_{f 1} ι_{f 3} μ + σ_{f} l_{1} l_{2} l_{3} l_{4} l_{5} ι_{f 3} + σ_{f} l_{1} l_{2} l_{8} j_{f 3} \\ + σ_{f} k_{3} ι_{f 1} ι_{f 3} (μ + b_{f} ψ_{f}) \\ + σ_{f} l_{1} ι_{f 1} ι_{f 3} (μ + k_{f} α_{f}) \\ + σ_{f} l_{1} l_{2} ι_{f 2} ι_{f 3} (μ + d_{u 1} g_{f}) + σ_{f} l_{1} l_{2} l_{3} (l_{4} + k_{6}) ι_{f 3} (μ + h_{u 1} y_{f}) \\ + σ_{f} l_{1} l_{2} {(j_{f 1} + j_{f 2} + j_{f 3}) μ + j_{f 1} d_{d 1} ζ_{1} + j_{f 2} h_{d 1} ζ_{2} + j_{f 3} q_{d 1} ζ_{3}} . \end{matrix}

To determine the sign of

Δ

, the system (D.2) is solved at the endemic equilibrium

{\tilde{E}}_{1 s}

as:

\begin{matrix} λ_{mf}^{* *} S_{f}^{* *} = k_{1} E_{f}^{* *}, \\ I_{f}^{* *} = \frac{σ_{f}}{k_{2}} E_{f}^{* *}, P_{f}^{* *} = \frac{σ_{f} l_{1}}{k_{2} k_{3}} E_{f}^{* *}, Q_{1 u}^{* *} = \frac{σ_{f} l_{1} l_{2} ι_{f 2}}{k_{2} k_{3} ι_{f 1}} E_{f}^{* *}, Q_{2 u}^{* *} = \frac{σ_{f} l_{1} l_{2} l_{3} k_{6}}{k_{2} k_{3} ι_{f 1}} E_{f}^{* *}, \\ Q_{3 u}^{* *} = \frac{σ_{f} l_{1} l_{2} l_{3} l_{4}}{k_{2} k_{3} ι_{f 1}} E_{f}^{* *}, Q_{1 d}^{* *} = \frac{σ_{f} l_{1} l_{2} j_{f 1}}{k_{2} k_{3} ι_{f 1} ι_{f 3}} E_{f}^{* *}, Q_{2 d}^{* *} = \frac{σ_{f} l_{1} l_{2} j_{f 2}}{k_{2} k_{3} ι_{f 1} ι_{f 3}} E_{f}^{* *}, \\ Q_{3 d}^{* *} = \frac{σ_{f} l_{1} l_{2} j_{f 3}}{k_{2} k_{3} ι_{f 1} ι_{f 3}} E_{f}^{* *}, \\ λ_{fm}^{* *} S_{m}^{* *} = e_{1} E_{m}^{* *}, \\ I_{m}^{* *} = \frac{σ_{m}}{e_{2}} E_{m}^{* *}, P_{m}^{* *} = \frac{σ_{m} (1 - b_{m}) ψ_{m}}{e_{2} e_{3}} E_{m}^{* *}, \\ J_{1 m}^{* *} = \frac{σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} ι_{m 2}}{e_{2} e_{3} ι_{m 1}} E_{m}^{* *}, \\ J_{2 m}^{* *} = \frac{σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} (1 - d_{m}) g_{m} e_{6}}{e_{2} e_{3} ι_{m 1}} E_{m}^{* *}, \\ J_{3 m}^{* *} = \frac{σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} (1 - d_{m}) g_{m} f_{1}}{e_{2} e_{3} ι_{m 1}} E_{m}^{* *} . \end{matrix}

Then, we observe that

\begin{matrix} λ_{mf}^{* *} S_{f}^{* *} & = a_{1} E_{m}^{* *} + a_{2} I_{m}^{* *} + a_{3} (P_{m}^{* *} + J_{1 m}^{* *}) + a_{4} (J_{2 m}^{* *} + J_{3 m}^{* *}) \\ = [a_{1} + a_{2} \frac{σ_{m}}{e_{2}} + a_{3} \frac{σ_{m} (1 - b_{m}) ψ_{m}}{e_{2} e_{3}} + a_{3} \frac{σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} ι_{m 2}}{e_{2} e_{3} ι_{m 1}}) \\ + (a_{4}, \frac{σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m} (1 - d_{m}) g_{m} (e_{6} + f_{1})}{e_{2} e_{3} ι_{m 1}}] E_{m}^{* *} \\ = e_{1} R_{m} \frac{S_{f}^{* *}}{S_{f}^{*}} E_{m}^{* *} = k_{1} E_{f}^{* *}, \\ λ_{fm}^{* *} S_{m}^{* *} & = b_{1} E_{f}^{* *} + b_{2} I_{f}^{* *} + b_{3} (P_{f}^{* *} + Q_{1 u}^{* *}) \\ + b_{4} (Q_{2 u}^{* *} + Q_{3 u}^{* *}) + b_{5} Q_{1 d}^{* *} + b_{6} (Q_{2 d}^{* *} + Q_{3 d}^{* *}) \\ = [b_{1} + b_{2} \frac{σ_{f}}{k_{2}} + b_{3} \frac{σ_{f} l_{1}}{k_{2} k_{3}} + b_{3} \frac{σ_{f} l_{1} l_{2} ι_{f 2}}{k_{2} k_{3} ι_{f 1}}) + b_{4} \frac{σ_{f} l_{1} l_{2} l_{3} (k_{6} + l_{4})}{k_{2} k_{3} ι_{f 1}} \\ + (b_{5} \frac{σ_{f} l_{1} l_{2} j_{f 1}}{k_{2} k_{3} ι_{f 1} ι_{f 3}} + b_{6} \frac{σ_{f} l_{1} l_{2} (j_{f 2} + j_{f 3})}{k_{2} k_{3} ι_{f 1} ι_{f 3}}] E_{f}^{* *} \\ = k_{1} R_{f} \frac{S_{m}^{* *}}{S_{m}^{*}} E_{f}^{* *} = e_{1} E_{m}^{* *} . \end{matrix}

Since

e_{1} E_{m}^{* *}, k_{1} E_{f}^{* *} > 0

and

{\tilde{R}}_{vs}^{2} = R_{m} R_{f}

, it follows that

1 - {\tilde{R}}_{vs}^{2} \frac{S_{m}^{* *} S_{f}^{* *}}{S_{m}^{*} S_{f}^{*}} = 0

, which means (D.19) can be simplified as

\begin{matrix} Δ = & [λ_{fm}^{* *} G_{m} k_{1} k_{2} k_{3} ι_{f 1} ι_{f 3} (μ + ξ_{f}) + λ_{mf}^{* *} G_{f} e_{1} e_{2} e_{3} ι_{m 1} (μ + ξ_{m})] k_{10} k_{11} e_{7} μ^{4} \\ + λ_{fm}^{* *} λ_{mf}^{* *} G_{m} G_{f} > 0 . \end{matrix}

Since the determinant (

Δ

) is positive, it follows that the system (D.5) has a unique solution, given by

Z = 0

, which contradicts the assumption that

Z

is not a zero vector.

Case 2: $ω \neq 0$ .

The objective here is to show that when

ω \neq 0

, then Re

(ω) \geq 0

gives a contradiction. It is convenient to define

\begin{matrix} F (ω) = min {| 1 + F_{i} (ω) |, i = 1, 2, \dots, 18} . \end{matrix}

Then, it follows, by taking absolute values in (D.16), that:

\begin{matrix} F (ω) | Y | \leq H | Y | . \end{matrix}

D.20

It can be shown (since

ω \neq 0

and Re

(ω) \geq 0

) that, in this case,

| 1 + F_{i} (ω) | > 1

for all

i = 1, 2, \dots, 18

, and therefore

F (ω) > 1

. It can be seen from (D.18) and (D.20) (and noting (D.17)), that

F (ω) | Y | \leq H | Y | \leq s H {\tilde{E}}_{1 s}^{p} = s {\tilde{E}}_{1 s}^{p}

, and hence

| Y | \leq \frac{s}{F (ω)} {\tilde{E}}_{1 s}^{p}

. Since

F (ω) > 1

, we obtain a contradiction to the minimality of s. Thus, Re

(ω) < 0

. This completes the proof.

□

Proofs of Theorems 3.1 and 3.3

Proof of Theorem 3.1

Proof

Consider the simplified model (3.1). Recall from Section 3 that the disease-free equilibrium (DFE) of the simplified model is given by:

\begin{matrix} E_{0 n} = (S_{f}^{*}, V_{f}^{*}, 0, 0, 0, 0, 0, 0, 0, S_{m}^{*}, V_{m}^{*}, 0, 0, 0, 0, 0, 0, 0), \end{matrix}

where, now,

\begin{matrix} S_{f}^{*} = \frac{(1 - q_{f}) π_{f}}{ξ_{f} + μ}, V_{f}^{*} = \frac{π_{f} (ξ_{f} + μ q_{f})}{μ (ξ_{f} + μ)}, \\ S_{m}^{*} = \frac{(1 - q_{m}) π_{m}}{ξ_{m} + μ}, V_{m}^{*} = \frac{π_{m} (ξ_{m} + μ q_{m})}{μ (ξ_{m} + μ)} . \end{matrix}

E.1

The local asymptotic stability of

E_{0 n}

is derived using the next generation operator method (Diekmann et al. 1990; Van den Driessche and Watmough 2002). The matrices F (for new infection terms near the DFE) and V (for linear transitions in the infected compartments) are given, respectively, by

F_{n} = [\begin{matrix} 0_{5 \times 5} & e_{1 n}^{T} f_{2 n} \\ e_{1 n}^{T} f_{1 n} & 0_{5 \times 5} \end{matrix}] and V_{n} = [\begin{matrix} V_{1 n} & 0_{5 \times 5} \\ 0_{5 \times 5} & V_{2 n} \end{matrix}],

where,

\begin{matrix} e_{1 n} & = (1, 0, 0, 0, 0), f_{1 n} = \frac{β_{fm} c_{mf} p_{m}^{*}}{N_{m}^{*}} (η_{f}, 1, {\tilde{θ}}_{f}, {\tilde{θ}}_{f} {\tilde{κ}}_{f}, 0), \\ f_{2 n} = \frac{β_{mf} c_{mf} p_{f}^{*}}{N_{m}^{*}} (η_{m}, 1, {\tilde{θ}}_{m}, {\tilde{θ}}_{m} {\tilde{κ}}_{m}, 0), \\ V_{1 n} & = [\begin{matrix} σ_{f} + μ & 0 & 0 & 0 & 0 \\ - σ_{f} & ψ_{f} + μ & 0 & 0 & 0 \\ 0 & - (1 - b_{f}) ψ_{f} & α_{f} + μ & 0 & 0 \\ 0 & 0 & - (1 - k_{f}) α_{f} & {\tilde{g}}_{f} + μ & 0 \\ 0 & 0 & 0 & (1 - d_{f}) {\tilde{g}}_{f} & r_{f} + \tilde{δ} + μ \end{matrix}], \\ V_{2 n} & = [\begin{matrix} σ_{m} + μ & 0 & 0 & 0 & 0 \\ - σ_{m} & ψ_{m} + μ & 0 & 0 & 0 \\ 0 & - (1 - b_{m}) ψ_{m} & α_{m} + μ & 0 & 0 \\ 0 & 0 & - (1 - k_{m}) α_{m} & {\tilde{g}}_{m} + μ & 0 \\ 0 & 0 & 0 & (1 - {\tilde{d}}_{m}) {\tilde{g}}_{m} & r_{m} + μ \end{matrix}], \end{matrix}

with,

p_{f}^{*} = S_{f}^{*} + (1 - ε_{v}) V_{f}^{*} = \frac{π_{f}}{μ} (1 - ε_{v} \frac{ξ_{f} + μ q_{f}}{ξ_{f} + μ})

and

p_{m}^{*} = S_{m}^{*} + (1 - ε_{v}) V_{m}^{*} = \frac{π_{m}}{μ} (1 - ε_{v} \frac{ξ_{m} + μ q_{m}}{ξ_{m} + μ}) .

It is convenient to define the quantity:

\begin{matrix} R_{vn} = ρ (F_{n} V_{n}^{- 1}) = \sqrt{R_{fn} R_{mn}}, \end{matrix}

where,

\begin{matrix} R_{fn} = \frac{β_{fm} c_{mf} μ p_{m}^{*}}{π_{m} (σ_{f} + μ)} \\ [η_{f} + \frac{σ_{f}}{ψ_{f} + μ} + \frac{{\tilde{θ}}_{f} σ_{f} (1 - b_{f}) ψ_{f}}{(ψ_{f} + μ) (α_{f} + μ)} + \frac{{\tilde{θ}}_{f} {\tilde{κ}}_{f} σ_{f} (1 - b_{f}) ψ_{f} (1 - k_{f}) α_{f}}{(ψ_{f} + μ) (α_{f} + μ) ({\tilde{g}}_{f} + μ)}], \\ R_{mn} = \frac{β_{mf} c_{mf} μ p_{f}^{*}}{π_{m} (σ_{m} + μ)} \\ [η_{m} + \frac{σ_{m}}{ψ_{m} + μ} + \frac{{\tilde{θ}}_{m} σ_{m} (1 - b_{m}) ψ_{m}}{(ψ_{m} + μ) (α_{m} + μ)} + \frac{{\tilde{θ}}_{m} {\tilde{κ}}_{m} σ_{m} (1 - b_{m}) ψ_{m} (1 - k_{m}) α_{m}}{(ψ_{m} + μ) (α_{m} + μ) ({\tilde{g}}_{m} + μ)}] . \end{matrix}

It follows from Theorem 2 of (Van den Driessche and Watmough 2002) that the DFE,

E_{0 n}

, of the simplified model (3.1) is locally-asymptotically stable if

R_{vn} < 1

, and unstable if

R_{vn} > 1

.

□

Proof of Theorem 3.3

Proof

Consider the special case of the simplified model (3.1) with no disease-induced mortality (i.e.,

\tilde{δ} = 0

); so that

N_{i} (t) \to \frac{π_{i}}{μ}

as

t \to \infty

for

i = f, m

. Let

R_{vn} < 1

. It is convenient to define

W_{i} (t) = \frac{S_{i} (t) + (1 - ε_{v}) V_{i} (t)}{π_{i} / μ}

and

W_{i}^{*} = \frac{S_{i}^{*} + (1 - ε_{v}) V_{i}^{*}}{π_{i} / μ},

for

i = f, m

. It can be seen from the bounds (3.5) and (3.6) that

\begin{matrix} W_{i} (t) \leq W_{i}^{*} for all t \geq 0 in D^{*} . \end{matrix}

E.2

Furthermore, let

x (t) = {(E_{f} (t), I_{f} (t), P_{f} (t), Q_{f} (t), C_{f} (t), E_{m} (t), I_{m} (t), P_{m} (t), Q_{m} (t), C_{m} (t))}^{T}

be the vector of the infected compartments of this special case of the simplified model (3.1). The equation for the rate of change of

x (t)

can then be re-written in terms of the next generation matrices

F_{n}

and

V_{n}^{*} : = V_{n} |_{\tilde{δ} = 0}

(with

F_{n}

and

V_{n}

matrices given in Section 3) as:

\begin{matrix} \dot{x} (t) = (F_{n} - V_{n}^{*}) x (t) - J_{n} x (t), \end{matrix}

E.3

where,

\begin{matrix} J_{n} = [\begin{matrix} 0_{5 \times 5} & e_{1 n}^{T} j_{2 n} \\ e_{1 n}^{T} j_{1 n} & 0_{5 \times 5} \end{matrix}], \end{matrix}

with,

e_{1 n} = (1, 0, 0, 0, 0), j_{1 n} = β_{fm} c_{mf} (W_{m}^{*} - W_{m}) (η_{f}, 1, {\tilde{θ}}_{f}, {\tilde{θ}}_{f} {\tilde{κ}}_{f}, 0),

and

j_{2 n} = β_{mf} c_{mf} (W_{f}^{*} - W_{f}) (η_{m}, 1, {\tilde{θ}}_{m}, {\tilde{θ}}_{m} {\tilde{κ}}_{m}, 0) .

Since

j_{1 n}

and

j_{2 n}

are non-negative vectors (from (E.2)), it follows that

J_{n}

is a non-negative matrix. Consequently, equation (E.3) can be written in terms of the following inequality:

\begin{matrix} \dot{x} (t) \leq (F_{n} - V_{n}^{*}) x (t) . \end{matrix}

E.4

Consider the linear ODE system given by the equality

\dot{x} (t) = (F_{n} - V_{n}^{*}) x (t)

in (E.4). If

R_{vn} < 1

, then

ρ (F_{n} {V_{n}^{*}}^{- 1}) < 1

(based on the local asymptotic stability result given in Theorem 3.1), which is equivalent to all eigenvalues of the Jacobian,

F_{n} - V_{n}^{*}

, of the linear system

\dot{x} (t) = (F_{n} - V_{n}^{*}) x (t)

, having negative real parts (Van den Driessche and Watmough 2002). Hence, the linearized differential inequality (E.4) is stable whenever

R_{vn} < 1

. Thus, it follows, by the standard comparison theorem (Lakshmikantham et al. 2016; Smith and Waltman 1995), that

x (t) \to 0

as

t \to \infty .

Substituting

x (t) = 0

into the equations of the special case of the simplified model (3.1) shows that

S_{i} (t) \to S_{i}^{*}

and

V_{i} (t) \to V_{i}^{*}

, as

t \to \infty

, for

i = f, m

. Hence, the disease-free equilibrium (

{\tilde{E}}_{0 n}

) of the special case of the simplified model (3.1) with negligible disease-induced mortality (i.e.,

\tilde{δ} = 0

) is globally-asymptotically stable in

D^{*}

, whenever

R_{vn} < 1

.

□

Monotonicity Verification for Sensitivity Analysis

Additional Simulations

Impact of vaccination of boys

Impact of vaccination of girls

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

ABG acknowledges the support, in part, of the National Science Foundation (Grant Number: DMS-2052363; transferred to DMS-2330801). The authors are grateful to the anonymous reviewers for their very constructive comments.

References

AlsalehAAGumelABAnalysis of Risk-Structured Vaccination Model for the Dynamics of Oncogenic and Warts-Causing HPV TypesIn: Bulletin of Mathematical Biology20147671670172610.1007/s11538-014-9972-425033777
ArbynMarcet al.Estimates of incidence and mortality of cervical cancer in 2018: a worldwide analysisIn: The Lancet Global Health202082e191e20310.1016/s2214-109x(19)30482-6PMC702515731812369
BaussanoIacopoet al.Impacts of human papillomavirus vaccination for different populations: A modeling studyIn: International Journal of Cancer201814351086109210.1002/ijc.31409PMC609933029603224
BlowerSally MDowlatabadiHadiSensitivity and Uncertainty Analysis of Complex Models of Disease Transmission: an HIV Model, as an ExampleIn: International Statistical Review199462229243
Bosch FX, de Sanjose S (2003) “Chapter 1: Human Papillomavirus and Cervical Cancer–Burden and Assessment of Causality”. In: JNCI Monographs 2003.31 , pp. 3–13. ISSN: 1745-6614. 10.1093/oxfordjournals.jncimonographs.a00347910.1093/oxfordjournals.jncimonographs.a00347912807939
BrauerFredBackward bifurcations in simple vaccination modelsIn: Journal of Mathematical Analysis and Applications2004298241843110.1016/j.jmaa.2004.05.045
Bray F, et al. “Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries”. In: CA: A Cancer Journal for Clinicians 74.3 (Apr. 2024), pp. 229–263. ISSN: 1542-4863. 10.3322/caac.2183410.3322/caac.2183438572751
Brewer DD, et al. “Prostitution and the Sex Discrepancy in Reported Number of Sexual Partners”. In: Proceedings of the National Academy of Sciences 97(22) (Oct. 24, 2000), pp. 12385–12388. 10.1073/pnas.210392097. (visited on 09/26/2025)10.1073/pnas.210392097PMC1735111027304
BulkmansNWJet al.High-risk HPV type-specific clearance rates in cervical screeningIn: British Journal of Cancer20079691419142410.1038/sj.bjc.6603653PMC236018317342094
CDC. CDC Vaccine Price List. Feb. 2024. https://www.cdc.gov/vaccines-for-children/php/awardees/current-cdc-vaccine-price-list.html
KCDC Center for Infectious Disease Control. Vaccination coverage for children nationwide in 2023. https://nip.kdca.go.kr/irhp/infm/goNatnVcntStatView.do. [Online; accessed 20-November-2024]. 2023
ChengLWangYDuJHuman Papillomavirus Vaccines: An Updated ReviewIn: Vaccines2020839110.3390/vaccines8030391PMC756529032708759
ChessonHarrell Wet al.The Estimated Lifetime Probability of Acquiring Human Papillomavirus in the United StatesIn: Sexually Transmitted Diseases2014411166066410.1097/olq.0000000000000193PMC674568825299412
Chin-Hong Peter V et al (2004) Age-Specific Prevalence of Anal Human Papillomavirus Infection in HIV-Negative Sexually Active Men Who Have Sex with Men: The EXPLORE Study. In: The Journal of Infectious Diseases 190(12):2070–2076. 10.1086/425906. (ISSN: 1537-6613.)10.1086/42590615551204
Choi JY, et al. “Human papillomavirus vaccine uptake in South Korea”. In: Clinical and Experimental Obstetrics & Gynecology 49(1) (Jan. 2022), p. 1. ISSN: 0390-6663. 10.31083/j.ceog4901022
ChoiWongyeongShimEunhaEvaluating the cost-effectiveness of HPV vaccination for adolescent girls in Japan: A comparison of 2-valent, 4-valent, and 9-valent HPV vaccines with consideration of cross-protectionIn: Preventive Medicine202417810774310.1016/j.ypmed.2023.10774337866695
Clifford G, et al. “Chapter 3: HPV type-distribution in women with and without cervical neoplastic diseases”. In: Vaccine 24 (Aug. 2006), S26–S34. ISSN: 0264-410X. 10.1016/j.vaccine.2006.05.02610.1016/j.vaccine.2006.05.02616950015
Cody P, et al. “Public health impact and cost effectiveness of routine and catch-up vaccination of girls and women with a nine-valent HPV vaccine in Japan: a model-based study”. In: BMC Infectious Diseases 21(1) (Jan. 2021). ISSN: 1471-2334. 10.1186/s12879-020-05632-010.1186/s12879-020-05632-0PMC778953933407188
DennyLAWrightTCHuman papillomavirus testing and screeningIn: Best Practice & Research Clinical Obstetrics & Gynaecology200519450151510.1016/j.bpobgyn.2005.02.00416150390
Devesa SS, et al. “Recent trends in cervix uteri cancer”. In: Cancer 64(10) (Nov. 1989), pp. 2184–2190. ISSN: 1097-0142. 10.1002/1097-0142(19891115)64:10<2184::aid-cncr2820641034>3.0.co;2-810.1002/1097-0142(19891115)64:10<2184::aid-cncr2820641034>3.0.co;2-82804908
Diekmann O, Heesterbeek JAP, Metz JAJ. “On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations”. In: Journal of Mathematical Biology 28(4) (June 1990). ISSN: 1432-1416. 10.1007/bf0017832410.1007/BF001783242117040
DushoffJonathanHuangWenzhangCastillo-ChavezCarlosBackwards bifurcations and catastrophe in simple models of fatal diseasesIn: Journal of Mathematical Biology199836322724810.1007/s0028500500999528119
Elbasha EH “Global Stability of Equilibria in a Two-Sex HPV Vaccination Model”. In: Bulletin of Mathematical Biology 70(3) (Nov. 2007). ISSN: 1522-9602. 10.1007/s11538-007-9283-010.1007/s11538-007-9283-017999117
ElbashaElamin HDasbachErik JInsingaRalph PModel for Assessing Human Papillomavirus Vaccination StrategiesIn: Emerging Infectious Diseases2007131284110.3201/eid1301.060438PMC272580117370513
ElbashaElamin HGalvaniAlison PVaccination against multiple HPV typesIn: Mathematical Biosciences200519718811710.1016/j.mbs.2005.05.00416095627
Elbasha EH, Gumel AB (2021) “Vaccination and herd immunity thresholds in heterogeneous populations”. In: Journal of Mathematical Biology 83.6–7 . ISSN: 1432-1416. 10.1007/s00285-021-01686-z10.1007/s00285-021-01686-zPMC865197934878609
Esteva L, Vargas C (2000) Influence of vertical and mechanical transmission on the dynamics of dengue disease. In: Mathematical Biosciences 167(1):51–64. 10.1016/s0025-5564(00)00024-9. (ISSN: 0025-5564.)10.1016/s0025-5564(00)00024-910942786
Fenigstein A, Preston M “The Desired Number of Sexual Partners as a Function of Gender, Sexual Risks, and the Meaning of “Ideal” ”. In: The Journal of Sex Research 44.1 (Apr. 1, 2007), pp. 89–95. ISSN: 0022-4499. 10.1080/00224490709336795. PMID: 17599267. (visited on 09/26/2025)10.1080/0022449070933679517599267
Gao S, et al. “A Dynamic Model to Assess Human Papillomavirus Vaccination Strategies in a Heterosexual Population Combined with Men Who have Sex with Men”. In: Bulletin of Mathematical Biology 83(1) (Jan. 2021). ISSN: 1522-9602. 10.1007/s11538-020-00830-y10.1007/s11538-020-00830-y33387083
GiulianoARet al.Age-Specific Prevalence, Incidence, and Duration of Human Papillomavirus Infections in a Cohort of 290 US MenIn: The Journal of Infectious Diseases2008198682783510.1086/59109518657037
GiulianoAnna Ret al.The Human Papillomavirus Infection in Men Study: Human Papillomavirus Prevalence and Type Distribution among Men Residing in Brazil, Mexico, and the United StatesIn: Cancer Epidemiology, Biomarkers & Prevention20081782036204310.1158/1055-9965.epi-08-0151PMC347177818708396
Global Strategy to Accelerate the Elimination of Cervical Cancer As a Public Health Problem. 1st ed. Description based on publisher supplied metadata and other sources. Geneva: World Health Organization, 2020. 156 pp. ISBN: 9789240014107
GoldieSue Jet al.Cost-effectiveness of HPV 16, 18 vaccination in BrazilIn: Vaccine200725336257627010.1016/j.vaccine.2007.05.05817606315
GoldieSue Jet al.Mathematical Models of Cervical Cancer Prevention in the Asia Pacific RegionIn: Vaccine200826M17M2910.1016/j.vaccine.2008.06.01818945411
Gumel AB, McCluskey CC, van den Driessche P. “Mathematical Study of a Staged-Progression HIV Model with Imperfect Vaccine”. In: Bulletin of Mathematical Biology 68(8) (July 2006), pp. 2105–2128. ISSN: 1522-9602. 10.1007/s11538-006-9095-710.1007/s11538-006-9095-716868850
GumelABCauses of backward bifurcations in some epidemiological modelsIn: Journal of Mathematical Analysis and Applications2012395135536510.1016/j.jmaa.2012.04.077
Ministry of Health and Welfare. Cancer Registration Statistics. 2021. https://kosis.kr/statHtml/statHtml.do?orgId=117&tblId=DT_117N_A00021 &conn_path=I2
HernandezBrenda Yet al.Transmission of Human Papillomavirus in Heterosexual CouplesIn: Emerging Infectious Diseases200814688889410.3201/eid1406.0706162PMC260029218507898
HethcoteHerbert WThe Mathematics of Infectious DiseasesIn: SIAM Review200042459965310.1137/s0036144500371907
HethcoteHerbert WThiemeHorst RStability of the endemic equilibrium in epidemic models with subpopulationsIn: Mathematical Biosciences198575220522710.1016/0025-5564(85)90038-0
Human papillomaviruses. This publication represents the views and expert opinions of an IARC Working Group on the Evaluation of Carcinogenic Risks to Humans, which met in Lyon, 6 - 13 June 1995. IARC monographs on the evaluation of carcinogenic risks to humans 64. Lyon: IARC, 1995. 551 pp. ISBN: 9283212649
JavameAliGumelAbba BQualitative Study of the Role of Pap Screening on HPV Transmission DynamicsIn: Journal of the Nigerian Mathematical Society201635580641
Kim SJ, et al (2018) Results from the School Entry Requirement Program for Immunization, 2018
Korea S. Total Survey Population by Administrative Area/Age/Gender. Population Census. (2000). https://kosis.kr/statHtml/statHtml.do?orgId=101&tblId=DT_1INOO &conn_path=I2
National Cancer Center Korea Central Cancer Registry. Annual report of cancer statistics in Korea in 2020. (2022)
Kwak K, Hwang S-s (2024) Predicted Cervical Cancer Prevention: Impact of National HPV Vaccination Program on Young Women in South Korea. In: Cancer Research and Treatment 56(3):898–908. 10.4143/crt.2023.981. (ISSN: 2005-9256.)10.4143/crt.2023.981PMC1126120738228083
KyoSet al.Detection Of High-Risk Human Papillomavirus In The Cervix And Semen Of Sex PartnersIn: Journal of Infectious Diseases1994170368268510.1093/infdis/170.3.6828077728
Bruni L, et al. Human Papillomavirus and Related Diseases in Republic of Korea. Summary Report. (2023)
Lakshmikantham V, Leela S, Martynyuk AA. Stability Analysis of Nonlinear Systems. 2nd ed. Systems and Control: Foundations and Applications Ser. Description based on publisher supplied metadata and other sources. Cham: Springer Basel AG, 2016. 1339 pp. ISBN: 9783319272009
Lee CH, et al. “Anal Human Papillomavirus Infection among HIV-Infected Men in Korea”. In: PLOS ONE 11.8 (Aug. 2016). Ed. by Maria Lina Tornesello, e0161460. ISSN: 1932-6203. 10.1371/journal.pone.016146010.1371/journal.pone.0161460PMC499336727548632
Lee H-K, et al. “Quality control program and its results of Korean Society for Cytopathologists”. In: The Korean Journal of Cytopathology 19(2) (2008), pp. 65–71
Lim S-C, Yoo CW (2019) Current Status of and Perspectives on Cervical Cancer Screening in Korea. In: Journal of Pathology and Translational Medicine 53(4):210–216. 10.4132/jptm.2019.04.11. (ISSN: 2383-7845.)10.4132/jptm.2019.04.11PMC663970131091859
Malik Tufail et al (2013) The impact of an imperfect vaccine and pap cytology screening on the transmission of human papillomavirus and occurrence of associated cervical dysplasia and cancer. In : Mathematical biosciences and engineering : MBE 10(4):1173–20510.3934/mbe.2013.10.117323906207
MarinoSimeoneet al.A methodology for performing global uncertainty and sensitivity analysis in systems biologyIn: Journal of Theoretical Biology2008254117819610.1016/j.jtbi.2008.04.011PMC257019118572196
Marra E, Lin C, Clifford GM. “Type-Specific Anal Human Papillomavirus Prevalence Among Men, According to Sexual Preference and HIV Status: A Systematic Literature Review and Meta-Analysis”. In: The Journal of Infectious Diseases 219 (4) (2018), pp. 590–598. ISSN: 0022-1899. 10.1093/infdis/jiy556. eprint: https://academic.oup.com/jid/article-pdf/219/4/590/27636038/jiy556.pdf10.1093/infdis/jiy55630239749
de MartelCatherineet al.Worldwide burden of cancer attributable to HPV by site, country and HPV typeIn: International Journal of Cancer2017141466467010.1002/ijc.30716PMC552022828369882
McLean AR, Blower SM. “Imperfect vaccines and herd immunity to HIV”. In: Proceedings of the Royal Society of London. Series B: Biological Sciences 253.1336 (July 1993), pp. 9–13. ISSN: 1471-2954. 10.1098/rspb.1993.007510.1098/rspb.1993.00758396781
MilwidRachael Met al.HPV Screening and Vaccination Strategies in an Unscreened Population: A Mathematical Modeling StudyIn: Bulletin of Mathematical Biology201881114313434210.1007/s11538-018-0425-329651669
MinKyung-Jinet al.The Korean guideline for cervical cancer screeningIn: Journal of Gynecologic Oncology201526323210.3802/jgo.2015.26.3.232PMC451034126197860
Moscicki A-B, et al. “Regression of low-grade squamous intra-epithelial lesions in young women”. In: The Lancet 364.9446 (Nov. 2004), pp. 1678–1683. ISSN: 0140-6736. 10.1016/s0140-6736(04)17354-610.1016/S0140-6736(04)17354-615530628
MoscickiAnna-Barbaraet al.The natural history of human papillomavirus infection as measured by repeated DNA testing in adolescent and young womenIn: The Journal of Pediatrics1998132227728410.1016/s0022-3476(98)70445-79506641
NakazawaHiroshiet al.Age at first sexual intercourse among young women with invasive cervical cancer: implications for routine vaccination against human papillomavirus in JapanIn: Japanese Journal of Clinical Oncology202353653053310.1093/jjco/hyad01736929021
World Health Organization. Human papillomavirus vaccines: WHO position paper. (2022)
PantBinodGumelAbba BMathematical assessment of the roles of age heterogeneity and vaccination on the dynamics and control of SARS-CoV-2In: Infectious Disease Modelling20249382887410.1016/j.idm.2024.04.007PMC1107946938725431
Podder C, Gumel A. “Transmission dynamics of a two-sex model for Herpes Simplex Virus type 2”. In: The Canadian Applied Mathematics Quarterly 17 (Jan. 2009)
ReidRichardet al.Should cervical cytologic testing be augmented by cervicography or human papillomavirus deoxyribonucleic acid detection?In: American Journal of Obstetrics and Gynecology199116461461147110.1016/0002-9378(91)91425-v1646568
SafdarSalmanNgonghalaCalistus NGumelAbba BMathematical assessment of the role of waning and boosting immunity against the BA.1 Omicron variant in the United StatesIn: Mathematical Biosciences and Engineering202220117921210.3934/mbe.202300936650762
Saltelli A, ed. Global Sensitivity Analysis: The Primer. Chichester, England Hoboken, NJ: John Wiley, (2008). 10.1002/9780470725184
Schiffman M, et al. “Human papillomavirus and cervical cancer”. In: The Lancet 370.9590 (Sept. 2007), pp. 890–907. ISSN: 0140-6736. 10.1016/s0140-6736(07)61416-010.1016/S0140-6736(07)61416-017826171
SharomiOet al.Role of incidence function in vaccine-induced backward bifurcation in some HIV modelsMath Biosci200721043646310.1016/j.mbs.2007.05.01217707441
ShinHai-Rimet al.Prevalence and Determinants of Genital Infection with Papillomavirus, in Female and Male University Students in Busan, South KoreaIn: The Journal of Infectious Diseases2004190346847610.1086/42127915243918
Shin HY, et al. “Trends in cervical cancer screening rates among Korean women: results of the Korean National Cancer Screening Survey, 2005–2020”. In: Journal of Gynecologic Oncology 33(4) (2022). ISSN: 2005-0399. 10.3802/jgo.2022.33.e3910.3802/jgo.2022.33.e39PMC925086335320888
Smith HL, Waltman P (1995) The Theory of the Chemostat: Dynamics of Microbial Competition. Cambridge University Press. 10.1017/cbo9780511530043
SmithJSet al.Human papillomavirus type distribution in invasive cervical cancer and high-grade cervical lesions: A meta-analysis updateIn: International Journal of Cancer2007121362163210.1002/ijc.2252717405118
Thieme HR. “Local Stability in Epidemic Models for Heterogeneous Populations”. In: Mathematics in Biology and Medicine. Springer Berlin Heidelberg, 1985, pp. 185–189. ISBN: 9783642932878. 10.1007/978-3-642-93287-8_26
Tyson G, et al. A First Look at User Activity on Tinder. July 7, 2016. 10.48550/arXiv.1607.01952. arXiv: 1607.01952[cs]. http://arxiv.org/abs/1607.01952 (visited on 09/26/2025). Pre-published
Van de Velde N, Brisson M, Boily M-C. “Modeling Human Papillomavirus Vaccine Effectiveness: Quantifying the Impact of Parameter Uncertainty”. In: American Journal of Epidemiology 165.7 (Feb. 2007), pp. 762–775. ISSN: 1476-6256. 10.1093/aje/kwk05910.1093/aje/kwk05917276976
Van den Driessche P, Watmough J. “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission”. In: Mathematical biosciences 180.1-2 (2002), pp. 29–4810.1016/s0025-5564(02)00108-612387915
Van DoornumGJet al.Regional distribution and incidence of human papillomavirus infections among heterosexual men and women with multiple sexual partners: a prospective studyIn: Sexually Transmitted Infections199470424024610.1136/sti.70.4.240PMC11952477959707
Walboomers JMM, et al. “Human papillomavirus is a necessary cause of invasive cervical cancer worldwide”. In: The Journal of Pathology 189 (1) (1999), pp. 12–19. 10.1002/(SICI)1096-9896(199909)189:1<12::AID-PATH431>3.0.CO;2-F10.1002/(SICI)1096-9896(199909)189:1<12::AID-PATH431>3.0.CO;2-F10451482

Author Contributions

Sigrid Bergström: Methodology

Declarations

Funding

Funding was provided by the European Research Council under the European Union’s Horizon Europe programme (grant agreement 101054321).

Competing Interests

The authors declare no competing interests.

Data Availability

The datasets generated and analyzed during the study are available from the corresponding author on reasonable request.

Demonstration content: the full text of this article is adapted from “Mathematical Assessment of the Roles of Vaccination and Pap Screening on the Burden of HPV and Related Cancers in Korea” by Park and colleagues (Bulletin of Mathematical Biology, 2025; https://doi.org/10.1007/s11538-025-01548-5), which is available under a Creative Commons license (https://creativecommons.org/licenses/by/4.0/).

Copyright and License

License: https://creativecommons.org/licenses/by-nc/4.0/

Download

PDF JATS XML

Abstract

Introduction

Formulation of Vaccination and Pap Screening Model

Basic qualitative properties

Lemma 2.1

Lemma 2.2

Proof

Existence and local asymptotic stability of equilibria of the model (2.5)

Disease-free equilibrium (DFE)

Theorem 2.3

Theorem 2.4

Existence and local asymptotic stability of endemic equilibria: special case

Theorem 2.5

Theorem 2.6

Calibration and Analysis of Simplified Model: No Screening

Theorem 3.1

Lemma 3.2

Proof

Theorem 3.3

Data fitting and parameter estimation for the simplified model

Sensitivity Analysis for the Model (2.5)

Numerical Simulation and Assessment of Control Strategies

Disease burden corresponding to each strategy

Impact of vaccination and pap screening on cancer cases and mortality

Spillover effect of vaccination-only strategy

Discussion and Conclusion

Appendices

Proof of Lemma 2.1

Proof

Non-negativity of the Control Reproduction Number (Rvs)

Coefficients of the Polynomial (2.20)

Proof of Theorem 2.6

Proof

Proofs of Theorems 3.1 and 3.3

Proof of Theorem 3.1

Proof

Proof of Theorem 3.3

Proof

Monotonicity Verification for Sensitivity Analysis

Additional Simulations

Impact of vaccination of boys

Impact of vaccination of girls

Acknowledgements

References

Author Contributions

Declarations

Funding

Competing Interests

Data Availability

Copyright and License

Non-negativity of the Control Reproduction Number ( $R_{vs}$ )