 Research article
 Open Access
 Published:
Modelling the trends of inpatient and outpatient rehabilitation for methamphetamine in the Western Cape province of South Africa
BMC Research Notesvolume 8, Article number: 797 (2015)
Abstract
Background
Dependence on methamphetamine remains one of the major health and social problem in the Western Cape province of South Africa. We consider a mathematical model that takes into account two forms of rehabilitation, namely; inpatient and outpatient. We examine the trends of these two types of rehabilitation. We also seek to investigate the global dynamics of the developed methamphetamine epidemic model.
Methods
The model is designed by likening the initiation process to an infection that spreads in a community through interactions between methamphetamine users and nonusers. We make use of Lyapunov functions obtained from a suitable combination of common quadratic and Volterratype functions to establish the global stability of the methamphetaminepersistent steady state. The least squares curve fit routine (lsqcurvefit) in Matlab with optimization is used to estimate the parameter values.
Results
The model analysis shows that the model has two equilibria, the methamphetamine free equilibrium and the methamphetamine persistent equilibrium, that are both globally stable when the threshold \(\mathcal {R}_a<1\) and \(\mathcal {R}_a>1\), respectively. Upon fitting the model to data on drug users under rehabilitation, parameter values that give the best fit were obtained. The projections carried out the long term trends of these forms of rehabilitation.
Conclusion
The results suggest that inpatient rehabilitation programs have an increased potential of enhancing the chances of recovery for methamphetamine addicts.
Background
The methamphetamine abuse problem and its consequences to communities in the Western Cape province present a complex scenario that drives sexually transmitted infections morbidity, mortality and heavy budgetary constraints. For example in South Africa over R20 billion is used annually, in treating drug users, cracking down drug traffickers and in prevention and media campaigns [1]. In the Western Cape Province of South Africa, the rising demand for substance abuse services and calls by communities for additional services has led to the Provincial Department of Social Development allocating additional resources to the prevention and treatment of substance use disorders (SUD’s). However, planning around the allocation of these resources has been hampered by several informational barriers. Decision makers within the Western Cape Department of Social Development (DOSD) do not have adequate information on the nature and extent of substance abuse in the province, the extent to which there is unmet substance abuse service needs, where these unmet needs are greatest, and which population subgroups have relatively high unmet needs [2]. There is therefore need for extensive research into methamphetamine abuse trends for the development of public policies and focused correctional services to decrease the gang populations.
Recently, many researchers have drawn a lot of scholarship from infectious diseases modelling with the aim of describing drug abuse spreading like an infectious disease [3–5]. Typical examples are models designed for heroin epidemics [12–15], peer influence [6–9] and the spread of alcoholism [10, 11]. Mathematical models can contribute to the understanding of the various aspects of methamphetamine use and are very useful in determining how prevalent drug use is. They can even help in designing and choosing proper interventions by providing a means of integrating data from different sources, describing a process to increase understanding and simulating policy experiments that are cumbersome in real life [16–18].
There are two primary forms of rehabilitation for methamphetamine addiction namely; inpatient and outpatient rehabilitation [19]. Rehabilitation helps individuals to improve their function, mobility, independence and quality of life. It helps individuals live fully regardless of impairment. It helps people who are living with various health conditions to maintain the functioning they have [20]. Inpatient methamphetamine rehabilitation provides an addict a place to live and at the same time providing 24 h support and treatment. Outpatient methamphetamine rehabilitation allows users to reside at home and come to the treatment center on a regular basis [19] (Figs. 1, 2).
Mathematical models on drug abuse that have been developed so far, have unified rehabilitation, see for instance [14, 15, 21, 31]. In this paper, we develop a mathematical model that takes into account both inpatient and outpatient forms of rehabilitation. Our aim is to examine the trends of these two types of rehabilitation. We also seek to investigate the global dynamics of the developed methamphetamine epidemic model. We analyze the global stability of the model equilibria and fit the model to data on these two forms of rehabilitation.
The paper is arranged as follows; in Sect. “Methods”, we formulate and establish the basic properties of the model. The model is analysed for stability in Sect. “Model analysis”. Parameter estimation and projection graphs are given in Sect. “Results and discussion”. Numerical results are also presented in this section. The paper is concluded in Sect. “Conclusions” (Figs. 3, 4).
Methods
Model formulation
The dynamics of inpatient and outpatient rehabilitation are modelled by considering the human population divided into four distinct compartments. The compartments comprise of S denoting the population at risk of being initiated into methamphetamine abuse, U those initiated into methamphetamine abuse, \(T_0\) those in rehabilitation as outpatients and \(T_i\) those in rehabilitation as inpatients. Here, all susceptible individuals range from age 15–64 years. The total human population is thus given by
We assume that individuals in each compartment are indistinguishable and there is homogeneous mixing. Susceptible humans enter the population through births or immigration at a rate \(\Pi\). Susceptible individuals are initiated into methamphetamine use following interaction with individuals using drugs. We assume an initiation function that is analogous to the force of infection for epidemic models. Thus, the per capita contact rate \(\beta\) is a product of the effective number of contacts c, between methamphetamine users and the susceptible population, and the probability \(\hat{\beta }\), that a contact results into initiation into methamphetamine use, that is, \(\beta =c\hat{\beta }\). A fraction U/N of the contacts is with those drug users not in rehabilitation and the average number of potential contacts resulting in susceptible individuals becoming methamphetamine abusers is \(\beta U/N\). Also, a fraction \(T_0/N\) of the contacts is with individuals under outpatient rehabilitation. The average number of potential contacts of each susceptible individual with individuals in outpatient rehabilitation is \(\beta \eta T_0/N\). The parameter \(\eta\) measures the relative ability for individuals in outpatient rehabilitation to initiate new methamphetamine users. Assuming that the rate at which individuals in outpatient rehabilitation recruit initiates is lower than that for drug users not in rehabilitation, we have, \(0<\eta <1\). This is due to the fact that individuals in outpatient rehabilitation will be receiving some form of counselling and therapy to aid them in quitting drug abuse. Such counselling may lead them to discourage susceptible individuals into becoming methamphetamine abusers. This done by raising awareness on the dangers associated with problematic drug use. Here, we make an assumption that individuals under inpatient rehabilitation cannot produce new initiates due to the fact that, they do not come in contact with the general population during the rehabilitation process. The total number of relevant contacts gives the initiation function,
Upon being initiated into methamphetamine use, a susceptible individual moves into the compartment U, of methamphetamine abusers. Upon realizing the demoralizing consequences related to methamphetamine abuse, individuals in the compartment U seek help via rehabilitation programs. We consider two primary types of rehabilitation programs, namely; inpatient and outpatient forms of rehabilitation. Here, we define the effectiveness of these rehabilitation programs to be, “the measure of the benefits and changes in the functioning of an individual acquired during the period when he/she was under treatment”. The effectiveness of rehabilitation programs is measured by the parameter \(\varepsilon\). If \(\varepsilon = 0\), then the rehabilitation programs are not effective, \(\varepsilon = 1\) corresponds to completely effective rehabilitation programs, while \(0<\varepsilon <1\) implies that rehabilitation programs will be effective to some degree. Amongst individuals in compartment U who are seeking help through rehabilitation, we have that a proportion p of these individuals are recruited into inpatient rehabilitation facilities and the complementary proportion \((1p)\) are recruited into outpatient rehabilitation. The rate at which methamphetamine users are recruited into rehabilitation (inpatient or outpatient) is given by \(\sigma\). Once an individual is in the rehabilitation phase, he/she can either quit permanently, relapse into methamphetamine use or die. Individuals experience natural death at a rate \(\mu\). The rate at which individuals under outpatient rehabilitation quit methamphetamine abuse permanently is represented by \(\delta _0\) and the rate at which those under inpatient rehabilitation quit methamphetamine abuse permanently is represented by \(\delta _1\). Individuals undergoing outpatient rehabilitation relapse into methamphetamine use at a rate \(\rho _1\) and individuals in inpatient rehabilitation facilities relapse into methamphetamine use at a rate \(\rho _2\). Substance abusers undergoing outpatient rehabilitation continue to keep contact with people and circumstances that trigger addiction whereas substance abusers undergoing inpatient rehabilitation are completely immersed in the program and separated from the lifestyle and habits that supported drug use, thus, relapsing may be more likely for outpatient rehabilitants as compared to inpatient rehabilitants. We can safely assume that \(\rho _2<\rho _1\). Individuals under inpatient rehabilitation programs can move to outpatient rehabilitation programs at a rate \(\gamma _2\). This movement may be as a result of some individuals failing to cope up with the higher costs usually associated with inpatient rehabilitation and thereby forcing them rather to seek help through outpatient rehabilitation which is generally cheaper. This movement might also be as a result of noted recovery to some individuals, who if discharged can still fully recover whilst residing at home. The rate at which individuals under outpatient rehabilitation move to inpatient rehabilitation facilities is given by \(\gamma _1\). This movement can be due to the fact that some individuals under outpatient rehabilitation might want to take full advantage of all the resources available in treatment, such as personal therapy, group therapy, educational classes on addiction as well as job skills and other related support structures in order to quicken their recovery process. So this in turn will lead them to seek help from inpatient rehabilitation facilities.
The description of the model formulation and the model diagram lead to the following set of nonlinear ordinary differential equations together with the initial conditions (Figs. 5, 6):
Model analysis
Model properties
Positivity of solutions
We now consider the positivity of model system Eqs. (2)–(5). We prove that all the state variables remain nonnegative and the solutions of model system Eqs. (2)–(5) with positive initial conditions will remain positive for all \(t > 0\). We thus state the following theorem (Figs. 7, 8).
Theorem 1
Given that the initial conditions of model system Eqs. (2)–(5) are \(S(0)>0\), \(U(0)>0\), \(T_0(0)>0\) and \(T_i(0)>0\). There exists \((S(t),U(t),T_0(t),T_i(t)): (0,\infty ) \rightarrow (0,\infty )\) which solve the model system Eqs. (2)–(5).
Proof
Assume that
Thus \(\hat{t}>0\), and it follows from the first equation of model system Eqs. (2)–(5) that
giving
From the second equation of model system Eqs. (2)–(5), we obtain
Similarly it can also be shown that \(T_0(t)>0\) and \(T_i(t)>0\) for all \(t>0\), and this completes the proof. \(\square\)
Invariant region
Theorem 2
The feasible region \(\mathcal {G}\) defined by
with initial conditions \(S_0\ge 0\), \(U_0\ge 0\), \(T_{00}\ge 0\) and \(T_{i0}\ge 0\) is positively invariant and attracting with respect to model system Eqs. (2)–(5) for all \(t>0\).
Proof
Summing up Eqs. (2)–(5), we obtain that the total population satisfies the differential equation
Applying a theorem by Birkhoff and Rota [22], on differential inequalities, we have
where N(0) represents the value of Eqs. (2)–(5) evaluated at the initial values of the respective variables. Taking the limit as \(t\rightarrow \infty\), we have that \(0\le N\le \dfrac{\Pi }{\mu }\). Thus, the state variables remain biologically meaningful in the set
for all positive initial conditions in \(\mathbb {R}^4_+\). Thus \(\mathcal {G}\) is a positively invariant region and all solutions of Eqs. (2)–(5), with \((S_0,U_0,T_{00},T_{i0})\in \mathbb {R}^4_+\) remain in \(\mathcal {G}\) for all \(t>0\). Therefore the \(\omega \text{ limit }\) set of solutions of Eqs. (2)–(5) in \(\mathcal {G}\) are contained in \(\mathcal {G}\). The uniqueness, existence and continuity results hold for Eqs. (2)–(5). The system is thus mathematically and epidemiologically wellposed, see also [4]. Our analysis will be based on the dynamics of the solutions generated in \(\mathcal {G}\). \(\square\)
Methamphetaminefree equilibrium and the abuse reproduction number
In this section, we carry out stability analysis of the methamphetaminefree equilibrium. The model has a methamphetaminefree equilibrium given by
a scenario depicting a methamphetaminefree state in the community or society. The abuse reproduction number of the model, denoted \(\mathcal {R}_a\), 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. The determination of \(\mathcal {R}_a\) through the next generation matrix [23] method has been explored in many papers (see for instance [24–29]). Using the same method we have
where
Here, \(\mathcal {R}_a\) is the sum of three subreproduction numbers representing the contributions of individuals in compartments U, \(T_0\) and \(T_i\) respectively. The term \(\Psi _1\) is the proportion of individuals who move back and forth between compartments \(T_0\) and \(T_i\), \(\Psi _2\) is the proportion of individuals who move back and forth between compartments \(T_0\) and U and \(\Psi _3\) is the proportion of individuals who move back and forth between compartments \(T_i\) and U. Interestingly, \(\Psi _4\) is the proportion of those individuals who were once inpatient rehabilitants and have relapsed after joining outpatient rehabilitation and lastly, \(\Psi _5\) is the proportion of those individuals who were once outpatient rehabilitants and have relapsed after joining inpatient rehabilitation facilities. The terms \(\Psi _4\) and \(\Psi _5\) play an essential role in evaluating the effectiveness of inpatient and outpatient rehabilitation programs.
The following result follows from [23].
Theorem 3
The methamphetaminefree equilibrium point \(\mathcal {G}_0\) of model Eqs. (2)–(5) is locally asymptotically stable for \(\mathcal {R}_a\le 1\) and unstable otherwise.
Global stability of the methamphetaminefree steady state
We shall now prove the global stability of the methamphetaminefree equilibrium point \(\mathcal {G}_0\) whenever the reproduction number is less than unity.
Theorem 4
The methamphetaminefree equilibrium point \(\mathcal {G}_0\) of model Eqs. (2)–(5) is globally asymptotically stable for \(\mathcal {R}_a\le 1\) and unstable if \(\mathcal {R}_a>1\).
Proof
Let \(V(U,T_0,T_i)\,=\,aU+bT_0+cT_i\) be a candidate Lyapunov function for some nonnegative constants a, b and c. The time derivative of V is given by
We now evaluate the coefficients of the suitable Lyapunov function such that the coefficients of \(T_0\) and \(T_i\) are equal to zero. We thus obtain
Using these coefficients, the time derivative of the Lyapunov function can be expressed as
We can deduce that \(\dfrac{dV}{dt}\le 0\) when \(\mathcal {R}_a\le 1\) with equality if \(\mathcal {R}_a\,=\,1\). Furthermore, \(\dfrac{dV}{dt}\,=\,0\) if and only if \(U\,=\,T_0\,=\,T_i\,=\,0\). Therefore, the largest compact invariant set in \((S,U,T_0,T_i)\in \mathcal {G}\) such that \(\dfrac{dV}{dt}\,=\,0\) when \(\mathcal {R}_a\le 1\) is the singleton \(\mathcal {G}_0\). By Lasalle invariance principle [30], this implies that \(\mathcal {G}_0\) is globally stable in \(\mathcal {G}\) if \(\mathcal {R}_a\le 1\). We observe that the Jacobian matrix evaluated at \(\mathcal {G}_0\) has a positive eigenvalue whenever \(\mathcal {R}_a>1\). Therefore the methamphetaminefree equilibrium is unstable if \(\mathcal {R}_a>1\). This completes the proof. \(\square\)
The methamphetamine persistent equilibrium point
In this section we determine the methamphetamine persistent equilibrium point by first setting the left hand side of model Eqs. (2)–(5) to zero. Some easy to follow algebraic manipulations give the methamphetaminepersistent equilibrium \(\mathcal {G}^*\,=\,(S^*,U^*,T^*_0,T^*_i)\) where
We thus have the following result:
Theorem 5
Model system Eqs. (2 3 4)–(5) has a unique methamphetaminepersistent equilibrium whenever \(\mathcal {R}_a>1\).
The results on persistence of the model are given in [31]. For more information about persistence of dynamical systems, we refer the reader to [32–36].
Global stability of the methamphetaminepersistent steady state
In this section, we prove the global stability of the methamphetaminepersistent steady state.
Theorem 6
Assume that \(\gamma _1\,=\,\theta \gamma _2\) where \(\theta \,=\,\dfrac{p\rho _1}{(1p)\rho _2}\). If \(\mathcal {R}_a>1\), then the unique methamphetaminepersistent equilibrium \(\mathcal {G}^*\) of model system Eqs. (2)–(5) is globally asymptotically stable in the interior of \(\mathcal {G}\).
Proof
The global asymptotic stability of the methamphetaminepersistent steady state is proved by constructing a global Lyapunov function following [37]. We propose the following Lyapunov function obtained from a suitable combination of common quadratic and Volterratype functions. Define \(\mathcal {V}:~\{(S,U,T_0,T_i)\in \mathcal {G}:~S,U,T_0,T_i>0\}\rightarrow \mathbb {R}\) by
This function is defined, continuous and positive definite for all \(S,U,T_0,T_i>0\). We also observe that \(\mathcal {V}(S,U,T_0,T_i)\,=\,0\) at the steady state \(\mathcal {G}^*\). Since \(\mathcal {G}^*\,=\,(S^*,U^*,T^*_0,T^*_i)\) is an endemic steady state point of model system Eqs. (2)–(5), we have
Computing the time derivative of \(\mathcal {V}(S,U,T_0,T_i)\) along the solution of model system Eqs. (2)–(5), we obtain
Using model system Eqs. (2)–(5) and (8), it can be easily shown that
Similarly, using model system Eqs. (2)–(5) and (8), we can also show that
We also have
and
Adding expressions (10), (11), (12) and (13) gives
where
We observe from expression (14) that (Figs. 9, 10)
only if \(S\,=\,S^*,~U\,=\,U^*\quad \text{ and }\quad T_0\,=\,T^*_0\). Thus we can deduce that, \(\dot{\mathcal {V}(S,U,T_0,T_i)}\) is negative definite if
and \(\dot{\mathcal {V}(S,U,T_0,T_i)}\,=\,0\) only if \(S\,=\,S^*,~U\,=\,U^*,~T_0\,=\,T^*_0 \quad \text{ and } \quad T_i\,=\,T^*_i\). Therefore, the largest compact invariant set in \(\{(S,U,T_0,T_i)\in \mathcal {G}:\) \(\dot{\mathcal {V}(S,U,T_0,T_i)}\,=\,0\}\) is the singleton \(\{{\mathbf {G}}^{*}\}\). By LaSalle’s Invariance Principle [30], we conclude that \({\mathbf {G}}^*\) is globally asymptotically stable in the interior of \(\mathcal {G}\). This completes the proof. \(\square\)
Results and discussion
Numerical simulations
Data
In this section we present an application of our model through fitting the model data on rehabilitation from the Medical Research Council’s (MRC’s), South African Community Epidemiology Network on Drug Use (SACENDU) project [40]. We fit the model system Eqs. (2)–(5) to data for individuals under inpatient and outpatient rehabilitation in Cape Town. We use data in Table 1 which was collected by the SACENDU for individuals who were admitted on inpatient and outpatient forms of rehabilitation in Cape Town.
While inpatient and outpatient rehabilitation programs may have similar components, they can have very different longterm outcomes and success rates. Statistics on success rate of substance abuse rehabilitation programs vary widely [39] and depend on a variety of factors, such as, the extent and nature of the patient’s problems, the appropriateness of treatment and related services used to address those problems, and the quality of interaction between the patient and his or her treatment providers.
Parameter estimation
We make use of Matlab programming language to estimate model parameters used in our numerical simulations and to analyze existing trends on inpatient and outpatient forms of rehabilitation. The model parameters that we use for numerical simulations are in Table 2. We make use of curve fitting, which is a process that allows us to quantitatively estimate the trend of the outcomes. The curve fitting process fits equations of approximating curves to the raw field data. Nevertheless, for a given set of data, the fitting curves of a given type are generally not unique. Thus, a curve with a minimal deviation from all data points is desired. This bestfitting curve can be obtained by the method of least squares. The least squares curve fit routine (lsqcurvefit) in Matlab with optimization is used to estimate the parameter values. Many parameters are known to lie within some intervals. During the estimation of parameter values, a Matlab code is used where unknown parameter values are given a lower and upper bound from which the set of parameter values that provide the best fit are obtained. The intervals used and a few parameters obtained from literature are given in Table 2.
Conclusions
In this paper, we designed a model that incorporates inpatient and outpatient rehabilitation of methamphetamine addicts to study the dynamics of methamphetamine abuse in Cape Town of South Africa. The reproduction number was derived and qualitatively used to investigate the stability of equilibrium states and the prevalence of methamphetamine abuse. The methamphetaminefree equilibrium point is shown to be globally asymptotically stable whenever the reproduction number is less than unity. Thus, methamphetamine abuse can be eliminated if control strategies are put in place and efforts are directed toward reducing the threshold number to a value less than unity. With the aid of a Lyapunov function obtained from a suitable combination of common quadratic and Volterratype functions, the methamphetaminepersistent steady state has been shown to be globally asymptotically stable whenever the reproduction number is greater than unity.
The least squares curve fit routine (lsqcurvefit) in Matlab with optimization has been used to fit the model to data on inpatient and outpatient rehabilitants with the objective of using the model parameters that give the best fit to obtain the incidence curve. We also used the parameter values that give the best fit to plot a graph of the ratio of inpatients to outpatients and as well predict future proportions of inpatient and outpatient rehabilitants. The results suggest that the proportion of patients admitted under inpatient rehabilitation facilities in Cape Town will continue to decrease for the next five years whereas that for outpatient rehabilitants will increase for the next five years. The estimated proportion of inpatient rehabilitants in specialist treatment centres of Cape Town was observed to be approximately \(31\,\%\) by the year 2018. Our estimated incidence of methamphetamine abuse related to data on inpatient and outpatient rehabilitants was observed to be generally decreasing over the years. However, it was noted that the estimated incidence for methamphetamine abuse related to data on inpatient rehabilitants had a sharp decrease as compared to that of outpatient rehabilitants, suggesting that inpatient rehabilitation programs have an increased potential of positively changing the lives of many methamphetamine addicts. The projections carried out show that the estimated incidence for methamphetamine abuse related to data on inpatient rehabilitants will have decreased down to below \(7\,\%\) by the year 2018 and that of outpatient rehabilitants will have decreased to below \(11.5\,\%\) by the year 2018. The study presented here is not exhaustive, it can be extended to include stratification of the population according to levels of methamphetamine use. Structuring the population in such a way would give some other helpful insights in studying the dynamics of methamphetamine abuse. In addition, the model did not take into account contextual dynamics, such as drug supply chains or changes in interdiction. Incorporating these processes will undoubtedly facilitate in the understanding of methamphetamine dynamics. Also, since the study focuses specifically on methamphetamine abuse, the model dynamics may also be affected by the availability and use of other drugs.
Abbreviations
 SUD:

Substance use disorder
 DOSD:

Department of Social Development
 MRC:

Medical Research Council
 SACENDU:

South African Community Epidemiology Network on Drug Use
 MA:

methamphetamine
References
 1.
The Naked Truth (TNT). 2010. http://www.tnt.org.za.
 2.
Myers B, Louw J, Fakier N. Alcohol and drug abuse: removing structural barriers to treatment for historically disadvantaged communities in Cape Town. Int J Soc Welf. 2008;17:156–65.
 3.
Anderson RM, May RM. Infectious diseases in humans: dynamics and control. Oxford: Oxford University Press; 1991.
 4.
Hethcote HW. The mathematics of infectious diseases’. Soc Ind Appl Maths Rev. 2000;42:599.
 5.
Brauer F, van den Driessche P, Wu J. Mathematical epidemiology. Lecture notes in mathematics. Mathematical biosciences subseries; 2008. pp. 1945.
 6.
Buonomo B, Lacitignola D. Modeling peer influence effects on the spread of high risk alcohol consumption behavior. Ricerche di Matematica. 2014;63:101–17.
 7.
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.
 8.
Manthey JL, Aidoob A, Ward KY. Campus drinking: an epidemiological model. J Biol Dyn. 2008;2:346356.
 9.
Do TS, Lee YS. A differential equation model for the dynamics of youth gambling. Osong Public Health Res Perspect. 2014;5:19.
 10.
Benedict B. Modeling alcoholism as a contagious disease: how infected drinking buddies spread problem drinking. SIAM News. 2007:40.
 11.
Snchez F, Wang X, CastilloChvez C, Gorman DM, Gruenewald PJ. Drinking as an epidemic: a simple mathematical model with recovery and relapse. Therapists guide to evidencebased relapse prevention. New York: Academic Press; 2007.
 12.
De Alarcon R. The spread of a heroin abuse in a community. Bull Narc. 1969;21:17–22.
 13.
Hunt LG, Chambers CD. The heroin epidemics. New York: Spectrum Publications Incorporated; 1976.
 14.
Mackintosh DR, Stewart GT. A mathematical model of a heroin epidemic: implications for control policies. J Epidemiol Commun Health. 1979;33:299–304.
 15.
Nyabadza F, HoveMusekwa SD. From heroin epidemics to methamphetamine epidemics: modelling substance abuse in a South African province. Math Biosci. 2010;225:132–40.
 16.
Behrens DA, Caulkins JP, Tragler G, Haunschmied JL, Feichtinge G. A dynamic model of drug initiation: implications for treatment and drug control. Math Biosci. 1999;159:1–20.
 17.
Behrens DA, Caulkins JP, Tragler G, Feichtinger G. Optimal control of drug epidemics: prevent and treatbut not at the same time? Manag Sci. 2000;46:333–47.
 18.
Billard L, Dayananda PWA. Drug addictionpusher generated from addicts. Biomed J. 1993;35:227–44.
 19.
Cross Roads Recovery Centre, Pretoria, South Africa. 2013. http://www.crossroadsrecovery.co.za.
 20.
Greater Toronto Area (GTA) Rehab Network. http://www.gtarehabnetwork.ca.
 21.
Nyabadza F, Njagarah JBH. Smith RJ? Modelling the dynamics of crystal meth (tik) abuse in the presence of drug supply chains in South Africa. Bull Math Biol. 2012. doi:10.1007/s1153801297905.
 22.
Birkhoff G, Rota G. Ordinary differential equations. Needham Heights Ginn. 1982;39:251–7.
 23.
Driessche P, Watmough J. Reproduction numbers and subthreshold endemic equilibria for the compartmental models of disease transmission. Math Biosci. 1999;180:29–48.
 24.
Kodaira JY, de Souza Passos JR. The basic reproduction number in SI staged progression model: a probabilistic approach. In: Dynamics days South America; International Conference on Chaos and Nonlinear Dynamics; 2010.
 25.
Driessche P, Zou X. Modeling relapse in infectious diseases. Math Biosci. 2007;207:89103.
 26.
Hsier YH, Wang YS. Basic reproduction number for HIV model incorporating commercial sex and behavior change. Bull Math Biol. 2006;68:551–75.
 27.
Mastroberardino T. Mathematical modeling of the HIV/AIDS epidemic in Cuba. AMS Eastern Sectional Meeting University of Maryland/Baltimore County; 2014.
 28.
Feng Z, CastilloChavez C. A model for tuberculosis with exogenous reinfection. Theor Popul Biol. 2000;57:235–47.
 29.
Capistrna MA, Morelesa MA, Larab B. Parameter estimation of some epidemic models. The case of recurrent epidemics caused by respiratory syncytial virus. Bull Math Biol. 2009. doi:10.1007/s1153800994293.
 30.
LaSalle JP. The stability of dynamical systems. Society for industrial and applied mathematics. Philadelphia; 1976.
 31.
Njagarah JBH, Nyabadza F. Modelling the impact of rehabilitation, amelioration and relapse on the prevalence of drug epidemics. J Biol Syst. 2013;21.
 32.
Magal P. Perturbation of a globally stable and uniform persistence. J Dyn Differ Equ. 2009;21:1–20.
 33.
Thieme HR. Uniform persistence and permanence for nonautonomous semiflows in population biology. Math Biosci. 2000;166:173–201.
 34.
Butler G, Waltman P. Persistence in dynamical systems. J Diff Equ. 1986;63:255–63.
 35.
Freedman HI, Ruan S. Uniform persistence functional differential equations. J Differ Equ. 1995;115:173–92.
 36.
Freedman HI, Ruan S, Tang M. Uniform persistence and flows near a closed positive set. J Dyn Differ Equ. 1994;6:583–600.
 37.
Cruz VDL. On the global stability of infectious disease model with relapse. Abst Appl. 2013;9:50–61.
 38.
Pluddemann A, Parry CDH, Cerff P, Bhana A, Sanca PE, Potgeiter H, Gerber W, Mohamed F, Petersen P, Carney T. The South African community epidemiology network on drug use (SACENDU). phase 22 SACENDU Research Brief. 2007;10.
 39.
Crystal meth addiction and holistic drug rehab. http://www.drugrehabadvisor.com.
 40.
The South African Community Epidemiology Network on Drug Use (SACENDU). http://www.mrc.ac.za/adarg/sacendu.html.
 41.
Jamison DT, Feachmen RG, Makgoba MW, Bos ER, Baingana FK, Hofman KJ, Rogo KO. Disease and mortality in subsaharan Africa, 2nd edn. Washington D.C.: World Bank; 2006.
Authors' contributions
JM participated in formulating the model, and carried out the stability analysis of the model steady states, and performed the numerical analysis and drafted the manuscript. FN conceived of the study, and participated in model formulation and helped to draft the manuscript. AGRS participated in the stability analysis of the endemic equilibria and numerical analysis. All authors read and approved the final manuscript.
Acknowledgements
The authors would like to thank an anonymous referee for many valuable comments and helpful suggestions. J. Mushanyu and A. G. R. Stewart authors acknowledge, with thanks, the support of the Department of Mathematics, University of Zimbabwe. The F. Nyabadza author acknowledges with gratitude the support from National Research Foundation and Stellenbosch University for the production of this manuscript.
Authors’ information
This work will be part of JM’s Doctor of Philosophy thesis.
Competing interests
The authors declare that they have no competing interests.
Author information
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
 Methamphetamine
 Reproduction number
 Inpatient rehabilitation
 Outpatient rehabilitation
 Least squares curve fitting