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.

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].

Table 1

Recorded prevalence of illicit drug use in Cape Town, South Africa.

Year	96b	97a	97b	98a	98b	99a	99b	00a	00b	01a	01b	02a
%	20	19	22.1	25	35.1	44.1	50.1	54.2	48.1	57.1	54.3	54.3
Year	02b	03a	03b	04a	04b	05a	05b	06a	06b	07a	07b	08a
%	52.8	58.7	62.9	61.1	66.3	65.6	74.9	69.9	73.2	70.5	70.5	69.4

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.

Table 4

Recorded prevalence of illicit drug use for Durban, South Africa.

Year	96b	97a	97b	98a	98b	99a	99b	00a	00b	01a	01b	02a
%	27	31	34	35	33	38	35	40	40	41	42	43
Year	02b	03a	03b	04a	04b	05a	05b	06a	06b	07a	07b	08a
%	42	43.7	44.6	42.1	46	48.9	45.4	48.6	46	50.2	55.2	50.5

Table 5

Recorded prevalence of illicit drug use for Gauteng, South Africa.

Year	96b	97a	97b	98a	98b	99a	99b	00a	00b	01a	01b	02a
%	—	—	—	31	32	33	37	40	40	46	48	46
Year	02b	03a	03b	04a	04b	05a	05b	06a	06b	07a	07b	08a
%	46	47.8	50.7	49.6	49	53.4	48.2	52.5	52.8	54.1	53	53

Table 6

Recorded prevalence of illicit drug use for Mpumalanga, South Africa.

Year	96b	97a	97b	98a	98b	99a	99b	00a	00b	01a	01b	02a
%	—	—	—	—	—	24	24	29	23	30	31	29
Year	02b	03a	03b	04a	04b	05a	05b	06a	06b	07a	07b	08a
%	32	30.9	38.9	36.2	39.2	44.4	45.7	45.5	52.7	56.3	56.7	60.4

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 u(t) 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(t) 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 N(t) is partitioned into subclasses of susceptible individuals S(t) (individuals who are not yet illicit drug users but interact with drug users), light or occasional drug users I(t), heavy drug users I (t), mentally ill population M(t) (individuals who suffer mental illness due to drug use), and detected illicit drug users R(t). Thus, N(t) = S(t) + I(t) + I (t) + M(t) + R(t). We assume a constant size population with recruitment and non-illicit-related death rate at time t 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 I and I on S; κ ≥ 1 is a modification factor which accounts for the increased likelihood of heavy illicit drug users I to influence more new drug users compared to light drug users I. 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 I, M, and I , respectively; σ and ϕ (σ ≤ ϕ) denote the rates at which light and heavy illicit drug users develop mental illness; ψ and d 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 ℰ 0 is globally asymptotically stable if ℛ ≤ 1. If ℛ > 1, 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 ℛ = 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.

Figure 1

Model system (2) fitted to data for individuals seeking treatment due to illicit drug use. The blue circles indicate the actual data and the solid line indicates the model fit to the data.

Table 2

Model parameters and their baseline values obtained after fitting data in Appendix A.

Parameter definition	Symbol	Range	Baseline value	Units
Recovery rate	ω	0.1–0.9	0.3	Per year
Recruitment rate	μ	0.02–0.03	0.02	Per year
Modification factor	κ	1–1.85	1.25	—
Transmission rate	β	0.31–0.36	0.35	Per year
Detection rate for light drug users	γ	0.01–0.14	0.1	Per year
Detection rate for heavy drug users	ρ	0.13–0.78	0.35	Per year
Detection rate for mentally ill population	ϵ	0.54–0.85	0.6	Per year
Escalation of a light user to a heavy user	α	0.01–0.75	0.01	Per year
Drug use-related death for light drug users	ψ	0.01–0.05	0.035	Per year
Drug use-related death for heavy drug users	δ	0.1–0.18	0.14	Per year
Drug use-related death for mentally ill population	d	0.077–0.23	0.2	Per year
Proportion of light users who develop mental illness	σ	0.10–0.33	0.05	Per year
Proportion of heavy users who develop mental illness	ϕ	0.08–0.09	0.09	Per year

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, u, that depends differentiably on a parameter, p, 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.

Table 3

Sensitivity indices of ℛ to parameters for model (2), evaluated at the baseline parameter values given in Table 2.

Parameter	Sensitivity index
β	+1
α	−0.067
κ	0.032
μ	−0.38
ρ	−0.013
d	−0.008
f	−0.0075
ψ	−0.19
σ	−0.059
γ	−0.19
ϕ	−0.0082

2.3. Optimal Control Problem and Its Analysis

In this section, we introduce two control functions, u(t) and v(t), 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, J, 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 u(t) and v(t) over the preselected time interval [0, T]. Mathematically, this corresponds to the minimization of the functional J over a set of feasible, u(t) and v(t), strategies on [0, T]. J functional is defined as follows:where A and B are balancing coefficients transforming the integral into dollars expended over a finite time period of T 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 ℝ 5 with the set of admissible control functions given byWe consider the optimal control problem of determining S (t), I (t), I (t), M (t), and R (t), associated with an admissible control pair (u (t), v (t)) ∈ Γ on the time interval [0, T], satisfying (7), given the initial conditions S(0), I(0), I (0), M(0), and R(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.

Theorem 3 .

Problem (7)–(10) with given initial conditions S(0), I(0), I (0), M(0), and R(0) and fixed final time T admits a unique optimal solution (S (t), I (t), I (t), M (t), R (t)) associated with an optimal control pair (u , v ) on [0, T].

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 S(0) = 0.97, I(0) = 0.02, I (0) = 0.01, and M(0) = R(0) = 0. 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 T = 20 years is simulated with the following assumed initial population levels: S = 0.97, I = 0.02, I = 0.01, and M = R = 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.

Figure 2

Dynamics of system (C.3) showing the effects of optimal control strategies on eliminating or reducing illicit drug use in the community. The control functions u(t) and v(t) were set to 0.95 and A = 5 × 10−7, B = 3 × 10−7.

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 u(t) and v(t). 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 A associated with control function u(t) is greater than parameter B associated with control function v(t). This implies that the cost associated with the implementation of control u(t) is greater than the one associated with control v(t). For this case, we observe that Figure 3(a) suggests that a slightly higher effort should be devoted to control v(t), since it is more feasible compared to control u(t). In Figure 3(b), it is evident that if the cost associated with control v(t) is higher than the cost associated with control u(t), then all the optimal intervention strategies are highly feasible.

Figure 3

Control profiles for control functions u(t) = 0.95 and v(t) = 0.95. In (a) A = 5 × 10−7 and B = 1 × 10−7; in (b) A = 1 × 10−7 and B = 5 × 10−7.

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 E(t)where I (t) and I (t) denote the optimal solutions associated with the optimal control of the corresponding variable and I(0) and I (0) 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 (u, v), by comparing the number of active illicit drug users at time t with the initial conditions for which there are no controls implemented (u = v = 0). By construction, E(t)∈[0,1] for all time t. Thus, the upper bound of E(t) 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.

Figure 4

Time series plot demonstrating the efficacy of optimal intervention strategies over a period of 500 days.

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.