Modelling the trends of inpatient and outpatient rehabilitation for methamphetamine in the Western Cape province of South Africa

Background Dependence on methamphetamine remains one of the major health and social problem in the Western Cape province of South Africa. We consider a mathematical model that takes into account two forms of rehabilitation, namely; inpatient and outpatient. We examine the trends of these two types of rehabilitation. We also seek to investigate the global dynamics of the developed methamphetamine epidemic model. Methods The model is designed by likening the initiation process to an infection that spreads in a community through interactions between methamphetamine users and non-users. We make use of Lyapunov functions obtained from a suitable combination of common quadratic and Volterra-type functions to establish the global stability of the methamphetamine-persistent steady state. The least squares curve fit routine (lsqcurvefit) in Matlab with optimization is used to estimate the parameter values. Results The model analysis shows that the model has two equilibria, the methamphetamine free equilibrium and the methamphetamine persistent equilibrium, that are both globally stable when the threshold \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_a<1$$\end{document}Ra<1 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_a>1$$\end{document}Ra>1, respectively. Upon fitting the model to data on drug users under rehabilitation, parameter values that give the best fit were obtained. The projections carried out the long term trends of these forms of rehabilitation. Conclusion The results suggest that inpatient rehabilitation programs have an increased potential of enhancing the chances of recovery for methamphetamine addicts.


Background
The methamphetamine abuse problem and its consequences to communities in the Western Cape province present a complex scenario that drives sexually transmitted infections morbidity, mortality and heavy budgetary constraints. For example in South Africa over R20 billion is used annually, in treating drug users, cracking down drug traffickers and in prevention and media campaigns [1]. In the Western Cape Province of South Africa, the rising demand for substance abuse services and calls by communities for additional services has led to the Provincial Department of Social Development allocating additional resources to the prevention and treatment of substance use disorders (SUD's). However, planning around the allocation of these resources has been hampered by several informational barriers. Decision makers within the Western Cape Department of Social Development (DOSD) do not have adequate information on the nature and extent of substance abuse in the province, the extent to which there is unmet substance abuse service needs, where these unmet needs are greatest, and which population subgroups have relatively high unmet needs [2]. There is therefore need for extensive research into methamphetamine abuse trends for the development of public policies and focused correctional services to decrease the gang populations.
Recently, many researchers have drawn a lot of scholarship from infectious diseases modelling with the aim of describing drug abuse spreading like an infectious disease [3][4][5]. Typical examples are models designed for heroin epidemics [12][13][14][15], peer influence [6][7][8][9] and the spread of alcoholism [10,11]. Mathematical models can contribute to the understanding of the various aspects of methamphetamine use and are very useful in determining how prevalent drug use is. They can even help in designing and choosing proper interventions by providing a means of integrating data from different sources, describing a process to increase understanding and simulating policy experiments that are cumbersome in real life [16][17][18].
There are two primary forms of rehabilitation for methamphetamine addiction namely; inpatient and outpatient rehabilitation [19]. Rehabilitation helps individuals to improve their function, mobility, independence and quality of life. It helps individuals live fully regardless of impairment. It helps people who are living with various health conditions to maintain the functioning they have [20]. Inpatient methamphetamine rehabilitation provides an addict a place to live and at the same time providing 24 h support and treatment. Outpatient methamphetamine rehabilitation allows users to reside at home and come to the treatment center on a regular basis [19] (Figs. 1, 2).
Mathematical models on drug abuse that have been developed so far, have unified rehabilitation, see for instance [14,15,21,31]. In this paper, we develop a mathematical model that takes into account both inpatient and outpatient forms of rehabilitation. Our aim is to examine the trends of these two types of rehabilitation. We also seek to investigate the global dynamics of the developed methamphetamine epidemic model. We analyze the global stability of the model equilibria and fit the model to data on these two forms of rehabilitation.
The paper is arranged as follows; in Sect. "Methods", we formulate and establish the basic properties of the model. The model is analysed for stability in Sect. "Model analysis". Parameter estimation and projection graphs are given in Sect. "Results and discussion". Numerical results are also presented in this section. The paper is concluded in Sect. "Conclusions" (Figs. 3, 4).

Model formulation
The dynamics of inpatient and outpatient rehabilitation are modelled by considering the human population divided into four distinct compartments. The . This decline, might be a result of retarding economic situation or expensive costs associated with inpatient rehabilitation centers. Also, due to the fact that inpatient rehabilitation programs are formally 24 h programs, some people would prefer outpatient rehabilitation programs to allow them to work for their families whilst receiving treatment. It demonstrates a good fit to the data from Table 1 compartments comprise of S denoting the population at risk of being initiated into methamphetamine abuse, U those initiated into methamphetamine abuse, T 0 those in rehabilitation as out-patients and T i those in rehabilitation as in-patients. Here, all susceptible individuals range from age 15-64 years. The total human population is thus given by We assume that individuals in each compartment are indistinguishable and there is homogeneous mixing. Susceptible humans enter the population through births or immigration at a rate . Susceptible individuals are initiated into methamphetamine use following interaction with individuals using drugs. We assume an initiation function that is analogous to the force of infection for epidemic models. Thus, the per capita contact rate β is a product of the effective number of contacts c, between methamphetamine users and the susceptible population, and the probability β , that a contact results into initiation into methamphetamine use, that is, β = cβ. A fraction U/N of the contacts is with those drug users not in rehabilitation and the average number of potential contacts resulting in susceptible individuals becoming methamphetamine abusers is βU /N. Also, a fraction T 0 /N of the contacts is with individuals under outpatient rehabilitation. The average number of potential contacts of each susceptible individual with individuals in outpatient rehabilitation is βηT 0 /N. The parameter η measures the relative ability for individuals in outpatient rehabilitation to initiate new methamphetamine users. Assuming that the rate at which individuals in outpatient rehabilitation recruit initiates is lower than that for drug users not in rehabilitation, we have, 0 < η < 1. This is due to the fact that individuals in outpatient rehabilitation will be receiving some form of counselling and therapy to aid them in quitting drug abuse. Such counselling may lead them to discourage susceptible individuals into becoming methamphetamine abusers. This done by raising awareness on the dangers associated with problematic drug use. Here, we make an assumption that individuals under inpatient rehabilitation cannot produce new initiates due to the fact that, they do not come in contact with the general population during the rehabilitation process. The total number of relevant contacts gives the initiation function, Upon being initiated into methamphetamine use, a susceptible individual moves into the compartment U, of methamphetamine abusers. Upon realizing the demoralizing consequences related to methamphetamine abuse, individuals in the compartment U seek help via rehabilitation programs. We consider two primary types of rehabilitation programs, namely; inpatient and outpatient forms of rehabilitation. Here, we define the effectiveness of these rehabilitation programs to be, "the measure of the benefits and changes in the functioning of an individual acquired during the period when he/ she was under treatment". The effectiveness of rehabilitation programs is measured by the parameter ε. If ε = 0, then the rehabilitation programs are not effective,   Table 1 ε = 1 corresponds to completely effective rehabilitation programs, while 0 < ε < 1 implies that rehabilitation programs will be effective to some degree. Amongst individuals in compartment U who are seeking help through rehabilitation, we have that a proportion p of these individuals are recruited into inpatient rehabilitation facilities and the complementary proportion (1 − p) are recruited into outpatient rehabilitation. The rate at which methamphetamine users are recruited into rehabilitation (inpatient or outpatient) is given by σ. Once an individual is in the rehabilitation phase, he/she can either quit permanently, relapse into methamphetamine use or die. Individuals experience natural death at a rate µ. The rate at which individuals under outpatient rehabilitation quit methamphetamine abuse permanently is represented by δ 0 and the rate at which those under inpatient rehabilitation quit methamphetamine abuse permanently is represented by δ 1 . Individuals undergoing outpatient rehabilitation relapse into methamphetamine use at a rate ρ 1 and individuals in inpatient rehabilitation facilities relapse into methamphetamine use at a rate ρ 2 . Substance abusers undergoing outpatient rehabilitation continue to keep contact with people and circumstances that trigger addiction whereas substance abusers undergoing inpatient rehabilitation are completely immersed in the program and separated from the lifestyle and habits that supported drug use, thus, relapsing may be more likely for outpatient rehabilitants as compared to inpatient rehabilitants. We can safely assume that ρ 2 < ρ 1 . Individuals under inpatient rehabilitation programs can move to outpatient rehabilitation programs at a rate γ 2 . This movement may be as a result of some individuals failing to cope up with the higher costs usually associated with inpatient rehabilitation and thereby forcing them rather to seek help through outpatient rehabilitation which is generally cheaper. This movement might also be as a result of noted recovery to some individuals, who if discharged can still fully recover whilst residing at home. The rate at which individuals under outpatient rehabilitation move to inpatient rehabilitation facilities is given by γ 1 . This movement can be due to the fact that some individuals under outpatient rehabilitation might want to take full advantage of all the resources available in treatment, such as personal therapy, group therapy, educational classes on addiction as well as job skills and other related support structures in order to quicken their recovery process. So this in turn will lead them to seek help from inpatient rehabilitation facilities.
The description of the model formulation and the model diagram lead to the following set of nonlinear ordinary differential equations together with the initial conditions (Figs. 5, 6): Time in years Ratio of inpatients to outpatients Fig. 6 Estimation of the proportions of inpatient rehabilitants relative to outpatient rehabilitants in Cape Town projected over 5 more years. It demonstrates the changes of the ratio of inpatient rehabilitants to outpatient rehabilitants over the modelling time and projects these changes over a period of 5 years. The graph shows a steady decrease in the proportion inpatient rehabilitants coupled with an increase in the proportion of outpatient rehabilitants. The reason for the decrease is similar to the one used to for Figs. 2 and 3

Positivity of solutions
We now consider the positivity of model system Eqs.
(2)- (5). We prove that all the state variables remain nonnegative and the solutions of model system Eqs. (2)-(5) with positive initial conditions will remain positive for all t > 0. We thus state the following theorem (Figs. 7, 8).
Proof Assume that Thus t > 0, and it follows from the first equation of model system Eqs. (2)-(5) that From the second equation of model system Eqs.
(2)- (5), we obtain Similarly it can also be shown that T 0 (t) > 0 and T i (t) > 0 for all t > 0, and this completes the proof.

Theorem 2
The feasible region G defined by with initial conditions S 0 ≥ 0, U 0 ≥ 0, T 00 ≥ 0 and T i0 ≥ 0 is positively invariant and attracting with respect to model system Eqs.
Proof Summing up Eqs. (2)-(5), we obtain that the total population satisfies the differential equation . This projected decrease in the proportion of inpatient rehabilitants can be an attribute of higher costs usually associated with inpatient rehabilitation programs and as a result, few individuals will afford receiving treatment at those facilities. As can be seen from the incidence curves related to data on inpatient rehabilitants, that inpatient rehabilitation programs have an increased potential of reducing methamphetamine incidence, there are fears that this decrease in the proportion of inpatient rehabilitants will to a greater extend result to higher prevalence rates of use observed for a long period of time  5), with (S 0 , U 0 , T 00 , T i0 ) ∈ R 4 + remain in G for all t > 0. Therefore the ω − limit set of solutions of Eqs. (2)-(5) in G are contained in G. The uniqueness, existence and continuity results hold for Eqs. (2)-(5). The system is thus mathematically and epidemiologically well-posed, see also [4]. Our analysis will be based on the dynamics of the solutions generated in G.

Methamphetamine-free equilibrium and the abuse reproduction number
In this section, we carry out stability analysis of the methamphetamine-free equilibrium. The model has a methamphetamine-free equilibrium given by a scenario depicting a methamphetamine-free state in the community or society. The abuse reproduction number of the model, denoted R a , is the average number of secondary cases generated by one drug user during his/her duration of drug use in a population of completely potential drug users. The determination of R a through the next generation matrix [23] method has been explored in many papers (see for instance [24][25][26][27][28][29]). Using the same method we have where Here, R a is the sum of three sub-reproduction numbers representing the contributions of individuals in compartments U, T 0 and T i respectively. The term 1 is the proportion of individuals who move back and forth between compartments T 0 and T i , 2 is the proportion of individuals who move back and forth between compartments T 0 and U and 3 is the proportion of individuals who move back and forth between compartments T i and U. Interestingly, 4 is the proportion of those individuals who were once inpatient rehabilitants and have relapsed after joining outpatient rehabilitation and lastly, 5 is the proportion of those individuals who were once outpatient rehabilitants and have relapsed after joining inpatient rehabilitation facilities. The terms 4 and 5 play an essential role in evaluating the effectiveness of inpatient and outpatient rehabilitation programs.
The following result follows from [23].

Theorem 3
The methamphetamine-free equilibrium point G 0 of model Eqs.
, Estimated incidence of methamphetamine abuse using data for inpatient rehabilitants in Cape Town and projected for 5 more years. We also note the decreasing trend for incidence of methamphetamine abuse. Our estimated incidence will have decreased to below 7 % by the year 2018

Global stability of the methamphetamine-free steady state
We shall now prove the global stability of the methamphetamine-free equilibrium point G 0 whenever the reproduction number is less than unity.

Theorem 4
The methamphetamine-free equilibrium point G 0 of model Eqs.
Proof Let V (U, T 0 , T i ) = aU + bT 0 + cT i be a candidate Lyapunov function for some non-negative constants a, b and c. The time derivative of V is given by We now evaluate the coefficients of the suitable Lyapunov function such that the coefficients of T 0 and T i are equal to zero. We thus obtain Using these coefficients, the time derivative of the Lyapunov function can be expressed as is the singleton G 0 . By Lasalle invariance principle [30], this implies that G 0 is globally stable in G if R a ≤ 1.
We observe that the Jacobian matrix evaluated at G 0 has a positive eigenvalue whenever R a > 1. Therefore the methamphetamine-free equilibrium is unstable if R a > 1 . This completes the proof.

The methamphetamine persistent equilibrium point
In this section we determine the methamphetamine persistent equilibrium point by first setting the left hand side of model Eqs.
The results on persistence of the model are given in [31]. For more information about persistence of dynamical systems, we refer the reader to [32][33][34][35][36].

Global stability of the methamphetamine-persistent steady state
In this section, we prove the global stability of the methamphetamine-persistent steady state.

Theorem
6 Assume that . If R a > 1, then the unique methamphetamine-persistent equilibrium G * of model system Eqs.
Proof The global asymptotic stability of the methamphetamine-persistent steady state is proved by constructing a global Lyapunov function following [37]. We propose the following Lyapunov function obtained from a suitable combination of common quadratic and Volterra-type functions. Define This function is defined, continuous and positive definite for all S, U , T 0 , T i > 0. We also observe that V(S, U , T 0 , T i ) = 0 at the steady state G * . Since (7) , , , Similarly, using model system Eqs.

Numerical simulations Data
In this section we present an application of our model through fitting the model data on rehabilitation from the Medical Research Council's (MRC's), South African Community Epidemiology Network on Drug Use (SACENDU) project [40]. We fit the model system Eqs.
(2)-(5) to data for individuals under inpatient and outpatient rehabilitation in Cape Town. We use data in Table 1 which was collected by the SACENDU for individuals who were admitted on inpatient and outpatient forms of rehabilitation in Cape Town. While inpatient and outpatient rehabilitation programs may have similar components, they can have very different long-term outcomes and success rates. Statistics on success rate of substance abuse rehabilitation programs vary widely [39] and depend on a variety of factors, such as, the extent and nature of the patient's problems, the appropriateness of treatment and related services used to address those problems, and the quality of interaction between the patient and his or her treatment providers.

Parameter estimation
We make use of Matlab programming language to estimate model parameters used in our numerical simulations and to analyze existing trends on inpatient and outpatient forms of rehabilitation. The model parameters that we use for numerical simulations are in Table 2. We make use of curve fitting, which is a process that allows us to quantitatively estimate the trend of the outcomes. The curve fitting process fits equations of approximating curves to the raw field data. Nevertheless, for a given set of data, the fitting curves of a given type are generally not unique. Thus, a curve with a minimal deviation from all data points is desired. This best-fitting curve can be obtained by the method of least squares. The least squares curve fit routine (lsqcurvefit) in Matlab with optimization is used to estimate the parameter values. Many parameters are known to lie within some intervals. During the estimation of parameter values, a Matlab code is used where unknown parameter values are given a lower and upper bound from which the set of parameter values that provide the best fit are obtained. The intervals used and a few parameters obtained from literature are given in Table 2.

Conclusions
In this paper, we designed a model that incorporates inpatient and outpatient rehabilitation of methamphetamine addicts to study the dynamics of methamphetamine abuse in Cape Town of South Africa. The reproduction number was derived and qualitatively used to investigate the stability of equilibrium states and the prevalence of methamphetamine abuse. The methamphetamine-free equilibrium point is shown to be globally asymptotically stable whenever the reproduction number is less than unity. Thus, methamphetamine abuse can be eliminated if control strategies are put in place and efforts are directed toward reducing the threshold number to a value less than unity. With the aid of a Lyapunov function obtained from a suitable combination of common  . 10 Estimated incidence of methamphetamine abuse using data for outpatient rehabilitants in Cape Town and projected for five more years. The estimated incidence of methamphetamine abuse is expected to gradually decrease for the next five years down to below 11.5 % by the year 2018 quadratic and Volterra-type functions, the methamphetamine-persistent steady state has been shown to be globally asymptotically stable whenever the reproduction number is greater than unity.
The least squares curve fit routine (lsqcurvefit) in Matlab with optimization has been used to fit the model to data on inpatient and outpatient rehabilitants with the objective of using the model parameters that give the best fit to obtain the incidence curve. We also used the parameter values that give the best fit to plot a graph of the ratio of inpatients to outpatients and as well predict future proportions of inpatient and outpatient rehabilitants. The results suggest that the proportion of patients admitted under inpatient rehabilitation facilities in Cape Town will continue to decrease for the next five years whereas that for outpatient rehabilitants will increase for the next five years. The estimated proportion of inpatient rehabilitants in specialist treatment centres of Cape Town was observed to be approximately 31 % by the year 2018. Our estimated incidence of methamphetamine abuse related to data on inpatient and outpatient rehabilitants was observed to be generally decreasing over the years. However, it was noted that the estimated incidence for methamphetamine abuse related to data on inpatient rehabilitants had a sharp decrease as compared to that of outpatient rehabilitants, suggesting that inpatient rehabilitation programs have an increased potential of positively changing the lives of many methamphetamine addicts.   [15] µ Natural death rate 0.019-0.021 0.020 year −1 [41] The projections carried out show that the estimated incidence for methamphetamine abuse related to data on inpatient rehabilitants will have decreased down to below 7 % by the year 2018 and that of outpatient rehabilitants will have decreased to below 11.5 % by the year 2018. The study presented here is not exhaustive, it can be extended to include stratification of the population according to levels of methamphetamine use. Structuring the population in such a way would give some other helpful insights in studying the dynamics of methamphetamine abuse.
In addition, the model did not take into account contextual dynamics, such as drug supply chains or changes in interdiction. Incorporating these processes will undoubtedly facilitate in the understanding of methamphetamine dynamics. Also, since the study focuses specifically on methamphetamine abuse, the model dynamics may also be affected by the availability and use of other drugs.