 Research note
 Open Access
Role of imitation and limited rehabilitation capacity on the spread of drug abuse
 Josiah Mushanyu^{1}Email author
 Received: 13 February 2018
 Accepted: 6 July 2018
 Published: 18 July 2018
Abstract
Objectives
We formulate a mathematical model for the spread of drug abuse using non linear ordinary differential equations. The model seeks to investigate both peer influence and limited rehabilitation effects on the dynamics of drug abuse. Peerinfluence is modelled through the mechanism of imitation and limited rehabilitation is described using a special treatment function. Center manifold theory is used to show that the model exhibits the phenomenon of backward bifurcation. Matlab has been used to carry out numerical simulations to support theoretical findings.
Results
The model analysis shows that the model has multiple equilibria. It has been shown that the classical \(\mathcal {R}_a\)—threshold is not the key to control drug abuse within a population. In fact drug abuse problems may persist in the population even with subthreshold values of \(\mathcal {R}_a\). This was shown to result, in particular when, \(\omega\), \(\eta _1\) and \(\eta _2\) are high enough such that \(\omega >\omega ^*\), \(\eta _1>\eta ^*_1\) and \(\eta _2>\eta ^*_2\). The results suggest the need for comprehensive and accessible substance abuse treatment services to curtail drug abuse.
Keywords
 Drug abuse
 Imitation
 Reproduction number
 Rehabilitation capacity
Introduction
Drug abuse has increased in recent years and is now an epidemic globally. The magnitude of the world drug problem becomes more apparent when considering that more than 1 out of 10 drug users is a problem drug user and the vast majority of these individuals continue to have no access to treatment [1]. There continues to be a large “treatment gap” for substance abuse problems as many countries have a large shortfall in the provision of services. According to the United Nations Office on Drugs and Crime [1], only one out of every six problem drug users in the world has access to treatment. Generally, the number of patients in need of rehabilitation often exceeds the carrying capacities of drug treatment facilities, especially those funded by the state.
We also include peer influence effects on the spread of drug abuse by assuming that the recruitment process happens through the mechanism of imitation. In this paper, we use the recruitment function given in [11]. Compared to previous drug epidemic models [2–9], a key novelty of our model is the inclusion of both imitation and limited rehabilitation on the dynamics of drug abuse.
The paper is arranged as follows; in “Model formulation” section, we formulate and establish the basic properties of the model. The model is analysed for stability in “Model analysis” section. In “Numerical simulations” section, we carry out some numerical simulations. Parameter estimation is also presented in this section. The paper is concluded in the “Conclusions” section.
Main text
Model formulation
Model analysis
Model properties
Invariant region
The drugfree equilibrium and the abuse reproduction number
Local stability of the drugfree steady state
Theorem 1
The drugfree equilibrium \(\mathcal {D}_0\) is locally asymptotically stable when \(\mathcal {R}_a<1\) and is unstable when \(\mathcal {R}_a>1\).
Proof
The drugpersistent equilibrium point
Backward bifurcation
Conditions for the existence of backward bifurcation follow from Theorem 4.1 proven in [24]. Let us make the following change of variables:
Theorem 2
If \(\omega >\omega ^*\), \(\eta _1>\eta ^*_1\) and \(\eta _2>\eta ^*_2\), then model system (2) has a backward bifurcation at \(\mathcal {R}_a=1\).
Results and discussion
Numerical simulations
Parameter estimation
Parameter values used in numerical simulations
Parameter  Range  Value  Source 

\(\beta _1\)  0.10–0.21  0.105  [7] 
\(\beta _2\)  0–0.10  0.063  [6] 
\(\omega\)  0–1  0.62  [5] 
\(\alpha\)  0–0.05024  0.02827  [4] 
p  0–1  0.352  [4] 
\(\eta _1\)  0–1  0.24  Assumed 
\(\eta _2\)  0–1  0.13  Assumed 
\(\delta _1\)  0.001–1  0.01  [4] 
\(\delta _2\)  0.01–1  0.3142  [4] 
\(\rho _1\)  0–0.054  0.0382  [4] 
\(\rho _2\)  0–0.0235  0.0020  [4] 
\(\gamma _1\)  0–0.06012  0.02961  [4] 
\(\gamma _2\)  0–0.008  0.003  [4] 
\(\Lambda\)  0.028–0.080  0.04  [7] 
\(\mu\)  0.019–0.021  0.020  [25] 
Parameter values used for numerical simulations are given in Table 1.
Numerical results
Conclusions
A mathematical model that incorporates imitation and limited rehabilitation has been formulated using nonlinear ordinary differential equations. It has been shown that the classical \(\mathcal {R}_a\)—threshold is not the key to control drug abuse within a population. In fact drug abuse problems may persist in the population even with subthreshold values of \(\mathcal {R}_a\). This was shown to result, in particular when \(\omega\), \(\eta _1\) and \(\eta _2\) are high enough such that \(\omega >\omega ^*\), \(\eta _1>\eta ^*_1\) and \(\eta _2>\eta ^*_2\). Considerable effort should be directed towards reducing \(\omega\), \(\eta _1\) and \(\eta _2\), this done by increasing the value of \(\omega ^*\), \(\eta ^*_1\) and \(\eta ^*_2\) so as to avoid backward bifurcation. Also, results from the model application show that increasing \(\omega\) and \(\eta _1\) lead to an increase in the prevalence of drug abuse. Thus, communities should have suitable capacity for the treatment of drug abusers and specific health and/or education programs may be employed to reduce the imitation coefficient \(\eta _1\).
Limitations

The model did not take into account contextual dynamics, such as drug supply chains or changes in interdiction.

Also, the study presented here ignored detailed social and economic characteristics.

Other initiation processes, not included in this work, for instance, initiation by selfconversion, drug supply chains etc. may form part of the author’s future research considerations.
Declarations
Authors' contributions
The author read and approved the final manuscript.
Acknowledgments
The author acknowledges, with thanks, the support of the Department of Mathematics, University of Zimbabwe for the production of this manuscript.
Competing interests
The author declares no competing interests.
Availability of data and materials
Estimation of parameters have been stated throughout the body of the paper and included in the reference section. The graphs were produced using the MATLAB software that is available from https://www.mathworks.com/products/matlab.html.
Consent to publish
Not applicable.
Ethics approval and consent to participate
No ethical approval was required for this project as this is secondary research.
Funding
Not applicable.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 United Nations Office on Drugs and Crime. World drug report. New York: United Nations publication; 2015.Google Scholar
 Mulone G, Straughan B. A note on heroin epidemics. Math Biosci. 2009;208:138–41.View ArticleGoogle Scholar
 Mushanyu J, Nyabadza F, Muchatibaya G, Stewart AGR. Modelling multiple relapses in drug epidemics. Ricerche di Matematica. 2015. https://doi.org/10.1007/s1158701502410.View ArticleGoogle Scholar
 Mushanyu J, Nyabadza F, Stewart AGR. Modelling the trends of inpatient and outpatient rehabilitation for methamphetamine in the Western Cape province of South Africa. BMC Res Notes. 2015b;8:797.View ArticlePubMedPubMed CentralGoogle Scholar
 Mushanyu J, Nyabadza F, Muchatibaya G, Stewart AGR. Modelling drug abuse epidemics in the presence of limited rehabilitation capacity. Bull Math Biol. 2016. https://doi.org/10.1007/s1153801602185.View ArticleGoogle Scholar
 Njagarah JBH, Nyabadza F. Modelling the impact of rehabilitation, amelioration and relapse on the prevalence of drug epidemics. J Biol Syst. 2013;21:1350001.View ArticleGoogle Scholar
 Nyabadza F, HoveMusekwa SD. From heroin epidemics to methamphetamine epidemics: modelling substance abuse in a South African province. Math Biosci. 2010;225:132–40.View ArticlePubMedGoogle Scholar
 Nyabadza F, Njagarah JBH, Smith RJ. Modelling the dynamics of crystal meth (Tik) abuse in the presence of drugsupply shains in South Africa. Bull Math Biol. 2012. https://doi.org/10.1007/s1153801297905.PubMedView ArticleGoogle Scholar
 White E, Comiskey C. Heroin epidemics, treatment and ODE modelling. Math Biosci. 2007;208:312–24.View ArticlePubMedGoogle Scholar
 Benedict B. Modeling alcoholism as a contagious disease: how infected drinking buddies spread problem drinking. SIAM News. 2007;40(3):11–3.Google Scholar
 Buonomo B, Lacitignola D. Modeling peer influence effects on the spread of highrisk alcohol consumption behavior. Ricerche di Matematica. 2013;63:101–17.View ArticleGoogle Scholar
 Mubayi A, Greenwood PE, CastilloChavez C, Gruenewald PJ, Gorman DM. The impact of relative residence times on the distribution of heavy drinkers in highly distinct environments. Socio Econ Plan Sci. 2010;44:4556.View ArticleGoogle Scholar
 Mulone G, Straughan B. Modelling binge drinking. Int J Biomath. 2012; 5:1250005.View ArticleGoogle Scholar
 Sánchez F, Wang X, CastilloChavez C, Gorman DM, Gruenewald PJ. Drinking as an epidemic: a simple mathematical model with recovery and relapse. In: Alan Marlatt G, Witkiewitz K, editors. Therapists guide to evidencebased relapse prevention. New York: Academic Press; 2007.Google Scholar
 Sharma S, Samanta GP. Analysis of a drinking epidemic model. Int J Dyn Control. 2015. https://doi.org/10.1007/s4043501501518.View ArticleGoogle Scholar
 Walters CE, Straughan B, Kendal JR. Modelling alcohol problems: total recovery. Ricerche Matematica. 2013;62:33–53.View ArticleGoogle Scholar
 Alkhudhari Z, AlSheikh S, AlTuwairqi S. Global dynamics of a mathematical model on smoking. Int Sch Res Not. 2014. https://doi.org/10.1155/2014/847075.View ArticleGoogle Scholar
 Bissell JJ, Caiado CCS, Goldstein M, Straughan B. Compartmental modelling of social dynamics with generalised peer incidence. Math Models Methods Appl Sci. 2014;24:719–50.View ArticleGoogle Scholar
 Sharomi O, Gumel AB. Curtailing smoking dynamics: a mathematical modeling approach. Appl Math Comput. 2008;195:475499.Google Scholar
 Hu Z, Ma W, Ruan S. Analysis of SIR epidemic models with nonlinear incidence rate and treatment. Math Biosci. 2012;238:12–20.View ArticlePubMedGoogle Scholar
 Zhang J, Jia J, Song X. Analysis of an SEIR epidemic model with saturated incidence and saturated treatment function. Sci World J. 2014. https://doi.org/10.1155/2014/910421.View ArticleGoogle Scholar
 Zhang X, Liu XN. Backward bifurcation of an epidemic model with saturated treatment function. J Math Anal Appl. 2008;348:433443.Google Scholar
 Zhou L, Fan M. Dynamics of an SIR epidemic model with limited medical resources revisited. Nonlin Anal Real World Appl. 2012;13:312–24.View ArticleGoogle Scholar
 CastilloChavez C, Song B. Dynamical models of tuberclosis and their applications. Math Biosci Eng. 2004;1(2):361–404.View ArticlePubMedGoogle Scholar
 Jamison DT, Feachmen RG, Makgoba MW, Bos ER, Baingana FK, Hofman KJ, Rogo KO. Disease and mortality in subsaharan Africa. 2nd ed. Washington DC: World Bank; 2006.Google Scholar