 Research note
 Open Access
 Published:
Role of imitation and limited rehabilitation capacity on the spread of drug abuse
BMC Research Notesvolume 11, Article number: 493 (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.
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.
Several mathematical models describing the spread of psychosocial ills in a community have been proposed, see for example, drug epidemics [2,3,4,5,6,7,8,9], alcoholism [10,11,12,13,14,15,16], smoking [17,18,19]. The basic assumption in most drug abuse models is that there is a direct proportional relationship between the number of drug users in need of treatment and the available health care resources present. In this paper, we develop a mathematical model that takes into account the possibility of the number of drug abusers in need of rehabilitation exceeding the capacity of rehabilitation centers. Recruitment into rehabilitation (inpatient or outpatient) is denoted by H(U) and defined as follows:
where U represents the proportion of individuals abusing drugs, \(\alpha\) is the maximum rehabilitation uptake per unit of time and \(\omega\) measures the extent of the effect of the problem of demand for treatment. Firstly, observe that for small U, \(H(U)\approx \alpha U\), that is, when the number of drug users is not too large, then the rate of entering treatment is proportional to the number of drug users present. Secondly, observe that for large U, \(H(U)\approx \alpha /\omega\), this means that the rate of entering rehabilitation takes a maximum value \(\alpha /\omega\). Finally, when \(\omega =0\), we again obtain the result as in the first case, \(H(U)=\alpha U\), that is, the function returns to a linear function mostly used in previous drug abuse models. Amongst drug abusers who are seeking help through rehabilitation, we have that a proportion p of these individuals are recruited into inpatient rehabilitation and the complementary proportion \((1p)\) are recruited into outpatient rehabilitation. It is also important to note that epidemic models including treatment functions of the form (1) are found in [20,21,22,23].
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,3,4,5,6,7,8,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
The model divides the population into four classes, S(t), U(t), \(R_{op}(t)\) and \(R_{ip}(t)\). The class S(t) represents the population at risk of being initiated into drug abuse. The class U(t) denotes those abusing drugs, \(R_{op}(t)\) denotes those in rehabilitation as outpatients and \(R_{ip}(t)\) denotes those in rehabilitation as inpatients. The total local population is thus given by
The general population enter the susceptible population at a rate \(\Lambda\), that is, the demographic process of individuals reaching age 15 in the modelling time period. Susceptible individuals become drug users upon contact with individuals in compartments U or \(R_{op}\). This results from the assumption that those in inpatient rehabilitation do not have contact with the population at risk. The per capita contact rate \(\beta _1\) is a product of the effective number of contacts \(c_1\), between drug users not in rehabilitation and the susceptible population, and the probability \(\hat{\beta _1}\), that a contact results into initiation into drug use, that is, \(\beta _1=c_1\hat{\beta _1}\). The per capita contact rate \(\beta _2\) is a product of the effective number of contacts \(c_2\), between drug users in outpatient rehabilitation and the susceptible population, and the probability \(\hat{\beta _2}\), that a contact results into initiation into drug use, that is, \(\beta _2=c_2\hat{\beta _2}\). Individuals under outpatient rehabilitation quit drug abuse permanently at a rate \(\delta _1\) and individuals under inpatient rehabilitation quit drug abuse permanently at a rate \(\delta _2\). The general population experience natural death at a rate \(\mu\). Drug users undergoing outpatient rehabilitation relapse into drug use at a rate \(\rho _1\) whereas those undergoing inpatient rehabilitation relapse at a rate \(\rho _2\). The relapse is thus assumed to be a voluntary process, that is not influenced by interaction with users. We allow the transfer from outpatient to inpatient rehabilitation, this happens at a rate \(\gamma _1\). We also allow the transfer from inpatient to outpatient rehabilitation, this rate is represented by \(\gamma _2\). We assume that individuals in each compartment are indistinguishable and there is homogeneous mixing. We have the following general set of nonlinear ordinary differential equations:
with the initial conditions:
where
Here \(\beta _2=\theta \beta _1\), with \(\theta =1\) signifying that the chance of initiating drug abuse habit upon contact with an individual in U or \(R_{op}\) is the same, \(\theta \in (0,1)\) signifying a reduced chance of initiating drug abuse habit upon contact with an individual in \(R_{op}\) as compared to an individual in U, \(\theta >1\) signifies an increased rate of initiating drug abuse habit upon contact with an individual in \(R_{op}\) as compared to an individual in U.
Model analysis
Model properties
Invariant region
It follows from system (2) that
Then, \(\displaystyle \limsup _{t\rightarrow \infty } N\le \dfrac{\Lambda }{\mu }\). Thus, the feasible region for system (2) is
It is easy to verify that the region \(\Omega\) is positively invariant with respect to system (2), see for instance [3,4,5].
The drugfree equilibrium and the abuse reproduction number
Model system (2) always has a drugfree equilibrium \(\mathcal {D}_0=\left( \dfrac{\Lambda }{\mu },\ 0,\ 0,\ 0\right)\). Denote the abuse reproduction number of model system (2) by
with
We can clearly note that \(\gamma _1\gamma _2\le h_1h_2\) and so \((1\Gamma _1)\ge 0\). Also, \(\gamma _1\gamma _2 +\gamma _2\rho _1+\rho _2 h_1\le h_1h_2\) and \(\gamma _1\gamma _2 +\gamma _1\rho _2 +\rho _1 h_2\le h_1h_2\). Therefore, \(\mathcal {R}_a\) is nonnegative. The abuse reproduction number \(\mathcal {R}_a\) of the model, is the average number of secondary cases generated by one drug user during his/her duration of drug use in a population of completely potential drug users.
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 Jacobian matrix of model system Eq. (2) at \(\mathcal {D}_0\) is given by
where \(h_1\) and \(h_2\) are defined as before and \(g_1=\frac{\Lambda }{\mu }\beta _1(\mu +\alpha )\), \(g_2=\frac{\Lambda }{\mu }\beta _2+\rho _1\). The local stability of the drugfree equilibrium is determined by the following submatrix of \(J(\mathcal {D}_0)\),
Since all offdiagonal elements are positive, we now consider matrix \(\bar{J}(\mathcal {D}_0)\). We claim that \(\bar{J}(\mathcal {D}_0)\) is an M—matrix. Multiplying matrix \(\bar{J}(\mathcal {D}_0)\) by the positive \(3\times 1\) matrix
we have
where \(W_2\) is a positive \(3\times 1\) matrix given by
Then, it follows from M—matrix theory that all eigenvalues of \(\bar{J}(\mathcal {D}_0)\) have negative real parts, which implies the local asymptotic stability of the drugfree equilibrium if \(\mathcal {R}_a<1\). On the other hand, it can be shown that the determinant of \(\bar{J}(\mathcal {D}_0)\) is given by
Thus, if \(\mathcal {R}_a<1\), then matrix \(\bar{J}(\mathcal {D}_0)\) has eigenvalues with negative real parts, which implies the stability of the drugfree equilibrium. This completes the proof.\(\square\)
The drugpersistent equilibrium point
The drugpersistent equilibrium \(\mathcal {D}^*=\left( S^*,U^*,R^*_{op},\;R^*_{ip}\right)\) always satisfies
From the last two equations of (5) we have that
where
Substituting expressions (6) into the first equation of (5), we get
Substituting expressions (6) and (7) into the second equation of (5) leads to the following sixth order polynomial equation
Solving (8) gives \(U^*=0\) which corresponds to the drugfree equilibrium or
where the coefficients \(\chi _i,~1\le i\le 5\) are in Additional file 1: Appendix S1. We can clearly note that, \(\chi _0>0\Leftrightarrow \mathcal {R}_a<1\) and \(\chi _0<0\Leftrightarrow \mathcal {R}_a>1\). The number of possible positive real roots of the polynomial (9) can be determined using the Descartes Rule of Signs. The number of positive roots are at most five.
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:
\(S=x_{1},~U=x_2~R_{op}=x_3,~R_{ip}=x_4\), so that \(\text{ N }=\displaystyle \sum _{n=1}^{4}{x_n}\). We now use the vector notation \(X=(x_{1},x_{2},x_{3},x_{4})^{T}\). System (2) can be written in the form \(\dfrac{dX}{dt}=F(t,x(t))=(f_{1},f_{2},f_{3},f_{4})^T\), where
with
Let \(\beta _1\) be the bifurcation parameter, \(\mathcal {R}_a=1\) corresponds to
The Jacobian matrix of system (2) at \(\mathcal {D}_0\) when \(\beta _1=\beta ^*_1\) is given by
where \(h_1\) and \(h_2\) are defined as before and \(g^*_1=\frac{\Lambda }{\mu }\beta ^*_1(\mu +\alpha )\), \(g^*_2=\frac{\Lambda }{\mu }\theta \beta ^*_1+\rho _1\).
System (10), with \(\beta _1=\beta ^*_1\) has a simple eigenvalue, hence the center manifold theory can be used to analyse the dynamics of system (2) near \(\beta _1=\beta ^*_1\). It can be shown that \(J^*(\mathcal {D}_0)\), has a right eigenvector given by \(w=(w_1,w_2,w_3,w_4)^{T}\), where
Further, the left eigenvector of \(J^*(\mathcal {D}_0)\), associated with the zero eigenvalue at \(\beta _1=\beta ^*_1\) is given by \(v=(v_1,v_2,v_3,v_4)^{T}\), where
The computations of a and b are necessary in order to apply Theorem 4.1 in CastilloChavez and Song [24]. For system (10), the associated nonzero partial derivatives of F at the drugfree equilibrium are in Additional file 1: Appendix S2. It thus follows that
where
with
Note that \(\omega ^*>0\), \(\eta ^*_1>0\) and \(\eta ^*_2>0\). Also note that if \(\omega >\omega ^*\), \(\eta _1>\eta ^*_1\) and \(\eta _2>\eta ^*_2\) then \(\text{\bf{a}}>0\) and \(\text{\bf{a}}<0\) if \(\omega <\omega ^*\), \(\eta _1<\eta ^*_1\) and \(\eta _2<\eta ^*_2\). Lastly,
We thus have the following result
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
Since we can rarely enumerate the incidence of drug users, data from treatment centers can be used as proxy for estimating parameters for drug related issues. We use data obtained from previous mathematical models with inpatient and outpatient rehabilitation [4, 5]. Some of the parameter values will be obtained from literature.
Parameter values used for numerical simulations are given in Table 1.
Numerical results
We carry out detailed numerical simulations using matlab to support our theoretical findings. The initial conditions used are: \(S(0)=0.95\), \(U(0)=0.05\), \(R_{op}(0)=0\), \(R_{ip}(0)=0\).
Figures 1 and 2 illustrate the effect of varying parameters \(\omega\) and \(\eta _1\) on the prevalence of drug abuse. Figures 1 and 2 demonstrate that increasing \(\omega\) and \(\eta _1\) results in an increase in the prevalence of drug abuse. This is a reflection that limited rehabilitation and imitation are of major concern in the fight against drug abuse.
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
Like in any model development, the model is not without 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.
References
 1.
United Nations Office on Drugs and Crime. World drug report. New York: United Nations publication; 2015.
 2.
Mulone G, Straughan B. A note on heroin epidemics. Math Biosci. 2009;208:138–41.
 3.
Mushanyu J, Nyabadza F, Muchatibaya G, Stewart AGR. Modelling multiple relapses in drug epidemics. Ricerche di Matematica. 2015. https://doi.org/10.1007/s1158701502410.
 4.
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.
 5.
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.
 6.
Njagarah JBH, Nyabadza F. Modelling the impact of rehabilitation, amelioration and relapse on the prevalence of drug epidemics. J Biol Syst. 2013;21:1350001.
 7.
Nyabadza F, HoveMusekwa SD. From heroin epidemics to methamphetamine epidemics: modelling substance abuse in a South African province. Math Biosci. 2010;225:132–40.
 8.
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.
 9.
White E, Comiskey C. Heroin epidemics, treatment and ODE modelling. Math Biosci. 2007;208:312–24.
 10.
Benedict B. Modeling alcoholism as a contagious disease: how infected drinking buddies spread problem drinking. SIAM News. 2007;40(3):11–3.
 11.
Buonomo B, Lacitignola D. Modeling peer influence effects on the spread of highrisk alcohol consumption behavior. Ricerche di Matematica. 2013;63:101–17.
 12.
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.
 13.
Mulone G, Straughan B. Modelling binge drinking. Int J Biomath. 2012; 5:1250005.
 14.
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.
 15.
Sharma S, Samanta GP. Analysis of a drinking epidemic model. Int J Dyn Control. 2015. https://doi.org/10.1007/s4043501501518.
 16.
Walters CE, Straughan B, Kendal JR. Modelling alcohol problems: total recovery. Ricerche Matematica. 2013;62:33–53.
 17.
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.
 18.
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.
 19.
Sharomi O, Gumel AB. Curtailing smoking dynamics: a mathematical modeling approach. Appl Math Comput. 2008;195:475499.
 20.
Hu Z, Ma W, Ruan S. Analysis of SIR epidemic models with nonlinear incidence rate and treatment. Math Biosci. 2012;238:12–20.
 21.
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.
 22.
Zhang X, Liu XN. Backward bifurcation of an epidemic model with saturated treatment function. J Math Anal Appl. 2008;348:433443.
 23.
Zhou L, Fan M. Dynamics of an SIR epidemic model with limited medical resources revisited. Nonlin Anal Real World Appl. 2012;13:312–24.
 24.
CastilloChavez C, Song B. Dynamical models of tuberclosis and their applications. Math Biosci Eng. 2004;1(2):361–404.
 25.
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.
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.
Author information
Additional file
13104_2018_3574_MOESM1_ESM.pdf
Rights and permissions
Open Access This 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.
About this article
Received
Accepted
Published
DOI
Keywords
 Drug abuse
 Imitation
 Reproduction number
 Rehabilitation capacity