Skip to main content

Role of imitation and limited rehabilitation capacity on the spread of drug abuse



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. Peer-influence 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.


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.


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 psycho-social 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:

$$\begin{aligned} H(U)=\frac{\alpha U}{1+\omega U} \end{aligned}$$

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 \((1-p)\) 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 out-patients and \(R_{ip}(t)\) denotes those in rehabilitation as in-patients. The total local population is thus given by

$$\begin{aligned} N(t)=S(t)+U(t)+R_{op}(t)+R_{ip}(t). \end{aligned}$$

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:

$$\begin{aligned} \left\{ \begin{array}{llll} \dfrac{dS}{dt}&{}=&{}\Lambda -f(S,U,R_{op})-\mu S, \\ \dfrac{dU}{dt}&{}=&{}f(S,U,R_{op})+\rho _1 R_{op}+\rho _2 R_{ip}-\mu U-\dfrac{\alpha U}{1+\omega U},\\ \dfrac{dR_{op}}{dt}&{}=&{}\gamma _2 R_{ip}-(\mu +\gamma _1+\rho _1+\delta _1)R_{op}+\dfrac{(1-p)\alpha U}{1+\omega U},\\ \dfrac{dR_{ip}}{dt}&{}=&{}\gamma _1 R_{op}-(\mu +\gamma _2+\rho _2+\delta _2)R_{ip}+\dfrac{p\alpha U}{1+\omega U}, \end{array} \right. \end{aligned}$$

with the initial conditions:

$$\begin{aligned} S(0)= & {} S_0>0,~U(0)=U_0\ge 0,~R_{op}(0)=R_{op0}\ge 0,~R_{ip}(0)=R_{ip0}\ge 0, \end{aligned}$$


$$\begin{aligned} f(S,U,R_{op})= & \; \beta _1 SU(1+\eta _1 U)+\beta _2 SR_{op} (1+\eta _2 R_{op})\\= & \; \beta _1 \left( SU(1+\eta _1 U)+\theta SR_{op} (1+\eta _2 R_{op})\right) . \end{aligned}$$

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

$$\begin{aligned} \frac{dN}{dt}\le \Lambda -\mu (S+U+R_{op}+R_{ip}). \end{aligned}$$

Then, \(\displaystyle \limsup _{t\rightarrow \infty } N\le \dfrac{\Lambda }{\mu }\). Thus, the feasible region for system (2) is

$$\begin{aligned} \Omega =\left\{ (S,U,R_{op},R_{ip})\in \mathbb {R}^{4}_{+}|~N\le \frac{\Lambda }{\mu }\right\} . \end{aligned}$$

It is easy to verify that the region \(\Omega\) is positively invariant with respect to system (2), see for instance [3,4,5].

The drug-free equilibrium and the abuse reproduction number

Model system (2) always has a drug-free equilibrium \(\mathcal {D}_0=\left( \dfrac{\Lambda }{\mu },\ 0,\ 0,\ 0\right)\). Denote the abuse reproduction number of model system (2) by

$$\begin{aligned} \mathcal {R}_a = & \; \mathcal {R}_{U}+\mathcal {R}_{R_{op}}~~~\text{ where }~~~\\ \mathcal {R}_{U}= & \; \left( \dfrac{\Lambda }{\mu }\right) \left[ \frac{\beta _1 (1-\Phi _1)}{\mu (1-\Phi _1)+\alpha p(1-\Phi _2) +\alpha (1-p)(1-\Phi _3)}\right] ~~~\text{ and }~~~\\ \mathcal {R}_{R_{op}}= & \; \left( \frac{\Lambda }{\mu h_1h_2}\right) \left[ \frac{\beta _2 \left((1-p)\alpha h_2+p\alpha \gamma _2 \right)}{\mu (1-\Phi _1)+\alpha p(1-\Phi _2) +\alpha (1-p)(1-\Phi _3)}\right] \end{aligned}$$


$$\begin{aligned} \Phi _1= & \; \frac{\gamma _1\gamma _2}{h_1h_2},~\Phi _2=\frac{\gamma _1\gamma _2 +\gamma _2\rho _1+\rho _2 h_1}{h_1h_2},~\Phi _3=\frac{\gamma _1\gamma _2 +\gamma _1\rho _2 +\rho _1 h_2}{h_1h_2},\\ h_1= & \; \mu +\gamma _1+\rho _1+\delta _1~~\text{ and }~~h_2=\mu +\gamma _2+\rho _2+\delta _2. \end{aligned}$$

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 non-negative. 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 drug-free steady state

Theorem 1

The drug-free equilibrium \(\mathcal {D}_0\) is locally asymptotically stable when \(\mathcal {R}_a<1\) and is unstable when \(\mathcal {R}_a>1\).


The Jacobian matrix of model system Eq. (2) at \(\mathcal {D}_0\) is given by

$$\begin{aligned} J(\mathcal {D}_0)=\begin{bmatrix} -\mu&-\frac{\Lambda }{\mu }\beta _1&\frac{\Lambda }{\mu }\beta _2&0\\ 0&g_1&g_2&\rho _2\\ 0&(1-p)\alpha&-h_1&\gamma _2\\ 0&p\alpha&\gamma _1&-h_2 \end{bmatrix} \end{aligned}$$

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 drug-free equilibrium is determined by the following submatrix of \(J(\mathcal {D}_0)\),

$$\begin{aligned} \bar{J}(\mathcal {D}_0)=\begin{bmatrix} g_1&g_2&\rho _2\\ (1-p)\alpha&-h_1&\gamma _2\\ p\alpha&\gamma _1&-h_2 \end{bmatrix}. \end{aligned}$$

Since all off-diagonal 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

$$\begin{aligned} W_1=\begin{bmatrix} h_1h_2(1-\Phi _1)\\\ p\alpha \gamma _2+(1-p)\alpha h_2\\ (1-p)\alpha \gamma _1 + p\alpha h_1 \end{bmatrix}, \end{aligned}$$

we have

$$\begin{aligned} -\bar{J}(\mathcal {D}_0)\cdot W_1=(1-\mathcal {R}_a)\cdot W_2 \end{aligned}$$

where \(W_2\) is a positive \(3\times 1\) matrix given by

$$\begin{aligned} W_2=\begin{bmatrix} h_1h_2\left[ \mu (1-\Phi _1)+\alpha p(1-\Phi _2)+\alpha (1-p)(1-\Phi _3)\right] \\ 0\\ 0 \end{bmatrix}. \end{aligned}$$

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 drug-free 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

$$\begin{aligned} \text{ det }~\bar{J}(\mathcal {D}_0)=h_1h_2\left[ \mu (1-\Phi _1)+\alpha p(1-\Phi _2)+\alpha (1-p)(1-\Phi _3)\right] (\mathcal {R}_a-1). \end{aligned}$$

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 drug-free equilibrium. This completes the proof.\(\square\)

The drug-persistent equilibrium point

The drug-persistent equilibrium \(\mathcal {D}^*=\left( S^*,U^*,R^*_{op},\;R^*_{ip}\right)\) always satisfies

$$\begin{aligned} \left\{ \begin{array}{llll} &{}&{}\Lambda -f\left( S^*,U^*,R^*_{op}\right) -\mu S^*=0, \\ &{}&{}f\left( S^*,U^*,R^*_{op}\right) +\rho _1 R^*_{op}+\rho _2 R^*_{ip}-\mu U^*-\dfrac{\alpha U^*}{1+\omega U^*}=0,\\ &{}&{}\gamma _2 R^*_{ip}-(\mu +\gamma _1+\rho _1+\delta _1)R^*_{op}+\dfrac{(1-p)\alpha U^*}{1+ \omega U^*}=0,\\ &{}&{}\gamma _1 R^*_{op}-(\mu +\gamma _2+\rho _2+\delta _2)R^*_{ip}+\dfrac{p\alpha U^*}{1+ \omega U^*}=0. \end{array} \right. \end{aligned}$$

From the last two equations of (5) we have that

$$\begin{aligned} R^*_{op}=\frac{\Psi _1 U^*}{1+\omega U^*}~~\text{ and }~~R^*_{ip}=\frac{\Psi _2 U^*}{1+\omega U^*} \end{aligned}$$


$$\begin{aligned} \Psi _1=\frac{\alpha p \gamma _2+\alpha (1-p)h_2}{h_1h_2\left( 1-\Phi _1\right) }~~\text{ and }~~\Psi _2=\frac{\alpha p h_1+\alpha (1-p)\gamma _1}{h_1h_2\left( 1-\Phi _1\right) }. \end{aligned}$$

Substituting expressions (6) into the first equation of (5), we get

$$\begin{aligned} S^*=\frac{\Lambda \left( 1+\omega U^*\right) ^2}{\left( \mu +\beta _1 U^*(1+\eta _1 U^*)\right) (1+\omega U^*)^2 +\beta _2 \Psi _1 U^*\left( 1+\omega U^* +\eta _2\Psi _1 U^*\right) }. \end{aligned}$$

Substituting expressions (6) and (7) into the second equation of (5) leads to the following sixth order polynomial equation

$$\begin{aligned} U^*\left( \chi _5 U^{*5}+\chi _4 U^{*4}+\chi _3 U^{*3}+\chi _2 U^{*2}+\chi _1 U^* +\chi _0\right) =0. \end{aligned}$$

Solving (8) gives \(U^*=0\) which corresponds to the drug-free equilibrium or

$$\begin{aligned} \chi _5 U^{*5}+\chi _4 U^{*4}+\chi _3 U^{*3}+\chi _2 U^{*2}+\chi _1 U^* +\chi _0=0, \end{aligned}$$

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

$$\begin{aligned} \left\{ \begin{array}{llll} x^{'}_{1}(t)&{}=&{}p\Lambda -h(x_1,x_2,x_3)-\mu x_1=f_1,\\ x^{'}_{2}(t)&{}=&{}h(x_1,x_2,x_3)+\rho _1 x_3 +\rho _2 x_4 -\mu x_2-\dfrac{\alpha x_2}{1+ \omega x_2}=f_{2},\\ x^{'}_{3}(t)&{}=&{}\gamma _2 x_4-h_1 x_3 +\dfrac{\alpha (1-p)x_2}{1+\omega x_2}=f_{3},\\ x^{'}_{4}(t)&{}=&{}\gamma _1 x_3-h_2 x_4 +\dfrac{\alpha p x_2}{1+\omega x_2}=f_{4}, \end{array} \right. \end{aligned}$$


$$\begin{aligned} h(x_1,x_2,x_3)=\beta _1\left( x_1x_2(1+\eta _1 x_2)+\theta x_1x_3(1+\eta _2 x_3)\right) . \end{aligned}$$

Let \(\beta _1\) be the bifurcation parameter, \(\mathcal {R}_a=1\) corresponds to

$$\begin{aligned} \beta _1=\beta ^*_1=\left( \frac{\mu }{\Lambda }\right) \left[ \frac{h_1h_2(\mu (1-\Phi _1)+\alpha p(1-\Phi _2) +\alpha (1-p)(1-\Phi _3))}{h_1h_2(1-\Phi _1)+\alpha p \theta \gamma _2 +\alpha (1-p)\theta h_2}\right] . \end{aligned}$$

The Jacobian matrix of system (2) at \(\mathcal {D}_0\) when \(\beta _1=\beta ^*_1\) is given by

$$\begin{aligned} J^*(\mathcal {D}_0)=\begin{bmatrix} -\mu&-\frac{\Lambda }{\mu }\beta ^*_1&\frac{\Lambda }{\mu }\theta \beta ^*_1&0\\\\ 0&g^*_1&g^*_2&\rho _2\\ 0&(1-p)\alpha&-h_1&\gamma _2\\ 0&p\alpha&\gamma _1&-h_2 \end{bmatrix} \end{aligned}$$

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

$$\begin{aligned} w_1= & \; -h_1h_2 \left(\mu (1-\Phi _1)+\alpha p (1-\Phi _2)+\alpha (1-p)(1-\Phi _3)\right),\\ w_2= & \; \mu h_1h_2(1-\Phi _1),~~w_3=\alpha \mu ((1-p)h_2+p\gamma _2),~~w_4=\alpha \mu (p h_1+(1-p)\gamma _1). \end{aligned}$$

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

$$\begin{aligned} v_1= & \; 0,~~v_2=h_1h_2\left( 1-\Phi _1\right) +\alpha (1-p)\theta h_2 +\alpha p\theta \gamma _2,\\ v_3= & \; h_2 \left( \theta (\alpha +\mu )+ \rho _1\right) +\rho _2 \left( \gamma _1 -\alpha \theta p\right) ,\\ v_4= & \; \rho _2 \left( h_1 +\alpha \theta (1-p)\right) +\gamma _2 \left( \theta (\alpha + \mu )+\rho _1\right) . \end{aligned}$$

The computations of a and b are necessary in order to apply Theorem 4.1 in Castillo-Chavez and Song [24]. For system (10), the associated non-zero partial derivatives of F at the drug-free equilibrium are in Additional file 1: Appendix S2. It thus follows that

$$\begin{aligned} \text{\bf{a}}= & \; v_1w_1w_2\frac{\partial ^2 f_1}{\partial x_1\partial x_2}+v_1w_1w_3\frac{\partial ^2 f_1}{\partial x_1\partial x_3}+v_1w^2_2\frac{\partial ^2 f_1}{\partial x^2_2}+v_1w^2_3\frac{\partial ^2 f_1}{\partial x^2_3}+v_2w_1w_2\frac{\partial ^2 f_2}{\partial x_1\partial x_2}\\&+v_2w_1w_3\frac{\partial ^2 f_2}{\partial x_1\partial x_3}+v_2w^2_2\frac{\partial ^2 f_2}{\partial x^2_2}+v_2w^2_3\frac{\partial ^2 f_2}{\partial x^2_3}+v_3w^2_2\frac{\partial ^2 f_3}{\partial x^2_2}+v_4w^2_2\frac{\partial ^2 f_4}{\partial x^2_2}\\= & \; 2\alpha \omega v_2w^2_2-2(1-p)\alpha \omega v_3 w^2_2 -2\alpha p\omega v_4w^2_2+\beta ^*_1v_2w_1w_2+\theta \beta ^*_1 v_2w_1w_3\\&+\frac{2\Lambda \beta ^*_1\eta _1 v_2w^2_2}{\mu }+\frac{2\theta \Lambda \beta ^*_1\eta _2 v_2 w^2_3}{\mu }\\= & \; \left[ A\omega -\mu ^2h_1h_2(1-\Phi _1)v^2_2\beta ^*_1\right] +\left[ B\eta _1-\mu \alpha p h_1h_2(1-\Phi _2)v^2_2\beta ^*_1\right] \\&+ \left[ C\eta _2 -\mu \alpha (1-p)h_1h_2(1-\Phi _3)v^2_2\beta ^*_1\right] \\= & {} A\left( \omega -\omega ^*\right) +B\left( \eta _1 -\eta ^*_1\right) +C\left( \eta _2 - \eta ^*_2\right) , \end{aligned}$$


$$\begin{aligned} \omega ^*=\frac{\mu ^2h_1h_2(1-\Phi _1)v^2_2\beta ^*_1}{A},~~\eta ^*_1=\frac{\mu \alpha p h_1h_2(1-\Phi _2)v^2_2\beta ^*_1}{B},~~\eta ^*_2=\frac{\mu \alpha (1-p)h_1h_2(1-\Phi _3)v^2_2\beta ^*_1}{C}, \end{aligned}$$


$$\begin{aligned} A= & \; 2\alpha \mu ^2h^2_1h^2_2\left( 1-\Phi _1\right) ^2\times \left[ ((1-\theta )\mu +\delta _1)((1-p)\rho _2+\gamma _2)+(\mu +\delta _2)((1-\theta )\mu +\mu \theta p+\gamma _1+\delta _1 +p\rho _1)\right] ,\\ B= & \; 2\Lambda \mu h^2_1h^2_2\left( 1-\Phi _1\right) ^2v_2\beta ^*_1~~\text{ and }~~C=2\Lambda \mu \theta \alpha ^2((1-p)h_2+p\gamma _2)^2v_2\beta ^*_1. \end{aligned}$$

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,

$$\begin{aligned} \text{\bf{b}}=\Lambda \left( h_2 \left( \alpha \theta (p-1)-h_1\right) +\gamma _2 \left( \gamma _1-\alpha \theta p\right) \right) {}^2>0. \end{aligned}$$

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.

Table 1 Parameter values used in numerical simulations

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\).

Fig. 1

Effects of varying \(\omega\) on the prevalence of drug abuse, starting from 0 up to 1.0 with a step size of 0.5

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.

Fig. 2

Effects of varying \(\eta _1\) on the prevalence of drug abuse, starting from 0 up to 1.0 with a step size of 0.5


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\).


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 self-conversion, drug supply chains etc. may form part of the author’s future research considerations.


  1. 1.

    United Nations Office on Drugs and Crime. World drug report. New York: United Nations publication; 2015.

    Google Scholar 

  2. 2.

    Mulone G, Straughan B. A note on heroin epidemics. Math Biosci. 2009;208:138–41.

    Article  Google Scholar 

  3. 3.

    Mushanyu J, Nyabadza F, Muchatibaya G, Stewart AGR. Modelling multiple relapses in drug epidemics. Ricerche di Matematica. 2015.

    Article  Google Scholar 

  4. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  5. 5.

    Mushanyu J, Nyabadza F, Muchatibaya G, Stewart AGR. Modelling drug abuse epidemics in the presence of limited rehabilitation capacity. Bull Math Biol. 2016.

    Article  Google Scholar 

  6. 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.

    Article  Google Scholar 

  7. 7.

    Nyabadza F, Hove-Musekwa SD. From heroin epidemics to methamphetamine epidemics: modelling substance abuse in a South African province. Math Biosci. 2010;225:132–40.

    Article  PubMed  Google Scholar 

  8. 8.

    Nyabadza F, Njagarah JBH, Smith RJ. Modelling the dynamics of crystal meth (Tik) abuse in the presence of drug-supply shains in South Africa. Bull Math Biol. 2012.

    PubMed  Article  Google Scholar 

  9. 9.

    White E, Comiskey C. Heroin epidemics, treatment and ODE modelling. Math Biosci. 2007;208:312–24.

    Article  PubMed  Google Scholar 

  10. 10.

    Benedict B. Modeling alcoholism as a contagious disease: how infected drinking buddies spread problem drinking. SIAM News. 2007;40(3):11–3.

    Google Scholar 

  11. 11.

    Buonomo B, Lacitignola D. Modeling peer influence effects on the spread of highrisk alcohol consumption behavior. Ricerche di Matematica. 2013;63:101–17.

    Article  Google Scholar 

  12. 12.

    Mubayi A, Greenwood PE, Castillo-Chavez 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.

    Article  Google Scholar 

  13. 13.

    Mulone G, Straughan B. Modelling binge drinking. Int J Biomath. 2012; 5:1250005.

    Article  Google Scholar 

  14. 14.

    Sánchez F, Wang X, Castillo-Chavez 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 evidence-based relapse prevention. New York: Academic Press; 2007.

    Google Scholar 

  15. 15.

    Sharma S, Samanta GP. Analysis of a drinking epidemic model. Int J Dyn Control. 2015.

    Article  Google Scholar 

  16. 16.

    Walters CE, Straughan B, Kendal JR. Modelling alcohol problems: total recovery. Ricerche Matematica. 2013;62:33–53.

    Article  Google Scholar 

  17. 17.

    Alkhudhari Z, Al-Sheikh S, Al-Tuwairqi S. Global dynamics of a mathematical model on smoking. Int Sch Res Not. 2014.

    Article  Google Scholar 

  18. 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.

    Article  Google Scholar 

  19. 19.

    Sharomi O, Gumel AB. Curtailing smoking dynamics: a mathematical modeling approach. Appl Math Comput. 2008;195:475499.

    Google Scholar 

  20. 20.

    Hu Z, Ma W, Ruan S. Analysis of SIR epidemic models with nonlinear incidence rate and treatment. Math Biosci. 2012;238:12–20.

    Article  PubMed  Google Scholar 

  21. 21.

    Zhang J, Jia J, Song X. Analysis of an SEIR epidemic model with saturated incidence and saturated treatment function. Sci World J. 2014.

    Article  Google Scholar 

  22. 22.

    Zhang X, Liu XN. Backward bifurcation of an epidemic model with saturated treatment function. J Math Anal Appl. 2008;348:433443.

    Google Scholar 

  23. 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.

    Article  Google Scholar 

  24. 24.

    Castillo-Chavez C, Song B. Dynamical models of tuberclosis and their applications. Math Biosci Eng. 2004;1(2):361–404.

    Article  PubMed  Google Scholar 

  25. 25.

    Jamison DT, Feachmen RG, Makgoba MW, Bos ER, Baingana FK, Hofman KJ, Rogo KO. Disease and mortality in sub-saharan Africa. 2nd ed. Washington DC: World Bank; 2006.

    Google Scholar 

Download references

Authors' contributions

The author read and approved the final manuscript.


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

Consent to publish

Not applicable.

Ethics approval and consent to participate

No ethical approval was required for this project as this is secondary research.


Not applicable.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information



Corresponding author

Correspondence to Josiah Mushanyu.

Additional file

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, 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 ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Mushanyu, J. Role of imitation and limited rehabilitation capacity on the spread of drug abuse. BMC Res Notes 11, 493 (2018).

Download citation


  • Drug abuse
  • Imitation
  • Reproduction number
  • Rehabilitation capacity