Curbing the Opioid Crisis: Optimal Dynamic Policies for Preventive and Mitigating Interventions

1. Introduction

The U.S. Department of Health and Human Services (HHS) declared a public health emergency in 2017 to tackle the nationwide opioid crisis (HHS 2017). This state of emergency remains in effect to this day. The opioid overdose crisis unfolded in successive waves: (1) fatalities related to prescription opioids, (2) deaths linked to heroin, (3) fatalities associated with synthetic opioids, notably fentanyl, and (4) an increase in deaths related to stimulant use (Ciccarone 2021, CDC 2022). Opioids are a class of drugs that reduce the intensity of pain signals and are commonly prescribed for pain management following an injury or for some health conditions such as cancer (CDC 2022). Well-known prescription opioids include oxycodone, hydrocodone, codeine, and morphine, which share the active ingredient with heroin. Misuse of prescription opioids can cause dependency on the pill or lead to secondary addictions to heroin or other synthetic opioids such as fentanyl. Medically known as opioid use disorder (OUD), this condition is defined as taking larger amounts of opioids than prescribed because of cravings, resulting in disruptions to work, school, or home obligations, and eventually leading to impairment and distress, according to the Diagnostic and Statistical Manual of Mental Disorders (American Psychiatric Association 2013). In addition, opioid misuse carries the potential for both fatal and nonfatal overdoses, posing significant risks to individuals’ lives and well-being.

Unprecedented promotion and pharmaceutical marketing campaigns since the 1990s significantly increased opioid prescriptions by healthcare providers (Van Zee 2009, Manchikanti et al. 2012, ASPA 2021), eventually contributing to widespread misuse and addiction in the population, because of both iatrogenic and pill diversion factors (Fan et al. 2019). OUD due to prescription misuse can be iatrogenic, leading to the patient’s addiction as a result of excessive prescription or misuse of medical care. OUD via pill diversion refers to the illicit use and misuse of opioid pills prescribed to someone else (Pitt et al. 2018). The Centers for Disease Control and Prevention (CDC) estimated the economic burden of prescription misuse, including the costs of healthcare, deaths, addiction treatment, and criminal justice involvement, to be $78.5 billion a year (Florence et al. 2016). Analyzing the Drug Enforcement Administration data, Joranson and Gilson (2005) show that nearly 6.5 million dosage units of opioid pills were diverted across 22 Eastern states. Moreover, overdose deaths remain a leading cause of injury-related death in the United States. Between 1999 and 2017, overdose deaths attributed to opioid prescriptions in the United States experienced a substantial rise from 3,442 to 17,029, making up nearly 40% of the total opioid-related overdose deaths (National Institute for Drug Abuse 2018). The significant individual and public health burden of this national crisis makes it a priority to carefully study and mitigate the devastating effects of OUD due to prescription misuse.

There have been many intervention efforts aimed at reducing the burden of the opioid pill diversion (HHS 2013). Although opioid prescription dispensing has significantly decreased with public health efforts, dispensing rates remain high in certain counties (CDC 2022). The increased emphasis on prescription drug monitoring programs (PDMPs) underscores the need to remain vigilant (Adalbert et al. 2023). PDMPs provide health professionals with an up-to-date history of controlled substance dispensing, to reduce doctor shopping, that is, visiting multiple physicians to obtain multiple prescriptions (Bao et al. 2016, Patrick et al. 2016). Prescription-limiting laws have been controversial as several medical societies have called on lawmakers to calibrate the trade-off between reducing misuse and still enabling proper access for patients with severe, chronic painful conditions (such as cancer) that require opioid pills for pain management (Pope 2018). Furthermore, Chen et al. (2019) argue that interventions such as PDMPs may ultimately have minimal effects on the number of opioid overdose deaths. Another intervention strategy is increasing the availability of naloxone (a medication designed to rapidly reverse opioid overdose) to reduce overdose mortality rates. However, there has been some controversy around the proper implementation of these types of interventions. For instance, Doleac and Mukherjee (2018) challenge the impact of promoting naloxone access and argue that the resulting net effect of such interventions is either zero or an increase in opioid misuse, mortality rates, and other unintended social harms.

Public health practitioners and policymakers are currently faced with the problem of selecting a combination of interventions to minimize the negative health outcomes of the opioid epidemic, under tight budget restrictions. Primary prevention efforts focus on intervening before OUD occurs, such as controlling the prescribing trends, whereas secondary or tertiary prevention (also referred to as mitigating interventions) focuses on early identification of OUD for supportive treatment, as well as harm reduction strategies such as naloxone distribution or similar overdose mortality-mitigating measures. Preventive and mitigating interventions are distinctive in terms of health benefits and effect sizes. Mitigating interventions usually require larger investments with immediate, significant, and shorter-term benefits (Coffin and Sullivan 2013). However, preventive interventions usually result in longer-term but less pronounced benefits with lower administration costs (Weiner et al. 2017). Nevertheless, both intervention types may or may not be perceived as cost-effective based on the prevailing level of opioid dependency in a population; see, for example, Perrone and Nelson (2012) and Mueller et al. (2015).

This study presents a decision-making framework built on the expected utility theory to support public health practitioners in effectively allocating their budget to the available intervention options, aiming to minimize adverse health effects in individuals with OUD. Considering the intricate dynamics of pill diversion and misuse, the proposed decision-making framework models the optimal selection of primary prevention and/or mitigating interventions to tackle the iatrogenic and diversion addiction pathways leading to OUD. To the best of our knowledge, this is the first mathematical formulation in the operations research and decision science literature that aims to inform the opioid crisis specifically focusing on pill diversion. The complexity of modeling OUD has been a significant factor contributing to the limited exploration in this area (Hyman 2005). Motivated by the seminal addiction literature by Freedman and Senay (1973), we develop a Markov decision process (MDP) model as a decision-making tool to identify the optimal selection of interventions to minimize the negative health outcomes of the epidemic at the current state given a restricted budget and in a certain population. The population dynamics is described via a susceptible-infected-recovered (SIR) compartmental model to represent the theoretical pathway from the state of opioid-naive, to prescription misuse or pill diversion, leading to OUD and eventually possible fatal overdose. The SIR compartmental model, inspired by the Vietnam War heroin addiction modeling, illustrates the spread of prescription opioid pills as a communicable disease (Green 2017), and it is utilized to dynamically assess OUD and overdose incidence, reflecting individuals’ transitions across compartments (Burke 2016, Pitt et al. 2018). We present an approximate solution algorithm to solve for large instances. We perform a comprehensive numerical analysis to illustrate the capabilities of our proposed model to allocate a budget to uncover underlying trade-offs between the effect sizes and the health costs associated with a set of intervention candidates for any stage of the epidemic. We also use simulation to evaluate the evolving consequences of implementing our proposed decision tool over time.

The remainder of this paper is organized as follows. In Section 2, we summarize the seminal literature on prescription pill diversion and misuse and, more broadly, the opioid epidemic. In Section 3, we describe the model of curbing opioid pill diversion, including the MDP model, the SIR compartmental model, and the approximation solution algorithm. We discuss data sources and parameter estimation in Section 4. In Section 5, we present our numerical analysis, including model validation, dynamic implementation of policies, and sensitivity analysis. We conclude and discuss the practical implications in Section 6.

2. Literature Review

Many efforts are being made to prevent and mitigate the negative consequences of the opioid epidemic. Preventive interventions aim to reduce the number of individuals who misuse opioid pills by (i) patient-level policies, for example, community-based prevention programs (Albert et al. 2011) or prescription drug take-back programs (Fass 2011, Stewart et al. 2015); (ii) prescriber-level policies, for example, prescription monitoring program (Baehren et al. 2010, Haffajee et al. 2015), primary/secondary/tertiary prevention policies (Kolodny et al. 2015), guidelines targeting prescriber behaviors (Chang et al. 2018), or peer-comparison opioid-prescribing feedback (Andereck et al. 2019); and (iii) regulator-level policies (Franklin et al. 2015, Rutkow et al. 2015). Mitigating interventions, on the other hand, aim to treat or alleviate the current cases of OUD by (i) expanding treatment options, for example, medication-assisted treatment (MAT) (Volkow et al. 2014, Connery 2015); (ii) directly reducing overdose death rates, for example, increasing the availability of naloxone (Maxwell et al. 2006, Straus et al. 2013); and (iii) reducing other substance-use-caused morbidity rates, for example, needle exchange programs to reduce infection among heroin users (Kaplan and O’Keefe 1993). More recently, Irvine et al. (2019) constructed a Bayesian hierarchical latent Markov process model to estimate monthly overdose and overdose-death risk along with the impact of interventions. They found the combined impact of three interventions (the take-home naloxone program, the introduction of overdose prevention service, and the uptake of opioid agonist therapy) is effective in averting overdose deaths during British Columbia’s opioid overdose crisis. All these studies focused on either primary preventive or mitigating interventions and synergies between interventions. However, the optimal allocation of a given budget to primary preventive or mitigating interventions has not been studied in this stream of literature.

In a review of recent harm-reduction solutions to the opioid epidemic, Hawk et al. (2015) recognize the important role of both prevention and treatment in reducing overdose fatalities. Reviewing studies that evaluate interventions targeting prevention, treatment, or management of opioid misuse and overdose, Barbosa et al. (2020) ask for future economic evaluations that consider interactions between interventions and examine combinations of interventions to inform the optimal policy response. In addition, Jalali et al. (2020) call for the need for more systems science approaches that can uncover the complexities of the opioid crisis, and help evaluate, analyze, and forecast the effectiveness of ongoing and new policy interventions.

There are a few papers in the literature that discuss the assessment of the cost-benefit of different interventions in the context of OUD, and more broadly substance use disorder (SUD). To evaluate the cost-benefit of alternative interventions, one needs to know how much is being spent on them and what benefits accrue from that spending. Rydell et al. (1996) present a model that estimates the cost-effectiveness of “supply control” and “demand control” programs to control the cocaine epidemic in the United States. They find the demand control programs (i.e., treating heavy users) to be the most cost-effective. Reviewing evidence-based interventions to prevent or reduce the damage that illicit drugs cause to public health, Strang et al. (2012) outline the potential effects of existing interventions if implemented thoroughly. They highlight the need for a guide to facilitate the selection of policies that maximize the public good. More recently, Homer and Wakeland (2021) develop a system dynamics simulation model to anticipate intended and unintended intervention effects. They find that no single intervention significantly reduces both the number of persons with OUD and overdose deaths, but a combination of interventions can do so. In a cost-effectiveness study, Fairley et al. (2021) show that MAT combined with contingency management, overdose education, and naloxone distribution to treat OUD is associated with significant health benefits and cost savings compared with no treatment, yielding savings of $25,000 to $105,000 in lifetime costs per person.

There are also a few studies that employ compartmental modeling of U.S. adults’ opioid use to estimate health outcomes of available interventions to mitigate the opioid epidemic. For example, using a compartmental model, Pitt et al. (2018) state that mitigating interventions focusing on individuals with OUD improve population health without harming any groups. They conclude that even though preventive interventions that reduce the prescription opioid supply could generate positive health benefits in the long term, they may increase heroin use and reduce quality of life in the short term. Extending and updating the compartmental model developed by Pitt et al. (2018), Rao et al. (2021) measure quality-adjusted life years (QALYs) and opioid-related deaths over 5 and 10 years to assess the effectiveness of interventions. Pitt et al. (2018) and Rao et al. (2021) call for a portfolio of interventions for eventual mitigation. Even though there are some similarities, our study is different from this stream of literature in several ways. First, whereas some seminal studies like Pitt et al. (2018) use quantitative approaches such as simulation models to evaluate the impact of one intervention on overall health outcomes, we focus on the stochastic dynamic optimization modeling to recommend the best course of action for a given stage of the epidemic. Our model takes healthcare costs and the impact of interventions quantified by public health practitioners as input parameters and recommends the optimal allocation of a budget to a set of interventions given the current state of the epidemic. The model can be updated to accommodate the changing nature of the epidemic to incorporate future waves of the opioid overdose crisis. Second, in our compartmental modeling, we narrow our focus to OUD and mortality that originated from prescriptions and that were caused by pill diversion.

MDPs, as a class of stochastic dynamic programs, have emerged as a powerful tool widely utilized in both operations research and decision sciences. They are particularly valuable for generating optimal policies in the realm of public health and medical interventions, where uncertainty and limited resources are prevalent. They are used to model illicit drug use in Australia (Caulkins et al. 2007), design and implement a dynamic health policy for the control of a novel strain of influenza virus (Yaesoubi and Cohen 2011a, b), address the decision of when to perform breast cancer biopsy (Chhatwal et al. 2010), reduce unnecessary mammography follow-ups (Alagoz et al. 2013), optimize the first response to sepsis (Rosenstrom et al. 2022), optimize resource allocation for intervening in human immunodeficiency virus (HIV) infection (Coşgun and Büyüktahtakın 2018), and model opioid prescribing trade-offs for pain management (Gökçınar et al. 2022), among many other healthcare and public health policy applications. In our paper, as the first milestone, we develop an MDP model of curbing prescription opioid pill diversion and misuse to specify the portfolio of interventions that minimizes the negative health outcomes of opioid pill diversion.

In sum, even though the existing literature has examined the impact of preventive or mitigating interventions on health outcomes, the optimal selection of a portfolio of such interventions is crucial to strategically address the ever-evolving opioid crisis and it is not carefully investigated. We model the opioid crisis caused by prescription pills and develop a decision-making framework that integrates compartmental modeling for the disease spread and the MDP to select the optimal portfolio of interventions to minimize the total health costs of the population of interest.

3. Model of Curbing Opioid Pill Diversion

We model the epidemic dynamics as a discrete-time Markovian process. Through a separate SIR model, we replicate access to prescription opioids and the ensuing pill diversion dynamic as a transmissible condition, where the population is classified into three compartments based on the disease stage: susceptible (S), infected (I), and removed (R). This approach allows us to account for the spread of opioids across individuals and populations, akin to the transmission of contagious diseases. By accounting for transmissibility in our study, we aim to gain valuable insights into the patterns and dynamics of opioid pill diversion, enabling us to develop targeted and tailored strategies to effectively curtail the propagation of the epidemic and address the broader public health crisis. Figure 1 demonstrates the high-level dynamics of our proposed framework and the interaction between the MDP and SIR models during each decision period.

Note that if the proposed decision framework is implemented dynamically, the optimal policy can inform the SIR parameters and alter the state of the population, indicating why there is an arrow from the optimal policies to field data, symbolizing the influence of policy decisions on population health outcomes. In Section 5.2, we further discuss the dynamic implementation of the optimal policies. Table 1 summarizes a complete list of notations used in this study.

Table 1. Summary of Notations

Within each decision period, for example, a one-year time frame, the public health practitioner is allocated a predetermined budget to manage the opioid crisis. This budget must be strategically allocated among a mix of primary preventive and mitigating interventions, taking into account the existing state of the epidemic. This budget assignment decision is the eventual output of the framework at each decision period through its MDP module. The required inputs to the MDP model include the fixed, predetermined budget for each period, the impact of each intervention per unit of budget invested, and the initial state of the epidemic.

Investing in the possible interventions and considering other external factors such as pill diversion and prescription trends are expected to change the state of the epidemic during the decision period and up to the next decision epoch, as replicated by the SIR model. Hence, the SIR model generates and updates the transition probabilities of the MDP model for the next decision period. This dynamic process enables the public health practitioner to incorporate the evolving nature of the epidemic and continued increases in overdose deaths into the decision-making framework.

The research question that this framework attempts to address is, given the dynamic state of the epidemic, and over a long planning horizon (e.g., 10 years), how to find the sequential optimal decisions of assigning the total fixed budget to the two main categories of interventions: (1) primary preventive interventions that would limit the prescription supply, and consequently reduce access to prescription opioids, and (2) mitigating interventions that would reduce the harm resulting from access to prescribed opioids by reducing the fatal overdose rate.

3.1. Key Underlying Assumptions

This paper exclusively addresses opioid pill diversion and misuse within the opioid epidemic, while excluding examination of illicit opioid access and its potential impact on the crisis. The foundation of our model is built upon the following key assumptions:

• In our MDP-based epidemic modeling, we assume complete observability of the epidemic state (e.g., susceptible and infected individuals) at discrete time intervals during the decision-making process. Moreover, we treat the variables defining the epidemic state as independent and mutually exclusive, leading to transition probabilities solely dependent on the immediate epidemic state and applied interventions. This assumption is well established in seminal literature (e.g., see Yaesoubi and Cohen 2011a, b). We posit that information about the number of new infections within a given period does not alter the probability distribution for the number of newly infected or removed individuals within the same period. This rationale is sound because the probability of opioid misuse is influenced by individual choices rather than solely relying on the overall prevalence of OUD within the population. Also, the independence assumption appropriately captures the dynamics of the opioid epidemic, setting it apart from other infectious diseases.

• The population for which the decision is being made is both fixed and known, similar to the seminal epidemic literature (Yaesoubi and Cohen 2011b, Tolles and Luong 2020). This assumption is justified by the fact that epidemics usually last over several years, whereas the rate of deaths from natural or other causes is generally small relative to the population size and assumed to equate to the birth rate. Also, the subpopulation in each compartment is considered homogeneous; that is, individuals in the same compartment are identical and mix uniformly among themselves (Brauer and Castillo-Chavez 2012).

• We consider a fixed budget to be assigned to two main categories of interventions introduced above. Each category of intervention may consist of several interventions as long as they impact the same transition rate. For example, primary preventive interventions could include prescription drug monitoring programs, drug disposal programs, and educational guidelines targeting prescriber behaviors because all alter the transition away from susceptibility.

• We assume the magnitude of the impact of interventions to be linear (see Section 3.2.3). That is, as the proportion of the budget assigned to an intervention category increases, the impact magnitude increases linearly. See case 2 in Section 5.3 for the discussion on nonlinear intervention impacts.

• We exclude a recovery pathway from the “infected” compartment to “susceptible,” as our susceptible compartment consists of opioid-naive individuals, who have not been regularly exposed to opioids and have no prior experience with these drugs (Stokowski 2010). The recovered individuals often exhibit different dynamics and incidence rates compared with opioid-naive individuals, potentially rendering them highly susceptible to opioid use once again. Furthermore, the recovery transition is infrequent, occurring in less than 10% of cases, and involves a prolonged process of mental and physical rehabilitation (Smyth et al. 2010). Because of our decision periods being relatively short compared with the time required for recovery, we assume that all individuals with OUD remain infected during each decision period. Note with this assumption that our modeling approach adopts a conservative stance, which is beneficial to public health practitioners seeking to implement impactful societal changes during a national emergency. However, this assumption does not disregard the possibility of self-recovery outside of the planning horizon; individuals who have fully recovered may not be counted as infected anymore in the next decision period. Focusing on the interventions designed to address the crisis, our model provides a cautious and pragmatic framework that prioritizes the effectiveness of external measures over relying on individual willpower. This approach ensures a comprehensive and proactive strategy to tackle the challenges posed by the national emergency and drive meaningful societal transformations.

3.2. MDP Model

The components of our MDP model are as follows:

State space. Let $X_S(t)$ denote the number of susceptible individuals, $X_I(t)$ denote the number of infected individuals, and $X_R(t)$ denote the number of removed individuals at time t. For a fixed population of size N, the state of the epidemic spread at any time t is ${\mathit{s}}_t=(X_S(t),X_I(t))$. For the convenience of notation, we drop the dependency of s on t. The state space is then $\mathcal{S}=\{(x_S,x_I)\in {\mathbb{N}}^2|x_S+x_I\le N\}$.

Decision epochs. The decisions are made at discrete time epochs $t=\{t_0,t_1,\dots ,T\}$, where T is the length of the decision horizon. Let $\mathrm{\Delta }t$ be the time interval between two consecutive decision epochs. To control the epidemic, a public health practitioner may continue to make decisions until the eventual containment. However, the stochastic nature of transitions makes it difficult to identify the time of this eventual ideal state. Hence, we consider an infinite decision horizon, that is, $t=\{t_0,t_1,\dots \}$.

Decision space. We consider two possible categories of interventions to control the epidemic: (1) primary preventive interventions and (2) mitigating interventions. We define $\mathcal{A}=\{({\text{a}}_P,{\text{a}}_M):{\text{a}}_P+{\text{a}}_M\le 1\}$ as the set of all possible combinations of budget proportions that can be allocated to primary preventive (${\text{a}}_P$) and mitigating (${\text{a}}_M$) interventions. The sum of these proportions cannot exceed one, indicating the fixed available budget for each decision period. At the end of this subsection, we provide a precise definition of the impact magnitude of interventions.

To formulate the optimality equation, we need to determine the transition probabilities within the MDP model. We employ SIR compartmental models, commonly used in epidemiology literature (Tolles and Luong 2020), to estimate the probabilities of transitioning between different states within the MDP. We discuss the proposed SIR model in Section 3.2.2.

3.2.1. Costs and Optimality Equation

Our objective is to determine the optimal intervention ratio to minimize negative health outcomes of the epidemic within the allocated budget. We assume that a fixed budget of B has been designated for allocation over a specified decision period of interest. There are two types of costs: (1) the cost associated with an infected individual with OUD (c_I) and (2) the cost of each fatal overdose (c_O). We represent the total cost when the epidemic is in state $\mathit{s}=(X_s(t),X_I(t))$ by $r(\mathit{s})$. Note the population of infected individuals $X_I(t)$ consists of individuals with access to opioid prescriptions, divided into two subcompartments: those with access to diverted opioids (p_d fraction of transition from susceptible to infected) and those with access to opioid prescriptions ($1-p_d$ fraction of transition from susceptible to infected). Within the subcompartment with access to diverted opioids, there is a possibility of misuse leading to OUD, denoted by the probability τ. Similarly, within the subcompartment with opioid prescriptions, there is a probability of misuse and subsequent development of OUD, denoted by p_m. Consequently, the cost c_I is attributed solely to $p_d\tau +(1-p_d)p_m$ portion of the infected individuals $X_I(t)$. Considering the removed individuals $X_R(t)$ as those in the total population N who are neither infected nor susceptible, the total cost $r(\mathit{s})$ is $c_I(p_d\tau +(1-p_d)p_m)X_I(t)+c_O(N-X_S(t)-X_I(t))$. See Section 3.2.2 for further elaboration on the dynamics described, which inform the necessary parameters for the MDP model. The optimality equation to minimize the expected total discounted cost when the epidemic is in state $\mathit{s}=(X_s(t),X_I(t))$ and decision $\mathit{a}=({\text{a}}_P,{\text{a}}_M)\in \mathcal{A}$ is made at time t is as follows:

where $\gamma \in [0,1]$ is the discount factor multiplied by the expected cost of the epidemic over the remaining time periods. The discount factor evaluates the importance of the accumulated future events in the present value. When the time span for the Markov chain is sufficiently large (like the epidemic studied in this paper), the time value of money should be taken into account in adding the costs of future periods to the costs of the current period.

3.2.2. SIR Compartmental Model

Our proposed SIR compartmental model for the opioid epidemic is composed of three umbrella compartments, susceptible, infected, and removed. Figure 2 depicts these compartments and their components, along with the transition pathways between them.

The “susceptible” compartment refers to the segment of the population that is opioid naive (box (a) in Figure 2). An individual in the susceptible compartment may transition to the “infected” compartment at transition rate η, representing the overall rate of obtaining prescription pills among opioid-naive individuals. This access to prescription pills is possible through either (1) the diversion pathway, with probability p_d (box (b) in Figure 2), or (2) the iatrogenic pathway, with probability $1-p_d$ (box (d) in Figure 2). Individuals with access to prescribed opioids can misuse them and develop OUD through two pathways. Firstly, people who obtain prescription pills via diversion might improperly use these pills despite not being intended recipients, with a probability of τ (box (c) in Figure 2). Secondly, those who are legitimately prescribed opioids for medical reasons might misuse their medication (Pitt et al. 2018), with a misuse probability of p_m, leading to iatrogenic addiction (box (e) in Figure 2). Therefore, the infected compartment comprises all individuals who have access to prescribed opioids, regardless of whether they misuse them.

Individuals with OUD may transition from the “infected” compartment to the “removed” compartment (i.e., fatal overdose) at transition rate μ (box (f) in Figure 2). Individuals who experience a nonfatal overdose stay in the infected compartment. Given we are modeling the development of access and pathways to OUD, individuals with prior opioid prescription access or those who have recovered from OUD are not identical to pain-free nonusers or opioid-naive persons in terms of relative susceptibility (Lemonick and Park 2007). In other words, the model focuses on the development of OUD from opioid prescription access, that is, incidence, and not on prevalence. Next, we introduce the driving events, discuss their probability distributions, and derive the transition probabilities between compartments.

3.2.2.1. Probability Distribution of Driving Events

There are two driving events in our SIR model during the period $[t,t+\mathrm{\Delta }t]$: (1) susceptible individuals obtain access to opioids and transition to the infected compartment, and (2) individuals in the infected compartment are removed because of fatal overdoses. Let S(t) be the number of susceptible individuals during the period $[t,t+\mathrm{\Delta }t]$, I(t) be the number of infection incidents during the period $[t,t+\mathrm{\Delta }t]$, and R(t) be the number of removed individuals from the infected compartment during the period $[t,t+\mathrm{\Delta }t]$. The random variables S(t), I(t), and R(t) are independently distributed and are determined only by the state of the epidemic at time t; see Section 3.1 for more details on our assumptions. Next, we discuss the probability distribution of the driving events I(t) and R(t), given the current state of the epidemic $(X_S(t),X_I(t))$.

Susceptible to infected. Two possible transitions are described below:

• (1) Pill diversion. The incidence of pill diversion depends on (i) the contact rate between susceptible and infected individuals and (ii) the size of the infected population. Similar to other communicable epidemics (Tolles and Luong 2020), the contact rate represents the average number of susceptible individuals with whom a person with access to the prescription makes sufficient contact to pass the infection (in our case, access to opioids) per unit time (Manfredi and D’Onofrio 2013). That is, the contact rate is the product of (1) the number of interactions between individuals in susceptible and infected compartments per person per unit time and (2) the probability of pill diversion once interacting. The product of (1) and (2) is $p_d\eta \mathrm{\Delta }t$. In addition to the contact rate, the incidence of opioid pill diversion depends on the size of the compartment with access to prescription, that is, $\beta (t)=X_I(t)/N$, which is the probability that the next interaction of a random susceptible is with an infected individual.

  We posit that during the interval $[t,t+\mathrm{\Delta }t]$, a randomly chosen susceptible individual will interact with n people and obtain diverted opioid pills, where the number n follows a Poisson distribution characterized by the rate $p_d\eta \beta (t)\mathrm{\Delta }t$. Let ${\rho }_d(t)$ denote the overall probability that a susceptible person becomes infected through pill diversion. The probability ${\rho }_d(t)$ is then calculated as follows:

  $${\rho }_d(t)=1-e^{-p_d\eta \beta (t)\mathrm{\Delta }t}.$$(2)

• (2) Medical cause for prescriptions. We consider the encounter rate between susceptible individuals and the healthcare system, analogous to contact rates in pill diversion, but focused on healthcare interactions. This encounter rate is quantified by $(1-p_d)\eta$, representing the average interactions per unit time a susceptible person has with healthcare providers, where each interaction has a high likelihood (effectively one) of resulting in an opioid prescription for medical reasons. Within a specific interval $[t,t+\mathrm{\Delta }t]$, a susceptible person is expected to engage in m such healthcare interactions, which could lead to being prescribed opioids medically. The quantity m follows a Poisson distribution with a rate of $(1-p_d)\eta \mathrm{\Delta }t$, capturing the dynamics of medical access to opioids and the prescribing habits within the healthcare system, offering a structured yet distinct pathway compared with the interpersonal transmission seen in pill diversion. Thus, the probability that a susceptible individual transitions to the infected compartment because of medical cause for prescriptions during the interval $\mathrm{\Delta }t$, denoted by ${\rho }_m(t)$, is given as follows:

  $${\rho }_m(t)=1-e^{-(1-p_d)\eta \mathrm{\Delta }t}.$$(3)

Hence, given state $(X_S(t),X_I(t))$ at time t, I(t) has a binomial distribution with $X_S(t)$ number of trials and probability of success $q_I(t)$. Then, assuming pill diversion and medical cause for prescriptions are independent and using the total law of probability, $q_I(t)=p_d{\rho }_d(t)+(1-p_d){\rho }_m(t)$. Thus, the probability that i individuals leave the susceptible compartment to the infected compartment, given the current state of the epidemic $(X_S(t),X_I(t))$, is

Infected to removed: We assume infected individuals are removed from the population because of fatal overdoses. We consider individuals with a nonfatal overdose to stay in the infected compartment unless they meet the criteria for “removal.” We assume the communicability period before fatal overdose is exponentially distributed with rate μ. Thus, the probability that an infected individual has a fatal overdose during the interval $\mathrm{\Delta }t$ can be calculated as ${\rho }_f(t)=1-e^{-\mu \mathrm{\Delta }t}$.

Remember τ refers to the probability that pill diversion results in addiction. Therefore, $p_d\tau$ represents the proportion of OUD after pill diversion. In addition to communicable OUD due to the pill diversion, we assume the non-OUD population with access to prescriptions ($1-p_d$ proportion of the infected or those with access to opioids) may misuse the pills and become addicted (i.e., iatrogenic addiction) with probability p_m. Therefore, the number of individuals with OUD at time t who transition to the removed compartment during the period $[t,t+\mathrm{\Delta }t]$ has a binomial distribution with $(p_d\tau +(1-p_d)p_m)X_I(t)$ number of trials and the success probability ${\rho }_f(t)$. Thus, the probability that r infected individuals transition to the R compartment during the interval $\mathrm{\Delta }t$, given the current state of the epidemic $(X_S(t),X_I(t))$, is

where $\lfloor x\rfloor$ means the greatest integer lesser than or equal to x.

Similar to the MDP literature for the epidemic modeling (see, for example, Yaesoubi and Cohen 2011a), we assume that a susceptible individual who is infected during $[t-\mathrm{\Delta }t,t]$ becomes infected at time t and interacts with the rest of the population during period $[t,t+\mathrm{\Delta }t]$. The individual is then removed from the population at the end of the period at time $t+\mathrm{\Delta }t$. This simplifying assumption allows us to apply the MDP to identify the optimal decision at each decision epoch, taking into account the worst possible outcome in a conservative manner.

3.2.2.2. SIR Dynamics and Transition Probabilities

Using the probability distribution of the driving events as discussed above, we can write the dynamics of the driving equations of the SIR model as follows:

We assume that I(t) and R(t) are independent. The dynamics feasibility constraints for the SIR model are then given as follows:

Note that only $(p_d\tau +(1-p_d)p_m)$ proportion of the infected individuals can transition to the removed compartment. Therefore, using inequalities in Equation (7), the set of possible transitions, that is, the probability support of the random variables $X_S(t+\mathrm{\Delta }t)$ and $X_I(t+\mathrm{\Delta }t)$ given the current state variables $X_S(t)$ and $X_I(t)$, is

Considering Equation (6) and the state space of transitions ${\mathcal{S}}_{(X_S(t),X_I(t))}$, we calculate the transition probabilities of the Markov chain $\{(X_S(t),X_I(t)):t=0,1,\dots \}$ as follows:

3.2.3. Impact Magnitude of Interventions

Implementing primary preventive interventions at time t reduces ${\rho }_d(t)$ in Equation (2) to ${\rho }_d(t)=1-e^{-(1-{\alpha }_P)p_d\eta \beta (t)\mathrm{\Delta }t}$ and ${\rho }_m(t)$ in Equation (3) to ${\rho }_m(t)=1-e^{-(1-{\alpha }_P)(1-p_d)\eta \mathrm{\Delta }t}$, where ${\alpha }_P\in [0,1]$ is the fractional reduction in the infection transmission rate and assumed to be a linear function of ${\text{a}}_P$. To ensure the robustness of our results, we examine the use of a nonlinear function, which will be discussed as case 2 in Section 5.3. We let ${\alpha }_M\in [0,1]$ denote the fractional reduction in the fatal overdose transmission rate and assume it to be a linear function of ${\text{a}}_M$. Similarly, implementing mitigating interventions at time t reduces probability ${\rho }_f(t)$ to ${\rho }_f(t)=1-e^{-(1-{\alpha }_M)\mu \mathrm{\Delta }t}$.

3.3. Solving the MDP

For a population of size N, the transition probability matrix of the Markov chain $\{(X_S(t),X_I(t)),t=1,2,\dots \}$ is of size $(N+1)(N+1)$, which quickly becomes computationally intractable for large populations. On the other hand, analyzing a Markov chain with absorbing states is complex (as the chain is not irreducible), and thus the characterization of the optimal policy is difficult. Therefore, we propose an approach, adopted from Yaesoubi and Cohen (2011b), for approximating the Markov chain $\{(X_S(t),X_I(t)),t\ge 0\}$ with the Markov chain $\{({\mathrm{\Theta }}_S(t),{\mathrm{\Theta }}_I(t)),t\ge 0\}$, where ${\mathrm{\Theta }}_S(t)$ is the proportion of susceptible population at time t and ${\mathrm{\Theta }}_I(t)$ is the proportion of infected population at time t. The proportions ${\mathrm{\Theta }}_S(t)$ and ${\mathrm{\Theta }}_I(t)$ can only take a limited number of values from the sets $\{{\theta }_1^S,{\theta }_2^S,\dots ,{\theta }_{d_S}^S\}$ and $\{{\theta }_1^I,{\theta }_2^I,\dots ,{\theta }_{d_I}^I\}$, respectively. We use a policy development algorithm, explained in Appendix B, to solve the MDP with the approximate state space.

To determine the set $\{{\theta }_1^S,{\theta }_2^S,\dots ,{\theta }_{d_S}^S\}$, we divide the interval $[0,1]$ into d_S regions with equal sizes of $1/d_S$. Let the boundary points be $\{b_0^S,b_1^S,\dots ,b_{d_S}^S\}$, where $0=b_0^S<b_1^S<\cdots <b_{d_S}^S=1$. Thus, the possible values for ${\mathrm{\Theta }}_S(t)$ are ${\theta }_i^S=(b_{i-1}^S+b_i^S)/2=(i-0.5)/d_S$, for $i=1,2,\dots ,d_S$. Similarly, dividing the interval $[0,1]$ into d_I equal regions and with boundary points $0=b_0^I<b_1^I<\cdots <b_{d_I}^I=1$, we can determine the set $\{{\theta }_1^I,{\theta }_2^I,\dots ,{\theta }_{d_I}^I\}$ by ${\theta }_i^I=(b_{i-1}^I+b_i^I)/2=(i-0.5)/d_I$, for $i=1,2,\dots ,d_I$.

We can then calculate the transition probabilities for the approximate Markov chain $\{({\mathrm{\Theta }}_S(t),{\mathrm{\Theta }}_I(t)),t\ge 0\}$ as follows:

Note that when a transition from $({\mathrm{\Theta }}_S(t),{\mathrm{\Theta }}_I(t))$ to $({\mathrm{\Theta }}_S(t+\mathrm{\Delta }t),{\mathrm{\Theta }}_I(t+\mathrm{\Delta }t))$ is not possible, that is, when $0\le {\mathrm{\Theta }}_S(t+\mathrm{\Delta }t)\le {\mathrm{\Theta }}_S(t),{\mathrm{\Theta }}_S(t)+(1-p_d\tau -(1-p_d)p_m){\mathrm{\Theta }}_I(t)\le {\mathrm{\Theta }}_S(t+\mathrm{\Delta }t)+{\mathrm{\Theta }}_I(t+\mathrm{\Delta }t)\le {\mathrm{\Theta }}_S(t)+{\mathrm{\Theta }}_I(t)$ does not hold, we should set the transition probability to zero. If the transition probability becomes greater than zero for an arbitrary state $(i,j),i<j$, an intuitive approximation is to add this probability to the corresponding diagonal state (i, i). The reader may refer to Appendix A for the details of the algorithm developed to approximate the transition probabilities for Markov chain $\{({\mathrm{\Theta }}_S(t),{\mathrm{\Theta }}_I(t)),t\ge 0\}$.

4. Data Sources and Parameter Estimation

In this section, using the publicly available opioid epidemic data and the existing literature, we estimate the parameters of the model.

4.1. Initial Compartment Sizes

Estimating the initial size of all compartments is necessary before calculating probabilities and transition rates. Our estimates rely on data from the 2020 U.S. population census and opioid-related statistics. According to the National Survey on Drug Use and Health (NSDUH), the total U.S. population aged 12 and older was 276.9 million in 2020 (Substance Abuse and Mental Health Services Administration 2021). To determine the initial size of each compartment at the beginning of the time period, we first calculate the infected and removed populations, and then approximate the susceptible population accordingly.

To gauge the initial size of the Infected compartment, we first estimate the number of individuals with iatrogenic addiction, representing OUD stemming from misusing prescribed opioids despite harm (box (e) in Figure 2). This misuse encompasses opioid usage contrary to directed or prescribed patterns, regardless of harm or adverse effects (Vowles et al. 2015). In 2020, approximately 3.4% (or 9.4 million) of individuals aged 12 or older in the United States misused opioids (excluding heroin) (Substance Abuse and Mental Health Services Administration 2021, table A.12B). Consistent with studies conducted in the United States and reported in Vowles et al. (2015), we assume the rate of addiction to be 11.5% among individuals who misuse opioid prescriptions, and estimate the population with iatrogenic addiction at 1.1 million (9.4 million people who misused × 11.5% addiction rate). Furthermore, with approximately 27% of those with prescriptions reported to misuse opioids (Naliboff et al. 2011), we estimate the population with access to opioid prescriptions (box (d) in Figure 2) at 34.8 million (9.4 million/0.27). For contagious addiction (box (c) in Figure 2), the NSDUH report (see table A.26B therein) indicates that 1.0% (or 2.8 million) of individuals aged 12 or older in the United States had OUD in 2020. Assuming OUD can arise from either contagious or iatrogenic addiction, we deduce the population with contagious addiction to be 1.7 million (2.8 million total OUD population − 1.1 million with iatrogenic addiction). Finally, by adopting a 6.8% probability of addiction for opioid-naive patients receiving diverted opioids (Pitt et al. 2018), we estimate the population with access to diverted opioids (box (b) in Figure 2) to be 25 million. In summary, we estimate the total size of the infected population at approximately 59.8 million individuals, comprising 34.8 million with access to opioid prescriptions (including iatrogenic addiction) and 25 million with access to diverted opioids (including contagious addiction).

Removed population includes individuals with OUD who have a fatal overdose (box (f) in Figure 2). According to the CDC (2022), prescribed opioids were involved in 49,860 overdose deaths in 2019 (70.6% of all drug overdose deaths). We, therefore, estimate the overdose deaths among people aged 12 or older in 2020 in the United States to be approximately 0.05 million.

Hence, the Susceptible population (including opioid-naive individuals including those who are pain free, those with acute pain, or those with chronic pain, represented in box (a) in Figure 2), is approximately 217.15 million individuals in the United States aged 12 and older, calculated as 276.9 million total model population minus 59.8 million Infected and 0.05 million Removed.

4.2. Transition Rates from S to I

We estimate the rate of infection, indicating access to opioids for opioid-naive patients, through two pathways: opioid prescriptions and diversion. To determine the rate of opioid prescriptions, we consider the proportion of people aged 12 or older in 2020 who are prescribed opioids as $34.8/276.9=0.126$. Therefore, we estimate the access rate to opioids through prescription as $(1-p_d)\eta =0.126$. For the diversion rate, we estimate the proportion of people aged 12 or older in 2020 with access to diverted prescribed opioids as $25/276.9=0.1$, and therefore we have $p_d\eta =0.1$. This aligns with the estimate provided by Pitt et al. (2018), which is derived from the product of “the number of interactions per person per month” and “the probability of interactions between the susceptible and infected population resulting in the diversion of opioids.” Their estimate indicates approximately 0.01 interactions per person per month resulting in the diversion of opioids to opioid-naive patients. With our study’s planning horizon set at 12 months (1 year), the access rate to opioids through diversion is estimated to be approximately 0.1 per year. Hence, we can now estimate the overall access rate to opioids (through prescription or diversion) as $\eta =p_d\eta +(1-p_d)\eta =0.1+0.126=0.226$ per year. Furthermore, we approximate the probability of misusing the prescription that leads to OUD for prescription-holding individuals to be p_m = 27% (misuse rate) × 11.5% (addiction rate) = 3% (Naliboff et al. 2011, Volkow et al. 2014). Finally, we estimate the probability of an individual transitioning from S to I through pill diversion as $p_d=0.1/0.226=44\%$.

4.3. Transition Rates from I to R

Infected individuals are removed because of fatal overdose. As discussed above, it is estimated that approximately 0.05 million individuals out of 2.8 million people with OUD suffer fatal overdoses annually in the United States, representing 2% of the OUD population (0.05 million/2.8). Thus, we estimate the fatal overdose rate as $\mu =0.02$ per year. The main input parameter values for the SIR compartmental model, as described in this section, are summarized in Figure 3.

5. Numerical Experiments

For a given state in the opioid epidemic, determining the most effective intervention portfolio relies not only on the impact of interventions but also on the public health practitioner’s perception of costs, specifically the marginal costs associated with individuals experiencing OUD and fatal overdoses in the population. To examine the performance of our proposed framework through numerical experimentation, we calibrate our MDP model for a population size normalized to 1,000 people (i.e., n = 1,000) for an annual decision epoch (i.e., $\mathrm{\Delta }t=1$). We choose $d_S=d_I=10$ in our MDP approximation procedure and assume a discount factor of $\gamma =0.9$.

This section begins with presenting a base case scenario, utilizing assumed estimates for addiction and fatal overdose marginal costs to validate the performance of our model. We then present the dynamic implementation of this decision framework and its impact on cost reduction through simulating the evolution of the epidemic under three scenarios: applying no intervention, applying a fixed policy, and applying the optimal intervention policy suggested by our framework. To conclude this section, we conduct a thorough sensitivity analysis on key model parameters, including the ratio of overdose and infection costs (i.e., $c_o/c_I$), transition rates, and probabilities, to assess the robustness of our model’s inputs and outputs. We illustrate the application of our framework under various scenarios to showcase the flexibility of the outputs and the optimal policies generated.

Before delving into the presentation of all the results, it is noteworthy to mention that we have incorporated a “smoothing strategy” across all numerical experiments in this section. The smoothing strategy is designed to improve the stability and consistency of our results by rectifying situations where the optimal policy for a particular epidemic state deviates from the optimal policies of all neighboring states, while all neighboring states share a uniform optimal policy. Recall the optimal policy for a given epidemic state is determined to minimize the expected total discounted cost of the epidemic. However, we observed cases in which the cost associated with the optimal policy of state (i, j) differs from (but is very close to) the cost associated with the optimal policy for all states around it, that is, $\{(i\prime ,j\prime )|i\prime \in \{i+\delta ,i-\delta \}\& j\prime \in \{j+\delta ,j-\delta \}\}$, where δ refers to the discretization level of states in the approximate Markov chain model. Therefore, when generating the final optimal policy plots in the subsequent sections, we implement the smoothing strategy to substitute the optimal policy with the alternative one, if it exists. This substitution assumes a predetermined threshold ϵ for the percentage change in total cost from optimality. The decision maker is open to alternative policies with slightly, albeit practically insignificantly, higher costs for a smoother intervention implementation, avoiding drastic changes in budget allocations as the epidemic evolves gradually. We consider $ϵ=1\text{E}-\text{}6$ in our numerical experiments.

5.1. Model Validation (Base Case)

We assume that public health practitioners always perceive the cost of a fatal overdose to be higher than that of an individual becoming addicted (i.e., c_O > c_I). Florence et al. (2021) estimated the annual value of reduced quality of life from OUD (i.e., c_I) to be around $183,185 and the annual value of life lost to opioid overdose (i.e., c_O) to be around $10.1 million, based on the data of 2017. Therefore, for the base case, we estimate that the public health practitioner perceives the negative consequences of one fatal overdose to be 55 times more taxing than having an addicted individual (i.e., $c_O/c_I=55$). Figure 4 shows the optimal policy for given epidemic states in the base case scenario. The dots on the graph represent the proportion of the budget allocated to the preventive intervention, selected from the set of 0.1, 0.5, and 0.9. Assuming that the entire budget must be spent on a combination of the two intervention types, the unallocated portion from the designated amount for preventive intervention signifies the budget allocated to mitigating interventions.

In Figure 4, as long as the susceptibility ratio is below 25%, it is optimal to allocate most of the budget to mitigation for minimizing overall costs by reducing overdose deaths. However, as the susceptible population increases, more of the budget should shift to prevention, and managing newly infected individuals. Optimal budget allocation involves a 50-50 split between prevention and mitigation when the susceptible population ranges from 25% to 35%, regardless of the infected population proportion. Exceptions arise for scenarios with 35% susceptibility and a 55% or higher infected rate, where it is optimal to prioritize preventive measures in budget allocation. A detailed cost analysis reveals that, in these instances, the cost difference between policy ${\text{a}}_p=0.9$ and policy ${\text{a}}_p=0.5$ is in fact less than $ϵ=1\text{E}-\text{}6$, a practically insignificant difference that could be smoothed out. Beyond 45% susceptibility, irrespective of the infected population proportion, it is optimal to allocate almost the entire budget to prevention.

Figure 5 illustrates the changing spectrum of optimal costs corresponding to the recommended optimal policy for the base case across various initial states. Figure 5 shows that initiating interventions only when a smaller proportion of the population is susceptible and infected would lead to higher costs attributed to increased opioid overdose deaths. We also observe that a 10% decrease in the initial proportion of susceptible individuals corresponds to approximately a $10 billion increase in total costs over the planning horizon. Next, we discuss the dynamic implementation of our proposed decision support tool spanning 10 consecutive years.

5.2. Dynamic Implementation of Policies

The findings presented in the previous section have depicted the optimal strategies for a planning horizon, specifically one year. These optimal policies are established for each given state of the epidemic, assuming a consistent pattern in epidemic dynamics where the transition rates from one compartment to the next remain unchanged. However, between consecutive planning horizons, the evolution of epidemic dynamics is contingent on various factors, including the interventions implemented during each specific planning period. In this section, we employ the modsim package in Python (Downey 2023) to construct a simulation model of the opioid epidemic. This model aims to investigate the costs associated with the unfolding epidemic, computed monthly, over 10 years, both with and without interventions. Utilizing the base case’s estimated model parameters (discussed in Section 4) with initial population fractions $({\mathrm{\Theta }}_S(0)=0.784,{\mathrm{\Theta }}_I(0)=0.215,{\mathrm{\Theta }}_R(0)=0.001)$, we recalculate monthly transition rates (${\eta }_i={\theta }_I^i/12$ and ${\mu }_i=X_R^i/12$) at the end of each year (for $i=1\cdots 10$). These updated rates are then applied to the subsequent year. The following differential equations, derived to capture the dynamics of the SIR model outlined in Section 3.2.2, are formulated for simulating the epidemic evolution over time:

where $\frac{d{\mathrm{\Theta }}_S(t)}{dt},\hspace{0.17em}\frac{d{\mathrm{\Theta }}_I(t)}{dt}$, and $\frac{d{\mathrm{\Theta }}_R(t)}{dt}$ refer to the rate of change in the size of the corresponding compartments per unit time. Consistent with our earlier assumption, no individuals are introduced to the susceptible group, as we disregard births and maintain a constant population throughout the simulation period.

First, we examine the changes in total cost, infection cost, and overdose cost under no intervention and compare them with the case where the optimal policy recommended by our model is applied. Figure 6 displays the changes in cost values, and Figure 7 depicts changes in the marginal costs (the difference in costs between the current and previous years).

According to Figure 6(a), a noticeable upward trend in total costs, including overdose and infection costs, emerges when allowing the epidemic to progress naturally over 10 years without interventions. Figure 7(a) shows that the marginal cost increases annually, reaching its peak quickly around year 3, with a continuous rise in overdose costs until the end of year 10. Although the marginal change in infection cost and total cost appears to decrease gradually, this does not indicate the natural containment of the epidemic. In fact, because of a substantial portion of the population being infected, the marginal changes are small, yet the total costs (comprising both overdose and infection costs) steadily increase. In contrast, looking at Figures 6(b) and 7(b) indicates that implementing optimal interventions each year over the 10 years results in a steady minimal overdose cost, with infection cost and total cost stabilizing shortly after the midpoint of the simulated time, around year 7. Notably, overdose cost constitutes the majority of the total cost in this scenario, whereas infection is effectively contained from the beginning.

Next, we assess how our optimal policy influences costs in contrast to a fixed policy, where the budget is consistently split equally between the two interventions, irrespective of the current state of the epidemic. For this comparison case, we define cost-saving measure $\text{CS}$ for policy π as ${\text{CS}}_{\pi }=(C_{\pi }-C^{*})/C_{\pi }$, where $C^{*}$ is the cost when the optimal budget allocation is used and $C_{\pi }$ is the cost when allocation π is used. Table 2 presents a summary of the comparison of cost savings between the optimal policy and a no-allocation policy, as well as the optimal policy versus a fixed 50-50 policy across a simulated epidemic spanning 10 years. As evident from the initial three columns of the table, the optimal policy consistently outperforms the no-allocation scenario, showcasing a progressively enhanced rate of overall cost reduction. Compared with the no-allocation case, implementing our proposed optimal policy yields an average improvement of 62% in infection costs, 13% in overdose costs, and 29% in total costs over 10 years.

Table 2. Cost Savings Comparison: Optimal vs. No Allocation and Optimal vs. Fixed Allocation

While enhancements are anticipated when allocating some budget toward mitigating the epidemic compared with a scenario with no allocation, it is intriguing to contrast the optimal policy, which adapts intelligently to the epidemic’s state at each decision point, with a fixed approach that remains oblivious to the epidemic’s evolution over time. These simplistic decisions are practical and straightforward, requiring no decision support tool or mathematical computations. However, as we observe in the second half of Table 2, in contrast to our proposed optimal policy, the fixed 50-50 policy leads to an average 49% rise in infection costs and a 12% increase in total costs, on average. The only case where the fixed policy presents some benefits is in the overdose cost with a 6% improvement. This is because the fixed policy guarantees a 50% allocation to mitigating measures every year and hence dramatically decreases the chance of overdose, even if that much budget allocated to mitigation is unnecessary for the current state of the epidemic, whereas the optimal policy adjusts to the state of the epidemic and favors preventive measures in many cases, with the objective of overall cost reduction over the years. As we can see in this experimental case, the optimal policy improves the total costs through increased investments in prevention. Figure C.1 in Appendix C shows the detailed cost trends over a 10-year simulation time for the fixed 50-50 allocation policy.

5.2.1. Impact of the Initial State of the Epidemic on Its Spread

We also investigate how the initial population proportions affect the epidemic’s trajectory to assess the effectiveness of prompt interventions over a 10-year simulation. To explore this, we conduct the simulation with two extreme sets for initial population proportions: (i) a high-susceptible, low-infected scenario (i.e., ${\mathrm{\Theta }}_S(0)=0.9,{\mathrm{\Theta }}_I(0)=0.099,{\mathrm{\Theta }}_R(0)=0.001$), and (ii) a low-susceptible, high-infected scenario (i.e., ${\mathrm{\Theta }}_S(0)=0.099,{\mathrm{\Theta }}_I(0)=0.9,{\mathrm{\Theta }}_R(0)=0.001$). In both cases, the proportion of removed individuals is assumed to remain constant at 0.1%. Figures C.2 and C.3 in Appendix C illustrate the total and marginal cost variations over 10 years for the initial proportion set (i). Similarly, Figures C.4 and C.5 in Appendix C depict the total and marginal cost changes over 10 years for the initial proportion set (ii). Contrasting situations with and without intervention highlights consistent patterns that underscore the importance of implementing optimal interventions. It demonstrates how the costs are controlled and eventually mitigated compared with allowing the epidemic to unfold naturally with no interventions, regardless of the initial state of the population. In particular, the patterns observed for the high-susceptible/low-infected set (i) in Figures C.2 and C.3 closely resemble the base case. However, in comparison with the base case, the low-susceptible/high-infected set (ii) in Figures C.4 and C.5 exhibits a drastic increasing trend in total costs, primarily driven by overdose costs in the absence of intervention. Furthermore, Figures C.4(b) and C.5(b) demonstrate that optimal interventions, even if introduced relatively late in the epidemic (given the initial population state in set (ii)), can effectively mitigate the consequences. It is noteworthy to see that in the initial population with high susceptibility and low infection, costs stabilize around the same time (after the midpoint of the simulated time) as the base case. However, in the initial population with low susceptibility and high infection, costs stabilize much earlier (within one year) than the base case, owing to the significant mitigating actions taken to address the consequences of a highly infected population. These observations indicate that the long-run impacts of implementing optimal interventions are not highly sensitive to the initial state of the population, provided they are not introduced excessively late.

5.3. Sensitivity Analysis and Robustness Checks

We consider the following cases to study the impact of various parameters and their magnitudes on the shape of the optimal policy:

• Case 1. Ratio of fatal overdose to addiction marginal costs

• Case 2. Impact magnitudes of interventions (nonlinear function)

• Case 3. Impact of parameters

Each of these cases is contrasted with the base case. To facilitate visual comparison, we overlay the two shaded areas initially depicted in Figure 4 to juxtapose each case with the optimal interventions suggested for the base case.

5.3.1. Case 1—Ratio of Fatal Overdose to Addiction Marginal Costs

In this study, we utilize the values provided by Florence et al. (2021) for the number of individuals affected by OUD and estimate a range for the ratio of fatal overdose to addiction marginal costs (i.e., $c_O/c_I$) to be between 48 and 62. Our findings reveal that the optimal interventions suggested by our model remain largely consistent with the base case across varying ratios of $c_O/c_I$ within and even beyond this range. However, considering the increasing trend in these marginal costs (CDC 2022), we also explore the optimal interventions recommended by our model for higher ratios of $c_O/c_I$, as depicted in Figure 8.

We observe in Figure 8 that as the $c_O/c_I$ ratio becomes much higher than the base case (i.e., 10⁴ and 10⁵), it is optimal to allocate a larger portion of our budget to mitigating interventions to control the epidemic by reducing the number of fatal overdoses as the cost of fatal overdose becomes proportionally higher.

5.3.2. Case 2—Impact Magnitudes of Interventions (Nonlinear Function)

In Section 3, we used the linear function to model the impact magnitude of implementing interventions. To examine the robustness of our results, in this case, we applied a square root function. The square root function, commonly used in utility theory and economic decision making (Samuelson 1983), exhibits diminishing marginal returns by gradually flattening the slope as the input increases. In our context, this represents a situation where the effectiveness of interventions is expected to level off after a certain threshold. We observed that the recommended optimal interventions remain almost the same as the base case. This implies that if we assume diminishing marginal returns on investments, the model’s outcomes appear unaffected by the particular mathematical function representing the decreasing trend.

5.3.3. Case 3—Impact of Parameters

We also perform sensitivity analysis on the main input variables: (1) transition rates from S to I (η) and from I to R (μ), and (2) probabilities τ, p_d, and p_m. Figures 9 and 10 demonstrate the impact of changing transition rates. We examine the impact of a 50% increase and decrease in the parameter values on the optimal policy. Figure 9 shows that as the transition rate from S to I decreases (fewer individuals become addicted), it is optimal to allocate a larger proportion of the budget to mitigating interventions to reduce the number of fatal overdoses. When the transition rate from S to I increases, it is optimal to allocate a larger portion of the budget to preventive interventions for most states, including those where the susceptible portion in the population is as low as 25%. Interestingly, Figure 10 shows that as the transition rate from I to R changes (the number of fatal overdoses changes), it is optimal to allocate a large proportion of the budget to preventive interventions. This recommendation may contradict intuition that suggests allocating a larger proportion of the budget to mitigating interventions to directly reduce the number of fatal overdoses. This is an interesting recommendation by the model to verify the saying that “an ounce of prevention is worth a pound of cure” as preventive interventions effectively reduce the number of infected individuals, which in turn decreases the fatal overdose cases in the long run. Note, based on supplementary experiments (omitted here for brevity), in situations where fatal overdoses are hypothetically exceptionally prevalent in the population, indicated by an overly high rate of transition from I to R (i.e., μ), it is beneficial to allocate a larger proportion of the budget to mitigating interventions even when the proportion of susceptible individuals exceeds 45% of the total population.

Figures C.6–C.8, in Appendix C, demonstrate the impact of changing probabilities. As shown in Figure C.6, as the probability of contagious addiction for an opioid-naive patient who receives diverted opioids (i.e., τ) changes, it is optimal to allocate a larger proportion of the budget to preventive interventions in states with lower proportions of susceptible individuals in the population. By implementing preventive interventions, we can reduce the number of infected individuals (consequently lowering the chance of pill diversion) as well as the number of fatal overdoses. Nonetheless, we do not observe substantial alterations to the optimal policy, when τ changes. As shown in Figure C.7, as the probability of an individual transitioning from S to I through pill diversion (i.e., p_d) decreases, it is optimal to allocate a larger proportion of the budget to preventive interventions to better manage the diversion of pills in the population. Similarly, this is because preventive interventions not only reduce addiction costs by curbing the number of newly infected individuals but also indirectly lessen the cost of fatal overdose. As shown in Figure C.8, as the probability of iatrogenic addiction for prescription-holding individuals (i.e., p_m) changes, it is optimal to allocate a larger proportion of the budget to preventive interventions to better manage the prescription of opioids to the population. Nevertheless, we do not observe significant changes to the optimal policy when p_m changes.

6. Discussion and Conclusion

Easy access to opioid prescriptions could potentially lead individuals to misuse the substance, develop OUD, and have an increased risk of fatal overdoses. It can also cause pill diversion to nonusers, resulting in more cases of OUD and a rise in fatal overdoses. Excessive opioid prescription and subsequent pill diversion eventually impose significant health costs on society. To minimize this taxing risk and enhance health outcomes, a handful of interventions (preventive or mitigating) are proposed (and implemented) in practice. However, the optimal allocation of a limited healthcare budget among mitigating and preventive interventions to minimize the opioid crisis health risks is not thoroughly investigated in the literature, to the best of our knowledge.

In response to the pressing opioid crisis, we introduced a robust decision framework that recommends optimal budget allocation strategies based on the current state of the epidemic and the predicted rates of addiction and fatal overdose that minimize its negative health impacts. By optimally allocating the budget to the available choices of interventions, our framework seeks to minimize the adverse consequences of the crisis including fatal overdoses and social and medical costs incurred by individuals with OUD. As input to the models within this framework, we utilized publicly available data and existing sources within the scientific literature, ensuring the accurate estimation of the key parameters. Through a comprehensive numerical analysis, we rigorously validated the framework’s effectiveness in generating results, demonstrating its capability to provide practical recommendations. Furthermore, we explored the results’ sensitivity to changes in the estimated parameters, revealing the adaptability and robustness of the framework under various scenarios.

We observed that when the ratio of the infected population is small (and therefore the susceptible population is large), the optimal approach is to focus solely on preventive interventions, aiming to manage the number of newly infected individuals and curb the epidemic. However, as the proportion of infected individuals increases, it becomes more effective to exclusively implement mitigating interventions, concentrating on reducing the number of fatal overdoses to control the epidemic. Our observations indicated that it is advantageous to allocate a considerable portion of the budget to preventive interventions, particularly when the rate of transition from susceptible to infected individuals rises. Interestingly, we found that even with an increase in the rate of transition from infected to removed individuals, it remains optimal to allocate a substantial budget to preventive interventions. Our numerical results also showed that it is optimal to allocate a larger portion of our budget to mitigating interventions when the ratio of fatal overdose to addiction marginal costs increases significantly.

Furthermore, we implemented a simulation model to assess the evolving trajectory of the epidemic under different budget allocation strategies during a simulated 10-year epidemic. Although the optimal policy dynamically adjusts allocations based on the epidemic’s state, scenarios of fixed 50-50 allocation policies and no allocation are also examined. The findings reveal that the optimal policy consistently surpasses fixed strategies, achieving notable cost savings of 29% compared with scenarios with no allocation and 12% compared with a fixed 50-50 allocation policy.

In summary, we believe that this research paves the way for informed and evidence-based policy-making, offering a vital tool for decision makers to effectively combat the opioid crisis and promote public health.

6.1. Limitations and Future Research

Although the literature includes a few quantitative approaches (such as simulation modeling) to evaluate the impact of a particular intervention on the opioid epidemic, to the best of our knowledge, no decision framework has been developed in this area. The proposed framework in this paper is the first milestone to mathematically formulate the budget allocation problem to address the opioid crisis as it relates to prescriptions. However, similar to any mathematical modeling approach, we had to make several simplifying assumptions (e.g., Poisson contacts and exponential communicability periods) to derive mathematically optimal solutions. In particular, we require that the number of new infections during each period is observable by the public health practitioner, and we assume that an infected individual interacts with the rest of the population only during the next period and then is effectively removed from the population. These assumptions can be relaxed to characterize optimal health policies using a partially observable Markov decision process (POMDP) (Sondik 1978, Molani et al. 2019). However, as the number of states required to model disease spread increases, MDP and POMDP models rapidly become computationally intractable, and even approximate dynamic programming methods may not efficiently solve such models. Nevertheless, such an extension is an attractive direction for future research. Also, it is essential to acknowledge that although assuming static transition rates and probabilities over one planning period simplifies modeling, the complex and dynamic nature of the opioid epidemic may involve fluctuations that this modeling approach does not fully capture. Future research should explore data-driven methods to better adapt to the dynamic nature of opioid epidemic dynamics over varying time intervals. Moreover, in this paper, we model the impact of interventions on the probability of transitions between compartments. Another possible direction for future research is to use policies that directly affect the number of infected or removed individuals. Furthermore, simplifying the prevalence of opioid pill diversion as a communicable disease by overlooking the heterogeneity of individuals in society concerning all socioeconomic factors may lead to inaccurate disease dynamics within the population. Future research may incorporate various subgroups (e.g., a combination of age groups and other demographics) under each compartment in the SIR model following the approach found in the literature on infectious diseases; see, for example, Enayati and Özaltın (2020).

As a final note, the numerical analysis performed in this paper is designed to illustrate the performance and capabilities of the proposed framework, and therefore the recommendation may not be directly applicable to policy-making. To successfully implement proposed dynamic health policies as a decision-aid tool in practice, the model’s parameters must be perpetually updated (which can be done using data analytics methods and machine learning techniques) and the epidemic state must be closely monitored. Even though dynamic optimization techniques are capable of handling noisy observations, inaccuracies in the parameter estimations and state monitoring processes may result in suboptimal or impractical policies for curbing the epidemic. Also, the notion of the “optimal” solution is ambitious given all the limitations in our analysis and the undeniable complexity of an issue like the opioid epidemic. In reality, no model can claim that an optimal solution exists for intricate public health issues such as the opioid epidemic. All these factors need to be considered when using the results of public health studies, including the current paper. Future research may translate these findings to real-world scenarios for use with real data to better understand the complexity of the opioid crisis and generate explicit state-dependent policies.