Modeling the Dynamics of Heroin and Illicit Opioid Use Disorder, Treatment, and Recovery.

Opioid use disorder (OUD) has become a serious leading health issue in the USA leading to addiction, disability, or death by overdose. Research has shown that OUD can lead to a chronic lifelong disorder with greater risk for relapse and accidental overdose deaths. While the prescription opioid epidemic is a relatively new phenomenon, illicit opioid use via heroin has been around for decades. Recently, additional illicit opioids such as fentanyl have become increasingly available and problematic. We propose a mathematical model that focuses on illicit OUD and includes a class for recovered users but allows for individuals to either remain in or relapse back to the illicit OUD class. Therefore, in our model, individuals may cycle in and out of three different classes: illicit OUD, treatment, and recovered. We additionally include a treatment function with saturation, as it has been shown there is limited accessibility to specialty treatment facilities. We used 2002-2019 SAMHSA and CDC data for the US population, scaled to a medium-sized city, to obtain parameter estimates for the specific case of heroin. We found that the overdose death rate has been increasing linearly since around 2011, likely due to the increased presence of fentanyl in the heroin supply. Extrapolation of this overdose death rate, together with the obtained parameter estimates, predict that by 2038 no endemic equilibrium will exist and the only stable equilibrium will correspond to the absence of heroin use disorder in the population. There is a range of parameter values that will give rise to a backward bifurcation above a critical saturation of treatment availability. We show this for a range of overdose death rate values, thus illustrating the critical role played by the availability of specialty treatment facilities. Sensitivity analysis consistently shows the significant role of people entering treatment on their own accord, which suggests the importance of removing two of the most prevalent SAMHSA-determined reasons that individuals do not enter treatment: financial constraints and the stigma of seeking treatment for heroin use disorder.

Introduction

A national crisis has emerged regarding opioid use disorder (OUD) (Vivolo-Kantor et al. 2018). Opioid overdose rates are on the rise and opioids are the primary cause of overdose deaths in the USA (Vivolo-Kantor et al. 2018; Jalal et al. 2018). In 2009, more than 20,000 people died in the USA by overdosing on opioids, including prescription opioids, heroin, and illicitly manufactured fentanyl; in 2019, the number of yearly opioid overdose deaths increased to nearly 50,000 according to the National Institute on Drug Abuse (Centers for Disease Control and Prevention 2020; National Institutes of Health 2019). Prescription opioid overdose, abuse, and dependence accounts for a total cost of 78.5 billion dollars a year reported by the Centers for Disease Control and Prevention (CDC). These costs are the result of elevated health care, drug abuse treatment, criminal justice, and loss of productivity expenditures (National Institutes of Health 2019; Florence et al. 2016). Other consequences of opioid abuse and dependence are exposure to sexually transmitted diseases, bacterial infections, and Neonatal abstinence syndrome (Hartnett et al. 2019; Centers for Disease Control and Prevention 2017; Haight et al. 2018; Volkow 2018). In addition, drug abuse is being closely linked to major depressive disorders and suicide attempts, which is now one of the increasing causes of death in the USA, according to the CDC (Dragisic et al. 2015; Center for Disease Control and Prevention 2016; Brook et al. 2002; Cole et al. 2019).

Opioids were commonly known in the past as naturally occurring substances derived from the opium poppy plant. They were thought to reduce the suffering of pain safely and effectively. Today, opioids now include the semisynthetic and fully synthetic drugs which invoke more intense, longer-lasting feelings of euphoria. Whether natural or synthetic, once in the bloodstream and traveled to the brain, they bind to -opioid receptors. This triggers the same reward system of pleasure and pain relief as do our body’s naturally occurring opioids called endorphins. The opioids activate the mid part of our brain generating feelings of pleasure from the discharge of dopamine in another part of our brain. This is known as the mesolimbic reward system (Kosten and George 2002; Lyden and Binswanger 2019; Incze and Steiger 2019; Veilleux et al. 2010).

Simultaneously, another part of the brain is remembering those good feelings of pleasure, specifically the details surrounding the event. Later, when faced with a similar situation, cravings for the drug taken are encountered. This is termed conditioned associations and it makes it very difficult for the user to not seek out that past feeling of pleasure. This leads to repeated use especially in the early stages. However, over time, repeated use switches from invoking those feelings of pleasure to avoiding the bad feelings of withdrawal. Another consequence of consistent opioid use is tolerance, which occurs when higher doses are needed to create the same previous sought after effects. The brain adjusted and now the individual feels right-minded when opioids are present, but abnormal when they are not, making withdrawal symptoms and cravings an issue. The implication of these complex brain processes then lead to the underlying causes for continuing use where a vicious cycle of repeated drug use has begun (Kosten and George 2002).

In recent years, opioid analgesics have been overprescribed and given their effect on the brain, this has resulted in an increased risk of OUD. This has also influenced an increase of heroin use where multiple users (4 out of 5 reported) have switched over from opioid pain reliever prescriptions because of lower cost and accessibility issues (Kolodny et al. 2015; Volkow 2018; Schuckit 2016; Connery 2015).

Fentanyl and other potent synthetic opioids on the black market have also fueled this problem. They are less expensive, more potent, and less costly to import and are either used to adulterate the heroin or replace it. The adulterated outcome of heroin mixed with fentanyl or other synthetics is unpredictable and dangerous (Volkow 2018; Williams et al. 2017; Spencer et al. 2019; Lyden and Binswanger 2019).

Many articles report a great need for OUD treatment, largely unmet, that signals a serious, widespread public health concern in the USA. For example, only about half of those individuals with heroin use disorder in the USA received treatment as stated in a 2014–15 study. Reasons mentioned include treatment not easily accessible, shortages of trained healthcare staff, insurance coverage issues, limited policy changes, limited financing of care, and limited means of quality care (Ghitza and Tai 2014; Mojtabai et al. 2019; Kolodny et al. 2015; Volkow 2018; Williams et al. 2019; Connery 2015; National Institutes of Health 2021).

Other than the previously mentioned obstacles, another barrier for treatment and limited access to care may include that the public’s view of drug abuse and dependence is stigmatized as opposed to being viewed as a chronic life-threatening disease in need of assistance. As a result of the stigma, in the past, the focus was on an abstinence-based treatment plan. Currently, there are three medications approved by the Food and Drug Administration (FDA) proven to reduce future overdoses and illicit drug use when combined with counseling and behavioral therapies. However, there still exists some reluctance on using these medications to treat OUD. The three medications are methadone, buprenorphine, and extended-release naltrexone. The combination of medication with counseling and behavioral therapies is called medication-assisted treatment (MAT) (Mojtabai et al. 2019; Volkow 2018; Coffa and Snyder 2019; Williams et al. 2019; Lyden and Binswanger 2019; SAMHSA 2021; National Institutes of Health 2021).

These medications remain underused where only a minority receives any treatment (including non-medication routes) and even a smaller amount receive MAT. Among treatment programs in the private sector, less than fifty percent offer opioid based medications and of these programs only thirty-three percent of patients are prescribed them. Therefore, many of the 2.4 million in the USA with OUD do not receive any MAT. To diminish the US OUD overdose epidemic, these barriers and misunderstandings for using these treatment steps must be tackled. OUD treatment is important to decrease the mortality of millions of Americans at risk of opioid-related overdoses. As a result, public health authorities are increasing efforts to integrate such treatment (Mojtabai et al. 2019; Kolodny et al. 2015; Volkow 2018; Williams et al. 2017; Lyden and Binswanger 2019; National Institutes of Health 2021)

A tool that can be used for understanding the complex issues surrounding OUD and illicit OUD, its treatment options, and methods for decreasing relapse, is a mathematical model. Mathematical models are very important to gain understanding of disease epidemiology. Using the spread of OUD viewed as the potential contagion, we can then use a mathematical model to describe the spread of OUD and the dynamics underlying those patterns that can best inform and assist policy makers in targeting prevention and treatment resources for maximum effectiveness (Bailey 1975; Anderson and May 1992; Murray 2007; Brauer and Castillo-Chavez 2012).

Studies of mathematical models on drug use have been previously conducted. White and Comiskey (2007) divided the population into susceptible, current, and in-treatment drug users for heroin addiction. A basic reproductive number, representing how many new users is produced per each current user, was found. A sensitivity analysis pertaining to control efforts was performed, which found that decreasing the transmission term of the contagion showed higher significance than increasing the proportion of users who enter treatment. The authors also found a condition where a backward bifurcation exists, which means that an endemic equilibrium may exist even when the reproductive number is less than one. Therefore, extra efforts would be needed to drive down the epidemic. Also noted in their model is the inclusion of enhanced death rates for the current users and users-in-treatment classes.

Some model studies of the White and Comiskey article were considered by other authors Mulone and Straughan (2009), Wang et al. (2011), Muroya et al. (2014), Ma et al. (2017) including Wangari and Stone.

Wangari and Stone (2017) had the added compartment class of individuals who left treatment but are not using. They also added a saturation term to deal with the shortcomings of the healthcare system when too many people seek treatment at the same time. They found when this saturation parameter was above a critical threshold, backward bifurcation existed. Their sensitivity analysis concluded that this parameter was of high importance in feeding the epidemic. The effective contact rate and relapse rate from treatment are other parameters they found with high sensitivity.

Additional models branched off of the White and Comiskey as well, including the distributed time delay (Liu and Zhang 2011; Liu and Wang 2016; Fang et al. 2014; Huang and Liu 2013; Samanta 2011) and the age structured models (Fang et al. 2015b, a; Wang et al. 2019).

Caldwell et al. (2019) implemented and analyzed a Vicodin epidemic model that focused only on the population of people who were prescribed Vicodin. They also included a global sensitivity analysis to show that preventative measures over treatment efforts are more successful for reduction of misuse.

Battista et al. (2019) proposed a model that added an opioid prescription drug user class where a potential user can become addicted through either the use of prescriptions, legally or illicitly, or through contact with another addicted person; they included a treatment class as well. Mathematical analysis was performed, showing that an addiction-free state cannot be attained without controls over prescriptions. Their sensitivity analysis showed that prevention, followed by vigorous treatment, may result in a low status of endemic misuse.

We propose an “illicit opioid use disorder” (IOUD) model to describe the role that black market opioids such as heroin, fentanyl, and other synthetic opioids play in the current opioid epidemic. Our model does not include a prescription class but will be extended to do so in future work. Thus, our proposed IOUD model can be viewed as what might happen if opioids were outlawed or, perhaps, severely restricted.

Novel to our model is the inclusion of a recovered class that does not allow for a past user to ever be considered as a nonusing susceptible individual in the future. Therefore, we must allow for relapse from both the recovered and treatment classes. According to Kosten and George (2002), repeated and prolonged drug use modifies physiological brain functions. Moreover, alternating between abstinence and withdrawal creates a “changed set point” model. Within this model, healthy dopamine (DA) transmitter activity is permanently altered by use of opioids. This effectively changes the natural baseline of DA tolerance in addicted individuals. Another model called the “cognitive deficits model of drug addiction” explains that damage to the prefrontal cortex may result due to habitual use. This further reduces judgment capacity and impulse constraint. The challenges arising from this neurobiological deterioration permanently increases the risk of relapse. Since chronic opioid use results in these brain transformations, cravings may be produced, causing a recovered individual who is no longer opioid dependent to relapse, following months or years of their abstinence (Kosten and George 2002; Kolodny et al. 2015).

Our IOUD model also considers a treatment class with a saturation term that slows down the rate at which people receive treatment, due to the previously mentioned barriers. We will see that both of these extensions play a role in the dynamics of the system.

Model Formulation and Basic Properties

Our proposed model assumes a homogeneous mixing of the human population. The total population at time t is denoted by N(t) and is divided into four mutually exclusive compartments as follows: susceptibles S(t), individuals with illicit OUD I(t), individuals in a treatment facility T(t), and recovered individuals R(t). Thus, ; see Fig. 1 for how individuals can move between compartments.

Fig. 1

Compartmental flow diagram of the illicit opioid use disorder (IOUD) model. S represents susceptible individuals, I represents individuals with illicit OUD, T represents those in specialty treatment facilities, and R represents recovered individuals. R is considered distinct from S due to an increased potential for relapse. The factor models the decreased rate of entrance into the T class due to limited access of care in specialty treatment facilities

The susceptible (potential individuals with illict OUD) class describes the number of the population who either have never used opioids or have used illicit opioids but never been considered to have illicit OUD. The susceptible population is increased by the constant recruitment rate, . A constant for recruitment was chosen because it will lead to an asymptotically constant population size as opposed to a linear one which might possibly lead to exponential growth or decay.

The IOUD class describes the individuals who have illicit OUD. OUD according to the Diagnostic and Statistical Manual of Mental Disorders, 5th Ed. is defined as the use of opioids leading to a precarious situation of repeated use and as a result, at least two destructive symptoms occur within a year period. These include problems such as strong, persistent cravings, failure to perform societal and personal obligations, increased physical endangerment, and an increased tolerance to opioids. A full list can be found in the manual (Edition et al 2013).

Someone who takes opioids illicitly a few times, in some kind of social circumstance (a few parties, music festivals, etc), but never has the kind of constant use that would result in the patterns discussed above would not be considered as having illicit OUD. Thus, this individual would not be considered in the IOUD class but will remain in the susceptible class.

This population class is considered infectious and as a consequence of interacting with individuals with illicit OUD, a susceptible individual may develop tendencies that could lead to illicit OUD. The value is the transmission rate of that interaction resulting in a change of class from S to I. In this way, the susceptible population may flow to the IOUD class.

There are multiple ways that individuals transition out of the IOUD.

The treatment class describes individuals with illicit OUD who are in a specialty treatment facility. Individuals with illicit OUD may decide to leave for the treatment class on their own at a rate of , or through the influence of a recovered individual or someone from the susceptible population; the last two interaction rates are , , respectively. These individuals may relapse back to the IOUD class by relapse rate or at the end of their treatment they may flow to the recovered class at rate .

From the yearly statistics from the Substance Abuse and Mental Health Services Administration (SAMHSA) published in the annual National Survey on Drug Use and Health (NSDUH), specialty treatment facilities (our T class) include hospitals (inpatient only), rehabilitation facilities (inpatient or outpatient), or mental health centers. In contrast, non-specialty treatment facilities include emergency room, private doctor’s office, self-help group, and prison/jail (Center for Behavioral Health Statistics and Quality 2020).

The recovered class describes all the individuals who either completed specialty treatment (i.e., went from T(t) to R(t)), or those with illicit OUD who quit on their own or with the help of a non-specialty treatment facility (i.e., went from I(t) to R(t) in either case).

Since illicit opioid use is a chronic condition (Kosten and George 2002), individuals remain in the recovered state unless they relapse which may be on their own at a rate of or as a consequence of interacting with an individual in the IOUD class at a rate of .

There is a removal from each class as the natural death rate , whereas the IOUD class, I(t), has an additional removal rate of . With this, the added component due to illicit OUD overdose death (Seth et al. 2018), an overall computed death rate for the individuals with IOUD would be .

Model Equations

The model is given by the following deterministic system of nonlinear differential equations:where and all parameters are nonnegative.

We use a saturation treatment function b(T) to modify the flow of individuals with illicit OUD to treatment, where the parameter models a saturation of availability of specialty treatment facilities. This limits the amount of individuals with illicit OUD that can go into specialty treatment facilities due to the limited access of care discussed previously in the introduction.

A description of variable and parameter values are listed in Table 1.

Table 1

Description of variables and parameters of the model

Variable	Description
S(t)	The total number of people who are susceptible at time t
I(t)	The total number of individuals with illicit OUD (for the first time and from relapse) not in specialty treatment or recovered at time t
T(t)	The total number of individuals in specialty treatment at time t
R(t)	The total number of individuals who have either completed specialty or non-specialty treatment or “quit cold turkey” at time t

The basic properties of the IOUD model were explored and those results can be found in Appendix.

Data Explanation and Parameter Estimation

Our model considers illicit OUD, treatment, and recovery, as well as overdose deaths. CDC data exists for overdose deaths due to synthetic opioids (primarily fentanyl) as well as heroin (sometimes in combination with synthetic opioids). However, the data from SAMHSA on illicit OUD and treatment is limited to heroin likely because the presence of synthetic opioids is a relatively recent phenomenon. Thus, for the purpose of comparing our model to data, we consider only heroin use or heroin use with synthetic opioids (both considered by SAMHSA) but do not include additional synthetic-opioid-only use. We consider a generic US city of approximately 200,000 people and scale the national data from the number of individuals with heroin use disorder (HUD), the number of individuals with treatment, and the number of overdose deaths to a city of this size by taking into account the increasing yearly US population. This allows us to consider a nearly constant population size as we analyze the dynamics of the model. For example, heroin use disorder (HUD) in 2002 in the top right graph of Fig. 2 is 148.97 = (214,000/287.3E+06) 200,000; see Table 2 for similar yearly numbers. We found CDC data with the number of deaths due to heroin or heroin mixed with synthetic opioids (third column of Table 2), which is dominated by fentanyl, as well as death from synthetic opioids alone (second column). SAMHSA data was found for the number of individuals with HUD, with the NSDUH counting those with HUD in the past year (fifth column, relates to our variable I). SAMHSA data is available for treatment in a specialty facility in the past year (sixth column, relates to our variable T) and also in a general treatment center in the past year (not presented). Data for "in specialty treatment for HUD" was only presented in the 2014–2017 SAMHSA survey results. In order to scale the other years to give an approximate number in specialty treatment facilities coming from HUD (our variable T), we looked at the ratio of specialty treatment from HUD to the specialty treatment data in the 4 years when it was available. The factor 0.6874 is the average of the ratio. The specialty treatment  0.6874 is labeled in the sixth column with an asterisk, and also given without error bars in the graph. The treatment data, similar to the HUD data, counted individuals in a specialty facility in the last year. This is the data presented in our graphs with the raw data given in Table 2. The error bars in the graphs represent the standard error given in the SAMHSA data (not presented in the table).

Fig. 2

Model output compared to data scaled to a population of 200,000 by taking into account the yearly US population values. (Top Left): CDC data for overdose deaths in HUD class due to heroin, obtained as 0.8 (total overdose deaths due to heroin), presented as red curve with diamonds compared with model output as blue curve with circles. (Top right): SAMHSA data for in “HUD in past year,” with error bars when given. Model approximation is the blue curve with circles, calculated with instantaneous model variable I (solid, cyan curve immediately below) averaged over each year and added to the “correction” for those that left and also possibly returned to I (see text). (Bottom right): SAMHSA data for in “specialty treatment in past year coming from I,” with error bars when given. Model approximation is the blue curve with circles, calculated with instantaneous model variable T (solid, cyan curve immediately below) averaged over each year and added to the “correction” for those that left and also possibly returned to T (see text). The bottom 2 curves in the right panels signify those who left I and T over the year presented with dash–dot curves and the corrected quantities of those who left the respective classes are presented with dotted curves. These last two quantities sum to give the solid curve with circles that we compare with the SAMHSA data. (Bottom left): Data-derived and least squares fit for . Asterisks and x-marks are calculated from data (see text and equation (3)) with blue x-marks used to obtain the horizontal (constant) line and black asterisks used to obtain the nonzero sloped line; both lines are calculated with a least squares fit (Color figure online)

Table 2

Data for USA, 2002–2020. The number of overdose deaths for 2002–2020 are from the CDC (Centers for Disease Control and Prevention 2020). US population comes from United Nations (2019). Use disorder and specialty treatment data come from SAMHSA’s NSDUH (Center for Behavioral Health Statistics and Quality 2020, 2018, 2016, 2015, 2014; Lipari and Hughes 2015; Center for Behavioral Health Statistics and Quality 2013; Substance Abuse and Mental Health Services Administration 2011, 2010, 2008, 2006). The derivation of values in the column -data are given in (3) where we used (HUD class data in year)0.903 to estimate average number with HUD during the year. The values in the column -fit are obtained from (4)

Year	Deaths due to overdose		US population	HUD in last year	Specialty treatment in last year from I	δ-data	δ-fit
	Synthetics	Heroin
2002	1295	2089	287.3E+06	214,000	Not available	0.008648	0.008089
2003	1400	2080	289.8E+06	189,000	Not available	0.009750	0.008089
2004	1664	1878	292.4E+06	270,000	107,200^∗	0.006162	0.008089
2005	1742	2009	295.0E+06	227,000	130,600^∗	0.007841	0.008089
2006	2707	2088	297.8E+06	324,000	259,100^∗	0.005709	0.008089
2007	2213	2399	300.6E+06	214,000	138,200^∗	0.009932	0.008089
2008	2306	3041	303.5E+06	283,000	156,000^∗	0.009520	0.008089
2009	2946	3278	306.3E+06	369,000	221,300^∗	0.007870	0.008089
2010	3007	3036	309.0E+06	361,000	188,300^∗	0.007451	0.008089
2011	2666	4397	311.6E+06	426,000	200,700^∗	0.009144	0.009255
2012	2628	5925	314.0E+06	467,000	201,400^∗	0.011240	0.01156
2013	3105	8257	316.4E+06	517,000	246,800^∗	0.014149	0.01387
2014	5544	10,574	318.7E+06	586,000	270,000	0.015986	0.01618
2015	9580	12,989	320.9E+06	591,000	242,000	0.019471	0.01848
2016	19,413	15,469	323.0E+06	626,000	235,000	0.021892	0.02079
2017	28,466	15,482	325.1E+06	652,000	358,000	0.021037	0.02310
2018	31,335	14,996	327.1E+06	526,000	291,500^∗	0.025258	0.02540
2019	36,259	14,019	329.1 E+06	438,000	321,000^∗	0.028356	0.02771
2020	56,883	13,058	331.0E+06	Not available	Not available	Not available	0.03002

Specialty treatment0.6874 because specialty treatment from HUD only asked in 2014–2017 SAMHSA surveys. The factor 0.6874 is the average of the ratio of specialty treatment from HUD to specialty treatment in the 4 years when data is available

Our variables I and T are instantaneous in time, whereas the SAMHSA data gives those in the respective classes in the past year. In the case of comparing our model output with the SAMHSA treatment data, individuals that were in treatment in the past year could (i) be currently in T, (ii) have relapsed and went from T back to I, or (iii) have successfully completed treatment and moved from T to R in the past year. In the case of comparing our model output with the SAMHSA HUD data, individuals with HUD could (i) be currently in I, (ii) have moved to treatment (I to T), or (iii) moved directly from I to R in the past year. Thus, we additionally need to keep track of the number of individuals who left each of these classes each year. We further correct our model output with a small discount for those that went back again (after having left and thus should not have been discounted). In the data matching plot, we present the model output, those that left I and T over the year, and a correction of those who left I and T over the year but then went back (estimated with , and , respectively). Those who left I and T over the year are presented with dash–dot curves, the corrected quantities of those who left the respective classes are presented with dotted curves, and the variable output of the class is a solid curve with no circles. These last two quantities sum to give the solid curve with circles that we compare with the SAMHSA data.

We were able to come up with reasonable estimates for many of the parameters based on the literature. We used from Wangari (2017), Wangari and Stone (2017) where it was assumed that the average person’s lifespan is 80 years old and thus . From Battista et al. (2019), Battista et al. (2019), we obtained an approximate range of as 0.1 to 0.4 Weiss and Rao (2017). We estimated to be in the range 0.4 to 0.9 from Smyth et al. (2010), Bailey et al. (2013), Weiss and Rao (2017). We set so that the population in the heroin-free model reaches 200,000 for the assumed and for . For parameters , , and , the entry to treatment rates, we used the range from Battista et al. (2019) and Wangari and Stone (2017), Zhang and Liu (2008). Both models had only a linear term from their addicted class to treatment, whereas our model has one linear term and two nonlinear terms between the comparable classes. Considering , we set estimates for , , and . Similarly, we found a rate from recovery back to HUD from the literature in a study by Gossop et al. (1989). We estimated to be in the range 0.1 to 1/3 with and significantly bigger than . We used , . The parameter for going directly from I to R, either “quitting cold turkey” or quitting through a general (non-specialty) treatment facility, was estimated to be in the range .01 to 0.2 (Wangari and Stone 2017).

The parameters and were difficult to determine, so we did parameter estimation with them as well as for and (where we used the above range from the literature for the latter two) (Banks et al. 2013; Banks and Bihari 2001; Cintrón-Arias et al. 2009). While we were able to approximately match the data for I and T, we were not able to come close to matching the overdose death data for a fixed , which increased significantly from 2010 through 2016, even allowing for a possible change in parameters in 2010 (see derivation below and Table 2).

We now consider in more detail. By definition, we have thatFor the numerator, the CDC data gives total yearly overdose deaths due to heroin, irrespective of whether an individual was with HUD or not (Centers for Disease Control and Prevention 2017). We note that the paper by Battista et al. on prescription opioids estimates a discount factor from the literature on what portion of opioid deaths were from someone addicted to opioids to address this analogous problem (Battista et al. 2019). In our current discussion that focuses on HUD, we did not find any comparable statement in the literature regarding the percentage of individuals who die from an overdose of heroin that were in the HUD class (in contrast to those who die from a heroin overdose but are “casual users”). We estimate that 80% of the heroin overdose deaths are from individuals with HUD as a first approximation that can be corrected if data becomes available. For the denominator, we need to estimate the average number of individuals with HUD during the year for a given year since the SAMHSA data gives the cumulative number of those with HUD in the past year (Center for Behavioral Health Statistics and Quality 2014, 2015, 2016, 2018, 2020; Lipari and Hughes 2015; Center for Behavioral Health Statistics and Quality 2013; Substance Abuse and Mental Health Services Administration 2006, 2008, 2010, 2011). In comparing the model output variable I with the model calculation to give the number in I in the past year (both with results from the parameter estimation), we observed the graphs were shaped similarly (solid cyan curve immediately underneath solid blue curve with circles in the top right graph of Fig. 2). We thus calculated the ratio of the average of the model output I over the past year to the model output I in the past year (described above) for each year and found its average value to be 0.903. In our calculation of , we thus estimated the average number in the HUD class over the year as the SAMHSA data for those individuals with HUD in the past year multiplied by 0.903. Thus, we calculate the yearly values asIn examining the data-derived yearly values of , we observed a significant year over year increase starting in 2012 through 2019; see the bottom left subgraph in Fig. 2 and Table 2. The 2020 overdose deaths were published recently by the CDC. During the revision of this manuscript, SAMHSA published the 2020 data for HUD; however, they changed the criteria for classifying an individual as HUD, thus making the 2020 data that were released not obviously compatible with the data from 2019 and earlier. Thus, we are not able to include the 2020 value in our parameter estimates. When plotted versus time, the -values follow a piecewise linear function (2002–2019), as shown in the bottom left subgraph in Fig. 2. We use the corresponding piecewise function obtained with a least squares fit (last column of Table 2) in our model calculations:where we present this number of significant digits to have agreement to six significant digits when the function switches branches. (In our computations, additional digits are kept.) Incorporating this piecewise function for into our parameter estimation, our baseline values areWe choose our initial conditions to approximately match the scaled data from : . While our model is for illicit opioid use and not just heroin, only heroin data is available for the I and T classes, and that is the data we use to fit our model. The data match is provided in Fig. 2.

Given the myriad of ways in which we varied parameters to try to match the data, we conclude that we cannot fix at a constant but must vary it according to the yearly data if we are to obtain agreement of model output with data. This increase in over time corresponds to a higher overdose death rate per individual in the HUD class. Given the agreement of data and model output, necessitated by an increasing , we interpret this deadlier as follows: the increase in number of heroin overdose deaths was driven by the increase in the prevalence of fentanyl in the heroin supply, as no increase in the HUD class (either seen in the SAMHSA data or our model output) could account for such a drastic increase with a fixed . Fentanyl first appeared in 2007, is 100 times more potent than heroin, and its prevalence is well known to have been getting greater in the heroin supply and illicit opioid use (Worth and House 2018; United States Drug Enforcement Administration 2021).

We note that the number of overdose deaths and the number in the HUD class both have decreased over the last 3 years, whereas our model continues to increase. (Interestingly, the data-derived -value for 2018, 2019 still increases in spite of this decrease.) As shown in the second column of Table 2, deaths due to synthetic opioids have skyrocketed. Numerous articles suggest that heroin users may be switching to synthetic opioids, but SAMHSA data does not keep track of synthetic opioid use explicitly and only in the last few years has considered illicit opioid use. A recent article from the RAND corporation states that “Cheap, accessible, and mass-produced synthetic opioids could very well displace heroin, generating important and hard-to-predict consequences” (Pardo et al. 2019).

Steady-State Analysis

Traditional epidemiological language uses the phrase “disease-free equilibria (DFE)” to describe the absence of the given disease. Our I class consists of those active users with illicit OUD. Thus, we will consider an “illicit opioid use disorder-free equilibria” (IOUDFE) that we will shorten to “disorder-free equilibria” (DFE) for convenience. We are interested in the DFE and its stability for (1). With , gives . Hence, the DFE of our IOUD model is

For the ensuing analysis, we consider a fixed so that the death rate due to overdose remains constant at some level (e.g., at its 2020 value). We tried to analyze the equilibria of the system using the local stability analysis and the Routh–Hurwitz criterion (Wirkus et al. 2017; Edelstein-Keshet 2005), but due to the complexity of the expressions we were not able to obtain any useful information.

Calculating the Basic Reproductive Number

The basic reproductive number, , is a quantity that represents the expected number of new infections produced per infected individual during their infectious period when a disease is introduced into a susceptible population. In the context of our model, it determines the additional number of new individuals with illicit OUD that each individual with illicit OUD will produce before entering treatment or recovery.

We will find the for our model (1) by using the next generation method as presented in Van den Driessche and Watmough (2002) and also by considering a heuristic derivation (see, e.g., Van den Driessche and Watmough 2002; Wangari and Stone 2017); both agree.

We restate (1) with the b(T) saturation term explicitly in the equations:For the heuristic derivation of , as presented by Van den Driessche and Watmough (2002) we observe that we can cycle in and out of the IOUD class, either through treatment or recovery or both due to relapse of individuals in the treatment or recovery classes (Van den Driessche and Watmough 2002; Wangari and Stone 2017). The average time an individual spends time as an opioid user in I without treatment is

The fraction of surviving I and moving to treatment is and the fraction of surviving I and moving to recovered is Now, the fraction of the surviving opioid users in T returning to I is seen to be , while the fraction of the surviving opioid users in R returning to I directly is We set , which defines going from IOUD to treatment and back to IOUD; and , which defines going from the IOUD class to recovered and back to the IOUD class.

In addition, we now have the fraction of surviving opioid users moving to treatment, then to R, and then to I:Our new expression for all possible combinations of multiple passes will now beAs this is a geometric sequence, we can write its sum as . Substitution for , , , and multiplication by gives uswhich can also be rearranged asIn this latter form, we see that entering treatment, either of ones own accord, , or through the interaction with a susceptible individual, , as well as recovering on one’s own, , (all increasing) will result in a lower value of . Decreasing the transmission rate, , or increasing the illicit OUD overdose death rate, , would also result in decreasing . The cycling that can occur between the I, R, and T classes makes the remaining parameters less obvious.

This expression for is the same as that obtained via the next generation method as we now show (Van den Driessche and Watmough 2002). The method requires that we identify “new infections” and “infected” compartments. We note that changes of the individual from T to I and R to I are not considered to be new infections, but rather the movement of an infected individual through the different compartments. According to the definitions of and , we computeandClearly I is an infected compartment as it holds those individuals with IOUD. Due to the structure of the equations and the mathematical method, T and R must also be considered as infected compartments because individuals can go from R or T into I without interaction because of the non-contact rates between them and the I class. In terms of the biological justification, T and R are infected compartments since opioid use may result in brain transformation with cravings that may be invoked, leading to relapse of an individual in treatment or a recovered individual (Kosten and George 2002). Thus, our infected compartments are I, T, and R giving .

According to the definitions of F and V and using our previously calculated DFE, we obtainandThe calculation of results in only one nonzero eigenvalue that contains only nonnegative parameter values. This maximum eigenvalue of gives us the same expression for as from (8).

One interesting observation is the absence of , , and from the expression. Since the interpretation of is often stated as one infected introduced into an entirely susceptible population, this would suggest that limited access to special facilities (modeled via ) will not play a role initially and the size of will be too small for or to have any effect.

Endemic Equilibria

We will now determine the existence of non-trivial endemic equilibria of the system. We will be particularly interested in the situation of a backward bifurcation, which is characterized by a stable endemic equilibria existing even when . In the region of bi-stability, both the endemic equilibria (EE) and the DFE exist and are stable. We begin by considering the case of no saturation, , so that the situation of limited availability in specialty treatment facilities does not occur:We will show that this case does not permit the existence of a backward bifurcation for but does permit one for large enough . We can obtain an equation in only the variable as follows. We set and solve for :We set and solve for ; see (17):We set and and solve for and :where and are defined as above.

We substitute into . After simplification, we obtain an equation of the form where B is a complicated expression involving the parameters as well as and andWe observe from (8) that does not depend on or since those factors are not present in . Thus, the presence of the factor ( suggests that altering or may affect the sign of the denominator. We first set , so that the denominator is always positive, and thus, we focus only on roots of the numerator.

From inspection, we observe that B is a cubic expression in without a constant term. Thus, our cubic expression for the roots of has the formwhereandThus, this c term from (10) is positive when and it is negative when . We will use this information to interpret whether or not it is possible to have a backward bifurcation when by using Descartes’ Rule of Signs. We know that when our c term must be negative. We also know that our a term in the quadratic in (10) must always be negative. According to Descartes’ Rule of Signs, there can be two or no positive real roots if and no positive real roots if . Using our baseline parameters discussed later with the modification that and , we observe that and the roots of (10) are positive. This is confirmed in the full system, and thus, we conclude that we can have a backward bifurcation for for sufficiently large (approximately for the given parameter values). We note before proceeding that the value for is 2 times its current estimated value; in contrast, is the value that fit the data, and thus, the value of nonlinear relapse rate needed for a backward bifurcation is at least 120,000 times greater than this and thus likely unrealistic.

We now keep and consider with . The denominator may become negative for sufficiently large . Trial and error shows that we can find roots of B/C that are positive. However, substituting these values into the full system yield negative values for some of the other variables. Following Battista et al. (2019), Castillo-Chavez and Song (2004) as shown in appendix, we show that this case cannot have a backward bifurcation. Thus, without saturation, we can have a backward bifurcation for an unrealistically large , the nonlinear relapse from R to I, but cannot have a backward bifurcation when the nonlinear relapse rate is zero.

Let us now look to analyze the equilibria when , i.e., we will include the saturation term. We will show that a critical value exists above which a backward bifurcation is permitted. Of particular note is that this critical value for is within a reasonable range and the value of could be 0 or its baseline value.

Proceeding in a more straightforward manner complicates things immediately due to the large algebraic expressions. We tried to simplify the saturation term through a Taylor series expansion for small but that approach did not work. Instead, we allow the system, whose total population is governed by , to reach its steady-state population level, , given by when . The resulting limiting system is as follows:where

We again try to obtain an equation involving only parameters and the variable . We proceed as before by solving for by setting and then plugging that result into and . Next, we solve for and in terms of simultaneously by setting and . We plug and into . The resulting equation, which we will refer to as (), is in terms of the variable . This is all done using Maple and is not presented here due to its length.

Obtaining a general expression for when a backward bifurcation occurred yielded pages of expressions that were too complicated to analyze. We thus chose to focus on three parameters, , , and the parameters addressing overdose death, saturation, and transmission, respectively. We extrapolate the -values from Fig. 2 as well as calculate the corresponding effective reproductive number to determine when the DFE and EE will be stable; see Fig. 3. The results that we now present use realistic parameter values based on data through 2019 and presented in (5) to give stability curves in terms of the overdose death, saturation, and transmission parameters. We can observe regions in the –– parameter space that correspond to the EE stable (only), both DFE and EE stable (bi-stability), and the DFE stable (only). In this latter situation, the EE no longer exists biologically with only the DFE persisting and stable. While this is clearly not a desirable situation, the increase in fentanyl in the heroin supply makes this scenario a potentially realistic one that needs consideration.

Fig. 3

(Left): Extrapolated -values. The blue x-marks and black asterisks are from the overdose data and are the same as in the bottom left panel of Fig. 2. The line obtained with a least squares fit of the data from 2011–2019 and given in (4) is extended to 2038. The labeled -values in 2020, 2029, and 2038 are from extrapolation (16) using the best fit line. (Right): The effective reproductive number, , is plotted as the solid black curve using the baseline values of the parameters from (5) and the extrapolated -values from the best fit line. Just above the curve, is plotted as a dashed blue curve; this close approximation is expected given that (Color figure online)

We leave , , and as parameters and substitute the remaining parameter values from (5) into () to obtain an equation in the parameters , , and and the variable :where the coefficients are given in appendix and the subscript refers to the power of . We eliminate the variable by simultaneously solving (15) and the derivative of it, thus requiring the condition for a saddlenode bifurcation. This results in a new equation in terms of , , and that is pages of output in Maple. However, we can plot this implicit equation numerically and present this three-dimensional –– surface with five cross-sectional subplots; see Fig. 4. The years given in Fig. 4 correspond with those shown in Fig. 3. For large , we have the situation where b(T) is very small, which is not allowing people to go into treatment due to a lack of availability in specialty treatment facilities.

Fig. 4

Regions of stability for equilibria. Top panel (left and middle): in the – plane with fixed at .09, the solid blue horizontal line corresponds to the constant for which =1 and below this line only the EE is stable; above this line for large enough is the curve that separates the region of bi-stability from where only the DFE is stable. Top panel (right): in the – plane with fixed at .0313, the two lines separate the regions of (i) EE stable (only), (ii) bi-stability of EE and DFE, and (iii) DFE stable (only); the upper line corresponds to =1. Right panel (middle and bottom): In the – plane with fixed at .0577 (its extrapolated 2032 value), the solid red horizontal line corresponds to the constant for which =1 and above this line only the EE is stable; below this line for large enough is the curve that separates the region of bi-stability from where only the DFE is stable. Bottom left panel: the previously described curves are put together in the three-dimensional –– space. The dots with years correspond to -values from the extrapolated -curve with all other parameters fixed at their baseline values from (5) with the color magenta corresponding to EE stable (only), blue corresponding to the region of bi-stability, and black corresponding to DFE stable (only) (Color figure online)

The presence of a backward bifurcation yields a region of bi-stability when . This means that we will have two asymptotically stable equilibria, the EE and the DFE, and which one a solution approaches simply depends on the initial conditions. Above the plane in the 3-d subplot, only the EE is stable. Below this plane, there is a range of parameter values where we may either have bi-stability or have the DFE as the only stable equilibrium.

For a given set of parameters, there is a critical that is required for bi-stability and a backward bifurcation. There is an inverse relationship between the saturation term, , and an availability of specialty treatment facilities. Thus, a lack of availability of specialty treatment facilities that occurs when can give rise to a situation in which the epidemic persists even though the conditions are such that . See Fig. 5.

Fig. 5

Backward bifurcation plots. The blue curves correspond to stable biologically relevant equilibria and the red curves correspond to unstable biologically relevant equilibria. This demonstrates the difficulty there may be in getting rid of the epidemic once it has taken hold. Top panel: is fixed at .0531, its extrapolated 2030 value; is varied to change . Bottom panel is fixed at .09; is varied to change with (2029 on its extrapolated curve) corresponding to . All other parameter values are from (5). The middle column differs from the first column in scale only (Color figure online)

Sensitivity Analysis

For our sensitivity analysis, we run the model from 2002 to 2020 using the parameters in (4)–(5) and then use the resulting 2020 model output values as our initial conditions. We use the baseline parameters given in (5) that generated this data match and consider two scenarios for : (i) assume that for , which we interpret as the fentanyl levels being kept at their 2020 levels, and (ii) assume that is defined by extrapolation based on its least squares fit line given in (16), which we interpret as the fentanyl levels in the heroin supply increasing. In both cases, we consider what happens in 2030 for the sensitivity. In the first scenario, is a constant and will be a parameter in our sensitivity analysis. For the second scenario, we explicitly rewrite in (4) in a form that allows for a vertical shift at 2020 as well as a potential shift in the slope of the line by at 2020. This is accomplished via the extension of the least squares fit line in (4), shown in the left panel of Fig. 3, with and written aswhere a percent change in b changes through a vertical shift by the same percent change of its 2020 value and a percent change in m changes the slope of by the same percent change. With baseline values of and from (4), we will examine these 2 additional parameters in our sensitivity analysis for the second scenario (McLeod et al. 2006).

In order to determine the sensitivity of the system to the input parameters, we perform a sensitivity analysis using partial rank correlation coefficient (PRCC) method (Marino et al. 2008). The PRCC method only applies when there is a monotonic relationship between the model parameters and the output values against which sensitivity is measured. We performed monotonicity checks for all our parameter and initial values and concluded that there is a monotonic relationship.

For our system, we consider the parameter values obtained through parameter estimation and given in (5) as the baseline parameter values. When we consider the extrapolated function for , we observe from Fig. 4 that the value puts the system in the region of bi-stability.

We let the parameters and initial conditions vary from their baseline values in 2020.

Discussion of the PRCC Values

We present the sensitivity of our variables S, I, T, and R to the parameters of the system in plots and tables in appendix and focus here on variables that may be of more interest to healthcare professionals and policy makers: number of those entering I (HUD) for the first time (yearly new HUD), the yearly number of relapses from T, the yearly number of relapses from R and heroin-related deaths.

While none of these are the variables in our original system, all can be calculated by keeping track of components that contribute to changes in our model variables.

We consider two graphs for each case corresponding to the sensitivities in 2030 for the constant death rate () vs. the variable death rate (16).

In describing the sensitivity results we will refer to a PRCC value of 0.85 or higher as “highly significant,” a PRCC value of 0.70–0.84 as “significant,” values of 0.55–0.69 “somewhat significant,” values of 0.45–0.54 as “slightly significant,” values of .40–0.44 as “borderline significant,” and under .40 as “not significant.”

As can be seen in the tables, some of the initial conditions may show up as significant or highly significant. We fit the 2002–2019 data to baseline parameters with the model output final (year 2020) values forming the initial conditions for our PRCC analysis. While S(0), I(0), T(0), and R(0) cannot really be changed, having somewhat different data (e.g., more accurate data) could represent the importance of their significance. Additionally, for the parameter (the natural death rate of the general population), regardless of its significance, it is not a parameter that can be altered since it is the natural death rate. Therefore, we would not focus on it either because it is something we do not have control over. For the following, only the parameters that we have control over will be discussed.

For the following variables’ PRCC results that will be discussed, it can be seen that the graphs at the end time of 2030 are similar for the constant death rate () vs. the variable death rate (16). However, the PRCC values of the parameters are equal to or lower in magnitude at the end time of 2030 for the constant death rate than for the variable death rate. This could be due to the fact that the variable death rate results in higher number of deaths, which has the effect of lowering the HUD class. We always want to lower the number of individuals in the HUD class in beneficial ways. However, with the higher death rate it becomes more crucial for individuals to exit out of the HUD class quicker due to the increased risk of heroin-related overdose. If the treatment rates and/or recovery rates could be increased and more users leave the HUD class and enter treatment, they would be protected from those resulting dangers that could lead to a heroin-related overdose death. It is vital at the higher death rates to get individuals out of the HUD class quicker than for the lower death rate.

Table 3

PRCC results for movement into I (HUD) using baseline parameters (5) and either constant or the variable in (16). The initial conditions for in 2020 were generated using (4), (5), 2002 values of , , , , and running the system until 2020 (as previously described to obtain Fig. 2). The PRCC values at 2030 are given here with the columns labeled “constant” corresponding to the constant death rate of (its extrapolated 2020 value) and the columns “variable” corresponding to the variable death rate defined in (16)

IC/ param	Yearly new I		Yearly relapse from T		Yearly relapse from R		Yearly deaths
	Constant	Variable	Constant	Variable	Constant	Variable	Constant	Variable
S(0)	–	–	–	–	–	–	–	–
I(0)	0.93	0.97	0.85	0.89	0.94	0.95	0.95	0.96
T(0)	0.72	0.81	0.50	0.62	0.76	0.78	0.69	0.79
R(0)	0.58	0.71	–	0.50	0.78	0.79	0.62	0.68
Λ∗	–	–	–	–	–	–	–	–
μ	–	- 0.49	–	–	- 0.45	- 0.54	- 0.40	- 0.52
β	0.99	0.99	0.79	0.85	0.85	0.88	0.93	0.95
η_1	- 0.58	- 0.72	0.78	0.81	–	–	- 0.65	- 0.75
η_2∗	–	–	–	–	–	–	–	–
η_3	–	–	0.39	0.40	–	–	- 0.29	–
ρ	–	–	- 0.47	- 0.56	0.89	0.90	–	–
κ	0.55	0.71	0.94	0.96	–	- 0.41	0.60	0.72
α_1	0.65	0.78	–	0.54	0.89	0.90	0.67	0.76
α_2∗	–	–	–	–	–	–	–	–
δ	- 0.63	–	- 0.45	–	- 0.53	–	0.96	–
m	–	- 0.42	–	–	–	–	–	0.85
b	–	- 0.74	–	- 0.52	–	- 0.49	–	0.87
ω	- 0.45	- 0.50	–	–	0.90	0.92	- 0.51	- 0.58
ϵ∗	0.58	0.69	- 0.75	- 0.82	–	–	0.65	0.64

All table entries without a value are not significant. The notation of * denotes that the parameter does not appear in the formula for . The corresponding graphs for this table are given in Figs. 6, 7, 8, 9

Fig. 6

PRCC results over time for those who are entering I (HUD) for the first time, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 3. Top: constant death rate of , its extrapolated 2020 value. Bottom: variable death rate defined in (16) (Color figure online)

Fig. 7

PRCC results over time for those who are entering I (HUD) by relapsing from T, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 3. Top: constant death rate of , its extrapolated 2020 value. Bottom: variable death rate defined in (16) (Color figure online)

Fig. 8

PRCC results over time for those who are entering HUD by relapsing from R, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 3. Top: constant death rate of , its extrapolated 2020 value. Bottom: variable death rate defined in (16) (Color figure online)

Fig. 9

PRCC results over time for the yearly HUD overdose deaths, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 3. Top: constant death rate of , its extrapolated 2020 value. Bottom: variable death rate defined in (16) (Color figure online)

The yearly new I variable keeps track of the number of individuals from the S class who are entering the I (HUD) class; see Fig. 6 and Table 3. The comparisons of the PRCC values for the yearly deaths due to overdose at the year-end time of 2030 were very similar for both death rates. What follows discusses both death rates unless otherwise noted. The parameter with the highest significance (ranked highly significant for both death rates) to focus on would be (the transmission rate of becoming an individual with HUD through interaction with others in the HUD class). Since this parameter is positively correlated, a decrease of the transmission rate would result in a decrease of the yearly new to HUD counts, as expected. Although not as significant as the transmission rate, but ranked somewhat significant to significant, other parameters to consider for focus would be (the rate of individuals in the recovered state relapsing back to the HUD class on their own), (saturation term for entering treatment), (the rate of individuals leaving treatment and returning to the HUD class), (the rate of individuals in HUD who enter a specialty treatment facility on their own), and (death rate of individuals in the HUD class due to overdose). Thus, decreasing the relapse rate from treatment and the recovered class, increasing availability of specialty treatment, or increasing the rate of access for someone to enter treatment on their own would all decrease the yearly new to HUD counts. Although an increase in the HUD death rate would decrease the counts, as expected since less individuals remaining in HUD would result in less of the S class moving to HUD, ethically, we do not want the counts to decrease because of less interactions due to the high death rate. Therefore, the previously mentioned parameters, other than the HUD death rate, are best to focus on.

The Yearly relapse T variable keeps track of the number of individuals who were in treatment and have relapsed back into the HUD class; see Fig. 7 and Table 3. The comparisons of the PRCC values for the yearly deaths due to overdose at the year-end time of 2030 were very similar for both death rates. What follows discusses both death rates unless otherwise noted. The parameter with the highest significance (ranked highly significant) to focus on would be (the rate of individuals leaving treatment and returning to the HUD class). Since is positively correlated, a decrease in the relapse rate of treatment decreases the yearly relapse T counts, as expected. Other parameters that followed in significance (ranked somewhat significant to significant) are (the transmission rate of becoming an individual with HUD through interaction with others in the HUD class), (the rate of individuals in HUD who enter a specialty treatment facility on their own), and (saturation term for entering treatment). Thus, decreasing the transmission rate, increasing the rate of individuals entering treatment by one’s own accord, and increasing availability of treatment would all decrease the yearly relapse T counts.

The yearly relapse R variable keeps track of the individuals who were in the R class and have relapsed back into the I (HUD) class whether on their own or by being in contact with someone from the HUD class; see Fig. 8 and Table 3. The comparisons of the PRCC values for the yearly deaths due to overdose at the year-end time of 2030 were very similar for both death rates. What follows discusses both death rates unless otherwise noted. The parameter with the highest significance (ranked highly significant) is (the rate of individuals in HUD who enter the recovered class by either treatment in non-specialty facilities and/or “quitting cold turkey”). Since is positively correlated, a decrease in the number of people entering the recovered class decreases the number of yearly relapse R counts; however, we do not want the count to decrease by lowering the rate individuals go into recovery. Hence, we look at the next most significant parameters (ranked highly significant) which are (the rate of individuals leaving treatment and entering a “recovered” state), (the rate of individuals in the recovered state relapsing back to the A class on their own), and (the transmission rate of becoming an individual with HUD through interaction with others in the HUD class). Similar to the analysis for , we ignore decreasing to decrease the counts, and thus, the most significant parameters to focus on would be and . Hence, decreasing the transmission rate and/or decreasing the relapse rate of individuals from R to HUD would decrease the yearly relapse R counts.

Yearly deaths:

The yearly deaths variable accounts for the number of yearly deaths due to overdose by HUD individuals; see Fig. 9 and Table 3. The comparisons of the PRCC values for the yearly deaths due to overdose at the year-end time of 2030 were very similar for both death rates. What follows discusses both death rates unless otherwise noted. The two most significant parameters (ranked highly significant) are the death rates, (death rate of individuals in the HUD class due to overdose)(scenario (i)), b and m (scenario (ii)), and (the transmission rate of becoming an individual with HUD through interaction with others in the HUD class). Hence, lowering the HUD death rate and transmission rate would decrease the yearly death counts, as expected. Other parameters (ranked somewhat significant to significant) are (the rate of individuals in the recovered state relapsing back to the HUD class on their own), (saturation term for entering treatment), (the rate of individuals leaving treatment and returning to the HUD class), (the rate of individuals in HUD who enter a specialty treatment facility on their own), and (the rate of individuals in HUD who enter the recovered class by either treatment in non-specialty facilities and/or “quitting cold turkey”). Therefore, lowering the relapse rates from treatment and the recovered class, increasing availability for treatment, increasing the rate of number of individuals entering treatment on their own, and/or increasing the rate of individuals entering the recovered class would all result in decreasing the yearly deaths.

Conclusion

Our paper presents a deterministic model for the dynamics of an illicit opioid use disorder (IOUD) model. Besides a traditional susceptible class and a class of individuals with illicit OUD, our model includes a treatment class for individuals in specialty treatment facilities. It further includes a recovered population class that holds individuals who have either completed treatment (specialty or non-specialty) or “quit cold turkey.” Here, they may remain or relapse back to the IOUD class. Our model also includes a saturation treatment function, which slows down the rate of entry into treatment due to the lack of availability of specialty treatment. Realistic parameter estimates are obtained from the literature and via parameter estimation to match the available SAMHSA data from 2002–2019. The overdose death rate for those in the IOUD class is seen to have been increasing at a linear rate since around 2011. In addition, since our model approaches a constant population , scaling the SAMHSA data to a population of 200,000 allows us to better see the dynamics of this heroin epidemic.

For the parameter estimates we found, the -value extrapolated to 2030 results in a situation where the effective reproductive number, , is less than one, yet a region of bi-stability exists in the –– space in which both EE and DFE are stable. There is a backward bifurcation that occurs just below as is varied (for fixed ) and just below as is varied (for fixed ) illustrating an additional difficulty of eradicating HUD. This region of bi-stability predicts a minimum -value below which we will not have bi-stability. Thus, ensuring adequate access to specialty treatment facilities is important. In addition, while our model has a backward bifurcation for no saturation, it requires an unrealistically large nonlinear relapse rate ; in contrast, with saturation, a backward bifurcation exists above the minimum -value for a realistic nonlinear relapse rate (including a value of ).

A surprising discovery in our analysis was that if the growth of the illicit OUD overdose death rate continues on its path of the last 10 years, by 2038 the DFE will be the only stable biologically relevant equilibrium. While we do want this epidemic to end, we do not want it to end because of overdose deaths from illicit opioid use. Law enforcement intervention, policies, and/or strategies can be taken to either slow the increase of , keep the rate constant, or possibly reduce it.

While many of the results of our sensitivity analysis were expected, one result stood out—the consistent importance of , which is the parameter quantifying the rate at which someone in HUD enters T on their own accord. Out of the three variables to move into treatment, was more important than , entering treatment because of interaction with a susceptible, or , entering treatment because of an interaction with a recovered person. It would seem beneficial in the short term to increase efforts for ways that make it easier for an individual to enter treatment if needed. This could be through things such as financial support for treatment or perhaps lowering the stigma to increase willingness to seek out help on their own as well.

Future work could include extensions to the model such as incorporating a prescription class, a “casual user” class, or a second treatment class for non-specialty. Finally, parameter estimation revealed the necessity of additional statistics that could be calculated and questions that could be asked by SAMHSA in their NSDUH that would allow for a better comparison of model output to data, including calculating heroin use disorder (HUD) in the last month and determining synthetic opioids use with all the time frames given including “in the last month.” Keeping track of whether those individuals in treatment came from the I class or from a “casual user” class would also help in estimating parameters.

Final notes: During the revisions of this paper, the SAMHSA data for 2020 were released. We observe that there was a change in the definition of individuals with substance use disorder (SUD), including HUD, due to the switch in criteria for classifying these individuals. “Beginning with the 2020 NSDUH, SUD estimates for alcohol and illicit drugs were based on criteria in the Diagnostic and Statistical Manual of Mental Disorders, 5th edition,” where previously the 4th edition was used (Center for Behavioral Health Statistics and Quality 2021). Due to the different definition for classifying HUD, we cannot directly incorporate the new data into our model and leave it to future work to determine how to incorporate it.

Table 4

PRCC results for those that completed treatment (yearly completed treatment), those that went to treatment (yearly treatment), and those who are in those who left the HUD class either quitting on their own or with the help of a non-specialty treatment facility (yearly I to R), using baseline parameters (5) and either constant or the variable in (16). The initial conditions for in 2020 were generated using (4)–(5), 2002 values of , , , , and running the system until 2020 (as previously described to obtain Fig. 2). The PRCC values at 2030 are given here with the columns labeled “constant” corresponding to the constant death rate of (its extrapolated 2020 value) and the columns “variable” corresponding to the variable death rate defined in (16)

IC/param	Yearly completed treatment		Yearly treatment		Yearly I to R
IC/ param	Constant	Variable	Constant	Variable	Constant	Variable
S(0)	–	–	–	–	–	–
I(0)	0.91	0.87	0.89	0.90	0.92	0.94
T(0)	0.65	0.61	0.58	0.65	0.68	0.73
R(0)	0.51	0.41	0.42	0.51	0.55	0.59
Λ∗	–	–	–	–	–	–
μ	- 0.40	–	–	–	–	- 0.48
β	0.88	0.85	0.89	0.91	0.90	0.93
η_1	0.84	0.78	0.82	0.82	- 0.52	- 0.63
η_2∗	–	–	–	–	–	–
η_3	0.44	0.48	0.45	0.48	–	-
ρ	0.98	0.97	–	–	–	–
κ	- 0.84	- 0.85	0.90	0.92	0.57	0.63
α_1	0.59	0.52	0.52	0.59	0.58	0.66
α_2∗	–	–	–	–	–	–
δ	- 0.49	–	- 0.55	–	- 0.63	–
m	–	–	–	–	–	–
b	–	- 0.40	–	- 0.59	–	- 0.65
ω	–	–	–	–	0.94	0.96
ϵ∗	- 0.85	- 0.82	- 0.80	- 0.83	0.54	0.55

All table entries without a value are not significant. The notation of * denotes that the parameter does not appear in the formula for . The corresponding graphs for this table are given in Figs. 10, 11, and 12

Table 5

PRCC results for those that are susceptible (model output S), those that are in the HUD class (model output I), those who are in treatment (model output T), and those that have recovered (model output R), using baseline parameters (5) and either constant or the variable in (16). The initial conditions for in 2020 were generated using (4)–(5), 2002 values of , , , , and running the system until 2020 (as previously described to obtain Fig. 2). The PRCC values at 2030 are given here with the columns labeled “constant” corresponding to the constant death rate of (its extrapolated 2020 value) and the columns “variable” corresponding to the variable death rate defined in (16)

IC/ param	S		I		T		R
	Constant	Variable	Constant	Variable	Constant	Variable	Constant	Variable
S(0)	0.99	0.99	–	–	–	–	–	–
I(0)	–	–	0.94	0.98	0.93	0.89	0.91	0.94
T(0)	–	–	0.72	0.85	0.72	0.63	0.63	0.73
R(0)	–	–	0.58	0.76	0.55	–	0.68	0.74
Λ∗	0.93	0.91	–	–	–	–	–	–
μ	- 0.93	- 0.91	- 0.40	- 0.63	- 0.41	–	–	- 0.52
β	–	–	0.93	0.97	0.90	0.84	0.81	0.84
η_1	–	–	- 0.56	- 0.80	0.87	0.78	–	–
η_2∗	–	–	–	–	–	–	–	–
η_3	–	–	- 0.23	- 0.35	0.49	0.44	–	–
ρ	–	–	–	- 0.42	- 0.66	- 0.55	0.81	0.86
κ	–	–	0.58	0.78	- 0.87	- 0.81	–	–
α_1	–	–	0.66	0.82	0.58	0.48	- 0.94	- 0.96
α_2∗	–	–	–	–	–	–	–	–
δ	–	–	- 0.66	–	- 0.56	–	–	–
m	–	–	–	0.49	–	–	–	–
b	–	–	–	- 0.81	–	- 0.50	–	- 0.44
ω	–	–	- 0.46	- 0.61	–	–	0.83	0.89
ϵ∗	–	–	0.61	0.75	- 0.87	- 0.79	–	–

All table entries without a value are not significant. The notation of * denotes that the parameter does not appear in the formula for . The corresponding graphs for this table are given in Figs. 14, 13, and 15