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 [ 35 ] . 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 ] . 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 . in 2007 , the canadian federal government provided the provinces and territories with $ 300 million to spend over 3 years for hpv immunization programs . by 2009 , every province and territory in canada had an hpv vaccination program in place . 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 ( 4986% within the first two years of the program ) . 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 . this is measured through the morbidity and mortality rates of hpv . 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 . 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 . 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 . 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 . 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 . 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 : ( 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 ? this paper is organized as follows . 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 . 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 ( c 
 iu , where 4  i  10 ) and vaccinated female children ( c 
 iv , where 4  i  10 ) . in the adult portion of the model , women are broken up into unvaccinated uninfected susceptible women ( a 
 u ) , vaccinated uninfected susceptible women ( a 
 v ) , unvaccinated and infected women ( i 
 u ) , or vaccinated infected women ( i 
 v ) . girls in grade 4 ( approximately 9 years old ) are described as ( 2)dc4dt=w1+cc4.for girls in grade i , where 5  i  10 , we have ( 3)dc(i+1)udt=1piciu1+cci+1u , dc(i+1)vdt=piciu+civ1+cc(i+1)v.uninfected adult women are described as ( 4)daudt=(1u)c10u+uiuf(wpw)au  waunaau+av , davdt=(1v)c10v+viv+f(wpw)au  (1)wavnaavav.infected adult women are described as ( 5)diudt=uc10u+waunuiuaiu+iv , divdt=vc10v+(1)wavnvivaiviv.uninfected men are described as ( 6)dmdt=mmiummivm+mnam.infected men are described as ( 7)dndt=mium+mivmmnan . a girl enters the model unvaccinated in grade 4 at a constant rate  
 w. at some grade between 4 and 10 , a proportion of the girls become vaccinated at rate p 
 i ( where p 
 i is the proportion of girls given the vaccine in grade i and  represents the vaccine efficacy for girls ) . when the girls enter grade 11 , a proportion of unvaccinated girls ( 1 u ) grow up to become uninfected unvaccinated women ( a 
 u ) , while the remaining unvaccinated girls (  
 u ) are categorized as infected unvaccinated women ( i 
 u ) to account for early childhood sexual activity . a proportion of vaccinated girls ( 1   
 v ) grow up to become uninfected vaccinated women ( a 
 v ) , while the rest of the vaccinated girls (  
 v ) are categorized as infected vaccinated women ( i 
 v ) 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 ) unvaccinated adult women ( a 
 u ) can become vaccinated at rate f( 
 w 
 p 
 w ) , where  
 w is the efficacy of the vaccine for adult women and p 
 w is the proportion of women who get the vaccine . the vaccine wanes at rate . unvaccinated adult women become infected ( i 
 u ) at rate  
 w when coming into contact with an infected man ( n ) . vaccinated adult women ( a 
 v ) can also become infected ( i 
 v ) at rate ( 1  ) 
 w , where  represents the ability of the vaccine to protect against all targeted types . men enter the model as susceptible ( m ) through a constant rate  
 m and stay in the model for approximately 10 years . unvaccinated susceptible men ( m ) can become infected after sexual activity with an infected woman ( i 
 u or i 
 v ) with transmission rate  
 m. infection clears at rate  
 u for unvaccinated women , at rate  
 v for vaccinated women , and  
 m for men . each compartment also includes a natural leaving rate ,  
 c or  
 a , depending on child or adult status . parameter descriptions , ranges , and sample values can be found in table 2 . for this model , it is assumed that hpv is only heterosexually transmitted . 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 . 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 . 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 . 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   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 ] . the disease - free equilibrium is(8)c4,c5u,c5v,c6u,c6v,c7u,c7v,c8u,c8v,c9u,c9v ,  c10u,c10v,au,av,iu,iv,m,n,where ( 9)c4=w1+c ,  and  we  define  c4u=c4 ,  c4v=0.then , for 5  i  10 , we have ( 10)ciu=1pi1c(i1)u1+c , civ=p(i1)c(i1)u+c(i1)v1+c , au=1uc10ufw  pw+a , av=fw  pwau+1vc10va , iu=0 ,  iv=0,m=ma ,  n=0 . 
  the jacobian matrix for this model evaluated at the disease - free equilibrium is j 
 dfe = [ j 
 dfe | j 
 dfe 
  | j 
 dfe 
 ] , where(11)the characteristic polynomial of the jacobian is ( 12)detji=1c13adethdetl , where ( 13)h=fafa,l=uawau0va(1)wavmmmmma. 
  solving det(h ) = 0 , we get(14)2+f+2a++af+a+=0,which has only eigenvalues with negative real part . solving det(l ) = 0 , we get ( 15)3+2++=0,where ( noting that m= ) ( 16)=3a+u+v+m+,=3a2+uv++um+v+m  + 2au+v+m+1wm  av  wmau,=a3+a2u+v+m+  + auv++um+v+m  + u(v+)m  1wm  ava+u+  wmaua+v+. 
  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 ) 
  > 0;(2 ) 
  > 0;(3 ) 
  > 0;(4 ) 
    > . 
  the first condition is satisfied , since all of the parameters are positive . the second condition is our threshold . if  < 0 , then we are guaranteed to have at least one positive real root of the characteristic polynomial . if  = 0 , 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  = 0 =  > 0 ( where  is a small positive perturbation ) . we will show that  > 0 and  > 0 in order to satisfy the routh  hurwitz criterion . we find ( 18)(1)wmav  = 1a+u+a3+a2+a  = 1a+u+wmaua+v+  = 1a+u++  uv+ma3 . 
  substituting this into  , we have ( 19)=3a2 + 2a+  1a+u+  a3+a2+  awmaua+v+  +  u(v+)mwmauwmau. 
  to require  > 0 , we need ( 20)(a+u+)3a2 + 2a+    + wmaua+v++  > a3+a2+  a+wmaua+u+  + uv+m,3a3 + 2a2+  a+3a2u+2au+  u+3a2  + 2a+  +(a+v+)wmau+  > a3+a2+  a+(a+u+)wmau  + uv+m,2a3+a2(+3u+3)+2au+  uu(v+)m  + 2a+  +(vu)wmau+>0.note that ( 21)uuv+m  = uuv++um+v+m  uvmum>0.it follows that  > 0 if ( 22)vu , 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  = 0 , the routh hurwitz conditions are satisfied , all roots have a negative real part , and our system is stable at the disease - free equilibrium . since  > 0 and  > 0 when  = 0 , then  = 0 is our threshold of stability . we used the product  =  
 w 
  
 m to numerically examine the routh and .   figure 2(a ) shows the scenario before vaccination while figure 2(b ) shows the scenario with 100% children and 100% adult vaccination . the region of stability does not significantly change depending on the efficacies (  
 c and  
 w ) or the probability of protection (  ) . in figure 2 , is represented by the solid red line while  is represented by the dashed blue line . when  > 0 ,  > 0 . the disease - free equilibrium is stable only if  > 0 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 . if  <  
 , then  > 0 and  > . this shows that we can use  as a ( local ) threshold of stability the basic reproductive number r 
 0 is a threshold that represents the average number of secondary infections caused by one infectious person in a completely susceptible population . from the analysis in section 3.1 , the r 
 0 threshold is ( 23)r0=wm((1)(a+u+)av+(a+v+)au)[a3+a2+  a+u(v+)],where au and av 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 r 
 0 , which in turn will give us insight into the most effective intervention strategies ( section 5.5 ) . there are critical vaccine efficacies for both children (  
 ) and women (  
 w 
 ) and there is a critical probability of protection (  ) that serves as a threshold where , even with 100% vaccination , the targeted types of hpv can not be eradicated in the population . these values are determined by first setting r 
 0 = 1 and rearranging for the desired parameter . in order to simplify the expression for r 
 0 , we rewrite the equilibrium values for our population in a general form . note here that k 
  represents the grade of childhood vaccination . for 4  k  10 , we find ( 24)cku=w1+ck3  for  kk,cku=w1pk11+ck3  for  k > k,ckv=0  for  kk,ckv=w1+ck3  for  k > k,au=wfpww+a1c7,av=wffpww+a1c7 . 
  to determine the expression for  
 , we look at only childhood vaccination ( i.e. , no adult vaccination ) , which we get by setting p 
 w = 0 . since childhood vaccination only occurs during one year , then(25)pk=1,when  k = k,0,when  kk. 
  rearranging the expression r 
 0 = 1 and solving for  
 , we find(26)=1+c7a2wm(a+u+)(1u)wa noticing that   0 ( i.e. ,  is not constant ) , we construct the following inequality from the definition of . we know that(27)wamc,w=wcwa=g , where w is the size of the female population when there are only women and g is the size of the female population when there are only girls . <  
 ( w ) , then 100% coverage rate of girls will not be sufficient for eradication , since r 
 0 > 1 . if  >  
 ( g ) , 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 can not 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 find(29)w=(1+)wmw(a+u+)(1u)a  ad = ld(ca)wmw(a+u+)=lllll1u1ca1,where d = ( 1 +  
 c) 
 a( 
 a 
  +  
 a 
 +   
 a 
  +  
 u( 
 v +  ) ) . w 
  is bounded as in ( 30)www()wg.the interpretation of  
 w 
  is similar to that of  
 . if the transmission in canada is below  
  = 1.02 , 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  
  = 2.95 however , if the transmission is above 2.95 , then even 100% vaccination can not eradicate hpv . these values allow us to estimate a range of likely transmission values , which are currently unknown . 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% ) . 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 30%   , giving a range of 0    0.3 . the recovery rate was found by determining the average infectious period for both high- and low - risk types of hpv ( 1/ ) and solving for  [ 26 , 27 ] . 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 ( 5096% ) and three ( 7188% ) 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 . 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 r 
 0 , while the box indicates the location of the upper and lower quartiles . 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 r 
 0 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 r 
 0 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 . comparing figures 3(a ) and 3(b ) , we notice that the r 
 0 values are generally smaller in the case of three doses , which is expected since the vaccine is less effective on a two - dose schedule . this can also be validated by figure 5 , which shows a low correlation between  or  
 w and r 
 0 . figure 4 shows the difference in the mean r 
 0 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 . since the true value of each parameter may fluctuate , we explore the sensitivity of r 
 0 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 r 
 0 , be it a positive or negative influence . here r 
 0 is most sensitive to the transmission of hpv from men to women ,  
 w , and from women to men ,  
 m. the next influential parameter is the probability of protection , . notice that the coverage rates ( p 
 i , where 4  i  10 , and p 
 w ) and vaccine efficacy (  and  
 w ) do not significantly affect the value of r 
 0 .   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 r 
 0 . 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 r 
 0 < 1 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 r 
 0 = 1 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 . 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 r 
 0 is always close to one when children are vaccinated , while the median value of r 
 0 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 r 
 0 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 r 
 0 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  
 w ) 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  
 w affects the value of r 
 0 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 . 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 . 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 ] . the model does not differentiate in hpv 6 , 11 , 16 , or 18 or include any other hpv types . we have also not ruled out the possibility of a backward bifurcation , meaning the disease may persist even for r 
 0 < 1 , 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 . note , however , that we have previously shown that this cohort assumption is negligible . 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 . 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 .