Analysis of a Heroin Epidemic Model with Saturated Treatment Function

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 amore complete representation of themodel 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.


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 by  0 , is one of the most important concepts, 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).
given its ability to predict the course of an epidemic.It will also prove invaluable in our study of heroin dynamics in society. 0 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  0 < 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  0 < 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.
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][7][8][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][12][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][15][16][17][18][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  ≥ 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.

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 (), individuals successfully treated from heroin use ( 3 ), heroin users undergoing treatment ( 2 ), and heroin users not in treatment ( 1 ).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 Figure 2, which can be written in terms of the following set of equations: In brief, the susceptible subpopulation () is generated at a constant rate through immigration and birth at rate Λ.Some susceptible individuals who come into contact with heroin users  1 () may begin to use heroin.Hence, the susceptible population is diminished due to contact with heroin users at rate  1 , while heroin users increase at the same rate.Heroin users also increase when those undergoing treatment relapse at rate  2 and return to their heroin using lifestyle.Heroin users reduce in number as a result of treatment which is represented by the treatment function ( 1 ).Moreover, the user subpopulation is reduced by heroin-induced death at rate  1  1 as well as a result of the self-decision to cease using heroin (also referred to as "self-cure") at rate  1 .Individuals undergoing treatment are diminished through the following processes: relapse to heroin using at rate  2 , heroin-induced death at rate  2  2 , and successful treatment at rate  2 .Finally, the recovered/successfully treated subpopulation  3 () 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 .
Heroin epidemic models studied to date [6][7][8][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  1 () receive treatment based on the following more general saturated treatment function: 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  1 the treatment function reduces to ( 1 ) ≈  1 while for large  1 it reduces to ( 1 ) ≈ / which actually characterizes the saturated phenomena of the treatment.Further, if  = 0, the treatment function becomes ( 1 ) =  1 which is the usual linear treatment rate.1/(1 +  1 ) 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.

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 variables ,  3 ,  2 ,  1 are nonnegative for all time  > 0.
Hence, the following result.
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 Considering that 0 <  1 () < (), 0 < 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 ,  3 ,  2 , and  1 are any solution of system (1) supplied with nonnegative initial conditions.Then it is straightforward to note that the total population  satisfies the inequality From (11) it follows that / ≤ Λ −  which implies / ≤ 0 if  ≥ Λ/.The standard comparison theorem [31] can be used to deduce that () ≤ (0) − + (Λ/)(1 −  − ).In particular () ≤ Λ/ if (0) ≤ Λ/ for all  > 0. Thus, under the flow induced by system (1), the region ℧ is positively invariant.Furthermore, for (0) > Λ/ the trajectory solutions () 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  0 = ( 0 , 0, 0, 0), a scenario representing a heroin-free state in the community. 0 = Λ/ represents the number of susceptibles when no one is using heroin.The basic reproduction number denoted by  0 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  0 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 is  0 = 1/( +  1 +  + ) and the probability of surviving this compartment and moving to the treatment compartment is  1 = /(+ 1 ++).Now the probability of surviving heroin users in treatment class and then returning to the heroin users class not in treatment is  2 = /( +  +  2 + ).Thus, the total average time spent by the heroin users not in the treatment compartment on multiple passes can be obtained as 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 Multiplying ( 13) with the effective contact rate  and the average recruitment rate Λ/ we obtain heroin basic reproduction number as It is easy to observe that  0 is inversely proportional to treatment , which implies that if treatment rate is maintained sufficiently high it can control a heroin epidemic (by reducing  0 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 provided  0 < 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  0 < 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.
Substituting ( 15) into the second equation of (1) yields where 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 where so that model (1) has the unique solution which is nonnegative if and only if  0 > 1.Hence, if  = 0, model (1) has a unique endemic equilibrium whenever  0 > 1 and this equilibrium approaches zero as  0 tends to one ( 0 → 1+) because  → 0. But there are no positive endemic equilibria if  0 < 1.These results are summarized in the following lemma.In what follows, we investigate the global stability for both the HFE and the unique endemic equilibrium  * 1 for the case  = 0.

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 = (,  3 ) and Y = ( 2 ,  1 ) with X ∈ R 2 representing the number of individuals not using heroin and Y ∈ R 2 representing the number of individuals using heroin (i.e., heroin users in treatment and heroin users not in treatment).Now suppose where X  and Y  denote differentiation with respect to time.The HFE is now denoted by  0 = (X 0 , 0), where X 0 = ( 0 , 0).The following conditions (1) and (2) have to be met to guarantee a local asymptotic stability: ) and ℧ is the region where model ( 1) is biologically realistic.Then, Castillo-Chavez and Song [35] have shown that the following lemma is satisfied.
Now consider the following theorem.
The epidemiological implication of HFE being g.a.s is that heroin epidemic will be eliminated from the community if the threshold quantity  0 is decreased to (or maintained at) a value below unity.
Now for  > 0 we establish the following theorem.
The theorem may be proved as follows.It is obvious to note that in quadratic equation ( 16)  is always positive and  is either positive or negative depending on whether the basic reproduction number is less than or greater than one, respectively.
In Case (iii) where  < 0 there is a nonnegative endemic equilibrium at  0 = 1.However, because ( 16) is quadratic and since the equilibrium is continuously determined by  0 then there must be an interval to the left of  0 = 1 on which two nonnegative equilibria coexist.That is, For Case (iv) where  > 0 and  > 0 or  2 < 4, (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  0 = 1 if and only if  < 0 (i.e.,  >   ).
Proof.Consider (16), ( * 1 ) =  * 2 1 +  * 1 +  = 0. Note that, at  0 = 1,  = 0 implies that the graph ( * 1 ) passes through the origin.If  < 0 it follows that ( * 1 ) = 0 has a nonnegative root.Since ( * 1 ) is a continuous function of , if we increase  such that  > 0, there is some open interval of  say (0, ) on which ( * 1 ) = 0 has two nonnegative roots.That is, there exist two nonnegative endemic equilibria when  0 < 1.This is indeed true since Case (iv) of Theorem 8 has already shown that for  ≥ 0 model (1) does not have positive real roots when  0 < 1.Note that, at  0 = 1,  = 0 the following equality holds: This together with condition  < 0 implies that Thus, the phenomenon of backward bifurcation (referring to Case (iii), a situation where there are two endemic equilibria) occurs at the left of  0 = 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.,  >   ).However, if  <   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  >   the bifurcation coefficient  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  0 should be reduced further below a certain threshold which we shall denote by   0 (see Figure 3(b)).

Computation of New
Threshold for Heroin Eradication   0 .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 is   0 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  > 0 (i.e.,  0 < 1) and  < 0,  2 > 4.It follows that if  = −2 √ , (16) has one nonnegative root −/2.Supposing there is backward bifurcation at  0 = 1, then there are two endemic equilibria for an interval of values of basic reproduction number starting from a threshold   0 defined by  = −2 √  to a point where  0 = 1.To obtain this threshold   0 which is often referred to as the critical value of the basic reproduction number we replace the values of , , and  into the equality  2 = 4 to obtain a quadratic equation in terms of Λ.For mathematical tractability we redefine coefficients  and  as  =  1 −  2 Λ and  =  1 −  2 Λ, where Note that  > 0 remains as previously defined in (16).Now the quadratic equation in terms of Λ can be obtained as Since we are considering the scenario where  < 0 and  0 < 1, we have  1  2 >  2  1 and thus just the single solution Thus, the critical value of basic reproduction number (i.e., the new threshold for heroin eradication),   0 , is given as Consequently, from the above analysis of computation of threshold for heroin eradication we can deduce the following lemma.(c) If  0 <   0 < 1, then model ( 1) has only the heroinfree equilibrium point  0 and in this case heroin users will disappear.
Figure 3 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,  = 0.467,  = 0.1, and 1 ≤ Λ ≤ 3. Figure 3(a) represents the forward bifurcation scenario where if  0 < 1 the heroin-free equilibrium is globally asymptotically stable while when  0 > 1 the heroin epidemic can persist.However, as we note from Figure 3(b), increasing parameter  from  = 0.11 to  = 0.25 such that  >   , a heroin epidemic can persist once established for a range of  0 values that are below unity which indicates the occurrence of backward bifurcation.This implies that reducing  0 below one will not necessarily be sufficient for eradication of heroin usage from the community.If  0 is sufficiently decreased such that  0 <   0 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   0 ≤  0 ≤ 1 there are a stable endemic equilibrium, an unstable endemic equilibrium, and a stable heroin-free equilibrium.When  0 > 1 there is only one stable endemic equilibrium.Figure 3(c) shows the effect of the saturation parameter  on   0 ; namely, increasing  decreases   0 .Figure 3(d) shows that increasing the treatment rate  increases   0 which epidemiologically implies that high cure rates of heroin users can lead to shrinking of the backward bifurcation regime.

Global Stability
According to Theorem 8, model (1) may have multiple equilibria when  0 < 1 and a unique endemic equilibrium whenever  0 > 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.

Global Stability of Endemic Equilibria 𝐸 *
1 Using Lyapunov Method (Special Case  =  =  = 0).Lyapunov functions have previously been used in proving global stability of epidemic models; for instance, see [37][38][39][40] and the references therein.Now consider the following theorem.

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  > 2 ordinary differential equations.This criterion is sometimes referred to as a geometric approach to global stability [14,42].
Suppose that the map   → () is a  1 function for  in an open subset  ⊂ R  and consider the following autonomous dynamical system: Let (,  0 ) be the solution to (39) satisfying (0,  0 ) =  0 .Now we make the following basic assumptions: (3)  is simply connected.
Now under the stated assumptions (3)-(5),  * is said to be globally stable in  if it is locally stable and all trajectories in  converge to the same equilibrium  * .That is, system (39) has no nonconstant periodic solutions.It is important to mention that global stability can be tested by Bendixson criteria.For  ≥ 2 a Bendixson criterion refers to a condition satisfied by field  which precludes the existence of nonconstant periodic solutions of (39).When  = 2, (i.e., the planar case) the classical results (Poincaré-Bendixson theorem and Dulac criteria; see [43]) adequately provide such global conditions.For  ≥ 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 (3)-(5) hold and differential equation ( 39) fulfils a Bendixson criterion that is robust under   ∈ }) of  at all nonequilibrium nonwandering (a point  0 ∈  is said to be nonwandering for system (39) if for any neighbourhood  of  0 in  there exists arbitrary large  such that  ∩ (, ) ̸ = 0.As an example, any equilibrium, alpha limit point, or omega limit point is nonwandering) points for system (39), then  * is globally stable in  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 (3)-(5).Let   → () (  2 )×(  2 ) be a matrixvalued function which is  1 for  ∈  and consider where   is the directional derivative of  in the direction of the vector field  in system (39) and it is defined as and  [2] represents the second additive compound matrix  (where () = ()).In [44] the relation of compound matrices to differential equations is established.It is shown that, for an arbitrary  ×  matrix  = ( , ),  [2] is an (  2 ) × (  2 ) matrix.Now define the following quantity: where () is the Lozinskiȋ measure of  with respect to vector norm | ⋅ | in R  ,  = (  2 ), and () is defined as (see [45,46]).In paper [13] it is proved that if conditions (3) and (4) are satisfied then  2 < 0, indicating that there are no orbits giving rise to simple closed rectifiable curve in  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 (3)-(5), quantity  2 < 0 implies the local stability of equilibrium point  * .As a result the following theorem is true.
Observe that whenever  0 > 1, there exists a unique and positive endemic equilibrium  * (see Lemma 10) for model system (1).The method outlined above requires that (i) the endemic equilibrium  * is unique in the interior of ℧ (i.e., condition (5) holds) and (ii) in the interior of ℧ there exists an absorbing compact set (condition (4) holds).The heroin model studied here with the assumption that  0 > 1 fulfils conditions (4)-(5).It is easy to prove that when  0 > 1, the heroin-free equilibrium  0 is unstable (see Theorem 3).The instability of the heroin-free equilibrium  0 combined with  0 ∈ ℧ signals uniform persistence [47].That is, there exists a positive constant  0 > 0 such that for every solution ((),  1 (),  2 (),  3 ()) of system ( 1) with ((0),  1 (0),  2 (0),  3 (0)) in the interior of biologically feasible region ℧ satisfies lim 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 (4) is satisfied.Also it is shown that whenever  0 > 1 the model system (1) has only one equilibrium  * in the interior of ℧, so that condition (5) is verified.Now for the heroin model system (1) the task involves verifying the Bendixson criterion (65).Note that the variable  3 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: The Jacobian matrix of subsystem ( 45) is found to be In working with Theorem 12 one needs to make use of additive compound matrices.For an arbitrary 3 × 3 matrix , the second additive compound matrix  [2] is defined as ) . (47) Thus, the second additive compound matrix of Jacobian matrix  of system ( 45) is given as where For the model system (45) a suitable vector norm | ⋅ | in R 3 and a 3 × 3 matrix-valued function () are given by ) , ) . ( Note that upper prime (  ) denotes differentiation with respect to time, ./.Thus,  =    −1 +  [2]  −1 can be obtained as It is helpful to write matrix  in block form as where Following [13], let (, V, ) represent the vectors in and let  represent the Lozinskiȋ measure with respect to this norm.Applying the method of approximating () as given in [46] leads to where Thus,  1 and  2 are, respectively, given as Now from second and third equation of (45) it is easy to obtain the following: Substituting (60) into (58) and ( 61) into (59), respectively, leads to Now based on the definition of the method of approximating Lozinskiȋ measure of , () as given in [46], we now approximate the supremum of both  1 and  2 .Hence, (63) Thus we get the following inequality: Now the next step involves substituting () in and deducing whether  2 < 0. And if the inequality  2 < 0 does not hold we will need to establish a condition that leads to  2 < 0 being fulfilled.Considering uniform persistence, there exist  0 > 0 and  > 0 such that, for  > , the following is implied: Now by letting Γ 1 =  0 /(1 +  0 ) 2 and Γ 2 =  the following claim is made: if then it follows that where Now, for  > , it can be deduced that and, thus, the Bendixson criterion given by ( 42) is verified.However, it is important to observe that  2 < 0 if and only if condition (67) holds true.Thus, the following theorem is established.
Theorem 13.Provided that  0 > 1, if Γ 1 < Γ 2 , then system (1) has a unique endemic equilibrium  * which is globally asymptotically stable with respect to solutions of (1) originating in the interior of ℧.
The validity of Theorem 13 will be shortly verified numerically.

Numerical Examples
In this section numerical simulations of the heroin epidemic model are presented to support theoretical findings.Figure 4 which shows backward bifurcation is obtained by plotting heroin users equilibrium as a function of  0 .The figures The heroin eradication threshold (also referred to as critical reproduction number)   0 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  0 < 1.
Figures 5(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)  0 = 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(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  0 < 1.Indeed, parameter  has to exceed a certain fixed threshold   , hence supporting our theoretical findings that a nonzero equilibrium when  0 < 1 can only be maintained when  is greater than   ; see (28). Figure 5(d) shows the effect of treatment  on heroin users.High treatment leads to a steady decline of heroin users.
Figure 6 presents a scenario where  0 > 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  0 , nevertheless it does affect the model dynamics.The impact of  when  0 > 1 is different from the case where  0 < 1.For  0 < 1 it plays a key role in inducing bistability while for  0 > 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.

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  value < 0.05. Figure 8 displays the sensitivity analysis results for the heroin users not in treatment  1 ().It is straightforward to see that recruitment rate Λ, effective contact rate , relapsing rate of heroin users in treatment to heroin users not in treatment , 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 , , 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  suggesting that at initial stage of a heroin epidemic effective contact between heroin users  1 () 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 , and the extent of saturation of heroin users rate  parameters.

Concluding Remarks
In this study we formulated a heroin epidemic model with bilinear and saturated treatment function.The threshold parameter  0 usually referred to as the basic reproduction number plays a key role in the prediction of disease persistence or extinction.Epidemiologically, when  0 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  0 is decreased below one.In the present heroin epidemic model the analytical results indicate that  0 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  0 below another threshold referred to as the critical value of the basic reproduction number   0 .That is, heroin users can be eradicated if  0 <   0 < 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   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 Figures 7(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 , and the extent of saturation of heroin users parameter  are the ones which contribute to persistence of heroin epidemic within the community.such that the model equations ( 1) can be rewritten as / = (), where  = ( 1 ,  2 ,  3 ,  4 ).Hence, it follows that where  * 1 = Λ/ =  0 .With  =  * the transformed model equations (B.1) have a simple eigenvalue with zero real part and all other eigenvalues are negative (i.e., they have a hyperbolic equilibrium point).Thus, we can apply center manifold theory to investigate the local dynamics of the transformed system (B.1)near  =  * .Now we proceed to obtain the Jacobian matrix of the transformed system evaluated at heroin-free equilibrium HFE as

Figure 1 :
Figure1: 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).

Figure 2 :
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.

Figure 3 :
Figure 3: (a) and (b) represent bifurcations where drug users,  * 1 , not in treatment equilibrium are plotted as a function of  0 .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 <   = 0.1156,  = 0.01,  = 0.001,  1 = 0.002,  2 = 0.001,  = 0.015,  = 0.9,  = 0.467,  = 0.1, and Λ ∈ {1, 3} while (b) represents backward bifurcation and parameters used are the same as in (a) except  = 0.25 >   = 0.1156.(c) depicts the critical value   0 as a function of saturation parameter .The black dotted line represents the threshold   which if exceeded gives rise to backward bifurcation.(d) depicts the critical value   0 as a function of treatment rate .

Figure 4 :
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.,   0 ) corresponding to (a)-(d) are now 0.5204, 0.6038, 0.7314, and 0.9131, respectively.

Figure 6 :
Figure 6: (a) and (b) show the effect of increasing the parameter  that accounts for the extent of saturation of treatment when  0 > 1. Parameters used are the same as in Figure 5(a) except  = 0.46 and Λ = 3 corresponding to  0 = 1.4855 > 1.

Figure 7 :
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.

Table 2 :
Parameter baseline values and ranges used in sensitivity analysis.