Abstract
Opioid addiction has become a global epidemic and a national health crisis in recent years, with the number of opioid overdose fatalities steadily increasing since the 1990s. In contrast to the dynamics of a typical illicit drug or disease epidemic, opioid addiction has its roots in legal, prescription medication—a fact which greatly increases the exposed population and provides additional drug accessibility for addicts. In this paper, we present a mathematical model for prescription drug addiction and treatment with parameters and validation based on data from the opioid epidemic. Key dynamics considered include addiction through prescription, addiction from illicit sources, and treatment. Through mathematical analysis, we show that no addictionfree equilibrium can exist without stringent control over how opioids are administered and prescribed, in which case we estimate that the epidemic would cease to be selfsustaining. Numerical sensitivity analysis suggests that relatively low states of endemic addiction can be obtained by primarily focusing on medical prevention followed by aggressive treatment of remaining cases—even when the probability of relapse from treatment remains high. Further empirical study focused on understanding the rate of illicit drug dependence versus overdose risk, along with the current and changing rates of opioid prescription and treatment, would shed significant light on optimal control efforts and feasible outcomes for this epidemic and drug epidemics in general.
This is a preview of subscription content, access via your institution.
References
Abdurahman X, Zhang L, Teng Z (2014) Global dynamics of a discretized heroin epidemic model with time delay. Abstr Appl Anal 2014:742385
Anderson RM, May RM (1979) Population biology of infectious diseases: part I. Nature 280:361–367
Bailey GL, Herman DS, Stein MD (2013) Perceived relapse risk and desire for medication assisted treatment among persons seeking inpatient opiate detoxification. J Subst Abuse Treat 45(3):302–305
Battista NA (2009) A comparison of heroin epidemic models. https://arxiv.org/pdf/1510.04581.pdf. Accessed 24 Jan 2017
Bicket MC, Long JJ, Pronovost PJ, Alexander GC, Wu CL (2017) Prescription opioid analgesics commonly unused after surgery: a systematic review. JAMA Surg 152(11):1066–1071
Bin F, Xuezhi L, Maia M, Liming C (2015) Global stability for a heroin model with agedependent susceptibility. J Syst Sci Complex 28:1243–1257
Boscarino JA, Rukstalis M, Hoffman SN, Han JJ, Erlich PM, Gerhard GS, Stewart WF (2010) Risk factors for drug dependence among outpatients on opioid therapy in a large us healthcare system. Addiction 105:1776–1782
Boudreau D, Von Korff M, Rutter CM, Sanders K, Ray GT, Sullivan MD, Campbell CI, Merrill JO, Silverberg MJ, BantaGreen C, Weisner C (2009) Trends in longterm opioid therapy for chronic noncancer pain. Pharmacoepidemiol Drug Saf 18:1166–1175
Campbell CI, Weisner C, LeResche L, Ray GT, Saunders K, Sullivan MD, BantaGreen CJ, Merrill JO, Silverberg MJ, Bondreau D, Satre DD, Korff MV (2010) Age and gender trends in longterm opioid analgesic use for noncancer pain. Am J Public Health 100(12):2541–2547
CastilloChavez C, Song B (2004) Dynamical models of tuberculosis and their applications. Math Biosci Eng 1(2):361–404
Centers for Disease Control and Prevention (CDC) (2014) Vital signs: opioid painkiller prescribing, where you live makes a difference. http://www.cdc.gov/vitalsigns/opioidprescribing/. Accessed 21 Aug 2017
Centers for Disease Control and Prevention (CDC) (2017) Annual surveillance report of drugrelated risks and outcomes—united states, 2017. Technical Report Surveillance Special Report 1, Department of health and Human Services, August 31, 2017
Davis JH (2017) Trump declares opioid crisis a ‘health emergency’ but requests no funds. NY Times. https://www.nytimes.com/2017/10/26/us/politics/trumpopioidcrisis.html. Accessed 27 Oct 2017
Diekmann O, Heesterbeek JAP, Metz JAJ (1990) On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases. J Math Biol 35:503–522
Diekmann O, Heesterbeek JAP, Roberts MG (2010) The construction of nextgeneration matrices for compartmental epidemic models. J R Soc Interface 7(47):873–885. https://doi.org/10.1098/rsif.2009.0386
Gladden RM, Martinez P, Seth P (2016) Fentanyl law enforcement submissions and increases in synthetic opioidinvolved overdose deaths—27 States, 2013–2014. MMWR Morb Mortal Wkly Rep 65:837–843
Gossop Michael, Bradley Brendan, Phillips Grania T (1987) An investigation of withdrawal symptoms shown by opiate addicts during and subsequent to a 21day inpatient methadone detoxification procedure. Addict Behav 12:1–6
Gwira Baumblatt JA, Wiedeman C, Dunn JR, Schaffner W, Paulozzi LJ, Jones TF (2014) Highrisk use by patients prescribed opioids for pain and its role in overdose deaths. JAMA Int Med 174(5):796–801
Han B, Compton WM, Blanco C, Crane E, Lee J, Jones CM (2017) Prescription opioid use, misuse, and use disorders in U.S. adults: 2015 national survey on drug use and health. Ann Intern Med 167(5):293–302
Hedegaard H, Warner M, Mini AM (2017) no. Drug overdose deaths in the united states, 1999–2016. National Center for Health Statistics, NCHS Data Brief no. 294:1–8
Heffernan JM, Smith RJ, Wahl LM (2005) Perspectives on the basic reproductive ratio. J R Soc Interface 2:281–293. https://doi.org/10.1098/rsif.2005.0042
Huang G, Liu A (2013) A note on global stability for a heroin epidemic model with distributed delay. Appl Math Lett 26(7):687–691
Hughes PH, Barker NW, Crawford GA, Jaffe JH (1972) The natural history of a heroin epidemic. Am J Public Health 62(7):995–1001
Hughes A, Williams MR, Lipari RN, Bose J, Copello EAP, Kroutil LA (2016) Prescription drug use and misuse in the United States: results from the 2015 national survey on drug use and health. NSDUH Data Review. Retrieved from http://www.samhsa.gov/data
Jones CM (2013) Heroin use and heroin use risk behaviors among nonmedical users of prescription opioid pain relievers–United States. Drug Alcohol Depend 132(1–2):95–100
Kermack WO, McKendrick AG (1927) A contribution to the mathematical theory of epidemics. Proc R Soc Lond A 115:700–721
Keyes KM, Cerda M, Brady JE, Havens JR, Galea S (2014) Understanding the ruralurban differences in nonmedical prescription opioid use and abuse in the United States. Am J Public Health 104(2):e52–e59
Kochanek KD, Murphy SL, Xu J, Arias E (2017) Mortality in the United States, 2016. National Center for Health Statistics, NCHS Data Brief no. 293:1–8
Lankenau SE, Teti M, Silva K, Jackson Bloom J, Harocopos A, Treese M (2012) Initiation into prescription opioid misuse amongst young injection drug users. Int J Drug Policy 23(1):37–44
Ma M, Liu S, Li J (2017) Bifurcation of a heroin model with nonlinear incidence rate. Nonlinear Dyn 88:555–565
Mandell BF (2016) The fifth vital sign: a complex story of politics and patient care. Clevel Clin J Med 83(6):400–401
Muhuri Pradip K, Gfroerer Joseph C, Christine Davies M (2013) Associations of nonmedical pain reliever use and initiation of heroin use in the United States. Technical report, Center for Behavioral Health Statistics and Quality (CBHSQ): Data Review
Njagarah HJB, Nyabadza F (2013) Modeling the impact of rehabilitation, amelioration and relapse on the prevalence of drug epidemics. J Biol Syst 21(1):1350001
Nyabadza F, HoveMusekwa SD (2010) From heroin epidemics to methamphetamine epidemics: modelling substance abuse in a South African province. Math Biosci 225(2):132–140
O’Donnell Julie K, Matthew Gladden R, Seth Puja (2017) Trends in deaths involving heroin and synthetic opioids excluding methadone, and law enforcement drug product reports, by census region United States, 2006–2015. MMWR Morb Mortal Wkly Rep 66:897–903
Office of the Surgeon General. Facing addiction in America: the surgeon general’s report on alcohol, drugs, and health. Technical report, U.S. Department of Health and Human Services (HHS), Washington, DC. Accessed 23 Nov 2017
Perry S, Heidrich G (1982) Management of pain during debridement: a survey of U.S. burn units. Pain 13(3):267–280
Peterson Alexis B, Matthew Gladden R, Delcher Chris, Spies Erica, GarciaWilliams Amanda, Wang Yanning, Halpin John, Zibbell Jon, Lullo McCarty Carolyn, DeFioreHyrmer Jolene, DiOrio Mary, Goldberger Bruce A (2016) Increases in fentanylrelated overdose deaths—Florida and Ohio, 2013–2015. MMWR Morb Mortal Wkly Rep 65:844–849
Pezalla EJ, Rosen D, Erensen JG, Haddox JD, Mayne TJ (2017) Secular trends in opioid prescribing in the USA. J Pain Res 10:383–387
Porter J, Jick H (1980) Addiction rare in patients treated with narcotics. N Engl J Med 302(2):123
Saltelli A (2002) Making best use of model evaluations to compute sensitivity indices. Comput Phys. Commun 145(2):280–297
Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S (2010) Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Comput Phys Commun 18(2):259–270
Samanta G (2011) Dynamic behaviour for a non autonomous heroin epidemic with time delay. J Appl Math Comput 35:161–178
Schug SA, Zech D, Grond S, Jung H, Meuser T, Stobbe B (1992) A longterm survey of morphine in cancer pain patients. J Pain Symptom Manag 7(5):259–266
Scully RE, Schoenfeld AJ, Jiang W, Lipsitz S, Chaudhary MA, Learn PA, Koehlmoos T, Haider AH, Nguyen LL (2018) Defining optimal length of opioid pain medication prescription after common surgical procedures. JAMA Surg 153(3):37–43
Seth P, Rudd RA, Noonan RK, Haegerich TM (2018) Quantifying the epidemic of prescription opioid overdose deaths. Am J Public Health 108(4):500–502
Shah Anuj, Hayes Corey J, Martin Bradley C (2017) Characteristics of initial prescription episodes and likelihood of longterm opioid use—United States, 2006–2015. Cent Disease Control Prev (CDC): Morb Mortal Week Rep 66(10):265–269
Smyth BP, Barry J, Keenan E, Ducray K (2010) Lapse and relapse following inpatient treatment of opiate dependence. Ir Med J 106(6):176–179
Sobol IM (2001) Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comput Simul 55(1–3):271–280
Stewart G, Mackintosh D (1979) A mathematical model of a heroin epidemic: implications for control policies. J Epidemiol Community Health 33:299–304
Substance Abuse and Mental Health Services Administration, Center for Behavioral Health Statistics and Quality. Treatment episode data set (TEDS): 2004–2014. National Admissions to Substance Abuse Treatment Services. Technical report, BHSIS Series S84, HHS Publication No. (SMA) 164986. Rockville, MD: Substance Abuse and Mental Health Services Administration, 2016
Twombly Eric C, Holtz Kristen D (2008) Teens and the misuse of prescription drugs: evidencebased recommendations to curb a growing societal problem. J Prim Prev 29:503–516
U.S. Census Bureau: International Database. Population, total. https://www.census.gov/. Accessed 08 May 2018
U.S. Food and Drug Administration. FDA analysis of longterm trends in prescription opioid analgesic products: quantity, sales, and price trends. Technical report, U.S. Food and Drug Administration, 03, 2018. https://www.fda.gov/downloads/AboutFDA/ReportsManualsForms/Reports/UCM598899.pdf
van den Driessche P, Watmough J (2002) Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math Biosci 180:29–48
Van Zee A (2009) The promotion and marketing of oxycontin: commercial triumph. Am J Public Health 99(2):221–227
Volkow ND, McLellan AT (2016) Opioid abuse in chronic pain—misconceptions and mitigation strategies. N Engl J Med 374:1253–1263
Vowles KE, McEntee ML, Julnes PS, Frohe T, Ney JP, van der Goes DN (2015) Rates of opioid misuse, abuse, and addiction in chronic pain: a systematic review and data synthesis. Pain 156(4):569–576
Watkins KE, Paddock SM, Hudson TJ, Ounpraseuth S, Schrader AM, Hepner KA, Stein BD (2017) Association between process measures and mortality in individuals with opioid use disorders. Drug Alcohol Depend 177:307–314
Weiss RD, Rao V (2017) The prescription opioid addiction treatment study: what have we learned. Drug Alcohol Dependence 173:S48–S54
White E, Comiskey C (2007) Heroin epidemics, treatment and ode modelling. Math Biosci 208(1):312–324
Acknowledgements
The authors would like to thank Christina Battista, Robert Booth, Namdi Brandon, Kathleen Carroll, Jana Gevertz, Anne Ho, Shanda Kamien, Grace McLaughlin, Gianni Migliaccio, Matthew Mizuhara, and Laura Miller for comments, suggestions, and informative conversations. NAB would like to thank Patricia Clark of RIT, whose mathematical biology course gave the original motivation for this project in 2009. We would also like to thank the anonymous reviewers and the associate editor of Bulletin of Mathematical Biology for their helpful comments.
Author information
Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix A
Appendix A
Here, we present supplemental material to support our findings including additional model analysis and validation, numerical stability analysis, and simulation data. We also provide details for the calculation determining a condition for backward bifurcation, the explicit Jacobian used in our stability analysis, and simulation results illustrating system sensitivity to the prescription addiction rate (\(\gamma \)), treatment success rate (\(\delta \)), and prescription rate (\(\alpha \)). Finally, we explore the relationship between prescription rate (\(\alpha \)) and prescription addiction rate (\(\gamma \)).
A.1 Initial Conditions for Validation
We estimated the initial prescribed population, \(P_0\), based off of the percentage of US population to whom were prescribed opioids at any given week in 2009 (\(2\%\)) (Boudreau et al. 2009). Since there were more prescriptions given in 2009 than 1999 (Shah et al. 2017), we estimated that roughly \(0.40\times 2\%\) of the population were prescribed opioids at any time in 1999; hence, \(P_0=0.008\). Note we estimated the coefficient of 0.40 by using the ratio of total opioids MME sold in 1999 to 2009 (U.S. Food and Drug Administration 2018).
We backed out the initial addicted population from the number of prescription opioid deaths in 1999 (2749) (Hedegaard et al. 2017), and normalized it by the fraction of deaths attributed to addicted persons (\(54.6\%\)) (Gwira Baumblatt et al. 2014) and the predicted number of deaths from our model with the ageadjusted US population in 1999 (\(259\times 10^6\)) (U.S. Census Bureau: International Database 2018), e.g., \(A_0 = \frac{(0.546)(2749)}{(259\times 10^6)(\mu ^*\mu )}=0.00136\). We then assumed \(R_0=0.1A_0\) (Office of the Surgeon General 2016) (fraction of population in treatment), making \(S_0=0.990504\).
A.2 Analysis of the AddictionFree Equilibrium
Here, we derive conditions on the existence of an addictionfree equilibrium (AFE) within the system defined by Eqs. 1–4. To begin, we set each equation to zero and require that \(A=0\). Equation 3 becomes \(0=(\delta +\sigma +\mu )R\), and since \(\mu >0\) as a natural death rate, this implies that \(R=0\) at any AFE (conversely, \(R=0\) requires that either \(A=0\) or \(\zeta =0\), which may apply at the beginning of an epidemic). We are left with the system
\(P^*\ne 0\) since otherwise the only solution is \(S^*=P^*=A^*=R^*=0\) and we require that \(S+P+A+R=1\). Then, \(0 = \gamma +\beta _P S\). Since all our parameters and dependent variables are nonnegative by definition, \(\gamma =\beta _P=0\). In this case, opioids are available only through the presence of current addicts (e.g., on the black market due to illicit demand) and not through currently prescribed users. We can now use our assumption that \(1=S+P+A+R\) to find that
A.3 Calculating the Basic Reproduction Number, \(R_0\)
Assuming that \(\gamma =\beta _P=0\), the necessary and sufficient conditions for the AFE to exist, Eqs. 3 and 4 reduce to
Using the next generation method (Diekmann et al. 1990; van den Driessche and Watmough 2002; Heffernan et al. 2005; Diekmann et al. 2010) with both A and R treated as “infected,” we compute the matrices F and V as
Then, \(R_0\) is given by the spectral radius of \(FV^{1}\),
Prevalence of opioid addicts will rise when \(R_0>1\) and fall when \(R_0<1\).
A.4 Jacobian Analysis and Alternative Relapse Models
Consider an alternative form of the model with the addition of two relapse rates \(\nu _1SP\) and \(\nu _2SA\),
The AFE for this system remains the same (with the same conditions for existence) as in Eq. 5. Calculating the basic reproduction number \(\mathcal {R}_0\) using the next generation method, we arrive at
so the addition of \(\nu _1RP\) contributes to \(\mathcal {R}_0\) in a way similar to \(\sigma \) (but scaled by \(P^*\)), while the addition of \(\nu _2RA\) does not contribute to \(\mathcal {R}_0\). We will now conduct further analysis on this model which, as a direct extension of our model given in Eqs. (1)–(4), will include it as a subcase.
Reducing the system to three equations for S, A, R using \(P=1SAR\) gives us
The Jacobian, J, of this system is
Evaluated at the AFE given by Eq. 5 with \(\gamma =\beta _P=0\), the Jacobian \(J(x_0)\) is
Following CastilloChavez and Song (2004), we now take \(\beta _A\) to be the bifurcation parameter (given the form of \(\mathcal {R}_0\)) and conduct analysis around
to analyze the bifurcation of this system when \(R_0=1\) and determine the bifurcation’s direction (CastilloChavez and Song 2004). First, we define the matrix \(\mathcal {A}\) as in CastilloChavez and Song (2004) but, via a change in coordinates, taking \(x_0\) to be the AFE and the bifurcation parameter to be \(\beta _A\). Writing our system of differential equations (including nonlinear relapse terms) as \(dx/dt=f(x,\beta _A)\), we have
It is easy to check that zero is a simple eigenvalue of \(\mathcal {A}\) and that all other eigenvalues of \(\mathcal {A}\) have negative real parts. \(\mathcal {A}\) has right eigenvector \(\mathbf {x}={(S^*(1+\widetilde{\varGamma }),1,\widetilde{\varGamma })^T}\) and left eigenvector \(\mathbf {y} = {(0,1,1\widetilde{\varLambda })}\) where \(\widetilde{\varGamma }\) is given by
and once again
The first component of \(\mathbf {x}\) is negative, but since \(S^*>0\) the analysis still applies (CastilloChavez and Song 2004). We now let \(f_k\) be the kth component of f and set
The nonzero derivatives are
Now,
To make \(a>0\), we therefore need
If this condition is satisfied, there will be a backward bifurcation at \(R_0=1\). Of course, for the model given in Eqs. (1)–(4) where \(\nu _2\) functionally equals zero, it is not possible for a backward bifurcation to occur.
A.5 AddictionFree Equilibrium Numerical Analysis
To examine the sensitivity of the model’s addictionfree equilibrium (AFE) to its parameters, we first ran simulations to see how the AFE changes when either \(\gamma \) or \(\beta _P\) shifts away from zero. Parameter values were chosen as in Table 1 with \(\epsilon =3\) and \(\zeta =0.25\). Our results show that for our estimated parameters resulting in \(R_0\approx 0.022\), shifting \(\beta _P\) away from zero has little noticeable effect, while shifting \(\gamma \) away from zero strongly moves the equilibrium away from the addictionfree state (see Fig. 7). This suggests that in a nearly addictionfree population, prescriptioninduced addiction remains far more important than securing prescriptions away from nonprescribed users. Note that in the exact case of an AFE, it is always stable when \(\gamma =\beta _P=0\) for a parameter space centered around the other parameters listed in Table 1.
Further analysis of the model parameter space when \(\gamma =\beta _P=0\) was conducted using the Sobol method (Sobol 2001). We chose \(N=800{,}000\) and generated \(N(2D+2)\) parameter sets (where \(D=9\) is the dimension of the parameter space) via Saltelli’s extension of the Sobol sequence (Saltelli 2002; Saltelli et al. 2010) for a total of 16 million samples. We then ran the model to 10,000 years for each set of parameters to arrive at an equilibrium. We subsequently conducted Sobol analysis (Sobol 2001) on the values for S, P, A, and R after the final year. Initial conditions for each simulation were \(S(0)=0.9435\), \(P(0)=0.05\), \(A(0)=0.0062\), and \(R(0)=0.0003\) (Fig. 8).
A.6: Further Numerical Exploration of Parameter Space
In this section, we expand our parameter space exploration for \(\{\epsilon ,\zeta \}\in [0.8,8.0]\times [0.2,2.0]\) by examining parameter sensitivity for each of S, P, A, and R instead of only the addicted class. More specifically, we examine the associated effects of \(\epsilon \) and \(\zeta \) on the predicted populations for 10 years into the future for each of the following cases:

1.
Prescription Addiction Rate (\(\gamma \)),

2.
Treatment Success Rate (\(\delta \)),

3.
Prescription Rate (\(\alpha \)),

4.
Prescription Rate vs. PrescriptionInduced Addiction (\(\alpha \) vs. \(\gamma \)).
Figures 5 and 9 show that as \(\gamma \) increases the addicted population grows. In particular, if \(\gamma \) doubles from its estimated value, there exists \((\epsilon ,\zeta )\) for which \(2\%\) of the population becomes addicted to opioids, which is approximately three times the number of addicts in 2016. Moreover, as \(\gamma \) increases, so does the rehabilitation class. Interestingly, for values of \((\epsilon ,\zeta )\) that make the addicted class roughly \(2\%\) of the population, the rehabilitation class makes up approximately \(1\%\). On the other hand, when the rehabilitation class composes roughly \(1.5\%\) of the population, the addicted class makes up roughly the same percentage. When \(\delta \) increases the rehabilitation class, population decreases near zero. The population of the addicted class decreases toward zero as well, while the populations of the susceptible class and prescribed class appear unaffected (Fig. 10).
Figure 11 shows that if the prescription rate \(\alpha \) is small enough, the entire population almost remains in the susceptible class. However, for certain values of \((\epsilon ,\zeta )\) roughly \(0.5\%\) of the population can still remain in the addicted population. Moreover, for all cases of \(\alpha \) and small \(\zeta \), the rehabilitation class’ population remains near zero for almost all values of \(\epsilon \).
Finally, we explore the relationship between prescriptioninduced addiction (\(\gamma \)) and completing the prescription and heading back into the susceptible class (\(\epsilon \)). Situations in which these two parameters do not add to one could be used to model long or shortterm opioid prescription use. The data are presented in Fig. 12. It is clear that a decrease in \(\epsilon \) corresponds to an increase in the number of addicts as might be expected for more chronic opioid prescription use. For large \(\gamma \), those differences are more subtle, as increasing \(\gamma \) leads to a profound escalation in the addicted population regardless of \(\epsilon \).
Rights and permissions
About this article
Cite this article
Battista, N.A., Pearcy, L.B. & Strickland, W.C. Modeling the Prescription Opioid Epidemic. Bull Math Biol 81, 2258–2289 (2019). https://doi.org/10.1007/s11538019006050
Received:
Accepted:
Published:
Issue Date:
Keywords
 Population biology
 Dynamical systems
 Epidemiology
 Compartmental model
 Mathematical biology
 Prescription drug addiction