#### Abstract

Human papillomavirus (HPV) infection is the most common sexually transmitted infection, which is linked to several cancers and genital warts. Depending on the Canadian province, the quadrivalent vaccine is given to girls in grades 4 through 10 with either a two- or three-dose schedule. We use a mathematical model to address the following research questions: (1) Does the grade at which the girls are vaccinated significantly affect the outcome of the program? (2) What coverage rate must the provinces reach in order to reduce the impact of HPV on the Canadian population? (3) What are the implications of vaccinating with two versus three doses? The model suggests the grade of vaccination and the number of doses do not make a significant difference to the outcome of the public vaccination program. The most significant factor is the coverage rate of children and adults. We recommend that provinces vaccinate as early as possible to avoid vaccine failure due to previous infection. We also recommend that the main focus of the program should be on obtaining a large enough coverage rate for children and/or adults in order to achieve the desired outcome with either two or three doses of the vaccine.

#### 1. Introduction

Human papillomavirus (HPV) is the most common sexually transmitted virus, infecting 75% of the sexually active Canadian population [1, 2]. It can cause cancer of the cervix, vulva, oropharynx, penis, and anus, as well as genital warts [3–5]. Until recently, there were no medical interventions against HPV infections except monitoring the virus’s progress through the Papanicolaou test (or “Pap test”) and treating the ailments that developed [6, 7].

Since 2006, two vaccines have been approved for use in Canada against HPV. One (Cervarix) protects against types 16 and 18 while the other (Gardasil) also protects against HPV types 6 and 11. These vaccines offer protection for at least eight years [8]. In 2007, the Canadian federal government provided the provinces and territories with $300 million to spend over 3 years for HPV immunization programs [9]. By 2009, every province and territory in Canada had an HPV vaccination program in place [10]. The creation and implementation of the programs are the responsibility of the individual provinces and territories. This results in eight distinct strategies throughout Canada, which can be seen in Table 1.

These strategies vary in the number of doses of the quadrivalent vaccine given (only two doses in British Columbia and Quebec, with three doses administered elsewhere), the grade of the girls who are given the vaccine (ranging from 4 to 10), and the coverage rate (49–86% within the first two years of the program) [10]. Although the implementation of these programs differs between provinces, they all share a common goal of reducing the impact of HPV on the Canadian population [11]. This is measured through the morbidity and mortality rates of HPV [11].

Since the vaccine’s introduction, there have been a number of mathematical models assessing its impact on several populations. In terms of the cost-effectiveness of the vaccine, Brown and White used an optimal control model for vaccinating adolescent females and males in the United Kingdom [12]. Brisson et al. used a cohort model to estimate how many girls need to be vaccinated to prevent cervical cancer and genital warts in Canada [13]. In terms of the disease impact on a population, Barnabas et al. developed a transmission model to measure the impact of vaccinating against HPV-16 in Finland for both women and men [14]. Although it does consider the effects caused by HPV later in life, this model neglects to consider herd immunity. There are current clinical trials in British Columbia looking at the effect of taking two or three doses of the quadrivalent vaccine [15]. This will provide valuable information about the levels of protection provided by two or three doses. Llamazares and Smith? developed an epidemiological model that includes a widespread childhood vaccination program with voluntary adult vaccination [4]. The model is used as a preliminary investigation to determine whether provincial healthcare programs in Canada should fund adult HPV vaccination. The authors determined a threshold of eradication for the targeted types of HPV and critical efficacy and immunogenicity levels such that eradication is not possible at any coverage rate.

Here we develop a mathematical model to explore the current vaccination strategies across Canada as well as potential alternative strategies. These strategies are defined by the number of doses given, the grade of the girls the vaccine is given to, and the coverage rate achieved. Our model will provide similar insight into the effectiveness of two and three doses of the quadrivalent vaccine. This model also provides insight into what the programs should look like in order to achieve a desired outcome. This model provides a unique perspective in its evaluation of the provincial programs based on the grade of vaccination, number of doses given, and the necessary coverage rate. We address the following research questions: Does the grade at which the girls are vaccinated significantly affect the outcome of the program? What coverage rate must the provinces reach in order to reduce the impact of HPV on the Canadian population? What are the implications of vaccinating with two versus three doses?

This paper is organized as follows. Section 2 gives the model description and equations. Section 3 details the analysis of the model. Section 4 gives the expressions for critical thresholds of efficacy and probability of protection. Section 5 explores how the parameters were estimated, the results on the infection when varying the grade, dosage and coverage rate, and the sensitivity analysis of the model. Section 6 discusses the implications.

#### 2. The Model

Our model is composed of 19 equations, 13 that describe the childhood vaccination strategy and six that describe the disease propagation through adults. The female children are broken up into seven grades (4 through 10). Within each appropriate grade class, there are unvaccinated (, where ) and vaccinated female children (, where ). In the adult portion of the model, women are broken up into unvaccinated uninfected susceptible women (), vaccinated uninfected susceptible women (), unvaccinated and infected women (), or vaccinated infected women (). Men are considered either uninfected susceptible () or infected ().

Define

Girls in grade 4 (approximately 9 years old) are described as For girls in grade , where , we have Uninfected adult women are described as Infected adult women are described as Uninfected men are described as Infected men are described as

A girl enters the model unvaccinated in grade 4 at a constant rate . At some grade between 4 and 10, a proportion of the girls become vaccinated at rate (where is the proportion of girls given the vaccine in grade and represents the vaccine efficacy for girls). When the girls enter grade 11, a proportion of unvaccinated girls () grow up to become uninfected unvaccinated women (), while the remaining unvaccinated girls () are categorized as infected unvaccinated women () to account for early childhood sexual activity [29]. A proportion of vaccinated girls () grow up to become uninfected vaccinated women (), while the rest of the vaccinated girls () are categorized as infected vaccinated women () to account for a proportion of girls who have contracted at least one of the targeted types before receiving the vaccine. Once in the adult stage (grade 11 to around the age of 26), all of the disease propagation takes place. Once the girl reaches grade 11, she is assumed to be sexually active.

Unvaccinated adult women () can become vaccinated at rate , where is the efficacy of the vaccine for adult women and is the proportion of women who get the vaccine. The vaccine wanes at rate . Unvaccinated adult women become infected () at rate when coming into contact with an infected man (). Vaccinated adult women () can also become infected () at rate , where represents the ability of the vaccine to protect against all targeted types. Men enter the model as susceptible () through a constant rate and stay in the model for approximately 10 years. Unvaccinated susceptible men () can become infected after sexual activity with an infected woman ( or ) with transmission rate . Infection clears at rate for unvaccinated women, at rate for vaccinated women, and for men. Each compartment also includes a natural leaving rate, or , depending on child or adult status. The flow diagram is illustrated in Figure 1. Parameter descriptions, ranges, and sample values can be found in Table 2.

For this model, it is assumed that HPV is only heterosexually transmitted. Sexual partnerships are not explicitly modelled [30]. Female children are in grades 4 through 10, which are the years childhood vaccination can take place. It is assumed that vaccination only occurs in one year, that male vaccination is negligible, and that the proportion of female children vaccinated during other years is negligible. The rate at which the HPV infection is cleared is independent of previous infection status. The transmission from women to men is higher than the transmission from men to women [2]. Both women and men are active in the adult model for approximately 10 years because the vaccine is not recommended for women over the age of 26 [7]. We keep men in the model for the same length of time as women, on the grounds that, while there may be an age difference between young women and older men, such men are “aging out” with their sexual cohorts [4]. That is, they may choose new partners within this cohort (e.g., friends of existing partners), but we assume they do not revert to an even younger cohort after such a cohort ages out. (However, we have explored this assumption in more detail elsewhere [31] and shown that relaxing it has a negligible effect on the epidemic, so it is only included here for analytical convenience.)

The assumptions made about the vaccine are that the vaccine does not wane in children, that it may not protect 100% (this is based on the probability of protection and the efficacy), and that the vaccine does not protect someone who is already infected with the virus [32, 33]. The difference between a two- and three-dose schedule is the efficacy of the vaccine. Lastly, we do not consider disease-induced death since it does not play a role in removing sexually active individuals from the pool that we are considering (complications due to HPV, such as cervical cancer, are generally seen much later in life) [7, 34, 35].

#### 3. Analysis

##### 3.1. Stability of the Disease-Free Equilibrium

The disease-free equilibrium iswhere Then, for , we have

The Jacobian matrix for this model evaluated at the disease-free equilibrium is , whereThe characteristic polynomial of the Jacobian is where

Solving , we getwhich has only eigenvalues with negative real part.

Solving , we get where (noting that )

In order to determine the stability of the disease-free equilibrium in the linearized system, we use the Routh–Hurwitz stability criterion. For our cubic characteristic polynomial (15), we have four Routh–Hurwitz conditions that must all be satisfied in order for the disease-free equilibrium to be locally stable:(1);(2);(3);(4).

The first condition is satisfied, since all of the parameters are positive.

The second condition is our threshold. If , then we are guaranteed to have at least one positive real root of the characteristic polynomial. If , then we have a nonhyperbolic fixed point and must use stability theory to determine the stability of the disease-free equilibrium. We determine stability for the case when (where is a small positive perturbation). We will show that and in order to satisfy the Routh–Hurwitz criterion.

Define

Rearranging the original expression of , we find

Substituting this into , we have

To require , we need Note that It follows that if that is, if the recovery rate from infection is faster for vaccinated individuals (or at least not worse), which we expect to be the case.

When , the Routh–Hurwitz conditions are satisfied, all roots have a negative real part, and our system is stable at the disease-free equilibrium. Since and when , then is our threshold of stability.

We used the product to numerically examine the Routh–Hurwitz coefficients and . Figure 2(a) shows the scenario before vaccination while Figure 2(b) shows the scenario with 100% children and 100% adult vaccination. Note that the grade of childhood vaccination does not change the results. The region of stability does not significantly change depending on the efficacies ( and ) or the probability of protection (). In Figure 2, is represented by the solid red line while is represented by the dashed blue line. When , . Since , then . This satisfies the first three Routh–Hurwitz stability criteria.

**(a)**

**(b)**

The disease-free equilibrium is stable only if and . The value of indicated on the inset graph of Figure 2 shows the threshold of transmission between a stable and unstable disease-free equilibrium. These values change depending on the coverage rate. However, if , then and . This shows that we can use as a (local) threshold of stability. This threshold is discussed in greater detail in Section 3.2.

##### 3.2. Basic Reproductive Number

The basic reproductive number is a threshold that represents the average number of secondary infections caused by one infectious person in a completely susceptible population [36]. From the analysis in Section 3.1, the threshold is where and are the values at the disease-free equilibrium.

Knowing this threshold is significant, since it will tell us which parameters are involved in shifting the disease from persistence to eradication. This threshold will be used to determine which parameters have the greatest influence on , which in turn will give us insight into the most effective intervention strategies (Section 5.5).

#### 4. Critical Thresholds

There are critical vaccine efficacies for both children () and women () and there is a critical probability of protection () that serves as a threshold where, even with 100% vaccination, the targeted types of HPV cannot be eradicated in the population. These values are determined by first setting and rearranging for the desired parameter. In order to simplify the expression for , we rewrite the equilibrium values for our population in a general form. Note here that represents the grade of childhood vaccination.

For , we find

To determine the expression for , we look at only childhood vaccination (i.e., no adult vaccination), which we get by setting . Since childhood vaccination only occurs during one year, then

Rearranging the expression and solving for , we find

Noticing that (i.e., is not constant), we construct the following inequality from the definition of . We know thatwhere is the size of the female population when there are only women and is the size of the female population when there are only girls. This allows us to bound between If , then 100% coverage rate of girls will not be sufficient for eradication, since . If , then—with a sufficient coverage rate—it is possible for female-only vaccination to eradicate the targeted types of HPV. If is between the bounds, we cannot accurately predict the outcome of the vaccine intervention.

The critical adult efficacy occurs when there is no childhood vaccination but there is 100% adult vaccination. Using a similar approach to the critical childhood efficacy, we findwhere . is bounded as in The interpretation of is similar to that of .

#### 5. Numerical Simulations

##### 5.1. Transmission

Figure 2(a) shows the region of stability when no vaccination occurs. If the transmission in Canada is below , then the disease would be eradicated without any intervention. Since HPV is endemic throughout the world and has not become eradicated, we conclude that the actual transmission parameters must be greater than 1.02 (assuming all other parameters used in the simulation are correct).

Figure 2(b) shows the region of stability for 100% childhood and 100% adult vaccination. As long as the transmission is lower than , then the introduction of the vaccine has the ability to eradicate HPV. However, if the transmission is above 2.95, then even 100% vaccination cannot eradicate HPV. These values allow us to estimate a range of likely transmission values, which are currently unknown.

##### 5.2. Estimation of Parameters

The probability of protection, , was estimated by including the transmission of the targeted HPV types and an estimate for the proportion of girls who become sexually active before grade 11 (30%) [25]. Assuming all of the sexually active girls came into contact with someone able to transmit one of the targeted HPV types, the infection for this grade class would spread at a rate of , giving a range of .

The recovery rate was found by determining the average infectious period for both high- and low-risk types of HPV () and solving for [26, 27].

##### 5.3. Varying Grade

We looked at the significance of vaccinating girls at different ages by constructing a box plot as seen in Figure 3. For a specific grade, the coverage rate was varied between 0 and 100%. In Figures 3(a) and 3(b), the vaccine efficacy in adults was used to represent two (50–96%) and three (71–88%) doses, respectively. Latin Hypercube Sampling was used to sample parameter values from the given ranges in Table 2 using a uniform distribution. Latin Hypercube Sampling is a statistical sampling method in which parameters are assigned a range of values, and a distribution of plausible collections of these parameter values are created [37]. was calculated 1000 times (diamonds) per grade. Note that only one grade (or adult) is considered during each run; for example, a box plot representing grade 4 only takes into account vaccinating girls in grade 4 and does not include vaccinating any other grade or adults. The thick red horizontal lines represent the median value of , while the box indicates the location of the upper and lower quartiles.

**(a)**

**(b)**

Comparing the median values in Figures 3(a) and 3(b), we notice that when girls are between grades 4 and 10, the values are all around 1, whereas the median is always above one for adult vaccination regardless of the number of doses given. The upper quartile values and lower quartile values also follow this trend, where the childhood vaccination values all have a similar range that is noticeably smaller than when adult vaccination occurs.

Figure 4 shows the mean values, rather than the median. While there is more variation between grades, vaccinating with 3 doses is always superior to vaccinating with 2 doses. Due to transmission rates (see Figure 5), neither form of vaccination will lead to eradication. However, childhood vaccination is always superior to adult vaccination.

##### 5.4. Varying Dosage

Comparing Figures 3(a) and 3(b), we notice that the values are generally smaller in the case of three doses, which is expected since the vaccine is less effective on a two-dose schedule. However the difference is minimal. This can also be validated by Figure 5, which shows a low correlation between or and .

Figure 4 shows the difference in the mean when using a two-dose versus a three-dose regime. Once again, having three doses is clearly superior to two, but the difference is not significant.

##### 5.5. Sensitivity to Variations

Since the true value of each parameter may fluctuate, we explore the sensitivity of to the parameter values in Table 2 using Latin Hypercube Sampling. This in turn allows us to use partial rank correlation coefficients to rank the parameters in terms of their influence on , be it a positive or negative influence. Figure 5 shows the partial rank correlation coefficient sensitivity analysis on . The parameters were varied for 1000 runs. Here is most sensitive to the transmission of HPV from men to women, , and from women to men, . The next influential parameter is the probability of protection, . Notice that the coverage rates (, where , and ) and vaccine efficacy ( and ) do not significantly affect the value of . Figure 6 shows the output of the Latin Hypercube Sampling for the three most influential parameters as well as some other parameters of interest: the recovery rates and the waning rate. Note the extremely weak correlation between the latter parameters and .

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

##### 5.6. Varying Coverage Rate

The coverage rate required to eliminate the targeted high-risk HPV types can be seen in Figure 7. Above the curves is the region where and the targeted types are theoretically eradicated, whereas below the curves is the region where the disease persists in the population. The curves indicate when for either two doses (lower red curves) or three doses (upper black curves) for different values of the waning rate. As expected, the curve for two doses is higher than three, since the efficacy for two doses is slightly lower. The targeted types can be eradicated if childhood vaccination is supplemented with significant adult vaccination. However, as the waning rate of the vaccine increases, the window for eradication shrinks, requiring a significant amount of both childhood and adult vaccination.

A childhood-only vaccination program can theoretically eradicate the targeted HPV types with either two or three doses, as long as the appropriate coverage rate is achieved (Figure 7). For example, if 80% of children are vaccinated with a 3-dose vaccine that does not wane (the lowest of the curves), then at least 40% of adults must be vaccinated to achieve eradication of targeted types. These requirements become harsher as the number of doses decreases or as the vaccine wanes. Note that these curves do not change with the grade of childhood vaccination.

#### 6. Discussion

Based on the results of the model, we conclude that, with sufficient childhood and adult vaccination, it is theoretically possible to eradicate targeted HPV types. We also determine that the grade of vaccination before sexual debut does not significantly affect the prevalence of the targeted HPV types for the Canadian population. Comparing Figures 3(a) and 3(b), we observe that the median value of is always close to one when children are vaccinated, while the median value of is above one when only adult women are vaccinated, independent of the number of doses given. This suggests that the targeted types are candidates for eradication if significant childhood vaccination can be achieved, whereas eradication is less likely from adult-only vaccination.

Looking only at childhood-only vaccination, there does not seem to be a large difference between the grades. From a mathematical standpoint, it is easy to understand why the grade of vaccination does not matter in the model, since there is no impact in terms of the disease dynamics, whether girls are vaccinated in grade 4 or 10. The large jump is due to the change in vaccine effectiveness for those previously exposed to the considered HPV types, as well as the possibility of overvaccinating women, since they remain within the vaccination cohort for 10 years, whereas children are only vaccinated within a single year. This is important to determine because it shows there are likely no underlying grade-dependent trends. However, from an epidemiological standpoint, we know that the vaccine has little to no impact on those individuals who have a previous infection with a targeted type [32, 33]. Therefore public vaccination strategies should choose the grade of vaccination based on epidemiological and vaccination program limitations.

Assuming that the number of doses only changes the efficacy of the vaccine, the evidence from this model suggests that the number of doses does not significantly change the outcome of the vaccine strategy. Figure 7 shows that the possibility of eradication of targeted HPV types exists with either two or three doses, independent of the grade the childhood vaccine is given in. This evidence follows the findings by several clinical studies that suggest that two doses of the bivalent vaccine may protect just as well as three doses [5, 15, 22]. The higher the coverage rate, the higher the success of the vaccine program (Figure 7). Based on this evidence, we suggest each province chooses the number of doses based on cost or ease of program implementation, as long as an appropriate coverage rate is matched to achieve the desired level of infection in the community.

The value of is mainly affected by the ability of the targeted HPV types to transmit and the ability of the vaccine to prevent the infection. Thus the most effective way to decrease is to decrease the transmission of the targeted types of HPV. This could be done through increasing condom use, reducing the number of sexual partners, or increasing the heterogeneity of the sexual contact network of the population [16, 38]. Since HPV has been infecting humans for millions of years and has not been eradicated, the real transmission rate in the Canadian population must be greater than 1.02. By studying the critical transmission value for 100% childhood and adult vaccination, we know that if the actual transmission rate is greater than 2.95, then, even with 100% vaccine coverage, the targeted HPV types could not be eradicated. Looking at the stability of the disease-free equilibrium based on the transmission is reasonable since the transmission parameters are the only two parameters that significantly affect the value of the basic reproductive number, as shown by the sensitivity analysis.

We chose the estimates for the efficacies ( and ) based on clinical evidence of both the bivalent and the quadrivalent vaccine. Obviously, a single vaccine program will only use either the bivalent or the quadrivalent vaccine. Since this model does not differentiate between different HPV types and the efficacies for either vaccine do not differ greatly, the numerical results hold for either vaccine. Through sensitivity analysis, using larger ranges of the efficacy than observed for either vaccine, we see in Figure 5 that it is not important to narrow down the range of efficacies since neither nor affects the value of significantly. In order to estimate the efficacy of a two-dose schedule of the quadrivalent vaccine, we used the efficacy of the bivalent vaccine from the Costa Rica clinical trial [22]. This model should be updated once better data is available for the efficacy and effectiveness of a two-dose quadrivalent HPV vaccine.

It is important to note that our model has several limitations. In terms of vaccination programs, the model does not include catch-up programs unless included in the initial population conditions. Our model does not allow the vaccine program to differ amongst individuals. Each person receiving the vaccine with either two or three doses must complete the regimen within one year. The possibility of male vaccination is not included, although it is now approved by Health Canada [1, 5, 11]. We limit transmission to that of heterosexual couples. The model does not differentiate in HPV 6, 11, 16, or 18 or include any other HPV types. This model does not include immigrating (infectious) adults. We have also not ruled out the possibility of a backward bifurcation, meaning the disease may persist even for , complicating eradication efforts.

From the base adult model, we assume that men who have sexual relations with women in the sexually active cohort (i.e., women who are eligible for adult vaccination) do not continue to find new partners in this age group as time goes on. The sexually active cohorts of men and women are linked only for the time (approximately 10 years) that adult women are sexually active and eligible for vaccination [4]. Note, however, that we have previously shown that this cohort assumption is negligible [31].

Since the provincial vaccination programs are already in place and each province has chosen their dosing schedule (although it is possible they may change), we recommend that all of the provinces focus on increasing their coverage rate for both the public vaccination of girls and private vaccination of women to at least 80%. This number is realistic for many provinces, since several have already reached this goal (Quebec, Nova Scotia, Newfoundland, and PEI). For the other provinces, increasing their coverage rate may be accomplished by creating public-education campaigns or may involve making protocols similar to vaccines necessary to attend school such as those for Hepatitis B [17, 39].

Our model suggests that the grade of vaccination does not affect the outcome of the vaccination program. Therefore we suggest provinces vaccinate girls as early as possible to avoid vaccine failure due to previous infection. We also recommend that the number of doses should be chosen for optimal uptake. The main focus should be on obtaining large enough coverage rates for children and/or adults in order to achieve the desired outcome: using the vaccine to reduce the prevalence of HPV types 6, 11, 16, and 18 in the population.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The authors thank Julien Arino, David Sankoff, and Lucy Campbell for technical discussions. Carley Rogers received a MITACS-funded internship at the Public Health Agency of Canada where the majority of this work was completed. Robert J. Smith? was funded by an Early Researcher Award from the Ontario Ministry of Research and Innovation. For citation purposes, note that the question mark in “Smith?” is part of his name.