Model for Assessing Human Papillomavirus Vaccination Strategies

A prophylactic quadrivalent vaccine can cost-effectively reduce the incidence of cervical cancer, cervical intraepithelial neoplasia, and genital warts.


Introduction
This technical report includes the mathematical model that was constructed to assess the impact of HPV vaccination strategies. It provides a detailed description of various model components as they relate to HPV infection, disease progression, vaccine characteristics, vaccination strategies, and the impact of HPV vaccination on epidemiologic and economic outcomes. The model allows for aggregating costs of vaccination, screening, and treatment of the population over time, compares them with total health outcomes as measured, for example, by quality adjusted life years (QALYs), and calculates incremental cost-e¤ectiveness ratios for various vaccination strategies. In constructing this model, we reviewed other relevant previous models and incorporated some of their structures and inputs. These included cervical cancer screening cohort models [17,18,21,58,61,50], HPV vaccination cohort models [76,71,53,31,32], and HPV vaccination dynamic models [39,28,5,19]. This model di¤ers from its predecessors in several ways. First, the approach is more comprehensive in the sense that it incorporates the epidemiology of HPV infection, disease, and economics into a single dynamic model. Besides capturing the direct and indirect 'herd immunity'bene…ts and costs of vaccination for the population over time, the added advantage of this latter approach is its transparency, making critical review of the model and reproducibility of results [81] feasible without needing to review the actual source code used to generate the results. In particular, publication of the model includes the mathematical equations that summarize in their entirety the actual workings of the model. These equations can then be entered into any standard mathematical software package such as Mathematica R (Wolfram Research, Champaign, IL) or MatLab R (MathWorks, Natick, MA) to reproduce the results. Second, we also convened an expert panel that reviewed model assumptions and provided guidance on some aspects of the natural history of disease where there was little or no clinical evidence. Finally, key inputs in this model are based on data from recent studies that were not available when previous models were constructed.
For ease of exposition, the model is divided into two major components. The …rst part, which is presented in section 2, is a description of the demographic aspects of the model. This component of the model is intended to mimic the current age structure of the US population. Section 3 includes the second part which consists of the epidemiologic model that describes HPV transmission, and progression to cervical intraepithelial neoplasia (CIN), cervical cancer, and genital warts. Because females who undergo hysterectomies for benign conditions are no longer at risk of developing CIN and cervical cancer but can contribute to the transmission of HPV, another submodule for benign hysterectomy is created. Descriptions of the forces of infection, mixing preferences, and estimates of the epidemiologic model completes section 3. In sections 4 and 5, we describe how the epidemiologic and economic impact of screening and vaccination strategies are assessed.

The demographic model 2.1 Demographic model structure
The demographic model is a modi…ed version of the initial-boundary-value problem for age-dependent population growth described in more details in [36]. The population is divided into n age groups de…ned by the age intervals [a i 1 ; a i ], where a 1 <a 2 < : : : <a n = 1 (all the symbols used to describe variables and parameters are de…ned in Table1 and 2). The number of individuals N i (t) at time t in the age interval [a i 1 ; a i ] is the integral of the age distribution function from a i 1 to a i . Assuming that the population distribution has reached a steady state with exponential growth or decay of the form e qt , Hethcote [36] derived a system of n ordinary di¤erential equations (ODEs) for the sizes of the n age groups. The simple demographic model used here divides the population into 2 gender (k = f; m) groups, and 17 age (i = 1; 2; : : : 17) groups (12-14, 15-17, 18-19, 20-24, 25-29, 30-34, 35-39,40-44, 45-49, 50-54, 55-59, 60-64, 65-69, 70-74, 75-79, 80-84, and over 85). This age grouping is chosen to accurately account for patterns of HPV transmission among sexually active groups, cervical cancer screening patterns, and risk of cervical cancer development among females, and genital wart occurrence among both males and females. Similar age groupings have been used by other sexually transmitted diseases models [24,25]. However, these models assumed an age of sexual debut of 15 years. By setting the age of sexual debut to 12 years, our model captures HPV transmission and disease that occurs before age 15. Recent data suggest age of …rst sexual intercourse is younger than 15 for some teenagers and adolescents. For example, according to data from the National Survey of Family Growth, 19% of female teenagers had had sex before age 15 in 1995, compared with 21% of male teenagers [1].
The sexually active population is further strati…ed into L sexual activity groups (l = 1; : : : ; L), de…ned according to the gender-, sexual activity-, and age-speci…c rate of sex partner change per unit time c kli . The number of sexual activity groups considered here is 3 (L = 3). New additions to the sexually active population enter gender k, sexual activity l, and cervical screening category b (b = 1; 2) at rate of B klb . Because males do not participate in cervical screening, throughout the model the subscript b does not apply to them. For example, B mlb = B ml . Individuals die of non-cervical cancer related causes at an ageand gender-speci…c per capita death rate ki per year and females with cervical cancer (categories CC s and DCC s ) also have an additional age-and stagedependent mortality rate si (s = L; R; D). It is assumed that being in any CIN or genital warts state does not pose an additional risk of death. Individuals are transferred between successive age groups at an age-and gender-speci…c per capita rate d ki per year given by [36] where band i is the number of years within age group i. The annual growth rate q of this demographic model should also satisfy a modi…ed age-group form of the Lotka characteristic equations [36] where % b denotes the fraction of females entering cervical screening category b, with % 1 + % 2 = 1.
After taking into account cervical cancer-induced mortality and replacing fertility rates in Hethcote's model [36] by recruitment rates into the sexually active population B klb , the demographic model is given by the following system of 102 (= 17 2 3) ODEs: i 2, s = L; R; D, where d k17 = 0. All variables, parameters, and subscripts are de…ned in Tables 1 and 2 and the text.

Estimates of the demographic model parameters
Death rates for males and females without cervical cancer are obtained from Vital Statistics data on gender-and age-speci…c mortality rates, all races, 2002 [51]. Cancer mortality data are obtained from Surveillance, Epidemiology, and End Results (SEER) Cancer Statistics age-speci…c mortality rates, 1997-2002 [75]. Because the U.S. population grew at a decennial rate of 13.2% between 1990 and 2000, the annual population growth rate was 1.23%. With recruitment rates into the sexually active population of 1.9% of the male active population and 1.7% of the female population, the largest annual growth rate q that satis…es the solution of the Lotka characteristic equation was 0.5%. Therefore, the annual growth rate q of this demographic model was set to zero, and B klb was chosen to satisfy the Lotka characteristic equation. This will also ensure that variation in the results across strategies is mainly due to epidemiologic and program features rather than peculiar characteristics of the demographic model [36]. The sensitivity of the results to this assumption will be tested using an annual population growth rate of 1.23%. The initial population size is set to 100,000, divided equally between males and females. With the proportion of adults in sexually activity class l given by ! l , the total number of individuals in sexual activity group l is given by    By using d ki 1 N kli 1 ( ki + d ki )N kli with the above equation, we obtain the initial number of individuals in the youngest age group (12-14 years) of each gender and sexual activity category as The initial numbers of other age groups are given by l = 1; 2; 3; i = 2; 3; : : : ; 17.
Note that the size of the male population in the model is always at a steadystate given by =2 = 50; 000. However, the size of the female population is not constant during the transient dynamics following vaccination because females are subject to additional cervical cancer-induced mortality.
The structure of the over 12-year old US population with 0% and 1.23% annual growth rates, together with data from the 2000 population census are plotted in Figures 1 and 2, respectively. The model …ts well for early age groups, underestimates around age 40, and overestimates the number of people over age 40 years. It should be noted that the current model does not capture special characteristics of the US population such as the "Baby Boom" and migration.

The epidemiologic model
The epidemiologic model can be thought of as comprising three components: HPV transmission, cervical cancer development, and genital warts occurrence.

HPV transmission
To simplify the analysis, only three ( types 16/18 = 1, types 6/11 = 2, and coinfection =12) HPV type groupings are modeled. The sexually active host population of size at time t is divided into distinct epidemiologic classes, depending on the host's susceptibility to infection or the host's status with respect to infection, disease, screening, and treatment. The HPV component consists of 17 epidemiologic classes (X; V; Y; W; U; P; Z; Q), with each class further strat-i…ed by gender (2 groups), age (17 groups), and sexual activity (3 groups). The female population has two additional strati…cations distinguishing females that are regularly screened from those who are never screened, and females who had hysterectomies from those with intact cervices. A schematical representation of the HPV transmission model is shown in Figures 3 and 4.

Susceptible individuals X
New additions to the sexually active population, at a rate of B klb , enter into the uninfected (susceptible) category of gender k, sexual activity group l, and screening category b. A fraction of them is vaccinated at rate kl0b and move to category V and the remaining fraction enter category X of susceptible individuals. The model also assumes that a proportion of individuals in other age groups and epidemiologic classes is vaccinated at rate klib and move into the vaccination classes V , W , P , or Q. It is assumed that the vaccine does not confer any therapeutic bene…ts to individuals already infected. Individuals in class X acquire HPV infection with type h at a gender, sexual activity group, age, and time dependent rate h kli ; where h = 1; 2; 12. In this notation, 1 kli denotes infection with types in group 1 (HPV 16/18) and 12 kli infection with types in both groups (HPV 16/18 and HPV 6/11). The number of people in category X klib is reduced by infection h kli , vaccination klib , benign hysterectomy ki , death from other causes ki , and aging d ki . The ODEs for category X are i = 2; : : : ; 17; l = 1; 2; 3; k = f; m; b = 1; 2.

Infected individuals Y
When transmission occurs, the unvaccinated X and vaccinated S susceptible individuals enter the Y class of infected individuals. Individuals enter class Y after they recover from genital warts at rate gk but are still infected with probability 1 h gki Females enter class Y if their CIN spontaneously regress at rate f s but are still infected with probability 1 h ki . Individuals leave this class and enter the Z class of recovered people with immunity when the infectious period for HPV ends. Unvaccinated infected individuals in the Y class resolve infection at an age-, gender-, and type-speci…c per capita rate of h ki . Individuals develop CIN and genital warts at rate h ks and h gk , respectively. The ODEs for category Y h are The ODEs for coinfection are given by

Partially immune individuals Z
Individuals enter class Z when recovered from CIN or genital warts and having resolved infection. It is assumed that immunity derived from natural infection can be temporary, and individuals in the Z category can eventually move to the susceptible class X at rate h zki . Individuals in the Z class who are susceptible to one type can be infected with that type and move to class U . The ODEs for category Z h are The ODEs for the fully immune individuals Z 12 are The vaccine provides incomplete protection against the high-risk and low-risk types at rates 1 ' 1 ki and 1 ' 2 ki , respectively. A vaccinated person moves into compartment W upon infection with any type. Upon clearance of infection at rate h ki faster than natural infection, the person moves to compartment Q. The vaccine-induced immunity wanes at rate ki .

Infected individuals with partial immunity U
The number of people in category U is reduced by vaccination klib , resolution of infection h ki , and onset of disease. The ODEs for category U are

Vaccinated individuals V
When 12-year olds are o¤ered the vaccine, a fraction of them kl0 are vaccinated and move into the vaccination class V . Also, individuals in class X are vaccinated at rate kl1b and enter category V . The vaccine-induced immunity of those in the vaccinated class V wanes, so that people eventually move to the susceptible class S at an age-and gender-dependent rate ki . It is assumed that when an individual loses vaccine-derived immunity, the individual becomes susceptible to infection with any of the types. Vaccinated individuals can also experience a break-through infection and enter the class W of infective people at per capita rate ' h k h kli . The ODEs for category V are 3.1.6 Vaccinated individuals with waned immunity S Individuals in this class can get infected at the same rate as those in the susceptible class X. The ODEs for class S are

Infectious vaccinated individuals W
Individuals infected with one type and susceptible to the other move category W when vaccinated. Vaccinated individuals are infected at an age-and gender-speci…c rate ' h k times slower, and recover from infection at a rate h ki faster than unvaccinated infected individuals and move to class Q. They also progress to disease at a di¤erent rate ( h wks or h gkw ) compared with that of infected unvaccinated individuals. The ODEs for category W are The ODEs for coinfection W 12 are Infected vaccinated individuals (category W ) recovering from infection and individuals with natural immunity to one type (category Z) receiving the vaccine move to category Q. Individuals in this class who are susceptible to one type can be infected with that type and move to class P . The ODEs for category Q are 3.1.9 Vaccinated, infected individuals with partial immunity P Coinfected vaccinated individuals recovering from one infection (category W 12 ), vaccinated individuals (category Q) getiing infected, and individuals infected with one type (category Z) receiving the vaccine move to category P . The ODEs for category P are Note that for males, mi = h wms = h ms = h ms = h cms = 0.

Cervical intraepithelial neoplasia
Infected females (whether vaccinated or not) can develop CIN and move to the CIN segment of the model. There are several states that represent the true histological health status of a female: infected with a normal cervix, CIN grade 1 (CIN 1), CIN grade 2 (CIN 2), and CIN grade 3 (CIN 3). Females in the CIN and cancer stages are further classi…ed into unknown, detected, or treated classes. There are also two additional absorbing states where only females who are no longer at risk of developing cervical cancer enter. These are benign hysterectomy for reasons other than cervical cancer (at an age-speci…c rate f i ) and treated and cured CIN at stage-speci…c rate (1 s ) s . Females in these two states are considered to be at no risk of developing cervical cancer [61]. However, females with hysterectomies for benign conditions can be infected and are at risk of developing genital warts [9]. Further, to take into account the fact that treatment of CIN does not completely eliminate the virus, another category of women with treated CIN who remain infected after treatment (ICIN ) was created. Females enter this category from the detected state at rate s s and stay there until their CIN recurs at rate h r s or they clear infection. An infected female with a normal cervix can only directly progress to CIN h s (at rate h f s if unvaccinated or h wf s if vaccinated), die due to causes other than cervical cancer, or remain infected without progressing to CIN ( Figure 5). The respective progression rates given coinfection are h cf s and h cwf s . For the base case, it is assumed that cases with coinfection progress to CIN according to the rate of high-risk HPV types. That is, 1 cf s = 1 f s , 2 cf s = 0, 1 cwf s = 1 wf s , and 2 cwf s = 0. It is assumed that infected females classi…ed as CIN can progress only to higher CIN states (CIN1 to CIN2, CIN2 to CIN3), or cancer (CIN3 to cervical carcinoma in situ, CIS) at rate h si , regress to normal at rate h s or CIN state g at rate h sg , die from other causes, be detected at rate sib and be treated and cured at rate s , or remain in that CIN state. Coinfection of females in CIN and cervical cancer states is not modeled. It is assumed that regression from CIN states does not necessarily imply recovery from HPV infection. A female whose CIN regresses to normal but is still infected moves to the infected category Y h f i at an age-and stage-speci…c rate h s (1 h f i ) regardless of her vaccination status. Only mutual regression from both HPV and CIN confers immunity against that type. Females regressing from CIN, whose HPV infection clears, move into class Z at an age-and state-speci…c rate h where s = 1; 2; 3.

Treated CIN T CIN s
It is assumed that treatment does not completely eliminate infection. A fraction of treated females s will remain infectious after treatment and move to the category treated but infectious ICIN s . Equations for treated CIN are where s = 1; 2; 3.

Cervical carcinoma in situ
It is assumed that females classi…ed as CIN can progress to carcinoma in situ (CIS). Because females spend, on average, a long time in CIS, two CIS states are modeled (CIS 1 and CIS 2). It is assumed that regression from CIS states is not possible. CIS is further divided into several epidemiologic classes (CIS s , DCIS s , T CIS s , ICIS s ; s = 1; 2), with each class further subdivided into age (= 17), sexual activity (= 3), and screening (= 2) groups.

Undetected CIS CIS s
The number of females with undetected CIS increases as they progress from CIN 3 (severe dysplasia) or fail treatment. Screening 3+sib and progression to higher disease grades h 3+si reduce the number of females in this category, s = 1; 2. Equations for undetected CIS are

Detected CIS DCIS s
Detection of CIS occurs only as result of screening at rate 3+sib . If it is not treated and cured at rate 3+s , CIS can progress to a higher grade h 3+si or cancer. Equations for detected CIS are where s = 1; 2.

Treated CIS T CIS s
It is assumed that treatment does not completely eliminate infection. A fraction of treated females 3+s will remain infectious after treatment and move to the category treated but infectious ICIS s . Equations for treated CIS are where s = 1; 2.

Treated CIS but infectious ICIS s
It is assumed that CIS recurs at rate h r3+s and women with recuring CIS move to category CIS s . Infection can also resolve and individuals enter category Z h . Equations for treated but infectious CIS are where s = 1; 2.

Cervical cancer
There are several states that represent the health status of a female with cervical cancer: localized cervical cancer (LCC), regional cervical cancer (RCC), distant cervical cancer (DCC), and cancer survivors who are free from cancer ( Figure  6). Females in cancer stages are further classi…ed into unknown, detected, or treated classes. A female with an invasive cancer can progress only to the next higher cancer state CC h s (LCC to RCC, RCC to DCC) at rate si (s = L; R), her cervical cancer is detected at rate sib and successfully treated and move to the cancer survivors state at rate s , die from cancer at rate si , or stay in that undetected cancer state. Regression from invasive cancer to normal is not allowed. It is assumed that females who were successfully treated for invasive cancer are no longer infectious.

Undetected cervical cancer CC
where i 2.

Detected cervical cancer DCC s
Detected cancer cases are treated and cured at rate s and move to the cancer survivors category SCC. Equations for detected CC are where s = L; R; D.

Cervical cancer survivors SCC
Equations for cancer survivors are

Genital warts GW
Individuals (whether vaccinated or not) infected with HPV 6/11 can develop genital warts at rate 2 gwk and move to the genital warts class GW . Of those, a proportion gs will remain asymptomatic and will not be treated whereas the rest will be recognized and treated. Individuals recovering from genital warts at rate 2 kg move to class Z. It is assumed that only infection with HPV 6/11 can cause genital warts whereas infection with HPV 16/18 does not lead to genital warts [84]. The asymptomatic genital warts class consists of the following ODEs The symptomatic genital warts class consists of the following ODEs

Susceptible individuals HX
The ODEs for category HX are

Infected individuals HY
The ODEs for category HY are The ODEs for HY 12 are

Partially immune individuals HZ
The ODEs for category HZ are The ODEs for HZ 12 are

Infected individuals with partial immunity HU
The ODEs for category HU are f li

Vaccinated individuals HV
The ODEs for category HV are

Vaccinated individuals with waned immunity HS
The ODEs for classes HS are

Infectious vaccinated individuals HW
The ODEs for category HW are The ODEs for HW 12 are 3.6.9 Vaccinated, infected individuals with partial immunity HP The ODEs for category HP are

Genital warts GW
The genital warts class consists of the following di¤erential equations

Forces of HPV infection
The rate at which susceptible individuals acquire infection with type h (per capita force of infection) h kli is gender, sexual activity, age, and time dependent. The rate h kli at which individuals of gender k, sexual activity group l, age class i, at time t acquire infection with type h depends on the number of gender partnerships and the way they form partnerships with individuals of the opposite gender k 0 , the fraction of infected sex partners, and the transmission probability h k per partnership. The force of HPV infection h kli is given by 3.8 Mixing preferences

Mixing matrix
The way sex partnerships are formed is governed by the conditional probability matrix . Thus, klmij is the probability of someone of gender k, sexual activity group l, age class i having a partner from the opposite gender from sexual activity group m and age class j. This depends on the proportion of sex partners from the opposite gender from sexual activity group m and age class j, c k 0 mj N k 0 mj (0), in the total sexually active population. In generating the mixing matrix , the parameters 1 and 2 are used to depict the degree of assortative mixing between age and sexual activity groups, respectively. Thus, mixing is fully assortative ( is the identity matrix klmij = lm ij , where ij is the Kronecker delta) if 1 = 2 = 0 and proportionate when 1 = 2 = 1 [24,25,26,27]. The mixing matrix klmij is given by The model should satisfy the constraints balancing the supply of and demand for sexual partnerships: c klmij klmij N kli = c k 0 mlji k 0 mlji N k 0 mj . This is accomplished by specifying the mean rates of sex partner change as functions of the initial imbalance in the supply and demand of sex partnerships. Thus, The di¤erential e¤ects of cervical cancer-induced mortality are also likely to cause an imbalance between the demand for and supply of sex partnerships. There are few options for rectifying this. One option is to let the rates of sex partner change and mixing pattern of one gender vary over time so as to satisfy the above constraints. Another option is to …x the mixing patterns of both sexes and to let their rates of sex partner change vary over time so as to balance the supply of and demand for sex partnerships [25]. However, this latter option requires adding additional di¤erential equations that may considerably increase the size of the model. Because of this additional complexity only the former option is tried. Thus, In the sensitivity analysis, the gender that will be chosen …rst will be varied to test the robustness of the results.

Estimates of the mixing matrix
Even though the crucial role of the mixing matrix in the spread of many sexually transmitted infections has been repeatedly emphasized before [24,25,26,27], there are no adequate data to generate such a matrix. The current analysis follows previous work in this area by examining the range of patterns that are likely to arise in practice. This range is governed by the parameters 1 and 2 whose respective values are set to 0.6 and 0.7 in the baseline analysis and varied over a wide range in the sensitivity analysis. These estimates are obtained from the National Health and Social Life Survey (NHSLS) [55,63,64]. Higher values for 2 are reported for high-risk populations. For example, Garnett et al [26] estimated a value of 0.9 using data from a sample of patients with STD seen at the Harborview Medical Center. The baseline parameter values for the rate of sex partner change, strati…ed by gender, sexual activity, and age, are calculated from Table 3 using data from the NHSLS and the procedure outlined in Garnett and Anderson [24,25]. Brie ‡y, this procedure can be described as follows. Let the relative partner acquisition rate of sexual activity group l relative to the lowest group be pc l . Similarly, de…ne the relative partner acquisition rate of age group i relative to the lowest group as pa i . Therefore, the rate of sex partner change for people in age group 18-59 is where c 3 is the weighted mean rate of sex partner change rate.

Balancing population
To close the model, the total number of people in each gender category k, (k = f; m), age group i ( i = 1, 2, . . . ,17) and sexual activity group l (l = 1, 2, 3) must be equal to the sum of individuals in each epidemiologic class in the respective gender, age, and sexual activity groups. That is,

Estimates of other clinical parameters
The values of screening, diagnosis, and treatment parameters are reported in Tables 6.

Estimates of vaccine parameters
The e¢ cacy of the vaccine against incident infection (HPV 6/11 or 16/18) was assumed to be 90%. It was also assumed that infected vaccinated individuals do not progress to disease [52,78]. We assumed the vaccine does not a¤ect the natural course of disease. The duration of immunity conferred by vaccination is currently unknown. We assumed the duration of protection of HPV vaccination to be lifelong for the base case as was done in previous models [32] and examined a duration of 10 years in sensitivity analyses. Given HPV vaccination coverage is unknown, we assumed that 70% of adolescents will receive a 3-dose vaccine before they turn 12 similar to the coverage rates used in previous models [71,32].
Coverage was also assumed to increase linearly from 0% up to 70% during the …rst …ve years of the program and remain at 70% thereafter. We assumed that vaccine coverage for the catch-up program would increase linearly from 0% up to 50% during the …rst 5 years and then drop to 0% after 5 years.

Epidemiologic impact of screening and vaccination strategies
To

Years of life
The …rst …nal outcome measure is the total number of years spent alive by the active population. Thus, the discounted total number of years of life achieved using strategy a is given by where N kli is the size of the population of gender k, in sexual activity group l, and in age group i; is the discount rate; and T is the planning horizon.

Quality-adjusted life years
The second …nal measure of e¤ectiveness assigns quality of life weights to each health state and integrates the sum of all these quality-adjusted health states  Table 6: Cervical cytology screening and colposcopy characteristics and rates of cure and symptom recognition over the planning horizon (0; T ). Let qcin s , qcis s , qcc s , qccs, qgw k , and q ki denote the quality of life weights for an individual in the detected health state CIN stage s, CIS stage s, cervical cancer stage s, genital warts, and normal of gender k at age i; respectively. The discounted total number of quality-adjusted life years using strategy a over the planning horizon (0; T ) is given by ; dt: Note that the quality-adjusted years of life for females are reduced by time spent in diagnosed genital warts, CIN, and cancer states DCIN s , DCC s , DGW , and SCC. Males' quality of life deteriorates by spending time with detected genital warts. The probability of genital warts being recognized and treated is assumed to be 75%. It is assumed here that if a persons's health condition is not detected, the quality of life of that person will be the same as that of a person without the condition. This assumption biases the results against the vaccine. In the sensitivity analysis, the magnitude of the quality of life improvements for persons with undetected conditions prevented by the vaccine will be quanti…ed.

Estimates of quality of life weights
Women diagnosed with CIN1 and CIN2/3 were assumed to have quality weight of 0.91 and 0.87, respectively [62,54]. The quality weight for genital warts is assumed to be 0.91 [62]. Females with local and regional cancer are assumed to have a quality of life weight of 0.76 and 0.67, respectively [62]. A quality weight for invasive distant cancer of 0.48 was derived from Gold et al [30] using the 25th percentiles of female genital cancer weights. It is assumed that the quality of life for cervical cancer survivors after successful treatment will continue to be lower (at 0.76) than that of healthy females [4,83]. Undiagnosed HPV, genital warts, CIN, and cervical cancer states and successfully treated CIN states are assumed to have a quality of life weight similar to those of individuals without HPV disease. Gender-and age-speci…c quality weights for other health states were derived from Gold [30]. Similar values were reported from the Beaver Dam Health Outcomes study [23]. CIN and cancer health states were multiplied by the age-and gender-speci…c weights to re ‡ect the variation in quality of life by age and gender groups.

Economic consequences of screening and vaccination strategies
The total costs of each strategy includes costs of cytology screening per unit time, cost of vaccination, lifetime cost of treating detected genital warts, CIN and invasive cancer cases, and the cost of following false positive results of screening.

Screening costs
The cost of cytology screening per unit time is the product of the cost per test scn, the test compliance rate cover ib given the frequency of administering the test per unit time (e.g., every year), and the size of the population eligible for screening P For simplicity, it is assumed that females in the hysterectomy class are not screened. However, this may not be the case as suggested by recent studies [70]. The cost of following false positive results of the cytology test is the product of the cost of colposcopy colp of those females who do not have a repeat cytology test, one minus cytology speci…city papsp and the size of the screened population that is truly negative P

Total costs
Discounted total cost over the planning horizon (0; T ) of following strategy a is [Screen a (t) + T reat a (t) + V accinate a (t)] e t dt:

Cost-e¤ectiveness ratio
To compare mutually exclusive vaccination strategies a and a 0 , we calculate the incremental cost-e¤ectiveness ratio [82] Cost a Cost a 0 QALY a QALY a 0 : 6 Analysis using the model 6.1 Simulations with the baseline estimates of the parameters Mathematica R (Wolfram Research, Champaign, IL) version 5.2 was used to generate numerical solutions of the model. The NDSolve subroutine in Mathematica is a general numerical di¤erential equations solver. Since the model consists of non-sti¤ ODEs, the Explicit Runge Kutta methods, with adaptive embedded pairs of 2(1) through 9(8), provide accurate and less expensive solutions [85]. Other methods such as the Predictor-Corrector Adams method, with orders 1 through 12, produced the same results, but took longer to compute the solution.
The following strategy for simulations was followed. First, the baseline parameter estimates were used to solve the model for the pre-vaccination steady-state values of the variables. Second, the pre-vaccination data were used as initial values for the vaccination model and the model was solved for the entire time path of the variables until the system approached the steady state (approximately 100 years). The solution approximates the potential impact of various HPV vaccination programs, including routine vaccination of 12-years old individuals. Finally, once the solution is obtained the results can be presented for various outcomes in many di¤erent formats.

Model validation
The validity of a complex model like this cannot be established directly. Instead, its face validity may be judged by how reasonable model assumptions are [34,81]. In the process of building this model, we comprehensively reviewed previous relevant models and consulted experts on the natural history of HPV infection and HPV-related diseases. A comprehensive review of the literature was conducted to identify studies to inform model inputs. To facilitate independent review of the model and the ability to replicate its results, all model equations and inputs are made available. All model equations and inputs are programmed in Mathematica TM (Wolfram Research, Champaign, IL). A series of tests were performed to debug and establish the technical accuracy of the Mathematica programs. For example, the sum of the number of individuals of a given gender, age, and sexual activity group in each compartment is veri…ed to be equal to the total number of people N kli at each point in time (see section 3.9 on balancing population). Finally, the predictive validity of the model was evaluated by looking at age-speci…c HPV prevalence, CIN, genital warts, and cervical cancer incidence rates predicted by the model and comparing them with those reported in the literature [29,47,73,74,75,41,45]. The model predictions were well within the range of values found in the literature. For example, the predicted HPV 16/18 attributable cervical cancer incidence curve in the absence of screening had a shape and magnitude at peak (55.9 per 100,000 women years for age 50-54) similar to that estimated for unscreened populations [33,58].