Computational and Mathematical Methods in Medicine
Volume 2015 (2015), Article ID 383154, 11 pages
http://dx.doi.org/10.1155/2015/383154

Research Article

Modeling Illicit Drug Use Dynamics and Its Optimal Control Analysis

Steady Mushayabasa and Gift Tapedzesa

University of Zimbabwe, Department of Mathematics, P.O. Box MP 167, Harare, Zimbabwe

Received 20 April 2015; Accepted 13 October 2015

Academic Editor: Seiya Imoto

Copyright © 2015 Steady Mushayabasa and Gift Tapedzesa. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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.

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 rate

The 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 form

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)

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.

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

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

(a)

(b)

(c)

(d)

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.

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.

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

Parameter
Sensitivity index

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

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) becomes

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 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 by

We 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)

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 and were set to 0.95 and , .

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)

Figure 3: Control profiles for control functions and . In (a) and ; in (b) and .

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.

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.

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 by

Thus,

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) is

with 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 that

here, 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 Define

It 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 that

Assume that . It suffices to show that for all If not, then there exists such that or . For , we have

It 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 that

Thus, when is sufficiently large, we have , , , and

By 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 equation

Denote 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 to

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

Then, the second compound system defined along the periodic solution of system (C.13) is given by

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

From condition (2), we have that there exists a constant such that

for 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 have

Thus,

Therefore,

Denote the period of the periodic solutions by We have

Thus, 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 by

where , 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 by

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

References

United Nations Office on Drugs and Crime (UNODC), World Drug Report, 2014, 1: Analysis.UNODC, United Nations Office on Drugs and Crime (UNODC), 2014.

L. Degenhardt, W. Hall, M. Lynskey, and M. Warner-Smith, “Illicit drug use,” in Comparative Quantification of Health Risks: Global and Regional Burden of Disease Attributed to Selected Major Risk Factors, M. Ezzati, A. D. Lopez, A. Rodgers, and R. Murray, Eds., pp. 1109–1176, World Health Organization, Geneva, Switzerland, 2004.

Global Burden of Disease (GBD), Mental Disorders and Illicit Drug Use Expert Group, Global Burden of Disease (GBD), Geneva, Switzerland, 2005, http://www.gbd.unsw.edu.au/.

M. Frischer, S. T. Green, and D. Goldberg, Substance Abuse Related Mortality: A World Wide Review, United Nations International Drug Control Programme, Vienna, Austria, 1994.

M. Donoghoe, “Illicit drugs,” in Quantifying Global Health Risks: The Burden of Disease Attributable to Selected Risk Factors, C. J. L. Murray and A. D. Lopez, Eds., Harvard University Press, Cambridge, Mass, USA, 1996.

R. M. Murray, P. D. Morrison, C. Henquet, and M. Di Forti, “Cannabis, the mind and society: the hash realities,” Nature Reviews Neuroscience, vol. 8, no. 11, pp. 885–895, 2007.

Canadian Mental Health Association, Mental Illness and Substance Use Disorders: Key Issues, Canadian Mental Health Association, 2005.

A. Pluddemann, S. Dada, C. Parry et al., “Monitoring alcohol and drug abuse trends in South Africa, South Africa Community Epidemiology Network on drug use (SACENDU),” SACENDU Research Brief, vol. 11, no. 2, 2008.

J. L. Gerberding, “Occupational HIV infection,” AIDS, vol. 11, supplement, pp. S57–S60, 1997.

E. H. Kaplan and R. Heimer, “A model-based estimate of HIV infectivity via needle sharing,” Journal of Acquired Immune Deficiency Syndromes, vol. 5, no. 11, pp. 1116–1118, 1992.

F. Nyabadza, J. B. H. Njagarah, and R. J. Smith, “Modelling the dynamics of crystal meth (‘Tik’) abuse in the presence of drug-supply chains in South Africa,” Bulletin of Mathematical Biology, vol. 75, no. 1, pp. 24–48, 2013.

G. P. Samanta, “Dynamic behaviour for a nonautonomous heroin epidemic model with time delay,” Journal of Applied Mathematics and Computing, vol. 35, no. 1-2, pp. 161–178, 2011.

E. White and C. Comiskey, “Heroin epidemics, treatment and ODE modelling,” Mathematical Biosciences, vol. 208, no. 1, pp. 312–324, 2007.

H. A. Whiteford, L. Degenhardt, J. Rehm et al., “Global burden of disease attributable to mental and substance use disorders: findings from the Global Burden of Disease Study 2010,” The Lancet, vol. 382, no. 9904, pp. 1575–1586, 2013.

M. Marmor and P. Hartsock, “Self-destructing (non-reusable) syringes,” The Lancet, vol. 338, no. 8764, pp. 438–439, 1991.

S. Mushayabasa and C. P. Bhunu, “Mathematical analysis of hepatitis C model for intravenous drug misusers: impact of antiviral therapy, abstinence and relapse,” Simulation, vol. 90, no. 5, pp. 487–500, 2014.

S. Mushayabasa and C. P. Bhunu, “Hepatitis C virus and intravenous drug misuse: a modeling approach,” International Journal of Biomathematics, vol. 7, no. 1, Article ID 1450006, 2014.

S. Mushayabasa, “On the role of optimal counseling and antiviral therapy on controlling HCV among intravenous drug misusers,” Journal of Computational Medicine, vol. 2014, Article ID 347827, 7 pages, 2014.

S. Mushayabasa, “The role of optimal intervention strategies on controlling excessive alcohol drinking and its adverse health effects,” Journal of Applied Mathematics, vol. 2015, Article ID 238784, 11 pages, 2015.

S. Mushayabasa and C. P. Bhunu, “Epidemiological consequences of non-compliance to HCV therapy among intravenous drug users,” International Journal of Research and Reviews in Applied Sciences, vol. 8, no. 3, pp. 288–295, 2011.

United Nations Office on Drugs and Crime (UNODC), World Drug Report, 2012 1: Analysis, UNODC, 2012.

S. Mushayabasa, C. P. Bhunu, and R. J. Smith, “Assessing the impact of educational campaigns on controlling HCV among women in prison settings,” Communications in Nonlinear Science & Numerical Simulation, vol. 17, no. 4, pp. 1714–1724, 2012.

S. Lenhart and J. T. Workman, Optimal Control Applied to Biological Models, Chapman & Hall, London, UK, 2007.

R. Leander, S. Lenhart, and V. Protopopescu, “Using optimal control theory to identify network structures that foster synchrony,” Physica D: Nonlinear Phenomena, vol. 241, no. 5, pp. 574–582, 2012.

L. Arriola, “Sensitivity analysis for quantifying uncertainty in mathematical models,” in Mathematical and Statistical Estimation Approaches in Epidemiology, G. Chowell, J. M. Hyman, L. M. A. Bettencourt, and C. Castillo-Chavez, Eds., pp. 195–248, Springer, New York, NY, USA, 2009.

L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko, The Mathematical Theory of Optimal Processes, Wiley, Hoboken, NJ, USA, 1962.

P. Van den Driessche and J. Watmough, “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission,” Mathematical Biosciences, vol. 180, pp. 29–48, 2002.

J. LaSalle, The Stability of Dynamical Systems, vol. 25 of CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, Pa, USA, 1976.

H. R. Thieme, “Persistence under relaxed point-dissipativity with an application to an epidemic model,” SIAM Journal on Mathematical Analysis, vol. 24, no. 2, pp. 407–435, 1993.

M. Y. Li and L. Wang, “Global stability in some SEIR epidemic models,” in Mathematical Approaches for Emerging and Reemerging Infectious Diseases: Models, Methods, and Theory, C. Castillo-Chavez, S. Blower, P. Van den Driessche, D. Denise Kirschner, and A. Yakubu, Eds., vol. 126 of The IMA Volumes in Mathematics and Its Applications, pp. 295–311, Springer, New York, NY, USA, 2002.

M. Fiedler, “Additive compound matrices and inequality for eigenvalues of stochastic matrices,” Czechoslovak Mathematical Journal, vol. 99, pp. 392–402, 1974.

M. Y. Li and J. S. Muldowney, “A geometric approach to global-stability problems,” SIAM Journal on Mathematical Analysis, vol. 27, no. 4, pp. 1070–1083, 1996.