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

Abstract

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 β₁ and β₂ 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.

Keywords:

• alcohol
• methamphetamine
• multiple substance/drug abuse
• alcohol-methamphetamine co-abuse
• latin hypercube sampling
• partial rank correlation coefficients

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.

1. Introduction

Substance/drug abuse poses a significant threat to the health, and socio-economic fabric of individuals and nations (WHO, 2018). It is estimated that out of a quarter of a billion people who used drugs in 2015, about 29.5 million suffered substance/drug-related disorders (UNODC, 2018). South Africa (SA) has experienced a continuous rise in multiple substance/drug abuse, especially in the recent past (Pluddemann et al., 2012). According to a report by the South African Community Epidemiology Network on Drug and Use (SACENDU), alcohol still remains the primary and most preferred substance of abuse among patients seen at specialist treatment centres across the country (Siphokazi et al., 2017). Drugs such as methamphetamine (MA), cannabis, cocaine and heroin have increasingly become popular among alcohol and other drug users in SA (Pluddemann et al., 2012). MA, cannabis and alcohol remain the most popular substances of abuse in the Western Cape (WC) Province of South Africa (Siphokazi et al., 2017), together comprising 91% of all admissions in the province. In the period July–December 2016, 29% of patients admitted to rehabilitation centres in the WC were being treated for MA and cannabis abuse, respectively. In the same period and region, 21% of all admissions were patients of alcoholism (Mposo, 2018). These numbers do not only present a worrying trend but a serious public health quandary.

Alcohol is a psychoactive substance with dependence-producing properties (2018). The harmful use of alcohol results in 3.3 million deaths each year (2018). Alcohol is known to have deleterious effects on the body’s immune system (NIH, 2018). The short-term effects of abuse may feature slowed down the activity of the brain causing speech to slow and the irregular drift in body temperatures (NIH, 2018). The long-term effects of alcoholism may include the emergence of serious organ diseases and socio-economic burden in families. Moreover, addicts may experience extreme pain during feeding as the digestive organs become inflammable (Barclay, Barbour, Stewart, Day, & Gilvarry, 2008).

MA, on the other hand, is a highly addictive illicit drug with over 26 million users worldwide (UNODC, 2018). The demand, spread and abuse of this vicious stimulant have increased dramatically in the recent past. MA (aka “tik” in WC) is a bitter white crystalline powder that is made up of ingredients that are readily available (Pluddemann et al., 2012). The neurotoxic effects of MA use have been well documented in (Chang, Alicata, Ernst, & Volkow, 2007; Kalechstein, Newton, & Green, 2003; Nordahl, Salo, & Leamon, 2003). The tremendous growth in the population of MA users in the Western Cape Province is documented to have reached epidemic levels (Siphokazi et al., 2017).

Despite the numerous negative consequences to the hosts, the trend of multiple substance abuse has continued to grow in recent years. Many young South Africans tend to use and abuse drugs in combination (SANCA, 2018). According to the 2017 SACENDU report (Siphokazi et al., 2017), multi-substance abuse remained high in the country, with 48% of patients in the WC region reporting using more than one substance of abuse. Clinical research indicates that the habit of taking several substances/drugs degrades the process of detoxification during treatment programmes (Center for Substance Abuse Treatment, 2005b). Denial by substance/drug abuse patients to using multiple substances/drugs of abuse only complicates the process of treatment/rehabilitation.

Efforts to resolve multi-substance abuse is further complicated by the fact that most if not all, the treatment centres and medical facilities in SA only provide experienced services at treating single substance/drug addiction cases (Bower, 1985). Moreover, Data on individuals who abuse more than one substance are either not available at all or is poorly document in SA. It has been shown that recovering from substance/drug abuse is easier when only one drug is involved (Center 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).

The concept of co-infection has been captured in several infectious disease models, see for instance, HIV and TB coinfection (Bolarin & Omatola, 2006; Naresh & Tripathi, 2005; Pawlowski, Jansson, Skold, Rottenberg, & Kallenius, 2012; Roeger, Feng, & Castillo-Chavez, 2009), HIV and hepatitis C coinfection (Kalichman et al., 2015), Malaria and Lymphatic Filariasis (Slater, Gambhir, Parham, & Michael, 2013). The understanding of how co-infection affects the dynamics and control of substance/drug abuse has become increasingly significant. Such dynamics have in the past provided meaningful evidence in pathogenic within-host interactions (Alter, 2006; Brooker, Akhwale, Pullan, Estambale, & Clarke et al., 2007; Cattadori, Boag, & Hudson, 2008; Cox, 2001; Petney & Andrews, 1998). Just like diseases, multi-substance abuse phenomenon is likely to change the course of “infection” of “co-infected” individuals. However, the actual mechanisms by which new cases of substance abuse are aggravated by multi-substance abuse are poorly understood.

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 multi-substance/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, co-infection 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.

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

(1) $N(t)=U_a(t)+R_a(t)+U_t(t)+R_t(t)+U_{at}(t)+R_{at}(t).$(1)

2.1. Basic assumptions

• Individuals get recruited into the susceptible population through birth and immigration at a constant rate, $\mathrm{\Lambda }$.

• 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 ${\lambda }_1$ with

(2) ${\lambda }_1={\beta }_1\left(\frac{U_a+{\zeta }_1{R}_a+{\zeta }_2{U}_{at}+{\zeta }_3{R}_{at}}{N}\right),$(2)

with the parameter ${\beta }_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, ${\zeta }_1,{\zeta }_3<1,$ while ${\zeta }_2>1.$ Similarly, individuals get initiated into using methamphetamine at the rate given by

(3) ${\lambda }_2={\beta }_2\left(\frac{U_t+{ϵ}_1{R}_t+{ϵ}_2{U}_{at}+{ϵ}_3{R}_{at}}{N}\right),$(3)

where the parameter ${\beta }_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 ${\eta }_a{\lambda }_2$ and move into the compartment $U_{at}$ with ${\eta }_t\ge 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 ${\delta }_1$ or die naturally at a rate $\mu ,$ 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 ${\sigma }_1$ or altogether, quit alcohol drinking at a rate ${\rho }_3.$ Alcohol users, upon successful treatment in the compartment $R_a$, permanently quit alcohol at the rate ${\rho }_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 ${\lambda }_2(t)$ through contact with those using methamphetamine. Some individuals in $U_t$ acquire alcohol drinking habits at the rate ${\eta }_t{\lambda }_1$ and move into the class $U_{at}$ with ${\eta }_t\ge 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 ${\sigma }_3$ and may relapse back to $U_t$ at a rate ${\gamma }_3.$ Individuals in the class $U_t$ may die due to methamphetamine-related causes at a rate ${\delta }_3$ or permanently quit methamphetamine at a rate ${\rho }_7.$

Users of both alcohol and methamphetamine in the class, $U_{at}$ may revert to alcohol abuse only at a rate ${\rho }_1$ or to methamphetamine abuse only at a rate ${\rho }_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 ${\delta }_2.$ They progress into the treatment class $R_{at}$ at a rate ${\sigma }_2.$ Individuals under treatment for alcohol and methamphetamine abuse in $R_{at},$ may relapse into alcohol-methamphetamine abuse compartment at a rate ${\gamma }_2$ or quit drinking of alcohol and relapse to $U_t$, at a rate ${\gamma }_5$. Similarly, they may quit methamphetamine and revert to the class of alcohol users not under treatment $U_a$ at a rate ${\gamma }_4$. We also allow permanent quitting for individuals in the classes $R_{at}$ and $R_t$ at rates ${\rho }_5$ and ${\rho }_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.

Figure 1. A compartmental representation of the epidemic of alcohol-methamphetamine co-use.

A description of parameters used in the co-abuse model and its sub-models are defined in Table 1.

Table 1. Descriptions of co-abuse model parameters

Parameter	Description
$\mathrm{\Lambda }$	Recruitment rate for susceptible persons
$\mu$	Natural mortality rate
${\beta }_1$	Alcohol transmission rate
${\beta }_2$	Methamphetamine transmission rate
${\gamma }_1$	Relapse rate for alcohol abusers
${\gamma }_2$	Relapse rate for co-abuser
${\gamma }_3$	Relapse rate for methamphetamine abusers
${\gamma }_4$	Co-abusers’ relapse rate into alcohol only abuse
${\gamma }_5$	Co-abusers’ relapse rate into methamphetamine only abuse
${\delta }_1$	Alcohol-induced mortality rate
${\delta }_2$	Co-abuse-induced mortality rate
${\delta }_3$	Methamphetamine-induced mortality rate
${\sigma }_1$	Recovery rate for alcohol abusers
${\sigma }_2$	Recovery rate for co-abusers of alcohol and methamphetamine
${\sigma }_3$	Recovery rate for methamphetamine abusers
${\rho }_1$	Rate at which co-abusers revert to methamphetamine abuse only without treatment
${\rho }_2$	Rate at which co-abusers revert to alcoholism only without treatment
${\rho }_3$	Quitting rate for alcohol abusers without treatment
${\rho }_4$	Quitting rate for alcohol abusers upon success treatment
${\rho }_5$	Quitting rate for co-abusers upon successful treatment
${\rho }_6$	Quitting rate for methamphetamine abusers upon successful treatment
${\rho }_7$	Quitting rate for methamphetamine abusers without treatment
${\eta }_a$, ${\eta }_t$, ${ϵ}_1$, ${ϵ}_2$	Modification parameters
${\zeta }_1$, ${\zeta }_2$, ${\zeta }_3$, ${ϵ}_3$	Modification parameters

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

(4) $\left.\begin{array}{cc}\frac{\mathrm{d}S}{\mathrm{d}t}=& \mathrm{\Lambda }-(\mu +{\lambda }_1+{\lambda }_2)S,\\ \frac{\mathrm{d}{U}_a}{\mathrm{d}t}=& {\lambda }_{1}S+{\rho }_1{U}_{at}+{\gamma }_1{R}_a+{\gamma }_4{R}_{at}-{\eta }_a{\lambda }_2{U}_a-(\mu +{\sigma }_1+{\delta }_1+{\rho }_3)U_a,\\ \frac{\mathrm{d}{R}_a}{\mathrm{d}t}=& {\sigma }_1{U}_a-(\mu +{\gamma }_1+{\rho }_4)R_a,\\ \frac{\mathrm{d}{U}_t}{\mathrm{d}t}=& {\lambda }_{2}S+{\rho }_2{U}_{at}+{\gamma }_3{R}_t+{\gamma }_5{R}_{at}-{\eta }_t{\lambda }_1{U}_t-(\mu +{\sigma }_3+{\delta }_3+{\rho }_7)U_t,\\ \frac{\mathrm{d}{R}_t}{\mathrm{d}t}=& {\sigma }_3{U}_t-(\mu +{\gamma }_3+{\rho }_6)R_t,\\ \frac{\mathrm{d}{U}_{at}}{\mathrm{d}t}=& {\eta }_a{\lambda }_2{U}_a+{\eta }_t{\lambda }_1{U}_t+{\gamma }_2{R}_{at}-(\mu +{\rho }_1+{\rho }_2+{\delta }_2+{\sigma }_2)U_{at},\\ \frac{\mathrm{d}{R}_{at}}{\mathrm{d}t}=& {\sigma }_2{U}_{at}-(\mu +{\gamma }_2+{\gamma }_4+{\gamma }_5+{\rho }_5)R_{at},\end{array}\right\}$(4)

$S(0)=S_0,U_a(0)=U_{a0},R_a(0)=R_{a0},U_t(0)=U_{t0},R_t(0)=R_{t0},U_{at}(0)=U_{at0},R_{at}(0)=R_{at0}$

where

(5) ${\lambda }_1={\beta }_1\left(\frac{U_a+{\zeta }_1{R}_a+{\zeta }_2{U}_{at}+{\zeta }_3{R}_{at}}{N}\right)\mathrm{a}\mathrm{n}\mathrm{d}{\lambda }_2={\beta }_2\left(\frac{U_t+{ϵ}_1{R}_t+{ϵ}_2{U}_{at}+{ϵ}_3{R}_{at}}{N}\right).$(5)

3. Mathematical analysis of the co-abuse model

3.1. 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\ge 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,

(6) $\underset{t\to \mathrm{\infty }}{lim}N(t)\le \frac{\mathrm{\Lambda }}{\mu }.$(6)

Proof. Assume that $T$ is the maximum time for the epidemic. That is, assume, $T=\mathrm{s}\mathrm{u}\mathrm{p}\{t>0,S>0,U_a\ge 0,R_a\ge 0,U_t\ge 0,R_t\ge 0,U_{at}\ge 0and{R}_{at}\ge 0\}\in [0,t]$. Therefore, $T>0$ and from the first equation of model system (4), we obtain

$\frac{\mathrm{d}}{\mathrm{d}t}\left[S(t)exp\left\{\mu t+{\int }_0^t({\lambda }_1(s)+{\lambda }_2(s))\mathrm{d}s\right\}\right]=\mathrm{\Lambda }exp\left[\mu t+{\int }_0^t({\lambda }_1(s)+{\lambda }_2(2))\mathrm{d}s\right],$

So that

$S(T)\ge S(0)exp\left[-\left\{\mu T+{\int }_0^T({\lambda }_1(s)+{\lambda }_2(s))\mathrm{d}s\right\}\right]$

$exp\left[-\left\{\mu T+{\int }_0^T({\lambda }_1(s)+{\lambda }_2(s))\mathrm{d}s\right\}\right]\left({\int }_0^T\mathrm{\Lambda }exp\left[\mu T+{\int }_0^{\hat{T}}(({\lambda }_1(w)+{\lambda }_2(w))\mathrm{d}w\right]\mathrm{d}\hat{T}\right)>0.$

Hence, $S(T)\ge 0\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{\forall }T>0.$

Now, from Equation (4), we have;

$\frac{\mathrm{d}{U}_a}{\mathrm{d}t}={\lambda }_{1}S+{\rho }_1{U}_{at}+{\gamma }_1{R}_a+{\gamma }_1{R}_{at}-{\eta }_a{\lambda }_2{U}_a-(\mu +{\sigma }_1+{\delta }_1+{\rho }_3)U_a,$

$\ge -({\eta }_a{\lambda }_2-\mu +{\sigma }_1+{\delta }_1+{\rho }_3)U_a,$

(7) $U_a(t)\ge U_{a0}exp\left[-\{(\mu +{\sigma }_1+{\delta }_1+{\rho }_3)t+{\int }_0^t{\eta }_a{\lambda }_2(s)\mathrm{d}s\}\right]>0.$(7)

It can similarly be shown that $U_t\ge 0$, $R_t\ge 0$, $U_{at}\ge 0$ and $R_{at}\ge 0$. Thus the solutions of (4) remain positive for all $t\ge 0$.

The evolution change in the population is given by

$\frac{\mathrm{d}N}{\mathrm{d}t}=\mathrm{\Lambda }-\mu N-({\delta }_1{U}_a+{\delta }_2{U}_{at}+{\delta }_3{U}_t)-({\rho }_3{U}_a+{\rho }_4{R}_a+{\rho }_5{R}_{at}+{\rho }_6{R}_t+{\rho }_7{U}_t),$

$\le \mathrm{\Lambda }-\mu N.$

It can easily be seen that

(8) $N(t)\le \frac{\mathrm{\Lambda }}{\mu }+\left(N_0-\frac{\mathrm{\Lambda }}{\mu }\right)exp(-\mu t).$(8)

From Equation (8), we observe that as $t\to \mathrm{\infty }$, $N(t)\to \frac{\mathrm{\Lambda }}{\mu }$. So that if $N_0\le \frac{\mathrm{\Lambda }}{\mu }$ then ${lim}_{t\to \mathrm{\infty }}N(t)=\frac{\mathrm{\Lambda }}{\mu }$. Clearly $\frac{\mathrm{\Lambda }}{\mu }$ is the upper bound for of $N$. On the other hand, if $N_0>\frac{\mathrm{\Lambda }}{\mu }$, then $N(t)$ will decrease to $\frac{\mathrm{\Lambda }}{\mu }$ as $t\to \mathrm{\infty }$. This means that if $N_0>\frac{\mathrm{\Lambda }}{\mu }$, 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

(9) $D=\left\{(S,U_a,R_a,U_t,R_t,U_{at},R_{at})\in {\mathbb{R}}_{+}^7|N\le \frac{\mathrm{\Lambda }}{\mu }\right\}$(9)

The region $D$ is positively invariant.

3.2. Steady states and stability analysis of the co-abuse model

The substance abuse free equilibrium ${\mathcal{E}}_{at}^0$ of the co-infection model (4) exists when there are no substance users in the community. That is, $U_a=0$, $R_a=0$, $U_t=0$, $R_t=0$, $U_{at}=0$, $R_{at}=0$. Therefore,

(10) ${\mathcal{E}}_{at}^0=\left(\frac{\mathrm{\Lambda }}{\mu },0,0,0,0,0,0\right).$(10)

Following (van Den Driessche & Watmough, 2002), we have the following results on the local stability of ${\mathcal{E}}_{at}^0$.

Theorem 3.2. The substance abuse free equilibrium, ${\mathcal{E}}_{at}^0$, is locally asymptotically stable if ${\mathcal{R}}_{at}<1$ and unstable if ${\mathcal{R}}_{at}>1$.

Next, we calculate the co-abuse endemic steady state(s), ${\mathcal{E}}_{at}^{\ast }$. 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 ${\lambda }_1$ and ${\lambda }_2$) from the model system (4), we obtain

$S^{\ast }=\frac{\mathrm{\Lambda }}{\mu +{\lambda }_1^{\ast }+{\lambda }_2^{\ast }}$, $R_{at}^{\ast }=\frac{{\sigma }_2}{b_6}{U}_{at}^{\ast }$, $R_t^{\ast }=\frac{{\sigma }_3}{b_4}{U}_t^{\ast }$, $R_a^{\ast }=\frac{{\sigma }_1}{b_2}{U}_a^{\ast }$, $U_a^{\ast }={\mathrm{\Theta }}_1+{\mathrm{\Theta }}_2{U}_{at}^{\ast }$, $U_t^{\ast }={\mathrm{\Theta }}_3+{\mathrm{\Theta }}_4{U}_{at}^{\ast }$, $U_{at}^{\ast }=\frac{{\eta }_t{\lambda }_1{\mathrm{\Theta }}_3+{\eta }_a{\lambda }_2{\mathrm{\Theta }}_1}{b_5(1-{\mathrm{\Phi }}_3)-({\eta }_a{\lambda }_2{\mathrm{\Theta }}_2+{\eta }_t{\lambda }_1{\mathrm{\Theta }}_4)},$ where,

${\mathrm{\Theta }}_1=\frac{{\mathcal{J}}_1}{{\mathrm{\Psi }}_1},{\mathrm{\Theta }}_2=\frac{{\mathcal{K}}_1}{{\mathrm{\Psi }}_1},{\mathrm{\Theta }}_3=\frac{{\mathcal{J}}_2}{{\mathrm{\Psi }}_2},{\mathrm{\Theta }}_4=\frac{{\mathcal{K}}_2}{{\mathrm{\Psi }}_2},{\mathcal{J}}_1=\frac{\mathrm{\Lambda }{\lambda }_1{b}_6}{\mu +{\lambda }_1+{\lambda }_2},{\mathcal{J}}_2=\frac{\mathrm{\Lambda }{\lambda }_2{b}_6}{\mu +{\lambda }_1+{\lambda }_2},$

${\mathcal{K}}_1=b_6{\rho }_1+{\sigma }_2{\gamma }_4,{\mathcal{K}}_2=b_6{\rho }_2+{\sigma }_2{\gamma }_5,{\mathrm{\Psi }}_1=b_1{b}_6\left(1-{\mathrm{\Phi }}_1+\frac{{\eta }_a{\lambda }_2}{b_1}\right),{\mathrm{\Psi }}_2=b_3{b}_6\left(1-{\mathrm{\Phi }}_2+\frac{{\eta }_t{\lambda }_1}{b_3}\right),$ ${\mathrm{\Phi }}_1=\frac{{\sigma }_1{\gamma }_1}{b_1{b}_2},{\mathrm{\Phi }}_2=\frac{{\sigma }_3{\gamma }_3}{b_3{b}_4},{\mathrm{\Phi }}_3=\frac{{\sigma }_2{\gamma }_2}{b_5{b}_6},b_1=\mu +{\sigma }_1+{\delta }_1+{\rho }_3,b_2=\mu +{\gamma }_1+{\rho }_4,$$b_3=\mu +{\sigma }_3+{\delta }_3+{\rho }_7,b_4=\mu +{\gamma }_3+{\rho }_6,b_5=\mu +{\rho }_1+{\rho }_2+{\delta }_2+{\sigma }_2,b_6=\mu +{\gamma }_2+{\gamma }_4+{\gamma }_5+{\rho }_5.$

Given that

(11) $\frac{{\lambda }_1}{{\lambda }_2}=\frac{{\beta }_1(U_a+{\zeta }_1{R}_a+{\zeta }_2{U}_{at}+{\zeta }_3{R}_{at})}{{\beta }_2(U_t+{ϵ}_1{R}_t+{ϵ}_2{U}_{at}+{ϵ}_3{R}_{at})},$(11)

we have,

(12) $p({\lambda }_1,{\lambda }_2)={\lambda }_1({\alpha }_0+{\alpha }_1{\lambda }_1+{\alpha }_2{\lambda }_2+{\alpha }_3{\lambda }_1{\lambda }_2+{\alpha }_4{\lambda }_1^2+{\alpha }_5{\lambda }_2^2),$(12)

where

${\alpha }_0=b_1{b}_2{b}_3{b}_4{b}_5{b}_6(1-{\mathrm{\Phi }}_1)(1-{\mathrm{\Phi }}_2)(1-{\mathrm{\Phi }}_3)({\mathcal{R}}_{t0}-{\mathcal{R}}_{a0}),$

${\alpha }_1=b_4{\beta }_1{\eta }_t(b_2+{\zeta }_1{\sigma }_1)({\mathcal{K}}_1-b_5{b}_6(1-{\mathrm{\Phi }}_3))+\mathcal{M}{\beta }_2(b_6{ϵ}_2+{\sigma }_2{ϵ}_3)+{\mathcal{K}}_2{b}_2{\beta }_2{\eta }_a(b_4+{ϵ}_1{\sigma }_3),$

${\alpha }_2=-{\beta }_1(b_6{\zeta }_2+{\zeta }_3{\sigma }_2)\mathcal{M}-{\mathcal{K}}_1{b}_4{\beta }_1{\eta }_t(b_2+{\zeta }_1{\sigma }_1)+b_2{\beta }_2{\eta }_a(b_4+{ϵ}_1{\sigma }_3)(b_5{b}_6(1-{\mathrm{\Phi }}_3)-{\mathcal{K}}_1),$

${\alpha }_3={\alpha }_4+{\alpha }_5,$

${\alpha }_4=b_2{b}_4{\beta }_2{\eta }_a{\eta }_t(b_6{ϵ}_2+{ϵ}_3{\sigma }_2),\mathrm{a}\mathrm{n}\mathrm{d}{\alpha }_5=-b_2{b}_4{\beta }_1{\eta }_a{\eta }_t(b_6{\zeta }_2+{\zeta }_3{\sigma }_2)$

and $\mathcal{M}=b_2{b}_4(b_3{\eta }_a(1-{\mathrm{\Phi }}_2)+b_1{\eta }_t(1-{\mathrm{\Phi }}_1)).$

Note that if ${\lambda }_1=0$, then clearly ${\lambda }_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

(13) $p({\lambda }_1,{\lambda }_2)={\alpha }_0+{\alpha }_1{\lambda }_1+{\alpha }_2{\lambda }_2+{\alpha }_3{\lambda }_1{\lambda }_2+{\alpha }_4{\lambda }_1^2+{\alpha }_5{\lambda }_2^2.$(13)

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.

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

(14) ${\mathcal{R}}_{at}=max\left\{\frac{{\beta }_1[b_2+{\zeta }_1{\sigma }_1]}{b_1{b}_2[1-{\mathrm{\Phi }}_1]},\frac{{\beta }_2[b_4+{ϵ}_1{\sigma }_3]}{b_3{b}_4[1-{\mathrm{\Phi }}_2]}\right\},$(14)

where ${\mathrm{\Phi }}_1=\frac{{\gamma }_1{\sigma }_1}{b_1{b}_2}$ and ${\mathrm{\Phi }}_2=\frac{{\gamma }_3{\sigma }_3}{b_3{b}_4}.$

The term $\frac{{\beta }_1[b_2+{\zeta }_1{\sigma }_1]}{b_1{b}_2[1-{\mathrm{\Phi }}_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 $\frac{{\beta }_2[b_4+{ϵ}_1{\sigma }_3]}{b_3{b}_4[1-{\mathrm{\Phi }}_2]}$represents the contribution to new “infections” resulting from one user of methamphetamine in a completely susceptible population.

4. Alcohol abuse model

This model is given when $U_t=R_a=U_{at}=R_{at}=0$. The co-abuse model (4) reduces to the following system of equations:

(15) $\left.\begin{array}{c}\frac{\mathrm{d}S}{\mathrm{d}t}=\mathrm{\Lambda }-(\mu +{\bar{\lambda }}_1)S,\\ \frac{\mathrm{d}{U}_a}{\mathrm{d}t}={\bar{\lambda }}_{1}S+{\gamma }_1{R}_a-(\mu +{\sigma }_1+{\delta }_1+{\rho }_3)U_a,\\ \frac{\mathrm{d}{R}_a}{\mathrm{d}t}={\sigma }_1{U}_a-(\mu +{\gamma }_1+{\rho }_4)R_a,\end{array}\right\}$(15)

where

${\bar{\lambda }}_1=\frac{{\beta }_1(U_a+{\zeta }_1{R}_a)}{N}.$

The alcohol model (15) has an alcohol-free equilibrium given by ${\mathcal{E}}_0^a=\left(\frac{\mathrm{\Lambda }}{\mu },0,0\right)$.

4.1. Reproduction number for the alcohol epidemic model (${\mathcal{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 $\mathrm{F}$ and $\mathrm{V}$ respectively at the alcohol-free equilibrium are given as follows:

(16) $\mathrm{F}=\left(\begin{array}{cc}{\beta }_1& {\beta }_1{\zeta }_1\\ 0& 0\end{array}\right)\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{V}=\left(\begin{array}{cc}(\mu +{\delta }_1+{\rho }_3+{\sigma }_1)& -{\gamma }_1\\ -{\sigma }_1& (\mu +{\gamma }_1+{\rho }_4)\end{array}\right).$(16)

It, therefore, follows that basic reproduction number due to alcohol abuse is given by

(17) ${\mathcal{R}}_{a0}=\frac{{\beta }_1[b_2+{\zeta }_1{\sigma }_1]}{b_1{b}_2[1-{\mathrm{\Phi }}_1]}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathrm{\Phi }}_1=\frac{{\gamma }_1{\sigma }_1}{b_1{b}_2}.$(17)

We now present the stability of the equilibrium points established from the alcohol abuse model (15).

4.2. Global stability of the alcohol-free equilibrium

Theorem 4.1. ${\mathcal{E}}_0^a$ is globally asymptotically stable in $\bar{\mathrm{\Omega }}\subset \mathrm{\Omega }$ if ${\mathcal{R}}_{a0}<1.$

Proof. Let

(18) $\mathrm{V}=b_2{U}_a+({\beta }_1{\zeta }_1+{\gamma }_1)R_a,$(18)

be a candidate Lyapunov function. The derivative of $\mathrm{V}$ with respect to time $(t)$ is given by

(19) $\frac{\mathrm{d}\mathrm{V}}{\mathrm{d}t}=[b_2({\beta }_1-b_1)+({\beta }_1{\zeta }_1+{\gamma }_1){\sigma }_1]U_a,$(19)

(20) $=b_1{b}_2(1-{\mathrm{\Phi }}_1)[{\mathcal{R}}_{a0}-1]U_a.$(20)

Noting that all the model parameters are positive, it follows that $\dot{\mathrm{V}}\le 0$ for $R_{ao}<1$ with $\dot{\mathrm{V}}=0$ only if $U_a=0$. Hence, $\mathrm{V}$ is a Lyapunov function on $\bar{\mathrm{\Omega }}=(S,U_a,R_a)$. Since $\bar{\mathrm{\Omega }}$ is invariant and attracting, it follows that the largest possible invariant set in $\{(S,U_a,R_a)\in \bar{\mathrm{\Omega }}:\dot{\mathrm{V}}=0\}$ is the singleton $\{{\mathcal{E}}_0^a\}$. 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 $\bar{\mathrm{\Omega }}$ approaches ${\mathcal{E}}_0^a$ as time approaches infinity. That is, as $t\to \mathrm{\infty }$, $(U_a(t),R_a(t))\to (0,0)$. Substituting for $U_a=R_a=0$ into the model system (15), gives $S\to \frac{\mathrm{\Lambda }}{\mu }$ as $t\to \mathrm{\infty }$. Thus $(S(t),U_a(t),R_a(t))\to \left(\frac{\mathrm{\Lambda }}{\mu },0,0\right)$ as $t\to \mathrm{\infty }$ for ${\mathcal{R}}_{ao}<1$ so that ${\mathcal{E}}_0^a$ is globally asymptotically stable in $\bar{\mathrm{\Omega }}$ if ${\mathcal{R}}_{a0}<1$. □

4.3. Endemic equilibrium and its stability

Upon equating model system (15) to zero and solving for $S^{\ast }$, $U_a^{\ast }$ and $R_a^{\ast }$ in terms of the force of infection ${\lambda }_1^{\ast }$, we obtain the following expressions.

(21) $S^{\ast }=\frac{\mathrm{\Lambda }}{\mu +{\bar{\lambda }}_1^{\ast }},U_a^{\ast }=\frac{\mathrm{\Lambda }{\bar{\lambda }}_1^{\ast }{b}_2}{b_1{b}_2(\mu +{\bar{\lambda }}_1^{\ast })(1-{\mathrm{\Phi }}_1)},R_a^{\ast }=\frac{\mathrm{\Lambda }{\bar{\lambda }}_1^{\ast }{\sigma }_1}{b_1{b}_2(\mu +{\bar{\lambda }}_1^{\ast })(1-{\mathrm{\Phi }}_1)}.$(21)

Furthermore, on substituting equations in (21) into (16), we obtain the polynomial

(22) ${\bar{\lambda }}_1^{\ast }\left\{(1-{\mathcal{R}}_{a0})+{\bar{\lambda }}_1^{\ast }\left(\frac{b_2+{\sigma }_1}{b_1{b}_2[1-{\mathrm{\Phi }}_1]}\right)\right\}=0.$(22)

From the polynomial (22), ${\bar{\lambda }}_1^{\ast }=0$ corresponds to the alcohol free-equilibrium. A scenario in which there is no alcohol abuse in the community.

Similarly, ${\bar{\lambda }}_1^{\ast }=\left\{\frac{b_1{b}_2(1-{\varphi }_1)}{b_2+{\gamma }_1}({\mathcal{R}}_{a0}-1)\right\}$ which exists for ${\mathcal{R}}_{a0}>1$ corresponds to the endemic equilibrium, a state in which alcohol abuse persists in the community. Therefore, the alcohol-endemic equilibrium point(s) is given by ${\mathcal{E}}_1^a=(S^{\ast },U_a^{\ast },R_a^{\ast })$ where $S^{\ast },U_a^{\ast }$ and $R_a^{\ast }$ are as described in Equation (21) and exists only if $R_{a0}>1$.

4.4. Global stability of the alcohol-endemic equilibrium point, (Anthuvan)

We show that if the population is constant, the endemic steady state ${\mathcal{E}}_1^a$, is globally asymptotically stable. We thus assume that

(23) $\mathrm{\Lambda }=\mu N+({\delta }_1+{\rho }_3)U_a+{\rho }_4{R}_a.$(23)

The differential equations in system (15) become

(24) $\left.\begin{array}{cc}\frac{\mathrm{d}S}{\mathrm{d}t}=& (\mu +{\delta }_1+{\rho }_3)U_a+(\mu +{\rho }_4)R_a-{\lambda }_{1}S,\\ \frac{\mathrm{d}{U}_a}{\mathrm{d}t}=& {\lambda }_{1}S+{\gamma }_1{R}_a-(\mu +{\delta }_1+{\rho }_3+{\sigma }_1)U_a,\\ \frac{\mathrm{d}{R}_a}{\mathrm{d}t}=& {\sigma }_1{U}_a-(\mu +{\gamma }_1+{\rho }_4)R.\end{array}\right\}$(24)

Since $N=S+U_a+R_a$, is a constant, we let $S$, $U_a$ and $R_a$:

$s=\frac{S}{N},v=\frac{U_a}{N}\mathrm{a}\mathrm{n}\mathrm{d}w=\frac{R_a}{N}.$

We have a non-dimensionalized system given by

(25) $\left.\begin{array}{cc}\dot{s}=& (\mu +{\delta }_1+{\rho }_3)v+(\mu +{\rho }_4)w-{\beta }_1(v+{\zeta }_{1}w)s,\\ \dot{v}=& {\beta }_1(v+{\zeta }_{1}w)s+{\gamma }_{1}w-(\mu +{\delta }_1+{\rho }_3+{\sigma }_1)v,\\ \dot{w}=& {\sigma }_{1}v-({\gamma }_1+{\rho }_4+\mu )w.\end{array}\right\}$(25)

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 ${\mathcal{E}}_1^a$ is globally asymptotically stable whenever ${\mathcal{R}}_{a0}$ is greater than unity, and the if the population is constant within the modelling time.

Proof. We propose a suitable Lyapunov function $\mathcal{V}$ such that

(26) $\mathcal{V}=\left(s-s^{\ast }-s^{\ast }ln\frac{s}{s^{\ast }}\right)+A\left(v-v^{\ast }-v^{\ast }ln\frac{v}{v^{\ast }}\right)+B\left(w-w^{\ast }-w^{\ast }ln\frac{w}{w^{\ast }}\right).$(26)

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

(27) $\frac{\mathrm{\partial }\mathcal{V}}{\mathrm{\partial }s}=\left(1-\frac{s^{\ast }}{s}\right),\frac{\mathrm{\partial }\mathcal{V}}{\mathrm{\partial }v}=A\left(1-\frac{v^{\ast }}{v}\right),\frac{\mathrm{\partial }\mathcal{V}}{\mathrm{\partial }w}=B\left(1-\frac{w^{\ast }}{w}\right),$(27)

are all zero at the corresponding alcohol-endemic steady state. That is, at the endemic steady states, $s=s^{\ast }$, $v=v^{\ast }$ and $w=w^{\ast }$. In addition, the second partial derivatives of $\mathcal{V}$ with respect to any of the three state variables are given as follows;

(28) $\frac{{\mathrm{\partial }}^2\mathcal{V}}{\mathrm{\partial }{s}^2}=\frac{s^{\ast }}{s^2},\frac{{\mathrm{\partial }}^2\mathcal{V}}{\mathrm{\partial }{v}^2}=A\frac{v^{\ast }}{v^2}\mathrm{a}\mathrm{n}\mathrm{d}\frac{{\mathrm{\partial }}^2\mathcal{V}}{\mathrm{\partial }{w}^2}=B\frac{w^{\ast }}{w^2}.$(28)

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

(29) $\dot{\mathcal{V}}=\left(1-\frac{s^{\ast }}{s}\right)\dot{s}+A\left(1-\frac{v^{\ast }}{v}\right)\dot{v}+B\left(1-\frac{w^{\ast }}{w}\right)\dot{w}.$(29)

Substituting for $\dot{s}$, $\dot{v}$ and $\dot{w}$ from system (25), we have

(30) $\dot{\mathcal{V}}=\left(1-\frac{s^{\ast }}{s}\right)[av+(\mu +{\rho }_4)w-{\beta }_1(v+{\zeta }_{1}w)s]+\left(1-\frac{v^{\ast }}{v}\right)A[{\beta }_1(v+{\zeta }_{1}w)s+{\gamma }_{1}w-b_{1}v]+\left(1-\frac{w^{\ast }}{w}B[{\sigma }_{1}v-b_{2}w\right).$(30)

where $a=\mu +{\delta }_1+{\rho }_3$, $b_1=\mu +{\gamma }_1+{\rho }_4$ and $b_2=\mu +{\delta }_1+{\rho }_3+{\sigma }_1$.

We now use the system of Equation (25) at ${\mathcal{E}}_1^a$ to obtain:

(31) $a=\frac{{\beta }_1(v^{\ast }+{\zeta }_1{w}^{\ast })s^{\ast }-(\mu +{\rho }_4)w^{\ast }}{v^{\ast }},b_1=\frac{{\beta }_1(v^{\ast }+{\zeta }_1{w}^{\ast })s^{\ast }+{\gamma }_1{w}^{\ast }}{v^{\ast }}\mathrm{a}\mathrm{n}\mathrm{d}{b}_2=\frac{{\sigma }_1{v}^{\ast }}{w^{\ast }}.$(31)

We now substitute the terms in Equation (31) into Equation (30) so that

$\dot{\mathcal{V}}=\left(1-\frac{s^{\ast }}{s}\right)\left[{\beta }_1{v}^{\ast }{s}^{\ast }\left(\frac{v}{v^{\ast }}-\frac{vs}{v^{\ast }{s}^{\ast }}\right)+{\beta }_1{\zeta }_1{w}^{\ast }{s}^{\ast }\left(\frac{v}{v^{\ast }}-\frac{ws}{w^{\ast }{s}^{\ast }}\right)+(\mu +{\rho }_4)w^{\ast }\left(\frac{v}{v^{\ast }}+\frac{w}{w^{\ast }}\right)\right]$

$+A\left(1-\frac{v^{\ast }}{v}\right)\left[{\beta }_1{v}^{\ast }{s}^{\ast }\left(-\frac{v}{v^{\ast }}-\frac{vs}{v^{\ast }{s}^{\ast }}\right)+{\beta }_1{\zeta }_1{w}^{\ast }{s}^{\ast }\left(-\frac{v}{v^{\ast }}+\frac{ws}{w^{\ast }{s}^{\ast }}\right)+{\gamma }_1{w}^{\ast }\left(\frac{v}{v^{\ast }}+\frac{w}{w^{\ast }}\right)\right]$

(32) $+B\left(1-\frac{w^{\ast }}{w}\right)\left[{\sigma }_1{v}^{\ast }\left(\frac{v}{v^{\ast }}-\frac{w}{w^{\ast }}\right)\right].$(32)

Let $\frac{s}{s^{\ast }}=J$, $\frac{v}{v^{\ast }}=K$, and $\frac{w}{w^{\ast }}=L$.

Therefore, $\dot{V}$ becomes

$\dot{\mathcal{V}}={\beta }_1{v}^{\ast }{s}^{\ast }\left[\left(2K-JK-\frac{K}{J}\right)+A(1+JK-(K+J))\right]+(\mu +{\rho }_4)\left[(L+K)-\left(\frac{L}{J}+\frac{K}{J}\right)\right]w^{\ast }$

$+{\beta }_1{\zeta }_1{w}^{\ast }{s}^{\ast }\left[(K+L)-\left(JL+\frac{K}{J}\right)+A(1+JL)-\left(K+\frac{JL}{K}\right)\right]+A{\gamma }_1{w}^{\ast }\left[(L+K)-\left(1+\frac{L}{K}\right)\right]$

(33) $+B{\sigma }_1{v}^{\ast }\left[(1+K)-\left(L+\frac{K}{L}\right)\right].$(33)

Setting the coefficients of $K$, $JK$, $JL$, $\frac{K}{J}$ to zero, we obtain

(34) $A=1\mathrm{a}\mathrm{n}\mathrm{d}B=\frac{{\beta }_1{\zeta }_1{w}^{\ast }{s}^{\ast }+(\mu +{\rho }_4)w^{\ast }}{{\sigma }_1{v}^{\ast }}>0.$(34)

Upon substituting the expressions in Equation (34) back into (30) we obtain

(35) $\dot{\mathcal{V}}=\left(1-\frac{s^{\ast }}{s}\right)\dot{s}+\left(1-\frac{v^{\ast }}{v}\right)\dot{v}+\left(1-\frac{w^{\ast }}{w}\right)\left(\frac{{\beta }_1{\zeta }_1{w}^{\ast }{s}^{\ast }+(\mu +{\rho }_4)w^{\ast }}{{\sigma }_1{v}^{\ast }}\right)\dot{w}.$(35)

Also, on substituting the values of $A$ and $B$ in Equation (34) into Equation (33), and on further simplification, we obtain

$\dot{V}={\beta }_1{v}^{\ast }{s}^{\ast }\left[(K-J)\left(1-\frac{1}{J}\right)\right]+{\beta }_1{\zeta }_1{w}^{\ast }{s}^{\ast }\left[(2+K-\left(\frac{K}{J}+\frac{K}{L}+\frac{JL}{K}\right)\right]$

(36) $+(\mu +{\rho }_4)w^{\ast }\left[1+2K-\left(\frac{L}{J}+\frac{K}{J}+\frac{K}{J}\right)\right]+{\gamma }_1{w}^{\ast }\left(L-\frac{L}{K}+K-1\right).$(36)

From Equation (36), we observe that the expression $(K-J)\left(1-\frac{1}{J}\right)$ 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)\left(1-\frac{1}{K}\right)\le 0$ with the equality holding if and only if $K=1$ or $L=K$.

We can draw similar conclusions from the remaining expressions $\left(2+K-\left(\frac{K}{J}+\frac{K}{L}+\frac{JL}{K}\right)\right)$ and $\left(1+2K-\left(\frac{L}{J}+\frac{K}{J}+\frac{K}{J}\right)\right)$. In the two cases, equality holds if and only if $K=J=L$.

Therefore, $\dot{\mathcal{V}}\le 0$ with equality holding if and only if $K=J=L$. Since $\dot{\mathcal{V}}=0$ only when $K=J=L$, which corresponds to $s=s^{\ast }$, $v=v^{\ast }$, and $w=w^{\ast }$ and subsequently, to $S=S^{\ast }$, $U_a=U_a^{\ast }$ and $R_a=R_a^{\ast }$, the largest invariant set in $\{(S,U_a,R_a)\in \bar{\mathrm{\Omega }}:\dot{\mathcal{V}}=0$ is the singleton $\{{\mathcal{E}}_1^a\}$. By LaSalle’s invariance principle (La Salle, 1976), we, therefore, conclude that the endemic equilibrium ${\mathcal{E}}_1^a$ is globally asymptotically stable in the interior of $\bar{\mathrm{\Omega }}$. This shows that every solution in $\bar{\mathrm{\Omega }}$ or that intersects $\bar{\mathrm{\Omega }}$, would approach the endemic equilibrium ${\mathcal{E}}_1^a$. □

5. Methamphetamine-abuse model

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

(37) $\left.\begin{array}{c}\frac{\mathrm{d}S}{\mathrm{d}t}=\mathrm{\Lambda }-(\mu +{\bar{\lambda }}_2)S,\\ \frac{\mathrm{d}{U}_t}{\mathrm{d}t}={\bar{\lambda }}_{2}S+{\gamma }_3{R}_t-(\mu +{\sigma }_3+{\delta }_3+{\rho }_7)U_t,\\ \frac{\mathrm{d}{R}_t}{\mathrm{d}t}={\sigma }_3{U}_t-(\mu +{\gamma }_3+{\rho }_6)R_t,\end{array}\right\}$(37)

where ${\bar{\lambda }}_2=\frac{{\beta }_2(U_t+{ϵ}_1{R}_t)}{N}.$

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

(38) ${\bar{\lambda }}_2^{\ast \ast }\left[b_3{b}_4(1-{\mathrm{\Phi }}_2)(1-{\mathcal{R}}_{t0})+{\bar{\lambda }}_2^{\ast \ast }\left({\sigma }_3+b_4\right)\right]=0$(38)

(39) $\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathcal{R}}_{t0}=\frac{{\beta }_2[b_4+{\zeta }_2{\sigma }_3]}{b_3{b}_4[1-{\mathrm{\Phi }}_2]}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Phi }}_2=\frac{{\gamma }_3{\sigma }_3}{b_3{b}_4}.$(39)

From Equation (38), ${\bar{\lambda }}_2^{\ast \ast }=0$ corresponds to the methamphetamine-free equilibrium, denoted by ${\mathcal{E}}_0^t$ and this is a scenario in which the community is free of methamphetamine abuse.

${\mathcal{E}}_0^t=(S^{\ast \ast },U_t^{\ast \ast },R_t^{\ast \ast })=\left(\frac{\mathrm{\Lambda }}{\mu },0,0\right).$

Note that ${\mathcal{R}}_{t0}$ is the methamphetamine-abuse reproduction number.

5.2. Stability analysis

From polynomial equation given in (38),

${\bar{\lambda }}_2^{\ast \ast }=\left\{\left(\frac{b_3{b}_4(1-{\mathrm{\Phi }}_2)}{{\sigma }_3+b_4}\right)(R_{t0}-1)\right\}$

which exists only if ${\mathcal{R}}_{t0}>1$, and corresponds to the methamphetamine-endemic equilibrium, denoted by ${\mathcal{E}}_1^t$ such that ${\mathcal{E}}_1^t=(S^{\ast \ast },U_t^{\ast \ast },R_t^{\ast \ast })$ where

(40) $S^{\ast \ast }=\frac{\mathrm{\Lambda }}{\mu +{\bar{\lambda }}_2^{\ast \ast }},U_t^{\ast \ast }=\frac{\mathrm{\Lambda }{\bar{\lambda }}_2^{\ast \ast }{b}_4}{b_3{b}_4(\mu +{\bar{\lambda }}_2^{\ast \ast })(1-{\mathrm{\Phi }}_2)},R_t^{\ast \ast }=\frac{\mathrm{\Lambda }{\bar{\lambda }}_2^{\ast \ast }{\sigma }_3}{b_3{b}_4(\mu +{\bar{\lambda }}_2^{\ast \ast })(1-{\mathrm{\Phi }}_2)}.$(40)

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, ${\mathcal{E}}_0^t$ is locally asymptotically stable for ${\mathcal{R}}_{t0}<1$ and unstable for ${\mathcal{R}}_{t0}>1.$

Theorem 5.2. The unique endemic equilibrium, ${\mathcal{E}}_1^t$ is globally stable for ${\mathcal{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 ${\mathcal{R}}_{at}=max\{{\mathcal{R}}_{a0},{\mathcal{R}}_{t0}\}$, where ${\mathcal{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 ${\mathcal{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.

6. Numerical simulation

6.1. 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, $\mu =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 $\mathrm{\Lambda }$ 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 ${\zeta }_1=0.005$, ${\zeta }_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.

Table 2. Description of parameters of the co-abuse model, their values and their ranges

Symbol	Value	Range	Source	Symbol	Value	Range	Source
$\mathrm{\Lambda }$	0.03	0.02–0.04	(Day & Gray, 2008)	${\sigma }_1$	0.2		Estimated
$\mu$	0.02	0.02–0.03	(Farai, Musekwa, & Senelani, 2010)	${\sigma }_2$	0.3		Estimated
${\beta }_1$	0.25	0.1–0.4	(Bhunu & Mushayabasa, 2012b)	${\sigma }_3$	0.3	0.09–0.3	(Farai et al., 2010)
${\beta }_2$	0.2	0.1–0.4	(Farai et al., 2010)	${\rho }_1$	0.35	0–1	(Bhunu & Mushayabasa, 2012b)
${\gamma }_1$	0.01	0–0.1	(Bhunu & Mushayabasa, 2012b)	${\rho }_2$	0.5	0–1	(Sharomi & Gumel, 2008)
${\gamma }_2$	0.1	0.1–0.5	Estimated	${\rho }_3$	0.35	0.2–0.5	(Bhunu & Mushayabasa, 2012b)
${\gamma }_3$	0.0011	0–0.9	(Kalula & Nyabadza, 1998)	${\rho }_3$	0.25	0–1	(Bhunu & Mushayabasa, 2012b)
${\gamma }_4$	0.3	0.1–0.5	Estimated	${\rho }_5$	0.25		Estimated
${\gamma }_5$	0.3	0.1–0.5	Estimated	${\rho }_6$	0.8312	0–91	(Kalula & Nyabadza, 1998)
${\delta }_1$	0.03	0.01–0.05	(Farai et al., 2010)	${\rho }_7$	0.0095	0–0.3	(Kalula & Nyabadza, 1998)
${\delta }_2$	0.04	0.01–0.05	Estimated	${\eta }_a$,${\eta }_t$,${\zeta }_2$${ϵ}_2$	$1.05$	1–3	Estimated
${\delta }_3$	0.033	0.01–0.05	(Bhunu & Mushayabasa, 2012b)	${\zeta }_1$, ${\zeta }_3$, ${ϵ}_1$, ${ϵ}_3$	$0.7$	0–1	Estimated

Observe from Figure 2 that when the co-abuse model reproduction number ${\mathcal{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, $\frac{\mathrm{\Lambda }}{\mu }$. This shows that when ${\mathcal{R}}_{at}0<1$; that is, (${\mathcal{R}}_{a}0<1$ and ${\mathcal{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 ${\mathcal{R}}_{at}0>1$. This implies that the co-abuse system (4) stabilizes at an endemic equilibrium point ${\mathcal{R}}_{at}0>1$. The co-abuse epidemic would persist when ${\mathcal{R}}_{at}0>1$.

Figure 2. Graphs showing the population dynamics of $S$, $U_a$, $R_t$, $U_t$, $R_t$, $U_{at}$ and $R_{at}$ for ${\mathcal{R}}_{at}0=0.6303<1$. Used parameter values are available in Table 2.

Figure 3. Graphs showing the population dynamics of $S$, $U_a$, $R_t$, $U_t$, $R_t$, $U_{at}$ and $R_{at}$ for ${\mathcal{R}}_{at}0=1.1254>1$. Used parameter values are available in Table 2.

6.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, ${\beta }_1$ and ${\beta }_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, ${\rho }_6$ is the parameter with the greatest potential to make the epidemic better when increased.

Figure 4. PRCCs for the full range of parameters from Table 2.

Figure 5. Graphs (a), (b), (c) and (d) show the Monte Carlo simulations for the four parameters with the greatest PRCC magnitude in the alcohol model. Parameter values in 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 ${\beta }_1$, ${\beta }_2$, ${\gamma }_3$ and $\mathrm{\Lambda }$ are positively correlated with the population of co-users of alcohol and methamphetamine ($U_{at}$). The parameters ${\beta }_1$ and ${\beta }_2$ are highly significant whereas the parameters $\mathrm{\Lambda }$ and ${\gamma }_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 ${\beta }_1$ and the relapse rate ${\gamma }_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).

Figure 6. Graphs (a), (b), (c) and (d) show the Monte Carlo simulations for the four parameters with the greatest PRCC magnitude in the methamphetamine (MA) model. Used parameter values can be found in Table (2). One thousand simulations per run were used.

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

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	2000a	2000b	2001a
% Meth users use	0	0.1	0	0.1	0.1	0.1	0.2	0.1	0.1
Year	2001b	2002a	2002b	2003a	2003b	2004a	2004b	2005a	2005b
% Meth users	0.3	0.3	0.8	2.3	2.3	10.7	19.3	26.1	34.7
Year	2006a	2006b	2007a	2007b	2008a	2008b	2009a	2009b	2010a
% Meth users	37.2	42.3	40.7	36.1	35.8	35.1	40.6	35.5	33.6
Year	2010b	2011a	2011b	2012a	2012b	2013a	2013b	2014a	2014b
% Meth users	35.1	35.3	38.8	34	33	28	33	33	35
Year	2015a	2015b	2016a	2016b
% Meth users	35	37	32	29

Table 4. Primary Alcohol abuse: by 6-month period (from 1997a to 2016b in % (Source: (Siphokazi et al., 2017))

Year	1997a	1997b	1998a	1998b	1999a	1999b	2000a	2000b	2001a
% Alcohol users	82.0	78.0	74.0	64.0	56.0	50.0	48.0	51.0	46.0
Year	2001b	2002a	2002b	2003a	2003b	2004a	2004b	2005a	2005b
% Alcohol users	46.0	48.0	47.0	43.6	39.4	38.3	33.7	34.4	25.1
Year	2006a	2006b	2007a	2007b	2008a	2008b	2009a	2009b	2010a
% Alcohol users	30.2	26.4	29.5	29.7	30.0	27.6	26.8	29.4	29.8
Year	2010b	2011a	2011b	2012a	2012b	2013a	2013b	2014a	2014b
% Alcohol users	27.5	27.5	23.7	24	22	20	22	20	22
Year	2015a	2015b	2016a	2016b
% Alcohol users	21	20	22	21

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.

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 short-term 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.

Figure 7. Shows methamphetamine model (37) fitted to data for addicts seeking treatment for methamphetamine as a primary substance of abuse in the Western Cape Province of South Africa. Parameter values are in Table (2).

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 under-reporting 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). 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, Shuper, & Rehm, 2017; Vellios & Van Walbeek, 2018).

Figure 8. Shows alcohol-model (15) fitted to data for individuals seeking treatment for alcohol abuse in Cape town and Western Cape province of South Africa. Parameter values are in 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.

Figure 9. Graphs showing the projected population of methamphetamine (MA), (a) alcohol (b) abuse patients seeking treatment services in Western Cape province, South Africa.

7. 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, ${\mathcal{R}}_{a0}$ and ${\mathcal{R}}_{t0}$. By constructing a suitable lyapunov function, the analysis shows that the alcohol-free equilibrium of the model is globally asymptotically stable whenever ${\mathcal{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 co-use 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.

Acknowledgements

We are grateful to the University of Stellenbosch, the African Institute for Mathematical Sciences (AIMS) and Strathmore University for their support during the production of this manuscript.

Disclosure statement

The authors declare that there are no competing interests regarding the publication of this paper.

Additional information

Funding

The authors received no direct funding for this research.

Notes on contributors

T. O. Orwa

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

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