The Role of Family in Initiating Methamphetamine Abuse Treatment: Insights through a Mathematical Model

It has been known for many years that the treatment gap" is massivethat is, among those who need treatment for a substance use disorder, few receive it [1]. For instance, the National Institute on Drug Abuse (NIDA) reported that after providing a survey on drug use that was conducted in 2007, 23:2 million people in the United States at least 12 years of age required treatment for drug use. However, only 2:4 million of this population actually received treatment at a facility specializing in drug rehabilitation. This means that about 8:4% of the US population ages 12 years and up needed treatment for drug use but did not receive it [2]. The question arises, \how do we reduce this gap?". Strategies include increasing access to effective treatment, reducing stigma and raising awareness among both drug addicts and their contacts of the value of addiction treatment. Families and friends play an incredibly critical role in motivating individuals with drug problems to enter treatment [3]. There is no single, immutable definition of family. For more information about a family [4]. A private, non-confrontational and honest talk between a non-addicted family member and an addicted family member can implore the addicted family member to seek treatment. Family members may disperse around the world, but still be connected emotionally and able to contribute to the dynamics of family functionally. Nearly all addicted individuals believe at the onset of their drug using career that they can stop using drugs on their own, and most try to stop without treatment. Although some people are successful, many attempts result in failure to achieve long-term abstinence [3]. The goal of treatment is to return people to productive functioning in the family, workplace and community. Methamphetamine (MA) is among the main targeted drugs of abuse by the Western Cape Department of Social Development (DOSD) due to its treacherous characteristics. Non drug users’ especially young girls are usually lured by its weight loss properties. It is a cheap street drug which can be made in clandestine kitchen laboratories. Despite methamphetamine's illegality, recipes for making the drug are found easily on the Internet and passed from user to user. It is this aspect that has given law enforcement agents torrid times due to its hard trafficking nature. MA use has been linked to sexually transmitted infections and risky sexual behavior [5]. Morris and Parry [6], reported high rates of increase of HIV infections in communities with high MA use. Thus, MA use has immense public health implications. That is why mathematical models have been recently designed to gain insight into the spread of MA within communities [7,8]. One of the first ordinary differential equation (ODE) models of opiate addiction [9] designed using the principles of mathematical epidemiology became instrumental in the development of the first mathematical model for MA [7]. The authors [7] introduced a mathematical model given by non-linear ODE's, extending the model [9] to qualitatively study the dynamics of MA use in the Western Cape Province (WCP) of South Africa. The authors [7] divided the class of drug users not on treatment into two classes: firstly, a class of light drug users (this class considers drug users on the period of concealed drug use at the beginning) and secondly, a class of addicted drug users which was referred to as hard drug users. Their extension also caters for the class of recovered drug users. Their model predicts that public health measures such as increased behavior change can reduce the prevalence of MA users in the community and an application based on real data is given. In Cape Town of South Africa, the MA epidemic has to some extent been exacerbated by drug lords who target innocent school children in disguise by selling them sweets tainted by drugs. This thereby has led to an extension of the work in [7] to include recruitment to drug use by drug lords [8]. More recently,


Introduction
It has been known for many years that the treatment gap" is massivethat is, among those who need treatment for a substance use disorder, few receive it [1]. For instance, the National Institute on Drug Abuse (NIDA) reported that after providing a survey on drug use that was conducted in 2007, 23:2 million people in the United States at least 12 years of age required treatment for drug use. However, only 2:4 million of this population actually received treatment at a facility specializing in drug rehabilitation. This means that about 8:4% of the US population ages 12 years and up needed treatment for drug use but did not receive it [2]. The question arises, \how do we reduce this gap?". Strategies include increasing access to effective treatment, reducing stigma and raising awareness among both drug addicts and their contacts of the value of addiction treatment. Families and friends play an incredibly critical role in motivating individuals with drug problems to enter treatment [3]. There is no single, immutable definition of family. For more information about a family [4]. A private, non-confrontational and honest talk between a non-addicted family member and an addicted family member can implore the addicted family member to seek treatment. Family members may disperse around the world, but still be connected emotionally and able to contribute to the dynamics of family functionally. Nearly all addicted individuals believe at the onset of their drug using career that they can stop using drugs on their own, and most try to stop without treatment. Although some people are successful, many attempts result in failure to achieve long-term abstinence [3]. The goal of treatment is to return people to productive functioning in the family, workplace and community. Methamphetamine (MA) is among the main targeted drugs of abuse by the Western Cape Department of Social Development (DOSD) due to its treacherous characteristics. Non drug users' especially young girls are usually lured by its weight loss properties. It is a cheap street drug which can be made in clandestine kitchen laboratories. Despite methamphetamine's illegality, recipes for making the drug are found easily on the Internet and passed from user to user. It is this aspect that has given law enforcement agents torrid times due to its hard trafficking nature. MA use has been linked to sexually transmitted infections and risky sexual behavior [5]. Morris and Parry [6], reported high rates of increase of HIV infections in communities with high MA use. Thus, MA use has immense public health implications. That is why mathematical models have been recently designed to gain insight into the spread of MA within communities [7,8]. One of the first ordinary differential equation (ODE) models of opiate addiction [9] designed using the principles of mathematical epidemiology became instrumental in the development of the first mathematical model for MA [7]. The authors [7] introduced a mathematical model given by non-linear ODE's, extending the model [9] to qualitatively study the dynamics of MA use in the Western Cape Province (WCP) of South Africa. The authors [7] divided the class of drug users not on treatment into two classes: firstly, a class of light drug users (this class considers drug users on the period of concealed drug use at the beginning) and secondly, a class of addicted drug users which was referred to as hard drug users. Their extension also caters for the class of recovered drug users. Their model predicts that public health measures such as increased behavior change can reduce the prevalence of MA users in the community and an application based on real data is given. In Cape Town of South Africa, the MA epidemic has to some extent been exacerbated by drug lords who target innocent school children in disguise by selling them sweets tainted by drugs. This thereby has led to an extension of the work in [7] to include recruitment to drug use by drug lords [8]. More recently, another mathematical model was developed in [10] to model the impact of rehabilitation, amelioration and relapse on the prevalence of drug epidemics. Their results show that both prevention and treatment/ rehabilitation are necessary strategies for reduction of drug epidemics and some recommendations were made with regards to both treatment and preventive strategies, that they should be directed toward reducing the contact rate and treatment should be combined with psychotherapy to accelerate quitting and reduce re-initiation.
In spite of the above extensions, the modeling efforts directed towards providing useful insights on MA use have not captured the notion that patients under treatment can be `in-patients' or `outpatients', except possibly for some recent work done in [11]. However, some authors do acknowledge that this is worth an issue to look into using mathematical models, see for instance [7]. Methamphetamine remains the substance of choice among 20 patients in the Western Cape province of South Africa and the province remain heavily affected by the methamphetamine epidemic [12]. The specialist substance abuse treatment centers in the Western Cape have felt an increasing pressure to treat young methamphetamine using clients.
The other underlying motivation for the development of this paper was that looking at the age of most patients who seem to be under the custodian of parents, the reaction put forward by South African Community Epidemiology Network on Drug Use (SACENDU) data makes it reasonable why there is an increase in referrals by family/ self/friends as compared to other referrals. In this paper, we add new features that were not captured in previous MA mathematical models. Precisely, we introduce: (i) A term that describes the roles of families and friends in initiating MA abuse treatment; (ii) Two treatment options for drug users namely, in-patient rehabilitation and out-patient rehabilitation; (iii) An incidence term of newly admitted rehabilitants with probabilities β 2 and β 3 for joining inpatient and outpatient rehabilitation respectively.
After some positive talk between a non-addicted family member and an addicted family member, an addicted family member will be compelled to enter either an inpatient or an outpatient rehabilitation program. Each patient's needs and means are different, and outpatient and inpatient programs have varying benefits for patients and family. 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 [13]. Involvement in an outpatient rehabilitation program means that patients are not separated from their families and can continue substance abuse treatment for an extended amount of time. Inpatient (residential) rehabilitation allows patients to undergo an intensive detoxification and recovery program. They are immersed in the recovery process and do not have the ability to leave the substance abuse treatment campus [14].
We formulate and establish the basic properties of the model. We re-scale parameters to allow for a comprehensive analysis of the model. Numerical simulations, parameter estimation and sensitivity analysis are also presented.

The Model and its Parameters
The local population is divided into four classes, S(t), U(t), Rop(t) and Rip(t). The class S(t) represents the population at risk of being initiated into drug abuse where all individuals range from age 15-64 years. The class U(t) denotes those initiated into drug abuse together with relapsed drug users, Rop(t) denotes those in rehabilitation as outpatients and Rip(t) denotes those in rehabilitation as in-patients. The total local population is thus given by N(t)= S(t)+ U(t)+ Rop(t)+Rip(t).
The remaining parameters given in model system equations (1)(2)(3)(4) are as follows: Λ: The rate at which the general population enter the susceptible population, that is, the demographic process of individuals reaching age 15 in the modelling time period.      ρ2: The proportion of drug users undergoing inpatient rehabilitation relapsing into drug use. 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. γ 1 : The transfer rate of individuals from outpatient rehabilitation to inpatient rehabilitation. 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. γ 2 : The transfer rate of individuals from inpatient rehabilitation to outpatient rehabilitation. 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 model involves certain assumptions. These consist of the following: (i) Individuals in each compartment are indistinguishable and there is homogeneous mixing.
(ii) The initiation of people into drug use together with the initiation of drug users into rehabilitation are functions that are analogous to the force of infection for epidemic models.
(iii) Drug users in treatment are not infectious to susceptible.
(iv) The referrers come from non-drug users.
The schematic diagram below shows the movement of humans as their status with respect to drug use changes. The model presented in Figure 1 above may be represented by the set of nonlinear ordinary differential equations below together with the initial conditions:

Proof:
We obtain that the total population satisfies the differential equation Applying a theorem by Birkho and Rota [15], on differential inequalities, we have where N(0) represents the value of equations (1-4) evaluated at the initial values of the respective variables. Taking the limit as t → ∞, we have that 0 Thus, the state variables remain biologically meaningful in the set S U R R remain in  for all t > 0. Therefore the ω-limit set of solutions of equations (1)-(4) in  are contained in . The uniqueness, existence and continuity results hold for equations (1-4). The system is thus mathematically and epidemiologically well-posed [16]. Our analysis will be based on the dynamics of the solutions generated in .

Drug-free equilibrium and the abuse reproduction number
Presented here are detailed analysis of the abuse reproduction number  a for the model and stability analysis of the drug-free equilibrium. The model has a drug-free equilibrium given by a scenario depicting a drug-free state in the community or society. The abuse reproduction number  a . of the model is defined herein in the drug-using context as the average number of people that each drug user will initiate to drug use during the drug-using career in a population of completely potential drug users. Usually,  a < 1 implies that drug abuse will die out, whereas  a > 1 implies that drug abuse will persist within a community and  a = 1 requires further investigation. The determination of a is done using the next generation matrix approach [17]. This method has been explored in many papers [18][19][20][21][22][23]. Using this method we have

Local Stability of the drug-free steady state
We shall now prove the local stability of the drug-free equilibrium point 0 whenever the abuse reproduction number is less than unity.

Theorem 2:
The drug-free equilibrium point 0 of model system

Proof:
The Jacobian matrix of model system equations (1)-(4) at  0 is given by where h, h1 and h2 are de_ned as before. The local stability of the drugfree equilibrium is determined by the following submatrix of J( 0 ), Since all o_-diagonal elements are positive, we now consider matrix (1 ) where X2 is a positive 3 × 1 matrix given by Then, it follows from M-matrix theory that all eigenvalues of J 1 ( 0 ) have negative real part, which implies the local asymptotic stability of the drug-free equilibrium if  a < 1. On the other hand, it can be shown that the determinant of J 1 ( 0 ) is given by Thus, if  a < 1, then matrix J 1 ( 0 ) has eigenvalues with negative real parts, which implies the stability of the drug-free equilibrium. This completes the proof.

Global stability of the drug-free steady state
We shall now prove the global stability of the drug-free equilibrium point  0 whenever the reproduction number is less than unity.
where 1 2 Then, following the definition and properties of M-matrices, the following conditions hold for the matrix V: 1. V is an M-matrix.
2. V-1 exists and is a non-negative matrix.
Multiplying a row vector (a 1 , a 2 , a 3 ) to the above equation (6) where a 1 , a 2 and a 3 are constants, we obtain Choose Thus, since V-1 is non-negative from the properties listed above, we have a 1 ≥ 0, a 2 ≥ 0 and a 3 ≥ 0. In particular, Substituting (8) to the right hand side of (7), we obtain For the choice of ai, i = 1, 2, 3 in (8), define a Lyapunov function Computing the time derivative of ν along solutions of model system equations (1)-(4), we obtain We also have the following equations and * 2 3 1 upon substituting expression (12) into the third and fourth equation of model system equations (1-4) respectively. Substituting equations (14) and (15) into (13) leads to Solving (17) gives U* = 0 which corresponds to the drug-free equilibrium or gives Substituting (19) into expressions (14), (15) and (16), we obtain the drug persistent equilibrium point

Numerical Simulations Data
We carry out detailed numerical simulations using matlab programming language to support our theoretical findings. We also present an application of our study through fitting the model to data from Tables 1 and 2, which was collected by the Medical Research Council's (MRC's), South African Community Epidemiology Network on Drug Use (SACENDU) [24].

Parameter estimation
In order to illustrate the model and check the feasibility of our analysis regarding existence conditions and stability of various equilibria, we conducted some numerical simulations by choosing the values of parameters in model system equations (1-4) as given in Table 3. We make use of Matlab programming language to estimate model parameters used in our numerical simulations and to identify parameters of interest in the methamphetamine using career. The collection of data relevant to such parameters by epidemiologists and treatment providers would be of significant use from an implementation and policy perspective. 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 feet routine 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 _t are obtained. The intervals used and a few parameters obtained from literature are given in Table 3. Data from SACENDU clearly reflects the referral sources for patients on admission. Nearly half of the patients admitted for the past fifteen years have been referred to treatment centers by their families or friends ( Table 2).
Data from SACENDU also reflects the type of treatment received by each patient in treatment centers of Cape Town (Table 1). From the year 1999 up to the first half of the year 2008, the majority of patients were treated on an inpatient basis. Starting from the second half of the year 2008 up to the last review period, the majority of patients were now being treated on outpatient basis. There are many factors that influence the decisions an individual has to make when attending rehab. For instance, cost of the treatment program, length of the treatment program and the nature and extent of aftercare services provided, will influence someone's decision on whether to choose inpatient or outpatient rehab. This makes the task of ascertaining the exact parameter values related to the uptake into inpatient or outpatient rehabilitation challenging. Thus, due to limitation of data in this aspect, we shall resort to parameter estimation.  family or friends. However, the proportions for this source of referrers have remained fairly high as compared to other referral sources [24]. Figure 2 demonstrates a good fit to the data from Table 2. Figure 2 shows that the proportion of patients in treatment centers of Cape Town whose referrals come from family/friends, is likely to remain stable for the next five years after the last review period.
We make use of dimensionless variables so that quantities on the vertical axis are now proportions. We use the re-scaling formulae:- , to obtain the following re-scaled subsystem of (1) The numerical simulations are carried out using the following initial conditions: s(0) = 0:95, u(0) = 0:05, rop(0) = 0, rip(0) = 0. When  a is less than unity, the generation of methamphetamine users is lower than their predecessors. Thus, methamphetamine use is eventually eliminated from the population (Figure 3). When Ra is greater than unity, the subsequent generation of methamphetamine users is greater than their predecessors. Therefore, methamphetamine abuse will persist in the population (Figure 4).
Latin hypercube sampling: Sensitivity analysis assesses the amount and type of change inherent in the model as captured by the terms that define the reproductive number. If the reproductive number is very sensitive to a particular parameter, then a perturbation of the conditions that connect the dynamics to such a parameter may prove useful in identifying policies or intervention strategies that reduce epidemic prevalence. We use Latin Hypercube Sampling (LHS) technique to perform sensitivity analysis of the abuse reproductive number  a . Latin hypercube sampling is a statistical sampling method that allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [25]. Partial rank correlation coefficients (PRCCs) were calculated to estimate the correlation between values of the reproductive number and the associated model parameters across 1000 random draws from the empirical distribution of  a and its associated parameters. Here, we only demonstrate for a few important selected parameters constituting the abuse reproduction number. Figure 5 shows partial rank correlation coefficients (PRCCs), which were calculated to estimate the correlation between a and constituent    parameters. Parameters with positive PRCCs will increase  a when they are increased, whereas parameters with negative PRCCs will decrease  a when they are increased. Figure 5 suggests that increasing the proportion of referrers leads to a decrease in the value of  a . Also, an increase in β 1 (probability of becoming a drug user) leads to an increase in the value of  a . interestingly, the transition from inpatient rehabilitation to outpatient rehabilitation is resulting in an increase in the value of  a . This might be a result of the reduced chances of a methamphetamine user relapsing during inpatient rehabilitation because they will be completely immersed in the program and separated from the lifestyle and habits that supported drug use. So when they transfer to outpatient rehabilitation where they have contact with people and circumstances that trigger addiction, the chances of relapsing may be more likely as compared when they were inpatients. We also examined the dependence of  a on these parameters in more detail. We used Latin Hypercube Sampling and Monte Carlo Simulations to run 1000 simulations, where all parameters were simultaneously drawn across their ranges [26][27][28][29][30]. Figure 6 illustrates the effect of the transmissibility parameter on the abuse reproductive number  a . It can be seen that as β 1 increases to a level above 65%, the drug epidemic converges to the endemic equilibrium. Therefore efforts targeted at reducing this transmissibility parameter will be useful for the reduction of  a . Figure 6b shows that if the proportion of referrers is greater than 75% then the abuse reproduction number  a can be kept to a level below one. This is an indication that referrers play a critical role in the fight against methamphetamine abuse.

Effect of the referrers proportion on the number of drug addicts:
We use parameter values that give the best fit to plot graphs for the number of drug addicts attributed to each indicated proportion of referrers. Each proportion has a different impact on the spread of the drug epidemic. These fractions are useful in bringing up information that will help to guide intervention strategies to slow the epidemic. Figure 7 shows the effect of varying the proportion of referrers on the number of drug addicts. The figure demonstrates that as we increase the proportion of referrers, the number of drug addicts decreases as shown by the arrow pointing downwards. This indicates that family and friends' involvement in enabling an addict to attend rehab is very crucial. Thus, efforts targeted at educating family members to understand how supportive they can be to their addicted loved ones in attaining a life of sobriety would be useful in controlling the methamphetamine epidemic [31][32][33][34][35][36][37][38][39].

Conclusion
In this paper, we designed a model that captures the best-case scenario where an addicted individual will be compelled to choose either an inpatient or an outpatient substance abuse treatment program after the intervention of family/friends. A key novelty of our model is an incidence term of newly admitted rehabilitants with probabilities β 2 and β 3 for joining inpatient and outpatient rehabilitation respectively. The model captures the reality of positive family/friends' involvement in motivating a drug addict to enter treatment. The abuse reproduction number was derived and qualitatively used to investigate the stability of equilibrium states. The drug-free equilibrium point is shown to be globally asymptotically stable whenever the abuse reproduction number is less than unity. With the aid of a Lyapunov function obtained from a suitable combination of Volterra-type functions, the drugpersistent steady state has been shown to be globally asymptotically stable whenever the abuse reproduction number is greater than unity. Latin hypercube sampling and partial rank correlation coefficients demonstrate that the parameter v, proportion of referrers, has great impact on the outcome. Increasing the proportion of referrers leads to a decrease in the value of  a . The results indicate that family members of drug addicts have a critical role to play in the dynamics of methamphetamine abuse. This is also seen from the results which indicate that as we increase the proportion of referrers, the number of drug addicts decreases. This suggests that family and friends' involvement in enabling an addict to attend rehab is very crucial. Thus, specific health and/or education programs may be employed to educate family members on how to encourage their addicted loved ones to initiate methamphetamine treatment. This will facilitate in controlling the methamphetamine epidemic.