Prevention of Prescription Opioid Misuse and Projected Overdose Deaths in the United States

Key Points Question What is the projected effect of lowering incident nonmedical prescription opioid use on the future trajectory of the opioid overdose crisis in the United States? Findings In this system dynamics model study, under current conditions, the opioid overdose crisis is expected to worsen—with the annual number of opioid overdose deaths projected to reach nearly 82 000 by 2025, resulting in approximately 700 000 deaths from 2016 to 2025. Interventions focused on lowering the incidence of prescription opioid misuse were projected to result in a 3.0% to 5.3% decrease in opioid overdose deaths over this period. Meaning Prevention of prescription opioid misuse alone is projected to have a modest effect on lowering opioid overdose deaths in the near future, and multipronged approach is needed to dramatically change the course of the epidemic.


eAppendix 1. Details of System Dynamics Model
We developed a system dynamics model (also known as compartmental model) (1), the Opioid Policy Model (OPyM), to simulate the opioid overdose crisis in the United States (US) from 2002-2025, and used it to assess the effects of reducing the incidence of nonmedical opioid analgesic use on the number of opioid overdose deaths.
We estimated the values of model parameters directly from data, where possible, including the National Survey on Drug Use and Health (NSDUH) (2), and the multiple cause of death data from the Centers for Disease Control and Prevention (CDC) (3). For parameters that could not be informed by available data sources, we fitted them using a model calibration approach. Further details of parameter estimation and calibration are described in eAppendix 2.
In this section, we describe the definition, structure, and evaluation of the system dynamic model in further detail, specifically: 1) Definition of compartments and the mathematical system dynamics model; 2) Formulation of time-dependent model parameters using joinpoint regression analysis; 3) Model evaluation accounting for uncertainty in input parameters.

A1.1 System Dynamics Model Structure
The system dynamics model consists of three compartments that represent the major subgroups of the population that use opioids (prescription and illicit) non-medically. The definitions of these population groups were adapted from the available information in the NSDUH survey data. In particular, we categorized the population with non-medical opioid use into the following three groups (eFigure 1): 1) People with non-medical use of prescription opioids, denoted by compartment (N). The corresponding selection criteria for this group in the NSDUH data include all the respondents who claimed to use opioids non-medically in the past year, did not claim heroin use in the past year, and did not meet opioid use disorder DSM-IV (Diagnostic and Statistical Manual of Mental Disorders 4th Edition) diagnosis criteria; 2) People with prescription opioid use disorder (OUD), denoted by compartment (D). The corresponding selection criteria for this group in the NSDUH data include respondents who met DSM-IV criteria for opioid abuse or dependence and did not claim to use heroin in the past year; 3) People with illicit opioid use, denoted by compartment (I). The corresponding selection criteria for this group in the NSDUH data include respondents who claimed to use illicit heroin/fentanyl (regardless of the route of administration). a. It is worth noting that it is possible that an individual uses both prescription opioids and heroin at the same time in real-world situations. However, by our hierarchical definitions of compartments, if the person claimed to use heroin/fentanyl in the past-year, he/she is categorized into the illicit opioid use compartment (I).
b. To keep the model parsimonious, we did not separate illicit opioid use without vs. with OUD. Separating it into different compartments would have increased the number of model parameters whose estimates are not readily available (e.g., transition rate from illicit opioid use to illicit OUD), without allowing us to query additional policy implications that were relevant the scope of this study.
Time progresses continuously in the model. As shown in eFigure 1, new individuals enter the model using either prescription opioids non-medically (compartment N) or illicit opioids (compartment I), and then transition through different compartments, i.e., progressing through different states of opioid use. From non-medical use of prescription opioids, individuals can develop a prescription OUD (compartment D). Individuals who used prescription opioids, with or without prescription OUD (compartments N and D), can transition to illicit opioid use (compartment I). Individuals were subject to opioid overdose death with mortality rates dependent on their compartment. In addition, individuals could transition out of the model when they either stopped using opioids or died from other (i.e., non-opioid-related) causes.
Mathematically, we used ( ), ( ), ( ) to represent the prevalence for each compartment at time , and the system dynamics were represented by the following ordinary differential equations: Annual incidence of non-medical use of prescription opioids and of illicit opioid use from other sources (i.e., not from prescription opioids), respectively, at time t. The model differentiates two cases of incident illicit opioid use: those who progress to illicit opioid use from prescription opioid use and those who initiate opioid use with illicit opioids; ( ) refers to the latter case. The system dynamics model was implemented in R programming language (Version 3.4.0) using "deSolve" package to solve the ordinary differential equations.
For some model parameters, we considered time-dependent structure in their values, because a constant value over time may not be sufficient to capture the highly dynamic trends in the current opioid epidemic, such as rapidly increasing overdose deaths in the past few years, and the increasing proportion of opioid initiates from illicit opioids. Next, we describe the formulation of time-dependent parameters.

A1.2 Joinpoint Model for Time-Dependent Parameters
To model the time-dependent structure of model parameters, we applied joinpoint regression analysis (4,5), which is commonly used to describe trends over time and identify changes in a trend at certain time points. A joinpoint regression model for parameter ( ) consists of the following components: 1) 0 , the baseline value at the initial time point 0 ; 2) 1 , annual percentage of change (APC) before the joinpoint; 3) , a joinpoint time indicating a significant change of APC from time onwards (i.e., an inflection point of trend at time ); and 4) 2 , APC after the joinpoint . Mathematically, time-dependent parameter ( ) is determined as follows: If there is no joinpoint specified, ( ) increases/decreases with APC 1 indefinitely since time 0 . Joinpoint regression can also be generalized to include more than one joinpoint.
For our system dynamics model, we first estimated the annual incidence of prescription opioid misuse from 2002-2014 using the NSDUH data (see A2.1 for estimation details), and then applied the joinpoint regression analysis on the annual incidence to obtain the joinpoint equation for ( ) (eFigure 2A, eTable 1). Instead of directly using the raw annual incidence estimates from NSDUH data, we used the joinpoint regression equation in our system dynamics model to avoid unsmooth variations in the incidence values across years.
To estimate time-dependent values of the overdose mortality rates, ( ) and ( ), we performed the joinpoint regression analysis on the multiple cause of death data from the CDC Wide-Ranging Online Data for Epidemiologic Research (WONDER) for 2002-2015, and identified inflection points in the numbers of annual overdose deaths ( Figures  2B and 2C). In particular, we incorporated the joinpoint year of 2011 for overdose mortality rate from illicit opioids ( ), and a joinpoint year of 2007 for overdose mortality rate from prescription OUD ( ) (eTable 1). However, the baseline value, 1 , and 2 obtained from the joinpoint regression equations cannot be directly used for parameters ( ) and ( ), because the observed data from CDC WONDER were the counts of overdose deaths, whereas the parameters ( ) and ( ) represented the transition rates per year (i.e., annual risks) of overdose deaths. Thus, in the joinpoint regression equations of the mortality rate parameters ( ) and ( ), we implemented the joinpoint years informed by the joinpoint analysis of CDC WONDER data, but determined the baseline value, 1 , and 2 via model calibration (see eAppendix 2).
For the incident illicit opioid use ( ( )), transition rate from prescription OUD to illicit opioid use ( ( )), and exit rate from compartments D ( ( )), we considered a simpler time-dependent model structure without a joinpoint, i.e., a simple growth models with fixed 1 , as we found that such structures provided better fit to the observed data (i.e., model calibration) compared with assuming constant transition rates. Due to lack of information for inferring the joinpoint time, we kept the time-dependent model simplistic, and their baseline values and APC 1 were determined through model calibration (eAppendix 2).

eAppendix 2. Model Parameter Estimation and Calibration
In this section, we describe the estimation and calibration of parameter values in our system dynamics model. To account for model uncertainty, we also describe the sampling distributions of input parameter values. Different input parameter values were sampled for each model evaluation, which was repeated by 1000 times. From the 1000 replications of model evaluations, we presented the mean results (i.e., the average value of model outputs), and presented the 95% uncertainty interval using the 2.5 th -97.5 th percentile range of model outputs.

A2.1. Parameter Estimation
For the following parameters, we estimated their values directly from the NSDUH data: 1) Initial prevalence values, 0 , 0 , and 0 for each compartment at year 2002. We first performed joinpoint regression analysis on the prevalence of each compartment estimated from the NSDUH data. The fitted joinpoint regression equations provided the estimates of the expected counts at each year, and we used the estimates at the year 2002 for the initial prevalence values. To account for uncertainty in those values, we used the standard error estimates from the NSDUH data and sampled the initial prevalence values assuming normal distribution as follows: • Initial value of compartment N: 0 with mean 10,029,859 and standard deviation 329,807; • Initial value of compartment D: 0 with mean 1,369,218 and standard deviation 116,347; • Initial value of compartment I: 0 with mean 328,731 and standard deviation 60,446.
2) Incidence of prescription opioid misuse ( ) at time . We first estimated the annual incidence of non-medical use of prescription opioids using the NSDUH survey data from 2003-2014. To capture the annual incidences for the full year, we calculated the incidence for each calendar year using the survey data from the following year. For example, incidence estimate for 2002 was calculated from survey data in 2003. Specifically, we queried for respondents of the survey in 2003 who claimed having initiated opioids in the previous year, i.e., in 2002. This way, it allowed us to capture an entire 12-month period for incidence (otherwise, using the 2002 survey data would underestimate the incidence for the 2002 year, because the survey could take place anywhere during the year, leaving the new incidences after the survey excluded). The survey data from 2015 was not used to calculate the 2014 incidence because the survey questions changed and the variables available were no longer synonymous. We used survey weights to estimate the incidence of non-medical use of prescription opioids for the entire US population. As indicated in Section A1.2, we then conducted the joinpoint analysis on the above estimates of annual incidence of prescription opioid misuse from the NSDUH data, which provided the estimates and the standard error of baseline value, 1 and 2 of the joinpoint regression equation (see eTable 1). We sampled these three values from a normal distribution based on the estimates of mean and standard deviation to determine the equation of ( ).

A2.2. Model Calibration
Because some model parameters are not directly observable from existing datasets, we used calibration methods (6)(7)(8) to estimate their values and ranges such that the model outcomes would closely match existing data. In this section, we describe the calibration parameters, targets, goodness-of-fit measure, search algorithm, and calibration results.

Calibration parameters
A total of 17 model parameters were determined by the model calibration, which are defined below: • : overdose mortality from compartment N (1 parameter);

Calibration targets
We compared the following model-predicted outcomes with observed data: prevalence of prescription opioid non-medical use, OUD, and illicit opioid use based on NSDUH data; (2) the number of overdoses from all opioids and illicit opioid use from on CDC WONDER data; (3) and the percentage of illicit opioid among initiating opioids based on a published study (eTable 2) (9).

Goodness-of-fit metric
For any given set of model parameters, we assessed how well the model was fitted to the calibration targets. That is, we compared the model outputs with corresponding calibration targets (eTable 3), and calculated the goodness-of-fit measure. For each calibration target, model error is defined by the sum of squared errors between model outputs and calibration targets over all years, and then normalized by the initial value of this target to preserve comparable magnitude across different calibration targets. The overall goodness-of-fit measure, the total error, was defined as the summation of model errors over all calibration targets (with equal weights):

Calibration procedure
To calibrate the model, we used a directed search method to identify the sets of parameters that lead to "good" model fits. An optimization algorithm was applied to iteratively search for the parameter sets with the lowest possible total model errors. We systematically tested multiple optimization functions in R (L-BFGS algorithm with package "lbfgsb3", Nelder-Mead algorithm with package "dfoptim", Differential Evolution Optimization with package "DEoptim", and Generalized Simulated Annealing with package "GenSA"), and selected package "GenSA", which consistently outperformed other algorithms with good solution qualities within reasonable computation time.
Search for the "optimal" values of calibration parameters was bounded within the predetermined range for each parameter. Because the search results from the optimization function (i.e., the calibration results) depend on the starting values of the search, we sampled the starting values for all parameters uniformly within their pre-determined ranges using Latin hypercube sampling approach and repeated the directed search for 1,000 replications. The parameter ranges were initially determined based on expert opinions and set conservatively wide.
For some parameters, we observed that the calibrated values were concentrated at the lower or upper bounds of their range. When this happened, we expanded the ranges and repeated the calibration process. We repeated this process until we observed that all calibrated parameters represented a "full" distribution of values (not a "truncated" distribution) over their search ranges. In addition, we made sure that the total calibration errors were reduced after each iteration of adjusting the search ranges. In this way, we believe that the search has adequately explored the promising areas with "good" values of model parameters that result in good model fitting.

Calibration results
The optimal solutions identified from 1,000 replications of model calibrations were collected and considered as the set of calibrated model parameters. To account for uncertainty in model input parameters, we did not average these parameter values into aggregated estimates for the model evaluation; instead, we evaluated the model outcomes using each of the 1,000 parameter sets and estimated 95% uncertainty intervals from the 1,000 replications of model evaluations. eFigure 3 shows the distribution of calibrated parameter values over the search ranges for the calibration, and eTable 4 summarizes values of the calibrated parameters. Model outcomes with calibrated parameters and comparisons with calibration targets can be found in Figures 2 and 3 of the main manuscript.

A3.1. Parameters for Projection Scenarios
To project model outcomes after year 2015, we considered three projection scenarios where the overdose mortality ( ) and annual incidences ( ) of illicit opioid use will continue to increase after 2015 but stabilize at different time points in future. In particular, we simulated: 1) a base-case scenario, which assumes that the overdose mortality rate (i.e., lethality) of illicit opioids and the annual rate of incident illicit opioid use would stabilize by 2020; and 2) a pessimistic scenario which assumes that both rates would stabilize by 2025.
To implement such a structure, for any given APC value (i.e., 2 for ( ) and 1 for ( ) that was obtained from model calibration), we decreased the value by a fixed amount every year after 2015 such that the APC decreases to 0 at the targeted stabilizing year (i.e., the year 2020 and 2025 in the base-case and pessimistic scenario, respectively); and the values of ( ) and ( ) remain constant afterwards. eFigure 4 illustrates the time-dependent structure of these two parameters for the base-case projection scenario.

A3.2 Simulated Strategies to Reduce the Incidence of Non-Medical Prescription Opioid Use
To evaluate the impact of reducing the incidence of non-medical opioid analgesic use on the number of opioid overdose deaths, we simulated the following four incidence trends representing different effects of strategies to reduce the incidence of non-medical prescription opioid use: (1) no change in the annual incidence of prescription opioid misuse from 2015 onward; (2) decreasing incidence of prescription opioid misuse from 2016-2025, at a rate of 7.5% per year as observed between 2008-2015 (i.e., 2 in the joinpoint model for ( ) before 2015), which may be achieved by continued implementation and success of prevention programs; (3) decreasing incidence of prescription opioid misuse at a rate that is 50% faster than strategy 2 (i.e., a 11.3% decrease per year), which assumes greater success of prevention programs; and (4) a hypothetical setting of no new incidence of prescription opioid misuse after 2015, which was included to assess the maximum possible benefits of prevention interventions under ideal conditions. eFigure 5 illustrates the retrospective incidence data prior to 2015 estimated from the NSDUH data, and four incidence trends from 2016 onwards that correspond with the above four prevention strategies (shaded areas represent uncertainty intervals).

eAppendix 4. Additional Sensitivity Analysis for NSDUH Estimates
One of our model limitations is that we used NSDUH data for model calibration and parameterization, which does not capture homeless and incarcerated population and is based on self-reported outcomes. Thus, the NSDUH data could represent underestimates of the actual epidemic situation. Since no other data are readily available for such underrepresented population at national level, we conducted sensitivity analysis on our model parameters estimated from NSDUH data to assess the impact of the underestimates on our model results.
To determine the possible variability to the estimates obtained from NSDUH data and to adjust the parameter values used in our analysis, we searched the literature to compare the differences between the estimates from NSDUH and other independent data sources. We found that the NSUDH data estimated the prevalence of opioid abuse and dependence in Massachusetts to be 1.17% between 2009-2011 (10), while the estimate was 2.72-2.87% in 2011-2012 based on the Chapter 55 dataset that linked the Massachusetts All-Payer Claims Database and other state databases (11). Therefore, we doubled our estimates of prevalence for each compartment and incidence of prescription opioid non-medical use from NSDUH data as an approximation for the adjustment. Then we re-calibrated our model to the adjusted prevalence and incidence estimates (eFigure 6) and repeated our analysis.
We found that, with the doubled prevalence and incidence estimates, the projected number of overdose deaths in the base case was increased by only 5.5% (739,100 compared with 700,400 in the base case results, eTable 5, eFigure 7). Such a small change is mostly because the model was re-calibrated to the overdose deaths data from CDC WONDER. Furthermore, the relative reduction in cumulative overdose deaths from 2015-2025 attributable to reducing incidence of prescription opioid misuse became 3.7%-5.2% (eTable 5), which is comparable to the 3.8%-5.3% reduction in the base case results. We acknowledge that our adjustment to the NSDUH estimates is not exact, but this sensitivity analysis demonstrated that our model results are robust, and the limitation of underestimation by the NSDUH data would not substantially change our conclusions.  Lines represent the average of 1000 outcomes from the model. Error bars represent 95% confidence intervals of the observed outcomes from the NSDUH data, and blue shaded regions represent the bootstrapped 95% uncertainty intervals of the model outcomes.
Abbreviations: NSDUH, National Survey on Drug Use and Health. Cicero (2017) (9) refers to the source of calibration targets.  ). * The joinpoint analysis provides SE for the slope, where APC = exp(slope)-1. Accordingly, slope for 1 is 0.00165, and slope for 2 is -0.0776. † No time-dependent structure was assumed for (the overdose mortality rate from compartment N). Considering the population without the diagnosis of opioid use disorder in compartment N, the rate is expected to be low and the temporal trend is not needed. For model simplicity, we considered as a constant parameter with respect to time. And thus, we assumed that the inflection point of the trends in overdose deaths from all opioids at year 2007 was attributed to the change in overdose mortality from compartment D, i.e., ( ).  Prevalence of prescription OUD ( ) 3 Prevalence of illicit opioid use ( ) 4 Overdose death from illicit opioid use ( ) ( ) 5 Overdose death from all opioids (from all compartments) Percentage of illicit opioid use at initial regular use of opioid ( ) + ( ) ( ) ( ) + ( ) + ( )