Global Dynamics of a Discretized Heroin Epidemic Model with Time Delay

Abstract

We derive a discretized heroin epidemic model with delay by applying a nonstandard finite difference scheme. We obtain positivity of the solution and existence of the unique endemic equilibrium. We show that heroin-using free equilibrium is globally asymptotically stable when the basic reproduction number [subscript]R0[/subscript] <1 , and the heroin-using is permanent when the basic reproduction number [subscript]R0[/subscript] >1 .

(ProQuest: ... denotes non-US-ASCII text omitted.)

Academic Editor:Wanbiao Ma

College of Mathematics and System Sciences, Xinjiang University, Urumqi, Xinjiang 830046, China

Received 24 June 2014; Accepted 15 September 2014; 20 October 2014

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

1. Introduction

As we all know, the use of heroin and other drugs in Europe and more specifically in Ireland and the resulting prevalence are well documented [1-3]. It shows that the use of heroin is very popular and causes many preventable deaths. Heroin is so soluble in the fat cells that it crosses the blood-brain barrier within 15-20 seconds, rapidly achieving a high level syndrome in the brain and central nervous system which causes both the "rush" experience by users and the toxicity. Heroin-related deaths are associated with the use of alcohol or other drugs [4]. Treatment of heroin users is a huge burden on the health system of any country.

We often study infectious diseases with mathematical and statistical techniques; see, for example, [5-11]; however, little has been done to apply this method to the heroin epidemics. In 1979, Mackintosh and Stewart [9] considered an exponential model which is simplified from infectious disease model of Kermack and McKendrick to illustrate how the heroin-using spreads in epidemic fashion. They arranged a numerical simulation to show how the dynamics of spread are influenced by parameters in the model. White and Comiskey [5] attempted to extend dynamic disease modeling to the drug-using career and formulated an ordinary differential equation. They divided the whole population into three classes, namely, susceptible, heroin users, and heroin users undergoing treatment. Their model allows a steady state (constant) solution which represents an equilibrium between the number of susceptible, heroin users, and heroin users in treatment. Furthermore, this ODE model was revisited by Mulone and Straughan [12]; the authors proved that this equilibrium solution is stable both linearly and nonlinearly under the realistic condition in which relapse rate of those in treatment returning to untreated drug use is greater than the prevalence rate of susceptible becoming drug users. Recently, the study of the global properties and permanence of continuous heroin epidemic models attracted the researchers and have some very good results; see [13-16]. Specially, Samanta [15] considered a model with time-dependent coefficients and with different removal rates for three different classes, introduced some new threshold values _R* and ^R* , and obtained the permanence of heroin-using career.

Motivated by Samanta [15] and Zhang and Teng [8], we alter a nonautonomous heroin epidemic model with time delay to an autonomous heroin epidemic model. For convenience, we replace _U1 and _U2 by U and V , respectively. Thus, we obtain the following continuous heroin epidemic model with a distributed time delay: [figure omitted; refer to PDF] where S(t) , U(t) , and V(t) represent the number of susceptible, heroin users not in treatment, and heroin users in treatment, respectively. We assume that the time taken to become heroin user is s . The function η(s):[0,h][arrow right][0,∞) is nondecreasing and has bounded variation such that ^∫0h ...η(s)ds=η(h)-η(0)=1 .

For understanding more realistic phenomenon of heroin, a little complicated epidemic model is helpful. By applying Micken's nonstandard discretization method [17] to continuous heroin epidemics model with time delay (1), we derive the following discretized heroin epidemic model with a distributed time delay: [figure omitted; refer to PDF] where _Sn is the susceptible class, _Un is the class of heroin users not in treatment, and _Vn is the class of heroin users in treatment at n th step. Since the sufficient condition can be obtained, independently of the choice of a time step-size, we let the time step-size be one for the sake of simplicity. The nonnegative constants _μ1 , _μ2 , and _μ3 denote the death rate of the susceptible, heroin users not in treatment, and heroin users in treatment class, respectively. Throughout the paper, it is biologically natural to assume that _μ1 ...4;min...{_μ2 ,_μ3 } . The constant λ>0 denotes the recruitment rate of susceptible population from the general population. Constant P>0 is the proportion of heroin users who enter the treatment class. The individuals in treatment who stop using heroin are susceptible at a constant rate _ξ2 ...5;0 . Constant _β3 represents the transmission rate from heroin users in treatment to untreated heroin users. _β1 (_Un ) is the probability per unit time and the transmission is used with the form _β1 (_Un )_Sn+1^∑k=0h ..._Un-kηk , which includes various delays. By a natural biological meaning, we assume that _β1 (U) is a positive function and that there exists a constant _Uβ >0 such that _β1 (U) is nondecreasing on the interval [0,_Uβ ] . The integer h...5;0 is the time delay. The sequence _ηk :-∞<_ηk <+∞ (k=0,1,...,h) is nondecreasing and has bounded.

The initial conditions of the system (2) are given by [figure omitted; refer to PDF] where ^ψn(i) ...5;0 (n=-h,-h+1,...,0, i=1,2,3) . Again, by biological meaning, we further assume that ^ψ0(i) >0 for all i=1,2,3 .

The paper is organized as follows. In Section 2, we prove the positivity and boundedness of the solution of system (2). In Section 3, we deal with the global asymptotic stability of the heroin-using free equilibrium. In Section 4, we consider the permanence of the discrete epidemic model applying Wang's technique. In the discretized epidemic model, sufficient condition for global asymptotic stability and permanence are the same as for the original continuous epidemic model. We give some numerical examples and conclusion in Sections 5 and 6.

2. Basic Properties

For system (2), the heroin-using free equilibrium is given by [figure omitted; refer to PDF] Define a positive constant A...1;^∑k=0h ..._ηk . The stability of ^E0 is studied by using the next generation method in [7]. The associated matrix F (of the new heroin-using terms) and the M-matrix V (of the remaining transfer terms) are given as follows, respectively: [figure omitted; refer to PDF] Clearly, F is nonnegative, V is a nonsingular M-matrix, and V-F has Z sign pattern. The associated basic reproduction number, denoted by _R0 , is then given by _R0 =ρ(F^V-1 ) , where ρ is the spectral radius of F^V-1 . It follows that [figure omitted; refer to PDF]

Now, we will consider the positivity and boundedness of solution to system (2). For most continuous epidemic models, positivity of the solution is clear, but, for system (2), the positivity of the sequences _Sn ,_Un , and _Vn holds in some condition.

Lemma 1.

Let (_Sn ,_Un ,_Vn ) be any solution of system (2) with initial condition (3); then (_Sn ,_Un ,_Vn ) is positive for any n∈N and _V0 <P(1+_μ3 +_ξ2 )/_β3 .

Proof.

Let (_Sn ,_Un ,_Vn ) be any solution of system (2) with initial condition (3). It is evident that system (2) is equivalent to the following iteration system: [figure omitted; refer to PDF] In the following, we will use the induction to prove the positivity of solution. When n=0 , we have [figure omitted; refer to PDF] From (8)-(10), we see that, as long as _U1 is obtained, _V1 and _S1 will be obtained too.

If _U1 >0 , from (9), we directly obtain _V1 >0 and, from (10), we further obtain that _S1 >0 . Furthermore, we also have _N1 =_S1 +_U1 +_V1 >0 .

Let x=_U1 ; then, from (8)-(10), we see that x satisfies the following equation: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Substituting _V1 in _W2 (x) gets [figure omitted; refer to PDF] Since _V0 <P(1+_μ3 +_ξ2 )/_β3 , we have [figure omitted; refer to PDF] and because of _W1 (0)=λ+_S0 +_ξ2V0 /(1+_μ3 +_ξ2 )>0 , thus [figure omitted; refer to PDF] Substituting _W1 and _W2 in ψ(x) yields [figure omitted; refer to PDF] Here, constants a , b , c , m , and n are as follows: [figure omitted; refer to PDF] We take the limit on both sides of the above equation: [figure omitted; refer to PDF] this means that ψ(x)=0 has at least one positive solution x-∈(0,+∞) . So, we have _U1 =x->0 . Therefore, the positivity of _S1 >0 , _U1 >0 , and _V1 >0 is finally obtained. When n=2 , we have [figure omitted; refer to PDF] a similar argument as in the above for _U1 , _V1 , and _S1 ; we also can obtain that _U2 >0 , _V2 >0 , and _S2 >0 . Lastly, by using the induction, we can finally obtain that _Sn >0 , _Un >0 , and _Vn >0 , for all n>0 .

Now, we define the total population as _Nn =_Sn +_Un +_Vn . Then, from system (2), we know that [figure omitted; refer to PDF] Notice the assumption that _μ1 ...4;min...(_μ2 ,_μ3 ) ; we obtain [figure omitted; refer to PDF] If λ/_μ1 ...5;_N0 , it is easy to see that _Nn ...4;λ/_μ1 =^S0 , for all large n . If λ/_μ1 <_N0 , from right hand side of system (2), we obtain [figure omitted; refer to PDF] Hence, we have _N1 <_N0 and there exists i∈N such that _Ni ...4;λ/_μ1 =^S0 . Then, we may use this _Ni as a starting value instead of _N0 . This argument leads to the following result.

Lemma 2.

For any solution (_Sn ,_Un ,_Vn ) of system (2), the total population _Nn =_Sn +_Un +_Vn satisfies [figure omitted; refer to PDF] thus (_Sn ,_Un ,_Vn ) is ultimately bounded.

Let Ω={(_Sn ,_Un ,_Vn ):_Sn ,_Un ,_Vn ...5;0, _Sn +_Un +_Vn ...4;λ/_μ1 } ; then Ω is the positive invariant set to the solution of system (2).

In the following, we will examine the existence of endemic equilibrium for a special case of system (2).

Lemma 3.

Assume that _β1 (U)=_β1 >0 is a constant. If _R0 >1 , system (2) admits a heroin-using equilibrium ^E* =(^S* ,^U* ,^V* ) when V<P/_β3 , where ^E* satisfies following equality: [figure omitted; refer to PDF]

Proof.

Consider the following equation: [figure omitted; refer to PDF] From the first equation and the second equation of the system (25), we have [figure omitted; refer to PDF] From the second equation and the third equation of the system (25), we obtain [figure omitted; refer to PDF] thus, [figure omitted; refer to PDF] Since U...0;0 , from the second equation of the system, we have [figure omitted; refer to PDF] Substituting S in (28), we obtain [figure omitted; refer to PDF] Substituting U and S in (26) yields a quadratic equation of V as follows: [figure omitted; refer to PDF] where the coefficients are given by [figure omitted; refer to PDF] Since _R0 =λA_β1 /_μ1 (_μ2 +P+_ξ1 )>1 , then it is easy to see that c<0 and b>0 . According to Descartes' rule of signs, if a...5;0 , then F(V)=0 has a positive solution; if a<0 , then F(V)=0 has two positive solutions. From the expression of S and U , we note that V<P/_β3 . Since [figure omitted; refer to PDF] This means that F(V)=0 has a unique positive solution ^V* ∈(0,P/_β3 ) . Therefore, there exists a unique positive solution (^S* ,^U* ,^V* ) of system (2).

For the local stability of the equilibria, we refer to Theorem 2 in [7] and have the following results.

Theorem 4.

Assume that _β1 (U)=_β1 , _β1 is a positive constant. The heroin-using free equilibrium ^E0 =(^S0 ,0,0) of system (2) is locally asymptotically stable if _R0 <1 and unstable if _R0 >1 .

3. Global Asymptotic Stability of the Heroin-Using Free Equilibrium

In this section, we still assume that _β1 (U)=_β1 >0 , and obtain a sufficient condition for global asymptotic stability of the heroin-using free equilibrium ^E0 of system (2).

Using a Lyapunov function similar to that in [11], we can easily prove the global asymptotic stability of the heroin-using free equilibrium ^E0 .

Theorem 5.

If _R0 <1 , the drug-using free equilibrium ^E0 of system (2) is globally asymptotically stable.

Proof.

Let us take the following Lyapunov function: [figure omitted; refer to PDF] where _ci >0 (i=1,2,3) are the constants to be defined later and ^S0 =λ/_μ1 . Using system (2), the difference of _Hn satisfies [figure omitted; refer to PDF] From _Sn ...4;_Nn ...4;^S0 , for all n...5;0 , we have [figure omitted; refer to PDF] Let us choose _ci >0 (i=1,2,3) such that these constants satisfy the following inequalities: [figure omitted; refer to PDF] From (37), we have _c3β1^Sn+12 +(_c3β1^S0 -_β1 )_Sn+1 +_c2 >0 ; since _Sn+1 >0 , then the following inequality is true: [figure omitted; refer to PDF] that is, [figure omitted; refer to PDF] Since _R0 <1 , which implies that _β1 A^S0 <_μ2 +P+_ξ1 , we can choose _c2 =_β1^S0 +... ; here, ... (0<...<(_μ2 +_ξ1 -A_β1^S0 )/A) is a sufficiently small positive number such that _β1 A^S0 +A...<_μ2 +(1-_c1 )P+_ξ1 . Since _β1^S0 -2_c2 <0 and (_β1^S0 -2_c2⁾² >(_β1^S0)2 , we can choose _c3 >0 to satisfy (41). We may further choose _c1 >1 to satisfy (38). Therefore, ΔH is negative definite and is equal to zero if and only if _Sn+1 =^S0 , _Un+1 =0 , and _Vn+1 =0 . The proof is complete.

4. Permanence of System (2)

The system (2) is said to be permanent if there are positive constants m and M such that [figure omitted; refer to PDF] hold for any sequence _Sn of the system (2), and the same inequalities hold for _Un and _Vn . For each class _Sn ,_Un , and _Vn , m and M are independent of initial conditions.

Following the method used by Wang in [6], we will prove the permanence of system (2) for the general case; that is, assume that _β1 (U) is related to U .

Theorem 6.

If _R0 >1 , then system (2) is permanent for any initial condition (3).

Proof.

Firstly, from system (2) and Lemmas 1 and 2, for any _...0 >0 , there exists sufficiently large _n0 >0 such that _Un ...4;λ/_μ1 +_...0 as n...5;_n0 +h . Then, we have [figure omitted; refer to PDF] Let ^β1M (_...0 )=_{max...U∈[0,λ/μ1 +...0 ]β1} (U) . Thus, we have [figure omitted; refer to PDF] Notice that _...0 can be arbitrarily small. Then, we have [figure omitted; refer to PDF] Next, let us consider the positive sequences _Sn and _Un of (2). According to these sequences, we define [figure omitted; refer to PDF] Then, for n...5;0 , we obtain [figure omitted; refer to PDF] Since _R0 =_β1 (0)Aλ/_μ1 (_μ2 +P+_ξ1 )>1 , there exist 0<α<_Uβ and ρ>0 such that [figure omitted; refer to PDF] note that [figure omitted; refer to PDF]

We claim that it is impossible that _Un ...4;α holds for all n...5;_n1 ...5;[ρh] . The function [x] gives the smallest integer not less than x . Suppose the contrary, for n...5;_n1 +h . Consider [figure omitted; refer to PDF] From Lemma 1, _Sn satisfies [figure omitted; refer to PDF] and we have that, for n...5;_n1 +h+[ρh] , we have [figure omitted; refer to PDF] Hence, for n...5;_n1 +h+[ρh] , we have [figure omitted; refer to PDF] Let ...=_min...θ {_{Un1 +[ρh]+h+θ} ;θ=-h,-h+1,...,0} . Now, we will show that _Un ...5;... for all n...5;_n1 +[ρh]+h . In fact, there is an integer n->0 such that [figure omitted; refer to PDF] However, for n=_n1 +[ρh]+h+n- , we have [figure omitted; refer to PDF] Which is a contradiction. Thus, _Un ...5;... for n...5;_n1 +[ρh]+h . Therefore, for n...5;_n1 +[ρh]+h , [figure omitted; refer to PDF] which implies that _Hn [arrow right]+∞ as n[arrow right]+∞ . But, from Lemma 2 and (46), there exists a sufficiently large integer ^{n1[variant prime]} >0 such that, for n>^{n1[variant prime]} , [figure omitted; refer to PDF] which is a contradiction. Hence, the claim is proved.

In the rest, we only need to consider the following two cases:

(i) _Un ...5;α for all large n .

(ii) _Un oscillates about α for all large n .

We show that _Un ...5;_mu for all large n , where 0<_mu ...4;α , is a constant which will be given later. Clearly, we only need to consider case (ii). Let positive integers _n1 and _n2 be sufficiently large that _Un1 ...5;α , _Un2 ...5;α , and _Un <α , for _n1 <n<_n2 .

If _n2 -_n1 <h+[ρh] , since [figure omitted; refer to PDF] we have [figure omitted; refer to PDF] Hence, _Un >_mu for n∈[_n1 ,_n2 ] .

If _n2 -_n1 >h+[ρh] , we can easily obtain that _Un >_mu for n∈[_n1 ,_n1 +h+[ρh]] . Assume that there exists an integer n^>0 such that [figure omitted; refer to PDF] However, for n=_n1 +h+[ρh]+n^ , [figure omitted; refer to PDF] This is a contradiction to the proposition that _Un+1 <_mu . Therefore, _Un ...5;_mu for n∈[_n1 ,_n2 ] . Since these positive integers _n1 and _n2 are chosen in an arbitrary way, we conclude that _Un ...5;_mu for all large n in case (ii). Hence, _{liminf...n[arrow right]∞Un} ...5;_mu .

Note that, from that third equation of system (2), we have [figure omitted; refer to PDF] From Lemma 2 and the discussion above, we have [figure omitted; refer to PDF] The proof is completed.

5. Numerical Example

In order to confirm the validity of our results, we consider the following heroin epidemic model with a discrete time delay: [figure omitted; refer to PDF] Now, we present a numerical example. For the sake of simplicity, we choose the parameters as _β1 =0.9 , _β3 =0.8 , λ=2 , _μ1 =0.1 , _μ2 =0.2 , _μ3 =0.1 , P=0.4 , _ξ1 =0.1 , and _ξ2 =0.2 ; we get _R0 =25.7143<1 . Figure 1 shows that the disease free equilibrium ^E0 of the system (64) is globally asymptotically stable when _R0 <1 . Figure 2 shows that the system (64) is permanent when _R0 >1 .

Figure 1: [figure omitted; refer to PDF]

Figure 2: [figure omitted; refer to PDF]

6. Conclusions

In this paper, we have modified the Samanta heroin epidemic model into an autonomous heroin epidemic model with distributed time delay. Further, we established a discretized heroin epidemic model with time delay, sufficient conditions have been obtained to ensure the global asymptotic stability of heroin-using free equilibrium when _R0 ...4;1 and _β1 (U) is replaced by a positive constant. We also carried out some discussion about the heroin-using equilibrium, but our results are only restricted to the existence of this equilibrium for _β1 (U)=_β1 >0 , a special case of system (2). The stability of heroin-using equilibrium is yet to be studied. As a main result of this paper, we obtained the permanence of the system (2). From the expression of _R0 =_β1 (0)λA/_μ1 (_μ2 +P+_ξ1 ) , we see that a decrease in _β1 (transmission coefficient from susceptible population) will cause a decrease of the same proportion in _R0 . If the rate of migration or recruitment is restricted into susceptible community, the spread of the disease can also be kept under control by reducing λ and thereby decreasing _R0 . The spread of the heroin users can also be controlled by educators, epidemiologists, and treatment providers to increase the values of ξ (removal rate of heroin users not in treatment who stop using heroin but are susceptible) and P (proportion of heroin users who enter treatment) and thereby to decrease _R0 . This analysis tells us that prevention is better than cure; efforts to increase prevention are more effective in controlling the spread of habitual drug use than efforts to increase the numbers of individuals accessing treatment.

Acknowledgment

This work was supported by the National Natural Science Foundation of China (Grants nos. 11261056, 11261058, and 11271312).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] European Monitoring Centre for Dru gs and Drug Addiction (CMCDDA), "Annual Report, 2005," http://www.emcdda.europa.eu/publications/annual-report/2005

[2] C. Comiskey, "National prevalence of problematic opiate use in Ireland,", EM-CDDA, 1999.

[3] A. Kelly, M. Carvalho, C. Teljeur Prevalence of Opiate Use in Ireland 2002-2001. A 3-Source Capture Recapture Study. A Report to the National Advisory Committee on Drugs, Sub-committee o n Prevalence , Small Area Health Research Unit, Department of Public, 2003.

[4] K. A. Sporer, "Acute heroin overdose," Annals of Internal Medicine , vol. 130, no. 7, pp. 584-590, 1999.

[5] E. White, C. Comiskey, "Heroin epidemics, treatment and ODE modelling," Mathematical Biosciences , vol. 208, no. 1, pp. 312-324, 2007.

[6] W. Wang, "Global behavior of an SEIRS epidemic model with time delays," Applied Mathematics Letters , vol. 15, no. 4, pp. 423-428, 2002.

[7] P. van den Driessche, J. Watmough, "Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission," Mathematical Biosciences , vol. 180, pp. 29-48, 2002.

[8] T. Zhang, Z. Teng, "Global behavior and permanence of SIRS epidemic model with time delay," Nonlinear Analysis: Real World Applications , vol. 9, no. 4, pp. 1409-1424, 2008.

[9] D. R. Mackintosh, G. T. Stewart, "A mathematical model of a heroin epidemic: Implications for control policies," Journal of Epidemiology and Community Health , vol. 33, no. 4, pp. 299-304, 1979.

[10] H. Su, X. Ding, "Dynamics of a nonstandard finite-difference scheme for Mackey-Glass system," Journal of Mathematical Analysis and Applications , vol. 344, no. 2, pp. 932-941, 2008.

[11] M. Sekiguchi, E. Ishiwata, "Global dynamics of a discretized SIRS epidemic model with time delay," Journal of Mathematical Analysis and Applications , vol. 371, no. 1, pp. 195-202, 2010.

[12] G. Mulone, B. Straughan, "A note on heroin epidemics," Mathematical Biosciences , vol. 218, no. 2, pp. 138-141, 2009.

[13] J. Liu, T. Zhang, "Global behaviour of a heroin epidemic model with distributed delays," Applied Mathematics Letters , vol. 24, no. 10, pp. 1685-1692, 2011.

[14] G. Huang, A. Liu, "A note on global stability for a heroin epidemic model with distributed delay," Applied Mathematics Letters , vol. 26, no. 7, pp. 687-691, 2013.

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

[16] X. Wang, J. Yang, X. Li, "Dynamics of a Heroin epidemic model with very population," Applied Mathematics , vol. 2, no. 6, pp. 732-738, 2011.

[17] R. E. Mickens, "Discretizations of nonlinear differential equations using explicit nonstandard methods," Journal of Computational and Applied Mathematics , vol. 110, no. 1, pp. 181-185, 1999.