Abstract
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.
1. 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–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–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 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.
2. Mathematical Model
2.1. 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-illicit-related death rate at time given by Susceptible individuals acquire illicit drug use habits at rateThe parameter measures the strength of interactions between the susceptible individuals and illicit drug users, that is, the “influence” of and on ; is a modification factor which accounts for the increased likelihood of heavy illicit drug users to influence more new drug users compared to light drug users The model takes the formwhere 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 byThe 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)The threshold quantity 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.
Theorem 1. The drug-free equilibrium point is globally asymptotically stable if . If , then a unique drug persistence equilibrium point of system (2) exists and is globally asymptotically stable.
2.2. 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 , 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.
(a)
(b)
(c)
(d)
2.2.1. 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. 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 The results here clearly show that is the most sensitive parameter to . An increase in by 10% would increase 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.
2.3. Optimal Control Problem and Its Analysis
In this section, we introduce two control functions, and , 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) becomesThe 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 over the preselected time interval . Mathematically, this corresponds to the minimization of the functional over a set of feasible, and , strategies on . 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 with the set of admissible control functions given byWe consider the optimal control problem of determining , , , , and , associated with an admissible control pair on the time interval , satisfying (7), given the initial conditions , , , , and , 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.
Theorem 3. Problem (7)–(10) with given initial conditions , , , , and and fixed final time admits a unique optimal solution associated with an optimal control pair on .
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].
2.3.1. Optimal Control Numerical Results
State system (7) and the adjoint system of differential equations together with the control characterization above from the optimality system are solved numerically using the “forward-backward sweep method (FBSM).” For a detailed discussion on FBSM, we refer the reader to [23]. We set , , , and . Parameters used are in Table 2.
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 years is simulated with the following assumed initial population levels: , , , and . 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.
(a)
(b)
(c)
(d)
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.
Simulations in Figure 3 illustrate the feasibility of the control functions and 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 This implies that the cost associated with the implementation of control is greater than the one associated with control For this case, we observe that Figure 3(a) suggests that a slightly higher effort should be devoted to control , since it is more feasible compared to control . In Figure 3(b), it is evident that if the cost associated with control is higher than the cost associated with control , then all the optimal intervention strategies are highly feasible.
(a)
(b)
2.3.2. 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 and denote the corresponding initial condition. Function (11) measures the proportional decrease in the number of active illicit drug users imposed by the intervention with controls , by comparing the number of active illicit drug users at time with the initial conditions for which there are no controls implemented By construction, for all time . Thus, the upper bound of is one.
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.
3. 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.
Appendices
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, over-the-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 byThus,Therefore, spectral radius is
C. Equilibrium States and Their Stability
C.1. Global Stability of the Drug-Free Equilibrium
To explore the global stability of , we consider the Lyapunov function:Its Lyapunov derivative along the solutions to system (2) iswith and . Thus, if with if and only if Therefore, by LaSalle’s Invariance Principle [28], every solution to system (2), with initial conditions in , approaches as That is, as . Substituting in (2) gives as Thus, is globally asymptotically stable in whenever
C.2. Global Stability of the Drug Persistence Equilibrium
Solving system (2) gives
Theorem C.1. If , system (2) is uniformly persistent; namely, there exists a constant such thathere, is independent with initial data in .
Proof. When , from system (2), we have the following limiting system:Observe that we have omitted the last two equations in (2) since the first three equations of system (2) are independent of the variables and Further, for ease of notation, we will continue to denote the drug-free equilibrium point of system (C.5) as DefineIt then suffices to show that (C.5) is uniformly persistent with respect to Firstly, from (C.5), it can be verified that both and are positively invariant. Clearly, is relatively closed in and system (C.5) is point dissipative. Consider the following set using solutions of system (C.5):Now we prove thatAssume that . It suffices to show that for all If not, then there exists such that or . For , we haveIt follows that there is such that for Clearly, we can restrict small enough such that for This means that for , which contradicts the assumption that For other cases, it can also be easily verified that they contradict the assumption that Thus, (C.8) holds.
Note that is globally asymptotically stable in Moreover, is an isolated invariant set in , every orbit in converges to , and is acyclic in By Theorem in [29], we only need to show that if
Now, we prove that Then, there exists a positive solution with , such that as . Since , we choose small enough such thatThus, when is sufficiently large, we have , , , andBy the comparison principle, it follows that , as , which contradicts , as . This proves
Since , , is an isolated invariant set in , and is acyclic in , and by Theorem in [29], hence we are able to conclude that system (C.5) is uniformly persistent with respect to Then, system (2) is uniformly persistent.
We proceed to investigate the global stability of .
Theorem C.2 (see [30]). Let be a function for in an open set Consider the differential equationDenote by the solution to (C.12) such that . A set is said to be absorbing in for (C.12) if for each 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.(3)System (C.12) satisfies a Poincare-Bendixson Property.(4)Each periodic orbit of (C.12) in is asymptotically orbitally stable.
Theorem C.3. The trajectory of any nonconstant periodic solution of system (2), if it exists, is asymptotically orbitally stable with asymptotic phase.
On investigating the global stability of , the equations for the mentally ill individuals and the recovered population can be omitted. Thus, system (2) reduces toWe 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 byThen, the second compound system defined along the periodic solution of system (C.13) is given byBased 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 functionFrom condition (2), we have that there exists a constant such thatfor any and any periodic solution . Let us estimate the right derivative of along a solution of system (C.16):We then need to estimate the following two functions:From system (C.13), we haveThus,Therefore,Denote the period of the periodic solutions by We haveThus, system (C.13) is asymptotically stable. Then, the periodic solutions are asymptotically orbitally stable with asymptotic phase.
Proof of Theorem 3. System (2) is converted into an equivalent problem, namely, the problem of minimizing the Hamiltonian given bywhere , for , 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 and corresponding states (, , , , and ), there exist adjoint functions satisfying with transversality conditions, , for . Observe that the right-hand side of differential equation is , and similarly for the other adjoint functions. The final time boundary conditions (transversality conditions) are zero since there is no dependence on the states at the final time in the objective functional. Furthermore, the optimal controls are characterized byThe control characterization for is obtained on whenever and taking bounds into account, similar to the other controls.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.