Primary and secondary substance use in the Western Cape Province of South Africa: A mathematical modelling approach

Abstract

Substance abuse is an issue of considerable concern due to its devastating effects to families and communities in many countries. We formulate a mathematical model for multiple substance abuse using nonlinear ordinary differential equations. The model includes two classes of individuals in the substance using phase, namely, those using the primary substance of abuse only and those using the primary substance of abuse together with additional substance(s). An epidemic threshold value, ${\mathcal{R}}_{in}$, termed the initiation reproduction number, is proposed and defined herein in the substance-using context. It has been shown that the substance-free equilibrium is locally asymptotically stable when ${\mathcal{R}}_{in}<1$ and the system has a unique positive substance-persistent equilibrium if and only if ${\mathcal{R}}_{in}>1$. The least-squares curve fit routine (lsqcurvefit) in Matlab with optimization has been used to fit the model to data on primary and secondary substance abusers in treatment centres of Cape Town, South Africa. The results indicate a continued decrease in the proportion of primary substance abusers in Cape Town whereas secondary substance abusers shall continue to increase but steadily. Numerical results suggest that there is need to increase treatment efforts targeting secondary substance abusers as this will have more impact in the fight against substance abuse.

Keywords:

• multiple substance abuse
• initiation reproduction number
• least-squares curve fitting
• numerical simulations

PUBLIC INTEREST STATEMENT

The paper proposes a mathematical model describing the dynamics of substance abuse. The model builds upon prior work in mathematical epidemiology by Anderson and May (1991), and, Brauer and Castillo-Chavez (2001). The framework provides an extension of specific population level substance abuse models whose mathematical and numerical analysis provides additional insights into the dynamics of substance abuse with clearly identified set of assumptions. The model is innovative in three ways. First, the model explicitly captures the levels of substance abuse by classifying substance abuse as either primary substance abuse or secondary substance abuse, not considered in past substance abuse models. Second, the model is fitted to data on primary and secondary substance abusers in treatment centres of Cape Town, South Africa. Third, the model can be linked to resource-constrained communities in which rehabilitation processes are underfunded.

Competing interests

The authors declare no competing interests.

1. Introduction

Substance abuse is an issue of considerable concern due to the social and economic costs brought about by substance-related problems such as mental illness, road traffic accidents, work-related accidents, violence and crimes committed under the influence of substances (Van Heerden et al., 2009). Substance abuse as defined by Kaufman (Kaufman, 1976) refers to the use of a psychoactive substance resulting in serious interference of one’s physical health, social relationships or vocational functioning. Following (Gropper & Burke, 1974), “multi-substance abuse” shall herein be defined as the nonmedical abuse of more than one substance.

The major factor attributing to the use or abuse of drugs by man is to alter his mental state (Kaufman, 1976). The case of multiple substance abuse considered in this paper involves simultaneous intake of substances so close in time that when the second and later substances are taken, the effects of the primary substance will not have died out. We neglect the stepwise progression of substances which sometimes is the experience of substance abusers (Kaufman, 1976). There are many possible combinations of substances of abuse. Generally, substance users combine a preferred substance of abuse with other chosen substances in order to enhance the effects of the primary substance (American Addiction Centers, 2019). In this study, we shall consider only the critical pattern of two types of substance abuse, namely, the primary and secondary substance of abuse as captured by the South African Community Epidemiology Network on Drug Use (SACENDU, 2019).

Recently, several mathematical studies for the dynamics of drug abuse have been developed (see for instance, Kalula & Nyabadza, 2012; Mushanyu, Nyabadza, & Stewart, 2015; Mushanyu, Nyabadza, Muchatibaya et al., 2015; Mushanyu et al., 2016a, 2016b; Njagarah & Nyabadza, 2013; Nyabadza & Hove-Musekwa, 2010; Nyabadza et al., 2012; Rossi, 2002; White & Comiskey, 2007), but consequently, there has been no attempt to investigate the dynamics of multiple substance abuse. The model developed in this study is the first attempt to model the dynamics of multiple substance abuse with particular emphasis on primary and secondary substance abuse. The model developed in this paper is applied to a case study of the Western Cape Province of South Africa. Data used were obtained from SACENDU which monitors alcohol, tobacco and other drug abuse trends from specialist substance use treatment centres in all nine provinces that now comprise the South African Community Epidemiology Network on Drug Use. The provinces include the Western Cape (WC), KwaZulu-Natal (mostly Durban and Pietermaritzburg), Eastern Cape (Port Elizabeth and East London), Gauteng province, Mpumalanga and Limpopo provinces (now termed the Northern Region [NR]), and the Central Region (comprising of the Free State, Northern Cape and North-West provinces [CR]; SACENDU, 2019). Poly-substance use remains high in the Western Cape Province of South Africa as compared to other provinces, with majority of patients indicating the use of more than one substance of use.

The paper is arranged as follows: in Section $2$, we formulate and establish the basic properties of the model. The model is analysed for stability in Section $3$. In Section $4$, we carry out some numerical simulations and parameter estimation. Numerical results are also presented in this section. The paper is concluded in Section $5$.

2. Model formulation

We propose a mathematical model that takes into account the dynamics of multiple substance abuse. To this effect, we include two separate classes of individuals in the substance using phase, that is, those who abuse only one substance of abuse and those who use more than one substance of abuse. The human population is thus divided into seven distinct classes, namely, $S(t)$, $U_p(t)$, $T_p(t)$, $R_p(t)$, $U_s(t)$, $T_s(t)$ and $R_s(t)$. The class $S(t)$ represents the population at risk of being initiated into substance abuse. The class $U_p(t)$ denotes those initiated into substance abuse and only using the primary substance of abuse. The class $T_p(t)$ represents primary substance abusers in treatment. In this paper, we assume that treatment is on an “inpatient” basis and individuals are released from treatment centres when they have recovered. So we allow for the class of the recovered. Thus, the class of recovered primary substance abusers is denoted by $R_p(t)$. It is upon this stage that a relapse can only occur. Individuals in the recovered primary substance abusers’ class relapse back to the class $U_p(t)$ at a rate ${\rho }_p$. The class $U_s(t)$ denotes individuals at the secondary stage of substance abuse, that is, those abusing their primary substance of abuse together with additional substance(s). The class $T_s(t)$ represents secondary substance abusers in treatment. The class of recovered secondary substance abusers is denoted by $R_s(t)$. Individuals in the recovered secondary substance abusers’ class relapse back to the class $U_s(t)$ at a rate ${\rho }_s$. The total population is thus given by:

$N(t)=S(t)+U_p(t)+T_p(t)+R_p(t)+U_s(t)+T_s(t)+R_s(t).$

Susceptible humans enter the population through birth or immigration at a rate $\mathrm{\Lambda }$. Susceptible individuals are initiated into substance use following interaction with individuals using substances. Substance users in treatment cannot initiate non-substance users since they are completely immersed in the program and separated from the general population. An initiation function that is analogous to the force of infection for epidemic models is assumed. Thus, the per capita contact rate $\beta$ is a product of the effective number of contacts $c$, between substance users and the susceptible population, and the probability $\hat{\beta }$, that a contact results in initiation into substance use, that is, $\beta =c\hat{\beta }$. A fraction $U_p/N$ of the contacts is with the primary substance abusers and the average number of potential contacts resulting in susceptible individuals becoming substance abusers is $\beta U_p/N$. Also, a fraction $U_s/N$ of the contacts is with individuals at the secondary stage of substance abuse. The average number of potential contacts of each susceptible individual with individuals at the secondary stage of substance abuse is $\beta \eta U_s/N$. The parameter $\eta$ measures the relative ability for individuals in the secondary substance abuse phase to initiate new substance users. Assuming that the rate at which individuals in $U_p$ initiate new recruits is lower than that of individuals in $U_s$, it follows that $\eta >1$. This is due to the fact that individuals in class $U_s$ have an enhanced familiarity with substance abuse as compared to individuals in class $U_p$ and this gives them more potential to recruit new substance abusers. The total number of relevant contacts gives the initiation function,

(2.1) $\lambda =\beta \left(\frac{U_p+\eta U_s}{N}\right).$(2.1)

Upon being initiated into substance use, a susceptible individual moves into the compartment $U_p$, of primary substance abusers. We assume that there is a slow progression of individuals from the susceptible class to the substance using class. Thus, individuals in the secondary stage of substance abuse are as a result of primary substance abusers progressing to secondary substance abuse at a rate ${\alpha }_p$. The model developed in this paper also allows for an ameliorative process of individuals in class $U_s$. As such, individuals in the secondary substance use phase can revert back to the primary substance abuse class $U_p$ at a rate ${\alpha }_s$. Upon realizing the negative effects associated with substance abuse, individuals in compartments $U_p$ and $U_s$ can seek help via rehabilitation programmes at rates ${\sigma }_p$ and ${\sigma }_s$, respectively. Also, some individuals in compartments $U_p$ and $U_s$ can naturally recover at rates ${\gamma }_p$ and ${\gamma }_s$, respectively. Substance users in classes $U_p$ and $U_s$ experience substance-related death at rates given by ${\omega }_p$ and ${\omega }_s$, respectively. Once an individual is in the rehabilitation phase, he/she can successfully recover. Individuals in class $T_p$ recover successfully at a rate ${\delta }_p$ and individuals in class $T_s$ recover successfully at a rate ${\delta }_s$. Individuals experience natural death at a rate $\mu$. We assume that individuals in each compartment are indistinguishable and there is homogeneous mixing so that those at risk of substance use are equally susceptible. In practice, susceptibility to substance use varies. This is due to differences in behavioural, social and environmental factors. The schematic diagram (Figure 1) shows the movement of humans as their status with respect to substance use changes.

Figure 1. Model flow diagram

Combining the parameters, assumptions and the schematic diagram (Figure 1), we obtain the following set of nonlinear ordinary differential equations:

(2.2) $\left\{\begin{array}{c}\frac{dS}{dt}=\mathrm{\Lambda }-\lambda S-\mu S,\\ \frac{d{U}_p}{dt}=\lambda S+{\alpha }_s{U}_s+{\rho }_p{R}_p-(\mu +{\omega }_p+{\gamma }_p+{\sigma }_p+{\alpha }_p)U_p,\\ \frac{d{T}_p}{dt}={\sigma }_p{U}_p-(\mu +{\delta }_p)T_p,\\ \frac{d{R}_p}{dt}={\delta }_p{T}_p+{\gamma }_p{U}_p-(\mu +{\rho }_p)R_p,\\ \frac{d{U}_s}{dt}={\alpha }_p{U}_p+{\rho }_s{R}_s-(\mu +{\omega }_s+{\gamma }_s+{\sigma }_s+{\alpha }_s)U_s,\\ \frac{d{T}_s}{dt}={\sigma }_s{U}_s-(\mu +{\delta }_s)T_s,\\ \frac{d{R}_s}{dt}={\delta }_s{T}_s+{\gamma }_s{U}_s-(\mu +{\rho }_s)R_s,\end{array}\right.$(2.2)

with the initial conditions:

$S(0)=S_0>0,U_p(0)=U_{p0}\ge 0,T_p(0)=T_{p0}\ge 0,R_p(0)=R_{p0}\ge 0,U_s(0)=U_{s0}\ge 0,$

$U_s(0)=U_{s0}\ge 0,T_s(0)=T_{s0}\ge 0,R_s(0)=R_{s0}\ge 0,$

where we assume that all the model parameters are positive.

3. Model analysis

3.1. Positivity of solutions

We now consider the positivity of system (2.2). We prove that all the state variables remain non-negative and the solutions of system (2.2) with positive initial conditions will remain positive for all $t>0$. We thus state the following theorem.

Theorem 3.2. Given that the initial conditions of system (2.2) are $S(0)>0$, $U_p(0)>0$, $T_p(0)>0$, $R_p(0)>0$, $U_s(0)>0$, $T_s(0)>0$ and $R_s(0)>0$. There exists $(S(t),U_p(t),T_p(t),R_p(t),U_s(t),T_s(t),R_s(t)):(0,\mathrm{\infty })\to (0,\mathrm{\infty })$ which solve system (2.2).

Proof. Assume that

$\hat{t}=sup\{t>0:S>0,U_p>0,T_p>0,R_p>0,U_s>0,T_s>0,R_s>0\}\in [0,t].$

Thus $\hat{t}>0$, and it follows from the first equation of system (2.2) that

$\frac{dS}{dt}=\mathrm{\Lambda }-(\lambda +\mu )S\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\lambda \mathrm{g}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{i}\mathrm{n}(2.1).$

Thus, we can solve the above equation with the help of the integrating factor technique for solving first-order linear differential equations to obtain

$\frac{d}{dt}\left[S(t)exp\left\{\mu t+{\int }_0^t\lambda (s)ds\right\}\right]=\mathrm{\Lambda }exp\left[\mu t+{\int }_0^t\lambda (s)ds\right].$

So,

$S(\hat{t})exp\left[\mu \hat{t}+{\int }_0^{\hat{t}}\lambda (s)ds\right]-S(0)={\int }_0^{\hat{t}}\mathrm{\Lambda }exp\left[\mu \hat{t}+{\int }_0^{\hat{t}}\lambda (\nu )d\nu \right]d\hat{t},$

giving

$S(\hat{t})=S(0)exp\left[-\left(\mu \hat{t}+{\int }_0^{\hat{t}}\lambda (s)ds\right)\right]$

$+exp\left[-\left(\mu \hat{t}+{\int }_0^{\hat{t}}\lambda (s)ds\right)\right]\left[{\int }_0^{\hat{t}}\mathrm{\Lambda }exp\left[\mu \hat{t}+{\int }_0^{\hat{t}}\lambda (\nu )d\nu \right]d\hat{t}\right]>0.$

From the second equation of system (2.2), we obtain

$\frac{d{U}_p}{dt}=\lambda S+{\alpha }_s{U}_s+{\rho }_p{R}_p-(\mu +{\omega }_p+{\gamma }_p+{\sigma }_p+{\alpha }_p)U_p,$

$\Rightarrow U_p(\hat{t})\ge U_{p0}{e}^{-(\mu +{\omega }_p+{\gamma }_p+{\sigma }_p+{\alpha }_p)\hat{t}}>0.$

In a similar fashion, it can also be shown that $T_p(t)>0$, $R_p(t)>0$, $U_s(t)>0$, $T_s(t)>0$ and $R_s(t)>0$ for all $t>0$, and this completes the proof.

3.2. Invariant region

It follows from system (2.2) that

$\frac{dN}{dt}=\mathrm{\Lambda }-\mu N-{\omega }_p{U}_p-{\omega }_s{U}_s\le \mathrm{\Lambda }-\mu N.$

Solving leads to

$N(t)\le \frac{\mathrm{\Lambda }}{\mu }+\left(N(0)-\frac{\mathrm{\Lambda }}{\mu }\right)e^{-\mu t}.$

Thus,

$0\le N(t)\le \frac{\mathrm{\Lambda }}{\mu }+\left(N(0)-\frac{\mathrm{\Lambda }}{\mu }\right)e^{-\mu t}=\frac{\mathrm{\Lambda }}{\mu }+N(0)e^{-\mu t}-\frac{\mathrm{\Lambda }}{\mu }{e}^{-\mu t}\le \frac{\mathrm{\Lambda }}{\mu }+N(0)e^{-\mu t}.$

That is, $0\le N(t)\le \frac{\mathrm{\Lambda }}{\mu }+N(0)e^{-\mu t}$ where $N(0)$ represents the value of Equation (2.2(2.2) $\left\{\begin{array}{c}\frac{dS}{dt}=\mathrm{\Lambda }-\lambda S-\mu S,\\ \frac{d{U}_p}{dt}=\lambda S+{\alpha }_s{U}_s+{\rho }_p{R}_p-(\mu +{\omega }_p+{\gamma }_p+{\sigma }_p+{\alpha }_p)U_p,\\ \frac{d{T}_p}{dt}={\sigma }_p{U}_p-(\mu +{\delta }_p)T_p,\\ \frac{d{R}_p}{dt}={\delta }_p{T}_p+{\gamma }_p{U}_p-(\mu +{\rho }_p)R_p,\\ \frac{d{U}_s}{dt}={\alpha }_p{U}_p+{\rho }_s{R}_s-(\mu +{\omega }_s+{\gamma }_s+{\sigma }_s+{\alpha }_s)U_s,\\ \frac{d{T}_s}{dt}={\sigma }_s{U}_s-(\mu +{\delta }_s)T_s,\\ \frac{d{R}_s}{dt}={\delta }_s{T}_s+{\gamma }_s{U}_s-(\mu +{\rho }_s)R_s,\end{array}\right.$(2.2) ) evaluated at the initial values of the respective variables. Then, $\underset{t\to \mathrm{\infty }}{limsup}\mathrm{N}(\mathrm{t})\le \frac{\mathrm{\Lambda }}{\mu }$. Thus, the feasible region for system (2.2) is

(3.1) $\mathrm{\Omega }=\left\{(S,U_p,T_p,R_p,U_s,T_s,R_s)\in {\mathbb{R}}_{+}^7|N\le \frac{\mathrm{\Lambda }}{\mu }\right\}.$(3.1)

It is easy to verify that the region $\mathrm{\Omega }$ is positively invariant with respect to system (2.2).

3.3. Substance-free equilibrium and the initiation reproduction number

The model has a substance-free equilibrium given by:

${\mathcal{D}}^f=(S^f,U_p^f,T_p^f,R_p^f,U_s^f,T_s^f,R_s^f)=\left(\frac{\mathrm{\Lambda }}{\mu },0,0,0,0,0,0\right),$

a scenario depicting a substance-free state in the community or society. The initiation reproduction number ${\mathcal{R}}_{in}$ of the model is defined herein in the substance-using context as the average number of people that each substance user will initiate to substance use during their substance-using career in a population of completely potential substance users. Usually, ${\mathcal{R}}_{in}<1$ implies that substance abuse will die out, whereas ${\mathcal{R}}_{in}>1$ implies that substance abuse will persist within a community and ${\mathcal{R}}_{in}=1$ requires further investigation. The determination of ${\mathcal{R}}_{in}$ is done using the next-generation matrix approach (Van den Driessche & Watmough, 2002). Using this method, the reproduction number of system (2.2) is the spectral radius of the next- generation matrix $F{V}^{-1}$ where

$F=\left(\begin{array}{cccccc}\beta & 0& 0& \beta \eta & 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)\mathrm{a}\mathrm{n}\mathrm{d}V=\left(\begin{array}{cccccc}{g}_{p_1}& 0& -{\rho }_p& -{\alpha }_s& 0& 0\\ -{\sigma }_p& g_{p_2}& 0& 0& 0& 0\\ -{\gamma }_p& -{\delta }_p& g_{p_3}& 0& 0& 0\\ -{\alpha }_p& 0& 0& g_{s_1}& 0& -{\rho }_s\\ 0& 0& 0& -{\sigma }_s& g_{s_2}& 0\\ 0& 0& 0& -{\gamma }_s& -{\delta }_s& g_{s_3}\end{array}\right).$

Thus,

(3.2) $\left\{\begin{array}{c}{\mathcal{R}}_{in}={\mathcal{R}}_p+{\mathcal{R}}_s\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\\ {\mathcal{R}}_p=\left(\frac{\beta }{g_{p_1}}\right)\left(\frac{1-\sum _{i=1}^2{\mathrm{\Phi }}_i}{\left(1-\sum _{i=1}^5{\mathrm{\Psi }}_i\right)+\sum _{i=6}^9{\mathrm{\Psi }}_i}\right)\mathrm{a}\mathrm{n}\mathrm{d}\\ R_s=\left(\frac{\beta \eta {\alpha }_p}{g_{p_1}{g}_{s_1}}\right)\left(\frac{1-\sum _{i=1}^2{\mathrm{\Phi }}_i}{\left(1-\sum _{i=1}^5{\mathrm{\Psi }}_i\right)+\sum _{i=6}^9{\mathrm{\Psi }}_i}\right)\end{array}\right.$(3.2)

with

${\mathrm{\Phi }}_1=\frac{{\gamma }_s{\rho }_s}{g_{s_1}{g}_{s_2}},{\mathrm{\Phi }}_2=\frac{{\delta }_s{\rho }_s{\sigma }_s}{g_{s_1}{g}_{s_2}{g}_{s_3}},{\mathrm{\Psi }}_1=\frac{{\alpha }_p{\alpha }_s}{g_{p_1}{g}_{s_1}},{\mathrm{\Psi }}_2=\frac{{\gamma }_p{\rho }_p}{g_{p_1}{g}_{p_3}},$

${\mathrm{\Psi }}_3=\frac{{\gamma }_s{\rho }_s}{g_{s_1}{g}_{s_3}},{\mathrm{\Psi }}_4=\frac{{\delta }_p{\rho }_p{\sigma }_p}{g_{p_1}{g}_{p_2}{g}_{p_3}},{\mathrm{\Psi }}_5=\frac{{\delta }_s{\rho }_s{\sigma }_s}{g_{s_1}{g}_{s_2}{g}_{s_3}},{\mathrm{\Psi }}_6=\frac{{\gamma }_p{\gamma }_s{\rho }_p{\rho }_s}{g_{p_1}{g}_{p_3}{g}_{s_1}{g}_{s_3}},$

${\mathrm{\Psi }}_7=\frac{{\gamma }_s{\delta }_p{\rho }_p{\rho }_s{\sigma }_p}{g_{p_1}{g}_{p_2}{g}_{p_3}{g}_{s_1}{g}_{s_3}},{\mathrm{\Psi }}_8=\frac{{\gamma }_p{\delta }_s{\rho }_p{\rho }_s{\sigma }_s}{g_{p_1}{g}_{p_3}{g}_{s_1}{g}_{s_2}{g}_{s_3}},{\mathrm{\Psi }}_9=\frac{{\delta }_p{\delta }_s{\rho }_p{\rho }_s{\sigma }_p{\sigma }_s}{g_{p_1}{g}_{p_2}{g}_{p_3}{g}_{s_1}{g}_{s_2}{g}_{s_3}},$

$g_{p_1}=\mu +{\omega }_p+{\gamma }_p+{\sigma }_p+{\alpha }_p,g_{p_2}=\mu +{\gamma }_p,g_{p_3}=\mu +{\rho }_p,$

$g_{s_1}=\mu +{\omega }_s+{\gamma }_s+{\sigma }_s+{\alpha }_s,g_{s_2}=\mu +{\gamma }_s\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{g}}_{{\mathrm{s}}_3}=\mu +{\rho }_{\mathrm{s}}.$

We can clearly note that

$\sum _{i=1}^2\text{,}{\mathrm{\Phi }}_i\le 1\text{,}and\sum _{i=1}^5\text{,}{\mathrm{\Psi }}_i\le 1.$

Therefore, ${\mathcal{R}}_{in}$ is non-negative. Here, ${\mathcal{R}}_{in}$ is the sum of two sub-reproduction numbers representing the contributions of individuals in compartments $U_p$ and $U_s$, respectively.

We now show that substance abuse can be eliminated from the community (when ${\mathcal{R}}_{in}<1$) if the initial sizes of the sub-populations of the model are in the basin of attraction of the substance-free equilibrium. We thus have the following theorem (Theorem 3.5).

Theorem 3.5. The substance-free equilibrium ${\mathcal{D}}^f$ is locally asymptotically stable when ${\mathcal{R}}_{in}<1$ and is unstable when ${\mathcal{R}}_{in}>1$.

Proof. The Jacobian matrix of system (2.2) at ${\mathcal{D}}_f$ is given by

$J({\mathcal{D}}^f)=\left(\begin{array}{ccccccc}-\mu & -\beta & 0& 0& -\beta \eta & 0& 0\\ 0& \beta -g_{p_1}& 0& {\rho }_p& \beta \eta +{\alpha }_s& 0& 0\\ 0& {\sigma }_p& -g_{p_2}& 0& 0& 0& 0\\ 0& {\gamma }_p& {\delta }_p& -g_{p_3}& 0& 0& 0\\ 0& {\alpha }_p& 0& 0& -g_{s_1}& 0& {\rho }_s\\ 0& 0& 0& 0& {\sigma }_s& -g_{s_2}& 0\\ 0& 0& 0& 0& {\gamma }_s& {\delta }_s& -g_{s_3}\end{array}\right),$

where $g_{p_1}$, $g_{p_2}$, $g_{p_3}$, $g_{s_1}$, $g_{s_2}$ and $g_{s_3}$ are defined as before. The local stability of the substance-free equilibrium is determined by the following submatrix of $J({\mathcal{D}}^f)$,

$\bar{J}({\mathcal{D}}^f)=\left(\begin{array}{cccccc}\beta -g_{p_1}& 0& {\rho }_p& \beta \eta +{\alpha }_s& 0& 0\\ {\sigma }_p& -g_{p_2}& 0& 0& 0& 0\\ {\gamma }_p& {\delta }_p& -g_{p_3}& 0& 0& 0\\ {\alpha }_p& 0& 0& -g_{s_1}& 0& {\rho }_s\\ 0& 0& 0& {\sigma }_s& -g_{s_2}& 0\\ 0& 0& 0& {\gamma }_s& {\delta }_s& -g_{s_3}\end{array}\right).$

Since all off-diagonal elements are positive, we now consider matrix $-\bar{J}({\mathcal{D}}^f)$. We claim that $-\bar{J}({\mathcal{D}}^f)$ is an $M$-matrix. Multiplying matrix $-\bar{J}({\mathcal{D}}^f)$ by the positive $6\times 1$ matrix

$W_1={\left[g_{p_2},{\sigma }_p,\frac{g_{p_2}{\gamma }_p+{\delta }_p{\sigma }_p}{g_{p_3}},\frac{g_{p_2}{\alpha }_p}{g_{s_1}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)},\frac{g_{p_2}{\alpha }_p{\sigma }_s}{g_{s_1}{g}_{s_2}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)},\frac{g_{p_2}{\alpha }_p\left(g_{s_2}{\gamma }_s+{\delta }_s{\sigma }_s\right)}{g_{s_1}{g}_{s_2}{g}_{s_3}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)}\right]}^T$

we have

$-\bar{J}(D^f)W_1=(1-R_{in})\cdot W_2,$

where $W_2$ is a positive $6\times 1$ matrix given by

$W_2={\left[\frac{\left(1-\sum _{i=1}^5{\mathrm{\Psi }}_i\right)+\sum _{i=6}^9{\mathrm{\Psi }}_i}{g_{p_3}{g}_{s_1}{g}_{s_2}{g}_{s_3}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)},0,0,0,0,0\right]}^T.$

Then, it follows from $M$-matrix theory that all eigenvalues of $\bar{J}({\mathcal{D}}^f)$ have negative real parts, which implies the local asymptotic stability of the substance-free equilibrium if ${\mathcal{R}}_{in}<1$. On the other hand, it can be shown that the determinant of $\bar{J}({\mathcal{D}}^f)$ is given by

$\mathrm{d}\mathrm{e}\mathrm{t}\left(\bar{J}({\mathcal{D}}^f)\right)=\left[\left(1-\sum _{i=1}^5{\mathrm{\Psi }}_i\right)+\sum _{i=6}^9{\mathrm{\Psi }}_i\right]({\mathcal{R}}_{in}-1).$

Thus, if ${\mathcal{R}}_{in}<1$, then matrix $\bar{J}({\mathcal{D}}^f)$ has eigenvalues with negative real parts, which implies the stability of the substance-free equilibrium. This completes the proof.

3.4. Substance-persistent equilibrium

The substance-persistent equilibrium ${\mathcal{D}}^{\ast }=\left(S^{\ast },U_p^{\ast },T_p^{\ast },R_p^{\ast },U_s^{\ast },T_s^{\ast },R_s^{\ast }\right)$ always satisfies

(3.3) $\left\{\begin{array}{c}0=\mathrm{\Lambda }-{\lambda }^{\ast }{S}^{\ast }-\mu S^{\ast },\\ 0={\lambda }^{\ast }{S}^{\ast }+{\alpha }_s{U}_s^{\ast }+{\rho }_p{R}_p^{\ast }-g_{p_1}{U}_p^{\ast },\\ 0={\sigma }_p{U}_p^{\ast }-g_{p_2}{T}_p^{\ast },\\ 0={\delta }_p{T}_p^{\ast }+{\gamma }_p{U}_p^{\ast }-g_{p_3}{R}_p^{\ast },\\ 0={\alpha }_p{U}_p^{\ast }+{\rho }_s{R}_s^{\ast }-g_{s_1}{U}_s^{\ast },\\ 0={\sigma }_s{U}_s^{\ast }-g_{s_2}{T}_s^{\ast },\\ 0={\delta }_s{T}_s^{\ast }+{\gamma }_s{U}_s^{\ast }-g_{s_3}{R}_s^{\ast }.\end{array}\right.$(3.3)

From the third, fourth, fifth, sixth and seventh equation of (3.3), we get that

(3.4) $\left\{\begin{array}{c}{T}_p^{\ast }=\frac{{\sigma }_p{U}_p^{\ast }}{g_{p_2}},R_p^{\ast }=\frac{U_p^{\ast }\left(g_{p_2}{\gamma }_p+{\delta }_p{\sigma }_p\right)}{g_{p_2}{g}_{p_3}},\\ U_s^{\ast }=\frac{g_{s_2}{g}_{s_3}{\alpha }_p{U}_p^{\ast }}{g_{s_1}{g}_{s_2}{g}_{s_3}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)},\\ T_s^{\ast }=\frac{g_{s_3}{\alpha }_p{U}_p^{\ast }{\sigma }_s}{g_{s_1}{g}_{s_2}{g}_{s_3}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)},\\ R_s^{\ast }=\frac{{\alpha }_p{U}_p^{\ast }\left(g_{s_2}{\gamma }_s+{\delta }_s{\sigma }_s\right)}{g_{s_1}{g}_{s_2}{g}_{s_3}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)}.\end{array}\right.$(3.4)

Adding the first and second equations of (3.3) leads to

(3.5) $0=\mathrm{\Lambda }-\mu S^{\ast }+{\alpha }_s{U}_s^{\ast }+{\rho }_p{R}_p^{\ast }-g_{p_1}{U}_p^{\ast }.$(3.5)

Using (3.4) to solve (3.5) for $S^{\ast }$ in terms of $U_p^{\ast }$ leads to

(3.6) $S^{*}=\frac{\left(1-{\displaystyle \sum _{i=1}^2{\Phi }_i}\right)-\left(\left(1-{\displaystyle \sum _{i=1}^5{\Psi }_i}\right)+{\displaystyle \sum _{i=6}^9{\Psi }_i}\right)g_{p_1}{U}_p^{*}}{\left(1-{\displaystyle \sum _{i=1}^2{\Phi }_i}\right)}.$(3.6)

Substituting (3.4) and (3.6) into the first equation of (3.3) and solving for $U_p^{\ast }$ leads to $U_p^{\ast }=0$ which corresponds to the substance-free equilibrium or

(3.7) $U_p^{\ast }=\frac{g_{p_2}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)}{g_{s_2}{g}_{s_3}\left[\beta g_{s_1}\left(1-\sum _{i=1}^2{\mathrm{\Phi }}_i\right)+\beta \eta {\alpha }_p\right]}\left({\mathcal{R}}_{in}-1\right).$(3.7)

Substituting (3.7) into expressions (3.4) and (3.6) gives the substance-persistent equilibrium ${\mathcal{D}}^{\ast }=\left(S^{\ast },U_p^{\ast },T_p^{\ast },R_p^{\ast },U_s^{\ast },T_s^{\ast },R_s^{\ast }\right)$. We thus have the following result on the existence of the substance-persistent equilibrium.

Theorem 3.7. System (2.2) has a positive substance-persistent equilibrium if and only if ${\mathcal{R}}_{in}>1$.

4. Numerical simulations

4.1. A case study for the Western Cape Province of South Africa

We present an application of our model to a case study of the Western Cape Province of South Africa. We focus mainly on the trends of polysubstance use. Data collected by the Medical Research Council (MRC) from specialists at substance abuse treatment centres in Cape Town (as part of the SACENDU, 2019) show that poly-substance use remains high over the recent years. Based on data from SACENDU, patients who seek treatment for substance abuse-related problems are categorized as either primary substance abusers or those abusing two or more substances. However, a few reports have indicated the number of substances abused by patients seeking treatment with some patients indicating that they have abused up to four substances. Since most SACENDU reports over the years have combined those abusing more than one substance into one class of individuals, the model in this paper considers such a setting for purposes of illustrating trends of poly-substance abuse. The data used are for the period covering the first half of $2001$ up to the first half of $2018$.

4.1.1. Parameter estimation

We make use of curve fitting in Matlab to estimate model parameters used in our numerical simulations. Curve fitting is a process that allows us to quantitatively estimate the trend of the outcomes where equations of approximating curves are fit to the raw field data (Mushanyu, Nyabadza, & Stewart, 2015). Nevertheless, for a given set of data, the fitting curves of a given type are generally not unique. Thus, a curve with a minimal deviation from all data points is desired. This best-fitting curve can be obtained by the method of least squares (Mushanyu, Nyabadza, & Stewart, 2015). The least-squares curve fit routine (lsqcurvefit) in Matlab with optimization is used to estimate the parameter values. Many parameters are known to lie within some intervals. During the estimation of parameter values, unknown parameter values are given a lower and upper bound from which the set of parameter values that provide the best fit are obtained. The parameters obtained from curve fitting and intervals used are given in Table 2 where data for primary substance abusers were used. We present an application of system (2.2) through fitting the model to data from the MRC’s SACENDU project. We fit the model system (2.2) to data on poly-substance use in Cape Town of South Africa. Also, we estimate the incidence of primary substance abusers using data on poly-substance use in the Western Cape Province of South Africa. Our estimated incidence of primary substance abuse is evaluated using the initiation function $\lambda S$ given in (2.1). We use data in Table 1 which was collected by SACENDU (2019).

Table 1. Polysubstance use in the Western Cape Province of South Africa for the period $2001$a to $2018$a $(\mathrm{\%})$. Letter “a” represents the first 6 months of the year and “b” represents the last 6 months of the year

Year	2001a	2001b	2002a	2002b	2003a	2003b	2004a	2004b
Primary $(\mathrm{\%})$	$58$	$60$	$69$	$61$	$53$	$52$	$51$	$52$
Primary$+2^{\mathbf{n}\mathbf{d}}$$(\mathrm{\%})$	$42$	$40$	$31$	$39$	$47$	$48$	$49$	$48$
Year	$2005$a	$2005$b	$2006$a	$2006$b	$2007$a	$2007$b	$2008$a	$2008$b
Primary $(\mathrm{\%})$	$51$	$51$	$55$	$52$	$56$	$56$	$54$	$57$
Primary$+2^{\mathbf{n}\mathbf{d}}$$(\mathrm{\%})$	$49$	$49$	$45$	$48$	$44$	$44$	$46$	$43$
Year	$2009$a	$2009$b	$2010$a	$2010$b	$2011$a	$2011$b	$2012$a	$2012$b
Primary $(\mathrm{\%})$	$58$	$54$	$57$	$56$	$61$	$53$	$57$	$51$
Primary$+2^{\mathbf{n}\mathbf{d}}$$(\mathrm{\%})$	$42$	$46$	$43$	$44$	$39$	$47$	$43$	$49$
Year	$2013$a	$2013$b	$2014$a	$2014$b	$2015$a	$2015$b	$2016$a	$2016$b
Primary$(\mathrm{\%})$	$46$	$52$	$51$	$56$	$49$	$51$	$52$	$52$
Primary$+2^{\mathbf{n}\mathbf{d}}$$(\mathrm{\%})$	$54$	$48$	$49$	$44$	$51$	$49$	$48$	$48$
Year	$2017$a	$2017$b	$2018$a
Primary$(\mathrm{\%})$	$54$	$56$	$51$
Primary$+2^{\mathbf{n}\mathbf{d}}$$(\mathrm{\%})$	$46$	$44$	$49$

Table 2. Parameter values and ranges obtained from data fitting using data for primary substance abusers in Cape Town

Parameter	Range	Value	Source
$\eta$	$1-10$	$2.7207$	Estimated
$\beta$	$0-1$	$4.7781\times {10}^{-5}$	Estimated
${\alpha }_p$	$0-1$	$0.6410$	Estimated
${\sigma }_p$	$0-1$	$0.0025$	Estimated
${\gamma }_p$	$0-1$	$0.0171$	Estimated
${\omega }_p$	$0-1$	$0.8100$	Estimated
${\delta }_p$	$0-1$	$0.1939$	Estimated
${\rho }_p$	$0-1$	$0.1318$	Estimated
${\alpha }_s$	$0-1$	$0.2099$	Estimated
${\sigma }_s$	$0-1$	$0.8518$	Estimated
${\gamma }_s$	$0-1$	$0.7694$	Estimated
${\omega }_s$	$0-1$	$0.5479$	Estimated
${\delta }_s$	$0-1$	$0.1733$	Estimated
${\rho }_s$	$0-1$	$0.1613$	Estimated
$\mathrm{\Lambda }$	$0.01-0.04$	$0.0257$	Estimated
$\mu$	$0.01-0.03$	$0.0221$	Estimated

4.1.2. Results

Figure 2(a) shows the trend in the proportion of primary substance abusers in treatment centres of Cape Town. The model fits well with the data from Table 1. The results show that the proportion of primary substance abusers in Cape Town shall continue to steadily decrease over the years. However, the proportion of primary substance abusers continues to dominate that of secondary substance abusers. Parameter values estimated using data for primary substance abusers in Cape Town are shown in Table 2. This estimation assumes that the dynamics remain the same over the entire period. Figure 2(b) shows the estimated incidence of primary substance abusers using data for primary substance abusers in the Western Cape Province of South Africa. We note that there is a steady increase in the incidence of primary substance abuse between the years 2001 and 2003 after which the incidence begins to decline up to the year 2010. After the year 2010, we observe some slight increase in the incidence of primary substance abusers up to the year 2018. Figure 3(a) shows the trend in the proportion of primary substance abusers in treatment centres of Cape Town projected up to the year 2030. We note from the projected trend that the proportion of primary substance abusers will continue to decrease. Figure 3(b) shows the projected incidence of primary substance abuse in Cape Town up to the year 2030. We observe that the incidence of primary substance abuse shall continue to steadily increase up to the year 2030.

Figure 2. (a) System (2.2) fitted to data for primary substance abusers in Cape Town. The blue circles indicate the actual data and the solid red line indicates the model fit to the data. (b) Estimated incidence of primary substance abusers in Cape Town

Figure 3. (a) System (2.2) fitted to data for primary substance abusers in Cape Town and projected up to the year 2030. The blue circles indicate the actual data and the solid red line indicates the model fit to the data. (b) Estimated incidence of primary substance abusers in Cape Town and projected up to the year 2030

Figure 4. System (2.2) fitted to data for secondary substance abusers $(\mathit{a})$ and projected up to the year $\mathbf{2030}$ $(\text{b})$ in Cape Town. The blue circles indicate the actual data and the solid red line indicates the model fit to the data

Figure 4 $(\mathrm{a})$ shows the trend in the proportion of secondary substance abusers in treatment centres of Cape Town. Figure 4 $(\mathrm{a})$ shows that the model fits well with the data from Table 1. The results show that the proportion of secondary substance abusers in Cape Town shall steadily increase over the years. Figure 4 $(\mathrm{b})$ shows the projected trends in the proportion of secondary substance abusers in treatment centres of Cape Town up to the year $2030$. The figure shows that the projected proportion of secondary substance abusers will not differ much with the current observed proportion.

4.2. Numerical results

We carry out detailed numerical simulations using Matlab to support our theoretical findings. Numerical solutions of a model depend on the values of all its parameters. The initial conditions used are: $S(0)=0.95$, $U_p(0)=0.03$, $T_p(0)=0$, $R_p(0)=0$, $U_s(0)=0.02$, $T_s(0)=0$, $R_s(0)=0$.

Figure 5. General population trajectory of individuals in class $\mathit{S}(t)$, (a); class ${\mathit{U}}_{\mathit{p}}(t)$, (b) and class ${\mathit{U}}_{\mathit{s}}(t)$, (c)

Figure 6. General population trajectory of individuals in class ${\mathit{T}}_{\mathit{p}}(t)$, (a); class ${\mathit{T}}_{\mathit{s}}(t)$, (b); class ${\mathit{R}}_{\mathit{p}}(t)$, (c) and class ${\mathit{R}}_{\mathit{s}}(t)$, (d)

Figure 7. Effects of varying ${\mathit{\sigma }}_{\mathit{p}}$, $\mathbf{(}a)$, ${\mathit{\sigma }}_{\mathit{s}}$, $\mathbf{(}b)$, ${\mathit{\alpha }}_{\mathit{p}}$, $(c)$ and ${\mathit{\alpha }}_{\mathit{s}}$, $\mathbf{(}d)$ on the prevalence of substance abuse, starting from $\mathbf{0.1}$ up to $\mathbf{1.0}$ with a step size of $\mathbf{0.01}$ across all the parameters

Figure 5 and Figure 6 illustrate the general population trajectory of individuals in classes S(t), U_p(t), T_p(t), R_p(t), U_s(t), T_s(t) and R_s(t). Figure 7 illustrates the effect of varying parameters ${\sigma }_p$, ${\sigma }_s$, ${\alpha }_p$ and ${\alpha }_s$ on the prevalence of substance abuse. Figure 7 $(\mathrm{a})$ shows that increasing the parameter ${\sigma }_p$ leads to a decrease in the prevalence of both primary and secondary substance abusers. Also, Figure 7 $(\mathrm{b})$ shows that increasing the parameter ${\sigma }_s$ leads to a decrease in the prevalence of both primary and secondary substance abusers. However, we observe that the parameter ${\sigma }_s$ has more positive impact in reducing the prevalence of both primary and secondary substance abusers as compared to the parameter ${\sigma }_p$. Figure 7 $(\mathrm{c})$ shows that increasing the parameter ${\alpha }_p$ results in a decrease in the prevalence of primary substance abusers and an increase in the prevalence of secondary substance abusers. Lastly, Figure 7 $(\mathrm{d})$ shows that increasing the parameter ${\alpha }_s$ results in an increase in the prevalence of primary substance abusers and a decrease in the prevalence of secondary substance abusers. We note from Figure 7 that although the prevalence of substance abuse either increases or decreases when a parameter is varied, the prevalence change for secondary substance abuse always dominates that of primary substance abuse. Thus, more effort should be directed towards ensuring adequate treatment of secondary substance abusers as this will substantially contribute to the control of substance abuse in the general population.

5. Conclusion

The goal of this paper is to model the trends of primary and secondary substance abuse in the Western Cape Province of South Africa. As such, the model developed in this paper explicitly captures the levels of substance abuse by classifying substance abuse as either primary substance abuse or secondary substance abuse, not considered in past substance abuse models. Structuring the model in such a way enables the application of the model to data from SACENDU.

An epidemic threshold termed the initiation reproduction number, ${\mathcal{R}}_{in}$, has been derived and qualitatively used to analyze the model equilibria. The substance-free equilibrium point is shown to be locally asymptotically stable whenever the initiation reproduction number ${\mathcal{R}}_{in}$ is less than unity. It was also shown that the system has a positive unique substance-persistent equilibrium if and only if ${\mathcal{R}}_{in}>1$. The least-squares curve fit routine (lsqcurvefit) in Matlab with optimization has been used to fit the model to data on primary and secondary substance abusers in treatment centres of Cape Town of South Africa. The model was observed to fit well with these data and parameter values from model fitting were obtained. Performing a visual predictive check on the graphs from data fitting suggests that the proportion of primary substance abusers in Cape Town shall continue to decrease over the years whereas that of secondary substance abusers shall continue to increase but steadily. However, the proportion of primary substance abusers shall continue to dominate that of secondary substance abusers.

Numerical simulations have also been performed to support theoretical results. For illustrative purposes, we have chosen a few important parameters, namely, ${\sigma }_p$, ${\sigma }_s$, ${\alpha }_p$ and ${\alpha }_s$. The results obtained suggest that there is need to increase control efforts targeting secondary substance abusers as this will have more impact in the fight against substance abuse. This, however, should not be taken to mean that control measures should do away with primary substance abusers, rather a concerted effort on both substance abusing classes is perceived to achieve enhanced positive results.

The study presented here is not exhaustive. We want to acknowledge that restricting our model to the primary and secondary substance abuse levels only has some implications in limiting the results of this study. Thus, depending on the availability of data, some future work on the model can be done by extending it to include more levels of substance abuse, for instance, primary substance abusers, primary $+2\text{nd}$ substance abusers, primary $+2\text{nd}$ substance abusers $+3\text{rd}$ substance abusers and so on.

Acknowledgements

The authors acknowledge, with thanks, the support of their respective departments for the production of this manuscript.

Additional information

Funding

The authors received no direct funding for this research.

Notes on contributors

Josiah Mushanyu

Josiah Mushanyu (PhD) is a Senior Lecturer in the Department of Mathematics, Faculty of Science, University of Zimbabwe. F. Nyabadza (PhD) is a Professor in Applied Mathematics and Deputy HOD, Department of Mathematics and Applied Mathematics, Faculty of Science, University of Johannesburg. P. Mafuta (PhD) is a Senior Lecturer in the Department of Mathematics, Faculty of Science, University of Zimbabwe. The authors’ research interests include Mathematical Modelling and Application of Graph Theoretic Approaches to Mathematical Biology.