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.
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 by , 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. 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 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 , 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–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 . 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 (), individuals successfully treated from heroin use (), heroin users undergoing treatment (), and heroin users not in treatment (). 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 may begin to use heroin. Hence, the susceptible population is diminished due to contact with heroin users at rate , while heroin users increase at the same rate. Heroin users also increase when those undergoing treatment relapse at rate and return to their heroin using lifestyle. Heroin users reduce in number as a result of treatment which is represented by the treatment function Moreover, the user subpopulation is reduced by heroin-induced death at rate as well as a result of the self-decision to cease using heroin (also referred to as “self-cure”) at rate Individuals undergoing treatment are diminished through the following processes: relapse to heroin using at rate , heroin-induced death at rate , and successful treatment at rate Finally, the recovered/successfully treated subpopulation 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–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 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 the treatment function reduces to while for large it reduces to which actually characterizes the saturated phenomena of the treatment. Further, if , the treatment function becomes which is the usual linear treatment rate. 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 Table 1.
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 variables , , , are nonnegative for all time Hence, the following result.
Theorem 1. Let the initial conditions supplied to model (1) be such that , and Then the trajectories of the model (1), with positive initial conditions, will remain positive for all time
Proof. Let Now from the first equation of model (1), it follows that which can be written as Hence, so that Following a similar procedure we can show that , , and for all time 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) yieldsConsidering that , and letting , it follows from (7) that Therefore
Theorem 2. The closed set 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 , , , and are any solution of system (1) supplied with nonnegative initial conditions. Then it is straightforward to note that the total population satisfies the inequalityFrom (11) it follows that which implies if The standard comparison theorem [31] can be used to deduce that In particular if for all Thus, under the flow induced by system (1), the region is positively invariant. Furthermore, for 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 , a scenario representing a heroin-free state in the community. represents the number of susceptibles when no one is using heroin. The basic reproduction number denoted by 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 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 and the probability of surviving this compartment and moving to the treatment compartment is Now the probability of surviving heroin users in treatment class and then returning to the heroin users class not in treatment is Thus, the total average time spent by the heroin users not in the treatment compartment on multiple passes can be obtained asClearly, 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 asMultiplying (13) with the effective contact rate and the average recruitment rate we obtain heroin basic reproduction number asIt is easy to observe that is inversely proportional to treatment , which implies that if treatment rate is maintained sufficiently high it can control a heroin epidemic (by reducing 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 ; 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 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. Therefore hold. To obtain the endemic equilibria we set (1) to zero and solve for the equilibrium quantities , and in terms of That is,Substituting (15) into the second equation of (1) yieldswhere 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, , (16) reduces to a linear equation where so that model (1) has the unique solution which is nonnegative if and only if Hence, if , model (1) has a unique endemic equilibrium whenever and this equilibrium approaches zero as tends to one because But there are no positive endemic equilibria if These results are summarized in the following lemma.
Lemma 5. The epidemic model (1) when has a unique positive endemic equilibrium whenever and no positive endemic equilibrium otherwise.
In what follows, we investigate the global stability for both the HFE and the unique endemic equilibrium for the case
3.2. Global Stability for Heroin-Free Equilibrium When
To investigate global stability we apply the method presented by Castillo-Chavez et al. [34]. First let and with representing the number of individuals not using heroin and representing the number of individuals using heroin (i.e., heroin users in treatment and heroin users not in treatment). Now suppose where and denote differentiation with respect to time. The HFE is now denoted by , where The following conditions and have to be met to guarantee a local asymptotic stability: For , is globally asymptotically stable (g.a.s)., where , for .
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 point is g.a.s equilibrium of model (1) provided that (locally asymptotically stable) and that assumptions () and () hold.
Now consider the following theorem.
Theorem 7. Suppose Then the HFE is g.a.s.
Proof. Let and , and , where
Then we have It is straightforward to see that, at the heroin-free equilibrium (HFE) , Thus, Now, as , . Hence, is globally asymptotically stable (i.e., condition (H1) is satisfied).
Now consider so that Since the total population is bounded by , we have Thus, , which now implies that conditions and are satisfied. Consequently by Lemma 6 the fixed point is globally asymptotically stable when , which indicates nonexistence of multiple nontrivial equilibria when The epidemiological implication of HFE being g.a.s is that heroin epidemic will be eliminated from the community if the threshold quantity is decreased to (or maintained at) a value below unity.
Now for we establish the following theorem.
Theorem 8. For model (1) has (i)a unique positive endemic equilibrium if and either or ,(ii)a unique positive endemic equilibrium if (i.e., ) and ,(iii)two positive endemic equilibria if , and ,(iv)no positive endemic equilibrium if and either or
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.
For Case (i) where and (i.e., ) (16) becomes linear and has a unique nonzero solution which is positive if and negative if Referring to (15) we see that if is unique then so are , and
For Case (ii) where (i.e., ) and , (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) where there is a nonnegative endemic equilibrium at However, because (16) is quadratic and since the equilibrium is continuously determined by then there must be an interval to the left of on which two nonnegative equilibria coexist. That is,For Case (iv) where and or , (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 if and only if (i.e., ).
Proof. Consider (16), Note that, at , implies that the graph passes through the origin. If it follows that has a nonnegative root. Since is a continuous function of , if we increase such that , there is some open interval of say on which has two nonnegative roots. That is, there exist two nonnegative endemic equilibria when This is indeed true since Case (iv) of Theorem 8 has already shown that for model (1) does not have positive real roots when Note that, at , the following equality holds:This together with condition 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 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 should be reduced further below a certain threshold which we shall denote by (see Figure 3(b)).
(a)
(b)
(c)
(d)
3.3. Computation of New Threshold for Heroin Eradication
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 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 (i.e., ) and , It follows that if , (16) has one nonnegative root Supposing there is backward bifurcation at , then there are two endemic equilibria for an interval of values of basic reproduction number starting from a threshold defined by to a point where To obtain this threshold which is often referred to as the critical value of the basic reproduction number we replace the values of , , and into the equality to obtain a quadratic equation in terms of For mathematical tractability we redefine coefficients and as and , whereNote that remains as previously defined in (16). Now the quadratic equation in terms of can be obtained as Since we are considering the scenario where and , we have and thus just the single solution Thus, the critical value of basic reproduction number (i.e., the new threshold for heroin eradication), , is given as Consequently, from the above analysis of computation of threshold for heroin eradication we can deduce the following lemma.
Lemma 10. (a) If , then model (1) has a unique endemic equilibrium point . In this case heroin epidemic will persist in the community.
(b) If , then model (1) has two endemic equilibria and and signals that model (1) has backward bifurcation.
(c) If , then model (1) has only the heroin-free equilibrium point 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 , , , , , , , , , and . Figure 3(a) represents the forward bifurcation scenario where if the heroin-free equilibrium is globally asymptotically stable while when the heroin epidemic can persist. However, as we note from Figure 3(b), increasing parameter from to such that , a heroin epidemic can persist once established for a range of values that are below unity which indicates the occurrence of backward bifurcation. This implies that reducing below one will not necessarily be sufficient for eradication of heroin usage from the community. If is sufficiently decreased such that 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 there are a stable endemic equilibrium, an unstable endemic equilibrium, and a stable heroin-free equilibrium. When there is only one stable endemic equilibrium. Figure 3(c) shows the effect of the saturation parameter on ; namely, increasing decreases Figure 3(d) shows that increasing the treatment rate increases which epidemiologically implies that high cure rates of heroin users can lead to shrinking of the backward bifurcation regime.
4. Global Stability
According to Theorem 8, model (1) may have multiple equilibria when and a unique endemic equilibrium whenever First, global stability of the endemic equilibrium of model (1) is investigated for a special case, that is, when , using Lyapunov direct method, and later proven for the general model (i.e., ) using a geometric approach.
4.1. Global Stability of Endemic Equilibria Using Lyapunov Method (Special Case )
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 the unique endemic equilibrium of model (1) is globally asymptotically stable in the interior of if
Proof. Defining the following Lyapunov candidate function:Now computing the time derivative of along the solutions of system (1) results inBecause is an endemic steady point of model (1) when , then it follows thatUsing (35) in (34) yields Note that Replacing the above equality in (36) results in Hence, for all Hence, the heroin endemic equilibrium is stable and if and only if The largest compact invariant set when in is the singleton Therefore, by LaSalle’s invariance principle, the endemic steady state is globally asymptotically stable in the interior of
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 ordinary differential equations. This criterion is sometimes referred to as a geometric approach to global stability [14, 42].
To apply the geometric approach on model (1), consider the autonomous dynamical system , where and 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 map is a function for in an open subset and consider the following autonomous dynamical system:Let be the solution to (39) satisfying Now we make the following basic assumptions: is simply connected.There exists a compact absorbing set .Equation (39) has a unique equilibrium in
Now under the stated assumptions –, 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 a Bendixson criterion refers to a condition satisfied by field which precludes the existence of nonconstant periodic solutions of (39). When , (i.e., the planar case) the classical results (Poincaré-Bendixson theorem and Dulac criteria; see [43]) adequately provide such global conditions. For a remarkable approach for proving global stability is given by the work of Li and Muldowney [11, 13, 16]. They showed that if conditions – hold and differential equation (39) fulfils a Bendixson criterion that is robust under local -perturbations (a function is called a local -perturbation of at if there exists an open neighbourhood of in such that the support and , where ) of at all nonequilibrium nonwandering (a point is said to be nonwandering for system (39) if for any neighbourhood of in there exists arbitrary large such that 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 – Let be a matrix-valued function which is for and considerwhere is the directional derivative of in the direction of the vector field in system (39) and it is defined as
and 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 , is an matrix. Now define the following quantity:where is the Lozinskiĭ measure of with respect to vector norm in , , and is defined as (see [45, 46]). In paper [13] it is proved that if conditions and are satisfied then , 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 –, quantity implies the local stability of equilibrium point As a result the following theorem is true.
Theorem 12 (see [13]). Assuming that conditions – hold, then, the equilibrium point is globally asymptotically stable in if a function and a Lozinskiĭ measure exist such that quantity
Observe that whenever , 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 () holds) and (ii) in the interior of there exists an absorbing compact set (condition () holds). The heroin model studied here with the assumption that fulfils conditions ()-(). It is easy to prove that when , the heroin-free equilibrium is unstable (see Theorem 3). The instability of the heroin-free equilibrium combined with signals uniform persistence [47]. That is, there exists a positive constant such that for every solution of system (1) with in the interior of biologically feasible region satisfies 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 () is satisfied. Also it is shown that whenever the model system (1) has only one equilibrium in the interior of , so that condition () is verified. Now for the heroin model system (1) the task involves verifying the Bendixson criterion (65). Note that the variable 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 beIn working with Theorem 12 one needs to make use of additive compound matrices. For an arbitrary matrix , the second additive compound matrix is defined asThus, the second additive compound matrix of Jacobian matrix of system (45) is given aswhere For the model system (45) a suitable vector norm in and a matrix-valued function are given byNote that upper prime denotes differentiation with respect to time, Thus, can be obtained as It is helpful to write matrix in block form aswhereFollowing [13], let represent the vectors in Now for the norm in select and let represent the Lozinskiĭ measure with respect to this norm. Applying the method of approximating as given in [46] leads to whereHere and are operator norms of and with respect to the vector norm, where they are both regarded as mapping from to . represents the Lozinskiĭ measure of the matrix with respect to the norm in To obtain we sum the absolute value of the off-diagonal elements to the diagonal one in each column of and then take the maximum of two sums. Assuming that , it follows that Thus, and 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 and Hence, Thus we get the following inequality: Now the next step involves substituting inand deducing whether And if the inequality does not hold we will need to establish a condition that leads to being fulfilled.
Considering uniform persistence, there exist and such that, for , the following is implied:Now by letting and the following claim is made: ifthen it follows thatwhere Now, for , it can be deduced that and, thus, the Bendixson criterion given by (42) is verified. However, it is important to observe that if and only if condition (67) holds true. Thus, the following theorem is established.
Theorem 13. Provided that , if , 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.
5. 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 The figures present a scenario where is varied via parameter (i.e., ) 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) 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 .
(a)
(b)
(c)
(d)
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) 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.
(a)
(b)
(c)
(d)
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 . Indeed, parameter has to exceed a certain fixed threshold , hence supporting our theoretical findings that a nonzero equilibrium when 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 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 , nevertheless it does affect the model dynamics. The impact of when is different from the case where For it plays a key role in inducing bistability while for 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.
(a)
(b)
We now verify the global stability condition obtained using the geometric approach based on the following parameter values:(i) With these parameters the corresponding and 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 This leads to 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.
(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 , with a PRCC value close to 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
Figure 8 displays the sensitivity analysis results for the heroin users not in treatment 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 (), “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 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.
6. Concluding Remarks
In this study we formulated a heroin epidemic model with bilinear and saturated treatment function. The threshold parameter usually referred to as the basic reproduction number plays a key role in the prediction of disease persistence or extinction. Epidemiologically, when 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 is decreased below one. In the present heroin epidemic model the analytical results indicate that is indeed the threshold when the parameter However, when a saturated treatment function (i.e., ) 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 below another threshold referred to as the critical value of the basic reproduction number That is, heroin users can be eradicated if 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.
Appendix
A. Derivation of
Generally the geometric sequence is given as and the sum of a certain number of terms of the geometric sequence is given as , where is the sum of terms (th partial sum), is the first term, and is the common ratio. Now considering the geometric sequence from expression (12), we note that and Substituting and in (A.1) yieldsMultiplying (A.2) by gives the required expression which if multiplied by the effective contact rate and the average recruitment rate yields the basic reproduction number .
B. Proof of Existence of Backward Bifurcation
Proof. To prove existence of backward bifurcation in model equations (1), the center manifold approach as outlined by Castillo-Chavez and Song in [35] is used. First for clarity and understanding of the center manifold theory the model equations (1) variables are transformed as follows: , , , , and the total population Define ( denotes transpose), such that the model equations (1) can be rewritten as , where Hence, it follows thatNow let be the bifurcation parameter. Observe that, at , where . 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 It is easy to obtain the right eigenvectors of this Jacobian matrix as , where where Similarly, it is possible to obtain the left eigenvectors which we denote by as Now we proceed to obtain the bifurcation coefficients and as defined in Theorem 4.1 of [35].
Calculation of Coefficient . First the nonvanishing partial derivatives of the transformed model (B.1) evaluated at heroin-free equilibrium are obtained asso that
Calculation of Coefficient . The bifurcation coefficient is obtained as According to Theorem 4.1 of [35] if both bifurcation coefficients and are positive then model (1) will exhibit backward bifurcation. Observe that is always positive while if and only if Thus, if then model (1) will exhibit the phenomenon of backward bifurcation.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors are grateful to Royal Melbourne Institute of Technology (RMIT) University for supporting this research.