Comparing Projected Fatal Overdose Outcomes and Costs of Strategies to Expand Community-Based Distribution of Naloxone in Rhode Island

Key Points Question What community-based naloxone distribution strategies would be the most effective and efficient in preventing opioid overdose deaths? Findings In this decision analytical model study evaluating the distribution of 10 000 additional naloxone kits annually in Rhode Island, the strategy focusing on distribution of naloxone according to geographic need to people who inject drugs resulted in the best outcomes at the lowest cost, averting an estimated 25.3% of opioid overdose deaths at an incremental cost of $27 312 per opioid overdose death averted. Meaning This study suggests that expanding naloxone distribution to people at highest risk for opioid overdose should be prioritized and that redirecting spatial distribution of naloxone to areas with the greatest need will improve both effectiveness and efficiency and reduce geospatial health inequality.

This supplementary material has been provided by the authors to give readers additional information about their work.
The PROFOUND model is an individual-based model composed of a microsimulation and a decision tree model. The individual-based microsimulation allows the inclusion of stochastic variation as well as variation due to individual characteristics, capturing population heterogeneities. 1 It also overcomes the Markovian assumption required by Markov models by allowing for "memory" in the model's structure, such as the individual's history of overdose. The microsimulation model is used to track changes in health and drug use states and to forecast overdose events among simulated individuals. The decision tree is used to determine the potential pathway and consequence of each overdose event.

Study population
We initialized the simulated study population representing all individuals in Rhode Island age 12 years and older who are at risk for opioid overdose, including people with exclusive prescription opioid misuse (excluding those who also use illicit opioids), any non-injection illicit opioid use (excluding those who also inject illicit opioids), any illicit opioid use via injection, and people who exclusively use stimulants without intended opioid use (but whose substances may be contaminated by fentanyl, see eFigure 1). We also included people with an opioid use history but who are currently not actively using opioids. In addition, we characterized the simulated study population by sex (male/female), age (continuous), race/ethnicity (Black/African American [Black], Hispanic/Latino [Hispanic], non-Hispanic white [White], and other [Other]), city/town of residence (39 in total), overdose history (yes/no), and fentanyl exposure (yes/no). eFigure 1. Drug Use State Stratification * Whose drugs may be contaminated with fentanyl; ** exclusive prescription opioid misuse.
In this study, we defined city/town as the subregion in Rhode Island (eFigure 2). In the study population initialization for the microsimulation, we first determined the population size and size of each demographic group (race * age group * sex) in each city/town based on the American Community Survey (2010 Census) for Rhode Island. 2 Using the National Survey on Drug Use and Health (NSDUH) data (northeast region), 3 we also estimated the prevalence of opioid misuse (all types) and exclusive stimulant use (cocaine and methamphetamine) among each demographic group. We then combined results from these two datasets to estimate the opioid misuse populations within each demographic group in each city/town. We compared the resulting estimated total number of people with opioid misuse with the Rhode Island's overall estimate, based on a statewide estimate of the prevalence of opioid use disorder for Rhode Island (5.2%) 4 and the state population, and adjusted the demographic and jurisdiction-specific estimates by a multiplier that reflects the difference between the two statewide estimates. We determined the population who exclusively use stimulants (without intended opioid use) within each demographic group in each city/town using the NSDUH data without adjustment (eTable 1).

eFigure 2. Map of Cities/Towns in Rhode Island
We then used the NSDUH dataset to derive, among the people with opioid misuse of different sex, the proportion of each type of opioid use. The variable names we used in the estimation process with the NSDUH data are presented in eTable 2. We also used a statewide cross-sectional assessment of the cascade of care for opioid use disorder in Rhode Island 4 to determine the proportion of the simulated people with an opioid use history but who are currently not actively using opioids (defined as in recovery from opioid use disorder We then assigned the initial overdose history variable value (yes/no) to simulated study individuals, stratified by type of drug use based on evidence from literature (eTable 3). This variable is updated during the microsimulation when overdoses occur. The initial level of fentanyl exposure (intentional use of fentanyl or fentanyl-contaminated drugs) was also estimated from the literature. [5][6][7][8] We assumed a monthly increase of 0.5% [0% -1%] 9 to account for secular trends in which we assigned by increasing in each monthly time-step the proportion of individuals that were exposed to fentanyl. We assumed that fentanyl exposure is limited to people who use illicit opioids, stimulants and prescription opioids not from prescribers, whereas people who exclusively misuse prescription opioids sourced from prescribers (including those from friends/relatives who were prescribed these opioids) were not at risk for fentanyl exposure. Estimates for the source of prescription opioids were based on NSDUH data. eTable 3 presents the resulting list of parameters used in defining the study population.

Opioid overdose risk (microsimulation)
Naloxone is effective in reversing overdoses caused by opioids. We modeled opioid overdose only and assumed different monthly risks of overdose for different patterns of opioid use based on the literature (eTable 4). Overdose events in each monthly timestep were randomly drawn among simulated individuals according to these probabilities of experiencing overdose. In general, people who use illicit opioids face a higher risk of overdose than those with exclusive prescription opioid misuse and those who use injection illicit opioids have a higher risk of overdose than those who use noninjection illicit opioids. We accounted for several factors that can elevate risk of overdose: (1) fentanyl exposure; (2) overdose history, where risk for subsequent overdose was assumed to be higher than risk for the initial overdose; and (3) first month of opioid use after being in the inactive opioid use state (see Section 2.4 below), to account for elevated risk associated with decreased tolerance for opioids after a period of abstinence. We assumed that the risk of opioid overdose was zero for people in the inactive opioid use state and people who exclusively use stimulants but whose drugs were not contaminated by fentanyl. In the absence of clear evidence about the risk of opioid overdose among those whose stimulants were contaminated by fentanyl, we assumed this risk was equal to the risk for those whose non-injection illicit opioids contain fentanyl.
Simulated individuals who experienced opioid overdose in each monthly time-step were assigned subsequent outcomes using the decision tree (see Section 2.5 below). Relative risk of overdose for fentanyl exposure 6.07 3.63 -10.16 17,18 Relative risk of overdose during the first month of opioid relapse 4.3 3.6 -5.2 19 Relative risk of overdose for subsequent versus first overdose 3.5 1.9 -6.4 14

Transitions between health states (microsimulation)
Individuals' transitions between health states are determined using a random process drawn from a transition probability matrix (Table S5). Transition probabilities were estimated from the published literature.
If a simulated individual is in the exclusive prescription opioid misuse health state, in each monthly time-step this individual is subject to a monthly probability of transiting to non-injecting illicit opioid use; 20 individuals in the non-injection illicit opioid use health state can transition to illicit injection opioid use; 21 and individuals in the illicit injection opioid use health state can transition back to illicit non-injection opioid use 22 . Individuals in each of these three opioid use states can also transition between these opioid use states and an inactive opioid use state (due to treatment or recovery). [23][24][25] We also included as a separate health state for the first month of opioid use after being in the inactive state (relapse) 26 to account for the elevated risk of overdose, after which the simulated individuals would transit (with a transition probability of 1) to the same active opioid use state prior to transitioning to the inactive opioid use state. For simplicity, we did not account for potential transitions between opioid use states and stimulant use states, and we did not include inactive stimulant use state.
In addition, simulated individuals in all health states are subject to age-specific risk of death due to causes other than overdose ("background mortality"), which we calculated by converting published annual death rates (based on life tables 27 ) into monthly probabilities (eTable 6) and applying a multiplier to reflect higher risks among people who use drugs. 28  * Calculated using a multiplier for the mortality rate among people who use drugs compared to general population. Only the death rates among people who use drugs were used in the model.

Opioid overdose events (decision tree)
We used a decision tree module to determine the pathways and consequences for each overdose event. Model parameters for the decision tree are presented in eTable 7. In the decision tree, each node is associated with branches defined by a set of probabilities that add up to one. The probabilities vary depending on the characteristics of the simulated individuals and the outcomes from the previous nodes. These nodes include the setting of the overdose (public versus private/semi-private), whether the overdose is witnessed, whether naloxone is administered (among witnessed overdoses), whether emergency medical services (EMS) is dispatched as a result of a witness calling 911 (among witnessed overdoses), whether the overdosed individual is transported to the hospital for emergency department (ED) care (among those for whom EMS is called), and whether the individual dies or survives.
In the decision tree, the probabilities of an overdose being witnessed and of a witness calling for help from EMS depends on the setting of overdose (public versus private/semi-private). The probability of dying from an overdose depends on whether naloxone is administered and whether EMS is called. The baseline probability of dying from an overdose where no naloxone is administered nor EMS called is based on observational studies. 30,31 Unlike previous modeling studies, we did not assume the probability of dying when naloxone was administered to be zero; instead we derived this estimate from a systematic review of observational studies which used a pooled estimate of the proportion of opioid overdoses where naloxone kits were used that resulted in death. 32 This might result in a conservative estimate of the effectiveness of naloxone but this assumption is deemed warranted in the current era of fentanyl. We also included a lower relative risk of death from overdose when EMS is called but no naloxone is administered by a witness. 33,34 If the simulated individual survives, this individual then returns to the microsimulation either in the inactive opioid use health state 35

Naloxone availability algorithm
To estimate the probability of naloxone being available and administered during a witnessed overdose, we adopted a previously used approach by Irvine et al. 39 assuming the probability is a nonlinear function of the number of naloxone kits in circulation and the size of the population at risk for opioid overdose. These variables can vary by city/town: where denotes the probability that naloxone is available and administered during a witnessed overdose in region r at time t; denotes the cap on naloxone availability, which was assumed to be 0.99; denotes the adjustment factor for naloxone availability, whose value is determined in model calibration (see Section 3.2); denotes the number of naloxone kits from OEND (community-based) programs in circulation in region r at time t; denotes the number of naloxone kits from pharmacies in circulation in region r at time t; denotes the ratio of effectiveness of pharmacy programs in reaching at-risk population compared to OEND programs (0.371 [0.345, 0.4], based on an observational study for the correlates of naloxone recipient characteristics with distribution program characterisics 40 ); denotes the number of at-risk (for opioid overdose) individuals residing in region r at time t based on model estimates; denotes the monthly rate of naloxone withdrawn due to expiration/loss (1/15.5 months, i.e., 6.5%), based on an unpublished analysis of the New York City Department of Health and Mental Hygiene Naloxone Recipient Form data for the average circulation time for naloxone kits; denotes the number of naloxone kits from OEND (community-based) programs received by residents in region r in a given year; denotes the number of naloxone kits from pharmacies received by residents in region r in a given year. Naloxone data for the number of naloxone kits distributed by OEND programs and received by residents of each subregion (i.e., city/town) were derived from data collected by the Rhode Island Department of Health (RIDOH) using a statewide, standardized reporting form for each individual receiving naloxone kit(s) from any OEND programs. Because subregional level data for naloxone distribution from pharmacies were unavailable in Rhode Island, we used the aggregate annual number of naloxone kits distributed by pharmacies in the state and assumed the amount received was proportional to the size of at-risk population in each subregion.

Population dynamics
When simulated individuals leave the model due to age-stratified background mortality (eTable 6) or overdose death (eTable 7), the deceased individuals were replaced by new ones with the same initial characteristics except for overdose history (reset as 0). This allowed us to maintain a fixed size of the simulated population of people at risk of overdose in Rhode Island. For simplicity we did not include in-migration/out-migration at the state or city/town level.

Model calibration
Model calibration refers to the process of matching model outcomes with observed data by adjusting uncertain model parameters and establishing plausible ranges that provide the best fit to available data. 41 Calibrating model inputs to observed epidemiological endpoints ensures the credibility of model results and thus strengthens our confidence in model inferences. We used a two-step calibration procedure for the PROFOUND model. First, we conducted an initial calibration of a smaller set of model parameters associated with overdose setting that are only used in the decision tree model. We then conducted a formal calibration for the entire model.

Initial calibration to adjust for overdose setting in the decision tree
During data collection to inform our model, we identified substantial differences in the setting of fatal opioid overdoses comparing data from the State Unintentional Drug Overdose Reporting Surveillance 42 and Rhode Island Office of the State Medical Examiners 43 versus data of non-fatal opioid overdose collected by the Rhode Island Emergency Medical Services (EMS) Information System. 44 In the first source approximately 10% of fatal opioid overdoses were reported to occur in public settings, as compared to approximately 31% of non-fatal opioid overdose reported in EMS data. This difference is likely attributable to the overrepresentation of opioid overdoses occurring in public captured by the EMS system and/or higher survival rates from overdoses occurring in public as a result of higher likelihood of being intervened and rescued. To address these differences, we performed an initial calibration for three parameters: (1) proportion of overdoses occurring in public settings; (2) relative risk of overdose being witnessed in private/semi-private versus public settings; and (3) relative risk of witness(es) calling EMS in private/semi-private versus public settings. The calibration was compared against the two targets for opioid overdose setting derived from fatal overdose surveillance data and EMS data (for non-fatal overdoses). The prior values and distributions for the calibration parameters and targets are presented in eTable 8.
In this initial calibration, we used a Bayesian method with the incremental mixture importance sampling (IMIS) algorithm 45 (using R package IMIS 46 ). Given the challenges in calibrating a stochastic model (due to Monte Carlo error), 45 we first transformed the decision tree model from individual-based to cohort-based (i.e., from stochastic to deterministic) assuming a cohort of 10,000 individuals experiencing overdose. All other decision tree parameters were held fixed at their prior point estimates, since they were independently defined from the calibration targets and had low influence on the calibration targets. 45 We generated a posterior sample of one million parameter sets from the IMIS algorithm. The results, including the posterior distribution of calibrated parameters and model fit to calibration targets, are presented in eFigure 3.

Formal calibration for the entire model
In the formal calibration for the entire model, we used a random calibration approach to repeatedly sample from estimated uncertainty ranges for 16 key model parameters that had the greatest influence on target estimates (eTable 9). We compared model projections against three targets in each year between 2016 to 2019 in Rhode Island reported by the Rhode Island Department of Health: (1) the number of opioid overdose deaths (OODs); (2) the percentage of OODs involving fentanyl; and (3) the number of emergency department visits related to opioid overdose (eTable 9). We used a Latin hypercube sampling method 47 to draw one million random parameters sets from the parameters' uncertainty range, augmented them with the one million sets of overdose setting parameters from the initial calibration and one million random seeds (to ensure reproducibility of results), simulated the model with these sets of inputs and seeds, and compared the resulting model projections against the selected calibration targets. We identified 500 calibrated subsets (and seeds) providing the best goodness-of-fit statistics for subsequent analysis. The goodness-of-fit (GoF) was measured by the mean percentage deviation, as shown in the equation below: where N represents the number of calibration targets, represents the modelprojected result for the ℎ target, and represents the observed estimate for the ℎ target. Smaller values of the GoF indicate a better fit to the observed data.
This approach, together with the stochastic process embedded within the microsimulation design, allowed us to derive calibrated model parameters providing  44 good fit to observed targets while simultaneously exploring the uncertainty of model outcomes resulting from both parameter uncertainty (the uncertainty in estimation of the parameter of interest) and stochastic uncertainty (random variability in outcomes between identical simulated agents). 41 We present in eTable 9 the posterior value of parameters after this calibration process.

Model validation
Model validation refers to the process of evaluating a model's accuracy in making relevant projections. 48 In particular, external validation entails the comparison of model projections to external estimates of key clinical and epidemiological data not used in the model. Given that our model was calibrated to epidemiological targets at the state level but our model was built to replicate the opioid overdose epidemics in each city/town, we compared our projections for the annual number of OODs in each city/town from the 500 calibrated parameter sets with Rhode Island Department of Health surveillance data at the city/town level from 2019 49 to examine whether the observed numbers fall within the 95% simulation intervals of the model estimates in each jurisdiction (eFigure 4). The results show that the surveillance data for all cities/towns but one (Lincoln) fell within the 95% simulation intervals of the model estimates.

Costs
In this study, we included two components of costs associated with the naloxone distribution expansion strategies: cost of naloxone kit, and cost of naloxone distribution. Given the short timeframe of evaluation (three years), all cost estimates were undiscounted.

Cost of naloxone kit
The cost of naloxone kit refers to the purchase cost of each naloxone kit (two doses) by the health department. There currently two types of naloxone kits distributed in Rhode Island, intramuscular versus intranasal, each with a disparate price. We used the RIDOH data to estimate an average naloxone kit cost based on the proportions of the two types of naloxone kit distributed during January 2020 to September 2021 by all OEND programs and the unit purchase price to RIDOH of each kit type, assuming the same composition of naloxone kits in the expanded naloxone distribution (eTable 10).

Cost of naloxone distribution
We estimated the cost of providing overdose education and naloxone distribution (OEND) which includes costs related to staff time spent receiving OEND training and delivering OEND and excludes naloxone medication costs. These costs were derived from OEND costs estimated in New York City. 50 We adjusted these costs by applying Bureau of Statistics national wage rates in place of NYC wage rates and removing NYC-specific program costs that would not be relevant for Rhode Island OEND programs, including time to assemble New York state-provided OEND bags by filling them with naloxone kits, educational materials, and breathing masks. After these adjustments, OEND costs per naloxone kit distributed were estimated at $22 per kit for syringe services programs (SSPs), $162 for single site non-SSPs, and $40 for multisite non-SSPs.
We used these results to assign costs to each of the four naloxone expansion strategies based on their exemplary OEND programs operating in Rhode Island: SSPs, street outreach programs (using a weighted average of single-site and multi-site non-SSPs based on the number of kits distributed by each type of program), community events (single site non-SSPs), and primary healthcare settings (using a weighted average of single-site and multi-site non-SSPs based on the number of kits distributed by each type of program

Summary inequality measures
We used the Theil index and between-group variance statistic 51 as summary measures of geospatial health inequality using the differences in the simulated rates of OOD across cities/towns under each scenario. The calculation equations and characteristics of each measure are described in eTable 12.