An economic and disease transmission model of human papillomavirus and oropharyngeal cancer in Texas

In 2017, 46,157 and 3,127 new oropharyngeal cancer (OPC) cases were reported in the U.S. and Texas, respectively. About 70% of OPC were attributed to human papillomavirus (HPV). However, only 51% of U.S. and 43.5% of Texas adolescents have completed the HPV vaccine series. Therefore, modeling the demographic dynamics and transmission of HPV and OPC progression is needed for accurate estimation of the economic and epidemiological impacts of HPV vaccine in a geographic area. An age-structured population dynamic model was developed for the U.S. state of Texas. With Texas-specific model parameters calibrated, this model described the dynamics of HPV-associated OPC in Texas. Parameters for the Year 2010 were used as the initial values, and the prediction for Year 2012 was compared with the real age-specific incidence rates in 23 age groups for model validation. The validated model was applied to predict 100-year age-adjusted incidence rates. The public health benefits of HPV vaccine uptake were evaluated by computer simulation. Compared with current vaccination program, increasing vaccine uptake rates by 50% would decrease the cumulative cases by 4403, within 100 years. The incremental cost-effectiveness ratio of this strategy was $94,518 per quality-adjusted life year (QALY) gained. Increasing the vaccine uptake rate by 50% can: (i) reduce the incidence rates of OPC among both males and females; (ii) improve the quality-adjusted life years for both males and females; (iii) be cost-effective and has the potential to provide tremendous public health benefits in Texas.


Results
Model validation. The model consisted of a large system of ordinary differential equations (ODEs) similar to the previous work of Elbasha and Peng 18,19 , which provided a solid foundation for depicting the population dynamics, HPV transmission, and vaccination simultaneously. In this study, the epidemiologic model of HPVassociated OPC was incorporated and the predicted HPV-associated OPC incidence rates were used for model validation.
The parameters and initial variables collected from Year 2010 in Texas used to validate the model. The predicted age-specific incidence rates (per 100,000) compared with the actual data for Years 2011 and 2012 in Texas, and both the real incidence rates and the predicted incidence rates were calculated (Fig. 1). The two curves matched each other reasonably well and had similar trends for all age groups, although the match was better for 2012 than for 2011, especially for males older than 45 years (Fig. S4). For males, the predicted result was slightly larger than the actual data for the Texas population older than 55 years; for females, the difference observed only for the population older than 60 years. Prediction. Age-adjusted incidence rates predicted for 100 years and all age groups were included. While other modelers such as Chesson et al. 20 set the starting age group at 12 years, Graham et al. 21 listed the vaccine target age for starting HPV vaccination for females in the U.S. at 11 years. Therefore, age-adjusted incidence rate was predicted by covering age groups started from 11 years (Supplemental Fig. S5). Because disease transmission occurs mainly among adults, our age-adjusted incidence rate prediction figure starts at age 20 years (Fig. 2).
The initial parameter values were determined for 2015 and calibrated for prediction of age-adjusted incidence rates. The parameter values fixed for the 100-year study time horizon. The predicted incidence rate (per 100,000) increased until Year 2030 and then declined for the remainder of the time span. The predicted male incidence rates were much higher than for females, which clearly reflects the known differences in OPC between genders. The predicted incidence rates covered Years 2019 to Year 2115 (the results for years 2016 to 2018 found to be very sensitive to the initial variable values and thus may be artifacts).
Cumulative cases prevented in Texas are presented in Table 1. Although increasing the vaccination coverage rate by 50% for females only did not affect males directly, vaccinating females resulted in a notable decrease in OPC among males. For example, increasing the vaccine uptake rate of females decreased the respective cumulative number of OPC cases for males and females by 1271 and 1297, respectively, within 100 years. If vaccination uptake rates increased for both females and males, providing both direct and indirect protection, even more cases  www.nature.com/scientificreports/ of disease among men were prevented; for example, 2362 cases among men and 1450 cases among women were prevented after 100 years of vaccination. Additionally, more cases expected prevented with a longer follow-up time. For instance, compared with 10 years, the number of OPC cases prevented among men and women was around 25 times and 20 times higher, respectively, over 100 years when vaccination coverage rate was increased for males and females.
Policy assessment and sensitivity analysis of model parameters. Policy focuses on the health gains from investments to increase the vaccination uptake rate among the target groups in Texas. We assessed the impact of different vaccination strategies or scenarios on health and economic consequences, such as HPVrelated OPC incidence rates, years of life, QALYs, and costs. Incremental cost-effectiveness ratios (ICERs) estimated for different vaccination coverage rates, vaccination covered age groups, and whether male vaccination rates increased in addition to those of females. The effects of rate of sexual partner change and detection rate of the OPC parameters on health and economic outcomes are unknown or highly uncertain and therefore assessed with sensitivity analysis.
Vaccination coverage rate. Vaccination adherence is one of the most important factors in determining the level of disease prevention. In our model, age and gender specific vaccination uptake rates, including uptake rates of first dose and second dose, were collected from CDC 22 and the values were presented in Supplements Text S2 and S3. In addition to predicting the impact of vaccination with first dose of vaccine, we considered the adherence to the second vaccine dose. Results derived for increasing the first and second vaccine uptake rates by multiplying 1.5 to each value of vaccine uptake rate respectively or simultaneously. As the model predicted results since 2016, the intervention started in Year 2016.
First dose. As shown in Fig. 3a and b, when the first dose vaccine uptake rate increased 50% for females only, the incidence rate (per 100,000) decreased by 2.79% for males and 5.30% for females by Year 2031 compared with the baseline predictions. When the first dose vaccine uptake increased by 50% for both males and females, the incidence rate (per 100,000) decreased by 6.91% for males and 10.04% for females by Year 2031, again compared with the baseline predictions.
Second dose. By changing the proportion of those with uptake receiving two doses, the incidence rate changed less than the incidence rate changes when we varied the first dose uptake rate (Fig. 3c,d). Compared with the baseline results, the incidence rate decreased by 0.10% for males and 0.12% for females by Year 2031 when we increased the vaccination uptake rate of females only; the incidence rates decreased by 2.36% for males and 3.15% for females at Year 2031 if uptake rates of both females and males were increased.
First and second dose. We also performed sensitivity analysis by changing the vaccine uptake rate for the first and second doses at the same time. As shown in Fig. 3e and f, the incidence rate (per 100,000) decreased by 2.89% for males and 5.44% for females by Year 2031 when we increased the vaccine coverage for females only. The incidence rate (per 100,000) decreased by 9.31% for males and 13.05% for females by Year 2031 when we increased the vaccine coverage rate 50% for both males and females.
Vaccination starting age. The starting age of vaccine uptake has a notable effect on incidence rate. We predicted the incidence rates for three different scenarios: (1) the current vaccination coverage strategy, starting from age 13-14 years for both males and females, (2) assuming that people would receive vaccination two age groups earlier, that is, age 9-10 years for males and females, and (3) assuming that people would receive vaccination two age groups later, that is, age 18 years for males and females (Fig. 4). Compared with the baseline results, the incidence rates (per 100,000) increased by 3.11% for males and 4.98% for females by Year 2031 when the vaccination started two age groups later. When the vaccination started two age groups earlier, the predicted incidence rates (per 100,000) were almost the same as the baseline results for both males and females over 100 years, with a decrease of 0.54% and 0.85%, respectively.
Rate of sexual partner change. Sexual behavior is an important contributor to HPV-related OPC incidence rate 18 . We increased and decreased the rate of sexual partner change by 50% from baseline; Fig. 5 shows that the www.nature.com/scientificreports/    www.nature.com/scientificreports/ Detection rate of OPC. Because of the unavailability of detection rates for OPC, we performed sensitivity analyses to assess the impact of detection rates on the disease incidence predictions. Figure 6 shows that when the detection rate was set at 0.8, the incidence rates increased by 2.32 cases (per 100,000) for males and 1.30 cases (per 100,000) for females in Year 2031 compared with the detection rate of 0.6. When the detection rate was set at 1, the incidence rates increased by 5.94 cases (per 100,000) for males and 3.36 cases (per 100,000) for females in Year 2031 compared with the detection rate of 0.6.
Economics. Incremental cost-effectiveness ratio.   www.nature.com/scientificreports/ Table 2 shows the short-term to long-term vaccine costs, promotion costs, and treatment costs. There was no tremendous increase for vaccine costs (0.99% in 10 years and 0.77% in 100 years) if uptake rate was increased 50% for females only. However, the extra costs were notable if the uptake rate increased 50% for both males and females, with a 6.89% increase in 10 years and a 4.33% increase in 100 years. The promotion costs were the major part of vaccination cost when we increased vaccine uptake rate: $1.6 billion in 10 years for increasing female only rate and $3.3 billion for increasing both gender rates. However, with the increase of vaccine uptake rate, the treatment costs would decline by 2.27% for the female only increase strategy and 3.87% for the both gender increase strategy in 10 years; the treatment cost would decrease by 3.94% and 9.20% in 100 years. Table 3 provides results for the cost-effectiveness analysis, which compared three vaccination coverage levels and three vaccination costs. The criteria for a strategy being cost-effective followed the convention of $100,000 per QALY gained 23 . For each vaccination scenario, the strategies ordered from least resource intensive (i.e., current vaccine uptake rate) in the top row to the most resource intensive (i.e., increase vaccination uptake rate 50% for males and females) in the bottom row. We varied combinations of vaccine acquisition cost and vaccine promotion uptake costs in a sensitivity analysis of the ICERs. Base-case vaccination cost values were derived from CDC data ($178.14 for the vaccine; $50.64 and $20.26 for promoting uptake of the first and second dose, respectively) 24 . The best-case cost assumptions ($142.51 for vaccine; $20.26 for extra first dose and second dose costs) and the worst-case cost assumptions ($217.11 per dose for vaccine; $87.92 and $20.26 for promoting first dose and second dose) were calculated 25 . The base-case results showed that changing from a strategy of a 50% increase over current vaccination levels among females only to a strategy of a 50% increase among females and males (Age-specific base vaccine uptake rates were presented in presented in Supplements Text S2 and S3) reduced the ICER from $126,689.20/QALY to $94,517.75/QALY. The variation of vaccine costs also had an important impact on the ICERs. For example, the ICERs increased to $194,065.98/QALY (baseline vs. strategy 1) and $144,497.70/QALY (strategy 1 vs. strategy 2) when the vaccination cost was $217.11 and dropped to $71,734.45/QALY and $53,718.93/QALY when the vaccination cost was $142.51.  www.nature.com/scientificreports/ To assess the impact of uncertainty of some model parameters on the ICERs, other scenarios were considered, such as no discounting of costs and QALYs or no promotion costs for the expansion of vaccine uptake (Table 4). In almost all scenarios, a strategy of increasing vaccine uptake by 50% for females only was dominated by the strategy of increasing both female and male vaccine uptake by 50%, except in the setting with no promotion cost and vaccination cost of $217.11. The strategy of adding male vaccination uptake rate was cost-effective in most of scenarios, especially for the setting of no discounting and no promotion cost. Finally, we also computed long-term ICER results for Years of Life Saved and the short-term to long-term costs, QALYs, and ICERs (See Supplemental Tables S3, S4). Table S4 showed that the increase of vaccination uptake rate became cost-effective after 50 years with the threshold of $100,000 per QALY and strategy of adding female and male dominated strategy of increasing female only after 50 years.
Cost-effectiveness acceptability curve. The uncertainty in the results of the cost-effectiveness analysis were summarized in CEACs, which are provided in Supplemental Figs. S6 and S7. The curves indicated the percent of simulations for which the strategy was cost-effective for a determined monetary value that the policy-maker would be willing to pay for a QALY saved. For instance, the CEAC for changing the vaccination coverage rate showed that increasing uptake rate for females only was cost-effective in all simulations if the willingness to pay was higher than $135,000 per QALY. And all simulations were cost-effective when the uptake rate was extended to include males at this monetary threshold. With willingness to pay reached $100,000 per QALY, the strategy of increasing female and male vaccine uptake rates was cost-effective in all the simulations.
Relative to the commonly reported threshold of $100,000 to $150,000 per QALY saved 23 , the results suggest that substantially increasing vaccination against HPV was very cost-effective, especially when the vaccination uptake rate was increased among males and females. The efficiency of increasing vaccination of males was driven by their relatively high and increasing OPC incidence rates.

Discussion
In this study, we built an age-structured HPV infection model to estimate the impact of HPV vaccination on OPC progression in Texas. We obtained model parameter values from the literature or calibrated them with data from Year 2010 and Year 2015 in Texas. The age-specific incidence rates predicted by the model validated against OPC incidence for Years 2011 and 2012 in Texas, and they matched the real data well. The model then applied to predict the age-adjusted incidence rate for Texas until 2115. Multiple sensitivity analyses conducted on several of the parameters, including vaccine uptake rate, coverage age and sexual activity parameters, to determine the impact of HPV vaccination and sexual behavior on disease progression.
Our work quantified the trend of OPC incidence rates for Texas. The model showed an expected incremental increase on predicted incidence rate for the first 20 years and then a decline for the rest of prediction period (80 years). While estimated tends in OPC incidence rates and the cost-effectiveness of increasing vaccination rates for boys and girls were similar to previous estimates of U.S. trends, there were some differences in the Texas results. A Canadian study 21 examined the cost-effectiveness of immunization of 12-year-old boys with a Markov model that did not account for herd immunity. With vaccine efficacy at 99 percent and 70 percent uptake, 0.05 QALYs and $145 estimated saved per person in the target population, or $28 mil over the lifetime of the cohort of 192,940 individuals. Vaccination of boys continued to be cost saving when efficacy and uptake were assumed to be 50 percent, but total savings declined to $8 mil. in Canada.
This study has some limitations. First, the study focused on only one of several diseases related to HPV infection. We focused on OPC because its incidence is increasing and it has overtaken cervical cancer as the leading health problem likely to decrease due to HPV immunization. Researchers have demonstrated the value of preventing other cancers, genital warts, and recurrent respiratory papillomatosis with HPV vaccination 26,27 . Addition of these diagnoses to the current work would make the increased investment in vaccination even more cost-effective and possibly cost saving, as shown in the work by investigators at the CDC 17 . Second, we did not attempt model fitting because of poor data availability and the lack of efficient regression techniques for large-scale ordinary differential equation models. Therefore, the parameter values obtained from the literature or were calibrated based on the actual data and assumptions. A better match would be possible if model fitting becomes feasible in the future. Third, the vaccination uptake rate held constant during the 100-year prediction period, which may not reflect the real evolution of vaccine policy and vaccination behavior. Fourth, because of the lack of population data, the validation was limited to only data for 2011 and 2012 in Texas. With more data in the future, we will improve model validation. For example, the predicted age-specific incidence rate is a little higher than the actual incidence rate in people aged 60 to 79 years, and the difference becomes smaller in people aged 80 years and older. Fifth, there is no dependable literature on the rates of detection of OPC. In one study 28 , the detection rate for oral cancer was 86%. We estimated that the detection rates for local and regional stages of OPC were lower than 86% and assumed the detection rate for distant stage OPC as 1. Therefore, detection rates were assumed 0.68, 0.85, and 1 for local, regional, and distant OPC. Finally, the predicted incidence rate declined briefly in the beginning of the prediction period; this drop may have resulted from the initial conditions.

Conclusion
The paper adds to the literature on health/economic policy analysis of HPV immunization by adapting a comprehensive population disease model for a state-level assessment of the long-term health and economic consequences of increasing adherence to HPV immunization guidelines. We incorporated up-to-date treatment and vaccine cost estimates for OPC based on Texas health insurance claims, published reports, up-to-date disease and sexual behavior data and cost estimates for promoting vaccine uptake among males and females to achieve a 50% increase in the HPV vaccination rates. Thus, this work provides information that may be actionable at  Table 5; and Table 6 contains initial variables' notations, definitions, and the sources of initial values. Our complete model equations provided in Supplementary Text S1; additionally, model structure diagrams given in Supplementary Figs. S1-S3. While the basic structure of our model follows the published model 18,29 , the natural history of OPC was modeled and incorporated, and the entire framework was updated and adapted to the Texas population for making specific predictions for Texas. In the demographic model, people transfer from their current age group to the next age group, except for the oldest age group (≥ 85 years), at age-and gender-specific rates. For instance, the calibrated transfer rate for Texas suggested that about 12.4% of females from the second age group (age 1-8 years) move to the third age group (age 9-10 years) each year, compared to 13.4% of males. The population size for every age group depended on the probability of transferring to the next age group and the probability of transferring from the younger group, and on the OPC-related and non-OPC-related death rates. The population growth for the first age group (0-1 years) calculated based on the birth rate because there is no younger age group, and there was no transfer out for the oldest age group. In the epidemiological model, HPV transmission specified by gender, age, and sexual activity. The important subpopulations related to HPV transmission included the susceptible population, infected population, persistently infected population, vaccinated population, infectious vaccinated population, persistently infected vaccinated population, and recovered vaccinated population. Additionally, since the three-dose HPV vaccination www.nature.com/scientificreports/ strategy has been changed into a two-dose strategy recently, the force of HPV infection (λ) was redefined by 19 , which was determined by the number of sexual partners, relative infectivity of vaccine breakthrough cases, and probability of people being in different groups.
In the OPC natural history model, the OPC stage contained detected local, regional, and distant cases and survived cases. We updated parameters of published models to reflect current information on HPV transmission and OPC progression. For OPC progression, equations now reflect the epidemiologic variations among different stages of detected OPC. The population of individuals with detected local OPC determined by progression rate, detection rate, and various subpopulations, as given below: where l denotes sexual activity group, c denotes gender, and i denotes age group. Parameter p L denotes the detection rate of local OPC and p R and p D denote the detection the rates of regional and distant cancers, respectively. We multiplied all subpopulations by their detection rates to predict the prevalence of detected OPC. Here, d was the population transfer rate and i = 1, 2, 3, . . . , 23 . As no younger group exists when i = 1 (age group 0-1 years), there is a small difference between the two equations. When i was greater than 1, the individuals with detected local OPC were from the detected local OPC cases in the first age group, HPV-infected individuals, persistently infected vaccinated individuals, infected vaccinated individuals with waned immunity, persistently infected individuals, infected vaccinated individuals who had undergone tonsillectomy, and infected vaccinated individuals who had not undergone tonsillectomy. We subtracted the population that died of any cause except OPC, the population that transferred to the next age group, the population that died of OPC, the population cured of OPC, and the population with detected local OPC that developed regional OPC; θ tR represented the progression rate from local to regional status. Equations given below model how the disease develops from local to regional: For the detected regional disease, when the disease progressed from regional to distant, we used the following equations: + p L * θ * U l,1,c + p L * θ hw * Hw l,1,c + p L * θ hy * Hy l,1,c www.nature.com/scientificreports/ In the equations, θ tL , θ tR , and θ tD represent the progression rates for three stages for an individual with OPC. Specifically, θ tL denotes the progression rate from HPV infection to local OPC, θ tR denotes the progression rate from local to regional OPC, and θ tD denotes the progression rate from regional to distant OPC. The median time of disease progression in all age groups was around 8 months 30 . The progression rates were age-specific, as the disease usually progresses more quickly in the older groups, and the θ tR and θ tD were larger than the θ tL since the disease progression accelerates when the disease progresses to regional or distant status compared with local status.

Economics.
Applying the 2010 and 2015 Texas base vaccination rates [31][32][33] , the incremental cost-effectiveness ratios (ICERs) were calculated for 50% increases in the female vaccination rate only and in the female and male rates (Age-specific base vaccine uptake rates presented in Supplements Text S2 and S3). The 50% increase was an important and feasible target and was also examined in a Canadian analysis of HPV immunization expansion 21 . Fifty percent increase in immunization would move Texas close to the 2020 Health People goal for the U.S. 34 . ICERs were computed by dividing the net economic costs (vaccine plus vaccination promotion costs minus cancer treatment cost averted) by life years and quality-adjusted life years (QALY) gained. Sensitivity analyses altering the vaccination strategies and HPV transmission rates examined the impact of behavior and vaccination on the incidence of OPC. To estimate the ICERs for the vaccination strategies, we calculated the vaccination cost, treatment cost, and QALYs based on the following formulas. All costs valued in 2018 U.S. dollars.

Vaccination cost.
Total vaccination costs included the unit cost per vaccination and the number of people vaccinated. Let B l,c denote the newborn population; φ l,c denote the proportion of newborns vaccinated; φ l,i,c represent the vaccine uptake rate for the first dose; Φ 1 and Φ 2 represent the proportion of receiving only one dose or two dose vaccine respectively; and Y l,i,c , X l,i,c , U l,i,c, denote the population vaccinated 18 . The cost of the HPV vaccine was $178.14 per dose 25 . cost pro_1 and cost pro_2 represent the extra costs for motivating people to obtain the immunization and complete the second dose, and they were included in the vaccination cost calculation; the estimated costs for promoting the first and second doses were $50.64 and $20.26, respectively [35][36][37] . pop 1 and pop 2 denote the corresponding target population for promotion. Sensitivity analysis assessed the ICER estimates over a range of vaccination costs and efficacies.
Cost-effectiveness ratio. To compare vaccination strategies a and a ′ , the ICER was calculated as follows: Here Cost a = T 0 (Vaccinate a (t) + Treat a (t)) , and the corresponding value for vaccination strategy a' was Cost a ′ = T 0 (Vaccinate a ′ (t) + Treat a ′ (t)) . Commonly, the lost quality-adjusted life years of the baseline strategy was larger than the strategy of increasing female vaccination or increasing vaccination for both genders, and the cost of the baseline strategy was smaller than the other improvement strategies. In practice, QALY lost,,a ′ − QALY lost,a applied instead of QALY lost,a − QALY lost,a ′.
Probabilistic sensitivity analysis. We computed the cost-effectiveness acceptability curve (CEAC) by performing probabilistic sensitivity analysis. Parameters for vaccine uptake rate and quality of life for patients with OPC were included in probabilistic sensitivity analyses. For quality of life, beta distribution was applied to generate 1000 random samples as inputs (α = 68.3, β = 27.24). For vaccine uptake rate, the Latin Hypercube sampling technique based on multivariate normal distribution conducted to generate random samples as input values. Because of the computational complexity, only 100 random samples for each parameter were generated as inputs in the simulation.
Model initial variables, parameters, and calibration. All of the values of variables, parameters were listed in Supplementary Text S2 (Year 2010) and S3 (Year 2015). For parameters that were not available in the literature, we calibrated from reasonable data or knowledge of HPV infection and the epidemiology of OPC. Seroconversion following oral HPV infection appears to be rare among both males and females [39][40][41][42][43] . The HPV vaccines function primarily by inducing durable systemic virus-neutralizing antibody responses, thereby preventing infection and possibly guarding against reinfection and disease 13,21 . Therefore, we assumed the corresponding parameters (ls, θ tw1 , θ tw2 , θ tws , θ tp1 , θ tp2 , θ tps ) to be 0. Implementation and computation. The model structure was implemented in MATLAB (MathWorks, Natick, MA) and the ode15s solver was applied to obtain numeric solutions. Both the relative and the absolute error tolerances were set as 1.0 × 10 -7 in the ordinary differential equation solver; we also set maximum step size as 1.0 × 10 -2 and restricted the total number of persons (variable N) as non-negative.
Simulations were performed for the Texas population. The entire population was stratified into 23 age groups, two gender groups, and three sexual activity groups. The variable values of initial conditions were from Year 2010 in Texas and we calculated the age-specific incidence rates (per 100,000) and age-adjusted incidence rates (per 100,000) for validation, prediction and sensitivity analyses.
To check the validity of the model, we compared the predicted age-specific incidence rates (per 100,000) with the real OPC incidence rates (per 100,000) from the Surveillance, Epidemiology, and End Results Program (SEER) database for Texas and the Texas Cancer Registry database. Because of limitations on data availability, comparisons conducted for Texas in 2011 and 2012 only. One-way sensitivity analyses implemented by changing the value of some parameters according to pre-specified ranges. Cost a − Cost a ′ QALY lost,a − QALY lost,a ′