Modeling Illicit Drug Use Dynamics and Its Optimal Control Analysis

The global burden of death and disability attributable to illicit drug use, remains a significant threat to public health for both developed and developing nations. This paper presents a new mathematical modeling framework to investigate the effects of illicit drug use in the community. In our model the transmission process is captured as a social “contact” process between the susceptible individuals and illicit drug users. We conduct both epidemic and endemic analysis, with a focus on the threshold dynamics characterized by the basic reproduction number. Using our model, we present illustrative numerical results with a case study in Cape Town, Gauteng, Mpumalanga and Durban communities of South Africa. In addition, the basic model is extended to incorporate time dependent intervention strategies.


Introduction
Illicit drug use continues to exert significant toll, with valuable human lives and productive years of many people being lost. An estimated 183,000 drug-related deaths were reported in 2012 [1]. Globally, it is estimated that, in 2012, between 162 million and 324 million people aged between 15 and 64 had used an illicit drug [1]. Illicit drug use is defined as the nonmedical use of a variety of drugs that are prohibited by international law [2]. These drugs include amphetamine-type stimulants, cannabis, cocaine, heroin, and other opioids, and MDMA (ecstasy) [3]. The risk of premature mortality and morbidity from illicit drug use is dependent on dose, frequency, and route of administration [2]. Further, the mortality risks of illicit drug consumption increase with increasing frequency and quantity of consumption [4]. The first global burden of death and disability attributable to illicit drug use was first estimated by Donoghoe [5], as part of the global burden of disease project [6]. Donoghoe estimated that illicit drug use was responsible for 10,000 deaths globally in 1990 with about 62% of them in developing countries [6]. Besides causing deaths, illicit drug use has been associated with 50% of mental illness cases [7]. Table 1 illustrates the cumulative prevalence of illicit drug use in Cape Town, South Africa, from 1996 (June) to 2008 (June); the data was adopted from SACENDU [8].
Additional data on recorded cases of illicit drug use in some communities of South Africa, namely, Durban (Table 4), Gauteng (Table 5), and Mpumalanga (Table 6), is presented in Appendix A.
Despite many clinical and theoretical studies [9][10][11][12][13][14][15][16][17][18][19][20] and tremendous educational campaigns [1,7,21,22], illicit drug use remains a significant threat to public health all over the world. The use of mathematical models to explore the dynamics of drug use has been an interesting topic for a couple of researchers [11][12][13]. In 2007, White and Comiskey [13] proposed a mathematical model to evaluate the role of treatment and relapse in the dynamics of heroin. Their work revealed among others that relapse of individuals who would have quit heroin use has a significant impact on the generation of new or secondary heroin users. Most recently, Nyabadza et al. [11] constructed a mathematical model to examine the dynamics of crystal meth "Tik" abuse in the presence of drug-supply chains. Using the data from South Africa, their work suggests among others that programs aimed at encouraging light drug users to quit drug use can be more effective to control "Tik" abuse in South Africa.  Our objective in this paper is to formulate a mathematical model for illicit drug use that includes relevant biological and social aspects, accounts for case detection of drug users, and allows optimal control methods to be used. We stress that our model is not for a specific substance abuse. To begin with, we integrate the aforementioned essential components into one SIR-type (Susceptible-Infectious (individuals who are drug users)-Recovered (individuals in rehabilitation)) model to accommodate the diverse dynamics of illicit drug use determined by population specific parameters such as the rates of light drug user to heavy drug user, light drug user to mental illness, and heavy drug user to mental illness. Our study differs from those in the literature in the fact that it accounts for individuals who become mentally ill due to illicit drug use. It is worth noting that the stigma attached to substance abuse and mental disorders often hinders early detection, diagnosis, and proper treatment [7]. We are aware that mental illness can lead to drug abuse. However, in this paper, we consider that individuals become mentally ill due to illicit drug use.
After comprehensive mathematical analysis, we extend the basic model to incorporate two control functions: the goal of the first control function ( ) is to reduce the intensity of "social influence," that is, to weaken the intensity of social interaction between the susceptible population and illicit drug users, while the second V( ) attempts to increase the rate of detection and rehabilitation of illicit drug users. Optimal control theory has found wide-ranging applications in solving biological problems [23] and network structures [24]. Numerical results followed by a brief discussion round up the paper.

Model Framework.
A population of size ( ) is partitioned into subclasses of susceptible individuals ( ) (individuals who are not yet illicit drug users but interact with drug users), light or occasional drug users ( ), heavy drug users ( ), mentally ill population ( ) (individuals who suffer mental illness due to drug use), and detected illicit drug users ( ). Thus, ( ) = ( ) + ( ) + ( ) + ( ) + ( ). We assume a constant size population with recruitment and non-illicitrelated death rate at time given by . Susceptible individuals acquire illicit drug use habits at rate The parameter measures the strength of interactions between the susceptible individuals and illicit drug users, that is, the "influence" of and on ; ≥ 1 is a modification factor which accounts for the increased likelihood of heavy where the upper dot represents the derivative of the component with respect to time. The parameter is the rate at which light drug users become heavy drug users; , , and denote the rates of detection and rehabilitation of individuals in classes , , and , respectively; and ( ≤ ) denote the rates at which light and heavy illicit drug users develop mental illness; and model the permanent exit rates of light and heavy users, respectively, due to either cessation or drug use-related death. Individuals in rehabilitation recover at rate and are assumed to permanently exit the model. Mentally ill individuals permanently exit the model at rate due to drug use-related death. Further, we assume that the mentally ill population does not influence the susceptible individuals to become illicit drug users. We study system (2) in the closed set: where Ω is positively invariant with respect to (2). In the absence of drug abuse in the community system (2) admits a drug-free equilibrium point given by The basic reproductive number provides an invasion criterion for the initial spread of the drug misuse in a susceptible population. For this case, the reproductive number is given by (see computations in Appendix B) R = ( + + + + ) ( + + + + ) ( + + + ) .
The threshold quantity R gives the average number of new drug abusers generated by a typical light or heavy drug user during his or her lifetime as an illicit drug user. Theorem 1 follows from computations in Appendix C.

Numerical Results.
In this section, we estimate the model parameters used in our numerical simulations. In order to carry out the estimation, we made use of the data from SACENDU (see Appendix A). The data has also been used to estimate the best fit for our model. Figure 1 shows "best" fit to the real data, illustrated in Appendix A. Here, the fitting process involved the use of the least squares-curve fitting method. A Matlab code has been used together with unknown parameters assigned lower and upper bound from which a set of parameter values that produce the best fit have been obtained. With the exception of Gauteng, simulation results in Figure 1 suggest that the prevalence of illicit drug use increases rapidly over the defined time intervals, while for Gauteng, the increase is gradual. Using the lower bounds of parameter values in Table 2, direct calculation shows that R = 6.22, which shows that effective methods should be designed to reduce illicit drug use in these communities. Further, from these simulations, one can suggest that system (C.3) can be an essential tool in predicting future illicit drug use in the communities.

Sensitivity Analysis of the Reproductive Number.
Sensitivity analysis of model parameters is very important to design and control strategies as well as a direction to future research. Local sensitivity indices allow us to measure the relative change in a state variable when a parameter changes. Table 3: Sensitivity indices of R to parameters for model (2), evaluated at the baseline parameter values given in Table 2.

Parameter
Sensitivity index In computing the sensitivity analysis, we adopt the approach described by Arriola [25]. The normalized forward sensitivity index of a variable to a parameter is the ratio of the relative change in the variable to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives.
Definition 2. The normalized forward sensitivity index of a variable, , that depends differentiably on a parameter, , is defined as To determine the numerical output for our sensitivity indices, we used the lower bounds of parameters in Table 2. Table 3 presents the model parameters and their sensitivity indices. Model parameters whose sensitivity index values are near −1 or +1 suggest that a change in their magnitude has a significant impact on either increasing or decreasing the size of the reproductive number R . The results here clearly show that is the most sensitive parameter to R . An increase in by 10% would increase R by 10%. Thus, effective methods aimed at reducing illicit drug and its adverse effects in the community should be designed with the strong view of weakening strength of interactions between the susceptible individuals and illicit drug users.

Optimal Control Problem and Its Analysis.
In this section, we introduce two control functions, ( ) and V( ), to system (2). The goal of the first control is to reduce the intensity of "social influence" while the second models the effort on the detection of illicit drug users. Thus, system (2) becomeṡ The objective functional, , is used to formulate the relevant optimization problem: finding the most effective strategy that reduces or eliminates the levels of illicit drug use in the community at minimal cost. This minimization goal will be achieved through the implementation of controls ( ) and V( ) over the preselected time interval [0, ]. Mathematically, this corresponds to the minimization of the functional over a set of feasible, ( ) and V( ), strategies on [0, ]. functional is defined as follows: where and are balancing coefficients transforming the integral into dollars expended over a finite time period of years. The balancing coefficients account for the relative size and importance preassigned by the modelers to the contributing terms in the objective functional. We consider state system (7) of ordinary differential equations in R 5 with the set of admissible control functions given by We consider the optimal control problem of determining * ( ), * ( ), * ( ), * ( ), and * ( ), associated with an admissible control pair ( * ( ), V * ( )) ∈ Γ on the time interval [0, ], satisfying (7), given the initial conditions (0), (0), (0), (0), and (0), and minimizing the cost functional (8); that is, The existence of optimal controls follows from standard results in optimal control theory [26]. Theorem 3 follows from Appendix C. The optimal control pair predicted by (5) represents the optimal intervention strategy, given cost constraints, and can be found by the application of the Pontryagin maximum principle [26].
Numerical results in Figure 2 illustrate the role of optimal intervention strategies in controlling illicit drug use in the community. In Figure 2, the balancing coefficients are fixed within the integral expression (8) and the optimal schedule of the two controls over = 20 years is simulated with the following assumed initial population levels: = 0.97, = 0.02, = 0.01, and = = 0. System (7) is used to determine the resultant population dynamics. From the illustrations, it is clear that optimal intervention strategies provide a more effective approach on the elimination or reduction of illicit drug use and mental health cases. In Figure 2(a), we note that in the presence of two optimal intervention strategies it can require 14 years for the population of light illicit drug users to die off. Although interventions which are not time dependent lead to a decrease in cases of light illicit drug use, implemented for 14 years, they may not yield cessation in illicit drug use cases.
In Figure 2(b), we set the initial population level of heavy drug users at 0.01 (1%) of the total population. Although this might be a high figure in reality, it allows us to observe the level of effectiveness of optimal intervention strategies. It is evident in Figure 2 that the presence of the aforementioned optimal intervention strategies can lead to elimination of heavy drug users even if their population is considerably high.
In Figures 2(c) and 2(d), the initial densities of mentally ill and detected population were set to zero. Results here show that in the presence of time dependent interventions mental illness associated with illicit drug use will remain low and will eventual die off after 14 years of implementing the aforementioned optimal intervention strategies. Further, it is evident that more illicit drug users will be detected when time dependent interventions are implemented. Figure 3 illustrate the feasibility of the control functions ( ) and V( ). If the control profile starts at the upper bound and drops on the final time to the origin, then it is regarded as highly feasible, implying that more effort and resources should be devoted to such an intervention strategy for effective problem-solving or control. Figure 3(a) assumes that parameter associated with control function ( ) is greater than parameter associated with control function V( ). This implies that the cost associated with the implementation of control ( ) is greater than the one associated with control V( ). For this case, we observe that Figure 3(a) suggests that a slightly higher effort should be devoted to control V( ), since it is more feasible compared to control ( ). In Figure 3(b), it is evident that if the cost associated with control V( ) is higher than the cost associated with control ( ), then all the optimal intervention strategies are highly feasible.

Efficacy of Optimal Intervention Strategy.
The efficacy of an intervention strategy on controlling the generation of new illicit drug users reflects the strength of the strategy to effectively control the drug use problem in the community. In this section, we explore the effectiveness of the aforementioned optimal intervention methods on reducing cumulative population illicit drug users. We define the efficacy function ( ) where * ( ) and * ( ) denote the optimal solutions associated with the optimal control of the corresponding variable and  Figure 4 demonstrates that optimal interventions may be effective to eliminate illicit drug use in the community after 14 years of implementation. It is worth noting that, after 4 years of implementing these methods, the efficacy level would reach 80% mark; this demonstrates that even for a short time horizon optimal intervention strategies can have a significant impact on controlling illicit drug use in the community.

Concluding Remarks
About 5% (230 million) of the people of the world's adult population are estimated to have used an illicit drug at least once in 2010 [21]. In this paper, a simple mathematical model to assess the impact of illicit drug use and its adverse health effects is proposed and analyzed. The model is further extended to incorporate two optimal intervention strategies. The goal of the first optimal intervention is tied on reducing the "social influence" between the susceptible and illicit drug users, while the second aims to increase the detection and rehabilitation of illicit drug users. From the illustrations in this study it is clear that time dependent intervention strategies can lead to elimination of illicit drug use in the community. Our analysis also suggests that the implementation of these controls should be devoted to both controls since they are highly feasible for a longer period considered in this study. Although the model is simple and has a couple of simplifying assumptions, qualitative conclusions can be reached in rather broad terms from the simulations and analysis, the kind of conclusions that can generate or provide useful insights into value and effectiveness of time dependent control efforts, aimed at elimination or effective control of illicit drug use and its adverse health problems.
Our model is not exhaustive; a new model that incorporates illicit drug use and age or gender can be developed and used to forecast future trends on illicit drug use.

A. Recorded Data on Illicit Drug Use
In this section, we present cumulative prevalence of illicit drug use in South Africa, for the following areas: Durban, Gauteng, and Mpumalanga. The data was adopted from Pluddemann et al. (2008) [8]. We obtained cumulative prevalence after summing the prevalences for the following illicit drugs: cannabis, Mandrax, cocaine, heroin, ecstasy, cocaine, overthe-counter and prescription medicines (OTC/PRE), methamphetamine, and other substances/polysubstance abuse. Key: 96b represents data recorded in 1996 from July to December and 97a represents January to June 1997.

B. Computation of the Reproductive Number
Following Van den Driessche and Watmough [27] and using the notation defined therein, the matrices and for the new infection terms and the remaining transfer terms are, respectively, given by

C.2. Global Stability of the Drug Persistence Equilibrium.
Solving system (2) gives here, is independent with initial data in Ω.
We proceed to investigate the global stability of E * . Theorem C.2 (see [30]). Let → ( ) ∈ R be a 1 function for in an open set ⊂ R . Consider the differential equatioṅ = ( ) . (C.12) Denote by ( , 0 ) the solution to (C.12) such that (0, 0 ) = 0 . A set is said to be absorbing in for (C.12) if ( , 1 ) ⊂ for each 1 ⊂ and is sufficiently large. If the following conditions are satisfied, then the unique equilibrium is globally asymptotically stable: (1) There exists a compact absorbing set ⊂ and (C. 12) has a unique equilibrium in in .
(2) is locally asymptotically stable. (C. 13) Computational and Mathematical Methods in Medicine 9 We study system (C.13) in the closed set: where Φ is positively invariant with respect to (C.13). The corresponding associated second compound matrix [31,32] is given by [2] = ( ) . Based on Theorem C.2, if we can prove that system (C.13) is asymptotically stable, then the periodic solution is asymptotically orbitally stable with asymptotic phase. We define the Lyapunov function We then need to estimate the following two functions: where , for = 1, 2, . . . , 5, are the adjoint functions associated with their respective states. Note that, in (C.25), each adjoint function multiplies the right-hand side of the differential equation of its corresponding state function. The first terms in come from the integrand of the objective functional. Given an optimal control pair ( * , V * ) and corresponding states ( * , * , * , * , and * ), there exist adjoint functions satisfyinġ The control characterization for * is obtained on / = 0 whenever 0 < * ( ) < 1 and taking bounds into account, similar to the other controls.