Modeling Hepatitis C Elimination Among People Who Inject Drugs in New Hampshire

Key Points Question What improvements in the hepatitis C (HCV) care cascade are required to eliminate HCV among people who inject drugs (PWID)? Findings This decision analytic model of HCV transmission found that improved testing, treatment, and access to harm reduction were all associated with reductions in HCV prevalence and mortality among PWID. Improvements in both testing and treatment were associated with HCV prevalence of less than 2% by 2030. Meaning These findings suggest that HCV elimination may be possible among PWID by 2030 with improved testing and treatment; improved harm reduction may reduce the time and number of treatments required to achieve similar outcomes.

This supplemental material has been provided by the authors to give readers additional information about their work.

New Hampshire Context
New Hampshire (NH) provides an interesting setting for studying HCV elimination among PWID in the United States (US) because of its historically low rates of HCV testing and treatment and its high rates of injection drug use during the opioid crisis. Though the rate of injection drug use in NH is high compared to most other states, the lack of infrastructure for testing, treatment, and harm-reduction services appears to be common throughout the US. This section provides some background on the current state of injection drug use, HCV prevalence, HCV testing, HCV treatment, and harm-reduction efforts in NH and, when possible, how NH compares to other states.

Injection Drug Use in New Hampshire
Rates of injection drug use rose substantially in NH during the opioid crisis, with a clear increase in past-year heroin use from 2013-2018 observed in the National Surveys on Drug Use and Health. 1 This increased use of injection drug use is reflected in a rapid increase in overdose deaths during this period (eFigures 1a and b). Heroin and fentanyl are now responsible for the majority of drug overdoses in the state 2 (eFigure 1b). The majority of overdose deaths, ED visits, and EMS naloxone occur among individuals age 20-49, mirroring national trends in HCV incidence attributable to PWID. 3 Rising rates of complications related to injection drug use -such as infective endocarditis or outbreaks of HIV similar to that seen in Scott County, Indiana -have also been noted. 4,5 NH has consistently been among states with the highest overdose mortality since 2014. 6 A 2017 "hotspot" report prepared by the National Drug Early Warning System studied this issue in some detail, identifying reasons for the rapid growth in overdoses from 2010-2016. Reasons identified include high rates of substance abuse, high historical rates of opioid prescribing, and poor access to substance abuse treatment and harm reduction programs (more below). 7 HCV Prevalence, Testing, and Treatment: Recent work by Romo et al. suggests HCV prevalence among PWID in NH is high though consistent with previous studies in the US. 8,9 HCV was not a reportable disease in NH until November 3, 2016, and there is no state-reported HCV prevalence data before then. Moreover, little testing has been done in potentially high-contact settings for PWID such as syringe service programs (SSPs) or medication-assisted treatment (MAT) programs. 10,11 Medicaid restrictions on treatment were strict in NH before 2016: patients were required to have F3 or greater liver disease, six months of sobriety and prescribers had to be a specialist or a physician with continuing medical education in HCV treatment. 12 These restrictions were relaxed in 2016, but prescribers must still document addiction treatment and periodic screens for alcohol use. 12 Such restrictions have limited HCV testing and treatment in NH and impeded HCV elimination efforts. 8,13 Though comprehensive data on PWID HCV testing and treatment in the US is not available, many studies suggest that insufficient rates of testing and treatment are common throughout the US. 8,[14][15][16] Recent published editorials from HCV experts describe a paucity of testing and treatment in the US for PWID, 17,18 though they do not comment on state-by-state specifics. Taken together, this suggests the low rates of testing and treatment in NH are common in other states as well.

Harm-reduction
There were no legally operating SSPs in NH before 2017. 19 Since SSPs were legalized in NH in 2017, there has been continued growth in the number of SSPs and the number of total syringes distributed in NH. 11 Previous studies demonstrate significant heterogeneity state-by-state in SSP coverage, leaving many PWID with insufficient access to clean needles. 20,21 Access to SSPs has grown during the opioid crisis -in NH and across the country -but legal barriers remain in many places. Overall, the evidence suggests that poor access to SSPs is common in the US. 20

New Hampshire Surveillance Data
Precise data on PWID are lacking both in NH and across the US. Our analysis relies on state-specific data on the HCV burden among HCV collected from a variety of sources cited throughout the main paper and supplement. In this section, we provide further detail on these state-specific sources and their methodologies. We also used data from the National Surveys on Drug Use and Health (NSDUH) and National Survey of Substance Abuse Treatment Services (N-SSAT), both available at the Substance Abuse and Mental Health Services Administration (SAMHSA) website. 24

Reportable Communicable Diseases in New Hampshire
As mandated by state law, the NH Department of Health and Human Services regularly collects and publishes case counts of select communicable diseases in the state. 25 Acute and chronic HCV incidence has been included in this report since November 3, 2016. Chronic HCV case counts in this report represent new diagnoses of chronic HCV. Data are also collected on whether HCV cases (acute and chronic) are related to injection drug use. HCV incidence stratified by injection drug use status is available from the state upon request. Stratified data was used for model validation of HCV incidence as described in the Validation section below (data as of 2/3/2021). 26 The communicable diseases report series is based on information reported to the Department of Health and Human Services and may represent an underestimate of the true absolute number and incident rates of cases in the state.

Drug Injection Surveillance and Care Enhancement for Rural Northern New England (DISCERNNE) Study
The DISCERNNE study is a multisite, multidisciplinary study conducting epidemiological and policy assessment of opioid-related outcomes, prevention and care services, and laws in eleven contiguous counties in rural MA, VT, and NH. Data was collected using both quantitative and qualitative methods including HCV testing, quantitative and social network survey, and PWID interviews on drug use and treatment histories. Participants were recruited using respondent-driven sampling and included 199 NH PWID. The details of this study have been published previously and provide useful context on the opioid crisis in the rural northeast. 8,27 New Hampshire Drug Monitoring Initiative (DMI) The DMI's stated purpose is to provide awareness and combat drug distribution and abuse. The DMI obtains data from various sources (including Public Health, Law Enforcement, and EMS) and provides monthly reports for stakeholders as well as situational awareness releases as needed. These data include overdose deaths, ED visits, EMS naloxone administration, and treatment admissions stratified by age and location. 28 where each term above is described below.

Inflow of PWID (IN)
PWID enter the model at rate θ. They are assumed to be initiating injection drug use for the first time (n=1); they are known to not have HCV (d=1), do not have liver disease (k=1), and are not enrolled in MAT or SSP (i=0, j=0). A fraction ϕ of these entering PWID are comorbid stimulant users (m=1); the remainder are not (m=0). Thus The inflow rate θ varies over time, with one inflow rate before 2013, another for the time of the opioid epidemic in 2013-2018, and afrer 2018. See the Model Calibration section below for more details. We assume ϕ is 16%. 28 In the Monte Carlo simulation, ϕ is drawn from a uniform distribution ranging from 10%-22%.
Injector Duration Progression (ID) eFigure 2 depicts the structure of the injector duration submodel. PWID progress through PWID duration categories at rates unless they become abstinent; the cessation rate is . Abstinent PWID relapse at rate . PWID exit the model with mortality rates and for active and inactive injectors; we assume higher mortality rates for active injectors The injector duration (ID) terms in the model are thus:

Comorbid Stimulant Use (CS)
The model recognizes comorbid use of stimulants such as cocaine or amphetamines as leading to higher risk for HCV infection. 29 PWID go through periods of stimulant use, entering and exiting at rates σ and ζ respectively. Thus the comorbid stimulant (CS) terms in the model are: The entering rate is selected to match a specified steady-state fraction f2 in this high-risk category.
The assumed parameter values are as follows:

MAT/SSP Enrollment (IT)
eFigure 3 depicts the structure of the MAT/SSP submodel. We assume entering or leaving one form of harm reduction is independent of the other. Active PWID enter MAT and SSP at rates β and η, respectively. PWID exit MAT and SSP at rates γ and κ respectively. Thus for active injectors (n=1, 2, 3), the intervention terms (IT) in the model are Inactive injectors (n=4) do not enter SSP and exit immediately upon becoming abstinent. The IT terms for inactive injectors are as above but with = 0 and ≈ ∞.
The entering rates β and η varied across different policy scenarios and were chosen to match steady-state enrollment targets for the given scenario. The assumed leaving rates were as follows: Parameter The scenario-specific steady-state enrollment targets were set as follows:  Pre-2013: MAT Enrollment Target = 7.5% SSP Enrollment Target = 0% The pre-2013 MAT enrollment target is estimated from an observed 1.53x increase in MAT participation from 2013-2020. 10, 23 We assume continued growth in MAT coverage in the base case such that the 2022 target is twice the pre-2013 target. There were no SSPs in NH before 2013 and this is assumed to be 0%.
 Base case harm prevention scenarios (Scenarios 1, 3, 4, 5), also 2013-2022: MAT Enrollment Target = 15% 36 SSP Enrollment Target = 15% 11 Base case harm prevention targets are based on reported opioid use disorder treatment rates nationally (19.7%), adjusted for below average rates of buprenorphine-waivered providers and rates of admissions for substance use in NH. 7,36 Base case SSP coverage target is based on state-reported data for the number of clean needles distributed. 11  Increased harm prevention scenarios (Scenarios 2, 6): MAT Enrollment Target = 35% SSP Enrollment Target = 50% Increased harm prevention scenario targets were selected based on prior work by Fraser et al. who assume targets of 50% for both MAT and SSP; we assume a more modest scale up in MAT participation. 37  Intermediate harm prevention scenario (Policy Sensitivity Analysis of eFigure 9): MAT Enrollment Target = 25% SSP Enrollment Target = 32.5%

Infection (FOI)
Acquiring a chronic HCV Infection in the model corresponds to uninfected people (those with d=1,…, 5) moving from compartment , , , , to the corresponding infected compartment , , , , . The force of infection (FOI) terms in the model capture this process. For d=1,…, 5, the FOI terms are of the form , , , is a weight representing additional risk or protection associated with specific categories where  is the relative weight by injector category; here = 1 because weights are relative to long-term injectors and = 0 because inactive PWID do not face transmission risk.
 Ξ , Γ , and Π are relative adjustments for comorbid stimulant use, enrollment in MAT and enrollment in SSP, respectively. Ξ = Γ = Π = 1, defining the base of the relative adjustment.
 Υ represents a risk-weighted fraction of total PWID that are infectious and is given as Here the numerator is the risk-weighted total number of infected PWID (d = 6 … 12) and the denominator is the risk-weighted total number of people, whether they are infected or not. Note that the risk-weighting ensure that only active injectors are included in these totals, since , , = 0 for inactive injectors (n=4).
The intuition behind this model is that PWID interact (e.g., by sharing needles) with other PWID with the riskweighting terms reflecting the frequency of interaction and likelihood for the interaction to result in transmission of HCV.
The intrinsic risk is assumed to be different before and after 2013 to capture the change in PWID behavior with the onset of the opioid crisis. See Model Calibration below for more information.
The assumed parameter values are as follows. Here LN(μ, σ) denotes a log-normal distribution with mean μ and standard deviation σ.
Liver Disease Progression (DS) Figure 1A in the main paper depicts the structure of the liver disease progression submodel. The disease progression rates depend on HCV infection status. For infected people (d=6, …, 12), the disease progression (DS) terms in the model are Here the are disease progression rates and the are increased mortality rates associated with advanced stages of liver disease.
Uninfected people (d=1, …, 5) do not progress through the first four liver disease states and have different rates of progression through the 5th and 6th stages of liver disease ( and instead of and ). For these people, the DS terms in the model are Our disease progression rates are based on Kabiri et al. and are converted from the annual transition probabilities given there to transition rates as used here. 40 Given an annual transition probability p, the corresponding transition rate is = − log(1 − ). The exception to this is the liver transplant to post-liver transplant transition ρ8 which is modeled as a deterministic process by Kaberi taking one year, by definition, since post-liver transplant is defined to be one year post-transplant. Thus we take ρ8 = 1.

Testing and Treatment (TT)
The structure of the testing and treatment submodel is depicted Figure 1B in main paper.
PWID (both susceptible and infected) who are not tested "age" from one testing category to the next (e.g., from < 6 months to 6-12 months since the last negative test result) until they reach the final category (> 24 months) where they remain until tested. The aging rate ν = 2.0 reflects the fact the time intervals are 6 months (=1/2 year) in length. PWID who are not infected and are tested revert to the "< 6 months" (d=1) since last negative category. Those PWID who are infected and are tested become linked to care (d=11). The testing rates Δ , , depend on enrollment in MAT/SSP (the i and j indices), the liver disease state (k index), and the time since last tested negative (captured through p index). The specific rate assumptions will vary by scenarios.
PWID who are linked to care exit via treatment at rate or are lost to follow up at rate , with the rates potentially depending on the liver disease state (the k index) and the injector status (the n index), to capture restrictions related to disease or sobriety in some scenarios.
PWID who are lost to follow up are reengaged at rate Δ , , , corresponding to the testing rate for the given category for a person who tested negatively (p=1) in the last 6 months; the assumption is that people get reengaged with care through a process that is similar to the process of becoming linked to care. PWID receiving treatment achieve a sustained viral response (SVR) at rate .
Combining the testing and treatment processes, the testing and treatment (TT) terms in the model are Testing Parameters. The assumed treatment parameters vary by timeframe and across scenarios. In all cases, the testing rates Δ , , were decomposed into factor-specific terms the form  is a multiplier that captures the idea that people who are asymptomatic for liver disease (i.e., in METAVIR stages F0, F1, F2, F3 corresponding k=1,2,3, 4) are less likely to be tested for HCV and those with more advanced liver disease are more likely to be tested.


is a multiplier that captures dependence on how long it has been since the person was last tested (p) and whether they are enrolled in MAT/SSP (i,j). With these assumptions, the asymptomatic testing rate varies from that of the Base Case Testing scenario to that of the Improved Testing scenario. The terms in these scenarios are as in the Base Case Testing scenario when there is community testing only and as in the Increased Testing scenario when there was annual testing at the MAT and SSP.
Treatment Parameters. The assumed treatment parameters vary by timeframe and across scenarios. In the base case scenario, the treatment rate reflects the time required to successfully complete treatment. The loss rates were selected to fit published data on attrition through the HCV care cascade. 8 Given treatment and loss rates and , a fraction = /( + ) of those linked to care will actually be treated. We chose these fractions and solved for the loss-to-follow-up rate = (1 − ) / . (We suppress the dependence of the treatment rates and fractions on the liver disease state k and injector status n in this discussion.) The assumed values vary by scenario as follows:

Model Calibration
To capture the effects of the opioid epidemic, the time horizon was divided into three time periods -before 2013, 2013-2018, and after 2018 -with the injector inflow rate varying across these time frames. The infectivity rate is different before and after 2013.
In the deterministic models, we assumed:  Before 2013, the population is assumed to be in steady-state and and were selected to give a prevalence of 40% among all PWID 18 and a population of 4000 (=8000/2) active injectors in 2013. The prevalence estimate represents a US average for PWID RNA prevalence, as no NH-specific prevalence data is available in that time period. The population estimate is based on an estimate of 8000 active injectors in 2018 (see below) and a two-fold increase during the opioid crisis. 1,37,41 In the base case, this calibration process results in  172 (i.e., 172 new injectors entering per year) and  0.069.
 From 2013 to 2022, the inflow rates before and after 2018 were selected to match a population of 8000 active injectors in 2018 and a 10% reduction of active injectors in 2022, that is, 7200 active injectors in 2022. We chose 8000 active injectors based on NSDUH estimates of past-year heroin, cocaine, and meth use, 1 and adjusted these estimates taking into account rates of injection use for each drug. [42][43][44] This method produced a peak population size of between 4917-13741 active injectors. NSDUH past-year heroin use and NH DMI overdose mortality data suggest peak use occurred between 2016 (NSDUH) and 2018 (NH DMI; see eFigure 1b). 1,2 Because NSDUH data may not capture some changes in substance use behavior (such as a switch from heroin to fentanyl), we assumed a peak population in 2018. The assumed reduction in active injectors is based on declining rates of reported past-year heroin use and drug overdose in NH, 1,2,44 though we assume a more modest reduction than what is seen in this data, based on reports of rising substance use during the COVID-19 pandemic. 45 The infectivity rate was selected to ensure a prevalence of 45% among all PWID in 2019 based on Romo et al. HCV seropositivity rates. 8 This calibration process results in  0.298 and  1363 before 2018 and 448 after 2018, reflecting an increase in both new injectors and infectivity associated with the opioid crisis.  The 2019 prevalence was taken to be normally distributed with mean 45% and standard deviation 5% and the 2013 prevalence as taken to be normally distributed with mean equal to 2019 prevalence less 5% and standard deviation 5%.

Model Validation
To validate our model, we compared model outcomes to published estimates not used in parameter selection and the calibration process. 34 These outcomes included HCV prevalence, HCV incidence, new chronic HCV diagnoses, HCV treatments, and HCV mortality as shown in eTable 1. Because PWID-specific data is scarce, some of the validation measures are not exclusive to PWID; these include the HCV prevalence, treatments, and mortality. We include these observations as they may still be useful in assessing the validity of the model forecasts.

Policy Sensitivity Analysis
In the six illustrative scenarios in the main paper, we considered various combinations of four dimensions of possible interventions -increased testing in the community (eg, primary care or emergency departments), scheduled testing in MAT/SSP, increased treatment uptake, and increased harm reduction coverage. To better understand the interactions among these dimensions, eFigure 10 shows the forecasted 2045 HCV prevalence and total HCV treatments associated with all possible combinations of interventions. In eFigure 10:  The outer x-axis considers increasing access to/enrollment in harm reduction programs. The base case is as described above (15% of active PWID in SSP and 15% in MAT), the aggressive case is that in Scenario 2 above (50% in SSP and 35% in MAT) and the moderate case is halfway between (32.5% in SSP and 25% in MAT).
 The outer y-axis considers increased treatment uptake. The base case assumes 20% treated, the aggressive case is that considered in Scenarios 4, 5, and 6 (100% treated), and the moderate case is halfway between (60% treated). Detailed descriptions of the treatment model parameters in these three scenarios are in the Treatment Model discussion above.
 The x-axis of each bar chart shows the community testing rate. The red and blue bars show the 2045 prevalence with and without (respectively) annual testing at MAT/NSPs. In the scenarios without annual testing at MAT/NSPs, everybody is tested at the community rate. The community rate (1 per 20 years) is as in the base case and the highest (1 per 2 years) is as in Scenario 3. Detailed descriptions of these testing model parameters in these scenarios are in the Testing Model discussion above.
 The numbers above the prevalence bars show the total number of HCV treatments from 2022 to 2045 for each scenario.
In eFigure 10, we see that improvements in any of the four intervention dimensions were associated with a decreased 2045 prevalence. The lowest 2045 prevalence is achieved by pursuing all interventions aggressively, which corresponds to Scenario 6 in the main paper.
We also observe that the improvements associated with the interventions are synergistic. For example, the improvements associated with annual testing for those enrolled in harm reduction programs are indicated by the difference in heights for the blue and red bars in the graph. These differences increase (in relative terms) with increases in the treatment uptake, harm reduction, and testing rate. Thus the reduction in prevalence associated with an intervention (in this example, annual testing) is greater when pursued with other interventions.
Moreover, among the scenarios that achieve significant reductions in prevalence, the total number of treatments required decreases as the intensity of the intervention(s) increases. For example, the scenario with the maximum improvement on each dimension (Scenario 6) achieves near-zero prevalence with fewer treatments than any other combination of interventions associated with a 2045 prevalence to 30% or less. As discussed in the main paper, aggressive treatment early in the forecast period decreases the total number of treatments required over the remainder of the forecast period.

Monte Carlo Simulation Results
To test the robustness of our results, we conducted a Monte Carlo simulation analysis where we randomly generated 500 sets of model parameters according to the distributions described above. We generated model forecasts for the six scenarios described in the main paper for each set of model parameters. The results of the simulation are summarized in eFigures 11-16, with one figure for each of the six scenarios. The figures show the 10 th , 50 th , and 90 th percentiles for key model forecasts from 2022 to 2045: prevalence, incidence (HCV infections and reinfections), HCV treatments and liver deaths.
The results of the simulation demonstrate the robustness of the conclusions from the deterministic analysis discussed in the main paper. Though there is considerable uncertainty in some outcomes, the two scenarios with improved testing and treatment (Scenarios 5 and 6) are the only ones that are consistently associated with a near-zero prevalence across this broad range of scenarios.
Scenario 6 (with increased harm reduction, testing, and treatment) is particularly robust in that model forecasts consistently show dramatic, rapid reductions in prevalence and incidence of HCV. The final prevalence ranges from 0.0% to 0.3% (10 th and 90 th percentiles) in the scenario. There is some uncertainty in the number of treatments required in this scenario, reflecting the uncertainty about the number infected in 2022 (which in turn affects the estimated intrinsic rate of infection). The testing and treatment resources in this scenario are adequate to handle this broad range of possible values. As in the deterministic analysis discussed in the paper, the number of treatments associated with Scenario 6 is less than Scenario 5.
In all scenarios, there is considerable uncertainty about the number of infections, treatments and liver deaths: this uncertainty reflects uncertainty about the number of PWID and the number infected in 2022, and for deaths, their state of liver disease, as well as uncertainty about the death rates and transplant rates. Scenarios 1-4 have consistently higher and sometimes increasing prevalences and rates of infections and deaths than Scenarios 5 and 6.

Treatment Rate Sensitivity Analysis
We also performed sensitivity analyses on Scenario 6 where we limit the number of annual treatments. The treatment rate in Scenario 6 is briefly over 3000 treatments per year in the first year of the intervention and it is natural to wonder what would happen if such a treatment rate were unattainable, for budget, supply, or other capacity limitations. To study this, we limit the treatment rate, considering limits of 500, 1000, 1500, 2000, 2500, 3000, and Infinity (i.e., no limit) treatments per year. In this analysis, when the demand for treatment exceeds the limit, treatments were prorated among those who would have been treated if there were no limits.
The results of this sensitivity analysis are summarized in eFigure 16. Here we see that with all of these limits, the prevalence is eventually near zero, but scenarios with lower limits reach this point later and were associated with more treatments in total. Placing a cap on treatments creates a bottleneck in the linked-to-care state as infected PWID are rapidly diagnosed, but cannot be treated as fast as they are identified. The larger limits considered (2000-3000 treatments per year) have little effect on the number of treatments required and the number of infections because these limits are reached only briefly. The lower treatment limits have increasingly severe negative impacts on all of the metrics considered.