Global stability for a heroin model with two distributed delays

In this paper, we consider global stability for a heroin model with two distributed delays. The basic reproduction number of the heroin spread is obtained, which completely determines the stability of the equilibria. Using the direct Lyapunov method with Volterra type Lyapunov function, we show that the drug use-free equilibrium is globally asymptotically stable if the basic reproduction number is less than one, and the unique drug spread equilibrium is globally asymptotically stable if the basic reproduction number is greater than one.

1. Introduction.It is well-known that heroin is an opiate drug that is synthesized from morphine [1].It is more soluble in the fat cells so it crosses the blood-brain barrier within 15-20 seconds, rapidly achieving a high level syndrome in the brain and the central nervous system.Its rapid action causes both the 'rush' experience by users and the toxicity [2].Heroin users are at high risk for addiction.It is estimated that about 23 percent of individuals who use heroin become dependent on it.Over the past two decades, China has faced a dramatic increase in illicit drug abuse accompanying rapid economic reform and development [3].In addition to their deleterious somatic and psychological effects, heroin abuse and dependence 716 BIN FANG, XUEZHI LI, MAIA MARTCHEVA AND LIMING CAI constitute one of the most important modes of transmission of human immunodeficiency virus (HIV) and Hepatitis C virus (HCV) [4,5].This brings tremendous pressures and damages to social and public health.However, treatment of heroin users or users of other drugs such as cocaine not only is a costly procedure but also is a major burden on the health system of any country.
In fact, the spread of heroin habituation and addiction can be well modeled by epidemic-type models as "transmission" occurs in the form of peer pressure where established users recruit susceptible individuals into trying and using the drug [6,7,8].Mathematical modelling is a means to provide a general insight for how classes of drug takers behave, and as such, could hopefully become a useful device to aid specialist teams in devising treatment strategies.
Modeling heroin addiction and spread in epidemic fashion is not new [9].More recently, the spread of heroin addiction has already stimulated a number of articles using mathematical modeling ([10]- [13]).Heroin epidemic models typically divide the population into three classes, namely susceptibles, heroin drug users not in treatment, and heroin drug users undergoing treatment.The authors in [10,11] considered susceptibles, untreated users, and treated users model with standard incidence rate and showed that the steady states of the heroin epidemic model are stable.Wang et al in [12] considered the mass action incidence rate and proved that the drug use-free equilibrium and the unique endemic equilibrium are globally asymptotically stable by using the second compound matrix and under some conditions.Samanta [13] considered a non-autonomous heroin epidemic model, and got the global asymptotic stability of the system under some sufficient conditions using the method of Lyapunov function.
However, the models considered in these articles are ODE models and the effects of the time delay is not taken into account.In fact, as many infectious disease epidemics, the heroin epidemic also shows the effects of the delay.Recognizing the delay experienced by a treated heroin user who returns to use, [14,15] introduced distributed delay in the relapse term.[15] established the global stability for the heroin epidemic model by using a Lyapunov function.Reportedly, most users begin using heroin out of curiosity [16] which may peak and lead to use long after the first encounter with a user.Furthermore, heroin users rarely begin use with heroin, undergoing several stages of drug use before switching to heroin [17].To account for this delay in heroin use after the contact with a user, in this paper we incorporate time delay to describe the time needed for a susceptible individual to become an infectious heroin user.In particular, in this paper we present a heroin epidemic model with two distributed delays, based on the principles of mathematical epidemiology.First, we introduce the progression-to-use time delay to describe the time needed for a susceptible individual to become an infectious heroin user.For more realistic consideration, we assume that this delay is a distributed parameter over a finite interval.Second, we introduce the relapse delay to describe the time needed for a treated drug user to return to untreated drug user which varies according to the drug users' different temporal, social, and physical contexts.We also assume it to be a distributed parameter over a finite interval.We analyze the existence and stability of the equilibria of the model.It is shown that the existence, local and global asymptotical stability of the equilibria is completely determined by the basic reproduction number.Using a suitable Lyapunov function, we establish the global asymptotic stability of the heroin epidemic model with two distributed delays.
The paper is organized as follows.A heroin epidemic model with two distributed delays is formulated in Section 2. The existence of a unique drug spread equilibrium is also established in this section.The local asymptotic stability of the drug usefree equilibrium and the drug spread equilibrium is discussed in Section 3. The global asymptotic stability of the drug use-free equilibrium and the drug spread equilibrium is investigated in Section 4 by the use of a suitable Lyapnuov function.Finally in Section 5 we summarize our results.
2. The model.White and Comiskey [10] have produced an interesting model for the dynamics of heroin users.Mainly motivated by this work, in this paper we have introduced two distributed delays and considered the effects of distributed time delays in the following modified White-Comiskey mathematical model for the dynamics of heroin users.The model is formulated on the premise that drug use follows a process that can be modeled in a similar way to modeling of infectious disease spread ([18]- [22]): Here N (t) = S(t) + U 1 (t) + U 2 (t) denotes the total number of high-risk human population at time t.The meanings of all parameters in the above model are as following ( see Table 1 ).All parameters in model (2.1) are non-negative, Λ > 0, and µ > 0.
Usually, there are stages of involvement with drugs.Most users after exposure to a user begin using first alcohol, tobacco or inhalants.Then they progress through one or more less potent drugs before they try and start using stronger drugs like heroin [17].To account for this delay in the start of heroin use, we incorporate a distributed delay in the progression of a susceptible individual to a habitual use.To consider a scenario closer to reality, therefore, we incorporate the first time delay in our model (2.1) and we use it to describe the time that elapses for a susceptible individual to become an infectious heroin user.In reality, however, the delay period given by the stages of progression is not a fixed amount of time but varies from individual to individual, i.e., this delay is not the same for the whole population that uses heroin.Hence, we assume that this delay τ is a distributed parameter ( [22]- [24]) over the interval [0, h 1 ], where h 1 is the maximum of the delay and h 1 > 0. We assume that a susceptible individual at time t, S(t), contacts an individual who himself got into infectious contact t − τ units ago and is now infectious.As a result, the force of infection becomes βS(t) h1 0 f (τ )U 1 (t − τ )e −(µ+δ1+p)τ dτ , where the kernel function f (τ ) represents the distribution of the infectivity of heroin users to susceptible individuals where the time taken to become infectious heroin users is τ , that is, those infected at time t − τ become infectious τ (0 ≤ τ ≤ h 1 ) time units later with different probabilities.We assume that f (τ ) is non-negative and continuous, and satisfies h1 0 f (τ )dτ = 1.Furthermore, e −(µ+δ1+p)τ denotes the the removal rate that includes the drug-related deaths of users in treatment and a rate of successful "cure" that corresponds to recovery to a drug free life and immunity to drug addiction for the duration of the modeling time period k the probability of a drug user in treatment relapsing to not in treatment µ the natural death rate of the general population probability that a newly infected susceptible individual will survive the stages of progression to become an infectious heroin user.Similarly, drug users would return to untreated drug user class after cessation of a drug treatment programme after a period of time.Therefore, in model (2.1) the second time delay was incorporated and used to describe the time needed that a drug user returns to untreated drug user.This delay varies according to drug users' different characteristics and external influence.We also assumed it to be a distributed parameter over the interval [0, h 2 ], where h 2 is the maximum of the delay and h 2 > 0. We assume that g (τ ) is the distribution function of τ (0 ≤ τ ≤ h 2 ), which is non-negative and continuous, and satisfies h2 0 g (τ )dτ = 1.In the last term of the equations of model (2.1), we have At time t, the rate at which habitual heroin users have each acquired treatment at time η ∈ (t − h 2 , t) is pU 1 (η), and the probability that the individuals will survive from becoming treated at time η until relapsing to heroin users at time t is e −(µ+δ2)(t−η) .Therefore, when η changes from t − h 2 to t, the last term in the equations of model (2.1) represents all possible times at which heroin users in treatment might have relapsed to the users not in treatment.
The system is equipped with the following initial conditions: For biological reasons, we further assume that S 0 > 0, ϕ(0) > 0, U 0 2 > 0. To obtain continuity of the solutions to system (2.1), in this paper, we require By the third equation of system (2.1) and the initial conditions (2.3), we have We define the following space of functions The standard theory of functional differential equation [25], with initial conditions (2.2) that belong to the positive cone X and equation (2.4), implies that it can be verified that the system (2.1) has a unique solution (S(t), U 1 (t), U 2 (t)) which remains nonnegative for all t ≥ 0 [14].Moreover, we can show the solutions of system (2.1) are ultimately uniformly bounded in X.To see that fact, we add all equations of system (2.1) and we have It is easy to see that the set Ω is positively invariant for system (2.1).We note that U 2 (t) can be removed from the equations of system (2.1); it is sufficient to analyze the dynamical behavior of solutions to system (2.1) without the equation of U 2 (t).
For simplicity, we introduce the following notation: Furthermore, we impose the following assumptions: Assumption 2.1.Let h be the same definition in (2.2).We assume that: 1. f (τ ) and g(τ Obviously, we have 0 < a, b < 1. Next, we investigate the dynamics of the following system (2.5) Now we introduce the reproduction number of the heroin epidemic model, which is given by the following expression: To interpret formula (2.6) as a secondary number of heroin users produced by one heroin user, that is R 0 , we note that the average time in the drug users not in treatment class on the first pass is 1 µ+δ1+p and the probability of surviving this class is p µ+δ1+p .Since b is the probability of surviving the drug users in treatment class, the total average time in the drug users not in treatment class (on multiple passes) is Multiplying this by β Λ µ a gives R 0 , which is the average number of new drug users produced by a typical drug user not in treatment introduced into an entirely susceptible population [26].Thus, R 0 is the basic reproduction number, and acts as a threshold as is shown in the following result.System (2.5) always has the drug use-free equilibrium.The drug use-free equilibrium, in which the drug users are not present, is given by Now we show that there exists a drug spread equilibrium whose components are positive.For system (2.5), the drug spread equilibrium (S * , U * 1 ) needs to satisfy the equations: (2.8) Solving the last equation of (2.8), we get Substituting it into the first equation of (2.8), yields Therefore, there is a unique drug spread equilibrium E * (S * , U * 1 ) if R 0 > 1.So we have the following result.
In the next section, we investigate the local stability of the equilibria of system (2.5).

3.
Local stability of the equilibria.First, we investigate the local stability of the drug use-free equilibrium.
Theorem 3.1.The drug use-free equilibrium E 0 is locally asymptotically stable if R 0 < 1 and is unstable if R 0 > 1.
Proof.Taking the linearization of system (2.5) at the point E 0 , i.e., letting we get (3.1) To analyze the asymptotic behavior around E 0 , we look for solutions of the form X(t) = Xe λt , Y (t) = Ȳ e λt , where X and Ȳ are to be determined.Thus, we can consider the following eigenvalue problem: The characteristic equation of system (2.5) at E 0 is Hence, one characteristic root λ = −µ is a negative, and another characteristic root is determined by the following equation: Hence H(λ) has at least one positive root.Thus, if R 0 > 1, then E 0 is unstable.
Taking λ = x + yi (x, y ∈ R), and assuming that x ≥ 0 gives This is impossible, Thus, if R 0 < 1, then x < 0, and E 0 is locally asymptotically stable.Now we investigate the local stability of the drug spread equilibrium E * .We have the following result.
Proof.We linearize the system (2.5) about E * , by defining the perturbation variables: We have

.7)
To analyze the asymptotic behavior around E * , we look for solutions of the form s(t) = ŝe λt , u 1 (t) = û1 e λt , where ŝ and û1 are to be determined.Thus, we can consider the following eigenvalue problem: (3.12)where for the last equality we used the equilibrium equation (2.8).Consequently, we have |LHS| > |RHS| which is a contradiction.Hence, using Rouche's Theorem ( [27], Theorem 9.17.4),we show that equation (3.10) cannot have any roots with a nonnegative real part.Therefore, if R 0 > 1, the equilibrium E * is locally asymptotically stable.
4. Global stability of the equilibria.In the previous section we established the local stability of the equilibria, that is, given the required conditions on the parameters, if the initial conditions are close enough to the equilibrium, the solution will converge to that equilibrium.In this section our objective is to extend these results to global results.That is, given the conditions on the parameters, convergence to the equilibrium occurs independently of the initial conditions.Next, we address global asymptotic stability of system (2.5) by constructing appropriate Volterra type Lyapunov functions.
4.1.Global stability of the drug use-free equilibrium.As a first step, we establish the global stability of the drug use-free equilibrium.We will use Lyapunov function to approach the problem.Theorem 4.1.Assume R 0 < 1.Then the drug use-free equilibrium is globally asymptotically stable.
Proof.We will use a Lyapunov function.We adopt the Volterra type Lyapunov function used in [28]- [34].Define We note that H(x) ≥ 0 for all x > 0. H(x) achieves its global minimum at one, with H(1) = 0. Additionally, we also have Because of the complexity of the expressions, we define the Lyapunov function in components and take the derivative of each component separately.We consider the following Lyapunov function: where and We now show that d dt V (t) is non-positive.For clarity, we first find the derivatives of V 1 (t), V 2 (t) and V 3 (t), respectively, before combining.First, calculating the derivative of V 1 (t) in Equation (4.2) along the solutions of system (2.5), we obtain Here by using Λ = µS * 0 , we have (4.4) Next, calculating the time derivative of V 2 (t) along the solutions of system (2.5) and using integration by parts, we have Similarly, calculating the time derivative of V 3 (t) along the solutions of system (2.5) and using integration by parts, we have Combining the equations (4.4), (4.5) and (4.6) above, we get Therefore, if R 0 < 1, it follows that the positive-definite functional V (t) has nonpositive derivative dV dt .Let M be the largest invariant subset of {(S(t), U 1 (t))| dV dt = 0}.In order to have dV dt equal to zero it has necessary to have S(t) = S * 0 , U 1 (t) = 0. Hence, M = {E 0 }.By the LaSalle invariance principle [35,36], every solution of (2.5) tends to the drug use-free equilibrium E 0 , which is globally asymptotically stable.

4.2.
Global stability of the drug spread equilibrium.Now we are ready to establish the global stability of the drug spread equilibrium E * .To demonstrate that with a suitable Lyapunov functional W (t), we have to establish that W (t) ≤ 0 along the solution curves of system (2.5).The following Theorem summarizes the result.
Theorem 4.2.Assume R 0 > 1.Then, the drug spread equilibrium E * is globally asymptotically stable, that is, for any initial condition x 0 ∈ X, the solution of system (2.5) converges to E * .Proof.We use again a Lyapunov function.With H(x) = x − 1 − ln x (x ∈ R + ), we take the Lyapunov function as follows where and First, calculating the derivative of W 1 (t) in Equation (4.9) along the solutions of system (2.5), we obtain Here by using Λ = βS * U * 1 a + µS * , we have Second, calculating the time derivative of W 2 (t), we have And using integration by parts, we obtain (4.12) Similarly, calculating the time derivative of W 3 (t) and using integration by parts, we have Combining the equations (4.11), (4.12) and (4.13) above, we get (4.14) From the second equation of (2.8), we know that βS * a + pb − (µ + δ 1 + p) U 1 (t) = 0. (4.15) We also have where We want to show that the largest invariant set in Υ is the singleton {E * }.First, we notice that equality W (t) = 0 occurs if and only if S(t) = S * , U 1 (t − τ ) = U 1 (t) for all τ such that f, g > 0. Thus, at each point in Υ we have S(t) = S * , and therefore dS(t) dt = 0 in Υ.From the first equation of (2.5), we obtain This holds for all t, which implies that U 1 (t) = U * 1 for all t.Hence, we conclude that the largest invariant set in Υ is the singleton {E * }.By the LaSalle invariance principle [35,36], every solution of (2.5) tends to the endemic equilibrium E * , which is globally asymptotically stable.Therefore, we conclude that the drug spread equilibrium E * is globally stable.5. Discussion.Recently, several mathematical models (as mentioned in introduction) have been developed to describe the heroin epidemic.Most of these heroin epidemic models are ODE models and do not incorporate the effects of the delays.In this paper, we present a heroin epidemic model with two distributed delays, based on the principles of mathematical epidemiology.First, we introduce the progression-to-use delay into our model (2.1) to describe the time needed that a susceptible individual becomes an infectious heroin user.To be more realistic, we assume that this delay τ is a distributed parameter over the interval [0, h 1 ], where h 1 is the maximum of the delay.In addition, we introduced a relapse delay to describe the time needed by a drug user to return to drug use after treatment.This time also varies according to drug users' personal characteristics.We assumed it to be a distributed parameter over the interval [0, h 2 ], where h 2 is the maximum of the delay.We analyze the existence and stability of the equilibria of the model.We characterize the threshold conditions of the heroin epidemic model with an explicit formula for the reproduction number of heroin use.The reproduction number gives the number of secondary untreated users that one untreated user will cause in an entirely susceptible population.The reproduction number is the threshold which completely determines the stability of the equilibria.Using a suitable Lyapunov function, we show that the drug use-free equilibrium is globally stable if R 0 < 1.We also show that if R 0 > 1 the drug use-free equilibrium is unstable.In addition, there is a unique drug spread equilibrium which suggests that the heroin use persists in the population.For R 0 > 1 we show that the drug spread equilibrium is locally stable.Furthermore, using a suitable Lyapunov function, establish that the drug spread equilibrium is globally stable.
The reproduction number R 0 is an increasing function of transmission coefficient β which gives the rate of becoming a drug user, pb the probability of a drug user in treatment relapsing to untreated use, and a decreasing function of p, the rate of drug users who enter treatment.Our mathematical analysis suggest that the spread of the heroin use should be controlled through stringent screening measures to reduce the values of β, through educational campaigns at all social levels, and particularly to epidemiologists and treatment providers in order to increase the values of p. Furthermore, we have which signifies that as b increases, R 0 increases.Since b is the probability of leaving the treatment class and then entering the untreated class, then long time treatment is beneficial to control the spread of habitual drug use.For practical purposes, these results suggest that prevention is better than treatment.Efforts to increase prevention are more effective in controlling the spread of habitual heroin use than efforts to increase the number of individuals who have access to treatment.These results are provided with intention to inform and assist policy-makers in targeting prevention and treatment resources for maximum effectiveness.