Dynamics of the tumor-immune-virus interactions: Convergence conditions to tumor-only or tumor-free equilibrium points

Table 1.  Table of acronyms.

Acronym	Its Expansion
APA	American Psychological Association
CDC	Centers for Disease Control and Prevention
DFE	Disorder-Free Equilibrium
DSM-IV	Diagnostic and Statistical Manual of Mental Disorders, Fourth Edition
EE	Endemic Equilibrium
HUD	Heroin Use Disorder
IOUD	Illicit Opioid Use Disorder
LHS	Latin Hypercube Sampling
NIH	National Institutes of Health
NSDUH	National Survey on Drug Use and Health
OUD	Opioid Use Disorder
PRCC	Partial Rank Correlation Coefficient
SAMSHA	Substance Abuse and Mental Health Services Administration
SUD	Substance Use Disorder

---

 | Show Table

DownLoad: CSV

1.   Introduction

Opioid use disorder (OUD) is a monumental health concern in the United States and considered to be an epidemic, claiming over 100,000 lives due to overdose deaths from April 2020 to April 2021, according to the Centers for Disease Control and Prevention (CDC) ^[1]. Data collected by the Substance Abuse and Mental Health Services Administration (SAMSHA) showed that in 2020, opioids were misused by 9.5 million people (Substance Abuse and Mental Health Services Administration (2021)). Yet OUD is a constantly changing phenomenon that is in great need of not only improved treatment approaches, but also more refined tracking and prediction of opioid usage trends across distinct populations.

As a point of emphasis, several "waves" of the opioid surge in opioid abuse and opioid-related deaths have been observed over the past 30 years. This was first noted as an overprescription and overmarketing of opioid pain relievers by medical professionals and pharmaceutical industries. This led to implementation of a number of restrictions on opioid prescriptions, resulting in opioid-dependent individuals resorting to the use of less expensive and more easily accessible heroin. In more recent years, drug markets have been infiltrated by illicitly manufactured fentanyl, an extremely potent and potentially fatal opioid, into the drug supply. Currently, the opioid epidemic has been confounded by an increase in the concurrent use of psychostimulants such as cocaine and methamphetamine to counteract the highly sedating effects of opioids, as well as co-abuse of alcohol and benzodiazepines such as alprazolam (Xanax) to lessen the symptoms of opioid withdrawal ^[2].

1.1.   Brain science and addiction

The term "opioids" refers to all chemicals that act on opioid receptors, which are proteins in the brain and body that are selectively activated by these molecules. Various subcategories of opioids include naturally occurring morphine, codeine and thebaine that is found in the opium poppy plant Papaver somniferum, semisynthetic opioids such as heroin, and fully synthetic opioids such as fentanyl and oxycodone. In addition to these exogenous opioids, throughout the body is a complex endogenous opioid system that utilizes natural peptide messengers such as endorphins to activate opioid receptors. The primary subtypes of opioid receptors are mu, delta and kappa opioid receptors (MOR, DOR and KOR, respectively; also referred to as $\mu$, $\delta$, and $\kappa$ receptors), as well as opioid receptor like-1 (ORL1), which binds the endogenous opioid peptide nociceptin. Abused opioids act primarily, but not exclusively, at $\mu$ opioid receptors as agonists, activating the receptor to produce specific responses in the cells on which they reside (i.e., neurons, immune cells). Remedies for opioid overdose such as naloxone (Narcan) act by competitively blocking the binding of opioids to their receptors.

Addiction to opioids and other drugs of abuse are considered chronic relapsing conditions, and more recently have been described as diseases of the brain ^[3]. A disease can be defined as a circumstance that impairs the normal operations of a living organism or one of its components, and evidenced by specific characteristics and symptoms. The disease theory of addiction challenges prior notions that addiction is a moral failure or weakness of character, and has helped decrease the stigmatization of addiction and increase accessibility to treatment. However, the disease theory of addiction is not universally accepted, as it has been argued that it downplays individual accountability, the importance of psychosocial and environmental factors, underemphasizes the need for understanding brain recovery from addiction, and biases treatment approaches towards medical approaches over those that are more holistic or psychological in nature ^[4].

Regardless of which side of the disease model debate one is on, there is a wealth of evidence that chronic use of opioids fundamentally alters brain structure and function. For example, interactions and synchronization between specific neural circuits, specifically the prefrontal cortex which governs executive function and the limbic system which mediates emotions, are disrupted by chronic use of prescription opioids ^[5]. Chronic heroin use has been shown to reduce the volume of gray matter in the cerebral cortex and as well as that of the brain's reward circuitry ^[5]. On a more microscopic level, repeated use of morphine and other opioids can modify the fine delicate structures of neurons in various regions of the brain including those that comprise the brain's reward system ^[6]. Medication assisted therapies for opioid dependence, which utilize substitutes for the previously abused opioid, are administered in a supervised manner and are efficacious in reducing relapse, particularly when coupled with psychosocial support. Such medications include the long acting $\mu$ receptor agonist methadone and the partial $\mu$ receptor agonist buprenorphine (formulated as Suboxone). Future treatment approaches that are being considered include non-invasive neuromodulation techniques (i.e., transcranial magnetic stimulation, low intensity focused ultrasound), psychedelic assisted therapy, and digitally based interventions.

1.2.   Mathematical epidemiological models

Mathematical epidemiological models have been used to study and track health issues and patterns in a wide variety of contexts ^[7,8]. In the context of this manuscript, previous work has examined the heroin epidemic (e.g., ^[9,10]) and opioid epidemic (e.g., ^[11]), and recent work has looked at an age-structured model for a drug addiction forecasting method that assimilates observational data of addiction and overdose mortality ^[12,13]. Our work does not venture into age structure to forecast but instead examines additional epidemiological classes in the epidemics where data is available, providing a more traditional epidemiological approach that can still be tested against data. This present work models the spread of heroin as a social contagion and and then extends this framework to examine illicit opioid drug use and abuse. We do not intend to assert that this model applies to the more general opioid epidemic where prescriptions are involved and we do not utilize such data. We instead focus only on illicit opioid use with an emphasis on heroin use. 'Social Contagion' is defined by the American Psychological Association (APA) Dictionary of Psychology as "the spread of behaviors, attitudes, and affect through crowds and other types of social aggregates from one member to another" ^[14]. While behavior of an individual is what classifies individuals in our model, we recognize that there are likely specific external conditions that trigger these conditions together with internal conditions that may perpetuate them. We have found studies whose results suggest that drug abuse can be transmitted in an environmentally mediated manner^[15,16], although exactly how this happens for any drug is not yet known. To help shed additional light on this important question of how, the National Institutes of Health (NIH) has recently created a grant mechanism to undertake investigations in understanding the social contagion of behavior and substance abuse ^[17]. While having an interaction with someone using illicit opioids is neither a necessary nor a sufficient condition for someone to transition from the susceptible class to using illicit opioids (since illicit opioid use can be self-induced, although the question then arises about how the illicit opioids were obtained), the spread of illicit drug use disorder across the population, the waves of opioid deaths, and recent studies suggesting drug abuse can be transmitted in an environmentally mediated manner all combine to support such an approach ^[15,16]. In addition, authors have studied mathematical models of drug use before from the perspective of the spread as a social contagion. A summary of those studies may be found in Cole and Wirkus ^[10].

This paper presents a mathematical model of nonlinear ordinary differential equations (ODE) to understand better the complex issues surrounding illicit OUD, its treatment options, and methods for decreasing relapse. By describing the spread of IOUD as a potential contagion, we assert that the IOUD treatment-relapse cycle is modeled within the context of disease epidemiology. Moreover, the dynamics underlying those patterns can best inform control.

2.   Materials and methods

2.1.   Mathematical model

The population is divided into 6 classes: susceptible ($S(t)$), exposed ($E(t)$), treatment for the exposed ($T_E(t)$), IOUD ($I(t)$), specialty treatment for IOUD ($T(t))$ (as defined by SAMHSA), and recovered ($R(t)$) at time $t > 2002$. $N = S+E+T_E+I+T+R$ denotes the total population. The constant influx, $\Lambda$, into the susceptible population together with the remaining equations result in the overall population approaching a constant level. We refer to our model as the IOUD model with a casual user class. While our earlier model, found in Cole and Wirkus (2022) ^[10], focused on heroin use disorder (HUD) this present model extends the HUD model to an IOUD model due to availability of SAMHSA data and the similarities with HUD when only illicit use (and not prescription use) is considered. Additionally, we introduce two new compartments. First, the exposed class, $E(t)$, holds the number of individuals who are using illicit opioids but do not have IOUD as defined by the Diagnostic and Statistical Manual of Mental Disorders-Fourth Edition (DSM-IV) ^[18].^* The treatment class for the exposed (i.e., casual users), $T_E(t)$, holds the number of individuals who are considered exposed but are in specialty treatment as defined in SAMHSA ^[19]. The IOUD class, $I(t)$, holds the number of individuals who are using illicit opioids and have IOUD as defined by the DSM-IV ^[18]. The treatment class for IOUD, $T(t)$, holds the number of individuals who have IOUD. The recovered class, $R(t)$, contains the number of individuals who had IOUD and either completed treatment, quit "cold turkey", or completed non-specialty treatment as defined in SAMHSA ^[19]. Finally, the susceptible class, $S(t)$, holds the number of individuals susceptible to opioid use or misuse. These susceptible individuals may have passed through $E$ or $T_E$ but not $I$, $T$, or $R$. We show the compartments and the interactions among them in Figure 1.

^*While DSM-V was released in 2013 and provides updated guidance on various classifications and definitions, SAMHSA did not use it until their gathering and release of the 2020 data. This current manuscript concerns data up to 2019. Additionally, the phrase “dependence or abuse” was used up until 2014 in the SAMHSA data whereas the phrase “use disorder”, defined in DSM-V, was used beginning in 2015 with the note “Substance Use Disorder is defined as meeting criteria for illicit drug or alcohol dependence or abuse.”

Figure 1.  Flow diagram of the IOUD Model with a casual user class: arrows show the progression of the change of classes. $S$ represents the susceptible individuals, $E$ represents the exposed individuals, $T_E$ represents those exposed in specialty treatment, $I$ represents individuals with IOUD, $T$ represents those in specialty treatment for IOUD, and $R$ represents recovered users. Due to a greater potential for relapse, $R$ is considered distinct from $S$. Due to limited access to care, $b(T, T_E) = \frac{1}{1+\epsilon T+ \epsilon_E T_E}$ represents the reduced rate of entry into the $T_E$ and $T$ class.

DownLoad: Full-Size Img PowerPoint

The recruitment rate into the susceptible population is denoted by $\Lambda$. The natural death rate for all classes is $\mu$. The transmission rate from susceptibles to exposed due to an interaction with someone who has IOUD is denoted by $\beta$. The transmission rate from susceptibles to exposed due to an interaction with someone using illicit opioids but not considered having IOUD is denoted by $\beta_E$. There are three avenues for someone from the exposed class to enter a specialty treatment facility. First, the rate of someone from the exposed class entering specialty treatment on their own accord is denoted by $\psi_1$. Second, the rate of someone from the exposed class entering a specialty treatment facility due to an interaction with someone from the recovered class is denoted by $\psi_2$. Third, the rate of someone from the exposed class entering a specialty treatment facility due to an interaction with someone from the susceptible class is denoted by $\psi_3$. An exposed individual could also stop using illicit opioids on their own or in some non-specialty treatment facility and return to the susceptible class at a rate denoted by $\zeta$. A previously exposed individual may complete specialty treatment and cycle back to the susceptibles by rate $\rho_E$; however, they may relapse from specialty treatment to using illicit opioids by rate $\kappa_E$. Lastly, an exposed individual could develop IOUD by rate $\chi$ and transfer to the IOUD class.

We now consider the remaining three classes identical to those considered in Cole and Wirkus (2022) ^[10]. Once an individual is in $I$ (the IOUD class), there are three avenues that they could enter into a specialty treatment facility. First, the rate of someone from the IOUD class entering specialty treatment on their own accord is denoted by $\eta_1$. Second, the rate of someone from the IOUD class entering a specialty treatment facility due to an interaction with someone from the recovered class is denoted by $\eta_2$. Third, the rate of someone from the IOUD class entering a specialty treatment facility due to an interaction with someone from the susceptible class is denoted by $\eta_3$. An IOUD individual could also stop using illicit opioids on their own or through a non-specialty treatment facility and transfer to the recovered class at a rate denoted by $\omega$. An IOUD individual may complete specialty treatment and flow into the recovered class by rate $\rho$. However, they may relapse from specialty treatment to using illicit opioids at a rate $\kappa$. Finally, individuals in the recovered class may cycle back to the IOUD class. They either may relapse on their own accord at rate $\alpha_1$ or by the influence of someone in the IOUD class by rate $\alpha_2$.

There is an added death rate component due to overdose deaths. This added rate is $\delta_E$ for the exposed class, and for the IOUD class, the added rate is $\delta$. The parameters $\delta$ and $\delta_E$ will be piecewise functions of time and not differentiable at the breakpoint. As will be discussed in the parameter estimation sections containing (2.3)–(2.4) and (2.6)–(2.7), function $\delta$ is continuous but just not-differentiable at point $c$, whereas $\delta_E$ has a jump discontinuity at point $c$.

Due to the limited access to care, a saturation treatment function limits the flow into the specialty treatment facilities. Because individuals in both the IOUD and $E$ classes can enter treatment, the saturation term is a function of both: $b(T, T_E)$. Therefore, the saturation parameter corresponding to the casual users is $\epsilon_E$, whereas the saturation term correlating to those in the IOUD class is $\epsilon$.

We thus have the following deterministic system of nonlinear ordinary differential equations based on the previous assumptions, casual use of illicit opioids, illicit opioid use, treatment, and recovery:

where $b(T, T_E) = \dfrac{1}{1+\epsilon T + \epsilon_E T_E}$ and all parameters are nonnegative. Table 2 gives a description of the variables. Table 3 gives a description of the parameters.

Table 2.  Description of variables of the IOUD model with a casual user class.

	Description
$S(t)$	Susceptible population
$E(t)$	Casual users
$T_E(t)$	Casual users now in specialty treatment
$I(t)$	Users that have IOUD
$T(t)$	Users from the $I$ class now in specialty treatment
$R(t)$	The recovered class

---

 | Show Table

DownLoad: CSV

Table 3.  Description of parameters of the IOUD model with a casual user class.

Parameter	Description	Units
$\Lambda$	Recruitment into the susceptible population	people/year
$\mu$	Natural death rate	1/year
$\beta$	Transmission rate from susceptible to exposed through interaction with someone from the $I$ class	1/year
$\beta_E$	Transmission rate from susceptible to exposed through interaction with someone from the $E$ class	1/year
$\zeta$	The rate of individuals in the $E$ class returning to the $S$ class	1/year
$\chi$	The rate of individuals in the $E$ class that transition to the $I$ class	1/year
$\psi_1$	The rate of individuals in $E$ who enter specialty treatment on their own	1/year
$\psi_2$	The rate of individuals in $E$ who enter specialty treatment through interaction with a recovered individual	1/year
$\psi_3$	The rate of individuals in $E$ who enter specialty treatment through interaction with a susceptible individual	1/year
$\rho_E$	The rate of casual users leaving treatment and entering the $S$ class	1/year
$\kappa_E$	The rate of casual users leaving treatment and returning to the $E$ class	1/year
$\eta_1$	The rate of individuals in $I$ who enter specialty treatment on their own	1/year
$\eta_2$	The rate of individuals in $I$ who enter specialty treatment through interaction with a recovered individual	1/year
$\eta_3$	The rate of individuals in $I$ who enter specialty treatment through interaction with a susceptible individual	1/year
$\omega$	The rate of individuals in $I$ who enter the recovered class by either completing treatment in non-specialty facilities or "quitting cold turkey"	1/year
$\rho$	The rate of individuals leaving treatment and entering the recovered class	1/year
$\kappa$	The rate of individuals leaving treatment and returning to the $I$ class	1/year
$\alpha_1$	The rate of individuals in the recovered state relapsing back to the $I$ class on their own	1/year
$\alpha_2$	The rate of individuals in the recovered state relapsing back to the $I$ class through interaction with an individual in the $I$ class	1/year
$\delta$	Added overdose death rate for the $I$ class	1/year
$\delta_E$	Added overdose death rate for the $E$ class	1/year
$\epsilon$	Saturation term for entering a specialty treatment facility from the $I$ class	1/people
$\epsilon_E$	Saturation term for entering a specialty treatment facility from the $E$ class	1/people

---

 | Show Table

DownLoad: CSV

2.2.   Non-negativity and boundedness

The following will investigate the fundamental dynamical properties of the illicit opioid use disorder (IOUD) model with a casual user class.

Since $N(t) = S(t)+E(t)+T_E(t)+I(t)+T(t)+R(t)$, we have $\frac{dN(t)}{dt} = \frac{dS(t)}{dt}+\frac{dE(t)}{dt}+\frac{dT_E(t)}{dt}+\frac{dI(t)}{dt}+\frac{dT(t)}{dt}+\frac{dR(t)}{dt}$.

We add the five equations of (2.1), and find the total population dynamics are driven by the following equation:

Since the IOUD model with a casual user class tracks physical entities, all associated parameters are nonnegative. We state the following two theorems, noting that the proofs are fairly standard and thus are not included (see, e.g., ^[20,21]).

Theorem 1 Local solutions to the IOUD model with a casual user class with initial data in the region

$S(0) = S_0 > 0, E(0) = E_0 > 0, T_E(0) = T_{E_0} > 0, I(0) = I_0 > 0, T(0) = T_0 > 0, R(0) = R_0 > 0, N(0) = N_0 > 0,$ exist and are unique.

Theorem 2 Given $N > 0$ and nonnegative initial conditions and parameter values, solutions to the IOUD model with a casual user class are nonnegative on the interval of existence.

Therefore the coordinate planes and hence the positive octant $\Omega = \{(S, E, T_E, I, T, R)\in \mathbb{R}_+^6:0 < S, 0 < E, 0 < T_E, 0 < I, 0 < T, 0 < R\},$ are invariant under the local flow.

2.3.   Heroin only - data explanation and parameter estimation

To compare the model to data, we first consider using a heroin dataset (defined by the CDC as overdose death due to heroin-only or heroin mixed with synthetic opioids) that does not include non-heroin use. This section compares our model to data considering only the use of heroin or heroin mixed with synthetic opioids (e.g., fentanyl). Finally, we present the data used for the IOUD model with the casual user class for a heroin-only use dataset in Table 4.

Table 4.  Data for U.S., 2002–2020. The number of overdose deaths for 2002–2020 are from the CDC ^[22]. U.S. population comes from ^[23]. Use disorder and specialty treatment data come from SAMHSA's National Survey on Drug Use and Health (NSDUH) ^{[19,24,25,26,27,28,29,30,31,32,33]}. Data with an asterisk $*$ = Specialty treatment data $\times$ 0.6874 because specialty treatment from $I$ only asked in 2014–2017 SAMHSA surveys. The factor 0.6874 is the average of the ratio of specialty treatment from $I$ to specialty treatment in the 4 years when data is available.

		US		specialty	specialty
	heroin	population	HUD	treatment	treatment	initiation	use = E + I
	deaths	(in millions)	in last yr.	in last yr.	from I	in last yr.	in last yr.
2002	2089	287.3	214,000	NA	NA	117,000	404,000
2003	2080	289.8	189,000	NA	NA	92,000	314,000
2004	1878	292.4	270,000	156,000	107,200$^*$	118,000	398,000
2005	2009	295.0	227,000	190,000	130,600$^*$	108,000	379,000
2006	2088	297.8	324,000	377,000	259,100$^*$	90,000	560,000
2007	2399	300.6	214,000	201,000	138,200$^*$	106,000	373,000
2008	3041	303.5	283,000	227,000	156,000$^*$	116,000	455,000
2009	3278	306.3	369,000	322,000	221,300$^*$	187,000	582,000
2010	3036	309.0	361,000	274,000	188,300$^*$	142,000	621,000
2011	4397	311.6	426,000	292,000	200,700$^*$	178,000	620,000
2012	5925	314.0	467,000	293,000	201,400$^*$	156,000	669,000
2013	8257	316.4	517,000	359,000	246,800$^*$	169,000	681,000
2014	10,574	318.7	586,000	428,000	270,000	212,000	914,000
2015	12,989	320.9	591,000	398,000	242,000	135,000	828,000
2016	15,469	323.0	626,000	365,000	235,000	170,000	948,000
2017	15,482	325.1	652,000	413,000	358,000	81,000	886,000
2018	14,996	327.1	526,000	424,000	291,500$^*$	117,000	808,000
2019	14,019	329.1	438,000	467,000	321,000$^*$	50,000	745,000
2020	13,058	331.0	NA	NA	NA	103,000	NA

---

 | Show Table

DownLoad: CSV

Column 2 of Table 4 displays the yearly number of overdose deaths due to heroin as found by the CDC ^[22]. Column 4 gives the number of individuals who reported having heroin use disorder (HUD) within the past year as provided by the NSDUH. (See Table 4 for those references.) This data relates to the state variable $I$. Column 5 shows the number of individuals who reported to SAMHSA that they were in a specialty treatment facility due to heroin, regardless of whether they had HUD, within the past year. Column 6 gives a count of those individuals who reported to SAMHSA that they were in a specialty treatment facility due to HUD within the past year; this data was provided for 2014–2017. Hence, the data modified by the asterisk is the data in the previous column for 2004–2013 and 2018–2019, scaled by a factor of 0.6874. (This factor was found by averaging the data from the specialty treatment due to HUD divided by the specialty treatment due to heroin, regardless of HUD for the given data found.) This column relates to the state variable $T$. Subtracting column 6 from column 5 gives us data related to the state variable $T_E$. Column 7 gives us the yearly count of those individuals who reported to SAMHSA that they had initiated heroin use for the first time within the past year. Column 8 shows a count of those individuals who said to SAMHSA that they used heroin within the past year, regardless if they had HUD. Column 4, subtracted from this column, gives us data related to our state variable $E$. Lastly, column 9 shows a count of the number of individuals who reported to SAMHSA that they used heroin within the past month, regardless of whether they had HUD.

The state variables $E$, $T_E$, $I$, and $T$ are instantaneous in time, whereas the SAMHSA data is not. SAMHSA gives a count over the year of those respective classes. Therefore, we correct comparing the data to the variables. The details on how we do this may be found in the Supplementary Information. It is a detailed explanation of how we added a correction to the $I$ variable and the $T$ variable to approximate SAMHSA's yearly count. Similarly, we use this method to find the corrections for the $S$, $E$, and $T_E$ variables.

Using the U.S. population (Table 4 column 3), we scale the national data to a city population of 200,000 individuals for the simulations: for the number of individuals in the HUD class, the number of individuals in specialty treatment from HUD, the number of individuals in specialty treatment from the casual user class, and the number of individuals in the casual user class. For example, in 2002, the HUD data would be calculated as $(214,000/287.3 \mbox{ million})\times200,000 = 148.97$. We depict this value in the top middle graph of Figure 2. This rescaling provides for our analysis a nearly constant population to keep the focus on the dynamics of the problem. This is the data presented in our graphs of Figure 2 with the raw data given in Table 4. The error bars in the graphs represent the standard error given in the SAMHSA data (not presented in the table).

Figure 2.  Fitting model output to scaled data (with error bars when given) for heroin. The red squares are the SAMHSA data and the blue curves are the model output as described in the text. (Top Left): Heroin overdose deaths by year. (Top Middle) Individuals with heroin use disorder (HUD) in the last year. (Top Right) Those who entered specialty treatment from HUD in the past year. (Middle Left) Individuals who used heroin in the last year. (Middle Middle) Those who entered specialty treatment from the casual user class, $E$. (Middle Right) Those who reported using heroin for the first time (initiation from $S$ to $E$). (Bottom Left) The effective reproductive number, $\mathcal{R}_{\mbox{eff}}(t) = (\mathcal{R}_0\cdot S(t)/N(t))$. (Bottom Middle) Overdose death rate for HUD class: asterisks and X-marks are calculated from data (see text and Eq (2.6)). Both lines are calculated witqh a least squares fit. (Bottom Right) Overdose death rate for $E$ Class: asterisks and X-marks are calculated from data (see text and Eq (2.7)).

DownLoad: Full-Size Img PowerPoint

For the data fitting, we have taken some parameter values for our model from the literature and performed a parameter estimation for the remaining parameters using MATLAB and its ${\tt fmincon}$ function. While few parameters in reality remain constant over time, our analysis of the data together with previous published work on similar models indicates that most of the parameters can be considered constant as reasonable first approximations^[9,10,11]. Additional modeling could be done where some of the parameters could be considered time-dependent; however, with the exception of $\delta$ and $\delta_E$, this is beyond the scope of this current work. First, we estimated parameter values and ranges for the IOUD with a casual user class from the literature noting that the time units are per year in all the rates. We used the natural yearly death rate for $\mu = 1/80$ ^[9]. Next, we estimated our yearly recruitment rate, $\Lambda = 2500$, for a population of 200,000, given the natural death rate, ignoring the additional deaths due to overdose (i.e., with $\delta = 0$). We used the approximated ranges from ^[11,34] for the yearly completed treatment rate, $\rho$, from a specialty treatment facility to the recovered class as 0.25–0.6. We assumed the yearly completed treatment rate, $\rho_E$, from a specialty treatment facility to the $S$ class would be higher than $\rho$ because casual users do not have opioid use disorder and may have a quicker recovery time. Hence, we chose a range of 0.5–2. We approximated from Weiss and Rao (2017), Bailey et al. (2013), and Smyth et al. (2010) ^[36,37] our range for the yearly relapse rate, $\kappa$, from the specialty treatment class back to the IOUD class as 0.18–4.0. We assumed that the yearly relapse rate, $\kappa_E$ was 1–2 to fit in this range. We determine those that go to specialty treatment from the IOUD class by estimating $\eta$, the overall yearly rate, to be 0.1–2 ^[9,11], where we used the equation $\eta = \eta_1+\eta_2(R/N)+\eta_3(S/N)$. We set $\eta_2 = 0.7$ and chose the range of 0.8–1.1 for $\eta_1$ and 0.2–5 for $\eta_3$. For the yearly relapse rates from the $R$ class back to the $I$ class we used a study by Gossop et al. (1989), who estimated a range for $\alpha$ of 0.1–1/3 ^[38,39]. However, additional research suggests the relapse rate is significantly higher due to the changes in the brain. Thus, we use the range of 0.1–1 for our $\alpha_i$ ^[40,41]. We use a field of 0.05 to 0.3 for the parameter $\omega$ ^[9]. The ranges for $\beta$ and $\beta_E$, our yearly transmission rates, and $\epsilon$ and $\epsilon_E$, our yearly saturation treatment parameters, were determined from Cole and Wirkus (2022) ^[10] and then determined via parameters estimation. Utilizing the SAMHSA and CDC data, we obtained the rates of the rest of the parameters via parameter estimation: $\chi$, the yearly rate of individuals who go from being a casual user to an individual who now has OUD, $\zeta$, the yearly rate of casual users going back to susceptibles, and $\psi_1$, $\psi_2$, and $\psi_3$, yearly rates from $E$ to $T_E$.

We now examine $\delta$ and $\delta_E$ for the IOUD model with a casual user class using the heroin-only dataset. Battista et al. ^[11] first gave the definition of the death rate $\delta$ due to drug overdose; we refer the reader to Cole and Wirkus (2022) ^[10] for the modification of the definition's derivation for these parameters, and what we present here is in their final forms. We highlight the discussion in ^[10] whereby we concluded that $\delta$ could not be approximated as a constant as the bottom middle subfigure of Figure 2 illustrates. The data from SAMHSA is not perfectly geared for our model because the HUD and OUD are presented as "in the past year" yet the mathematical model gives variables continuous in time. For example, the mathematical model allows for an individual to leave HUD to go into treatment but then relapse back to HUD in the same year. This individual would not be able to overdose from heroin while in treatment but would be able to die from an overdose while using heroin. To address this, we assumed a simple inflow-outflow subdiagram (one for susceptible-casual—casual treatment, and another for HUD—treatment—recovered) presented in the Supplementary Information, that uses the SAMHSA data to estimate what percentage of the given class is actually in that class at any given time. The average ratio of the model output state variable $I$ over the model calculation for each year (keeping track of individuals' movement in and out of classes) is 0.88 while for the variable $E$ it is 0.41. In words, the multiplication factor of 0.88 in the denominator of $\delta$ for heroin use says that of the individuals that SAMHSA classifies as having HUD in the past year, 88% are in the HUD class at any given time with the remaining percentage either in treatment, having stopped using, or having died of an overdose. Finally, based on the data fit from parameter estimation, we estimate the number of HUD overdose deaths due to heroin as 0.89 for $\delta$, and we approximate the number of HUD overdose deaths to heroin as 1 - 0.89 for $\delta_E$. Thus, we have the following definitions.

and

Incorporating these piecewise functions for $\delta$ and $\delta_E$ into our parameter estimation, our baseline values are

Our initial conditions, chosen to approximately pair with the scaled data are $S_0 = 199,800, E_0 = 26, T_{E_0} = 16, I_0 = 132, T_0 = 45, R_0 = 96$. The data match is provided in Figure 2.

2.4.   All-illicit opioids dataset - data explanation and parameter estimation

This section compares the IOUD model with a casual user class to data using all-illicit opioids. Data for the all-illicit opioids dataset is given in Tables 5 and 6. Additionally, we fit the IOUD model with casual users to the new dataset. We use data found in SAMHSA for the IOUD model with the casual user class for an all-illicit opioid use dataset (CDC data includes heroin). The following discussion references the data presented in Table 5. Column 2 of Table 5 displays the yearly number of overdose deaths due to all opioids, as found by the CDC ^[42]. Column 3 gives a count of those individuals who reported having substance use disorder (SUD) within the past year due to the use of pain medications as given by the NSDUH. (See Table 5 for those references.) Column 4 shows a count of those individuals who reported having HUD within the past year given by the NSDUH. Column 5 counts those individuals who reported to SAMHSA, given only for 2015–2019, that they had OUD within the past year (whether due to heroin or pain medication). Given the five years of data in column 5, we found a formula to approximate the data for the missing years of column 5 using columns 3 and 4. This formula (Column 3 $\times$ 0.87 + Column 4) fills in the absent years of column 5 (marked with an asterisk), and that is the data used for our state variable $I$. We found the factor 0.87 gave the best fit and, for the known years, gives approximately 2364.1 (vs. known 2375), 2116.1 (vs. known 2144), 2111.9 (vs. known 2110), 1999.8 (vs. known 2028), and 1626.4 (vs. known 1622). Therefore, an interpretation of 0.87 is that 87% of those diagnosed with SUD for pain medication have SUD for opioids; SAMHSA defines all heroin use disorders as being in the OUD class.

Table 5.  Data for U.S., 2003–2019: overdose deaths, use disorder, and use in past year. Numbers in thousands for all data values. The number of overdose deaths for 2003–2020 are from the CDC ^[22]. Pain medication use disorder, heroin use disorder, opioid use disorder, and use in past year data come from SAMHSA's NSDUH ^{[19,24,25,26,27,28,29,30,31,32,33]}. Data with an asterisk (*) are missing data and are calculated from other columns as described in the text. See the text for discussions of these categories.

	all-illicit				use in	use in	use in
	opioids	SUD-		SUD -	past year	past year	past year
year	OD deaths	pain med	HUD	opioid	pain med	heroin	opioid
2003	12.940	943	189	1009.4*	11,671	314	NA
2004	13.756	1388	270	1477.6*	11,256	398	11,401.5*
2005	14.918	1546	227	1572.0*	11,815	379	11,091.2*
2006	17.545	1635	324	1746.5*	12,649	560	11,603.3
2007	18.516	1707	214	1699.1*	12,466	373	12,576.6*
2008	19.582	1716	283	1775.9*	11,885	455	12,215.7*
2009	20.422	1854	369	1982.0*	12,405	582	11,745.8*
2010	21.089	1923	361	2034.0*	12,242	621	12366.8*
2011	22.784	1768	426	1964.2*	11,143	620	12,250.9*
2012	23.166	2056	467	2255.7*	12,489	669	11,205.9*
2013	25.052	1879	517	2151.7*	11,082	681	12,533.6*
2014	28.647	1918	586	2254.7*	10,337	914	10,734.2*
2015	33.091	2038	591	2375	12,462	828	12,693
2016	42.249	1753	626	2144	11,517	948	11,824
2017	47.600	1678	652	2110	11,077	886	11,401
2018	46.802	1694	526	2028	9948	808	10,250
2019	49.860	1366	438	1622	9724	745	10,065

---

 | Show Table

DownLoad: CSV

Table 6.  Data for U.S., 2003–2019: initiation and specialty treatment. Numbers in thousands for all data values. Pain medication initiation, heroin initiation, and specialty treatment data come from SAMHSA's NSDUH ^{[19,24,25,26,27,28,29,30,31,32,33]}. Data with an asterisk (*) are missing data and are calculated from other columns as described in the text. See the text for discussions of these categories.

					specialty	specialty	specialty
			specialty	specialty	treatment	treatment	treatment
	initiation-	initiation-	treatment	treatment	from $I$	from $I$	from $I$
year	pain med	heroin	heroin	pain med	heroin	pain med	opioid
2003	2456	92	NA	199	NA	132	NA
2004	2422	118	156	226	107.2*	152	215.2*
2005	2193	108	190	259	130.6*	NA	NA
2006	2150	90	377	347	259.1*	238	428.1*
2007	2147	106	201	299	138.2*	195	276.6*
2008	2176	116	227	350	156.0*	201	298.7*
2009	2179	187	322	466	221.3*	320	448.5*
2010	2013	142	274	408	188.3*	271	380.8*
2011	1888	178	292	438	200.7*	335	438.6*
2012	1880	156	293	514	201.4*	427	504.6*
2013	1539	169	359	421	246.8*	345	491.7*
2014	1425	212	428	475	270	348	517.1*
2015	2126	135	398	470	242	371	505.4*
2016	2139	170	365	374	235	307	453
2017	2010	81	413	481	358	346	603
2018	1908	117	424	415	NA	NA	NA
2019	1607	50	467	425	NA	NA	NA

---

 | Show Table

DownLoad: CSV

Column 6 gives the counts by SAMHSA of those individuals who disclosed pain medication use within the past year. Column 7 shows the counts of those individuals who reported heroin use by SAMHSA within the past year. Column 8 counts those individuals who reported to SAMHSA, given for 2015–2019, illicit opioid use within the past year (whether heroin use or opioid pain medication use). We found a formula to approximate the data for the missing years of column 8 using columns 6 and 7. This formula (Column 6 $\times$ 0.95 + Column 7) fills in the absent years of column 8 (and is marked with an asterisk). Comparing this formula with the known data years gives 12,666.9 (vs. known 12,693), 11,889.2 (vs. known 11,824), 11,409.2 (vs. known 11,401), 10,258.6 (vs. known 10,250), and 9982.8 (vs. known 10,065). An interpretation of 0.95 is that 95% of those who misused pain medication in the past year specifically misused opioid pain medication; SAMHSA defines all heroin users as opioid misusers. With this column filled in, we subtract column 5 (approximated and existing values) from the values of column 8 (approximated and existing values) to find our state variable $E$.

The following discussion references the data presented in Table 6.

Column 2 gives us the yearly count of those individuals who reported to SAMHSA having initiated illicit pain medication use for the first time within the past year. Column 3 gives us the yearly count of those individuals who disclosed to SAMHSA having started heroin use for the first time within the past year. Unfortunately, SAMHSA provided no data for years of initiation of illicit opioid use. We thus "guess" that the same factor used to determine opioid misuse from known pain medication misuse and heroin use can give an approximation for initiation to illicit opioid use. Thus, we use the formula to find initiation as Column 2 $\times$ 0.95 + Column 3.

Column 4 gives us the yearly count of those individuals who reported to SAMHSA having entered a specialty treatment facility due to heroin use within the past year. Column 5 gives us the yearly count of those individuals who reported to SAMHSA having entered a specialty treatment facility due to pain medication use within the past year. Column 6 gives us the yearly count of those individuals with HUD due to heroin use who reported to SAMHSA having entered a specialty treatment facility within the past year. Column 7 gives us the yearly count of those individuals with OUD due to illicit pain medication use who reported to SAMHSA having entered a specialty treatment facility within the past year. Finally, column 8 gives us the yearly count of those individuals with OUD (due to heroin or pain medication) who reported to SAMHSA for 2016 and 2017, having entered a specialty treatment facility within the past year.

We use the data values in columns 4 through 8 to calculate data related to our state variables $T_E$ and $T$. The values for specialty treatment from $I$ for heroin are in column 6 of Table 2. We refer the reader to that section of the text to explain data calculation for the absent years. Using the values from column 6 of Table 4 and columns 7 and 8 of Table 6, we find the data for our state variable $T$. Given column 8 from Table 6; we found a formula to approximate the data for the missing years. This formula (Column 7 $\times$ 0.71 + Column 6) fills in the unaccounted-for values of column 8 (and is marked with an asterisk). For the only two years given by SAMHSA (2016 and 2017), the formula applied gives 453.0 (vs. 453) and 603.7 (vs. 603).

To find values for specialty treatment from opioids, we use the formula (Column 5$\times$ 0.83 + Column 4), where the factor of 0.83 is the average of the factor used to find the specialty treatment from $I$ data and the use in year data. (We have no data to compare; using 0.71 slightly altered the numbers (e.g., 316.5 vs. 343.6 in the year 2003), but the parameter estimation was primarily affected in the $\psi$ terms.) Then, with the values for the specialty treatment from opioids, we subtract the values of the $T$ data per year, and that result gives us the corresponding values for our state variable $T_E$.

The state variables $E$, $T_E$, $I$, and $T$ are instantaneous in time, whereas the SAMHSA data is not. SAMHSA gives a cumulative count of those respective classes. Therefore, we apply a correction when comparing the data to the variables (see Supplemental Information). We explain the computations for these corrections in the analogous heroin-only section. These explanations are precisely the same for comparing our model to the all-illicit opioids dataset.

As in heroin-only section, we consider a city with a population size of 200,000 individuals and scale the corresponding data for this set for the simulations. Parameter ranges used are the same as the values discussed in the analogous heroin-only section, except for $\omega$. We extended the range for this parameter from 0.05 to 1.3 because the SAMSHA data showed a higher number of individuals in the general treatment versus the specialty treatment facilities.

Analogous to the section for the heroin-only dataset, we define $\delta$ and $\delta_E$:

and

As in the heroin datasets, we note that the data used in the definition of the death rate suggest a piecewise function; see bottom middle subfigure of Figure 9. Incorporating these piecewise functions for $\delta$ and $\delta_E$ into our parameter estimation, our baseline values are

Our initial conditions, chosen to approximately pair with the scaled data are $S_0 = 199,600, E_0 = 3240, T_{E_0} = 48, I_0 = 769, T_0 = 71, R_0 = 325$. The data match is provided in Figure 9.

2.5.   Basic reproduction number

To better understand the dynamics of transmission, the basic reproduction number $\mathcal{R}_0$ is computed and analyzed. $\mathcal{R}_0$ is the number of secondary cases produced by one infectious individual introduced into a population of wholly susceptible individuals during their infectious period.

The $\mathcal{R}_0$ for the IOUD model with a casual user class (2.1) calculated using the next generation method as presented in Van den Driessche & Watmough (2002) ^[43] is as follows: $\mathcal{R}_0 = \mathcal{R}_1 + \mathcal{R}_2$, where

and

$\mathcal{R}_0^{{ SITR}}$ represents the basic reproduction number for the IOUD model without a casual user class in Cole and Wirkus (2022) ^[10] and duplicated here for convenience to the reader.

Shown in the flow diagram, Figure 1, is how the compartments $S$, $E$, and $T_E$ are coupled to $I$, $T$, and $R$, through $E$ going to $I$. We see this connection in the reproduction number. The terms $\mathcal{R}_1$ and $\mathcal{R}_2$ are analyzed analogously to "Reproduction numbers of infectious disease models" by Pauline van den Driessche (2017) ^[44]. If we introduce an individual into the $I$ class (i.e., one infected person into the population), then that individual's influence has two parts described as follows:

PART 1 ($\mathcal{R}_1$): The introduction of the infected person into the population may influence someone from $S$ into $E$. Therefore, this first part provides the contributions attributed to the $E$ class.

PART 2 ($\mathcal{R}_2$): This part provides the contributions attributed from the $I$ class as before in Cole and Wirkus (2022) ^[10], whereas the first factor describes the proportion of individuals from the $E$ class.

2.6.   Endemic equilibria

To proceed in determining the existence of non-trivial endemic equilibria of our IOUD model with a casual user class, the system is set to the steady-state population level, denoted by $N^*$. Since the total population is driven by $\frac{dN}{dt} = \Lambda -\mu N -\delta I -\delta_E E$, the steady-state population level is reached at $N^* = \frac{\Lambda-I\delta-\delta_E E}{\mu}$. This substitution results in the following system of equations:

where $b(T, T_E) = \dfrac{1}{1+\epsilon T + \epsilon_E T_E}.$

We numerically found parameter regimes of bistability, that is, where both Disorder-Free Equilibrium (DFE) and Endemic Equilibrium (EE) exist biologically and are stable for both heroin-only and illicit opioid datasets. We tried to obtain an explicit numerical approximation to the analytical curve separating the region of bistability from EE-only stability using the method from Cole and Wirkus (2022) ^[10] that successfully gave analytical curves for the regions of bistability: Using Maple, an equation is obtained by setting $\frac{d\tilde{S}}{dt} = 0$ and then solving for $S^*$; this result is substituted into $\frac{d\tilde{E}}{dt} = 0$, $\frac{d\tilde{T_E}}{dt} = 0$, $\frac{d\tilde{T}}{dt} = 0$ and $\frac{d\tilde{R}}{dt} = 0$; $E^*$, $T_E^*$, $T^*$ and $R^*$ are solved simultaneously. However, this step was not able to be solved by Maple. We use parameter values obtained from the heroin-only dataset, which one could refer to the baseline values (2.5), and then the all-illicit opioids dataset, which one could refer to the baseline values (2.8), for the ensuing investigation. We extrapolate the $\delta$ values and $\delta_E$ values for the two overdose death rates; see left panels of Figure 3 for the heroin-only dataset and Figure 4 for the all-illicit opioids dataset. Next, we determine the effective reproductive number $\mathcal{R}_{\mbox{eff}}(t) = (\mathcal{R}_0\cdot S(t)/N_0)$; see top right panel in Figures 3 and 4. For a range of overdose death rates, we numerically observed bistability. That is, with realistic parameter values and the $\delta$ and $\delta_E$ values extrapolated to future values; for example their 2026 values ($\delta\approx 0.0501$, $\delta_E\approx 0.0105$), we found both the DFE and an EE were stable. For the heroin-only dataset, we see that $\mathcal{R}_{\mbox{eff}}$ becomes less than 1 during 2024 and for the all-illicit opioids dataset $\mathcal{R}_{\mbox{eff}}$ becomes less than 1 in year 2044. Given the large number of parameters that have been estimated, which have natural temporal variability, we only extrapolate 10 years (to 2030) to lessen the chances of erroneous conclusions and generalizations. To this end, the bottom right subfigure of Figures 3 and 4 shows a plot of $\mathcal{R}_{\mbox{eff}}$ with an envelope of $\pm 10\%$ of the baseline parameter values where each significant parameter (determined by LHS-PRCC) was changed $+$ or $-$ 10% in the direction that it most influences $\mathcal{R}_{\mbox{eff}}$. Including a stochastic approach would be another way to take into account the randomness in future extensions of this work; see e.g., ^[45].

Figure 3.  Heroin death rate and effective reproductive number. (Top Left): Extrapolated $\delta$-values. The black X-marks and asterisks are from the heroin-only overdose data and are found in Figure 2; the extrapolated $\delta$-values and corresponding piecewise curve are represented in black. The magenta marks and piecewise curve represented are for the extrapolated $\delta$-values using the parameter values from the HUD model in Cole and Wirkus (2022) ^[10]. (Bottom Left): extrapolated $\delta_E$-values. The black X-marks and asterisks are from the heroin-only overdose data and are found in Figure 2; the extrapolated $\delta_E$-values and corresponding piecewise curve are represented in black. The magenta marks and piecewise curve represented are for the extrapolated $\delta_E$-values using the parameter values from the HUD model in Cole and Wirkus (2022) ^[10]. (Top Right): The effective reproductive number, $\mathcal{R}_{\mbox{eff}}(T) = (\mathcal{R}_0\cdot S(T)/N(T))$, is plotted as the solid black curve using the baseline values of the parameters of the heroin dataset, (2.5) and the extrapolated $\delta$-values are from the best fit line. Just above the $\mathcal{R}_{\mbox{eff}}$ curve, $\mathcal{R}_0$ is plotted as a dashed blue curve; also plotted in magenta is the curve from the original HUD model without the casual user class in Cole and Wirkus (2022) ^[10]. (Bottom Right): Parameters found for which $\mathcal{R}_0$ is most sensitive to are varied 10% in the direction of their greatest influence on $\mathcal{R}_0$ to give an envelope surrounding $\mathcal{R}_{\mbox{eff}}$.

DownLoad: Full-Size Img PowerPoint

Figure 4.  Illicit opioid death rate and effective reproductive number. (Top Left): Extrapolated $\delta$-values. The black X-marks and asterisks are from the all-illicit opioids overdose data and are found in Figure 9. The extrapolated $\delta$-values and curve are represented in magenta. (Bottom Left): Extrapolated $\delta_E$-values. The black X-marks and asterisks are from the all-illicit opioids overdose data and are found in Figure 9. (Top Right): The effective reproductive number, $\mathcal{R}_{\mbox{eff}}(t) = (\mathcal{R}_0\cdot S(t)/N(t))$, is plotted as the solid black curve using the baseline values of the parameters of the all-illicit opioids dataset, (2.8) and the extrapolated $\delta$-values from the best fit line. Just above the $\mathcal{R}_{\mbox{eff}}$ curve, $\mathcal{R}_0$ is plotted as a dashed blue curve. (Bottom Right): Parameters found for which $\mathcal{R}_0$ is most sensitive to are varied 10% in the direction of their greatest influence on $\mathcal{R}_0$ to give an envelope surrounding $\mathcal{R}_{\mbox{eff}}$.

DownLoad: Full-Size Img PowerPoint

3.   Sensitivity analysis

We execute a sensitivity analysis using the Partial Rank Correlation Coefficient (PRCC) methodology ^[46] to determine the input parameter's system's sensitivity and use Latin hypercube sampling (LHS) to allow for a $\pm 10\%$ uncertainty in the parameter inputs. For the study, we use the parameter values captured through the parameter estimation and the literature given in (2.5) for heroin and (2.8) for illicit opioids as our baseline values for 2020. We vary the parameters and initial conditions by $\pm 10\%$ from their baseline values.

3.1.   Heroin only analysis

Discussion of the PRCC values

There must be a monotonic relationship between the output values and model parameters when measuring sensitivity for the PRCC method ^[46]. Therefore, we performed monotonicity checks for all initial conditions and parameter values.

For the yearly number of casual users who enter the HUD class, a monotonic relationship for all variables and initial conditions was concluded from 2022 to 2040 for the constant and variable death rates.

For the yearly number of relapses from $T$ counts variable, plots were non-monotonic for several years for some of the parameters and initial conditions. For the constant and variable death rate, the yearly number of relapses from $T$ was not monotonic in the initial condition $T_E(0)$ from 2022 to 2024, and it was not monotonic in the parameter $\rho_E$ from 2028 to 2030; however, neither of these showed up as significant on the PRCC graphs for the relapse $T$ variable.

For the yearly number of relapses from $R$ counts variable, plots were non-monotonic for several years for some of the parameters and initial conditions. For both the constant and the variable death rate, the yearly number of relapses from $R$ was not monotonic in $\Lambda$ from 2026 to 2028; however, this parameter did not show up as significant on the PRCC graphs for the relapse $R$ variable. For the constant death rate, the yearly number of relapses from $R$ was not monotonic in $\epsilon$ from 2034 to 2036, but this parameter showed up as insignificant during this time period on the PRCC graph for the relapse $R$ variable. For the constant death rate, the yearly number of relapses from $R$ was not monotonic in $\rho_E$ from 2024 to 2026, but this parameter did not show up as significant on the PRCC graphs for the relapse $R$ variable. For the constant and variable death rate, the yearly number of relapses from $R$ was not monotonic in $S(0)$ from 2022 to 2024; on the other hand, this parameter did not show up as significant on the PRCC graphs for the relapse $R$ variable. For the constant death rate, the yearly number of relapses from $R$ was not monotonic in $\epsilon_E$ from 2034 to 2036; however, this parameter did not show up as significant on the PRCC graphs for the relapse $R$ variable. For the constant and variable death rate, the yearly number of relapses from $R$ was not monotonic in $T_E(0)$ from 2022 to 2024, but this parameter did not show up as significant on the PRCC graphs for the relapse $R$ variable. For the variable death rate, the yearly number of relapses from $R$ was not monotonic in the parameter $\epsilon$ from 2036 to 2038; on the other hand, this parameter did not show up as significant during this time period on the PRCC graph for the relapse $R$ variable. For the variable death rate, the yearly number of relapses from $R$ was not monotonic in the parameter $\rho_E$ from 2036 to 2038; however, this parameter did not show up as significant on the PRCC graph for the relapse $R$ variable. For the variable death rate, the yearly number of relapses from $R$ was not monotonic in the parameter $\epsilon_E$ from 2024 to 2026 for the relapse $R$ variable; however, this parameter did not show up as significant on the PRCC graphs for the relapse $R$ variable.

Finally, we checked the monotonicity results for the yearly number of opioid overdose deaths. For the constant and variable death rate, the yearly number of opioid overdose deaths was not monotonic in $\Lambda$ from 2022 to 2024, although this parameter did not show up as significant on the PRCC graphs for the yearly number of opioid overdose deaths. For the variable death rate, the yearly number of opioid overdose deaths was not monotonic in the parameter $b$ from 2036 to 2038, but this parameter showed up as not significant during this time period on the PRCC graph for the yearly number of opioid overdose deaths.

As the theory requires, we do not use results of significance for the parameters in the years where monotonicity fails. Similarly, in the ensuing discussion, we don't consider the variables for the parameters in the years when monotonicity fails.

This following discussion presents variables of interest to the healthcare industry and policymakers. The focus will be on the yearly number of casual users who enter the HUD class for the first time, the yearly number of individuals who relapse from the $T$ class, the yearly number of individuals who relapse from the $R$ class, and the number of yearly opioid overdose deaths due to heroin. Although these variables are not of the original system of equations, we calculate them by keeping track of their cumulative yearly totals.

Four graphs for each case correspond to the sensitivities for the constant death rates (at their 2020 values) in 2030 and 2040 ($\delta = 0.0343$ and $\delta_E$ = 0.0078) versus the variable death rate in 2030 and 2040.

As was done in Cole and Wirkus (2022) ^[10], but duplicated here for the reader's convenience, we will refer to the sensitivity results as being "highly significant" if it has a PRCC value of 0.85 or higher, "significant" if it has a PRCC value of 0.70–0.84, "somewhat significant" if it has a PRCC value of 0.55–0.69, "slightly significant" if it has a PRCC value of 0.45–0.54, "borderline significant" if it has a PRCC value of 0.40–0.44, and "not significant" if it has a PRCC value of under 0.40. The significance of the initial conditions will also not be discussed as reasoned in Cole and Wirkus (2022) ^[10]. Additionally, only parameters that may be changed due to external influence will be discussed.

Yearly new I from E:

The variable Yearly new $I$ from $E$ gives the count of the number of casual users from the $E$ class who entered the $I$ (HUD) class; see Figure 5 and Table 7. The comparisons of the PRCC values graphs are similar, and the discussion will be relevant for all four unless otherwise noted. The analysis ranked three parameters highly significant. $\beta$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $I$) is positively correlated. An increase in this transmission rate would cause an increase in the number of individuals who transition into the casual user $E$ class, as expected. $\chi$ (the rate of individuals in $E$ that transition to $I$) is also positively correlated. An increase in this rate would also cause an increase in the number of individuals entering the $I$ class, as expected. Not only are more individuals entering $I$, but there are also more individuals interacting with $S$ to influence them into $E$. $\zeta$ (The rate of individuals in $E$ returning to $S$) is negatively correlated. Hence, as expected, lowering this rate would decrease the number of individuals entering the $I$ class. Thus, it is recommended in the short and long term to reduce the rate of transmission of those casual users leading to IOUD and increase the rate that casual users stop using by entering back into the $S$ class. At the year-end of 20 years, two parameters ranked somewhat significant, $\delta$ (HUD overdose death rate) and $b$ (one of the parameters for the variable death rate). They are negatively correlated, so increasing these rates would decrease the number of individuals entering the $I$ class. However, we ethically would not want the individuals entering the $I$ class to drop in this manner; therefore, it would be beneficial to concentrate on the other ones.

Figure 5.  Heroin Data: PRCC results over time for those who are entering $I$ for the first time, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 7. The left figures have a final time of 2030 whereas the right figures have a final time of 2040. The top figures keep $\delta$ and $\delta_E$ constant at their 2020 values whereas the bottom figures use the extrapolation functions for $\delta$ and $\delta_E$.

DownLoad: Full-Size Img PowerPoint

Table 7.  Heroin Only Data: PRCC results for movement into $I$, relapse from $T$, relapse from $R$, and yearly deaths using the baseline parameters and initial conditions and using either the constant $\delta$ or the variable $\delta$. The PRCC values are given at year end time of 2030 and year end time of 2040. Table values without an entry are not significant or undefined (in the case of $m$ and $b$ for the constant death rate and $\delta$ and $\delta_E$ for the variable death rate). The corresponding graphs for this Table are given in Figures 5–8.

Param	Yearly new $I$ from $E$				Yearly relapse $T$				Yearly relapse $R$				Yearly Deaths
	Constant		Variable		Constant		Variable		Constant		Variable		Constant		Variable
	2030	2040	2030	2040	2030	2040	2030	2040	2030	2040	2030	2040	2030	2040	2030	2040
$\mu$	-	-	-	-0.46	-	-	-	-0.53	-0.45	-0.4	-0.49	-0.61	-0.43	-	-0.49	-0.63
$\beta$	0.97	0.96	0.97	0.98	0.81	0.93	0.87	0.96	0.94	0.94	0.87	0.97	0.94	0.93	0.94	0.98
$\delta$	-0.46	-0.59	-	-	-	-0.66	-		-0.64	-0.68	-	-	0.93	0.5	-	-
m	-	-	-	-0.52	-	-	-	-0.58	-	-	-	-0.67	-	-	0.74	0.44
b	-	-	-0.47	-0.69	-	-	-	-0.71	-	-	-0.57	-0.8	-	-	0.74	-
$\Lambda$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
k	-	-	-	-	0.92	0.89	0.95	0.91	-0.43	-	-	-0.41	0.49	-	0.47	-
$\rho$	-	-	-	-	-0.64	-0.63	-0.64	-0.59	0.89	0.52	0.8	0.8	-	-	-	-
$\eta_1$	-	-	-	-	0.78	0.65	0.79	0.73	0.42	-	0.41	-	-	-	-0.41	-
$\eta_2$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\eta_3$	-	-	-	-	0.53	0.53	0.62	0.49	-	-	-	-	-	-	-	-
$\alpha_1$	0.43	-	0.45	0.47	-	-	0.5	0.49	0.7	0.5	0.59	0.63	0.61	-	0.64	0.57
$\alpha_2$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\omega$	-	-	-	-	-	-	-	-	0.86	0.58	0.8	0.78	-	-	-0.42	-
$\epsilon$	0.43	-	0.51	0.55	-0.66	-0.53	-0.72	-0.52	-	-	-	-	0.52	0.42	0.56	0.65
$\beta_E$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\delta_E$	-	-	-	-	-	-	-	-	-	-		-	-	-	-	-
mE	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
bE	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
k_E	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\rho_E$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\psi_1$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\psi_2$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\psi_3$	-	-	-	-	-	-	-	-	-	-	-	-0.44	-	-	-	-
$\zeta$	-0.93	-0.9	-0.93	-0.95	-0.68	-0.86	-0.72	-0.9	-0.85	-0.86	-0.77	-0.94	-0.87	-0.85	-0.88	-0.94
$\chi$	0.95	0.93	0.94	0.97	0.8	0.91	0.83	0.94	0.9	0.9	0.85	0.96	0.9	0.88	0.9	0.96
$\epsilon_E$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$S(0)$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$E(0)$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$T_E(0)$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$I(0)$	0.86	0.65	0.85	0.85	0.84	0.79	0.9	0.87	0.94	0.79	0.91	0.92	0.94	0.74	0.94	0.91
$T(0)$	0.46	-	-	-	0.42	0.47	0.48	0.43	0.64	-	0.56	0.55	0.65	-	0.6	0.51
$R(0)$	-	-	-	-	-	-	-	-	0.61	-	-	0.41	0.53	-	0.53	-

---

 | Show Table

DownLoad: CSV

Yearly relapse T:

The variable yearly relapse $T$ gives the count of the individuals who relapsed from the $T$ class back to the $I$ class; see Figure 6 and Table 7. The graphs for the year end of 2030 for both death rates are similar. The parameter $\kappa$ (rate of individuals leaving treatment and returning to $I$) is ranked highly significant. Since it is positively correlated, increasing this rate will increase the number of individuals who relapse from $T$ back to $I$. The parameters $\beta$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $I$), $\chi$ (the rate of individuals in $E$ that transition to $I$), and $\eta_1$ (rate of individuals in $I$ who enter specialty treatment on their own) are ranked significant and are positively correlated. Hence, a decrease in these rates would cause a reduction in the number of individuals who relapse from $T$. Although we want to see a drop, we still wish for individuals to enter into treatment even if they may retreat to $I$; hence, we do not consider it beneficial to decrease $\eta_1$. The parameters $\zeta$ (rate of individuals in $E$ returning to $S$), $\epsilon$ (saturation term for entering a specialty treatment facility), and $\rho$ (rate of individuals leaving treatment and entering the recovered class) came up as somewhat significant and all are negatively correlated. Thus, an increase in these parameters would decrease the yearly relapse $T$ counts. Since an increase in $\epsilon$ would lower the limit of available treatment facilities, we would not focus on this parameter since we want individuals to get into treatment. Hence, in the short term of ten years, decreasing the relapse rate from $T$, decreasing the rate of casual users ending up with OUD, reducing the transmission rate of the HUD class, increasing the rate of casual users returning to the $S$ class, and increasing the rate of individuals who complete treatment are all beneficial avenues for decreasing the yearly relapse $T$ counts. The graphs for the year end of 2040 are similar with a few differences. The significance of parameters $\beta$ and $\zeta$ increased, and they were the most significant parameters. Although $\kappa$ decreased in relevance, it remained highly influential. As with the variable yearly new $I$ from $E$ variable, the parameters $\delta$ and $b$ were ranked as somewhat significant in the long term and negatively correlated. Similarly, we do not want death by overdose to decrease the counts.

Figure 6.  Heroin: PRCC results over time for the number of those individuals who relapsed from $T$ and went back to $I$, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 7. The left figures have a final time of 2030 whereas the right figures have a final time of 2040. The top figures keep $\delta$ and $\delta_E$ constant at their 2020 values whereas the bottom figures use the extrapolation functions for $\delta$ and $\delta_E$.

DownLoad: Full-Size Img PowerPoint

Yearly relapse R:

The variable yearly relapse $R$ gives the count of the individuals who relapsed from the $R$ class back to the $I$ class; see Figure 7 and Table 7. The graphs are similar for both death rates for 2030 and 2040. The highly significant parameters are $\beta$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $I$) and $\chi$ (rate of individuals in $E$ that transition to $I$); both are positively correlated. Hence, as expected, increasing these rates would increase the number of individuals who relapse to $I$ from the recovered class. Increasing those values would result in an overall increase in the number of individuals who would enter $I$ and then possibly flow into the recovered class either directly or indirectly through a specialty treatment facility. The parameter $\rho$ (rate of individuals leaving specialty treatment and entering the recovered class) ranked as significant to highly significant. As expected, since it is positively correlated, increasing this rate would increase the number of yearly relapse $R$ counts. The parameter $\omega$ (rate of individuals in $I$ who enter $R$ by either completing treatment in non-specialty facilities or "quitting cold turkey") ranked significant. Positively correlated, a decrease in those directly entering the $R$ class from $I$ would decrease the Yearly relapse $R$ counts. However, this reduction would reduce the overall number of individuals in $R$. The goal is always to move individuals out of the $I$ class in beneficial ways, even if the relapse count could be higher. Hence, we do not consider decreasing $\omega$. The parameter $\zeta$ (rate of individuals in $E$ returning to $S$) ranked as significant. Negatively correlated, increasing this parameter would reduce the yearly relapse $R$ counts. For the year-end of 20 years and the variable death rate, this parameter increased its significance over time to a ranking of highly significant. The parameter $\alpha_1$ (rate of individuals in $R$ relapsing to $I$ on their own accord) ranked as somewhat influential. Since this parameter is positively correlated, decreasing the rate reduces the yearly relapse $R$ counts. Hence, in the short and long term, it would be best to focus on reducing the HUD class transmission rate, lowering the rate of individuals entering $I$ from $E$, increasing the rate of return from $E$ to $S$, and decreasing the relapse rate of $R$ for those relapsing on their own accord.

Figure 7.  Heroin: PRCC results over time for the number of those individuals who relapsed from $R$ and went back to $I$, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 7. The left figures have a final time of 2030 whereas the right figures have a final time of 2040. The top figures keep $\delta$ and $\delta_E$ constant at their 2020 values whereas the bottom figures use the extrapolation functions for $\delta$ and $\delta_E$.

DownLoad: Full-Size Img PowerPoint

As with the previous variables, the parameters $\delta$ and $m$ showed somewhat significant in the long run and were negatively correlated. Additionally, $b$ ranked as sensitive. We again reiterate that we do not want death by overdose to be the reason for a decrease in the counts.

Yearly deaths:

The variable yearly deaths count the number of HUD overdose deaths from the $E$ class and the $I$ class; see Figure 8 and Table 7. The comparisons of the PRCC values graphs are similar, and the discussion will be relevant for all four plots unless otherwise noted. Three parameters ranked highly significant. First, $\beta$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $I$) is positively correlated. An increase in this transmission rate would cause an increase in the number of HUD overdose deaths, as expected. Second, $\chi$ (rate of individuals in $E$ that transition to $I$) is also positively correlated. An increase in this rate would also cause an increase in the number of HUD overdose deaths, as hypothesized. Finally, $\zeta$ (rate of individuals in $E$ returning to $S$) is negatively correlated. Hence, as one would predict, decreasing this rate would decrease the number of HUD overdose deaths. Therefore, it is recommended in the short and long term to reduce the HUD transmission rate, decrease the rate of casual users with OUD, and increase the rate that casual users stop using and enter back into the $S$ class.

Figure 8.  Heroin: PRCC results over time for the number of yearly deaths due to illicit opioid use, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 7. The left figures have a final time of 2030 whereas the right figures have a final time of 2040. The top figures keep $\delta$ and $\delta_E$ constant at their 2020 values whereas the bottom figures use the extrapolation functions for $\delta$ and $\delta_E$.

DownLoad: Full-Size Img PowerPoint

Figure 9.  Fitting model output to scaled data (with error bars when given) for all-illicit opioids. The red squares are the SAMHSA data and the blue curves are the model output as described in the text. (Top Left): Illicit opioid overdose deaths by year. (Top Middle) Individuals with illicit opioid use disorder (IOUD) in the last year. (Top Right) Those who entered specialty treatment from IOUD in the past year. (Middle Left) Individuals who used illicit opioids in the last year. (Middle Middle) Those who entered specialty treatment from the illicit opioid casual user class, $E$. (Middle Right) Those who reported using illicit opioids for the first time (initiation from $S$ to $E$). (Bottom Left) The effective reproductive number, $\mathcal{R}_{\mbox{eff}}(t) = (\mathcal{R}_0\cdot S(t)/N(t))$. (Bottom Middle) Overdose death rate for IOUD class: asterisks and X-marks are calculated from data (see text and Eq (2.6)). Both lines are calculated with a least squares fit. (Bottom Right) Overdose death rate for $E$ Class: asterisks and X-marks are calculated from data (see text and Eq (2.7)).

DownLoad: Full-Size Img PowerPoint

The parameter $\delta$ (HUD overdose death rate) ranked highly significant for the constant death rate at year-end of ten years. Since it was positively correlated, an increase in the death rate would increase the number of yearly death counts, as expected. However, at the year-end of 20 years, for the constant death rate, $\delta$ showed up as only slightly significant. Counter-intuitively, the significance of this parameter decreased as time went on. An interpretation would be that the consequence of a high death rate leading to a higher number of overdose deaths led to fewer users in the $I$ class over time. As a result, there are fewer individuals in $I$ to interact with susceptibles. Therefore, this parameter is less sensitive in the long run because the number of influenced individuals decreases, leading to an overall decrease in the $I$ population and fewer yearly deaths. This case advocates the urgency to expedite users out of the $I$ class into treatment to protect them from the high overdose death rate. This circumstance is the same for the variable death rate and the parameters $m$ and $b$. In the short term, these parameters were significant, although they were not in the long term. This discovery could also be why $\epsilon$ ranked as somewhat significant for the year-end of 20 years and the variable death rate. This significance increases from the 10-year mark. Positively correlated, as one would predict, increasing this parameter increases the yearly deaths. $\epsilon$ is inversely proportional to the availability of specialty treatment facilities, and increasing this parameter will decrease the availability of an individual to get care.

3.2.   All-illicit opioids analysis

Discussion of the PRCC values

There must be a monotonic relationship between the output values and model parameters when measuring sensitivity for the PRCC method. Therefore, we performed monotonicity checks for all initial conditions and parameter values.

For the yearly number of casual users who enter the IOUD class, a monotonic relationship for all variables and initial conditions was concluded from 2022 to 2040 for the constant and variable death rates.

For the yearly number of relapses from T counts variable, plots were non-monotonic for several years for some of the parameters and initial conditions. For the constant and variable death rates, the yearly number of relapses from $T$ was not monotonic in parameter $\beta_E$ from 2024 to 2025; however, this parameter did not show up as significant during this time period on the PRCC graphs for the yearly number of relapses from $T$. For the constant and variable death rates, the yearly number of relapses from $T$ was not monotonic in parameter $\zeta$, from 2024 to 2025; however, this parameter did not show up as significant during this time period on the PRCC graphs for the yearly number of relapses from $T$. For the constant death rate, the yearly number of relapses from $T$ was not monotonic in parameter $\beta$, from 2024 to 2025; however, this parameter did not show up as significant during this time period on the PRCC graphs for the yearly number of relapses from $T$. For the constant death rate, the yearly number of relapses from $T$ was not monotonic in the initial condition $T_E(0)$ from 2024 to 2026; however, this parameter did not show up as significant on the PRCC graphs for the yearly number of relapses from $T$. For the variable death rate, the yearly number of relapses from $T$ was not monotonic in the parameter $\rho_E$ from 2024 to 2026; however, this parameter did not show up as significant on the PRCC graphs for the yearly number of relapses from $T$.

For the yearly number of relapses from R counts variable, plots were non-monotonic for several years for some of the parameters and initial conditions. For the constant death rate, the yearly number of relapses from $R$ was not monotonic in the parameters $\kappa_E$, $\psi_1$, and $\psi_3$ from 2022 to 2023; however, these parameters did not show up as significant on the PRCC graphs for the yearly number of relapses from $R$. For the constant death rate, the yearly number of relapses from $R$ was not monotonic in parameter $\omega$ from 2036 to 2037; however, this parameter did not show up as significant during this time period on the PRCC graph for the yearly number of relapses from $R$. For the variable death rate, the yearly number of relapses from $R$ was not monotonic in parameters $\kappa_E$, $\rho_E$, $\psi_1$, $\psi_2$, and $\psi_3$ from 2022 to 2024; however, these parameters did not show up as significant on the PRCC graphs for the yearly number of relapses from $R$. For the constant death rate, the yearly number of relapses from $R$ was not monotonic in parameter $\omega$ in 2040. However, this parameter did not show up as significant during this time on the PRCC graph for the yearly number of relapses from $R$.

Finally, we checked the monotonicity results for the yearly opioid overdose deaths. For the constant death rate, the yearly number of opioid overdose deaths was not monotonic in parameter $\kappa_E$ from 2024 to 2025; however, this parameter did not show up as significant on the 2030 PRCC graph and showed up as insignificant during this time period on the 2040 PRCC graph for the yearly number of opioid overdose deaths. For the constant death rate, the yearly number of opioid overdose deaths was not monotonic in the parameters $\psi_1$, $\psi_2$, and $\psi_3$ from 2024 to 2025; however, these parameters did not show up as significant on the PRCC graphs for the yearly number of opioid overdose deaths. For the variable death rate, the yearly number of opioid overdose deaths was not monotonic in parameters $\kappa_E$, $\rho_E$, $\psi_1$, $\psi_2$, and $\psi_3$ from 2024 to 2026; however, this parameter did not show up as significant on the PRCC graphs for the yearly number of opioid overdose deaths.

As the theory requires, we do not use results of significance for the parameters in the years where monotonicity fails. Similarly, in the ensuing discussion, we don't consider the variables for the parameters in the years when monotonicity fails.

The following discussion presents variables of interest to the healthcare industry and policymakers. The focus will be on the yearly number of the casual users who enter the IOUD class for the first time, the yearly number of individuals who relapse from the $T$ class, the yearly number of individuals who relapse from the $R$ class, and the yearly opioid overdose deaths due to all-illicit opioids. Although these variables are not of the original system of equations, we calculate them by keeping track of their cumulative yearly totals.

Four graphs for each case correspond to the sensitivities for the constant death rates (at their 2020 values) in 2030 and 2040 ($\delta = 0.0343$ and $\delta_E$ = 0.0078) versus the variable death rate in 2030 and 2040.

We will refer to the sensitivity results as being "highly significant" if it has a PRCC value of 0.85 or higher, "significant" if it has a PRCC value of 0.70–0.84, "somewhat significant" if it has a PRCC value of 0.55–0.69, "slightly significant" if it has a PRCC value of 0.45–0.54, "borderline significant" if it has a PRCC value of 0.40–0.44, and "not significant" if it has a PRCC value of under 0.40. The significance of the initial conditions will also not be discussed as reasoned in Cole and Wirkus (2022) ^[10]. Additionally, only parameters that may be changed due to external influence will be discussed.

Yearly new I from E:

The variable yearly new $I$ from $E$ gives the count of the number of casual users from the $E$ class who entered the $I$ (all-illicit opioids) class; see Figure 10 and Table 8. The comparisons of the PRCC values graphs are similar, and the discussion will be relevant for all four plots unless otherwise noted. The two parameters $\beta_E$ and $\zeta$ ranked highly significant. Since the parameter $\beta_E$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $E$) was positively correlated, increasing this rate would increase the number of yearly new $I$ from $E$ counts. On the other hand, since the parameter $\zeta$ (rate of individuals in $E$ returning to $S$) is negatively correlated, increasing this rate would decrease the number of yearly new $I$ from $E$ counts. Finally, the parameter $\beta$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $I$) ranked as somewhat significant to significant. Since it is positively correlated, increasing this parameter would increase the number of yearly new $I$ from $E$ counts. Therefore, decreasing the transmission rates, with more focus on $\beta_E$, plus increasing the rate of individuals going back to the $S$ class from $E$ would be very beneficial to drive down the yearly new $I$ from $E$ counts.

Figure 10.  Illicit opioids: PRCC results over time for those who are entering $I$ for the first time, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 8. The left figures have a final time of 2030 whereas the right figures have a final time of 2040. The top figures keep $\delta$ and $\delta_E$ constant at their 2020 values whereas the bottom figures use the extrapolation functions for $\delta$ and $\delta_E$.

DownLoad: Full-Size Img PowerPoint

Table 8.  All illicit opioids data: PRCC results for movement into $I$, relapse from $T$, relapse from $R$, and yearly deaths using the baseline parameters and initial conditions and using either the constant $\delta$ or the variable $\delta$. The PRCC values are given at year end time of 2030 and year end time of 2040. Table values without an entry either are not significant or undefined (in the case of $m$ and $b$ for the constant death rate and $\delta$ and $\delta_E$ for the variable death rate). The corresponding graphs for this table are given in Figures 10–13.

Param	Yearly new $I$ from $E$				Yearly relapse $T$				Yearly relapse $R$				Yearly Deaths
	Constant		Variable		Constant		Variable		Constant		Variable		Constant		Variable
	2030	2040	2030	2040	2030	2040	2030	2040	2030	2040	2030	2040	2030	2040	2030	2040
$\mu$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\beta$	0.6	0.66	0.57	0.45	-	0.48	-	-	-	0.58	0.52	0.45	-	0.61	0.47	0.43
$\delta$	-	-	-	-	-	-	-	-	-	-	-	-	0.7	0.69	-	-
m	-	-	-	-	-	-	-	-	-	-	-	-	-	-	0.41	0.53
b	-	-	-	-	-	-	-	-	-	-	-	-	-	-	0.51	-
$\Lambda$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
k	-	-	-	-	0.73	0.62	0.71	0.45	-	-	-	-	-	-	-	-
$\rho$	-	-	-	-	-0.62	-	-0.65	-0.48	-	-	-	-	-	-	-	-
$\eta_1$	-	-	-	-	0.75	0.61	0.7	0.61	-	-	-	-	-	-	-	-
$\eta_2$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\eta_3$	-	-	-	-	0.68	0.5	0.67	0.45	-	-	-	-	-	-	-	-
$\alpha_1$	-	0.4	-	-	0.81	0.75	0.8	0.7	0.87	0.75	0.87	0.69	0.67	0.71	0.72	0.65
$\alpha_2$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\omega$	-	-0.48	-	-0.41	-0.81	-0.77	-0.82	-0.72	-	-	-	-	-0.7	-0.71	-0.75	-0.66
$\epsilon$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\beta_E$	0.98	0.98	0.97	0.97	0.88	0.96	0.88	0.95	0.96	0.97	0.95	0.96	0.95	0.97	0.95	0.97
$\delta_E$	-	-	-	-	-	-	-	-	-	-		-	-	-	-
$m_E$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$b_E$	-	-	-	-	-	-	-	-			-	-			-	-
$k_E$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\rho_E$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\psi_1$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\psi_2$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\psi_3$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$\zeta$	-0.98	-0.98	-0.97	-0.97	-0.86	-0.96	-0.88	-0.94	-0.95	-0.97	-0.96	-0.96	-0.95	-0.97	-0.95	-0.96
$\chi$	-	-	-	-	0.59	0.5	0.55	0.4	0.44	-	-	-	-	-	-	-
$\epsilon_E$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$S(0)$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$E(0)$	-	-	-	-	-	-	0.43	-	-	-	-	-	-	-	-	-
$T_E(0)$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$I(0)$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$T(0)$	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
$R(0)$	-	-	-	-	0.68	0.46	0.6	0.45	0.75	-	0.7	0.49	0.53	-	0.47	0.45

---

 | Show Table

DownLoad: CSV

Yearly relapse T:

The variable yearly relapse $T$ gives the count of the individuals who relapsed from the $T$ class back to the $I$ class; see Figure 11 and Table 8. The comparisons of the PRCC values graphs are similar, and the discussion will be relevant for all four plots unless otherwise noted. The two parameters $\beta_E$ and $\zeta$ were highly significant. Since the parameter $\beta_E$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $E$) is positively correlated, increasing this rate would increase the number of yearly relapse $T$ counts, as expected. Since the parameter $\zeta$ (rate of individuals in $E$ returning to $S$) is negatively correlated, an increase in this rate would decrease the number of yearly relapse $T$ counts, as expected. Two parameters, $\omega$ (rate of individuals in $I$ who enter $R$ by either completing treatment in non-specialty facilities or "quitting cold turkey") and $\alpha_1$ (rate of individuals in $R$ relapsing to $I$ on their own accord) came up as significant. Since $\omega$ is negatively correlated, increasing this rate decreases the yearly relapse $T$ counts. This decrease happens because fewer individuals would go to the $T$ class. However, we disregard this parameter because we want individuals to transition out of $I$ whether they go to $R$ or $T$. Since $\alpha_1$ is positively correlated, a decrease in this rate would decrease the yearly relapse $T$ counts, which would be agreeable. The parameters $\eta_1$ (rate of individuals in $I$ who enter specialty treatment on their own) and $\kappa$ (rate of individuals leaving treatment and returning to $I$) ranked as somewhat significant to significant. Since $\eta_1$ is positively correlated, decreasing this rate would reduce the number of yearly relapses from $T$ counts because more individuals would be in the $T$ class. Although we want to decrease the number of relapses, we still wish for individuals to go to treatment Therefore, we disregard this parameter. Since $\kappa$ is positively correlated, decreasing this rate would reduce the number of yearly relapses from $T$ counts, as expected and agreeable. The parameters $\eta_3$ (rate of individuals in $I$ who enter $T$ through interaction with a susceptible), $\rho$ (rate of individuals leaving specialty treatment and entering the recovered class), and $\chi$ (rate of individuals in $E$ that transition to $I$) showed up as somewhat significant. Since $\eta_3$ is positively correlated, increasing this rate would increase the yearly relapse $T$ counts. Following the same reasoning as $\eta_1$, we disregard this parameter. Since the parameter $\rho$ is negatively correlated, an increase in this rate would decrease the number of yearly relapse $T$ counts, as expected and agreeable. Since $\chi$ is positively correlated, lowering this rate would reduce the number of yearly relapse $T$ counts, as expected and agreeable. First, it is recommended to focus on starting with, most importantly, decreasing the casual user transmission rate and increasing the rate of users returning to the $S$ class from the $E$ class. Then, it is beneficial to decrease the relapse rates from the $R$ class followed by the $T$ class, increase the treatment completion rate from the $T$ class, and decrease the rate of those casual users entering the $I$ class.

Figure 11.  Illicit opioids: PRCC results over time for the number of those individuals who relapsed from $T$ and went back to $I$, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 8. The left figures have a final time of 2030 whereas the right figures have a final time of 2040. The top figures keep $\delta$ and $\delta_E$ constant at their 2020 values whereas the bottom figures use the extrapolation functions for $\delta$ and $\delta_E$.

DownLoad: Full-Size Img PowerPoint

Yearly relapse R:

The variable yearly relapse $R$ gives the count of the individuals who relapsed from the $R$ class back to the $I$ class; see Figure 12 and Table 8. The comparisons of the PRCC values graphs are similar, and the discussion will be relevant for all four plots unless otherwise noted. The two parameters $\beta_E$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $E$) and $\zeta$ (rate of individuals in $E$ returning to $S$) came up as highly significant. Since the parameter $\beta_E$ is positively correlated, increasing this rate would increase the number of yearly relapse $R$ counts. On the other hand, since the parameter $\zeta$ is negatively correlated, an increase in this rate would decrease the number of relapse $R$ counts. The parameter $\alpha_1$ (rate of individuals in $R$ relapsing to $I$ on their own accord) came up as significant in the short term but then decreased to somewhat significant in the long run for both death rates. Since $\alpha_1$ is positively correlated, lowering this parameter would reduce the number of relapse $R$ counts. Hence in the short term, it is recommended to focus on decreasing the casual user transmission rate, increasing the rate casual users go back to the susceptibles, and reducing the relapse rate of individuals from the $R$ class on their own. Finally, in the long run, the focus should be more concentrated on the casual user transmission rate and the rate of the casual users back to $S$.

Figure 12.  Illicit opioids: PRCC results over time for the number of those individuals who relapsed from $R$ and went back to $I$, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 8. The left figures have a final time of 2030 whereas the right figures have a final time of 2040. The top figures keep $\delta$ and $\delta_E$ constant at their 2020 values whereas the bottom figures use the extrapolation functions for $\delta$ and $\delta_E$.

DownLoad: Full-Size Img PowerPoint

Yearly deaths:

The variable yearly death gives the count of the number of opioid overdose deaths from the $E$ class and the $I$ class; see Figure 13 and Table 8. Since the comparisons of the PRCC values graphs are similar, the discussion will be relevant for all four plots unless otherwise noted. The two parameters $\beta_E$ (transmission rate of moving to $E$ from $S$ through interaction with someone from $E$) and $\zeta$ (rate of individuals in $E$ returning to $S$) came up as highly significant. Since the parameter $\beta_E$ is positively correlated, a decrease in this rate would decrease the number of yearly death counts. Since the parameter $\zeta$ is negatively correlated, an increase in this rate would reduce the number of yearly death counts. Two parameters, $\omega$ (rate of individuals in $I$ who enter $R$ by either completing treatment in non-specialty facilities or "quitting cold turkey") and $\alpha_1$ (rate of individuals in $R$ relapsing to $I$ on their own accord) came up as significant. Since $\omega$ is negatively correlated, increasing this rate will decrease the yearly death count. Since $\alpha_1$ is positively correlated, reducing this rate would lower the yearly death count. Therefore, we recommend focusing on starting with, most importantly, decreasing the casual user transmission rate and increasing the rate casual users go back to the susceptible class. Then, increasing the rate of individuals going to the recovered class from $I$ and decreasing the relapse rate of the $R$ class on one's own.

Figure 13.  Illicit opioids: PRCC results over time for the number of yearly deaths, with grayed region denoting a lack of significance. These results are summarized in the text and in Table 8. The left figures have a final time of 2030 whereas the right figures have a final time of 2040. The top figures keep $\delta$ and $\delta_E$ constant at their 2020 values whereas the bottom figures use the extrapolation functions for $\delta$ and $\delta_E$.

DownLoad: Full-Size Img PowerPoint

The sensitivity results for the overdose death rate parameters are similar to the heroin-only dataset. Although in the beginning, the parameters $\delta$ and $b$ display high significance as time goes on, this significance drops. A high death rate leading to a higher number of overdose deaths led to fewer users in the $I$ class over time. As a result, there are fewer individuals in $I$ to interact with susceptibles. Therefore, this parameter is less sensitive in the long run because the number of individuals being influenced decreases, leading to an overall reduction in the $I$ population and fewer yearly deaths. This case advocates the urgency to expedite users out of the $I$ class into treatment to be protected from the high overdose death rate.

4.   Summary and discussion

We presented a model for the heroin epidemic dynamics and illicit opioid use epidemic dynamics. We do not intend to assert that these results apply to the more general opioid epidemic, which was initially driven by legitimate prescriptions. In Cole and Wirkus (2022) ^[10], we illustrate one such model called the IOUD model, which features four classes: susceptibles, IOUD class, specialty treatment facilities, and recovered and use datasets from heroin use. Here we additionally extend this to include illicit opioid use (IOUD). There are multiple ways an individual could cycle among the classes. For example, one could cycle from IOUD to treatment to recovery and then back to the IOUD class. However, once an individual had IOUD, they could no longer cycle back to the susceptibles. This model also featured a saturation treatment function which limits the flow into the specialty treatment facilities from the IOUD class due to the limited availability of care. We extended the IOUD model in Cole and Wirkus (2022) ^[10] to include a casual user class and a corresponding specialty treatment facilities class. We referred to this new version as the IOUD model with a casual user class.

We found realistic parameter values through the literature and parameter estimation and matched them to the CDC and SAMHSA data. The IOUD model with a casual user class displayed linearly increasing overdose death rates. This increase started in 2011 for the HUD overdose death rates, $\delta$ and $\delta_E$ using the heroin-only dataset. However, this increase started in 2013 using the all-illicit opioids dataset. On the basis that the IOUD model with a casual user class approaches constant population $N^* = (\Lambda- \delta I^*- \delta_E E^*)/\mu$, we scaled the SAMHSA data to a population of 200,000 (ignoring the overdose death rate). Scaling in such a way permitted us to enhance our understanding of the heroin/illicit opioids epidemic dynamics.

With the parameter estimates, we determined that extrapolated $\delta$-values resulted in bistability. In this region, although the effective reproductive number, $\mathcal{R}_{\mbox{eff}}$, is less than one, we found both the DFE and EE stable. For the extended IOUD model, with realistic parameter values for heroin data and the $\delta$ and $\delta_E$ values extrapolated to future values; for example their 2026 values ($\delta\approx 0.0501, \delta_E \approx 0.0205$), we found both the DFE and an EE were stable; and for sufficiently large values of $\delta$ and $\delta_E$ (the late-2045 values of $0.1025$ and $0.042$, respectively), we found that only the DFE was stable. This discovery of backward bifurcation emphasizes a complication for eliminating HUD. For the IOUD model with the heroin-only dataset, there is a minimum threshold $\epsilon$ value below which we did not have bistability. In other words, increasing accessibility to specialty treatment facilities is vital to ending this epidemic. In addition, including the casual user class also appears to increase the region of bistability.

We discovered an alarming result concerning the overdose death rate for the PRCC results for yearly death counts in the IOUD model with the casual user class. The following applies to both heroin and all-illicit opioids datasets. The significance of the overdose death rate was initially high, as expected. However, its relevance decreased as time moved on, indicating the higher death rate reduced the population in the IOUD class to a degree where fewer individuals interacted with the susceptibles. This decreased interaction led to fewer individuals flowing into the IOUD class. Although we want fewer individuals individuals departing to go to the IOUD class, we do not wish for the reason to be higher overdose deaths. Therefore, there is an urgency to expedite users out of the IOUD class into treatment. This PRCC result concurs with our startling revelation discovered for our original HUD model. For the sake of considering potential scenarios, suppose the growth rate of overdose death rates continues while the other parameters remain at their current estimated values. In that case, the DFE will be the only stable, biologically relevant equilibrium predicted to happen by 2038 in the original HUD model and by 2046 for the HUD model ( = IOUD model with heroin dataset) with casual users. Again, this emphasizes the importance of driving down the future outlook of this epidemic ending with overdose deaths from heroin use. Strategies that could reduce this rate or keep it constant include increased police and law intervention, updated enforcement policies, and unprecedented procedures targeted at enforcement of laws.

Although one would intuitively predict many of our sensitivity analysis results, some interesting results emerged. A surprising result in the sensitivity analysis for the all-illicit opioids dataset was the importance of the casual user transmission rate over-ranking the IOUD transmission rate by far. Thus, it is essential not to overlook the casual users contributing to this epidemic.

Although the last parameter discussed showed significance for some of the variables in the extended model with both datasets, the parameter quantifying the rate of relapse from recovery on their own accord consistently had significance for most of the variables. Furthermore, this case was especially apparent in the opioid dataset and exemplified when comparing the heroin-only and all-illicit opioids dataset parameter values. We note that the parameter for individuals moving to the recovered class through non-specialty treatment facility means or quitting on their own is much higher for the all-illicit opioid epidemic than the heroin epidemic. We hypothesize that relapse may be more problematic for those individuals who go directly to the recovered class instead of going through specialty treatment facilities. Focus on recovered people that are relapsing came there by omega than not through specialty treatment route. Efforts to increase those going into specialty treatment and decrease the relapse rate for the individuals in recovery are exceptionally beneficial for the all-illicit opioid epidemic.

Comparing parameters between the datasets for the IOUD model with a casual user, we discover intriguing results. We note that the rate quantifying how many individuals go back to being susceptible from casual use is much higher for the heroin epidemic than the all-illicit opioids epidemic. Individuals remain longer in the casual user class of the all-illicit opioid epidemic. This development concurs with the PRCC results on the importance of not considering how influential the casual users are in driving the opioid epidemic. We also note that the going to treatment rates were significantly smaller for those in the opioid epidemic, whether casual users or the IOUD class. It signifies that illicit opioid pain medication users are more unlikely to seek treatment. It would be beneficial to raise awareness of this fact.

Future work extends the IOUD model with a casual user class by adding a general treatment class to distinguish between specialty treatment facilities and non-specialty treatment facilities. Furthermore, adding a prescription class to the IOUD model with a casual user class to incorporate those using opioid prescriptions by a physician's order is also another consideration forthcoming. As mentioned earlier, an age-structured model as well as allowing for parameters to change with time could be additional extensions ^[12,13]. Finally, additional data collection could be undertaken (e.g., by SAMHSA) that would give more frequent data points for what is presented in Tables 4–6. Besides likely requiring additional funding, this would necessitate further collaboration with modelers with the hope that the results would give additional insight for best ways to address this epidemic.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors gratefully thank the three reviewers for their important insights in the revisions of this manuscript.

Conflict of interest

The authors declare there is no conflict of interest.