Basic properties Invariant region Because the model monitors changes in the human population, the variables and the parameters are assumed to be positive for all t ≥ 0. [System 1] will therefore be analysed in a suitable feasible region G of biological interest. The following lemma applies to the region that [System 1] is restricted to: Lemma 1 The feasible region G defined by http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3458 with initial conditions http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3468 http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3469 is positively invariant and attracting with respect to [System 1] for all t > 0. Proof: Adding the equations of [System 1] we obtain http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3459 The solution NC(t) of the differential equation, [Eqn 5], has the following property http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3470 where NC(0) represents the sum of the initial values of the variables. http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3471 So if http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3472 thenhttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3473 This means that http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3474 is the upper bound of NC. On the other hand, if http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3475 then NC(t) will decrease to http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3476 as t → ∞. This means that if http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3477 , then the solution http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3478 enters G or approaches it asymptotically. Hence, G is positively invariant under the flow induced by [System 1]. Thus in G, [System 1] is well-posed mathematically. Hence, it is sufficient to study the dynamics of the model in G. Positivity of solutions For [System 1], it is important to prove that all the state variables remain non-negative so that the solutions of [System 1] with positive initial conditions will remain positive for all t > 0. We thus give Lemma 2.Lemma 2 Given that the initial conditions of [System 1] are: S(0) > 0, UL(0) > 0, UH(0) > 0, UT(0) > 0, Q(0) > 0, the solutions S(t), UL(t), UH(t), UT(t) and Q(t) are non-negative for all t > 0.Proof: Assume that http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3481. Thus t ˃ 0 and it follows from the first equation of [System 1] that http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3483 . We thus have http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3484 Hencehttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3485 From the second equation of [System 1], we have http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3486 Similarly, it can be shown that UH(t) > 0, UT(t) > 0 and Q(t) > 0 for all t > 0. Model equilibria and stability analysis Local stability of the drug-free equilibrium [System 1] has a drug-free equilibrium given by: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3488 Following van den Driessche and Watmough17, the linear stability of E0 can be established using the next generation matrix method in [System 1]. Using the notations in van den Driessche and Watmough17 for our system, the matrices for the new infection terms (F) and transition terms (V ) are, respectively, given by: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3489 where b1 = μ + σ + ρ2 and b2 = μ + γ + ψ + δ1. It follows then that the basic reproduction number R0 is given by the spectral radius of FV-1 where V-1 denotes the inverse of V . http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3460 wherehttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3490 R0 in this case represents the average number of secondary cases that one drug user can generate during his or her duration of drug use in a population of potential drug users. The expression of R0 is the sum of two terms representing the contribution of light drug users and hard drug users. Hence, using Theorem 2 of van den Driessche and Watmough17, we establish Theorem 1.Theorem 1 The drug-free equilibrium point, E0 , is locally asymptotically stable if R0 ˂ 1, and unstable if R0 ˂ 1. We now illustrate the above theorem numerically. We performed numerical simulations using a fourth-order Runge–Kutta scheme in Matlab.18 The aim was to verify the analytical results obtained on the stability of [System 1]. We first established the parameter values to be used in the simulations. For the purpose of these simulations, we considered hypothetical populations of one million individuals for the core group and four million individuals for the non-core group. We arbitrarily set the initial conditions for the system for illustrative purposes. We considered the case when R0 ˂ 1, in particular when R0 = 0.6541. For varying initial conditions when R0 = 0.6541, the dynamics of drug users is represented by Figure 2. These results show that the population of drug users declines to zero, that is, it approaches the drug-free equilibrium. The results also show that the drug-free equilibrium is locally asymptotically stable whenever R0 ˂ 1. These results support Theorem 1 on the stability of a drug-free equilibrium.Drug-persistent equilibrium In order to determine the drug-persistent equilibrium of [System 1], we set the equations equal to zero. Let E1 = (S*, UL*, UH*, UT*, Q) represent the drug-persistent equilibrium and let http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3461 be the force of infection at steady state E1. In terms of λ*, the components of E1 are http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3462 wherehttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3492 http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3493 and Np*= Np. By substituting [Eqn 8] into [Eqn 7], and simplifying, it can be shown, after some tedious algebraic manipulations, that the non-zero equilibria of the model satisfy the following quadratic equation in terms of λ*: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3463 wherehttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3500 Thus, the positive drug-persistent equilibria of [System 1] are obtained by solving for λ* from the quadratic equation, [Eqn 9], and substituting the results into the expressions in [Eqn 8]. Clearly, the coefficient a0 of [Eqn 9] is always negative and http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3495 We thus produce Theorem 2 on the existence of the drug-persistent equilibrium. Theorem 2 [System 1] has four cases: 1. a unique drug-persistent equilibrium if R0 ˃ 1 2. a unique drug-persistent equilibrium if b0 ˃ 0 and c0 = 0 or b02 – 4a0c03. two drug-persistent equilibria if b0 ˃ 0 and R0 ˂ 1 4. no drug-persistent equilibrium otherwise. It is clear from Theorem 2 Case 1 that the model has a unique drug-persistent equilibrium whenever R0 ˃ 1. Further, Case 3 suggests the possibility of a backward bifurcation. To check for this, we set the discriminant zero and the result solved for the critical value of R0 , giving http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3464 where R0c is a critical value of R0 , below which no drug-persistent equilibrium exists. (For an effective drug abuse control, the reproduction number should be brought below R0c. The condition R0 ˂ 1 is not sufficient for a complete reversal of the substance abuse epidemic described by [System 1].) Backward bifurcation The phenomenon of backward bifurcation has been observed in many epidemiological models such as models for tuberculosis with exogenous re-infection,19,20,21 vector disease models,22 susceptible-infected-susceptible models with saturation of recoveries,23,24 and in particular, models for drug abuse.13,15 The phenomenon has epidemiological significance whereby the classical requirement of R0 ˂ 1 is, although necessary, no longer sufficient to end the substance abuse epidemic. Theorem 2 can be illustrated in the bifurcation diagram shown as Figure 3. Figure 3 is reminiscent of a standard backward bifurcation diagram (see for instance Dushoff25). We emphasise here that the parameter values chosen are for illustrative purposes only and may not necessarily reflect a real substance abuse phenomenon.The simulation results depicted in Figure 3 show that [System 1] only has the drug-free equilibrium when R0 ˂ R0c ˂ 1, two drug-persistent equilibria when R0 ˂ R0c ˂ 1 and one drug-persistent equilibrium when R0 ˃ 1, as shown by Regions A, B and C, respectively. In Region A, the drug-free equilibrium is locally asymptotically stable, whilst in Region B one of the drug-persistent equilibria is stable and the other is unstable. This result clearly shows the coexistence of two stable equilibria when R0 ˂ R0c ˂ 1, confirming that [System 1] exhibits backward bifurcation. In Region C, the drug-persistent equilibrium is stable. The results shown in Figure 3 are summarised in Table 2. The simulations were in agreement with Theorem 2. The time series plots shown in Figure 4, for different initial conditions, also reflect the existence of multiple steady states. The parameter values are as given for Figure 3 with β values being within each of Regions A, B and C. It can be seen that, irrespective of the initial conditions, the force of infection stabilises to a drug-free equilibrium in Region A, one drug-persistent equilibrium and one drug-free equilibrium in Region B and to a drug-persistent equilibrium in Region C. Lemma 3 is thus established. Lemma 3 [System 1] undergoes backward bifurcation when Case 3 of Theorem 2 holds and R0C < R0 < 1. 3 The number of methamphetamine users (primary and secondary) in South Africa that sought treatment from 1996 to 2009. http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3453

The role of relapse One of the major problems relating to treatment for substance abuse is the relapse of those under treatment into hard drug use. We considered the situation in which there is no relapse to hard drug use for the sake of comparison with the case in which relapse occurs. In this situation we considered the case of r = 0, such that [System 1] reduces to http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3496 [System 2] has the same drug-free equilibrium point as [System 1]. The drug-persistent equilibrium can be obtained by considering quadratic equation, [Eqn 9], when r = 0. The coefficients a0, b0 and c0 in [Eqn 9] reduce to: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3497 In this case, the force of infection at the steady state is λ* = μ[R0 – 1], which is positive when R0 ˃ 1. Then one can show that the drug-persistent equilibriumhttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3498 exists and is unique. S*, UL*, UH*, UT * and Q* are given by: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3499 with Np*= Np and http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3501 Hence, in this case (with r = 0), no drug-persistent equilibrium exists whenever R0 < 1. It follows that, owing to the absence of multiple drug-persistent equilibria for [System 1] with r = 0 and R0 < 1, a backward bifurcation is unlikely for [System 1] with r = 0 and R0 < 1. Figure 5 shows the contribution of the relapse rate r on the prevalence of drug use. In the presence of relapse, the prevalence of drug use is higher. It is important to note that when r = 0, the ability of drug users not in treatment to recruit initiates from the susceptible population is the same as the ability to recruit from individuals in treatment.Global stability of the drug-free equilibrium The absence of multiple drug-persistent equilibria when r = 0, suggests that the drug-free equilibrium of [System 1] is globally asymptotically stable when R0 < 1. We thus produce Theorem 3.Theorem 3 Consider [System 2] with r = 0. The drug-free equilibrium is globally asymptotically stable in G whenever R0 < 1. Proof: Let us consider the following Lyapunov candidate function: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3502 where α1 and α2 are positive constants to be determined. Its time derivative along the trajectories of [System 2] satisfies http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3465 The constants α1 and α2 are chosen such that the coefficient of UH is equal to zero. Thus, one can easily show that α1 = b1 – βSθη and α2 = βη(1 – θ)S + ψ. Because S ≤ S*, after substituting α1 and α2 in [Eqn 11], we obtain V(t) ≤ b1b2(1 – q1)(1 – R0)UL. Thus, V(t) ≤ 0 when R0 ≤ 1. Furthermore V(t) = 0 when R0 = 1, that is, when UL = UH = UT = Q = 0. By LaSalle’s invariance principle, the largest invariant set in G, contained inhttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3503 is reduced to the drug-free equilibrium. This pro ves the global asymptotic stability of E0 in G (see Bhatia and Szegö26, Theorem 3.7.11, page 346).Local stability of the drug-persistent equilibrium We determined the stability of the drug-persistent equilibrium and further investigated the possibility of backward bifurcation as a result of the existence of multiple equilibria as indicated in Theorem 2 Case 3. The stability analysis of the drug-persistent equilibrium point required us to determine the eigenvalues of the Jacobian matrix evaluated at the drug-persistent equilibrium. As expressing drug-persistent equilibria explicitly is complicated for [System 1], the calculation of eigenvalues is mathematically cumbersome. So we used the centre manifold theory as presented by Castillo-Chavez and Song20. To apply this method, we first changed the variables of [System 1] such that S = x1, UL = x2, UH = x3, UT = x4 and Q = x5 with http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3504 [System 1] then becomes http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3505 We choose ϕ = β as the bifurcation parameter. We thus equate R0 = 1 and obtain http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3506 The Jacobian of [System 3] at drug-free equilibrium E0 when ϕ = β , is given by http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3507 We note that the Jacobian J(ϕ) of the linearised system has a simple zero eigenvalue. We can thus use the centre manifold theory to analyse the dynamics of [System 3]. The right eigenvector associated with zero eigenvalue is given by w = (w1, w2, w3, w4, w5)T, where http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3508 http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3509 and (.)T denotes a vector transpose. Further, J(ϕ) has a corresponding left eigenvector v = (v1, v2, v3, v4, v5)T, where http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3510 andhttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3511 We note that all the eigenvectors are positive except for w1 and the value of α is chosen such that v.w = 1. In order to establish the local stability of E1, we used Theorem 4 proven in Castillo-Chavez and Song20 and adopted the use of a and b as in Castillo-Chavez and Song20. In particular, because v1 = v4 = v5 = 0, http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3512 To compute the value of a and b, we first computed the non-zero second-order partial derivatives of [System 3] at drug-free equilibrium such that, http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3513 http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3514 http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3515 and http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3516 It thus follows that http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3517 , where and http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3518 Alsohttp://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3519 Hence the sign of a depends on the value of Γ and X, such that if Γ ˃ X then a ˃ 0 and if Γ ˂ X then a < 0 whilst b ˃ 0. We thus obtain Theorem 4. Theorem 4 If Γ ˃ X, then [System 1] has a backward bifurcation at R0 = 1. Alternatively, if Γ ˂ X, then [System 1] undergoes forward bifurcation and the drug-persistent equilibrium is locally asymptotically stable for R0 ˃ 1 but close to one. Further, using the same initial conditions when R0 = 1.7443, the population of drug users tends to a drug-persistent equilibrium in Figure 6. This pattern indicates that, irrespective of the initial conditions, the population of drug users eventually settles at the drug-persistent equilibrium with increasing time. This result is in agreement with Theorem 4. The role of key parameters It is also important to investigate how some key parameters jointly influence the epidemic. This investigation was performed using contour plots. In Figure 7, contours of R0 are plotted as a function of transition rates σ and ρ2 in Figure 7a, ρ2 and ψ in Figure 7b, γ and ρ2 in Figure 7c and γ and ψ in Figure 7d. Based on the parameter values used in the simulation, Figure 7 shows that increasing σ, ρ2 and γ reduces the model reproduction number, whilst increasing ψ increases R0. This pattern indicates that R0 is a decreasing function of σ, ρ2 and γ, and is an increasing function of ψ. These results can also be obtained by performing a sensitivity analysis on R0. According to the model, to decrease the reproduction number, it is thus necessary to increase the rate at which individuals become hard drug users, the rate at which they permanently quit and the rate at which they are rehabilitated. This result makes sense as increasing forward progression rates eventually leads to more individuals quitting. The significance of increasing σ to fight the epidemic is an outcome of the model formulation for two reasons. Firstly, hard drug users have been assumed to be less effective recruiters and secondly, the class of hard drug users is the entry point into treatment programmes. It is thus advantageous according to the model, for identification purposes, that an individual remains a light drug user for only a short time. In reality, this result remains debatable. Application of the model We applied the model to data on methamphetamine abuse in the Western Cape. [System 1] was fitted to the data for individuals who attended specialist treatment centres in the Western Cape. This data is collected every 6 months by the South Africa Community Epidemiology Network on Drug Use2 for individuals who attend specialist treatment centres in the Western Cape. The data on treatment demand trends was used to model the growth of individuals in the UT class in our model. The data for the growth of methamphetamine users in the Western Cape is given in Table 3. Table 3 includes all individuals who use methamphetamine as their primary and secondary substance of abuse. As the data is collected at 6 monthly intervals, the letter ‘a’ represents the first 6 months of the year (January to June) and ‘b’ represents the second 6 months (July to December). Because of the unavailability of data on transmission and progression rates, we estimated most of the parameters, which makes the setting of initial conditions difficult. Nevertheless, for the purpose of the simulations and illustrating the usefulness of the model, we assumed an initial population of one million for the population of individuals who are prone to become methamphetamine abusers. We set the natural death rate of 0.025.15 Many parameters are known to lie within limits. Only a few parameters are known exactly and it is thus important to estimate the others. The estimation process attempts to find the best accordance between the computed and observed data. The estimation can be carried out by ‘trial and error’ or by the use of software programs that are designed to find parameters that give the best fit. Here, the fitting process involved the use of the least squares curve fitting method. A Matlab18 code was used where unknown parameter values were given a lower and upper bound from which the set of parameter values that produced the best fit were obtained. The parameter values obtained from the fitting are shown in Table 4. Figure 8 is a graphical representation of the model fitted to the data for individuals seeking treatment for methamphetamine abuse. As can be seen in Figure 8, the model fits well with the data. For planning and management of interventions, it is important to project the prevalence of the methamphetamine epidemic over a number of years. In our case, we chose 5 years. The projected prevalence over 5 years to 2015 is shown in Figure 9. The model projection shows that there will be a gradual decrease in prevalence. For the given parameter values, the prevalence declines from a peak value of approximately 2.3×105 to 1.9×105 over a period of 5 years. This estimation, of course, assumes that the dynamics remain the same over the entire period. 2A numerical summary of the backward bifurcation shown in Figure 3 with the corresponding reproduction number (R0) and the local stability of equilibria for each of Regions A, B and C. http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3452

4Parameter values that give the best fit to the data in the model of substance abuse. http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/654/3454

DiscussionWe modified the compartmental deterministic model of Nyabadza and Hove-Musekwa15 to incorporate slow and fast progression of initiates and a cycle of light and hard drug use. We also included individuals who permanently stop using drugs and relapse for those under treatment. Relapse was considered synonymous to re-infection in epidemiological models. Our model was used to gain some insights into the dynamics of substance abuse. We established the local and global stability of the drug-free equilibrium. We noted that the drug-free equilibrium point is locally stable whenever R0 < 1. Also, the model has a unique drug-persistent equilibrium whenever R0 ˃ 1, which shows persistence of substance use in the community. For some specific conditions established in Theorem 2, the model exhibits backward bifurcation and some bifurcation diagrams are presented in Figures 3 and 4. The numerical results suggest that the spread of substance abuse can be controlled through a reduction in the relapse rate ψ, increasing interventions at the light drug users’ phase and increasing the uptake rates into treatment. The existence of backward bifurcation in our model is indicative of complex dynamics. It is not sufficient to reduce R0 below unit to control the substance abuse epidemic but rather the value of R0 should be reduced to below R0c. It was shown that backward bifurcation is caused by relapsing to hard drug use when individuals in treatment are lured back into substance abuse by light drug users. The process remains a subject of debate as individuals in treatment are more likely to revert to drug use without due influence. The model thus suggests that strengthening of treatment programmes to prevent relapse is vital. As with many models, the model presented in this article should be treated with circumspection because of the assumptions made and the difficulty in the estimation of the model parameters. As part of future work to improve the model in this article, the model considered here can be refined to incorporate drug users who start using drugs on their own without having contact with other drug users; the impact of behavioural changes induced by campaigns; age structure; and recruitment by drug lords. The model can also be refined for a specific substance of abuse and be fitted to data. Despite its shortcomings, the model provides useful insights into the possible impact of treatment and relapse in communities struggling with substance abuse. AcknowledgementsA.S.K. appreciates the support of SACEMA in the production of the manuscript. F.N. acknowledges with gratitude the Department of Mathematical Sciences for its support in the production of this manuscript. Competing interests We declare that we have no financial or personal relationships which may have inappropriately influenced us in writing this article. Authors’ contributions A.S.K. was the main author of the manuscript, which is part of an MSc thesis. A.S.K. designed the model, performed the analysis and wrote the manuscript. F.N. supervised the research, and helped with the simulations and revisions. 1.http://dx.doi.org/10.4065/81.1.77164384822.PlüddemannADadaSParryCet al.Monitoring alcohol and substance abuse trends in South Africa3.http://dx.doi.org/10.1007/s10461-008-9367-3183244704.http://dx.doi.org/10.1016/j.drugpo.2007.11.0181820772324352995.ParryCDHSubstance abuse intervention in South Africa6. http://dx.doi.org/10.1016/S0038-0121(03)00029-67.RossiCThe role of dynamic modelling in drug abuse epidemiology8.http://dx.doi.org/10.1016/j.mbs.2009.01.006195637399.De AlarconRThe spread of a heroin abuse in a communityBull Narc196921172210.HuntLGChambersCDThe heroin epidemicsNew York: Spectrum Publications Inc.197611.ttp://dx.doi.org/10.1136/jech.33.4.29912. http://dx.doi.org/10.1016/j.amc.2007.05.01213. http://dx.doi.org/10.1016/j.mbs.2006.10.0081717434614.http://dx.doi.org/10.1016/S0895-7177(98)00095-815.http://dx.doi.org/10.1016/j.mbs.2010.03.0022029870316. http://dx.doi.org/10.1016/0025-5564(94)00066-917. http://dx.doi.org/10.1016/S0025-5564(02)00108-618.Matlab. Version 7.01.Natick, MA: Mathworks200419.http://dx.doi.org/10.1016/j.jtbi.2008.06.0231864438620. Castillo-ChavezCSongBDynamical models of tuberculosis and their applicationsMath Biosci Eng2004136140421.http://dx.doi.org/10.1006/tpbi.2000.14511082821622.http://dx.doi.org/10.1016/j.mbs.2008.05.0021857350723.http://dx.doi.org/10.1016/j.jtbi.2008.05.0151858627724.http://dx.doi.org/10.1016/j.mbs.2007.05.0121770744125.http://dx.doi.org/10.1006/jtbi.1996.0094875952726. BhatiaNPSzegöGPStability theory of dynamical systemsBerlin: Springer-Verlag1970