# Analysis of a Heroin Epidemic Model with Saturated Treatment Function

Authors: Isaac Mwangi Wangari; Lewi Stone
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

A mathematical model is developed that examines how heroin addiction spreads in society. The model is formulated to take into account the treatment of heroin users by incorporating a realistic functional form that “saturates” representing the limited availability of treatment. Bifurcation analysis reveals that the model has an intrinsic backward bifurcation whenever the saturation parameter is larger than a fixed threshold. We are particularly interested in studying the model’s global stability. In the absence of backward bifurcations, Lyapunov functions can often be found and used to prove global stability. However, in the presence of backward bifurcations, such Lyapunov functions may not exist or may be difficult to construct. We make use of the geometric approach to global stability to derive a condition that ensures that the system is globally asymptotically stable. Numerical simulations are also presented to give a more complete representation of the model dynamics. Sensitivity analysis performed by Latin hypercube sampling (LHS) suggests that the effective contact rate in the population, the relapse rate of heroin users undergoing treatment, and the extent of saturation of heroin users are mechanisms fuelling heroin epidemic proliferation.


---

## Body

## 1. Introduction

In 1897, Germany’s Bayer pharmaceutical company synthesised heroin and soon after marketed the product as a nonaddictive miracle drug, for use as a cough syrup and pain reliever [1]. Cough medicine was in fact in high demand, since tuberculosis and pneumonia were fast-spreading diseases of the time. As such, the miracle drug heroin was rapidly disseminated across the globe. Fast forward to today, and we know that addiction to heroin is an extremely common phenomenon among heroin users; some 23% of individuals who consume the drug become dependent on it. Worldwide, many countries are affected by the heroin drug-trafficking industry and its growing number of users. America is currently in the midst of another heroin epidemic [2] with approximately 700,000 Americans using heroin in the past year [2]. The number of people using heroin for the first time is increasing at an alarming rate, with >150,000 Americans engaging in heroin use in 2012, which is almost double that recorded in 2006 [2]. Heroin also leads to other diseases and is considered a major pathway responsible for fuelling proliferation of human immunodeficiency virus (HIV) and Hepatitis B and Hepatitis C virus (HBV, HCV) [3, 4].The development of heroin habituation and addiction has similar characteristics to an epidemic, in terms of its disturbingly contagious spread through a susceptible population. In the last decades, a whole range of mathematical models have been developed to forecast how diseases spread in time and space and how they can be controlled. Recently, the same mathematical modelling techniques have been extended for the purpose of understanding and combating drug addiction problems. The aim of the present study is to propose a novel heroin epidemic model and make use of it to study issues arising with treatment and establish conditions that may signal heroin persistence within the community.The ultimate goal of mathematical epidemiology is to understand how to control and eliminate infectious diseases and these ideas have a place for also dealing with social problems. In epidemic theory the basic reproduction number, usually denoted byR0, is one of the most important concepts, given its ability to predict the course of an epidemic. It will also prove invaluable in our study of heroin dynamics in society. R0 is defined as the number of secondary infections that are likely to occur when a single infectious individual is introduced into an entirely susceptible population [5]. Until recently, it has been widely accepted that the condition R0<1 is an essential requirement for the eradication of a disease. However, this viewpoint has been recently challenged with a number of theoretical studies demonstrating that this criterion may not always be sufficient. Instead, the phenomenon of backward bifurcation offers a different interpretation since it shows that although the basic reproduction number is below unity and the infection free equilibrium is stable, there might still be another stable endemic equilibrium and unstable endemic equilibrium coexisting simultaneously. Thus even though R0<1, a population may still reside at an endemic equilibrium in which the disease persists indefinitely. In a scenario where multiple equilibria concurrently exist, the extinction or persistence of an epidemic is dependent on initially infected size of subpopulations. The qualitative features of backward bifurcation are illustrated in Figure 1.Figure 1
Illustration of the qualitative features of backward bifurcation. The red line represents the unstable equilibria (i.e., unstable endemic equilibria and unstable heroin-free equilibrium) while the blue line represents stable equilibria (i.e., stable endemic equilibria and stable heroin-free equilibrium).A variety of both behavioural and pharmacological medications can be administered to effectively treat heroin addiction. The side effects associated with quitting using heroin (such as pain, diarrhoea, nausea, and vomiting) are very severe and very often compel heroin addicts to relapse. To prevent such cases there are available medications that can be administered during the detoxification stage to relieve craving and physical symptoms. A number of studies have established that pharmacological therapy has positive impact in facilitating drug addicts to remain in treatment programs. Furthermore, it has been noted that during addiction treatment there is a decline in drug consumption, infectious disease transmission, and crime rates [2]. In our present study we shall incorporate a saturated treatment function and derive threshold conditions that indicate when heroin is able to persist within a community. Besides incorporation of a saturated treatment function, our model also included an extra class of individuals, namely, those who have been successfully treated from heroin using. This class has been neglected in previous heroin epidemic models [6–9]. Much of our work will be focused on exploring the conditions for global stability of the heroin model with treatment. Our work deals with global stability of a heroin model with bilinear incidence rate, self-cure, relapse, and saturated treatment function using the Bendixson criterion.With this in mind we will extend the SIR (Susceptible-Infected-Recovered) model by [8] to represent a heroin epidemic model and investigate global stability properties. To be precise we study the conditions of global stability for the nontrivial equilibrium states by using two distinct approaches: the Lyapunov direct method and the Li and Muldowney’s geometric approach to global stability. It is with no doubt that the famous Lyapunov direct method is a powerful tool for nonlinear stability analysis [10]. One of the main advantages of Lyapunov direct technique is that it is directly applicable to nonlinear systems [11]. However, the major challenge with using Lyapunov direct method is that it requires an auxiliary function which is often hard to construct. And this difficulty is exacerbated especially if the model exhibits backward bifurcation phenomena because Lyapunov functions for such models may not exist. To address these difficulties another powerful tool, the geometric technique due to Li and Muldowney, was developed in the middle of nineties [11–13]. Their method involves generalization of Bendixson’s criterion to systems of any finite dimensions and applies compound matrices. Presently this method has gained popularity due to its vast range of applications, in particular to mathematical models that are of biological interest. Although this method is mainly applied in epidemic models (e.g., see, [14–19]) its use can be found in other population dynamics contexts (see [20]). It has been shown in [16] that the geometric technique is more appropriate in mathematical models of SEIR-like structure since their analysis can be easily reduced to a three-dimensional system. Nevertheless, the method has been extended to four-dimensional systems that may be difficult to reduce. In the sequel applications to four-dimensional systems are rare because the procedure becomes mathematically involved when n≥4. Examples of four-dimensional systems can be traced in the work of Ballyk and coworkers who applied compound matrices to a four-dimensional population model [21] and also by Gumel and coworkers [22] who studied a SVEIR (Susceptible-Vaccinated-Exposed-Infected-Recovered) model of severe acute respiratory syndrome (SARS) epidemic spread.The four-dimensional model studied here can be reduced to a three-dimensional system. Both Lyapunov direct method and geometric approach are applied to investigate global properties of a four-dimensional heroin epidemic model. Lyapunov direct method will be applied in a special case in particular where the parameter that triggers bistability phenomena is switched off. On the other hand geometric approach will be applied in the general model where all parameters are present including the one that causes bistability. Here we follow the procedure in [11, 16] to obtain sufficient condition for global stability.


## 2. Model Formulation

In the spirit of the SIR (Susceptible-Infected-Recovered) model in the literature (i.e., [23]), we formulate a heroin epidemic model based on the assumption that heroin use follows a process that can be modelled similar to infectious diseases [24, 25]. The general population is stratified into four mutually exclusive classes, namely, susceptibles (S), individuals successfully treated from heroin use (U3), heroin users undergoing treatment (U2), and heroin users not in treatment (U1). The proposed heroin epidemic model is based on key assumptions which include the following:(i)
Uniform mixing: individuals in the above-mentioned classes freely interact with each other.(ii)
Individuals undergoing treatment are still often using drugs [26].(iii)
Heroin users in treatment relapse to heroin users not in treatment as a result of the self-decision to terminate treatment [27].(iv)
Heroin users in treatment do not infect susceptibles.Given these assumptions the heroin model may be described by the processes illustrated in Figure2, which can be written in terms of the following set of equations:(1)dSdt=Λ-βU1S-μS,dU1dt=βU1S+pU2-μ+δ1+ξU1-TU1,dU2dt=TU1-p+σ+δ2+μU2,dU3dt=σU2+ξU1-μU3.In brief, the susceptible subpopulation S(t) is generated at a constant rate through immigration and birth at rate Λ. Some susceptible individuals who come into contact with heroin users U1(t) may begin to use heroin. Hence, the susceptible population is diminished due to contact with heroin users at rate βU1S, while heroin users increase at the same rate. Heroin users also increase when those undergoing treatment relapse at rate pU2 and return to their heroin using lifestyle. Heroin users reduce in number as a result of treatment which is represented by the treatment function T(U1). Moreover, the user subpopulation is reduced by heroin-induced death at rate δ1U1 as well as a result of the self-decision to cease using heroin (also referred to as “self-cure”) at rate ξU1. Individuals undergoing treatment are diminished through the following processes: relapse to heroin using at rate pU2, heroin-induced death at rate δ2U2, and successful treatment at rate σU2. Finally, the recovered/successfully treated subpopulation U3(t) is generated when heroin users undergoing treatment are successfully cured and also through “self-cure.” All subpopulations are decreased by natural death via the background mortality parameter μ.Figure 2
A heroin epidemic model with bilinear incidence rate and saturated treatment function. The blue solid arrows represent deaths due to either heroin or natural causes while the black solid arrows represent change of status from one compartment to another.Heroin epidemic models studied to date [6–9] assume the classical view that the treatment rate of the infective population should be proportional to the number of infective individuals [28]. This viewpoint was criticised during the SARS (Severe Acute Respiratory Syndrome) outbreaks in 2003. The dramatic increase of SARS cases in Beijing challenged the normal public-health system because it was only possible to treat a limited number of SARS patients at a given time. The experience with SARS epidemic sparked a renewed interest among modellers to investigate the implication of the capacity of the healthcare system. Authors in [29] considered an SIR epidemic model and assumed a Heaviside treatment function while Wang [30] restudied the same SIR model but assumed a piecewise linear treatment function. Here we will assume that the heroin users U1(t) receive treatment based on the following more general saturated treatment function: (2)TU1≜αU11+ωU1, where α is positive and ω is nonnegative. In our present model the parameter ω accounts for the extent of saturation of heroin users. Note that for small U1 the treatment function reduces to T(U1)≈αU1 while for large U1 it reduces to T(U1)≈α/ω which actually characterizes the saturated phenomena of the treatment. Further, if ω=0, the treatment function becomes T(U1)=αU1 which is the usual linear treatment rate. 1/(1+ωU1) is a measure of inhibition due to a saturation of heroin users who are usually too many to be dealt with given the limited available treatment.A summary of the model variables and parameters is given in Table1.Table 1
Description of variables and parameters of model (1).
Variable
Description
S
Number of susceptible individuals at timet
U
3
Number of heroin users who have been successfully treated from heroin use, as well as individuals who have voluntarily stopped using the heroin (and have withdrawal symptoms) “self-cure" at timet
U
2
Number of heroin users undergoing treatment at timet
U
1
Number of drug users not undergoing treatment at timet, that is, the initial and relapsed heroin users
N
Total population at timet (N=S(t)+U3(t)+U2(t)+U1(t))
Parameter
Description
Λ
Recruitment rate of individuals in the general population entering the susceptible population
β
Effective contact rate
μ
Death rate due to natural causes
p
Rate at which individuals undergoing treatment relapse to heroin use
ξ
“Self-cure rate": rate at which heroin users stop using heroin and join the successfully cured class of individuals not taking heroin
δ
1
Heroin-related death rate of heroin users not in treatment
δ
2
Heroin-related death rate of individuals undergoing treatment
σ
Rate at which heroin users in treatment are successfully cured (i.e., completely detoxicated) from heroin use
α
Rate at which heroin users are treated
ω
The extent of saturation of heroin users within the community


## 3. Basic Properties and Basic Reproduction Number

Since we are studying a human population, the model must be able to ensure that all the associated parameters and the state variablesS, U3, U2, U1 are nonnegative for all time t>0. Hence, the following result.Theorem 1.
Let the initial conditions supplied to model (1) be such that S(0)>0,U3(0)>0,U2(0)>0, and U1(0)>0. Then the trajectories (S(t),U3(t),U2(t),U1(t)) of the model (1), with positive initial conditions, will remain positive for all time t.Proof.
Lett1=sup{t>0:S(t)>0,U3(t)>0,U2(t)>0,U1(t)>0}>0. Now from the first equation of model (1), it follows that (3)dSdt=Λ-ΦS-μS≥Λ-ΦS-μSwhereΦ=βU1,which can be written as (4)ddtStexpμt+∫0tΦτdτ≥Λexpμt+∫0tΦτdτ.Hence, (5)St1expμt1+∫0t1Φτdτ-S0≥∫0t1Λexpμx+∫0xΦτdτdx,so that (6)St1≥S0exp-μt1-∫0t1Φτdτ+exp-μt1-∫0t1Φτdτ×∫0t1Λexpμx+∫0xΦτdτdx>0. Following a similar procedure we can show that U3(t)>0, U2(t)>0, and U1(t)>0 for all time t>0. Thus, all trajectories of model (1) remain positive for all nonnegative initial conditions, as required.Now in what follows we establish the region where model (1) is considered to be biologically feasible. Summing all the equations of the basic model (1) yields(7)dNtdt=Λ-μNt-δ1U1t-δ2U2t.Considering that 0<U1(t)<N(t), 0<U2(t)<N(t) and letting δ¯=max{δ1,δ2}, it follows from (7) that (8)Λ-μ+2δ¯Nt≤dNtdt<Λ-μNt. Therefore (9)Λμ+2δ¯≤liminf⁡t→∞⁡Nt≤limsup⁡t→∞⁡Nt≤Λμsuchthatlimsup⁡t→∞⁡Nt≤Λμ.Theorem 2.
The closed set(10)℧=S,U3,U2,U1:0≤S,U3,U2,U1;S+U3+U2+U1≤Λμis positively invariant and absorbing with respect to the set of nonlinear differential equation (1).Proof.
Here we show that the feasible solutions of model (1) are uniformly bounded in the region ℧. Suppose S, U3, U2, and U1 are any solution of system (1) supplied with nonnegative initial conditions. Then it is straightforward to note that the total population N satisfies the inequality(11)dNdt=Λ-μN-δ1U1-δ2U2≤Λ-μN.From (11) it follows that dN/dt≤Λ-μN which implies dN/dt≤0 if N≥Λ/μ. The standard comparison theorem [31] can be used to deduce that N(t)≤N(0)e-μt+Λ/μ(1-e-μt). In particular N(t)≤Λ/μ if N(0)≤Λ/μ for all t>0. Thus, under the flow induced by system (1), the region ℧ is positively invariant. Furthermore, for N(0)>Λ/μ the trajectory solutions N(t) either enter in the region ℧ in finite time or asymptotically approach Λ/μ. Thus, in the region ℧, model (1) is said to be mathematically and epidemiologically well posed [32] and the solution of all the trajectories generated by model (1) is considered in a biologically feasible region ℧.Clearly system (1) has an intrinsic Heroin-free equilibrium (HFE) given by D0=(S0,0,0,0), a scenario representing a heroin-free state in the community. S0=Λ/μ represents the number of susceptibles when no one is using heroin. The basic reproduction number denoted by R0 is defined as the number of secondary infections that are likely to be triggered by a single infectious individual when introduced into a wholly susceptible population [32]. Here R0 is interpreted as the mean number of secondary cases of heroin users generated by a typical heroin user not in treatment during his/her duration of heroin use in a population of potential drug users.To obtain the basic reproduction number we observe that the average time an individual spends as a heroin user without treatment isT0=1/(μ+δ1+ξ+α) and the probability of surviving this compartment and moving to the treatment compartment is T1=α/(μ+δ1+ξ+α). Now the probability of surviving heroin users in treatment class and then returning to the heroin users class not in treatment is T2=p/(p+σ+δ2+μ). Thus, the total average time spent by the heroin users not in the treatment compartment on multiple passes can be obtained as(12)T=T01+T1·T2+T1·T22+⋯.Clearly, the terms inside the square brackets in (12) constitute a geometric sequence (see Appendix A for detailed derivation) and therefore expression (12) can be written as(13)T=p+σ+δ2+μμ+δ1+ξp+δ2+σ+μ+ασ+δ2+μ.Multiplying (13) with the effective contact rate β and the average recruitment rate Λ/μ we obtain heroin basic reproduction number as(14)R0=βΛp+δ2+μ+σμαμ+δ2+σ+μμ+δ1+ξp+δ2+μ+σ.It is easy to observe that R0 is inversely proportional to treatment α, which implies that if treatment rate is maintained sufficiently high it can control a heroin epidemic (by reducing R0 to less than one). However, as we will see later when parameter ω (representing the extent of saturation of heroin users) is accounted for, this control is no longer guaranteed.Theorem 3.
The HFE is locally asymptotically stable providedR0<1; otherwise it is unstable.This general result has been reviewed in [33] and thus not proved again here. The theorem implies that heroin users will disappear from the community when R0<1 if the initial sizes of the subpopulations of system (1) are in the basin of attraction of the heroin-free equilibrium.Remark 4.
It is instructive to note that the basic reproduction number does not include the parameterω that accounts for the extent of saturation of heroin users. In the next subsection we investigate endemic equilibria of the model where the parameter ω plays a key role in emergence of bistability.


### 3.1. Endemic Equilibria

Within the context of our heroin model the endemic equilibrium refers to a state when heroin addiction is maintained over long time scales in the population. ThereforeS∗,U3∗,U2∗,U1∗>0 hold. To obtain the endemic equilibria E∗=(S∗,U3∗,U2∗,U1∗) we set (1) to zero and solve for the equilibrium quantities S∗,U3∗, and U2∗ in terms of U1∗. That is,(15)S∗=ΛβU1∗+μ,U3∗=σαU1∗+ξp+σ+δ2+μ1+ωU1∗U1∗μp+σ+δ2+μ1+ωU1∗,U2∗=αU1∗p+σ+δ2+μ1+ωU1∗.Substituting (15) into the second equation of (1) yields(16)fU1∗=AU1∗2+BU1∗+C=0,where (17)A=μ+δ1+ξp+σ+δ2+μβω,B=μ+δ1+ξp+σ+δ2+μβ+μωμ+δ1+ξp+σ+δ2+μ+ασ+δ2+μβ-βΛωp+σ+δ2+μ,C=μ+δ1+ξp+σ+δ2+μμ+αμσ+δ2+μ1-R0.The quadratic equation (16) can be analysed to investigate the existence of multiple equilibria when the basic reproduction number is below unity.If the parameter that accounts for the extent of saturation of heroin users in model (1) is excluded, that is, ω=0, (16) reduces to a linear equation (18)U1∗B~+C=0,where (19)B~=μ+δ1+ξp+σ+δ2+μβ+ασ+δ2+μβso that model (1) has the unique solution (20)U1∗=-CB~which is nonnegative if and only if R0>1. Hence, if ω=0, model (1) has a unique endemic equilibrium whenever R0>1 and this equilibrium approaches zero as R0 tends to one (R0→1+) because C→0. But there are no positive endemic equilibria if R0<1. These results are summarized in the following lemma.Lemma 5.
The epidemic model (1) when ω=0 has a unique positive endemic equilibrium E1∗=(S∗,U3∗,U2∗,U1∗) whenever R0>1 and no positive endemic equilibrium otherwise.In what follows, we investigate the global stability for both the HFE and the unique endemic equilibriumE1∗ for the case ω=0.


### 3.2. Global Stability for Heroin-Free Equilibrium Whenω=0

To investigate global stability we apply the method presented by Castillo-Chavez et al. [34]. First let X=(S,U3) and Y=(U2,U1) with X∈R2 representing the number of individuals not using heroin and Y∈R2 representing the number of individuals using heroin (i.e., heroin users in treatment and heroin users not in treatment). Now suppose (21)X′=FX,Y,Y′=GX,Y,GX,0=0,where X′ and Y′ denote differentiation with respect to time. The HFE is now denoted by H0=(X0,0), where X0=(S0,0). The following conditions (H1) and (H2) have to be met to guarantee a local asymptotic stability: H1
ForX′=FX,0,  X0is globally asymptotically stable (g.a.s).(H2)
G(X,Y)=BY-G^(X,Y), where G^(X,Y)≥0, for (X,Y)∈℧.B
=
D
Y
G
(
X
,
0
) and ℧ is the region where model (1) is biologically realistic. Then, Castillo-Chavez and Song [35] have shown that the following lemma is satisfied.Lemma 6.
The fixed pointH0=(X0,0,0) is g.a.s equilibrium of model (1) provided that R0<1 (locally asymptotically stable) and that assumptions (H1) and (H2) hold.Now consider the following theorem.Theorem 7.
SupposeR0<1. Then the HFE D0 is g.a.s.Proof.
LetX=(S,U3) and Y=(U2,U1), and H0=(X0,0,0), where X0=Λ/μ,0.
Then we have(22)X′=FX,Y=Λ-βU1S-μSσU2+ξU1-μU3.It is straightforward to see that, at the heroin-free equilibrium (HFE) S=S0=Λ/μ, F(X,0)=00. Thus, (23)X′=FX,0=Λ-μS-μU3.
Now, ast→∞, X→X0. Hence,X0 is globally asymptotically stable (i.e., condition (H1) is satisfied).
Now consider(24)GX,Y=-p+σ+δ2+μαpβS0-μ+δ1+ξ+αU2U1-0βU1S0-S so that (25)B=-p+σ+δ2+μαpβS0-μ+δ1+ξ,G^X,Y=0βU1S0-S. Since the total population is bounded by N=S+U3+U2+U1≤Λ/μ, we have S≤N≤Λ/μ. Thus, G^(X,Y)≥0, which now implies that conditions (H1) and (H2) are satisfied. Consequently by Lemma 6 the fixed point H0 is globally asymptotically stable when R0<1, which indicates nonexistence of multiple nontrivial equilibria when ω=0. The epidemiological implication of HFE being g.a.s is that heroin epidemic will be eliminated from the community if the threshold quantity R0 is decreased to (or maintained at) a value below unity.Now forω>0 we establish the following theorem.Theorem 8.
Forω>0 model (1) has (i)
a unique positive endemic equilibrium ifB<0 and either C=0 or B2-4AC=0,(ii)
a unique positive endemic equilibrium ifC<0 (i.e., R0>1) and B<0,(iii)
two positive endemic equilibria ifC>0,B<0, and B2-4AC>0,(iv)
no positive endemic equilibrium ifB>0 and either C>0 or B2<4AC.The theorem may be proved as follows. It is obvious to note that in quadratic equation (16) A is always positive and C is either positive or negative depending on whether the basic reproduction number is less than or greater than one, respectively.For Case (i) whereB<0 and C=0 (i.e., R0=1) (16) becomes linear AU1∗+B=0 and has a unique nonzero solution U1∗=-B/A which is positive if B<0 and negative if B>0. Referring to (15) we see that if U1∗ is unique then so are S∗,U2∗, and U3∗.For Case (ii) whereC<0 (i.e., R0>1) and B<0, (16) is quadratic and according to Descartes Rule of Signs (see [36]), (16) has one change of signs indicating (16) has a unique positive root and therefore there is a unique endemic equilibrium.In Case (iii) whereB<0 there is a nonnegative endemic equilibrium at R0=1. However, because (16) is quadratic and since the equilibrium is continuously determined by R0 then there must be an interval to the left of R0=1 on which two nonnegative equilibria coexist. That is,(26)U1,1∗=-B-B2-4AC2A,U1,2∗=-B+B2-4AC2A.For Case (iv) where B>0 and C>0 or B2<4AC, (16) has no positive real root as can be seen in (26), implying nonexistence of a positive endemic equilibrium.Case (iii) suggests that model (1) exhibits the phenomenon of backward bifurcation since the classical requirement for the occurrence of the phenomenon of backward bifurcation is satisfied, that is, the existence of multiple equilibria when the basic reproduction number is less than one. Thus, we have the following lemma.Theorem 9.
Model (1) has backward bifurcation at R0=1 if and only if B<0 (i.e., ω>ωc).Proof.
Consider (16), f(U1∗)=AU1∗2+BU1∗+C=0. Note that, at R0=1, C=0 implies that the graph f(U1∗) passes through the origin. If B<0 it follows that f(U1∗)=0 has a nonnegative root. Since f(U1∗) is a continuous function of C, if we increase C such that C>0, there is some open interval of C say (0,ψ) on which f(U1∗)=0 has two nonnegative roots. That is, there exist two nonnegative endemic equilibria when R0<1. This is indeed true since Case (iv) of Theorem 8 has already shown that for B≥0 model (1) does not have positive real roots when R0<1. Note that, at R0=1, C=0 the following equality holds:(27)μ+δ1+ξp+σ+δ2+μμ+αμμ+δ2+σ=p+σ+δ2+μβΛ.This together with condition B<0 implies that(28)ω>μ+δ1+ξp+σ+δ2+μβ+ασ+δ2+μβμασ+δ2+μ≜ωc.Thus, the phenomenon of backward bifurcation (referring to Case (iii), a situation where there are two endemic equilibria) occurs at the left ofR0=1 if and only if condition (28) is satisfied. This suggests that backward bifurcation will only occur if the parameter ω that accounts for the extent of saturation of heroin users exceeds a certain threshold (i.e., ω>ωc). However, if ω<ωc backward bifurcation cannot occur. Thus, parameter ω plays a critical role in the formation of backward bifurcation for model (1). It is instructive to note that similar results as the one shown in inequality (28) can be obtained by center manifold theory (see Appendix B), where it is emphasized that if ω>ωc the bifurcation coefficient a is positive indicating that the model system (1) undergoes the phenomenon of backward bifurcation. The epidemiological implication of backward bifurcation is that although it is necessary to reduce the basic reproduction number below one it is not sufficient to eradicate a heroin epidemic; rather R0 should be reduced further below a certain threshold which we shall denote by R0C (see Figure 3(b)).Figure 3
(a) and (b) represent bifurcations where drug users,U1∗, not in treatment equilibrium are plotted as a function of R0. The blue solid line represents the stable equilibrium while the red solid line represents the unstable equilibrium. (a) represents forward bifurcation and the parameters used include ω=0.11<ωc=0.1156, μ=0.01, β=0.001, δ1=0.002, δ2=0.001, ξ=0.015, α=0.9, p=0.467, σ=0.1, and Λ∈{1,3} while (b) represents backward bifurcation and parameters used are the same as in (a) except ω=0.25>ωc=0.1156. (c) depicts the critical value R0C as a function of saturation parameter ω. The black dotted line represents the threshold ωc which if exceeded gives rise to backward bifurcation. (d) depicts the critical value R0C as a function of treatment rate α.
(a)
(b)
(c)
(d)


### 3.3. Computation of New Threshold for Heroin EradicationR0C

Here we compute the value of the basic reproduction number where the two nontrivial endemic equilibria (both stable and unstable) collide and annihilate each other leaving only the heroin-free equilibrium point as the stationary solution. This isR0C in Figure 3(b). We choose Λ as the parameter of backward bifurcation. As we have noted in Case (iii) of Theorem 8, (16) has nonnegative roots corresponding to two endemic equilibria if and only if C>0 (i.e., R0<1) and B<0, B2>4AC. It follows that if B=-2AC, (16) has one nonnegative root -B/2A. Supposing there is backward bifurcation at R0=1, then there are two endemic equilibria for an interval of values of basic reproduction number starting from a threshold R0C defined by B=-2AC to a point where R0=1. To obtain this threshold R0C which is often referred to as the critical value of the basic reproduction number we replace the values of A, B, and C into the equality B2=4AC to obtain a quadratic equation in terms of Λ. For mathematical tractability we redefine coefficients B and C as B=B1-B2Λ and C=C1-C2Λ, where(29)B1=μ+δ1+ξp+σ+δ2+μβ+μωμ+δ1+ξp+σ+δ2+μ+ασ+δ2+μβB2=βωp+σ+δ2+μ,C1=μμ+δ1+ξp+σ+δ2+μ+αμσ+δ2+μ,C2=βp+σ+δ2+μ,B1,B2,C1,C2>0.Note that A>0 remains as previously defined in (16). Now the quadratic equation in terms of Λ can be obtained as (30)B22Λ2+4AC2-2B1B2Λ+B12-4AC1=0 Since we are considering the scenario where B<0 and R0<1, we have C1B2>C2B1 and thus just the single solution (31)Λc=B1B2-2AC2+2A2C22+AB2C1B2-C2B1B22>0.Thus, the critical value of basic reproduction number (i.e., the new threshold for heroin eradication), R0C, is given as (32)R0C=βΛcp+δ2+μ+σμαμ+δ2+σ+μμ+δ1+ξp+δ2+μ+σ.Consequently, from the above analysis of computation of threshold for heroin eradication we can deduce the following lemma.Lemma 10.
(a) IfR0>1, then model (1) has a unique endemic equilibrium point E∗. In this case heroin epidemic will persist in the community.
(b) IfR0C<R0<1, then model (1) has two endemic equilibria E¯1 and E¯2 and signals that model (1) has backward bifurcation.
(c) IfR0<R0C<1, then model (1) has only the heroin-free equilibrium point D0 and in this case heroin users will disappear.Figure3 exhibits typical bifurcation diagrams for model (1). To obtain the graphs we vary recruitment rate Λ while other parameter values are held fixed. The parameters used for the numerical simulation that leads to Figure 3(a) include ω=0.11, μ=0.01, β=0.001, δ1=0.002, δ2=0.001, ξ=0.015, α=0.9, p=0.467, σ=0.1, and 1≤Λ≤3. Figure 3(a) represents the forward bifurcation scenario where if R0<1 the heroin-free equilibrium is globally asymptotically stable while when R0>1 the heroin epidemic can persist. However, as we note from Figure 3(b), increasing parameter ω from ω=0.11 to ω=0.25 such that ω>ωc, a heroin epidemic can persist once established for a range of R0 values that are below unity which indicates the occurrence of backward bifurcation. This implies that reducing R0 below one will not necessarily be sufficient for eradication of heroin usage from the community. If R0 is sufficiently decreased such that R0<R0C the positive equilibrium no longer exists and heroin usage will cease to thrive and will eventually fall from its relatively high endemic level to the heroin-free equilibrium. From Figure 3(b) we note that when R0C≤R0≤1 there are a stable endemic equilibrium, an unstable endemic equilibrium, and a stable heroin-free equilibrium. When R0>1 there is only one stable endemic equilibrium. Figure 3(c) shows the effect of the saturation parameter ω on R0C; namely, increasing ω decreases R0C. Figure 3(d) shows that increasing the treatment rate α increases R0C which epidemiologically implies that high cure rates of heroin users can lead to shrinking of the backward bifurcation regime.


## 3.1. Endemic Equilibria

Within the context of our heroin model the endemic equilibrium refers to a state when heroin addiction is maintained over long time scales in the population. ThereforeS∗,U3∗,U2∗,U1∗>0 hold. To obtain the endemic equilibria E∗=(S∗,U3∗,U2∗,U1∗) we set (1) to zero and solve for the equilibrium quantities S∗,U3∗, and U2∗ in terms of U1∗. That is,(15)S∗=ΛβU1∗+μ,U3∗=σαU1∗+ξp+σ+δ2+μ1+ωU1∗U1∗μp+σ+δ2+μ1+ωU1∗,U2∗=αU1∗p+σ+δ2+μ1+ωU1∗.Substituting (15) into the second equation of (1) yields(16)fU1∗=AU1∗2+BU1∗+C=0,where (17)A=μ+δ1+ξp+σ+δ2+μβω,B=μ+δ1+ξp+σ+δ2+μβ+μωμ+δ1+ξp+σ+δ2+μ+ασ+δ2+μβ-βΛωp+σ+δ2+μ,C=μ+δ1+ξp+σ+δ2+μμ+αμσ+δ2+μ1-R0.The quadratic equation (16) can be analysed to investigate the existence of multiple equilibria when the basic reproduction number is below unity.If the parameter that accounts for the extent of saturation of heroin users in model (1) is excluded, that is, ω=0, (16) reduces to a linear equation (18)U1∗B~+C=0,where (19)B~=μ+δ1+ξp+σ+δ2+μβ+ασ+δ2+μβso that model (1) has the unique solution (20)U1∗=-CB~which is nonnegative if and only if R0>1. Hence, if ω=0, model (1) has a unique endemic equilibrium whenever R0>1 and this equilibrium approaches zero as R0 tends to one (R0→1+) because C→0. But there are no positive endemic equilibria if R0<1. These results are summarized in the following lemma.Lemma 5.
The epidemic model (1) when ω=0 has a unique positive endemic equilibrium E1∗=(S∗,U3∗,U2∗,U1∗) whenever R0>1 and no positive endemic equilibrium otherwise.In what follows, we investigate the global stability for both the HFE and the unique endemic equilibriumE1∗ for the case ω=0.


## 3.2. Global Stability for Heroin-Free Equilibrium Whenω=0

To investigate global stability we apply the method presented by Castillo-Chavez et al. [34]. First let X=(S,U3) and Y=(U2,U1) with X∈R2 representing the number of individuals not using heroin and Y∈R2 representing the number of individuals using heroin (i.e., heroin users in treatment and heroin users not in treatment). Now suppose (21)X′=FX,Y,Y′=GX,Y,GX,0=0,where X′ and Y′ denote differentiation with respect to time. The HFE is now denoted by H0=(X0,0), where X0=(S0,0). The following conditions (H1) and (H2) have to be met to guarantee a local asymptotic stability: H1
ForX′=FX,0,  X0is globally asymptotically stable (g.a.s).(H2)
G(X,Y)=BY-G^(X,Y), where G^(X,Y)≥0, for (X,Y)∈℧.B
=
D
Y
G
(
X
,
0
) and ℧ is the region where model (1) is biologically realistic. Then, Castillo-Chavez and Song [35] have shown that the following lemma is satisfied.Lemma 6.
The fixed pointH0=(X0,0,0) is g.a.s equilibrium of model (1) provided that R0<1 (locally asymptotically stable) and that assumptions (H1) and (H2) hold.Now consider the following theorem.Theorem 7.
SupposeR0<1. Then the HFE D0 is g.a.s.Proof.
LetX=(S,U3) and Y=(U2,U1), and H0=(X0,0,0), where X0=Λ/μ,0.
Then we have(22)X′=FX,Y=Λ-βU1S-μSσU2+ξU1-μU3.It is straightforward to see that, at the heroin-free equilibrium (HFE) S=S0=Λ/μ, F(X,0)=00. Thus, (23)X′=FX,0=Λ-μS-μU3.
Now, ast→∞, X→X0. Hence,X0 is globally asymptotically stable (i.e., condition (H1) is satisfied).
Now consider(24)GX,Y=-p+σ+δ2+μαpβS0-μ+δ1+ξ+αU2U1-0βU1S0-S so that (25)B=-p+σ+δ2+μαpβS0-μ+δ1+ξ,G^X,Y=0βU1S0-S. Since the total population is bounded by N=S+U3+U2+U1≤Λ/μ, we have S≤N≤Λ/μ. Thus, G^(X,Y)≥0, which now implies that conditions (H1) and (H2) are satisfied. Consequently by Lemma 6 the fixed point H0 is globally asymptotically stable when R0<1, which indicates nonexistence of multiple nontrivial equilibria when ω=0. The epidemiological implication of HFE being g.a.s is that heroin epidemic will be eliminated from the community if the threshold quantity R0 is decreased to (or maintained at) a value below unity.Now forω>0 we establish the following theorem.Theorem 8.
Forω>0 model (1) has (i)
a unique positive endemic equilibrium ifB<0 and either C=0 or B2-4AC=0,(ii)
a unique positive endemic equilibrium ifC<0 (i.e., R0>1) and B<0,(iii)
two positive endemic equilibria ifC>0,B<0, and B2-4AC>0,(iv)
no positive endemic equilibrium ifB>0 and either C>0 or B2<4AC.The theorem may be proved as follows. It is obvious to note that in quadratic equation (16) A is always positive and C is either positive or negative depending on whether the basic reproduction number is less than or greater than one, respectively.For Case (i) whereB<0 and C=0 (i.e., R0=1) (16) becomes linear AU1∗+B=0 and has a unique nonzero solution U1∗=-B/A which is positive if B<0 and negative if B>0. Referring to (15) we see that if U1∗ is unique then so are S∗,U2∗, and U3∗.For Case (ii) whereC<0 (i.e., R0>1) and B<0, (16) is quadratic and according to Descartes Rule of Signs (see [36]), (16) has one change of signs indicating (16) has a unique positive root and therefore there is a unique endemic equilibrium.In Case (iii) whereB<0 there is a nonnegative endemic equilibrium at R0=1. However, because (16) is quadratic and since the equilibrium is continuously determined by R0 then there must be an interval to the left of R0=1 on which two nonnegative equilibria coexist. That is,(26)U1,1∗=-B-B2-4AC2A,U1,2∗=-B+B2-4AC2A.For Case (iv) where B>0 and C>0 or B2<4AC, (16) has no positive real root as can be seen in (26), implying nonexistence of a positive endemic equilibrium.Case (iii) suggests that model (1) exhibits the phenomenon of backward bifurcation since the classical requirement for the occurrence of the phenomenon of backward bifurcation is satisfied, that is, the existence of multiple equilibria when the basic reproduction number is less than one. Thus, we have the following lemma.Theorem 9.
Model (1) has backward bifurcation at R0=1 if and only if B<0 (i.e., ω>ωc).Proof.
Consider (16), f(U1∗)=AU1∗2+BU1∗+C=0. Note that, at R0=1, C=0 implies that the graph f(U1∗) passes through the origin. If B<0 it follows that f(U1∗)=0 has a nonnegative root. Since f(U1∗) is a continuous function of C, if we increase C such that C>0, there is some open interval of C say (0,ψ) on which f(U1∗)=0 has two nonnegative roots. That is, there exist two nonnegative endemic equilibria when R0<1. This is indeed true since Case (iv) of Theorem 8 has already shown that for B≥0 model (1) does not have positive real roots when R0<1. Note that, at R0=1, C=0 the following equality holds:(27)μ+δ1+ξp+σ+δ2+μμ+αμμ+δ2+σ=p+σ+δ2+μβΛ.This together with condition B<0 implies that(28)ω>μ+δ1+ξp+σ+δ2+μβ+ασ+δ2+μβμασ+δ2+μ≜ωc.Thus, the phenomenon of backward bifurcation (referring to Case (iii), a situation where there are two endemic equilibria) occurs at the left ofR0=1 if and only if condition (28) is satisfied. This suggests that backward bifurcation will only occur if the parameter ω that accounts for the extent of saturation of heroin users exceeds a certain threshold (i.e., ω>ωc). However, if ω<ωc backward bifurcation cannot occur. Thus, parameter ω plays a critical role in the formation of backward bifurcation for model (1). It is instructive to note that similar results as the one shown in inequality (28) can be obtained by center manifold theory (see Appendix B), where it is emphasized that if ω>ωc the bifurcation coefficient a is positive indicating that the model system (1) undergoes the phenomenon of backward bifurcation. The epidemiological implication of backward bifurcation is that although it is necessary to reduce the basic reproduction number below one it is not sufficient to eradicate a heroin epidemic; rather R0 should be reduced further below a certain threshold which we shall denote by R0C (see Figure 3(b)).Figure 3
(a) and (b) represent bifurcations where drug users,U1∗, not in treatment equilibrium are plotted as a function of R0. The blue solid line represents the stable equilibrium while the red solid line represents the unstable equilibrium. (a) represents forward bifurcation and the parameters used include ω=0.11<ωc=0.1156, μ=0.01, β=0.001, δ1=0.002, δ2=0.001, ξ=0.015, α=0.9, p=0.467, σ=0.1, and Λ∈{1,3} while (b) represents backward bifurcation and parameters used are the same as in (a) except ω=0.25>ωc=0.1156. (c) depicts the critical value R0C as a function of saturation parameter ω. The black dotted line represents the threshold ωc which if exceeded gives rise to backward bifurcation. (d) depicts the critical value R0C as a function of treatment rate α.
(a)
(b)
(c)
(d)


## 3.3. Computation of New Threshold for Heroin EradicationR0C

Here we compute the value of the basic reproduction number where the two nontrivial endemic equilibria (both stable and unstable) collide and annihilate each other leaving only the heroin-free equilibrium point as the stationary solution. This isR0C in Figure 3(b). We choose Λ as the parameter of backward bifurcation. As we have noted in Case (iii) of Theorem 8, (16) has nonnegative roots corresponding to two endemic equilibria if and only if C>0 (i.e., R0<1) and B<0, B2>4AC. It follows that if B=-2AC, (16) has one nonnegative root -B/2A. Supposing there is backward bifurcation at R0=1, then there are two endemic equilibria for an interval of values of basic reproduction number starting from a threshold R0C defined by B=-2AC to a point where R0=1. To obtain this threshold R0C which is often referred to as the critical value of the basic reproduction number we replace the values of A, B, and C into the equality B2=4AC to obtain a quadratic equation in terms of Λ. For mathematical tractability we redefine coefficients B and C as B=B1-B2Λ and C=C1-C2Λ, where(29)B1=μ+δ1+ξp+σ+δ2+μβ+μωμ+δ1+ξp+σ+δ2+μ+ασ+δ2+μβB2=βωp+σ+δ2+μ,C1=μμ+δ1+ξp+σ+δ2+μ+αμσ+δ2+μ,C2=βp+σ+δ2+μ,B1,B2,C1,C2>0.Note that A>0 remains as previously defined in (16). Now the quadratic equation in terms of Λ can be obtained as (30)B22Λ2+4AC2-2B1B2Λ+B12-4AC1=0 Since we are considering the scenario where B<0 and R0<1, we have C1B2>C2B1 and thus just the single solution (31)Λc=B1B2-2AC2+2A2C22+AB2C1B2-C2B1B22>0.Thus, the critical value of basic reproduction number (i.e., the new threshold for heroin eradication), R0C, is given as (32)R0C=βΛcp+δ2+μ+σμαμ+δ2+σ+μμ+δ1+ξp+δ2+μ+σ.Consequently, from the above analysis of computation of threshold for heroin eradication we can deduce the following lemma.Lemma 10.
(a) IfR0>1, then model (1) has a unique endemic equilibrium point E∗. In this case heroin epidemic will persist in the community.
(b) IfR0C<R0<1, then model (1) has two endemic equilibria E¯1 and E¯2 and signals that model (1) has backward bifurcation.
(c) IfR0<R0C<1, then model (1) has only the heroin-free equilibrium point D0 and in this case heroin users will disappear.Figure3 exhibits typical bifurcation diagrams for model (1). To obtain the graphs we vary recruitment rate Λ while other parameter values are held fixed. The parameters used for the numerical simulation that leads to Figure 3(a) include ω=0.11, μ=0.01, β=0.001, δ1=0.002, δ2=0.001, ξ=0.015, α=0.9, p=0.467, σ=0.1, and 1≤Λ≤3. Figure 3(a) represents the forward bifurcation scenario where if R0<1 the heroin-free equilibrium is globally asymptotically stable while when R0>1 the heroin epidemic can persist. However, as we note from Figure 3(b), increasing parameter ω from ω=0.11 to ω=0.25 such that ω>ωc, a heroin epidemic can persist once established for a range of R0 values that are below unity which indicates the occurrence of backward bifurcation. This implies that reducing R0 below one will not necessarily be sufficient for eradication of heroin usage from the community. If R0 is sufficiently decreased such that R0<R0C the positive equilibrium no longer exists and heroin usage will cease to thrive and will eventually fall from its relatively high endemic level to the heroin-free equilibrium. From Figure 3(b) we note that when R0C≤R0≤1 there are a stable endemic equilibrium, an unstable endemic equilibrium, and a stable heroin-free equilibrium. When R0>1 there is only one stable endemic equilibrium. Figure 3(c) shows the effect of the saturation parameter ω on R0C; namely, increasing ω decreases R0C. Figure 3(d) shows that increasing the treatment rate α increases R0C which epidemiologically implies that high cure rates of heroin users can lead to shrinking of the backward bifurcation regime.


## 4. Global Stability

According to Theorem8, model (1) may have multiple equilibria when R0<1 and a unique endemic equilibrium whenever R0>1. First, global stability of the endemic equilibrium of model (1) is investigated for a special case, that is, when ω=ξ=σ=0, using Lyapunov direct method, and later proven for the general model (i.e., ω,σ,ξ>0) using a geometric approach.


### 4.1. Global Stability of Endemic EquilibriaE1∗ Using Lyapunov Method (Special Case ω=ξ=σ=0)

Lyapunov functions have previously been used in proving global stability of epidemic models; for instance, see [37–40] and the references therein. Now consider the following theorem.Theorem 11.
Ifω=ξ=σ=0 the unique endemic equilibrium E1∗ of model (1) is globally asymptotically stable in the interior of Ω if R0>1.Proof.
Defining the following Lyapunov candidate function:(33)WS,U1,U2=S-S∗22S∗+U1-U1∗-U1∗ln⁡U1U1∗+pU2∗αU1∗U2-U2∗-U2∗ln⁡U2U2∗.Now computing the time derivative of W(S,U1,U2) along the solutions of system (1) results in(34)W′S,U1,U2=S-S∗S∗dSdt+U1-U1∗U1dU1dt+pU2∗αU1∗1-U2∗U2dU2dt=S-S∗S∗Λ-βU1S-μS+U1-U1∗U1βU1S+pU2-μ+δ1+αU1+pU2∗αU1∗1-U2∗U2αU1-p+δ2+μU2.Because (S∗,U2∗,U1∗) is an endemic steady point of model (1) when ω=ξ=σ=0, then it follows that(35)Λ=βU1∗S∗+μS∗,p+δ2+μ=αU1∗U2∗,μ+δ1+α=βS∗+pU2∗U1∗.Using (35) in (34) yields(36)W′S,U1,U2=-S-S∗S∗μS-S∗+βSU1-S∗U1∗+U1-U1∗βS-S∗+pU2U1-U2∗U1∗+pU2∗αU1∗1-U2∗U2αU1-αU2U1∗U2∗. Note that (37)SU1-S∗U1∗=S∗U1-U1∗+U1S-S∗.Replacing the above equality in (36) results in(38)W′S,U1,U2=-S-S∗2S∗μ+βU1+pU2∗2-U1∗U2U1U2∗-U2∗U1U2U1∗=-S-S∗2S∗μ+βU1-pU2∗U1∗U2U1U2∗-U2∗U1U2U1∗2.
Hence,W′(S,U1,U2)≤0 for all S,U1,U2>0. Hence, the heroin endemic equilibrium E1∗ is stable and W′(S,U1,U2)=0 if and only if S=S∗,U2=U2∗,U1=U1∗. The largest compact invariant set when ω=ξ=σ=0 in {(S∗,U1∗,U2∗)∈R+3} is the singleton {E1∗}. Therefore, by LaSalle’s invariance principle, the endemic steady state E1∗ is globally asymptotically stable in the interior of R+3.The previous global stability analysis was only relevant for very specific case. In the subsequent subsection we use the geometric approach by Li and Muldowney [11, 13, 41] to obtain sufficient condition that ensures that the unique endemic equilibrium is globally asymptotically stable for a wide range of parameter values.


### 4.2. A Geometric Approach to Global Stability

For the general model, global stability is investigated using the Li and Muldowney [11, 13, 41] generalizations of the Poincaré-Bendixson approach for systems of n>2 ordinary differential equations. This criterion is sometimes referred to as ageometric approach to global stability [14, 42].To apply the geometric approach on model (1), consider the autonomous dynamical system dy/dt=f(y), where f=(f1,f2,f3,f4)T and f1,f2,f3,f4 represent the right-hand side of system (1), respectively. We first briefly outline the general mathematical framework of the procedure developed in Li and Muldowney [13, 16].Suppose that the mapy↦f(y) is a C1 function for y in an open subset D⊂Rn and consider the following autonomous dynamical system:(39)y′=fy.Let y(t,y0) be the solution to (39) satisfying y(0,y0)=y0. Now we make the following basic assumptions: (H3)
D is simply connected.(H4)
There exists a compact absorbing setK⊂D.(H5)
Equation (39) has a unique equilibrium y∗ in D.Now under the stated assumptions(H3)–(H5), y∗ is said to be globally stable in D if it is locally stable and all trajectories in D converge to the same equilibrium y∗. That is, system (39) has no nonconstant periodic solutions. It is important to mention that global stability can be tested byBendixson criteria. For n≥2 aBendixson criterion refers to a condition satisfied by field f which precludes the existence of nonconstant periodic solutions of (39). When n=2, (i.e., the planar case) the classical results (Poincaré-Bendixson theorem and Dulac criteria; see [43]) adequately provide such global conditions. For n≥3 a remarkable approach for proving global stability is given by the work of Li and Muldowney [11, 13, 16]. They showed that if conditions (H3)–(H5) hold and differential equation (39) fulfils a Bendixson criterion that is robust under C1 local ε-perturbations (a function g∈C1(D→Rn) is called a C1 local ε-perturbation of f at y0∈D if there exists an open neighbourhood U of y0 in D such that the support supp(f-g)⊂U and |f-g|C1<ε, where |f-g|C1=sup{|f(y)-g(y)|+|fy(y)-gy(y)|:y∈D}) of f at all nonequilibrium nonwandering (a point y0∈D is said to be nonwandering for system (39) if for any neighbourhood U of y0 in D there exists arbitrary large t such that U∩y(t,U)≠∅. As an example, any equilibrium, alpha limit point, or omega limit point is nonwandering) points for system (39), then y∗ is globally stable in D provided that it is stable. We now state the new Bendixson criterion based on the use of the Lozinskiĭ measure as developed in [13]. Consider the differential equation (39) under the stated assumptions (H3)–(H5). Let y↦P(y)(n2)×(n2) be a matrix-valued function which is C1 for y∈D and consider(40)A=PfP-1+PJ2P-1,where Pf is the directional derivative of P in the direction of the vector field f in system (39) and it is defined as(41)pi,jyf=∂pi,jy∂xT·fy=∇pi,j·fy,andJ[2] represents the second additive compound matrix J (where J(y)=Df(y)). In [44] the relation of compound matrices to differential equations is established. It is shown that, for an arbitrary n×n matrix J=(Ji,j), J[2] is an (n2)×(n2) matrix. Now define the following quantity:(42)q¯2=limsupt→∞supy0∈℧1t∫0tρAs,y0ds,where ρ(A) is the Lozinskiĭ measure of A with respect to vector norm |·| in RN, N=(n2), and ρ(A) is defined as(43)ρA=limh→0+1+hA-1h (see [45, 46]). In paper [13] it is proved that if conditions (H3) and (H4) are satisfied then q¯2<0, indicating that there are no orbits giving rise to simple closed rectifiable curve in D that is invariant for system (39) (i.e., periodic orbits, homoclinic orbits, and heteroclinic cycles). Furthermore, it has been demonstrated in [13] that, under the stated assumptions (H3)–(H5), quantity q¯2<0 implies the local stability of equilibrium point y∗. As a result the following theorem is true.Theorem 12 (see [13]).
Assuming that conditions(H3)–(H5) hold, then, the equilibrium point y∗ is globally asymptotically stable in D if a function P(y) and a Lozinskiĭ measure ρ exist such that quantity q¯2<0.Observe that wheneverR0>1, there exists a unique and positive endemic equilibrium E∗ (see Lemma 10) for model system (1). The method outlined above requires that (i) the endemic equilibrium E∗ is unique in the interior of ℧ (i.e., condition (H5) holds) and (ii) in the interior of ℧ there exists an absorbing compact set (condition (H4) holds). The heroin model studied here with the assumption that R0>1 fulfils conditions (H4)-(H5). It is easy to prove that when R0>1, the heroin-free equilibrium D0 is unstable (see Theorem 3). The instability of the heroin-free equilibrium D0 combined with D0∈δ℧ signals uniform persistence [47]. That is, there exists a positive constant c0>0 such that for every solution (S(t),U1(t),U2(t),U3(t)) of system (1) with (S(0),U1(0),U2(0),U3(0)) in the interior of biologically feasible region ℧ satisfies (44)limt→∞infSt,U1t,U2t,U3t≥c0. Because of boundedness of the region ℧, uniform persistence is equivalent to the existence of a compact set in the interior of ℧ which is absorbing for (1) (see [48]). Hence, condition (H4) is satisfied. Also it is shown that whenever R0>1 the model system (1) has only one equilibrium E∗ in the interior of ℧, so that condition (H5) is verified. Now for the heroin model system (1) the task involves verifying the Bendixson criterion (65). Note that the variable U3 does not affect first, second, and third equation of system (1). Thus, the fourth equation can be dropped from the analysis, and we only need to consider the following subsystem:(45)dSdt=Λ-βU1S-μS,dU1dt=βU1S+pU2-μ+δ1+ξU1-αU11+ωU1,dU2dt=αU11+ωU1-p+σ+δ2+μU2.The Jacobian matrix of subsystem (45) is found to be(46)J=-βU1+μ-βS0βU1βS-μ+δ1+ξ-α1+ωU12p0α1+ωU12-p+σ+δ2+μ.In working with Theorem 12 one needs to make use of additive compound matrices. For an arbitrary 3×3 matrix B, the second additive compound matrix B[2] is defined as(47)B=b11b12b13b21b22b23b31b32b33,B2=b11+b22b23-b13b32b11+b33b12-b31b21b22+b33.Thus, the second additive compound matrix of Jacobian matrix J of system (45) is given as(48)J2=J11p0α1+ωU12J22-βS0βU1J33,where (49)J11=-βU1-2μ+δ1+ξ-α1+ωU12+βS,J22=-βU1-2μ+p+σ+δ2,J33=-α1+ωU12-2μ+p+σ+δ1+δ2+ξ+βS.For the model system (45) a suitable vector norm · in R3 and a 3×3 matrix-valued function P(y) are given by(50)PS,U1,U2=1000U1U2000U1U2,ThusPfP-1=0000U1′U1-U2′U2000U1′U1-U2′U2,PJ2P-1=J11pU2U10αU11+ωU12U2J22-βS0βU1J33.Note that upper prime (′) denotes differentiation with respect to time, d./dt. Thus, A=PfP-1+PJ[2]P-1 can be obtained as(51)A=J11pU2U10αU11+ωU12U2J22+U1′U1-U2′U2-βS0βU1J33+U1′U1-U2′U2. It is helpful to write matrix A in block form as(52)A=A11A12A21A22,where(53)A11=-βU1-2μ+δ1+ξ-α1+ωU12+βS,A12=pU2U10,A21=αU11+ωU12U20T,A22=J22+U1′U1-U2′U2-βSβU1J33+U1′U1-U2′U2.Following [13], let (u,v,w) represent the vectors in R3≅R(32). Now for the norm |·| in R3 select (54)u,v,w=maxu,v+w and let ρ represent the Lozinskiĭ measure with respect to this norm. Applying the method of approximating ρ(A) as given in [46] leads to (55)ρA≤supg1,g2,where(56)g1=ρ1A11+A12g2=A21+ρ1A22.Here |A12| and |A21| are operator norms of A12 and A21 with respect to the l1 vector norm, where they are both regarded as mapping from R2 to R.  ρ1(A22) represents the Lozinskiĭ measure of the 2×2 matrix A22 with respect to the l1 norm in R2. To obtain ρ1(A22) we sum the absolute value of the off-diagonal elements to the diagonal one in each column of A22 and then take the maximum of two sums. Assuming that (1/2)δ1+ξ+α/(1+ωU1)2>βS, it follows that (57)ρ1A11=-βU1-2μ+δ1+ξ-α1+ωU12+βS,A12=maxpU2U1,0=pU2U1,A21=maxαU11+ωU12U2,0T=αU11+ωU12U2,ρ1A22=maxU1′U1-U2′U2-2μ+p+σ+δ2,U1′U1-U2′U2-2μ+p+σ+δ2-δ1-ξ-α1+ωU12+2βS=U1′U1-U2′U2-2μ+p+σ+δ2. Thus, g1 and g2 are, respectively, given as (58)g1=ρ1A11+A12=βS+pU2U1-βU1-α1+ωU12-2μ+δ1+ξ(59)g2=A21+ρ1A22=αU11+ωU12U2+U1′U1-U2′U2-2μ+p+σ+δ2.Now from second and third equation of (45) it is easy to obtain the following:(60)U1′U1=βS+pU2U1-α1+ωU1-μ+δ1+ξ(61)U2′U2=αU11+ωU1U2-p+σ+δ2+μ.Substituting (60) into (58) and (61) into (59), respectively, leads to (62)g1=U1′U1-μ-βU1+αωU11+ωU12≤U1′U1-μ+αωU11+ωU12,g2=U1′U1-μ-αωU121+ωU12U2≤U1′U1-μ.Now based on the definition of the method of approximating Lozinskiĭ measure of A, ρ(A) as given in [46], we now approximate the supremum of both g1 and g2. Hence, (63)ρA≤supg1,g2=supU1′U1-μ+αωU11+ωU12,U1′U1-μ=U1′U1-μ+supαωU11+ωU12,0≤U1′U1-μ+αωU11+ωU12.Thus we get the following inequality: (64)ρA≤U1′U1-μ+αωU11+ωU12.Now the next step involves substituting ρ(A) in(65)q¯2=limsupt→∞supy0∈℧1t∫0tρAs,y0dsand deducing whether q¯2<0. And if the inequality q¯2<0 does not hold we will need to establish a condition that leads to q¯2<0 being fulfilled.Considering uniform persistence, there existc0>0 and T>0 such that, for t>T, the following is implied:(66)St≥c0,U1t≥c0,U2t≥c0,U3t≥c0.Now by letting Γ1=αωc0/(1+ωc0)2 and Γ2=μ the following claim is made: if(67)Γ1<Γ2,then it follows that(68)ρA≤U1′U1-V~,where (69)V~=μ-αωc01+ωc02>0.Now, for t>T, it can be deduced that (70)q¯2=1t∫0tρAds=1t∫0TρAds+1t∫TtρAds≤1tlogU1tU1T+1t∫0TρAds-V~t-Tt<0whenV~>0, and, thus, the Bendixson criterion given by (42) is verified. However, it is important to observe that q¯2<0 if and only if condition (67) holds true. Thus, the following theorem is established.Theorem 13.
Provided thatR0>1, if Γ1<Γ2, then system (1) has a unique endemic equilibrium E∗ which is globally asymptotically stable with respect to solutions of (1) originating in the interior of ℧.The validity of Theorem13 will be shortly verified numerically.


## 4.1. Global Stability of Endemic EquilibriaE1∗ Using Lyapunov Method (Special Case ω=ξ=σ=0)

Lyapunov functions have previously been used in proving global stability of epidemic models; for instance, see [37–40] and the references therein. Now consider the following theorem.Theorem 11.
Ifω=ξ=σ=0 the unique endemic equilibrium E1∗ of model (1) is globally asymptotically stable in the interior of Ω if R0>1.Proof.
Defining the following Lyapunov candidate function:(33)WS,U1,U2=S-S∗22S∗+U1-U1∗-U1∗ln⁡U1U1∗+pU2∗αU1∗U2-U2∗-U2∗ln⁡U2U2∗.Now computing the time derivative of W(S,U1,U2) along the solutions of system (1) results in(34)W′S,U1,U2=S-S∗S∗dSdt+U1-U1∗U1dU1dt+pU2∗αU1∗1-U2∗U2dU2dt=S-S∗S∗Λ-βU1S-μS+U1-U1∗U1βU1S+pU2-μ+δ1+αU1+pU2∗αU1∗1-U2∗U2αU1-p+δ2+μU2.Because (S∗,U2∗,U1∗) is an endemic steady point of model (1) when ω=ξ=σ=0, then it follows that(35)Λ=βU1∗S∗+μS∗,p+δ2+μ=αU1∗U2∗,μ+δ1+α=βS∗+pU2∗U1∗.Using (35) in (34) yields(36)W′S,U1,U2=-S-S∗S∗μS-S∗+βSU1-S∗U1∗+U1-U1∗βS-S∗+pU2U1-U2∗U1∗+pU2∗αU1∗1-U2∗U2αU1-αU2U1∗U2∗. Note that (37)SU1-S∗U1∗=S∗U1-U1∗+U1S-S∗.Replacing the above equality in (36) results in(38)W′S,U1,U2=-S-S∗2S∗μ+βU1+pU2∗2-U1∗U2U1U2∗-U2∗U1U2U1∗=-S-S∗2S∗μ+βU1-pU2∗U1∗U2U1U2∗-U2∗U1U2U1∗2.
Hence,W′(S,U1,U2)≤0 for all S,U1,U2>0. Hence, the heroin endemic equilibrium E1∗ is stable and W′(S,U1,U2)=0 if and only if S=S∗,U2=U2∗,U1=U1∗. The largest compact invariant set when ω=ξ=σ=0 in {(S∗,U1∗,U2∗)∈R+3} is the singleton {E1∗}. Therefore, by LaSalle’s invariance principle, the endemic steady state E1∗ is globally asymptotically stable in the interior of R+3.The previous global stability analysis was only relevant for very specific case. In the subsequent subsection we use the geometric approach by Li and Muldowney [11, 13, 41] to obtain sufficient condition that ensures that the unique endemic equilibrium is globally asymptotically stable for a wide range of parameter values.


## 4.2. A Geometric Approach to Global Stability

For the general model, global stability is investigated using the Li and Muldowney [11, 13, 41] generalizations of the Poincaré-Bendixson approach for systems of n>2 ordinary differential equations. This criterion is sometimes referred to as ageometric approach to global stability [14, 42].To apply the geometric approach on model (1), consider the autonomous dynamical system dy/dt=f(y), where f=(f1,f2,f3,f4)T and f1,f2,f3,f4 represent the right-hand side of system (1), respectively. We first briefly outline the general mathematical framework of the procedure developed in Li and Muldowney [13, 16].Suppose that the mapy↦f(y) is a C1 function for y in an open subset D⊂Rn and consider the following autonomous dynamical system:(39)y′=fy.Let y(t,y0) be the solution to (39) satisfying y(0,y0)=y0. Now we make the following basic assumptions: (H3)
D is simply connected.(H4)
There exists a compact absorbing setK⊂D.(H5)
Equation (39) has a unique equilibrium y∗ in D.Now under the stated assumptions(H3)–(H5), y∗ is said to be globally stable in D if it is locally stable and all trajectories in D converge to the same equilibrium y∗. That is, system (39) has no nonconstant periodic solutions. It is important to mention that global stability can be tested byBendixson criteria. For n≥2 aBendixson criterion refers to a condition satisfied by field f which precludes the existence of nonconstant periodic solutions of (39). When n=2, (i.e., the planar case) the classical results (Poincaré-Bendixson theorem and Dulac criteria; see [43]) adequately provide such global conditions. For n≥3 a remarkable approach for proving global stability is given by the work of Li and Muldowney [11, 13, 16]. They showed that if conditions (H3)–(H5) hold and differential equation (39) fulfils a Bendixson criterion that is robust under C1 local ε-perturbations (a function g∈C1(D→Rn) is called a C1 local ε-perturbation of f at y0∈D if there exists an open neighbourhood U of y0 in D such that the support supp(f-g)⊂U and |f-g|C1<ε, where |f-g|C1=sup{|f(y)-g(y)|+|fy(y)-gy(y)|:y∈D}) of f at all nonequilibrium nonwandering (a point y0∈D is said to be nonwandering for system (39) if for any neighbourhood U of y0 in D there exists arbitrary large t such that U∩y(t,U)≠∅. As an example, any equilibrium, alpha limit point, or omega limit point is nonwandering) points for system (39), then y∗ is globally stable in D provided that it is stable. We now state the new Bendixson criterion based on the use of the Lozinskiĭ measure as developed in [13]. Consider the differential equation (39) under the stated assumptions (H3)–(H5). Let y↦P(y)(n2)×(n2) be a matrix-valued function which is C1 for y∈D and consider(40)A=PfP-1+PJ2P-1,where Pf is the directional derivative of P in the direction of the vector field f in system (39) and it is defined as(41)pi,jyf=∂pi,jy∂xT·fy=∇pi,j·fy,andJ[2] represents the second additive compound matrix J (where J(y)=Df(y)). In [44] the relation of compound matrices to differential equations is established. It is shown that, for an arbitrary n×n matrix J=(Ji,j), J[2] is an (n2)×(n2) matrix. Now define the following quantity:(42)q¯2=limsupt→∞supy0∈℧1t∫0tρAs,y0ds,where ρ(A) is the Lozinskiĭ measure of A with respect to vector norm |·| in RN, N=(n2), and ρ(A) is defined as(43)ρA=limh→0+1+hA-1h (see [45, 46]). In paper [13] it is proved that if conditions (H3) and (H4) are satisfied then q¯2<0, indicating that there are no orbits giving rise to simple closed rectifiable curve in D that is invariant for system (39) (i.e., periodic orbits, homoclinic orbits, and heteroclinic cycles). Furthermore, it has been demonstrated in [13] that, under the stated assumptions (H3)–(H5), quantity q¯2<0 implies the local stability of equilibrium point y∗. As a result the following theorem is true.Theorem 12 (see [13]).
Assuming that conditions(H3)–(H5) hold, then, the equilibrium point y∗ is globally asymptotically stable in D if a function P(y) and a Lozinskiĭ measure ρ exist such that quantity q¯2<0.Observe that wheneverR0>1, there exists a unique and positive endemic equilibrium E∗ (see Lemma 10) for model system (1). The method outlined above requires that (i) the endemic equilibrium E∗ is unique in the interior of ℧ (i.e., condition (H5) holds) and (ii) in the interior of ℧ there exists an absorbing compact set (condition (H4) holds). The heroin model studied here with the assumption that R0>1 fulfils conditions (H4)-(H5). It is easy to prove that when R0>1, the heroin-free equilibrium D0 is unstable (see Theorem 3). The instability of the heroin-free equilibrium D0 combined with D0∈δ℧ signals uniform persistence [47]. That is, there exists a positive constant c0>0 such that for every solution (S(t),U1(t),U2(t),U3(t)) of system (1) with (S(0),U1(0),U2(0),U3(0)) in the interior of biologically feasible region ℧ satisfies (44)limt→∞infSt,U1t,U2t,U3t≥c0. Because of boundedness of the region ℧, uniform persistence is equivalent to the existence of a compact set in the interior of ℧ which is absorbing for (1) (see [48]). Hence, condition (H4) is satisfied. Also it is shown that whenever R0>1 the model system (1) has only one equilibrium E∗ in the interior of ℧, so that condition (H5) is verified. Now for the heroin model system (1) the task involves verifying the Bendixson criterion (65). Note that the variable U3 does not affect first, second, and third equation of system (1). Thus, the fourth equation can be dropped from the analysis, and we only need to consider the following subsystem:(45)dSdt=Λ-βU1S-μS,dU1dt=βU1S+pU2-μ+δ1+ξU1-αU11+ωU1,dU2dt=αU11+ωU1-p+σ+δ2+μU2.The Jacobian matrix of subsystem (45) is found to be(46)J=-βU1+μ-βS0βU1βS-μ+δ1+ξ-α1+ωU12p0α1+ωU12-p+σ+δ2+μ.In working with Theorem 12 one needs to make use of additive compound matrices. For an arbitrary 3×3 matrix B, the second additive compound matrix B[2] is defined as(47)B=b11b12b13b21b22b23b31b32b33,B2=b11+b22b23-b13b32b11+b33b12-b31b21b22+b33.Thus, the second additive compound matrix of Jacobian matrix J of system (45) is given as(48)J2=J11p0α1+ωU12J22-βS0βU1J33,where (49)J11=-βU1-2μ+δ1+ξ-α1+ωU12+βS,J22=-βU1-2μ+p+σ+δ2,J33=-α1+ωU12-2μ+p+σ+δ1+δ2+ξ+βS.For the model system (45) a suitable vector norm · in R3 and a 3×3 matrix-valued function P(y) are given by(50)PS,U1,U2=1000U1U2000U1U2,ThusPfP-1=0000U1′U1-U2′U2000U1′U1-U2′U2,PJ2P-1=J11pU2U10αU11+ωU12U2J22-βS0βU1J33.Note that upper prime (′) denotes differentiation with respect to time, d./dt. Thus, A=PfP-1+PJ[2]P-1 can be obtained as(51)A=J11pU2U10αU11+ωU12U2J22+U1′U1-U2′U2-βS0βU1J33+U1′U1-U2′U2. It is helpful to write matrix A in block form as(52)A=A11A12A21A22,where(53)A11=-βU1-2μ+δ1+ξ-α1+ωU12+βS,A12=pU2U10,A21=αU11+ωU12U20T,A22=J22+U1′U1-U2′U2-βSβU1J33+U1′U1-U2′U2.Following [13], let (u,v,w) represent the vectors in R3≅R(32). Now for the norm |·| in R3 select (54)u,v,w=maxu,v+w and let ρ represent the Lozinskiĭ measure with respect to this norm. Applying the method of approximating ρ(A) as given in [46] leads to (55)ρA≤supg1,g2,where(56)g1=ρ1A11+A12g2=A21+ρ1A22.Here |A12| and |A21| are operator norms of A12 and A21 with respect to the l1 vector norm, where they are both regarded as mapping from R2 to R.  ρ1(A22) represents the Lozinskiĭ measure of the 2×2 matrix A22 with respect to the l1 norm in R2. To obtain ρ1(A22) we sum the absolute value of the off-diagonal elements to the diagonal one in each column of A22 and then take the maximum of two sums. Assuming that (1/2)δ1+ξ+α/(1+ωU1)2>βS, it follows that (57)ρ1A11=-βU1-2μ+δ1+ξ-α1+ωU12+βS,A12=maxpU2U1,0=pU2U1,A21=maxαU11+ωU12U2,0T=αU11+ωU12U2,ρ1A22=maxU1′U1-U2′U2-2μ+p+σ+δ2,U1′U1-U2′U2-2μ+p+σ+δ2-δ1-ξ-α1+ωU12+2βS=U1′U1-U2′U2-2μ+p+σ+δ2. Thus, g1 and g2 are, respectively, given as (58)g1=ρ1A11+A12=βS+pU2U1-βU1-α1+ωU12-2μ+δ1+ξ(59)g2=A21+ρ1A22=αU11+ωU12U2+U1′U1-U2′U2-2μ+p+σ+δ2.Now from second and third equation of (45) it is easy to obtain the following:(60)U1′U1=βS+pU2U1-α1+ωU1-μ+δ1+ξ(61)U2′U2=αU11+ωU1U2-p+σ+δ2+μ.Substituting (60) into (58) and (61) into (59), respectively, leads to (62)g1=U1′U1-μ-βU1+αωU11+ωU12≤U1′U1-μ+αωU11+ωU12,g2=U1′U1-μ-αωU121+ωU12U2≤U1′U1-μ.Now based on the definition of the method of approximating Lozinskiĭ measure of A, ρ(A) as given in [46], we now approximate the supremum of both g1 and g2. Hence, (63)ρA≤supg1,g2=supU1′U1-μ+αωU11+ωU12,U1′U1-μ=U1′U1-μ+supαωU11+ωU12,0≤U1′U1-μ+αωU11+ωU12.Thus we get the following inequality: (64)ρA≤U1′U1-μ+αωU11+ωU12.Now the next step involves substituting ρ(A) in(65)q¯2=limsupt→∞supy0∈℧1t∫0tρAs,y0dsand deducing whether q¯2<0. And if the inequality q¯2<0 does not hold we will need to establish a condition that leads to q¯2<0 being fulfilled.Considering uniform persistence, there existc0>0 and T>0 such that, for t>T, the following is implied:(66)St≥c0,U1t≥c0,U2t≥c0,U3t≥c0.Now by letting Γ1=αωc0/(1+ωc0)2 and Γ2=μ the following claim is made: if(67)Γ1<Γ2,then it follows that(68)ρA≤U1′U1-V~,where (69)V~=μ-αωc01+ωc02>0.Now, for t>T, it can be deduced that (70)q¯2=1t∫0tρAds=1t∫0TρAds+1t∫TtρAds≤1tlogU1tU1T+1t∫0TρAds-V~t-Tt<0whenV~>0, and, thus, the Bendixson criterion given by (42) is verified. However, it is important to observe that q¯2<0 if and only if condition (67) holds true. Thus, the following theorem is established.Theorem 13.
Provided thatR0>1, if Γ1<Γ2, then system (1) has a unique endemic equilibrium E∗ which is globally asymptotically stable with respect to solutions of (1) originating in the interior of ℧.The validity of Theorem13 will be shortly verified numerically.


## 5. Numerical Examples

In this section numerical simulations of the heroin epidemic model are presented to support theoretical findings. Figure4 which shows backward bifurcation is obtained by plotting heroin users equilibrium as a function of R0. The figures present a scenario where R0 is varied via parameter β (i.e., 0.0005≤β≤0.0015) and other parameters are fixed. Figures 4(a)–4(d) show that increasing the parameter ω leads to the expansion of the region of bistability while decreasing ω results in contraction of the bistability region. The heroin eradication threshold (also referred to as critical reproduction number) R0C shifts from right to left when ω increases and vice versa when ω decreases. High ω value implies not enough treatment for a large population of heroin users, thus favouring a situation where there will always be heroin users within the community even though R0<1.Figure 4
Illustration of the effect of increasing parameterω that accounts for the saturation of heroin users. Here the parameters remain as in caption of Figure 3 except Λ=2 while ω is shown in the figure. The heroin eradication thresholds (i.e., R0C) corresponding to (a)–(d) are now 0.5204, 0.6038, 0.7314, and 0.9131, respectively.
(a)
(b)
(c)
(d)Figures5(a) and 5(b) exhibit the time course of the heroin endemic in a parameter regime where there is a backward bifurcation. In both Figures 5(a) and 5(b)  R0=0.7506<1. The figures show the dependence of heroin usage on the size of the initial conditions supplied to the system, which is a common characteristic of models that have a bistability region. If the model is supplied with initial conditions that are below the unstable curve (see the red solid line in Figure 3(b)) the solution trajectories are attracted to the heroin-free equilibrium while if initial conditions are chosen such that they are above the unstable curve, then the solution trajectories are attracted to a stable nontrivial equilibrium. Thus, in the case there is backward bifurcation, the initial number of people engaging in heroin use governs the course of the heroin epidemic.Figure 5
(a) and (b) illustrate the dependence of the heroin epidemic extinction or persistence on the initial states provided to the model. Parameters used includeω=0.15, β=0.001, δ1=0.002, δ2=0.001, ξ=0.015, μ=0.01, α=0.9, p=0.467, σ=0.1, and Λ=1.5. With these parameters the corresponding R0=0.7506. (c) The figure illustrates the impact of increasing parameter ω (the extent of saturation of treatment) on heroin users when R0<1. Parameters used include β=0.001, δ1=0.002, δ2=0.001, ξ=0.015, μ=0.01, α=0.9, p=0.46, σ=0.1, and Λ=1.5 which correspond to R0=0.7427<1. (d) The figure shows the effect of treatment α on heroin users U1 when all other parameters and initial conditions are fixed.
(a)
(b)
(c)
(d)Figure5(c) shows the time course of the heroin users when the parameter ω that accounts for the extent of saturation of heroin users is varied while the initial states and all other parameter values are fixed to constant values. It can be seen that not all values of ω will trigger rapid growth towards an endemic equilibrium when R0<1. Indeed, parameter ω has to exceed a certain fixed threshold ωc, hence supporting our theoretical findings that a nonzero equilibrium when R0<1 can only be maintained when ω is greater than ωc; see (28). Figure 5(d) shows the effect of treatment α on heroin users. High treatment leads to a steady decline of heroin users.Figure6 presents a scenario where R0>1. In this scenario we expect that when a heroin user enters a heroin-free community there will be rapid growth of heroin users until a globally stable equilibrium point is reached. Recalling that parameter ω does not appear in R0, nevertheless it does affect the model dynamics. The impact of ω when R0>1 is different from the case where R0<1. For R0<1 it plays a key role in inducing bistability while for R0>1 the parameter ω impacts the heroin dynamics by determining the time taken for an epidemic to occur. For relatively high ω values there is a sudden decrease of susceptible subpopulation while for relatively low values of ω there is a gradual decrease of susceptible subpopulation. Moreover, Figure 6(b) depicts that for any given value of ω the heroin users gradually approach a stable endemic equilibrium point. The only striking difference is the time taken to reach the heroin endemic equilibrium. At high values of parameter ω the heroin endemic will rapidly approach an equilibrium.Figure 6
(a) and (b) show the effect of increasing the parameterω that accounts for the extent of saturation of treatment when R0>1. Parameters used are the same as in Figure 5(a) except p=0.46 and Λ=3 corresponding to R0=1.4855>1.
(a)
(b)We now verify the global stability condition obtained using the geometric approach based on the following parameter values:(i)
β=0.001,δ1=0.002,δ2=0.001,ξ=0.015,μ=0.01,α=0.9,p=0.467,σ=0.1,ω=2,Λ=3,c0=50. With these parameters the corresponding R0=1.5012>1 and Γ1=0.0088<Γ2=0.0100. In this scenario condition (67) is satisfied and the model system (1) should be globally asymptotically stable. Figures 7(a), 7(b), 7(c), and 7(d) show existence of an apparently stable equilibrium.(ii)
We now use the same set of parameter values as Case (i) exceptω=0.05. This leads to Γ1=0.1837>Γ2=0.0100. In this case the asymptotic stability condition is not satisfied and unsurprisingly model system (1) has periodic solutions as shown in Figures 7(e), 7(f), 7(g), and 7(h). The epidemiological interpretation of this is that the heroin epidemic will fluctuate between low and high endemic levels. The cycles are induced by time delays of the transmission processes.Figure 7
Illustration of validity of global stability condition. (a), (b), (c), and (d) represent a scenario where condition (67) holds true and global stability is predicted. (e), (f), (g), and (h) represent a scenario where condition (67) does not hold and oscillations are expected.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


### 5.1. Uncertainty and Sensitivity Analysis

Here we conduct sensitivity analysis so as to identify critical inputs of our heroin epidemic model and gain insights on how input uncertainty influences model outcome [49]. To achieve this we make use of the Latin hypercube sampling (LHS) technique which provides a comprehensive method of assessing model sensitivity to parameters over multidimensional parameter space. One of the advantages of using the LHS technique is that it requires fewer samples of parameters than simple random sampling to achieve the same accuracy (see [49] and the references therein for in-depth discussion on LHS). In our heroin epidemic model the LHS technique is important due to the relatively large uncertainty of the model parameter estimates we have used. The technique works in combination with the partial rank correlation coefficient (PRCC) which estimates the sign and strength of the relationship that exists between each model parameter and any specified output variable [50, 51]. The PRCC values are bounded between 1 and -1, with a PRCC value close to 1 (-1) indicating very strong positive (or negative) correlation. The relative importance of the model parameters can be directly evaluated by comparing the values of the PRCC [51]. The uncertainty and sensitivity analysis using the LHS technique involves first selecting a baseline value and a range for each parameter of the heroin epidemic model (1) (see Table 2) and then performing multiple runs for a given outcome variable or response function. To enhance accuracy, 1500 random samples of parameter values were used for the sensitivity analysis and significant levels are set for pvalue<0.05.Table 2
Parameter baseline values and ranges used in sensitivity analysis.
Parameter
Baseline value
Range
Source
Λ
2
1–5
Assumed
β
0.001
0.0005–0.015
[7]
μ
1
80
[0.01125, 0.01375]
[52]
p
0.467
0.1–0.8
[7]
α
0.5
0.2–0.95
[23]
δ1
0.002
0.0008–0.0025
Assumed
δ2
0.001
0.00095–0.002
Assumed
ξ
0.5
0.05–0.5
Assumed
σ
0.1
0.1–0.7
Assumed
ω
0.01
0.008–0.25
[23]Figure8 displays the sensitivity analysis results for the heroin users not in treatment U1(t). It is straightforward to see that recruitment rate Λ, effective contact rate β, relapsing rate of heroin users in treatment to heroin users not in treatment p, and saturation parameter ω are positively correlated while natural death μ, treatment rate α, heroin-induced death rates (δ1,δ2), “self-cure” rate ξ, and successful recovery rate of heroin users in treatment σ are negatively correlated. Among the positively correlated PRCC values the parameters β,p, and ω are strongly positively correlated to heroin users not in treatment as evidenced by their high PRCC values. However, at time point year 15 the effective contact rate β has a slightly higher PRCC value than relapsing rate p suggesting that at initial stage of a heroin epidemic effective contact between heroin users U1(t) and susceptibles significantly contributes to the emergence of a heroin epidemic. On the other hand, at time point year 30 the situation observed at time point year 15 is reversed. That is, relapsing of heroin users in treatment has slightly higher PRCC values than the effective contact rate β. Hence, long-term relapsing of heroin users in treatment back to the heroin users not in treatment also plays a role in ensuring that there will always be heroin users within the community. Thus, attempting to control heroin usage within the community measures that ensure that heroin users undergoing treatment do not relapse should be of great importance. The extent of saturation of heroin users as a result of failure to treat heroin users promptly which is accounted by parameter ω also contributes to sustaining heroin epidemic. As suggested by the strongly negatively correlated PRCC value of parameter σ, ensuring that heroin users in treatment are successfully treated (i.e., they do not relapse) can substantially reduce the subpopulation of heroin users. In general the sensitivity analysis results suggest that, to combat heroin epidemic, policy makers and clinicians should target controlling effective contact rate β, relapsing rate p, and the extent of saturation of heroin users rate ω parameters.Figure 8
The PRCC output for heroin usersU1(t) not in treatment.


## 5.1. Uncertainty and Sensitivity Analysis

Here we conduct sensitivity analysis so as to identify critical inputs of our heroin epidemic model and gain insights on how input uncertainty influences model outcome [49]. To achieve this we make use of the Latin hypercube sampling (LHS) technique which provides a comprehensive method of assessing model sensitivity to parameters over multidimensional parameter space. One of the advantages of using the LHS technique is that it requires fewer samples of parameters than simple random sampling to achieve the same accuracy (see [49] and the references therein for in-depth discussion on LHS). In our heroin epidemic model the LHS technique is important due to the relatively large uncertainty of the model parameter estimates we have used. The technique works in combination with the partial rank correlation coefficient (PRCC) which estimates the sign and strength of the relationship that exists between each model parameter and any specified output variable [50, 51]. The PRCC values are bounded between 1 and -1, with a PRCC value close to 1 (-1) indicating very strong positive (or negative) correlation. The relative importance of the model parameters can be directly evaluated by comparing the values of the PRCC [51]. The uncertainty and sensitivity analysis using the LHS technique involves first selecting a baseline value and a range for each parameter of the heroin epidemic model (1) (see Table 2) and then performing multiple runs for a given outcome variable or response function. To enhance accuracy, 1500 random samples of parameter values were used for the sensitivity analysis and significant levels are set for pvalue<0.05.Table 2
Parameter baseline values and ranges used in sensitivity analysis.
Parameter
Baseline value
Range
Source
Λ
2
1–5
Assumed
β
0.001
0.0005–0.015
[7]
μ
1
80
[0.01125, 0.01375]
[52]
p
0.467
0.1–0.8
[7]
α
0.5
0.2–0.95
[23]
δ1
0.002
0.0008–0.0025
Assumed
δ2
0.001
0.00095–0.002
Assumed
ξ
0.5
0.05–0.5
Assumed
σ
0.1
0.1–0.7
Assumed
ω
0.01
0.008–0.25
[23]Figure8 displays the sensitivity analysis results for the heroin users not in treatment U1(t). It is straightforward to see that recruitment rate Λ, effective contact rate β, relapsing rate of heroin users in treatment to heroin users not in treatment p, and saturation parameter ω are positively correlated while natural death μ, treatment rate α, heroin-induced death rates (δ1,δ2), “self-cure” rate ξ, and successful recovery rate of heroin users in treatment σ are negatively correlated. Among the positively correlated PRCC values the parameters β,p, and ω are strongly positively correlated to heroin users not in treatment as evidenced by their high PRCC values. However, at time point year 15 the effective contact rate β has a slightly higher PRCC value than relapsing rate p suggesting that at initial stage of a heroin epidemic effective contact between heroin users U1(t) and susceptibles significantly contributes to the emergence of a heroin epidemic. On the other hand, at time point year 30 the situation observed at time point year 15 is reversed. That is, relapsing of heroin users in treatment has slightly higher PRCC values than the effective contact rate β. Hence, long-term relapsing of heroin users in treatment back to the heroin users not in treatment also plays a role in ensuring that there will always be heroin users within the community. Thus, attempting to control heroin usage within the community measures that ensure that heroin users undergoing treatment do not relapse should be of great importance. The extent of saturation of heroin users as a result of failure to treat heroin users promptly which is accounted by parameter ω also contributes to sustaining heroin epidemic. As suggested by the strongly negatively correlated PRCC value of parameter σ, ensuring that heroin users in treatment are successfully treated (i.e., they do not relapse) can substantially reduce the subpopulation of heroin users. In general the sensitivity analysis results suggest that, to combat heroin epidemic, policy makers and clinicians should target controlling effective contact rate β, relapsing rate p, and the extent of saturation of heroin users rate ω parameters.Figure 8
The PRCC output for heroin usersU1(t) not in treatment.


## 6. Concluding Remarks

In this study we formulated a heroin epidemic model with bilinear and saturated treatment function. The threshold parameterR0 usually referred to as the basic reproduction number plays a key role in the prediction of disease persistence or extinction. Epidemiologically, when R0 exceeds one, an epidemic persists and if it is below unity the disease will die out. This classical viewpoint has recently been challenged by many researchers since it is not always true that a disease will disappear if R0 is decreased below one. In the present heroin epidemic model the analytical results indicate that R0 is indeed the threshold when the parameter ω=0. However, when a saturated treatment function (i.e., ω>0) rather than a linear treatment rate is used, the heroin model exhibits the phenomenon of backward bifurcation where a heroin-free equilibrium and two nontrivial equilibria coexist even though the basic reproduction number is below unity (see Theorem 8—Case (iii)). The appearance of backward bifurcation indicates that it is not sufficient to decrease the basic reproduction number below unity for the eradication of heroin users within the community. Thus, to effectively control the spread of heroin users one has to reduce R0 below another threshold referred to as the critical value of the basic reproduction number R0C. That is, heroin users can be eradicated if R0<R0C<1. It is important to note that although the parameter ω might be present in the model, not every value of ω will lead to bistability. Instead ω has to be greater than a certain threshold ωc which is an aggregate of model parameters (see (28)). In general both analytical (see Appendix A for center manifold) and numerical results suggest that the saturation parameter ω is the one responsible for backward bifurcation. Failure to intervene before heroin users have accumulated in the community will lead to a situation where a heroin epidemic exists even though basic reproduction number is below one. Improvement of existing medical technology as well as channelling sufficient resources in medicines can significantly facilitate early intervention by ensuring that heroin users receive treatment promptly.In addition global stability properties using both the Lyapunov direct method and geometric approach by Li and Muldowney have been investigated. We note that, even for a four-dimensional model, the use of the two nonlinear stability techniques becomes nontrivial. In fact when all the parameters of the model are accounted for, it is difficult if not impossible to design a Lyapunov function. Using the geometric approach we establish a global condition that accounts for all parameters that if satisfied signals that heroin persistence within the community is globally stable. However, if the global condition is not satisfied heroin users can oscillate periodically in number (see Figures7(e), 7(f), 7(g), and 7(h)). Moreover, sensitivity and uncertainty analysis using LHS results indicate that effective contact between susceptibles and heroin users β, the relapsing of heroin users in treatment p, and the extent of saturation of heroin users parameter ω are the ones which contribute to persistence of heroin epidemic within the community.



---

*Source: 1953036-2017-08-31.xml*