THE IMPACT OF AN IMPERFECT VACCINE AND PAP CYTOLOGY SCREENING ON THE TRANSMISSION OF HUMAN PAPILLOMAVIRUS AND OCCURRENCE OF ASSOCIATED CERVICAL DYSPLASIA AND CANCER

A mathematical model for the natural history of human papillomavirus (HPV) is designed and used to assess the impact of a hypothetical anti-HPV vaccine and Pap cytology screening on the transmission dynamics of HPV in a population. Rigorous qualitative analysis of the model reveals that it undergoes the phenomenon of backward bifurcation. It is shown that the backward bifurcation is caused by the imperfect nature of the HPV vaccine or the HPV-induced and cancer-induced mortality in females. For the case when the disease-induced and cancer-induced mortality is negligible, it is shown that the disease-free equilibrium (i.e., equilibrium in the absence of HPV and associated dysplasia) is globally-asymptotically stable if the associated reproduction number is less than unity. The model has a unique endemic equilibrium when the reproduction threshold exceeds unity. The unique endemic equilibrium is globally-asymptotically stable for a special case, where the associated HPVinduced and cancer-induced mortality is negligible. Numerical simulations of the model, using a reasonable set of parameter values, support the recent recommendations by some medical agencies and organizations in the USA to offer Pap screening on a 3-year basis (rather than annually).

1. Introduction.Cervical cancer, the second most frequent form of cancer among women worldwide, is an important global public health problem [18].In 2002, an estimated 529,000 new cases of cervical cancer occurred worldwide, representing nearly 10% of all cases of cancer among women [18].An estimated 1,300 Canadian women were diagnosed with invasive cervical cancer in 2009, corresponding to an annual incidence rate of 7 per 100,000 women [6].
There is overwhelming epidemiologic and experimental evidence that cervical infection with certain "high-risk" types of the human papillomavirus (HPV) is responsible for virtually all cervical cancer cases [20,54].Genital HPV infection is the most common sexually-transmitted disease in the US and Canada [20].It is estimated that up to 75% of women may eventually become infected in their lifetime [29] (some other studies, such as those reported in [31,41,44], show that up to 70% of sexually-active men and women acquire HPV infection at least once in their lifetime).The HPV infection is prevalent among younger women, with evidence of infection found in 5-40% of asymptomatic women in the reproductive age group.Fortunately, most of these infections are transient (last less than a year) and asymptomatic.However, persistent infection with thirteen or so high-risk (oncogenic) HPV (HR-HPV) types could lead to cancer [37].
The natural history of cervical cancer is perhaps the best understood of all malignancies.Cervical carcinogenesis is a complex stepwise process characterized by slow progression over a continuum of increasingly more severe precancerous changes known collectively as cervical intraepithelial neoplasia (CIN) [20].Depending on the extent and severity of dysplastic features, the spectrum of CIN is traditionally divided into three histopathological categories, namely: CIN1, CIN2 and CIN3.In CIN1, cells with malignant changes are limited to the superficial layer of the cervical epithelium.Most CIN1 lesions are likely to disappear without treatment.However, a small percentage may progress to high-grade CIN (CIN2 and CIN3).This (high) grade CIN is characterized by more severe dysplastic changes and higher degree of epithelial basal cell involvement.The risk of progression to invasive cervical cancer increases significantly with worsening CIN grades [29].
In most developed countries, Pap cytology screening for the early detection of cervical neoplasia has been successful in reducing cervical cancer incidence and mortality, especially in jurisdictions with organized screening programs [10].However, Pap cytology, as a primary screening tool, has important limitations.A single Pap smear has both low sensitivity (approximately 51% [35]) and poor specificity which necessitates repeated screening at regular intervals throughout adulthood.In many countries, women are screened annually starting at the initiation of sexual activity [47].The present policy in the Canadian province of Manitoba requires all women 18 years or older to undergo annual screening with conventional Pap cytology if they have ever been sexually-active.After 3 consecutive normal smears, screening is continued at 2-year intervals [43].This policy could result in a woman undergoing 27 routine smears in her lifetime assuming full compliance and no abnormal results.Some countries offer Pap screening at a less frequent interval (see, for instance, [25,30] for the screening guidelines in Norway).Furthermore, some agencies in the USA (such as American Cancer Society, American Society for Colposcopy and Cervical Pathology, and American Society for Clinical Pathology [46]) have recently recommended Pap screening to be carried out every three years.
The recent introduction of effective vaccines against HPV is expected to reduce the burden of cervical cancer substantially [23].However, screening for cervical cancer will still be needed because about 30% of cervical cancers are caused by types other than the two high-risk HPV (HR-HPV) types targeted by the current vaccines (Types 16 and 18) [38], and because women who are already infected will continue to develop cervical cancer [19].In the future, the cost-effectiveness of Pap screening may be significantly reduced because as the prevalence of cervical neoplasia decreases, the positive predictive value of the Pap test will also decrease, and as a result, more women will be referred for unnecessary diagnostic procedures and follow-up [19].
Mathematical modeling is a useful tool for assessing the potential impact of intervention strategies against HPV infection (and the associated cervical cancer), such as mass vaccination and Pap screening.A number of authors have reported on the use of such modeling to evaluate the impact of HPV vaccine [1,2,4,9,13,16].In particular, Al-arydah [2] developed a two-sex, age-structured model to describe the vaccination program for an HPV vaccine in childhood (under 13 years) and adult stages.It is shown that vaccinating a single age cohort in one gender can result in eventual control of HPV across all age groups.Furthermore, Elbasha and Galvani [13] showed that for the case of synergistic interaction between HPV types, the use of mass vaccination may reduce the prevalence of HPV types that are not included in the currently available HPV vaccines.Models for the combined impact of HPV vaccination and Pap screening have also been published in the literature (see, for instance, [15,17,24,39,40]).Myers et al. [39] proposed a model for the natural history of HPV infection and cervical carcinogenesis.
The purpose of the current study is to extend some of the aforementioned studies by designing a new and comprehensive sex-structured (male and female) model for the natural history of cervical cancer, and use the model to assess the public health impact of mass vaccination and Pap screening of sexually-active females.Although there are many oncogenic HPV types, this theoretical study only considers the dynamics and control of the two vaccine-targeted oncogenic types (Types 16 and 18).Furthermore, the only modality of screening considered in this study is Pap screening (and not HPV DNA testing).The model is formulated in Section 2 and analysed in Section 3. The impact of vaccination and screening are assessed in Section 4.
2. Model formulation.The model is formulated as follows.The total sexuallyactive female population at time t, denoted by N f (t), is split into the sub-populations of unvaccinated susceptible (S f (t)), susceptible individuals vaccinated against the relevant vaccine-targeted types (V (t)), HPV-infected (I f (t)), HPV-infected with persistent infection (P (t)), HPV-infected individuals with undetected CIN (three grades, Q i (t); i = 1, 2, 3, representing low-, medium-and high-grade squamous intraepithelial lesions, respectively), infected individuals with detected CIN (three grades, Q id ; i = 1, 2, 3), those with undetected cancer (C(t)), those with detected cancer (C d (t)), those who recovered from infection or the CIN without developing cancer (R f (t)) and those who recovered from cancer (R c (t)), so that The total sexually-active male population, denoted by N m (t), is sub-divided into susceptible (S m (t)), infected (I m (t)) and recovered (R m (t)) males.Hence, The population of unvaccinated susceptible females (S f ) is generated by the recruitment of new sexually-active females (at a rate π f ; a fraction, φ, of which is vaccinated).Although this study does not explicitly model the various HPV types, it is assumed that the theoretical vaccine targets some HPV types.The vaccines currently in the market are GSK's bivalent Cervarix (which targets HPV Types 16 and 18) and Merck's quadrivalent Gardasil (which targets HPV Types 6,11,16,18).These vaccines, which have been approved for use in a number of countries, including Canada [44], have 90-100% efficacy against infection with Types 16 and 18 [22,53].
The unvaccinated susceptible female population is decreased by the acquisition of HPV infection, following effective contact with infected males, at a rate λ m , given by where β m is the proportion of potentially infectious contacts which result in the transmission of infection from males to females and c f is the average number of contacts per female per unit time.It is also decreased due to natural death (at a rate µ f ; females in all compartments suffer natural death at this rate, µ f ).The class of vaccinated susceptible females (V ) is populated by the vaccination of new sexually-active females (at the aforementioned rate π f φ).It is decreased by HPV infection (at a reduced rate (1 − v )λ m , where 0 < v ≤ 1 represents the vaccine efficacy against HPV infection) and natural death.It is assumed that the vaccine does not wane during the period under consideration [16].The population of infected females (I f ) is generated by the infection of unvaccinated and vaccinated susceptible females.It is assumed that a fraction, r 1 , of the infected female population recovers (at a rate σ f r 1 ), while the remaining fraction, 1 − r 1 , develops persistent infection (and moves to the class P at the rate σ f (1 − r 1 )).It is assumed that recovery confers permanent immunity against re-infection with HPV.It should be mentioned that although re-infection (especially with a different HPV strain) is possible in HPV dynamics, the model (3) does not incorporate re-infection.The reader is referred to [13] for a discussion of immunity against HPV.The class of infected females is further reduced by natural death and infection-induced death (at a rate δ f ).
The population of infected females with persistent infection (P ) is generated by the fraction (1 − r 1 ) of infected females in the I f class that develop persistent infection (at the rate (1 − r 1 )σ f ).It is decreased by progression to Grade 1 of undetected CIN (at a rate κ p ), recovery (at a rate ψ 0 ) and natural death.The population of infected females in Grade 1 of undetected CIN (CIN1) is populated by females with persistent infection (at the rate κ p ) and those regressing from CIN Grade 2 (at a rate κ q ), and is diminished by recovery (at a rate ψ 1 ), progression to undetected CIN Grade 2 (at a rate ζ 1 ), detection (at a rate α 1 ) and natural death.The class of infected females in Grade 2 of undetected CIN (CIN2) is generated by progression from CIN Grade 1 (at the rate ζ 1 ), and is decreased by remission to CIN Grade 1 (at a rate κ q ), recovery (at a rate ψ 2 ), progression to Grade 3 (at a rate ζ 2 ), detection (at a rate α 2 ) and natural death.Similarly, the population of those in undetected CIN Grade 3 (CIN3) is generated at the rate ζ 2 , and is diminished by recovery (at a rate ψ 3 ), development of cervical cancer (at a rate ζ 3 ), detection (at a rate α 3 ) and natural death.
The population of females with detected CIN Grade 1 (Q 1d ) is populated by the detection of undetected infected females in class Q 1 (at the rate α 1 ) and by the regression of those in Q 2d class (at a rate κ d ).It is decreased by progression to detected CIN Grade 2 (at the rate ζ 1 ), recovery (at a rate ρ 1 ) and natural death.The class of detected CIN Grade 2 (Q 2d ) is generated by the detection of infected individuals in the Q 2 class (at the rate α 2 ) and by progression of those in the Q 1d class (at the rate ζ 1 ).It is decreased by remission to Q 1d class (at the rate κ d ), progression to Grade 3 (at the rate ζ 2 ), recovery (at a rate ρ 2 ) and natural death.Similarly, the population of those in detected CIN Grade 3 class (Q 3d ) is generated at the rates α 3 (for those in Q 3 class) and ζ 2 (for those females progressing from Q 2d class), and is decreased by recovery (at a rate ρ 3 ) and natural death.
The detection rate, α i i = 1, 2, 3, is an aggregate parameter that depends on the proportion of sexually-active females screened for cervical dysplasia (χ -henceforth referred to as 'proportion screened'), the screening sensitivity of the Pap test in the respective CIN grade ( i , i = 1, 2, 3) and the screening frequency per year (∇).Thus for the present model, it is defined as Van de Velde et al. [52] assume the Pap test to be 60% sensitive for the detection of CIN Grade 1 and 78% sensitive for Grades 2 and 3.They assume the proportion screened (χ) to range from a low 7% (in the age group 70-100) to 32% (in the age group 20-29) per year.The proportions screened vary from 37% in the Canadian provinces of British Columbia and Ontario to 44% in Nova Scotia [26].As the base-case Pap test assumption, the current model assumes χ = 0.32.Due to the wide range of values reported in the literature, the sensitivity of the model to the proportion screened (χ) is examined in Section 4.4.
The population of individuals with cervical cancer (C) is generated by individuals in the Q 3 class who develop cancer (at the rate ζ 3 ).It is decreased by cancer detection (at a rate α c ).If undetected, the population of individuals with cervical cancer suffers an additional cancer-related mortality (at a rate δ c ).The population of females with detected cervical cancer (C d ) is populated by the detection of females in the cancer class (at the rate α c ) and decreased by recovery (at a rate γ), natural death, and cancer-related death (at a rate δ c d ).The population of infected females who recovered from cervical cancer (R c ) is generated by the recovery of females with detected cancer (at the rate γ) and is decreased by natural death.Similarly, the population of those who recovered from HPV infection and CIN (R f ) is generated by the recovery of infected individuals in the , and is decreased by natural death.
The population of susceptible males (S m ) is generated by the recruitment of new sexually-active males (at a rate π m ).It is decreased by infection, following effective contact with infected females (it is assumed that only infected females in the compartments I f and P can transmit infection to males) at a rate λ f , given by where β f is the proportion of potentially infectious contacts which result in the transmission of infection from females to males and c m is the average number of contacts per male per unit time.
The modification parameter θ, 0 < θ ≤ 1, accounts for the variability in the transmission probability of individuals in the P class in relation to those in the I f class.It is assumed that infected females with CIN i (i = 1, 2, 3) are not infectious.Males in all epidemiological compartments suffer natural death at a rate µ m .The population of infected males is generated at the rate λ f , and diminished by recovery (at a rate σ m ) and natural death.No persistent infection or disease-related mortality is assumed in males.
Furthermore, the following conservation law (for number of sexual contacts made by males balancing those made by females) [16] is assumed to hold: Male sexual partners are assumed to be abundant, so that females can always have enough number of sexual partners per unit time.Hence, it is assumed that c f is constant, and c m is calculated from the relation Nm .Using (1), the infection rates λ m and λ f are given by Based on the above formulations and assumptions, and using the conservation law (1), it follows that the model for the natural history of HPV, in the presence of an imperfect vaccine and Pap cytology screening, is given by the following deterministic system of non-linear differential equations (a flow diagram of the model is depicted in Figure 1, and the associated variables and parameters of the model are described in Tables 1 and 2, respectively): The model ( 3) is a deterministic model for the natural history of the cervical cancer, inspired by the Markov model in [39], incorporating sex structure (the dynamics of sexually-active males and females) as well as a mass vaccination program (against some vaccine-targeted HPV types).Furthermore, it extends the model in [13] by incorporating sex structure and by accounting for persistent HPV infections, cervical dysplasia and cancer, and Pap screening cytology.It also extends the models in [16,28] by incorporating the dynamics of individuals with persistent infection.
The model (3) will now be qualitatively analysed to gain insight into its dynamical features.It should be further emphasized that the model (3) only considers the vaccine-targeted HPV types.
3. Analysis of the model.
Population of unvaccinated susceptible females V Population of vaccinated susceptible females I f Population of infected females P Population of females with persistent HPV infection Total male population Table 1.Description of state variables of the model (3).
is positively-invariant and attracting with respect to the model (3).
Proof.Adding the first 14 equations of the model (3) gives Since N f (t) ≥ 0, it follows using a standard comparison theorem [33] that Similarly, adding the last three equations of the model (3) gives: so that, Thus, the positive invariance of D has been proved.
To prove that D is attracting, it is clear from ( 4) and ( 5) that (for i = m, f ), dN i dt < 0, whenever N i (t) > π i /µ i .Thus, either the solution enters D in finite time, or N i (t) approaches π i /µ i , and the variables denoting infected classes approach zero.Hence, D is attracting and all solutions in R 17 + eventually enter D.
Since the region D is is positively-invariant and attracting, the usual existence, uniqueness, continuation results hold for the system, and the system (3) is mathematically and epidemiologically well-posed in D [27].Therefore, it is sufficient to consider the dynamics of the flow generated by the model (3) in D. Table 2. Description of the parameters of the model (3).A: Assumed.D: Derived.
Lemma 3.2 implies that the combined use of Pap screening and an HPV vaccine can lead to the effective control of the disease in the population if the threshold quantity (R v ) can be brought to (and maintained at) a value less than unity.Mathematically-speaking, the result implies that the vaccine-targeted HPV types can be eliminated from the population (when R v < 1) if the initial sizes of the sub-populations of the model ( 3) are in the basin of attraction of the DFE (E 0 ).

Existence of endemic equilibrium point (EEP).
Let ) represents an arbitrary endemic equilibrium of the model (3) (i.e., equilibria for which at least one of the infected components of the model ( 3) is non-zero).Furthermore, let be the forces of infection for males and females at steady-state, respectively.Solving the equations of the model (3) at steady-state gives: Substituting ( 9) in (8), and simplifying, gives, respectively, and, It follows, by substituting ( 11) into (10), that the positive (endemic) equilibria of the model (3) satisfy the following quadratic (in terms of λ * * m ): where, ].The components of the positive endemic equilibria of the model ( 3) are then obtained by solving the quadratic (12) for λ * * m and substituting the results in ( 9) and (11).It should be observed that (11) gives λ * * f as a rational function of λ * * m , so that a fixed value of λ * * m corresponds to a unique value of λ * * f .Clearly, the coefficient a 0 of ( 12) is always positive (since 0 < v ≤ 1), and c 0 is positive (negative) if R v is less than (greater than) unity.Thus, the following result is established.It is clear from Case 1 of Theorem 3.3 that the model has a unique endemic equilibrium whenever R v > 1. Numerical simulation results, depicted in Figure 2, show convergence of solutions to this EEP when R v > 1 suggesting that the EEP is asymptotically-stable when it exists.Furthermore, Case 3 indicates the possibility of backward bifurcation (where the locally-asymptotically stable DFE co-exists with a locally-asymptotically stable endemic equilibrium when R v < 1 [3,14,32,48]).The existence of backward bifurcation in vaccination models, such as (3), has been established in a number of epidemiological settings (see, for instance, [3,12,14,21,36,48,49]).We claim the following (the proof, based on using Centre Manifold Theory [7,8,51], is given in Appendix A).
The epidemiological implication of the backward bifurcation phenomenon of the model ( 3) is that having R v < 1 is not sufficient (albeit necessary) to effectively control the spread of HPV and the associated dysplasia in the population.In other words, the combined use of mass vaccination and Pap cytology screening may fail to lead to effective control of HPV and the associated dysplasia even when R v < 1 (due to the phenomenon of backward bifurcation).In such a backward bifurcation scenario, effective disease control, when R v < 1, is dependent on the initial sizes of the populations of the model.Clearly, this (backward bifurcation) phenomenon makes effective control of HPV and the associated dysplasia difficult.It is worth stating that with the realistic parameter values in Table 2, the associated backward bifurcation parameter, a, is given by a = −1.14× 10 −5 < 0 (so that backward bifurcation cannot occur).In other words, this study shows that the phenomenon of backward bifurcation is not feasible in HPV dynamics using a realistic set of parameter values.
It is worth noting that for the case when the vaccine is perfect (i.e., v = 1), the coefficient a 0 of the quadratic ( 12) is zero, and b 0 > 0. Hence, the quadratic (12) has a unique root, given by λ * * m = −c 0 /b 0 , with c 0 ≤ 0. The case c 0 = 0 implies λ * * m = 0, which corresponds to the DFE.On the other hand, c 0 < 0 corresponds to Case 1 of Theorem 3.3, which shows the existence of a unique endemic equilibrium, only when R v > 1.
Corollary 1.The model (3) with perfect vaccine ( v = 1) does not undergo a backward bifurcation at R v = 1.
Thus, this study shows that the imperfect nature of the HPV vaccine (0 < v < 1) causes the phenomenon of backward bifurcation in the model (3).

3.4.
Global stability of equilibria: Special case.Consider the model (3) with the associated disease-induced and cancer-induced mortality set to zero (i.e., δ := δ f = δ c = δ c d = 0).When δ = 0, it follows from ( 4) and ( 5) that It is convenient to write Thus the system (3) can now be re-written as: where, g1 = g 1 | δ=0 = σ f + µ f .Since the variable W does not feature in any of the other equations of the model (13), it is removed from the analysis.The invariant region of the reduced system (13) and the associated DFE is given by Ẽ0 = (S * f , V * , 0, 0, 0, S * m , 0, 0), where, S * f , V * and S * m are given by (6).Furthermore, the reproduction number of the reduced model ( 13) is It is worth noting that the associated backward bifurcation parameter, a, for the reduced model ( 13) is given by so that the reduced model (13) will not undergo backward bifurcation (by Theorem 4.1 of [8]).To further confirm the absence of backward bifurcation in the reduced model ( 13), a global-asymptotic stability proof is given below for the DFE of the model (13).
Proof.Consider the non-linear Lyapunov function with the Lyapunov derivative given by (where the prime denotes differentiation with respect to t) In the above calculations, the following DFE relations (obtained from ( 13) at E 0 ) have been used: Since the arithmetic mean exceeds the geometric mean (that is, Therefore, F ≤ 0 if Rv ≤ 1.The equality F = 0 holds only (a) at the DFE Ẽ0 or (b) when Rv = 1 and is the singleton { Ẽ0 }.It follows from the LaSalle's Invariance Principle [34] that every solution of the system (13) with initial conditions in D converges to the DFE Ẽ0 as t → ∞.That is, (I f (t), P (t), I m (t)) → (0, 0, 0) as t → ∞.Substituting I f = P = I m = 0 in the equations for the susceptible females and males of the model (13) gives S f (t) → S * f and S m (t) → S * m as t → ∞.Similarly, it is easy to see that (when (I f , P, I m ) = (0, 0, 0)) the variable W (t) → 0 as t → ∞.Thus, (S f (t), V (t), I f (t), P (t), W (t), I m (t)) → Ẽ0 , as t → ∞ for Rv < 1.Hence, the DFE Ẽ0 , of the model ( 13), is GAS in D if Rv < 1.
The above analyses clearly show that the backward bifurcation property of the model ( 3) can be removed if the associated disease-induced and cancer-induced mortality is set to zero.The epidemiological implication of Theorem 3.5 is that for this special case (of the model ( 3) with δ = 0), the vaccine-targeted HPV types and associated dysplasia will be eliminated from the community if the combined use of (cohort) vaccination and Pap screening can make Rv < 1 (i.e., Rv < 1 is necessary and sufficient for the elimination of the relevant HPV types).Figure 3 depicts solution profiles of the model ( 3) with δ = 0, for various initial conditions, showing convergence to the DFE ( Ẽ0 ) for Rv < 1 (in line with Theorem 3.5).
The following result can be shown using the approach in Section 3.3.
Theorem 3.6.The reduced model ( 13) has an EEP of the form whenever Rv > 1, where

3.4.2.
Global stability of EEP.In this section, the global stability of the EEP of the reduced model ( 13) is established.Define: Theorem 3.7.The EEP, Ẽ1 , of the reduced model ( 13), is GAS in D \D 0 whenever Rv > 1.
Proof.Let Rv > 1, so that the unique EEP of the model ( 13), Ẽ1 , exists (Theorem 3.6).Here, too, the variable W of the model ( 13) is not included in the analysis.Consider the non-linear Lyapunov function The above can be simplified using the following endemic equilibrium conditions

This gives:
It follows, using Proposition 1 (Appendix B) in (15), that Since the arithmetic mean exceeds the geometric mean, the following inequalities hold: Hence F ≤ 0, with equality if and only if equality holds in each of the inequalities in (16).Thus, (S f , V, . Substituting these (endemic equilibrium) values in the equation for W in (13) , so that the EEP ( Ẽ1 ) of ( 13) is GAS in D \ D 0 .4. Assessment of vaccine and screening impact.In this section, the community-wide impact of mass vaccination and Pap cytology screening will be assessed using the model (3) with δ = 0 (or, equivalently, (13)).4.1.Vaccine impact.Let x = V * /N * f represents the fraction of females vaccinated at the disease-free steady state.It follows from (14) that where, is the basic reproduction number (the average number of new HPV cases generated by a typical infected individual in a completely susceptible population) associated with the model (3).Considering Rv as a function of x (i.e., Rv = Rv (x)), it can be shown that Since 0 < v < 1, it follows that ∂ Rv ∂x < 0 for 0 ≤ x < 1.Thus, Rv is a decreasing function of x.Since, in general, a reduction in the reproduction number signifies a reduction in disease burden (measured in terms of number of infections, HPVrelated morbidity, hospitalization, mortality etc.), the above analysis shows that an imperfect HPV vaccine will have a positive impact (in reducing HPV burden) for any x > 0 and v > 0 (that is, for a given vaccine efficacy, v > 0, vaccinating any fraction of susceptible females at steady-state will result in a decrease in HPV burden).Furthermore, there is a unique x c such that Rv (x) = 1, given by Lemma 4.1.The DFE of the reduced model ( 13), Ẽ0 , is GAS in Proof.Since Rv is a decreasing function of x, it follows that Rv < 1 whenever x > x c .Thus, it follows from Theorem 3.5 that the DFE of the model ( 13) is GAS in D in this case (with x > x c ).
The above result implies that if the fraction of vaccinated susceptible females exceeds the threshold level x c , then the DFE of the model ( 3) with δ = 0 is GAS (and the vaccine-targeted types will be eliminated from the community).A contour plot of Rv as a function of vaccine efficacy and coverage rate, depicted in Figure 4, shows that with the assumed vaccine efficacy of 90%, the relevant HPV types can be eliminated from the community if at least 88% of the susceptible females are vaccinated at steady-state (it should be recalled from Theorem 3.5 that the disease dies out when Rv < 1).
It should be noted that the above analyses also hold for the case where only (cohort) vaccination is used (i.e., for the case where Pap screening is not also used in the community).In other words, the singular use of a vaccination program can lead to the elimination of the relevant (vaccine-targeted) HPV types from the community if the vaccine efficacy and coverage rate are high enough.Figure 5 shows that the prevalence of the vaccine-targeted HPV types can be reduced by about 97.14% in 50 years.4.2.Screening impact.To assess the singular impact of Pap cytology screening, the model ( 3) is simulated in the absence of vaccination (φ = 0) and no diseaseinduced and cancer-induced mortality (δ = 0).Thus, this case models the situation where Pap cytology screening is the only intervention strategy adopted.Figure 6 depicts the cumulative number of females detected with cervical cancer or CIN (in the absence of vaccination) using various screening intervals (1, 2, and 3 years).It is evident from Figure 6 that the use of annual screening results in the detection of about 534 cases of females with cancer or CIN over 10 years.Furthermore, this figure shows that extending the screening interval from 2 to 3 years has marginal impact on the cumulative number of cases detected (while 2-year screening detects about 320 cases over 10 years, 3-year screening leads to the detection of about 228 cases over the same time period).  2 with δ = 0.It is shown that, unlike in the case of the vaccine-only control strategy (which reduces the prevalence of the vaccinetargeted HPV types and dysplasia by 97.14% in 50 years) (Figure 5), the combined vaccination-screening strategy (with 2-year or 3-year screening interval) achieved the same reduction in about 35 years (Figure 7).It is worth noting from Figure 7 that the 2-year and 3-year screening strategies resulted in essentially the same reduction (over the same time period).4.4.Model sensitivity to the proportion screened.As noted, a wide range of values of the proportion screened (χ) are cited in the literature (7% to 44%).Therefore the sensitivity of the model to this parameter is analyzed in this section.Under the base-case Pap test assumptions (with χ = 32%), the model predicts a total of 534, 320 and 228 cases in 10 years with the annual, 2-year and 3-year screening (Figure 6).Low proportion (7%) screened results in 157, 82 and 56 cases, respectively (Figure 6 top inset) and the high proportion (44%) screened results in 650, 410 and 298 cases detected, respectively (Figure 6 bottom inset).Whereas χ is an influential parameter in determining the number of detected cumulative cervical cancer and CIN cases, a comparison of the cumulative number of cases under annual and the 2-year screening respectively, relative to the 3-year screening reveals that the model is not much sensitive to the proportion screened; while the base-case value of the proportion screened (χ) results in the detection of 2.4 and 1.4 times more cases respectively relative to the 3-year screening, low proportion screened results in the detection of 2.8 and 1.5 times more cases with annual and 2-year screening respectively and the high proportion screened results in the detection of 2.2 and 1.4 times more cases.The comparison also signifies that (i) 3-year (rather than 2-year) screening is reasonable, and (ii) with low proportions screened, a higher screening frequency is more desirable.
Comparing the reduction in the prevalence of HPV infection and associated dysplasia leads to similar conclusions.The base-case value of the proportion screened results in 97.14% reduction in 35 years (Figure 7).The same reduction takes 41 years with a low proportion screened (Figure 7 top inset) and 33 years with a high proportion screened (Figure 7 bottom inset).As under the base-case scenario, the annual, 2-year and 3-year screening with low or high proportions screened makes no significant difference in providing the same amount of reduction in prevalence (hence the solid, dashed and dotted lines are indistinguishable in the inset figures).
Conclusions.A two-sex deterministic model is designed and used to study the transmission dynamics of HPV in the presence of Pap cytology screening and mass vaccination against some vaccine-targeted HPV types.The main theoretical results shown are summarized below: (i) The model exhibits the phenomenon of backward bifurcation where the stable disease-free equilibrium coexists with two endemic equilibria (one of which is suggested to be stable, by numerical simulations) when the associated reproduction number (R v ) is less than unity.The epidemiological implication of this result is that having the reproduction number less than unity, while necessary, is no longer sufficient for eliminating the HPV-infection and associated dysplasia from the community; (ii) The backward bifurcation phenomenon of the model is caused by any of the following mechanisms: (a) imperfect nature of the vaccine in preventing infection (0 < v < 1) (b) HPV-and cancer-induced mortality in the female population (δ f > 0, δ c > 0 or δ c d > 0); (iii) For the case when the HPV-induced and cancer-induced mortality is negligible (δ = 0), it is shown that the disease-free (HPV-infection-and dysplasia-free) equilibrium of the model is globally-asymptotically stable when the associated reproduction number ( Rv ) is less than unity; (iv) The model has a unique endemic equilibrium when the reproduction number (R v ) exceeds unity.This equilibrium is shown to be globally-asymptotically stable for the special case when the HPV-induced and cancer-induced mortality is negligible.Numerical simulations of the model suggest the following: (v) Pap screening every three years is reasonable.This result is consistent with the empirical findings [45,46]; (vi) The use of mass vaccination alone (with a vaccine efficacy of 90%) can lead to the elimination of the vaccine-targeted HPV types and associated dysplasia from the community if at least 88% of the susceptible female population is vaccinated at steady-state; (vii) Mass vaccination alone (without Pap screening) could result in 97.14% reduction in the prevalence of the vaccine-targeted HPV types and the associated dysplasia in about 50 years; (viii) If Pap screening (2-year or 3-year) is implemented in addition to mass vaccination, the time required to achieve the same amount of reduction in the prevalence of HPV infection and the associated dysplasia is reduced to about 35 years.
Appendix A: Proof of Theorem 3.4.
Proof.It is convenient to use the change of variables: denote the vector field of (3) in the notation (19), so that the model ( 3), with the conservation law (1), is re-written in the form: The Jacobian of the system (20), at the DFE E 0 , is given by It can be shown, from the Jacobian J(E 0 ), that (as in 7)), Consider the case when R v = 1.Suppose, further, that β f is chosen as a bifurcation parameter.Solving for β f from R v = 1 in (21) gives It is easy to verify that the transformed system (20), with β f = β * has a simple eigenvalue with zero real part, and all other eigenvalues have negative real parts.Hence, the centre manifold theory can be used to analyse the dynamics of (20) near β f = β * .In particular, Theorem 4.1 in [8] will be used.The application of this theorem entails the following computations.
Eigenvectors of J(E 0 )| β f =β * .The Jacobian of (20) at β f = β * , denoted by J β * , has a right eigenvector (corresponding to the zero eigenvalue) given by w = [w 1 , • • • , w 17 ] T , where where, g 3 g 4 − κ q ζ 1 > 0, and g 6 g 7 − κ d ζ 1 > 0. Further, the Jacobian J β * has a left eigenvector (associated with the zero eigenvalue) given by v where Computations of bifurcation coefficients a and b: For the system (20), the associated non-zero partial derivatives of f (at the DFE E 0 ) are given by The bifurcation coefficient b is similarly given by b Thus, it follows from Theorem 4.1 of [8] that the model (3) undergoes backward bifurcation whenever the bifurcation coefficient, a, given in (23), is positive.The inequality a > 0 can be expressed in terms of the model parameters as below.( S * * which is a contradiction.Hence, − ( S * *   3) showing the number of infected individuals (males and females) as a function of time, using various initial conditions.Parameter values used are as given in Table 2, with φ = 0.3 and c f = 4 (so that, R v = 3.7897 > 1).  2 with δ = 0, φ = 0.9 and c f = 1 (so that, Rv = 0.4833 < 1).2, with δ = 0.     2, with δ = φ = 0.The inset plots show the same corresponding to low screening rate due to low proportion screened (χ = 7%) and high screening rate due to high proportion screened (χ = 44%).

3. 1 .
Basic properties.The basic dynamical features of the model (3), subject to the conservation law (1), will now be explored.Lemma 3.1.The closed set

4. 3 .
Combined impact of vaccination and screening.The combined impact of Pap screening and mass vaccination is assessed by simulating the model (3) using the parameter values in Table

Figure 2 .
Figure 2. Simulations of the model (3) showing the number of infected individuals (males and females) as a function of time, using various initial conditions.Parameter values used are as given in Table2, with φ = 0.3 and c f = 4 (so that, R v = 3.7897 > 1).

Figure 3 .
Figure3.Simulations of the model(3) showing the number of infected individuals (males and females) as a function of time, using various initial conditions.Parameter values used are as given in Table2with δ = 0, φ = 0.9 and c f = 1 (so that, Rv = 0.4833 < 1).

Figure 4 .
Figure 4. Simulation of the model (3) showing the contour plot of Rv as a function of vaccine efficacy ( v ) and the fraction of susceptible females vaccinated at steady-state (V * /N * f ).Parameter values used are as in Table2, with δ = 0.
Years since the start of vaccination Percent reduction in prevalence of HPV and associated dysplasia

Figure 5 .
Figure 5. Simulation of the model (3) showing the percentage reduction in prevalence of the HPV infection and the associated dysplasia as a function of time, in the absence of Pap screening.Parameter values used are the same as those used to generate Figure 4.

Figure 6 .
Figure 6.Cumulative number of detected cervical cancer and CIN cases as a function of time.The dotted line represents 3-year screening, solid line represents 2-year screening, and dashed line represents annual screening.Parameter values used are as in Table2, with δ = φ = 0.The inset plots show the same corresponding to low screening rate due to low proportion screened (χ = 7%) and high screening rate due to high proportion screened (χ = 44%).

Figure 7 .
Figure 7. Simulation of the model (3) with δ = 0 showing the percentage reduction in prevalence of the HPV infection and the associated dysplasia as a function of time, with Pap screening.The solid curve represents annual screening, dashed line represents 2year screening and the dotted line represents 3-year screening.Parameter values used are the same as those used to generate Figure 4.The inset plots show the same corresponding to low screening rate due to low proportion screened (χ = 7%) and high screening rate due to high proportion screened (χ = 44%).
3 Population of infected females in Grade i of detected CIN 32T.MALIK, J. REIMER, A. GUMEL, E. ELBASHA AND S. MAHMUD