Mathematical modelling and analysis of alcohol-methamphetamine co-abuse in the Western Cape Province of South Africa

Substance/drug abuse poses a significant threat to the health and socio-economic fabric of individuals and nations. The combined abuse of alcohol and the highly addictive methamphetamine has worsened the drug epidemic in South Africa, especially in the Western Cape province. In this paper, a mathematical model is formulated to model the dynamics of alcohol and methamphetamine co-abuse. We prove that the equilibria of the submodels are locally and globally asymptotically stable when the sub-model threshold parameters are less than unity. The basic reproduction number due to co-abuse is shown to be the maximum of the two sub-model reproduction numbers. Sensitivity analysis reveals that the most sensitive parameters in the co-abuse epidemic are the alcohol and methamphetamine recruitment rates β1 and β2 respectively. The prevalence curve is indicative of a persistent drug problem in the region. Hence, the need to promote social programs that raise awareness of the dangers posed by multiple substance abuse, through educational campaigns in learning institutions, social media and health institutions. Transmission control must focus on enhancing the quitting process while promoting support services to drug users during and after treatment to minimize cases of relapse. Subjects: Toxicology; Epidemiology; Advanced Mathematics; Applied Mathematics T. O. Orwa ABOUT THE AUTHORS T. O. Orwa is currently a doctoral fellow pursuing his PhD in Biomathematics at the Institute of Mathematical Sciences, Strathmore University (Kenya). He holds a Masters of Science (MSc) degree in Mathematical Sciences from the African Institute for Mathematical Sciences (AIMS), South Africa and Masters of science (MSc) in Mathematics from University of Stellenbosch in South Africa. He has a range of interest in Mathematical Biology with specific interest in substance abuse modelling and infectious disease modelling. F. Nyabadza is a professor in the Department of Mathematical Sciences of Stellenbosch University, South Africa. He has many research articles and has supervised several students at masters level, PhD level and even at post-docs levels. His area of interest is Mathematical Epidemiology (specifically, infectious disease modelling). PUBLIC INTEREST STATEMENT Substance/drug abuse poses a significant threat to the health and socio-economic fabric of individuals and nations. The habit of taking several substances/drugs only degrades the process of detoxification during treatment programmes. The presented model provides the required insights to understanding co-abuse in substance abuse epidemic. The results are relevant in designing intervention strategies aimed at combating multiple substance abuse in the Western Cape Province of South Africa. Orwa & Nyabadza, Cogent Mathematics & Statistics (2019), 6: 1641175 https://doi.org/10.1080/25742558.2019.1641175 © 2019 The Author(s). This open access article is distributed under a Creative Commons Attribution (CC-BY) 4.0 license. Received: 03 March 2018 Accepted: 02 July 2019 First Published: 10 July 2019 *Corresponding author: T. O. Orwa, Institute of Mathematical Sciences, Strathmore University, P.O Box 59857-00200, Nairobi, Kenya E-mail: torwa@strathmore.edu Reviewing editor: J. Alberto Conejero, Universitat Politecnica de Valencia, Spain Additional information is available at the end of the article

for Substance Abuse Treatment, 2005a). The abuse of more than one substance will certainly increase the cost of treatment and relapses (Center for Substance Abuse Treatment, 2005a).
While there could be numerous reasons for using more than one substance, some of the attractions included the following: desire to increase, to balance and to maintain the effects of the primary substance of choice, to experience greater euphoria, and even longer hours of sleep (Bower, 1985). Research work in (Kirkpatrick, Gunderson, R Frances, Foltin, & Hart, 2012), revealed that in comparison to single drug effects, alcohol-MA combination produces a greater elevation of heartbeat and elation, which is arguably a motivation for its popularity and use in the WC region (Mendelson, Jones, Upton, & Jacob, 1995).
One of the purposes of modelling substance abuse epidemics is to provide a rational basis for policies designed to control the spread of substance/drug abuse. A sundry of research works on drug epidemic models have focussed on the abuse of single substances/drugs, see for instance (Burattini, Massad, & Coutinho, 1998;Mulone & Straughan, 2009;Nyabadza & Hove-Musekwa, 2010;Samanta, 2011;Samanta & Sharma, 2012;Sharma & Samanta, 2015;White & Comiskey, 2007a, 2007b. Studies in (Sharma & Samanta, 2015) reveal the possible existence of backward bifurcation in an alcohol drinking epidemic with treatment control. Moreover, a heroin epidemic model in (White & Comiskey, 2007b) is modified in (Samanta, 2011) to identify parameters that drive the epidemic and hence policy formulations that target optimal treatment and prevention measures against the epidemic. Very little work has been done on the mathematical modelling of multi-abuse of substances/drugs. This could be attributed to a lack of data on multiple substance/ drug use. This paper is drawn from a master's thesis that was done by the first author while at Stellenbosch University in South Africa (Orwa, 2014). Treatment geared towards and campaign against multi-substance abuse remain the only viable form of control against the multisubstance/drug abuse epidemics and its spread in communities. We focus on mathematical modelling of co-use of alcohol and methamphetamine. To the author's best knowledge, coinfection modelling has not been applied to substance/drug abuse and transmission.
The rest of the paper is organized as follows: The co-abuse model of alcohol and methamphetamine is formulated in Sections 2. Mathematical analysis of the co-abuse model is performed in Section 3. In Sections 4 and 5, the individual submodels are formulated and analysed. Numerical simulation including parameter estimation, sensitivity analysis and data fitting is performed and discussed in Section 6. The paper is concluded in Section 7 with relevant discussions and recommendations on the co-abuse of alcohol and methamphetamine.

Co-abuse model formulation
We formulate a simple mathematical model that captures the co-use of alcohol and methamphetamine in the population. The model subdivides the human population NðtÞ into seven mutually exclusive compartments, namely those individuals that have never taken alcohol and/ or used methamphetamine, but are at risk of using either substance S (they are also termed susceptibles); alcohol users not under treatment U a ; alcohol users under treatment R a ; methamphetamine users who are not under treatment U t ; methamphetamine users under treatment R t ; users of both alcohol and methamphetamine who are not under treatment U at and users of both alcohol and methamphetamine under treatment R at : Therefore, the total human population at any time t is given by NðtÞ ¼ U a ðtÞ þ R a ðtÞ þ U t ðtÞ þ R t ðtÞ þ U at ðtÞ þ R at ðtÞ: (1)

Basic assumptions
• Individuals get recruited into the susceptible population through birth and immigration at a constant rate, Λ.
• Susceptibles become alcohol users through interactive activities with drinking friends/colleagues and drinking family members.
• Individuals may become methamphetamine users through contact with other methamphetamine users who may include friends, family members and even colleagues in work places.
• Entry into substance abuse is mainly driven by peer pressure. We thus assume that substance abuse spreads like a disease that spreads through contact or proximity to pathogen carriers.
• Alcohol users begin to use methamphetamine through interactive activities with methamphetamine users in places of substance abuse such as clubs and bars.
• Mathemaphetamine users begin to take alcohol through contact with alcohol users.
The class of substance abusers is often divided into light/moderate users and heavy users, see for instance (Nyabadza & Hove-Musekwa, 2010;Nyabadza, Njagarah, & Smith, 2013). For simplicity, we combine these two classes because our primary interest is in those who abuse both alcohol and methamphetamine. Assuming homogeneous mixing of populations, non-alcohol drinkers acquire alcohol drinking habits at rate λ 1 with with the parameter β 1 denoting the effective contact rate (i.e. the contact with an alcohol drinker that will result in one taking alcohol). To account for the decreased chances of becoming an alcohol drinker after being in contact with alcoholics and co-users in rehabilitation, it is assumed here that those under treatment tend to have lower recruitment aptness relative to addicts. Those under the influence of alcohol are more likely to be co-users when they get into contact with those who abuse both substances and not under rehabilitation. Therefore, ζ 1 ; ζ 3 < 1; while ζ 2 > 1: Similarly, individuals get initiated into using methamphetamine at the rate given by where the parameter β 2 is the effective contact rate for initiation into methamphetamine abuse. The modification parameters 1 and 3 are both assumed to be less than unity. This would account for the reduced tendency to initiate new users into methamphetamine abuse by methamphetamine users under treatment. Like in the case of alcohol, we assume that 2 > 1.
Individuals in the compartment U a begin to use methamphetamine at a rate η a λ 2 and move into the compartment U at with η t ! 1 accounting for the increased chances of using methamphetamine by alcohol users when compared to non-alcohol users. Individuals not under treatment may die due to alcohol-related causes at a rate δ 1 or die naturally at a rate μ; a rate that is also assumed for the rest of the population. Furthermore, those that abuse alcohol may seek treatment for alcohol dependence at a rate σ 1 or altogether, quit alcohol drinking at a rate ρ 3 : Alcohol users, upon successful treatment in the compartment R a , permanently quit alcohol at the rate ρ 4 : We assume that those that quit do so permanently. Ideally, quitters should be allowed to relapse. In this model, we only assume relapse in the form of rehabilitation failure. This is a plausible assumption in the Western Cape scenario where drop-out rates of patients at facilities may be as high as 40% (Ramlagan, Peltzer, & Matseke, 2010).
Similarly, susceptible individuals in S, get recruited into the class U t at a rate λ 2 ðtÞ through contact with those using methamphetamine. Some individuals in U t acquire alcohol drinking habits at the rate η t λ 1 and move into the class U at with η t ! 1 accounting for the increased chances of drinking alcohol for methamphetamine users when compared to those not using methamphetamine. Others may enter into rehabilitation in compartment R t at a rate σ 3 and may relapse back to U t at a rate γ 3 : Individuals in the class U t may die due to methamphetamine-related causes at a rate δ 3 or permanently quit methamphetamine at a rate ρ 7 : Users of both alcohol and methamphetamine in the class, U at may revert to alcohol abuse only at a rate ρ 1 or to methamphetamine abuse only at a rate ρ 2 upon quitting either substances of abuse. In addition, individuals in this class may die as a result of death related to abuse of both alcohol and methamphetamine such as liver cirrhosis and pancreatitis at a rate δ 2 : They progress into the treatment class R at at a rate σ 2 : Individuals under treatment for alcohol and methamphetamine abuse in R at ; may relapse into alcohol-methamphetamine abuse compartment at a rate γ 2 or quit drinking of alcohol and relapse to U t , at a rate γ 5 . Similarly, they may quit methamphetamine and revert to the class of alcohol users not under treatment U a at a rate γ 4 . We also allow permanent quitting for individuals in the classes R at and R t at rates ρ 5 and ρ 6 respectively. We assume here that permanent quitting results from effective rehabilitation programs. Clearly, the presented model has additional assumptions. Individuals under treatment take alcohol and/or methamphetamine as the treatment process is taken to be outpatient and hence can initiate others. Although quitting both substances at the same time is unlikely, in this model, we allow permanent quitting of individuals in rehabilitation for both alcohol and methamphetamine.
The flow of individuals between classes is given in Figure 1.
A description of parameters used in the co-abuse model and its sub-models are defined in Table 1.

Model equations
From the above assumptions, parameter definitions and variables, we have the following model system of differential equations, with non-negative initial conditions that describe the dynamics of the co-use alcohol and methamphetamine epidemic. The compartment Q is considered to be superfluous.
3. Mathematical analysis of the co-abuse model

Positivity and boundedness of solutions
For every dynamical system, it is important to establish the long-term behaviour of its solutions. The formulated model system model monitors changes in the human population. We, therefore, assume that the variables and the parameters used are all non-negative for all time, t ! 0.
Theorem 3.1. If Sð0Þ, U a ð0Þ, R a ð0Þ, U t ð0Þ, R t ð0Þ, U at ð0Þ, R at ð0Þ, are non-negative, then so are SðtÞ, U a ðtÞ, R a ðtÞ, U t ðtÞ, R t ðtÞ, U at ðtÞ, R at ðtÞ for all time t > 0. Moreover, Proof. Assume that T is the maximum time for the epidemic. That is, assume, Therefore, T > 0 and from the first equation of model system (4), we obtain So that Hence, SðTÞ ! 0 for "T > 0: Now, from Equation (4), we have; It can similarly be shown that U t ! 0, R t ! 0, U at ! 0 and R at ! 0. Thus the solutions of (4) remain positive for all t ! 0.
The evolution change in the population is given by Quitting rate for methamphetamine abusers without treatment η a , η t , 1 , 2 Modification parameters ζ 1 , ζ 2 , ζ 3 , 3 Modification parameters Λ À μN: It can easily be seen that From Equation (8), we observe that as t ! 1, On the other hand, if N 0 > Λ μ , then NðtÞ will decrease to Λ μ as t ! 1. This means that if N 0 > Λ μ , then the solution ðSðtÞ; U a ðtÞ; R a ðtÞ; U t ðtÞ; R t ðtÞ; U at ðtÞ; R at ðtÞÞ enters the feasible region D or approaches it asymptotically. □ The feasible region, D for system (4) is therefore given as The region D is positively invariant.

Steady states and stability analysis of the co-abuse model
The substance abuse free equilibrium E 0 at of the co-infection model (4) exists when there are no substance users in the community. That is, Following (van Den Driessche & Watmough, 2002), we have the following results on the local stability of E 0 at .
Theorem 3.2. The substance abuse free equilibrium, E 0 at , is locally asymptotically stable if R at < 1 and unstable if R at > 1.
Next, we calculate the co-abuse endemic steady state(s), E Ã at . This is a scenario in which there are some substance users in the society. Upon equating to zero and solving the derived system (in terms of the forces of infections λ 1 and λ 2 ) from the model system (4), we obtain Given that we have, where Note that if λ 1 ¼ 0, then clearly λ 2 ¼ 0. This gives the co-use-free equilibrium derived in (10). The solutions to the remaining part of the polynomial (12), described by Equation (13) define the possible endemic states of the model system (4). A scenario in which the combined abuse of alcohol and methamphetamine persists in the society. We thus have The existence of the endemic stable states for the co-use epidemic model depends on the solutions of (13). The solutions must, however, be real and positive. Determining the explicit existence of the endemic equilibrium in clearly mathematically intractable. We thus comprehensively analyse the sub-models of the co-abuse model and resort to numerical simulations for further applications of the co-abuse model.

Basic reproduction number of the co-abuse model
The basic reproduction number of the co-abuse model system (4), denoted by R at , is calculated using the next generation matrix method (van Den Driessche & Watmough, 2002). It is given by where The term β 1 ½b 2 þζ 1 σ 1 b 1 b 2 ½1ÀΦ 1 represent contribution to new cases of co-abuse of alcohol and methamphetamine resulting from an existing alcohol user. On the other hand, the term β 2 ½b 4 þ 1 σ 3 b 3 b 4 ½1ÀΦ 2 represents the contribution to new "infections" resulting from one user of methamphetamine in a completely susceptible population.

Alcohol abuse model
This model is given The co-abuse model (4) reduces to the following system of equations: where The alcohol model (15) has an alcohol-free equilibrium given by E a 0 ¼ Λ μ ; 0; 0 .

Reproduction number for the alcohol epidemic model (R a0 )
We establish a threshold parameter resulting from the alcohol model. We are interested in establishing the possible number of new cases that would be produced by a single alcohol user, assuming the absence of the other drugs or substances of abuse, in a completely susceptible population. Using the approach and notations described in (van Den Driessche & Watmough, 2002), the matrices for the new infections terms and the transfer terms denoted by F and V respectively at the alcohol-free equilibrium are given as follows: It, therefore, follows that basic reproduction number due to alcohol abuse is given by We now present the stability of the equilibrium points established from the alcohol abuse model (15).

Global stability of the alcohol-free equilibrium
be a candidate Lyapunov function. The derivative of V with respect to time ðtÞ is given by Noting that all the model parameters are positive, it follows that Since Ω is invariant and attracting, it follows that the largest possible invariant set in fðS; U a ; R a Þ 2 Ω : _ V ¼ 0g is the singleton fE a 0 g. Therefore, by the La-Salle's invariance principle (La Salle, 1976), every solution to the equation in the alcohol-only model in system (15) with initial conditions in Ω approaches E a 0 as time approaches infinity. That is, as t ! 1, ðU a ðtÞ; R a ðtÞÞ ! ð0; 0Þ. Substituting for U a ¼ R a ¼ 0 into the model system (15), gives S ! Λ μ as t ! 1. Thus ðSðtÞ; U a ðtÞ; R a ðtÞÞ ! Λ μ ; 0; 0 as t ! 1 for R ao < 1 so that E a 0 is globally asymptotically stable in Ω if R a0 < 1. □

Endemic equilibrium and its stability
Upon equating model system (15) to zero and solving for S Ã , U Ã a and R Ã a in terms of the force of infection λ Ã 1 , we obtain the following expressions.
Furthermore, on substituting equations in (21) into (16), we obtain the polynomial From the polynomial (22), λ Ã 1 ¼ 0 corresponds to the alcohol free-equilibrium. A scenario in which there is no alcohol abuse in the community.
n o which exists for R a0 > 1 corresponds to the endemic equilibrium, a state in which alcohol abuse persists in the community. Therefore, the alcoholendemic equilibrium point(s) is given by E a 1 ¼ ðS Ã ; U Ã a ; R Ã a Þ where S Ã ; U Ã a and R Ã a are as described in Equation (21) and exists only if R a0 > 1.

Global stability of the alcohol-endemic equilibrium point, (Anthuvan)
We show that if the population is constant, the endemic steady state E a 1 , is globally asymptotically stable. We thus assume that The differential equations in system (15) Since N ¼ S þ U a þ R a , is a constant, we let S, U a and R a : We have a non-dimensionalized system given by We argue here that system (25) has a unique endemic equilibrium following the analysis of the original model. We thus have the following theorem.
Theorem 4.2. The alcohol-endemic equilibrium E a 1 is globally asymptotically stable whenever R a0 is greater than unity, and the if the population is constant within the modelling time.
Proof. We propose a suitable Lyapunov function V such that The positive constants A and B are to be determined. We observe that from the proposed Lyapunov function in Equation (26), the first partial derivatives with respect to any of the state variables is given by are all zero at the corresponding alcohol-endemic steady state. That is, at the endemic steady states, s ¼ s Ã , v ¼ v Ã and w ¼ w Ã . In addition, the second partial derivatives of V with respect to any of the three state variables are given as follows; We observe from Equation (28) that all the second partial derivatives are positive. This indicates that the alcohol endemic equilibrium is the minimum of each of the state variables.
The time derivative of the Lyapunov function in (26) is given by Substituting for _ s, _ v and _ w from system (25), we have We now use the system of Equation (25) at E a 1 to obtain: We now substitute the terms in Equation (31) into Equation (30) so that Setting the coefficients of K, JK, JL, K J to zero, we obtain Upon substituting the expressions in Equation (34) back into (30) we obtain From Equation (36), we observe that the expression ðK À JÞ 1 À 1 J À Á is less than or equal to zero with the equality holding if and only if J ¼ 1 or K ¼ J.
Also, the expression ðL À KÞ 1 À 1 K À Á 0 with the equality holding if and only if K ¼ 1 or L ¼ K.
We can draw similar conclusions from the remaining expressions 2 þ K À In the two cases, equality holds if and only if K ¼ J ¼ L.
Therefore, _ V 0 with equality holding if and only if Salle, 1976), we, therefore, conclude that the endemic equilibrium E a 1 is globally asymptotically stable in the interior of Ω. This shows that every solution in Ω or that intersects Ω, would approach the endemic equilibrium E a 1 . □

Methamphetamine-abuse model
The co-use model reduces to the following model when only methamphetamine is being abused.

Model equilibria
Setting the right-hand side of equations in model system (37) to zero, and solving for S, U t and R t , we obtain the polynomial From Equation (38), λ ÃÃ 2 ¼ 0 corresponds to the methamphetamine-free equilibrium, denoted by E t 0 and this is a scenario in which the community is free of methamphetamine abuse.
Note that R t0 is the methamphetamine-abuse reproduction number.

Stability analysis
From polynomial equation given in (38), which exists only if R t0 > 1, and corresponds to the methamphetamine-endemic equilibrium, Due to the similarity between the structures of the two sub-models, we state the following theorems on the stability of the steady states.
Theorem 5.1. The unique methamphetamine-free equilibrium point, E t 0 is locally asymptotically stable for R t0 < 1 and unstable for R t0 > 1: Theorem 5.2. The unique endemic equilibrium, E t 1 is globally stable for R t0 > 1 if the population is constant over the modelling time.

Remark:
We notice from Equation (14) that the reproduction number due to the combined abuse of alcohol and methamphetamine is basically the maximum of either of the sub-model reproduction numbers. Therefore R at ¼ maxfR a0 ; R t0 g, where R a0 is the average number of new alcohol drinkers produced as a result of associating with an individual who drinks alcohol during his or her entire drinking life and R t0 is the average number of new methamphetamine users produced as a result of associating with an individual who uses methamphetamine, during his or her entire methamphetamine-using career.

Parameter estimation
In this section, we estimate the co-use model parameters to be used in the simulations. While most parameters are obtained from the fitting of the sub-models to alcohol and methamphetamine data, some parameters are simply estimated based on available literature, see (Bhunu & Mushayabasa, 2012a;Burattini et al., 1998;Nyabadza & Hove-Musekwa, 2010;Nyabadza et al., 2013;UNODC, World Drug Report, 2013). On average, the life expectancy in Sub-Sahara Africa and that of South Africa at the beginning of the modelling time, i.e. 1997 was about 50 years (Jamison et al., 2006). This, therefore, corresponds to a mortality rate of 0.02 per annum. We thus have the natural mortality rate, μ ¼ 0:02. According to Gray in (Day & Gray, 2008), the average birth rate in South Africa is approximated 0.028 per annum. To account for impact immigration into the Western Cape Province of South Africa, we shall set our recruitment rate Λ to be greater than 0.028 per annum. The actual mortality rate due to substance abuse is complex to estimate. This could be due to variations in risky behaviours by individuals while under the influence of substances of abuse (Burattini et al., 1998). This variation also extends within and among a population. The mortality rate, for example, among injecting crack-cocaine users in (Burattini et al., 1998) was 0.018 per year. According to a report in (NIDA, National Institute on Drug Abuse, 2018), a smokers life expectancy is increased by about 14% if he/she quits smoking at about age 35. Since treatment impacts positively on the quality of life, we assume that it reduces mortality rate related to substance abuse. Thus, we choose ζ 1 ¼ 0:005, ζ 3 ¼ 0:05, 1 ¼ 0:05 and 3 ¼ 0:00001 per year. The observed treatment demand for methamphetamine users was 17% in (Nyabadza & Hove-Musekwa, 2010). In this paper, we choose the average treatment demand of 30% as the corresponding treatment rate of 0.3. The summary of parameter values used in the numerical simulations is given in Table 2.
Observe from Figure 2 that when the co-abuse model reproduction number R at 0 < 1, the populations of alcohol and methamphetamine abusers with and without treatment decline to near zero in numbers. However, the susceptibles population approaches a constant, Λ μ . This shows that when R at 0 < 1; that is, (R a 0 < 1 and R t 0 < 1), the alcohol-methamphetamine co-abuse epidemic dies out and only the susceptible population remains. The given system approaches the substance-free equilibrium (SFE) which is consistent with Theorem 3.2. On the other hand, Figure 3 shows that the populations of susceptibles, methamphetamine and alcohol abusers in treatment and those without treatment initially increase and later level off at different heights when R at 0 > 1.
Graphs showing the population dynamics of S, U a , R t , U t , R t , U at and R at for R at 0 ¼ 0:6303 < 1. Used parameter values are available in Table 2.
This implies that the co-abuse system (4) stabilizes at an endemic equilibrium point R at 0 > 1. The co-abuse epidemic would persist when R at 0 > 1. Figure 3. Graphs showing the population dynamics of S, U a , R t , U t , R t , U at and R at for R at 0 ¼ 1:1254 > 1. Used parameter values are available in Table 2.

Sensitivity analysis
We use Latin Hypercube Sampling (LHS) to determine the Partial Rank Correlation Coefficients (PRCCs) with 1000 simulations per run to establish the sensitivity of the model parameters to the population of alcohol-methamphetamine abusers in the compartment U at . We observe from Figure 4, that the parameters with the greatest potential to worsen the co-use epidemic of alcohol and methamphetamine abuse are the effective person to person contact rates, β 1 and β 2 . Similarly, the relapse rate from the treatment to the class of methamphetamine users is observed to be highly significant, and may easily worsen the epidemic. Moreover, the quitting parameter, ρ 6 is the parameter with the greatest potential to make the epidemic better when increased.  Table 2.  Table 2 and 1,000 simulations per run were used.
Using the parameter values defined in Table 2, we illustrated the Monte Carlo simulations for the four parameters whose PRCC magnitudes are as shown in Figure 4. We observe from Figure 4, that the parameters β 1 , β 2 , γ 3 and Λ are positively correlated with the population of co-users of alcohol and methamphetamine (U at ). The parameters β 1 and β 2 are highly significant whereas the parameters Λ and γ 3 are less significant in the correlation. This means that an increase in these parameters will worsen the substance abuse epidemic. Those that show some negative correlation decrease the epidemic when they are increased. Typically, this is observed in the permanent quitting rates.
The PRCC scatter plots clearly show that the alcohol reproduction number is mainly influenced by the transmission parameter β 1 and the relapse rate γ 1 . Control measures against alcohol abuse should, therefore, focus on reducing new cases and ensuring that all the treatment facilities provide effective treatment services to minimize the rampant cases of relapse. The above trend is also observed in the methamphetamine model (see Figure 6).

Numerical results
In this subsection, we fit the model systems (15) and (37) to the data of individuals seeking treatment for alcohol and methamphetamine as their primary substance of abuse at specialised treatment centres in the Western Cape Province. The data were collected from 1996 to 2016 on a 6-month interval by SACENDU (Siphokazi et al., 2017) and is given in Tables 3 and 4. The least squares curve fit is used to fit the models to data. A Matlab code is used in which, the unknown parameter values are given a lower bound and an upper bound from which the set of parameter values that produce the best fit are obtained. The data on the demand for treatment are used to model growth in the R t class in the methamphetamine sub-model.  Table  (2). One thousand simulations per run were used. Figure 7 shows the graphical representation of model (15) fitted to data for persons seeking treatment for methamphetamine abuse. The circles and the solid line each represent the actual data points and the model fit to the data, respectively. We observe that the model fits well with the data. Furthermore, we notice that there was no population of methamphetamine users before 1996. The unavailability of such data could be attributed to a lack of data collection before the year 1996 (Nyabadza & Hove-Musekwa, 2010). Individuals may have been using the drugs but never really sort treatment services. It is equally important to notice that the population methamphetamine users under treatment peaked in the year 2006. The results showed a shortterm and fast-growing methamphetamine epidemic in which there was a significant increase in the number of users between 2002 and 2005, followed by a significant slowdown in the generation of new cases. The data shows an epidemic that is stabilizing at about 35% of the rehabilitants. The model also shows a steady state solution close to this value.
The data showing the demand for treatment for alcohol addiction are shown in Table 4. As with methamphetamine, the alcohol data were similarly collected by SACENDU (Siphokazi et al., 2017) from 1996 to 2016 in time intervals of 6 months.
From Figure 8, we observe that the proportion of individuals seeking treatment for alcohol abuse has been on a steady decline. The decline could be associated with heavy underreporting of alcohol abuse. Given the stigma associated with alcohol drinking in many communities in SA, people tend to under-report their consumption of alcohol, as it is perceived to be socially undesirable (Livingston & Callinan, 2015;Parry et al., 2005;Stockwell et al., 2004). Table 3. Primary Methamphetamine-abuse for the period 1997a to 2016b in % (Source: (Siphokazi et al., 2017)). The letters a and b represent the first 6 months (January-June) and the second 6 months (July-December) of the year, respectively Year 1997a 1997b 1998a 1998b 1999a 1999b   Moreover, it is also thought that some alcohol users may have found refuge in other substances/drug of abuse, thereby quitting alcohol or reducing their alcohol intake. The presented prevalence trends from data collected in WC by SACENDU may therefore not be representative of the actual prevalence and intensity of alcohol intake and abuse in the region. Further research work is thus necessary to establish why very few people only seem to recognize other substances/drugs and not alcohol as being part of their substance abuse problem. One way to address the current under-reporting in alcohol abuse is to use data on binge drinking to assist policymakers to monitor trends in alcohol use, especially risky drinking (Probst,  Time ( Table (2).
The projected population of alcohol and methamphetamine patients under treatment in the Western Cape Province until the year 2027 are given in Figure 9. We observe that the population of methamphetamine users would remain constant for the 9-year period. This is supported by the data that seems to fluctuate around this steady state. The population of alcohol users under treatment is however shown to be on the decline. The decline once again could be attributed to the growing popularity of other substances/drugs of abuse such as cannabis and methamphetamine in the region coupled with under-reporting.

Discussion and conclusion
In this paper, a deterministic model of co-use of alcohol and methamphetamine is presented and analysed. The model is decomposed into two sub-models; the alcohol submodel and methamphetamine submodel. We establish that the co-use model reproduction number is the maximum of the sub-model reproduction numbers, R a0 and R t0 . By constructing a suitable lyapunov function, the analysis shows that the alcohol-free equilibrium of the model is globally asymptotically stable whenever R a0 is less than unity. The health implication of this observation is that keeping the reproduction numbers below unity may be necessary for the effort to control the growth of the alcohol epidemic and subsequently the co-use epidemic of alcohol and methamphetamine.
The results from the sensitivity analysis suggest that the control of the combined abuse of methamphetamine and alcohol pivots around social intervention programs aimed at new users. Limiting new cases through preventive measures such as educational campaigns and other activities that encourage new users to quit. Secondly, cases of relapse should be minimised. Furthermore, there is the need to adopt in-patient treatment methods as opposed to the usual out-patient method, as this would limit contact with other drug users and hence lowering new relapse cases.
During treatment, alcohol-methamphetamine co-abuse could be managed by ensuring that treatment providers and rehabilitation centres treat patients for their concurrent substance abuse aggressively or refer them appropriately to other medical facilities. The providers should try to understand and address the underlying causes of concurrent substance use among visiting patients. This is vital for future policy formulations.
The model presented has emphasised that data collection on substance abuse should not be restricted to the primary substances of abuse but should also extend to multiple abuses of substances/drugs. The different substance combinations should be well documented. Owing to the numerous modelling assumptions made in this paper, the presented co-use substance abuse model may not provide an exact and a true representation of the epidemic resulting from the couse of alcohol and methamphetamine in the Western Cape Province of South Africa. The dynamics  are much more complex than described here. We argue that the presented model provides the required insights to understanding co-use in substance abuse epidemic. The results are relevant in designing intervention strategies aimed at combating multiple substance abuse in the Western Cape Province of South Africa.