# Stability Analysis of a Reaction-Diffusion Heroin Epidemic Model

Authors: Liang Zhang; Yifan Xing
Journal: Complexity
 (2020)

---

## Abstract

A reaction-diffusion (R-D) heroin epidemic model with relapse and permanent immunization is formulated. We use the basic reproduction numberR0 to determine the global dynamics of the models. For both the ordinary differential equation (ODE) model and the R-D model, it is shown that the drug-free equilibrium is globally asymptotically stable if R0≤1, and if R0>1, the drug-addiction equilibrium is globally asymptotically stable. Some numerical simulations are also carried out to illustrate our analytical results.


---

## Body

## 1. Introduction

Heroin comes from opioids, commonly known as “opium,” derived from the poppy plant. Pure heroin is a white powdery substance or a white crystalline powder. It is well known that long-term consumption and injection of heroin can cause personality disintegration, psychological metamorphosis, and lifespan to be reduced. There are more and more registered drug users, and this number is growing continuously [1, 2]. Heroin abuse and independence has been bringing tremendous pressures on the social and public health systems due to its prevalence all over the world [3, 4].Since the spread of heroin is as contagious as infectious diseases, it is a new trend to study heroin transmission from the perspective of infectious disease dynamics [5–16]. In 2007, White and Comiskey established an ODE model for heroin infectious diseases [14]. They studied the dynamics using a threshold R0 and showed that prevention is better than treatment. In 2009, this model was reconsidered by Mulone and Straughan [11], and the stability of the positive equilibrium point of the model is obtained by the authors by using the eigenvalue equation and Poincare–Bendixson theory. In 2011, Wang et al. [15] used the bilinear law incidence function instead of standard incidence, and they also analyzed the dynamic behavior of the heroin model. In order to better study the dynamics of heroin infectious diseases, many other different epidemic models have been formulated and studied in various ways [5–10, 12, 13, 16].As known to all that many people who are detoxifying are not really willing to go, and lots of people still keep in touch with friends who use drugs after they have a successful drug treatment, thus they have a high probability of reusing drugs. Hence, some scholars considered relapse in the drug model and analyzed its stability [12, 17]. However, there are some drug users who experienced the harmful effect of drugs and realized the happiness of being an ordinary person after they have a successful drug treatment, such persons will be far away from drugs, so we are optimistic that they will continue this good habit for a lifetime. In addition, there are some people who have been well educated since childhood and have a healthy living environment and strong willpower. So, they do not take drugs from the start to the finish. We call these two types of people permanently immunized against drugs.In recent years, it has been well recognized that spatial diffusion and environmental heterogeneity have important effects on the persistence and extinction of contagious diseases [18–38]. In the case of the heroin epidemic model, the spatial distribution of the susceptible or infected person is uneven, and the density may change at any time and place, thus it is more reasonable to use the R-D equations to describe the spread of drug abusers. Moreover, to the best of our knowledge, the heroin infectious disease models in which population density depends on both time and space variables are rarely studied.In this work, we first formulate an R-D heroin model with relapse and permanent immunization, and then study its global dynamics. We organize this paper as follows: in Section2, the model is derived and the positive property of the solutions for the model is proved. In Section 3, the model without diffusion is analyzed, and we obtain the basic reproduction number and show the stability of all the equilibria. In Section 4, the stability of the R-D model is obtained. In Section 5, we illustrate our results by some numerical simulations. Finally, we finish this paper with a concluding discussion.


## 2. The Model

### 2.1. Model Formulation

We divide the total population into five compartments:S,U1,U2,Q,andR. Here, S represents the number of susceptible individuals who have never used heroin; U1 represents the number of heroin users; U2 represents the number of heroin users undergoing treatment; Q represents the number of people who have used drugs and are not taking drugs at this stage, but may be taking drugs in the future; and R represents the number of people who never use drugs or these successful detoxification people do not take drugs anymore. We assume that drug users are not able to heal themselves through self-control, and if they want to abstain from drugs, they have to go through treatment. We also assume that not all people can be cured completely. If the drug users are successfully cured, individuals in compartment U2 will enter into the compartment of Q. If the treatment is terminated or failure, the people who failed the treatment will still take drugs. Some of these successful detoxification people will redrug because they cannot resist the temptation of drugs, and some will never take drugs because they know the harm of drug abuse. Thus, the total population is given by(1)N=S+U1+U2+Q+R.We give the transfer diagram of the model in Figure1. According to Figure 1, we obtain the following R-D heroin epidemic model with relapse and permanent immunization:(2)∂S∂t=d1ΔS+Λ−β1SU1−μ+α1S,l∈Ω,t>0,∂U1∂t=d2ΔU1+β1SU1+k2Q+β2U2−μ+δ1+kU1,l∈Ω,t>0,∂U2∂t=d3ΔU2+kU1−μ+δ2+k1+β2U2,l∈Ω,t>0,∂Q∂t=d4ΔQ+k1U2−μ+k2+α2Q,l∈Ω,t>0,∂R∂t=α1S+α2Q−μRl∈Ω,t>0,∂S∂n=∂U1∂n=∂U2∂n=∂Q∂n=∂R∂n=0,l∈Ω,t>0,Sl,0=S0l>0,U1l,0=U10l>0,U2l,0=U20l>0,Ql,0=Q0l>0,Rl,0=R0l>0,l∈Ω.Figure 1
Transfer diagram of the heroin epidemic model with relapse and permanent immunization.We assume that the bounded domainΩ⊂ℝ5 has a smooth boundary ∂Ω. As shown in (2), the homogeneous Neumann boundary conditions indicate that the population movements will not cross the border. Δ is the usual Laplacian operator on ℝ5. di>0i=1,2,3,4 are the diffusion coefficients. As mentioned above, R is a compartment of complete rehabilitation. The individuals in R are permanently immunized, and we can see that the R equation is uncoupled with the other equations of (2). Hence, the diffusion of R is not considered. We assume that all parameters in the model are positive constants, and the meaning of the parameters is described in Table 1.Table 1
Description of parameters.
ParameterDescriptionData estimatedData sourcesSl,tNumber of susceptible people at locationl and time tU1l,tNumber of heroin users at locationl and time tU2l,tNumber of heroin users undergoing treatment at locationl and time tQl,tNumber of people who have used drugs at locationl and time tRl,tNumber of people who never use drugs at locationl and time tΛRecruitment rate of the population1[39]μNatural death rate0.02[40]β1Addition rate fromS to abusersVariableβ2The proportion of failure treatment0.0011[40]δ1The heroin-related death rate ofU10.01Estimateδ2The heroin-related death rate of being treated0.005EstimatekProgression rate toU2 from U10.0095[40]k1The proportion of successful treatmentVariablek2Addition rate fromQ to abusersVariableα1The permanent withdrawal rate fromS to RVariableEstimateα2The permanent withdrawal rate fromQ to R0.0001EstimateTo investigate the global dynamic behavior of system (2), we first study its ODE counterpart version as follows:(3)S˙=Λ−β1SU1−μ+α1S,U˙1=β1SU1+k2Q+β2U2−μ+δ1+kU1,U˙2=kU1−μ+δ2+k1+β2U2,Q˙=k1U2−μ+k2+α2Q,R˙=α1S+α2Q−μR.


### 2.2. The Basic Properties of Model (3)

Lemma 1.
If the initial valuesS0,U10,U20,Q0,andR0 are positive, then model (3) has positive solutions St, U1t, U2t, Qt, and Rt for all t>0.Since the proof of the above lemma is direct, we omit it.Lemma 2.
All feasible solutions of the model (3) are bounded and enter the following region:(4)Ω=St,U1t,U2t,Qt,Rt∈R+50≤S+U1+U2+Q+R≤Λμ.Proof.
IfS,U1,U2,Q,R is a solution of model (3) with nonnegative initial conditions, adding the five equations yields(5)S˙+U˙1+U˙2+Q˙+R˙=Λ−μS−μU1−μU2−μQ−μR−δ1U1−δ2U2=Λ−μS+U1+U2+Q+R−δ1U1+δ2U2≤Λ−μS+U1+U2+Q+R=Λ−μNt,where(6)Nt=St+U1t+U2t+Qt+Rt,which indicates that(7)0≤Nt≤Λμ+N0e−μt,where N0 is the initial value. Thus, 0≤Nt≤Λ/μ, as t⟶∞. This completes the proof.


## 2.1. Model Formulation

We divide the total population into five compartments:S,U1,U2,Q,andR. Here, S represents the number of susceptible individuals who have never used heroin; U1 represents the number of heroin users; U2 represents the number of heroin users undergoing treatment; Q represents the number of people who have used drugs and are not taking drugs at this stage, but may be taking drugs in the future; and R represents the number of people who never use drugs or these successful detoxification people do not take drugs anymore. We assume that drug users are not able to heal themselves through self-control, and if they want to abstain from drugs, they have to go through treatment. We also assume that not all people can be cured completely. If the drug users are successfully cured, individuals in compartment U2 will enter into the compartment of Q. If the treatment is terminated or failure, the people who failed the treatment will still take drugs. Some of these successful detoxification people will redrug because they cannot resist the temptation of drugs, and some will never take drugs because they know the harm of drug abuse. Thus, the total population is given by(1)N=S+U1+U2+Q+R.We give the transfer diagram of the model in Figure1. According to Figure 1, we obtain the following R-D heroin epidemic model with relapse and permanent immunization:(2)∂S∂t=d1ΔS+Λ−β1SU1−μ+α1S,l∈Ω,t>0,∂U1∂t=d2ΔU1+β1SU1+k2Q+β2U2−μ+δ1+kU1,l∈Ω,t>0,∂U2∂t=d3ΔU2+kU1−μ+δ2+k1+β2U2,l∈Ω,t>0,∂Q∂t=d4ΔQ+k1U2−μ+k2+α2Q,l∈Ω,t>0,∂R∂t=α1S+α2Q−μRl∈Ω,t>0,∂S∂n=∂U1∂n=∂U2∂n=∂Q∂n=∂R∂n=0,l∈Ω,t>0,Sl,0=S0l>0,U1l,0=U10l>0,U2l,0=U20l>0,Ql,0=Q0l>0,Rl,0=R0l>0,l∈Ω.Figure 1
Transfer diagram of the heroin epidemic model with relapse and permanent immunization.We assume that the bounded domainΩ⊂ℝ5 has a smooth boundary ∂Ω. As shown in (2), the homogeneous Neumann boundary conditions indicate that the population movements will not cross the border. Δ is the usual Laplacian operator on ℝ5. di>0i=1,2,3,4 are the diffusion coefficients. As mentioned above, R is a compartment of complete rehabilitation. The individuals in R are permanently immunized, and we can see that the R equation is uncoupled with the other equations of (2). Hence, the diffusion of R is not considered. We assume that all parameters in the model are positive constants, and the meaning of the parameters is described in Table 1.Table 1
Description of parameters.
ParameterDescriptionData estimatedData sourcesSl,tNumber of susceptible people at locationl and time tU1l,tNumber of heroin users at locationl and time tU2l,tNumber of heroin users undergoing treatment at locationl and time tQl,tNumber of people who have used drugs at locationl and time tRl,tNumber of people who never use drugs at locationl and time tΛRecruitment rate of the population1[39]μNatural death rate0.02[40]β1Addition rate fromS to abusersVariableβ2The proportion of failure treatment0.0011[40]δ1The heroin-related death rate ofU10.01Estimateδ2The heroin-related death rate of being treated0.005EstimatekProgression rate toU2 from U10.0095[40]k1The proportion of successful treatmentVariablek2Addition rate fromQ to abusersVariableα1The permanent withdrawal rate fromS to RVariableEstimateα2The permanent withdrawal rate fromQ to R0.0001EstimateTo investigate the global dynamic behavior of system (2), we first study its ODE counterpart version as follows:(3)S˙=Λ−β1SU1−μ+α1S,U˙1=β1SU1+k2Q+β2U2−μ+δ1+kU1,U˙2=kU1−μ+δ2+k1+β2U2,Q˙=k1U2−μ+k2+α2Q,R˙=α1S+α2Q−μR.


## 2.2. The Basic Properties of Model (3)

Lemma 1.
If the initial valuesS0,U10,U20,Q0,andR0 are positive, then model (3) has positive solutions St, U1t, U2t, Qt, and Rt for all t>0.Since the proof of the above lemma is direct, we omit it.Lemma 2.
All feasible solutions of the model (3) are bounded and enter the following region:(4)Ω=St,U1t,U2t,Qt,Rt∈R+50≤S+U1+U2+Q+R≤Λμ.Proof.
IfS,U1,U2,Q,R is a solution of model (3) with nonnegative initial conditions, adding the five equations yields(5)S˙+U˙1+U˙2+Q˙+R˙=Λ−μS−μU1−μU2−μQ−μR−δ1U1−δ2U2=Λ−μS+U1+U2+Q+R−δ1U1+δ2U2≤Λ−μS+U1+U2+Q+R=Λ−μNt,where(6)Nt=St+U1t+U2t+Qt+Rt,which indicates that(7)0≤Nt≤Λμ+N0e−μt,where N0 is the initial value. Thus, 0≤Nt≤Λ/μ, as t⟶∞. This completes the proof.


## 3. Global Dynamics of the ODE Model (3)

### 3.1. The Basic Reproduction Number and Stability of Drug-Free Equilibrium

It is easy to get the drug-free equilibrium of system (3):(8)E0=Λμ+α1,0,0,0,α1Λμμ+α1.We now use the next-generation matrix method formulated in [41] to derive the basic reproduction number R0 of model (3).Letl=U1,U2,Q,S,RT, then system (3) can be written as(9)dldt=Pl−Wl,where(10)Pl=β1SU10000,Wl=μ+δ1+kU1−k2Q−β2U2μ+δ2+k1+β2U2−kU1μ+k2+α2Q−k1U2μ+α1S+β1SU1−ΛμR−α1S−α2Q.The Jacobian matrices ofPl and Wl at the drug-free equilibrium E0 are as follows:(11)DPE0=P3×300000000,DWE0=W3×300β1Λμ+α100μ+α1000−α2−α1μ,where(12)P3×3=β1Λμ+α100000000,W3×3=μ+δ1+k−β2−k2−kμ+δ2+k1+β200−k1μ+k2+α2.Then, the basic reproduction numberR0 is(13)R0=ρPW−1=β1Λμ+δ2+k1+β2μ+k2+α2μ+α1μ+δ1+kμ+δ2+k1+β2μ+k2+α2−kβ2μ+k2+α2−kk1k2.The following result on the local stability ofE0 can be obtained directly by Theorem 2 in [41], and we thus omit its proof.Theorem 1.
The drug-free equilibriumE0 is locally asymptotically stable for R0<1 and unstable for R0>1.Theorem 2.
IfR0≤1, the drug-free equilibrium E0 of system (3) is globally asymptotically stable.Proof.
Inspired by [28], we introduce the following Lyapunov function:(14)V=S−S0−S0lnSS0+U1+U2+Q.
It follows that the derivative ofV is(15)V˙=S˙+U˙1+U˙2+Q˙−S0SS˙=S0μ+α1−Sμ+α1−μ+δ2U2−μ+α2Q−S0SS0μ+α1−β1S0U1+μ+α1S0=−μ+α1S0S0S+SS0−2−μ+δ2U2−μ+α2Q−μ+δ11−β1Λμ+α1μ+δ1U1≤−μ+α1S0S0S+SS0−2−μ+δ2U2−μ+α2Q−μ+δ11−R0U1.
Here, we sued the equalitiesΛ=S0μ+α1. Since S0/S+S/S0−2≥0, if R0≤1, we have V˙≤0. Clearly, V˙≤0 if and only if S=S0andU1=U2=Q=0. Substituting S=S0andU1=U2=Q=0 into system (3), we get R⟶α1Λ/μμ+α1 as t⟶∞. According to LaSalle’s invariance principle [42], we obtain that the drug-free equilibrium E0 of system (3) is globally asymptotically stable in Ω if R0≤1. This completes the proof.


### 3.2. Existence and Stability of the Drug-Addiction Equilibrium

Let the right side of model (3) be equal to zero. Then we get(16)Λ−β1S∗U1∗−μ+α1S∗=0,β1S∗U1∗+k2Q∗+β2U2∗−μ+δ1+kU1∗=0,kU1∗−μ+δ2+k1+β2U2∗=0,k1U2∗−μ+k2+α2Q∗=0,α1S∗+α2Q∗−μR∗=0.Then, the unique drug-addiction equilibriumE∗=S∗,U1∗,U2∗,Q∗,R∗ exists if R0>1, where(17)S∗=Λβ1U1∗+μ+α1,U1∗=μ+α1R0−1β1,U2∗=kU1∗μ+δ2+k1+β2,Q∗=k1kU1∗μ+k2+α2μ+δ2+k1+β2,R∗=α1Λμβ1U1∗+μ+α1+α2k1kU1∗μμ+k2+α2μ+δ2+k1+β2.Theorem 3.
IfR0>1, then the unique drug-addiction equilibrium E∗ of system (3) is global asymptotically stable.Proof.
IfR0>1, then there exists a unique drug-addiction equilibrium E∗. We now introduce a Lyapunov function V as follows:(18)V=S−S∗−S∗lnSS∗+AU1−U1∗−U1∗lnU1U1∗+BU2−U2∗−U2∗lnU2U2∗+CQ−Q∗−Q∗lnQQ∗,where A, B, and C are the positive constants to be determined later. It follows that the derivative of V is(19)V˙=S˙1−S∗S+AU˙11−U1∗U1+BU˙21−U2∗U2+CQ˙1−Q∗Q=1−S∗SΛ−β1SU1−μ+α1S+A1−U1∗U1β1SU1+k2Q+β2U2−μ+δ1+kU1+B1−U2∗U2kU1−μ+δ2+k1+β2U2+C1−Q∗Qk1U2−μ+k2+α2Q=1−S∗Sβ1S∗U1∗+μ+α1S∗−β1SU1−μ+α1S+A1−U1∗U1β1SU1+k2Q+β2U2−β1S∗U1∗+k2Q∗+β2U2∗U1∗U1+B1−U2∗U2kU1−kU1∗U2∗U2+C1−Q∗Qk1U2−k1U2∗Q∗Q.
LetS/S∗=x, U1/U1∗=y, U2/U2∗=z, and Q/Q∗=u. Then, we have(20)V˙=1−1xβ1S∗U1∗+μ+α1S∗−β1S∗U1∗xy−μ+α1S∗x+A1−1yβ1S∗U1∗xy+k2Q∗u+β2U2∗z−β1S∗U1∗y−k2Q∗y−β2U2∗y+B1−1zkU1∗y−kzU1∗+C1−1uk1U2∗z−k1U2∗u=−μ+α1S∗1−x2x+β1S∗U1∗1−xy−1x+y+Aβ1S∗U1∗xy−y−x+1+Ak2Q∗u−y−uy+1+Aβ2U2∗z−y−zy+1+BkU1∗y−z−yz+1+Ck1U2∗z−u−zu+1=−μ+α1S∗1−x2x+β1S∗U1∗+Aβ1S∗U1∗+Ak2Q∗+Aβ2U2∗+BkU1∗+Ck1U2∗+xy−β1S∗U1∗+Aβ1S∗U1∗+yβ1S∗U1∗−Aβ1S∗U1∗−Ak2Q∗−Aβ2U2∗+BkU1∗+uAk2Q∗−Ck1U2∗+zAβ2U2∗−BkU1∗+Ck1U2∗−xAβ1S∗U1∗−uyAk2Q∗−zyAβ2U2∗−yzBkU1∗−zuCk1U2∗−1xβ1S∗U1∗.
The variables with nonnegative coefficients in (20) are xy, y, u, and z. If all the coefficients are positive, then V˙ is positive. If all the coefficients of xy, y, u, and z are equal to zero, then we get(21)Aβ1S∗U1∗−β1S∗U1∗=0,β1S∗U1∗−Aβ1S∗U1∗−Ak2Q∗−Aβ2U2∗+BkU1∗=0,Ak2Q∗−Ck1U2∗=0,Aβ2U2∗−BkU1∗+Ck1U2∗=0.
By (21), we obtain(22)A=1,B=β2U2∗+k2Q∗kU1∗,C=k2Q∗k1U2∗.
Hence, we have(23)V˙=−μ+α1S∗1−x2x+β1S∗U1∗2−x−1x+k2Q∗3−uy−yz−zu+β2U2∗2−zy−yz.
It is clear that−μ+α1S∗1−x2/x≤0 if x>0 and −μ+α1S∗1−x2/x=0 if only if x=1. By the relationships between the arithmetic mean and the geometric mean, we get 2−x−1/x≤0 if x>0 and 2−x−1/x=0 if and only if x=1; 3−u/y−y/z−z/u≤0 for x>0,y>0,u>0 and 3−u/y−y/z−z/u=0 if and only if y=z=u; 2−z/y−y/z≤0 for z>0,y>0 and 2−z/y−y/z=0 if and only if z=y. Therefore, V˙≤0 if x,y,z,u>0 and V˙=0 if and only if x=1andy=z=u. Substituting S=S∗andU1/U1∗=U2/U2∗=Q/Q∗ into the first equation of system (3), we get 0=Λ−β1S∗U1−μ+α1S∗. It follows from the first equation of (16) that U1=U1∗. Therefore, the maximum invariant set of system (2) on set x,y,z,u:V˙=0 is the singleton 1,1,1,1. This implies that the largest invariant set where V˙=0 is the singleton S∗,U1∗,U2∗,Q∗,R∗. Thus, by LaSalle’s invariance principle in [42], the drug-addiction equilibrium E∗ of model (3) is globally asymptotically stable when R0>1. This completes the proof.


## 3.1. The Basic Reproduction Number and Stability of Drug-Free Equilibrium

It is easy to get the drug-free equilibrium of system (3):(8)E0=Λμ+α1,0,0,0,α1Λμμ+α1.We now use the next-generation matrix method formulated in [41] to derive the basic reproduction number R0 of model (3).Letl=U1,U2,Q,S,RT, then system (3) can be written as(9)dldt=Pl−Wl,where(10)Pl=β1SU10000,Wl=μ+δ1+kU1−k2Q−β2U2μ+δ2+k1+β2U2−kU1μ+k2+α2Q−k1U2μ+α1S+β1SU1−ΛμR−α1S−α2Q.The Jacobian matrices ofPl and Wl at the drug-free equilibrium E0 are as follows:(11)DPE0=P3×300000000,DWE0=W3×300β1Λμ+α100μ+α1000−α2−α1μ,where(12)P3×3=β1Λμ+α100000000,W3×3=μ+δ1+k−β2−k2−kμ+δ2+k1+β200−k1μ+k2+α2.Then, the basic reproduction numberR0 is(13)R0=ρPW−1=β1Λμ+δ2+k1+β2μ+k2+α2μ+α1μ+δ1+kμ+δ2+k1+β2μ+k2+α2−kβ2μ+k2+α2−kk1k2.The following result on the local stability ofE0 can be obtained directly by Theorem 2 in [41], and we thus omit its proof.Theorem 1.
The drug-free equilibriumE0 is locally asymptotically stable for R0<1 and unstable for R0>1.Theorem 2.
IfR0≤1, the drug-free equilibrium E0 of system (3) is globally asymptotically stable.Proof.
Inspired by [28], we introduce the following Lyapunov function:(14)V=S−S0−S0lnSS0+U1+U2+Q.
It follows that the derivative ofV is(15)V˙=S˙+U˙1+U˙2+Q˙−S0SS˙=S0μ+α1−Sμ+α1−μ+δ2U2−μ+α2Q−S0SS0μ+α1−β1S0U1+μ+α1S0=−μ+α1S0S0S+SS0−2−μ+δ2U2−μ+α2Q−μ+δ11−β1Λμ+α1μ+δ1U1≤−μ+α1S0S0S+SS0−2−μ+δ2U2−μ+α2Q−μ+δ11−R0U1.
Here, we sued the equalitiesΛ=S0μ+α1. Since S0/S+S/S0−2≥0, if R0≤1, we have V˙≤0. Clearly, V˙≤0 if and only if S=S0andU1=U2=Q=0. Substituting S=S0andU1=U2=Q=0 into system (3), we get R⟶α1Λ/μμ+α1 as t⟶∞. According to LaSalle’s invariance principle [42], we obtain that the drug-free equilibrium E0 of system (3) is globally asymptotically stable in Ω if R0≤1. This completes the proof.


## 3.2. Existence and Stability of the Drug-Addiction Equilibrium

Let the right side of model (3) be equal to zero. Then we get(16)Λ−β1S∗U1∗−μ+α1S∗=0,β1S∗U1∗+k2Q∗+β2U2∗−μ+δ1+kU1∗=0,kU1∗−μ+δ2+k1+β2U2∗=0,k1U2∗−μ+k2+α2Q∗=0,α1S∗+α2Q∗−μR∗=0.Then, the unique drug-addiction equilibriumE∗=S∗,U1∗,U2∗,Q∗,R∗ exists if R0>1, where(17)S∗=Λβ1U1∗+μ+α1,U1∗=μ+α1R0−1β1,U2∗=kU1∗μ+δ2+k1+β2,Q∗=k1kU1∗μ+k2+α2μ+δ2+k1+β2,R∗=α1Λμβ1U1∗+μ+α1+α2k1kU1∗μμ+k2+α2μ+δ2+k1+β2.Theorem 3.
IfR0>1, then the unique drug-addiction equilibrium E∗ of system (3) is global asymptotically stable.Proof.
IfR0>1, then there exists a unique drug-addiction equilibrium E∗. We now introduce a Lyapunov function V as follows:(18)V=S−S∗−S∗lnSS∗+AU1−U1∗−U1∗lnU1U1∗+BU2−U2∗−U2∗lnU2U2∗+CQ−Q∗−Q∗lnQQ∗,where A, B, and C are the positive constants to be determined later. It follows that the derivative of V is(19)V˙=S˙1−S∗S+AU˙11−U1∗U1+BU˙21−U2∗U2+CQ˙1−Q∗Q=1−S∗SΛ−β1SU1−μ+α1S+A1−U1∗U1β1SU1+k2Q+β2U2−μ+δ1+kU1+B1−U2∗U2kU1−μ+δ2+k1+β2U2+C1−Q∗Qk1U2−μ+k2+α2Q=1−S∗Sβ1S∗U1∗+μ+α1S∗−β1SU1−μ+α1S+A1−U1∗U1β1SU1+k2Q+β2U2−β1S∗U1∗+k2Q∗+β2U2∗U1∗U1+B1−U2∗U2kU1−kU1∗U2∗U2+C1−Q∗Qk1U2−k1U2∗Q∗Q.
LetS/S∗=x, U1/U1∗=y, U2/U2∗=z, and Q/Q∗=u. Then, we have(20)V˙=1−1xβ1S∗U1∗+μ+α1S∗−β1S∗U1∗xy−μ+α1S∗x+A1−1yβ1S∗U1∗xy+k2Q∗u+β2U2∗z−β1S∗U1∗y−k2Q∗y−β2U2∗y+B1−1zkU1∗y−kzU1∗+C1−1uk1U2∗z−k1U2∗u=−μ+α1S∗1−x2x+β1S∗U1∗1−xy−1x+y+Aβ1S∗U1∗xy−y−x+1+Ak2Q∗u−y−uy+1+Aβ2U2∗z−y−zy+1+BkU1∗y−z−yz+1+Ck1U2∗z−u−zu+1=−μ+α1S∗1−x2x+β1S∗U1∗+Aβ1S∗U1∗+Ak2Q∗+Aβ2U2∗+BkU1∗+Ck1U2∗+xy−β1S∗U1∗+Aβ1S∗U1∗+yβ1S∗U1∗−Aβ1S∗U1∗−Ak2Q∗−Aβ2U2∗+BkU1∗+uAk2Q∗−Ck1U2∗+zAβ2U2∗−BkU1∗+Ck1U2∗−xAβ1S∗U1∗−uyAk2Q∗−zyAβ2U2∗−yzBkU1∗−zuCk1U2∗−1xβ1S∗U1∗.
The variables with nonnegative coefficients in (20) are xy, y, u, and z. If all the coefficients are positive, then V˙ is positive. If all the coefficients of xy, y, u, and z are equal to zero, then we get(21)Aβ1S∗U1∗−β1S∗U1∗=0,β1S∗U1∗−Aβ1S∗U1∗−Ak2Q∗−Aβ2U2∗+BkU1∗=0,Ak2Q∗−Ck1U2∗=0,Aβ2U2∗−BkU1∗+Ck1U2∗=0.
By (21), we obtain(22)A=1,B=β2U2∗+k2Q∗kU1∗,C=k2Q∗k1U2∗.
Hence, we have(23)V˙=−μ+α1S∗1−x2x+β1S∗U1∗2−x−1x+k2Q∗3−uy−yz−zu+β2U2∗2−zy−yz.
It is clear that−μ+α1S∗1−x2/x≤0 if x>0 and −μ+α1S∗1−x2/x=0 if only if x=1. By the relationships between the arithmetic mean and the geometric mean, we get 2−x−1/x≤0 if x>0 and 2−x−1/x=0 if and only if x=1; 3−u/y−y/z−z/u≤0 for x>0,y>0,u>0 and 3−u/y−y/z−z/u=0 if and only if y=z=u; 2−z/y−y/z≤0 for z>0,y>0 and 2−z/y−y/z=0 if and only if z=y. Therefore, V˙≤0 if x,y,z,u>0 and V˙=0 if and only if x=1andy=z=u. Substituting S=S∗andU1/U1∗=U2/U2∗=Q/Q∗ into the first equation of system (3), we get 0=Λ−β1S∗U1−μ+α1S∗. It follows from the first equation of (16) that U1=U1∗. Therefore, the maximum invariant set of system (2) on set x,y,z,u:V˙=0 is the singleton 1,1,1,1. This implies that the largest invariant set where V˙=0 is the singleton S∗,U1∗,U2∗,Q∗,R∗. Thus, by LaSalle’s invariance principle in [42], the drug-addiction equilibrium E∗ of model (3) is globally asymptotically stable when R0>1. This completes the proof.


## 4. Global Dynamics of the R-D Model (2)

### 4.1. Positivity and Boundedness of the Solutions

Theorem 4.
LetSl,t,U1l,t,U2l,t,Ql,t,Rl,t be a solution of system (2) and Sl,t,U1l,t,U2l,t,Ql,t,Rl,t∈CΩ¯×0,T∩C2,1Ω×0,T, where T is the maximal existing time. If Sl,0>0,U1l,0>0,U2l,0>0,Ql,0>0,andRl,0>0, then Sl,t>0,U1l,t>0,U2l,t>0,Ql,t>0,andRl,t>0, for all l,t∈Ω×0,T.Proof.
By the first equation of model (2), we get(24)∂S∂t≥d1ΔS−β1SU1−μ+α1S,that is,(25)∂S∂t−d1ΔS+β1U1+μ+α1S≥0.
SinceSl,0>0≡0, by Lemma 2.4.1 in [33], we get Sl,t>0.
Let(26)gU1=β1SU1+k2Q+β2U2−μ+δ1+kU1,gU2=kU1−μ+δ2+k1+β2U2,gQ=k1U2−μ+k2+α2Q.
Then, the second to the fourth equations of system (2) can be rewritten as(27)∂U1∂t−d2ΔU1−gU1l,t,U1,U2,Q−gU1l,t,0,U2,Q+gU1l,t,0,U2,Q−gU1l,t,0,0,Q+gU1l,t,0,0,Q−gU1l,t,0,0,0=∂U1∂t−d2ΔU1−β1S−μ−δ1−kU1+β2U2+k2Q=0,∂U2∂t−d3ΔU2−gU2l,t,U1,U2,Q−gU2l,t,0,U2,Q+gU2l,t,0,U2,Q−gU2l,t,0,0,Q+gU2l,t,0,0,Q−gU2l,t,0,0,0=∂U2∂t−d3ΔU2−kU1−μ+δ2+k1+β2U2=0,∂Q∂t−d4ΔQ−gQl,t,U1,U2,Q−gQl,t,0,U2,Q+gQl,t,0,U2,Q−gQl,t,0,0,Q+gQl,t,0,0,Q−gQl,t,0,0,0=∂Q∂t−d4ΔQ−k1U2−μ+k2+α2Q=0.
Define operatorℒi (i=1,2,3) on Ω×0,T as follows:(28)ℒ1U1=∂U1∂t−d2ΔU1,ℒ2U2=∂U2∂t−d3ΔU2,ℒ3Q=∂Q∂t−d4ΔQ.
We now consider the following parabolic system:(29)ℒ1U1+h11U1+h12U2+h13Q=0,ℒ2U2+h21U1+h22U2+h23Q=0,ℒ3Q+h31U1+h32U2+h33Q=0,U1l,0>0,U2l,0>0,Ql,0>0,∂U1∂t=∂U2∂t=∂Q∂t=0,where(30)h11=−β1S−μ−δ1−k,h12=−β2,h13=−k2,h21=−k,h22=−μ+δ2+k1+β2,h23=0,h31=0,h32=−k1,h33=μ+k2+α2.
Applying Theorem 4.2.4 in [33] to model (29), we get U1l,t≥0,U2l,t≥0,andQl,t≥0. So, by the second equation of model (2), we get(31)∂U1∂t≥d2ΔU1−μ+δ1+kU1.
That is,(32)∂U1∂t−d2ΔU1+μ+δ1+kU1≥0.
SinceU1l,0>0≡0, according to Lemma 2.4.1 in [33], U1l,t>0. Analogously,(33)∂U2∂t≥d3ΔU2−μ+δ2+k1+β2U2,∂Q∂t≥d4ΔQ−μ+k2+α2Q.
That is,(34)∂U2∂t−d3ΔU2+μ+δ2+k1+β2U2≥0,∂Q∂t−d4ΔQ+μ+k2+α2Q≥0.
Similarly, we obtainU2l,t>0andQl,t>0. If Rl,t>0 does not hold, then there exist l1∈Ω,t1∈0,T such that Rl1,t1=0,R˙l1,t1≤0. By the fifth equation of model (2), we get(35)R˙l1,t1=α1Sl1,t1+α2Ql1,t1.
Since∀l∈Ω,t∈0,T,Sl,t>0,Ql,t>0, R˙l1,t1=α1Sl1,t1+α2Ql1,t1>0. This is a contradiction. Therefore, Rl,t>0 for all l∈Ω,t∈0,T. This completes the proof.Theorem 5.
LetSl,t,U1l,t,U2l,t,Ql,t,Rl,t be a solution of model (2) and Sl,t,U1l,t,U2l,t,Ql,t,Rl,t∈CΩ¯×0,T∩C2,1Ω×0,T, where T is the maximum time of existence. If Sl,0>0,U1l,0>0,U2l,0>0,Ql,0>0,andRl,0>0, then the solution Sl,t,U1l,t,U2l,t,Ql,t,Rl,t is bounded on Ω¯.Proof.
By model (2), we get(36)∂S∂t+∂U1∂t+∂U2∂t+∂Q∂t−d1ΔS−d2ΔU1−d3ΔU2−d4ΔQ=Λ−μS+U1+U2+Q−α1S+δ1U1+δ2U2+α2Q.
It then follows that(37)∫Ω∂S∂t+∂U1∂t+∂U2∂t+∂Q∂t−d1ΔS−d2ΔU1−d3ΔU2−d4ΔQdl=∫ΩΛ−μS+U1+U2+Q−α1S+δ1U1+δ2U2+α2Qdl.
By Green’s formulas and Newman boundary∂S/∂t=∂U1/∂t=∂U2/∂t=∂Q/∂t=0,l∈∂Ω,t>0, we have(38)d1∫ΩΔSdl=d1∫∂Ω∂S∂nds=0,d2∫ΩΔU1dl=d2∫∂Ω∂U1∂nds=0,d3∫ΩΔU2dl=d3∫∂Ω∂U2∂nds=0,d4∫ΩΔQdl=d4∫∂Ω∂Q∂nds=0.
Hence,(39)∫Ω∂S∂t+∂U1∂t+∂U1∂t+∂Q∂tdl=∫ΩΛ−μS+U1+U2+Q−α1S+δ1U1+δ2U2+α2Qdl≤∫ΩΛ−μS+U1+U2+Qdl=ΛΩ−μ∫ΩS+U1+U2+Qdl.
Let∫ΩS+U1+U2+Qdl=Ft. Then, (39) becomes(40)dFtdt≤ΛΩ−μFt.
Which indicates that0≤Ft≤Λ/μΩ+F0e−μt, here(41)F0=∫ΩSl,0+U1l,0+U2l,0+Ql,0dl≤∫ΩSl,0+U1l,0+U2l,0+Ql,0∞dl=Sl,0+U1l,0+U2l,0+Ql,0∞Ω.
This shows thatFt=∫ΩS+U1+U2+Qdl is bounded. Let Z=Λ/μΩ+F0, then(42)Ft=∫ΩS+U1+U2+Qdl≤Z.
By Theorem 3.1 in [43], there is a positive constant Z∗ depending on Z so that(43)Sl,t+U1l,t+U2l,t+Ql,tL∞Ω≤Z∗.
Hence, we obtain thatSl,t,U1l,t,U2l,t,andQl,t are uniformly bounded on Ω¯.
For the last equation of model (2), let Sl,tandQl,t be bounded by S¯andQ¯. Then,(44)∂R∂t≤α1S¯+α2Q¯−μR.
Hence,Rt≤α1S¯+α2Q¯/μ+Rl,0∞ on Ω¯. This completes the proof.LetL≔CΩ¯;ℝ be a Banach space with a supremum norm:(45)w∞≔supwl,∀w∈CΩ¯;ℝ.DefineB:L5⟶L5, where L5=L×L×L×L×L, and let Bii=1,2,3,4 be a linear operator on L defined by(46)Biwl≔diΔwl,DBi≔w∈L:Δw∈L,∂w∂n=0on∂Ω.Then, by [23], we obtain that Bi are the infinitesimal generators of a strongly continuous semigroup etBit≥0 in L. For ∀φ≔φ1,φ2,φ3,φ4,φ5∈L5, there are(47)Bφl≔B1φ1l,B2φ2l,B3φ3l,B4φ4l,0T,DB≔DB1×DB2×DB3×DB4×L.They are also the infinitesimal generator of a strongly continuous semigroupetBt≥0 in Y≔L5, where(48)etB=etB1,etB2,etB3,etB4,1T,etBφ=etB1φ1,etB2φ2,etB3φ3,etB4φ4,φ5T,since Y is a Banach space with norm(49)φ1,φ2,φ3,φ4,φ5TY≔φ1L+φ2L+φ3L+φ4L+φ5L.LetG be a nonlinear operator on Y defined by(50)Gφ≔g1φ,g2φ,g3φ,g4φ,g5φ,where(51)g1φ=Λ−β1φ1φ2−μ+α1φ1,g2φ=β1φ1φ2+k2φ4+β2φ3−μ+δ1+kφ2,g3φ=kφ2−μ+δ2+k1+β2φ3,g4φ=k1φ3−μ+k2+α2φ4,g5φ=α1φ1+α2φ4−μφ5.Then, model (2) can be rewritten into a more abstract form in Y as follows:(52)dφtdt=Bφt+Gφt,φt≔S⋅,t,U1⋅,t,U2⋅,t,Q⋅,t,R⋅,tT,φ0≔S0⋅,U10⋅,U20⋅,Q0⋅,t,R0⋅,tT.By Proposition 4.16 in [44], it can be obtained that the unique continuously differentiable solution φ:0,Tmax⟶Y of the above equations has a maximum interval of existence 0,Tmax so that(53)φt=etBφ0+∫0tet−sBGφsds,and either Tmax⟶∞ or limsupt⟶Tmax−0φtY⟶+∞. By Theorem 5, we obtain Tmax=+∞ holds. Therefore, φt is a global solution.


### 4.2. Global Stability of the R-D Model

Since theR equation of model (2) is uncoupled with the other equations, model (2) can be reduced by ignoring R. It is clear that the R-D heroin epidemic model (2) has a drug-free equilibrium Er=Λ/μ+α1,0,0,0 and a unique positive drug-addiction equilibrium Er∗=S∗,U1∗,U2∗,Q∗ if R0>1.Theorem 6.
IfR0≤1, the drug-free equilibrium Er is globally asymptotically stable.Proof.
We give a Lyapunov function as follows:(54)W1t=∫ΩVφl,tdl,where V is given by (14) and φl,t=Sl,t,U1l,t,U2l,t,Ql,t. Direct computations show that(55)dW1dt=∫ΩgradφV⋅∂φ∂tdl=∫Ω1−S0S,1,1,1⋅S˙+d1ΔS,U˙1+d2ΔU1,U˙2+d3ΔU2,Q˙+d4ΔQdl=∫Ω1−S0SS˙+U1.+U2.+Q˙dl+∫Ω1−S0Sd1ΔSdl+∫Ωd2ΔU1dl+∫Ωd3ΔU2dl+∫Ωd4ΔQdl=∫ΩdVdtdl+∫Ω1−S0Sd1ΔSdl.
It follows by Green’s identity that(56)∫Ω1−S0Sd1ΔSdl=∫∂Ω1−S0Sd1∂S∂nds−∫Ωd1∇1−S0S∇Sdl=−∫Ωd1S0∇S2S2dl≤0.
Taking (56) into (55) and by Theorem 2, we obtain that if R0≤1, dV/dt≤0,∫ΩdV/dtdl≤0. Furthermore, if R0≤1,dW1/dt<0. By LaSalle’s invariance principle in [42], the drug-free equilibrium Er is globally asymptotically stable. This completes the proof.Theorem 7.
IfR0>1, then the drug-addiction equilibrium Er∗ is globally asymptotically stable.Proof.
We give the following Lyapunov function:(57)W2t=∫ΩV2φl,tdl,where φl,t=Sl,t,U1l,t,U2l,t,Ql,t and(58)V2=S−S∗−S∗lnSS∗+U1−U1∗−U1∗lnU1U1∗+β2U2∗+k2Q∗kU1∗U2−U2∗−U2∗lnU2U2∗+k2Q∗k1U2∗Q−Q∗−Q∗lnQQ∗.
Direct calculation yields(59)dW2dt=∫ΩgradφV2⋅∂φ∂tdl=∫Ω1−S∗S,1−U1∗U1,β2U2∗+k2Q∗kU1∗1−U2∗U2,k2Q∗k1U2∗1−Q∗Q⋅S˙+d1ΔS,U˙1+d2ΔU1,U˙2+d3ΔU2,Q˙+d4ΔQdl=∫Ω1−S∗SS˙+1−U1∗U1U˙1+β2U2∗+k2Q∗kU1∗1−U2∗U2U˙2+k2Q∗k1U2∗1−Q∗QQ˙dl+∫Ω1−S∗Sd1ΔSdl+∫Ω1−U1∗U1d2ΔU1dl+∫Ωβ2U2∗+k2Q∗kU1∗1−U2∗U2d3ΔU2dl+∫Ωk2Q∗k1U2∗1−Q∗Qd4ΔQdl=∫ΩdV2dtdl−d1S∗∫Ω∇S2S2dl−d2U1∗∫Ω∇U12U12dl−d3U2∗β2U2∗+k2Q∗kU1∗∫Ω∇U22U22dl−d4Q∗k2Q∗k1U2∗∫Ω∇Q2Q2dl.
By Theorem3, we know that if R0>1, then dV1/dt≤0,∫ΩdV2/dtdl≤0. Therefore, if R0>1, then dW2/dt<0. By LaSalle’s invariance principle in [42], the drug-addiction equilibrium Er∗ is globally asymptotically stable. This completes the proof.


## 4.1. Positivity and Boundedness of the Solutions

Theorem 4.
LetSl,t,U1l,t,U2l,t,Ql,t,Rl,t be a solution of system (2) and Sl,t,U1l,t,U2l,t,Ql,t,Rl,t∈CΩ¯×0,T∩C2,1Ω×0,T, where T is the maximal existing time. If Sl,0>0,U1l,0>0,U2l,0>0,Ql,0>0,andRl,0>0, then Sl,t>0,U1l,t>0,U2l,t>0,Ql,t>0,andRl,t>0, for all l,t∈Ω×0,T.Proof.
By the first equation of model (2), we get(24)∂S∂t≥d1ΔS−β1SU1−μ+α1S,that is,(25)∂S∂t−d1ΔS+β1U1+μ+α1S≥0.
SinceSl,0>0≡0, by Lemma 2.4.1 in [33], we get Sl,t>0.
Let(26)gU1=β1SU1+k2Q+β2U2−μ+δ1+kU1,gU2=kU1−μ+δ2+k1+β2U2,gQ=k1U2−μ+k2+α2Q.
Then, the second to the fourth equations of system (2) can be rewritten as(27)∂U1∂t−d2ΔU1−gU1l,t,U1,U2,Q−gU1l,t,0,U2,Q+gU1l,t,0,U2,Q−gU1l,t,0,0,Q+gU1l,t,0,0,Q−gU1l,t,0,0,0=∂U1∂t−d2ΔU1−β1S−μ−δ1−kU1+β2U2+k2Q=0,∂U2∂t−d3ΔU2−gU2l,t,U1,U2,Q−gU2l,t,0,U2,Q+gU2l,t,0,U2,Q−gU2l,t,0,0,Q+gU2l,t,0,0,Q−gU2l,t,0,0,0=∂U2∂t−d3ΔU2−kU1−μ+δ2+k1+β2U2=0,∂Q∂t−d4ΔQ−gQl,t,U1,U2,Q−gQl,t,0,U2,Q+gQl,t,0,U2,Q−gQl,t,0,0,Q+gQl,t,0,0,Q−gQl,t,0,0,0=∂Q∂t−d4ΔQ−k1U2−μ+k2+α2Q=0.
Define operatorℒi (i=1,2,3) on Ω×0,T as follows:(28)ℒ1U1=∂U1∂t−d2ΔU1,ℒ2U2=∂U2∂t−d3ΔU2,ℒ3Q=∂Q∂t−d4ΔQ.
We now consider the following parabolic system:(29)ℒ1U1+h11U1+h12U2+h13Q=0,ℒ2U2+h21U1+h22U2+h23Q=0,ℒ3Q+h31U1+h32U2+h33Q=0,U1l,0>0,U2l,0>0,Ql,0>0,∂U1∂t=∂U2∂t=∂Q∂t=0,where(30)h11=−β1S−μ−δ1−k,h12=−β2,h13=−k2,h21=−k,h22=−μ+δ2+k1+β2,h23=0,h31=0,h32=−k1,h33=μ+k2+α2.
Applying Theorem 4.2.4 in [33] to model (29), we get U1l,t≥0,U2l,t≥0,andQl,t≥0. So, by the second equation of model (2), we get(31)∂U1∂t≥d2ΔU1−μ+δ1+kU1.
That is,(32)∂U1∂t−d2ΔU1+μ+δ1+kU1≥0.
SinceU1l,0>0≡0, according to Lemma 2.4.1 in [33], U1l,t>0. Analogously,(33)∂U2∂t≥d3ΔU2−μ+δ2+k1+β2U2,∂Q∂t≥d4ΔQ−μ+k2+α2Q.
That is,(34)∂U2∂t−d3ΔU2+μ+δ2+k1+β2U2≥0,∂Q∂t−d4ΔQ+μ+k2+α2Q≥0.
Similarly, we obtainU2l,t>0andQl,t>0. If Rl,t>0 does not hold, then there exist l1∈Ω,t1∈0,T such that Rl1,t1=0,R˙l1,t1≤0. By the fifth equation of model (2), we get(35)R˙l1,t1=α1Sl1,t1+α2Ql1,t1.
Since∀l∈Ω,t∈0,T,Sl,t>0,Ql,t>0, R˙l1,t1=α1Sl1,t1+α2Ql1,t1>0. This is a contradiction. Therefore, Rl,t>0 for all l∈Ω,t∈0,T. This completes the proof.Theorem 5.
LetSl,t,U1l,t,U2l,t,Ql,t,Rl,t be a solution of model (2) and Sl,t,U1l,t,U2l,t,Ql,t,Rl,t∈CΩ¯×0,T∩C2,1Ω×0,T, where T is the maximum time of existence. If Sl,0>0,U1l,0>0,U2l,0>0,Ql,0>0,andRl,0>0, then the solution Sl,t,U1l,t,U2l,t,Ql,t,Rl,t is bounded on Ω¯.Proof.
By model (2), we get(36)∂S∂t+∂U1∂t+∂U2∂t+∂Q∂t−d1ΔS−d2ΔU1−d3ΔU2−d4ΔQ=Λ−μS+U1+U2+Q−α1S+δ1U1+δ2U2+α2Q.
It then follows that(37)∫Ω∂S∂t+∂U1∂t+∂U2∂t+∂Q∂t−d1ΔS−d2ΔU1−d3ΔU2−d4ΔQdl=∫ΩΛ−μS+U1+U2+Q−α1S+δ1U1+δ2U2+α2Qdl.
By Green’s formulas and Newman boundary∂S/∂t=∂U1/∂t=∂U2/∂t=∂Q/∂t=0,l∈∂Ω,t>0, we have(38)d1∫ΩΔSdl=d1∫∂Ω∂S∂nds=0,d2∫ΩΔU1dl=d2∫∂Ω∂U1∂nds=0,d3∫ΩΔU2dl=d3∫∂Ω∂U2∂nds=0,d4∫ΩΔQdl=d4∫∂Ω∂Q∂nds=0.
Hence,(39)∫Ω∂S∂t+∂U1∂t+∂U1∂t+∂Q∂tdl=∫ΩΛ−μS+U1+U2+Q−α1S+δ1U1+δ2U2+α2Qdl≤∫ΩΛ−μS+U1+U2+Qdl=ΛΩ−μ∫ΩS+U1+U2+Qdl.
Let∫ΩS+U1+U2+Qdl=Ft. Then, (39) becomes(40)dFtdt≤ΛΩ−μFt.
Which indicates that0≤Ft≤Λ/μΩ+F0e−μt, here(41)F0=∫ΩSl,0+U1l,0+U2l,0+Ql,0dl≤∫ΩSl,0+U1l,0+U2l,0+Ql,0∞dl=Sl,0+U1l,0+U2l,0+Ql,0∞Ω.
This shows thatFt=∫ΩS+U1+U2+Qdl is bounded. Let Z=Λ/μΩ+F0, then(42)Ft=∫ΩS+U1+U2+Qdl≤Z.
By Theorem 3.1 in [43], there is a positive constant Z∗ depending on Z so that(43)Sl,t+U1l,t+U2l,t+Ql,tL∞Ω≤Z∗.
Hence, we obtain thatSl,t,U1l,t,U2l,t,andQl,t are uniformly bounded on Ω¯.
For the last equation of model (2), let Sl,tandQl,t be bounded by S¯andQ¯. Then,(44)∂R∂t≤α1S¯+α2Q¯−μR.
Hence,Rt≤α1S¯+α2Q¯/μ+Rl,0∞ on Ω¯. This completes the proof.LetL≔CΩ¯;ℝ be a Banach space with a supremum norm:(45)w∞≔supwl,∀w∈CΩ¯;ℝ.DefineB:L5⟶L5, where L5=L×L×L×L×L, and let Bii=1,2,3,4 be a linear operator on L defined by(46)Biwl≔diΔwl,DBi≔w∈L:Δw∈L,∂w∂n=0on∂Ω.Then, by [23], we obtain that Bi are the infinitesimal generators of a strongly continuous semigroup etBit≥0 in L. For ∀φ≔φ1,φ2,φ3,φ4,φ5∈L5, there are(47)Bφl≔B1φ1l,B2φ2l,B3φ3l,B4φ4l,0T,DB≔DB1×DB2×DB3×DB4×L.They are also the infinitesimal generator of a strongly continuous semigroupetBt≥0 in Y≔L5, where(48)etB=etB1,etB2,etB3,etB4,1T,etBφ=etB1φ1,etB2φ2,etB3φ3,etB4φ4,φ5T,since Y is a Banach space with norm(49)φ1,φ2,φ3,φ4,φ5TY≔φ1L+φ2L+φ3L+φ4L+φ5L.LetG be a nonlinear operator on Y defined by(50)Gφ≔g1φ,g2φ,g3φ,g4φ,g5φ,where(51)g1φ=Λ−β1φ1φ2−μ+α1φ1,g2φ=β1φ1φ2+k2φ4+β2φ3−μ+δ1+kφ2,g3φ=kφ2−μ+δ2+k1+β2φ3,g4φ=k1φ3−μ+k2+α2φ4,g5φ=α1φ1+α2φ4−μφ5.Then, model (2) can be rewritten into a more abstract form in Y as follows:(52)dφtdt=Bφt+Gφt,φt≔S⋅,t,U1⋅,t,U2⋅,t,Q⋅,t,R⋅,tT,φ0≔S0⋅,U10⋅,U20⋅,Q0⋅,t,R0⋅,tT.By Proposition 4.16 in [44], it can be obtained that the unique continuously differentiable solution φ:0,Tmax⟶Y of the above equations has a maximum interval of existence 0,Tmax so that(53)φt=etBφ0+∫0tet−sBGφsds,and either Tmax⟶∞ or limsupt⟶Tmax−0φtY⟶+∞. By Theorem 5, we obtain Tmax=+∞ holds. Therefore, φt is a global solution.


## 4.2. Global Stability of the R-D Model

Since theR equation of model (2) is uncoupled with the other equations, model (2) can be reduced by ignoring R. It is clear that the R-D heroin epidemic model (2) has a drug-free equilibrium Er=Λ/μ+α1,0,0,0 and a unique positive drug-addiction equilibrium Er∗=S∗,U1∗,U2∗,Q∗ if R0>1.Theorem 6.
IfR0≤1, the drug-free equilibrium Er is globally asymptotically stable.Proof.
We give a Lyapunov function as follows:(54)W1t=∫ΩVφl,tdl,where V is given by (14) and φl,t=Sl,t,U1l,t,U2l,t,Ql,t. Direct computations show that(55)dW1dt=∫ΩgradφV⋅∂φ∂tdl=∫Ω1−S0S,1,1,1⋅S˙+d1ΔS,U˙1+d2ΔU1,U˙2+d3ΔU2,Q˙+d4ΔQdl=∫Ω1−S0SS˙+U1.+U2.+Q˙dl+∫Ω1−S0Sd1ΔSdl+∫Ωd2ΔU1dl+∫Ωd3ΔU2dl+∫Ωd4ΔQdl=∫ΩdVdtdl+∫Ω1−S0Sd1ΔSdl.
It follows by Green’s identity that(56)∫Ω1−S0Sd1ΔSdl=∫∂Ω1−S0Sd1∂S∂nds−∫Ωd1∇1−S0S∇Sdl=−∫Ωd1S0∇S2S2dl≤0.
Taking (56) into (55) and by Theorem 2, we obtain that if R0≤1, dV/dt≤0,∫ΩdV/dtdl≤0. Furthermore, if R0≤1,dW1/dt<0. By LaSalle’s invariance principle in [42], the drug-free equilibrium Er is globally asymptotically stable. This completes the proof.Theorem 7.
IfR0>1, then the drug-addiction equilibrium Er∗ is globally asymptotically stable.Proof.
We give the following Lyapunov function:(57)W2t=∫ΩV2φl,tdl,where φl,t=Sl,t,U1l,t,U2l,t,Ql,t and(58)V2=S−S∗−S∗lnSS∗+U1−U1∗−U1∗lnU1U1∗+β2U2∗+k2Q∗kU1∗U2−U2∗−U2∗lnU2U2∗+k2Q∗k1U2∗Q−Q∗−Q∗lnQQ∗.
Direct calculation yields(59)dW2dt=∫ΩgradφV2⋅∂φ∂tdl=∫Ω1−S∗S,1−U1∗U1,β2U2∗+k2Q∗kU1∗1−U2∗U2,k2Q∗k1U2∗1−Q∗Q⋅S˙+d1ΔS,U˙1+d2ΔU1,U˙2+d3ΔU2,Q˙+d4ΔQdl=∫Ω1−S∗SS˙+1−U1∗U1U˙1+β2U2∗+k2Q∗kU1∗1−U2∗U2U˙2+k2Q∗k1U2∗1−Q∗QQ˙dl+∫Ω1−S∗Sd1ΔSdl+∫Ω1−U1∗U1d2ΔU1dl+∫Ωβ2U2∗+k2Q∗kU1∗1−U2∗U2d3ΔU2dl+∫Ωk2Q∗k1U2∗1−Q∗Qd4ΔQdl=∫ΩdV2dtdl−d1S∗∫Ω∇S2S2dl−d2U1∗∫Ω∇U12U12dl−d3U2∗β2U2∗+k2Q∗kU1∗∫Ω∇U22U22dl−d4Q∗k2Q∗k1U2∗∫Ω∇Q2Q2dl.
By Theorem3, we know that if R0>1, then dV1/dt≤0,∫ΩdV2/dtdl≤0. Therefore, if R0>1, then dW2/dt<0. By LaSalle’s invariance principle in [42], the drug-addiction equilibrium Er∗ is globally asymptotically stable. This completes the proof.


## 5. Numerical Simulations

In this section, we shall carry some numerical simulations to illustrate our analytic results by using the parameter values in Table1. We fix Λ=1,μ=0.02,β2=0.0011,k=0.0095,α2=0.0001,δ1=0.01,andδ2=0.005.If we chooseβ1=0.0002,k2=0.00008,k1=0.008,andα1=0.02, then R0≈0.13<1 and E0=Λ/μ+α1,0,0,0,α1Λ/μμ+α1. Give different initial values I1=5,5,5,5,5andI2=5,10,3,6,1, respectively, and we can see that the drug-free equilibrium E0 is globally asymptotically stable (Figure 2). When α1=μ, we get Λ/μ+α1=α1Λ/μμ+α1, this is verified by the figure. We can see all solutions of the system converge to the drug-free equilibrium 25,0,0,0,25. If we keep β1,k2,andk1 unchanged and let α1=0.08>μ=0.02, then R0=0.05<1. According to the above discussion, E0 is asymptotically stable, and when α1>μ, Λ/μ+α1<α1Λ/μμ+α1. Figure 3 not only illustrates the stability of E0, but also shows the number of people who are permanently immunized against drugs, (Rt) is greater than the number of susceptible people (St) in the equilibrium E0 when α1>μ. By Figures 2 and 3, we can not only clearly see that U1t declined sharply and got to zero finally, but also see that all solutions of the system infinitely close to the drug-free equilibrium E0. It verifies the existence of E0.Figure 2
The stability ofE0 with different initial values when α1=μ. (a) The stability of E0 with an initial value I1. (b) The stability of E0 with an initial value I2.
(a)(b)Figure 3
The stability ofE0 with different initial values when α1≠μ. (a) The stability of E0 with an initial value I1. (b) The stability of E0 with an initial value I2.
(a)(b)If we choose parametersβ1=0.01,k2=0.008,k1=0.01,α1=0.02, and initial values I3=25,25,25,25,25 and I4=25,5,16,12,3, then R0≈6.5>1. Hence, the endemic equilibrium E∗ is global asymptotically stable (Figure 4). This also verifies the existence of E∗.Figure 4
The stability ofE∗ with different initial values. (a) The stability of E∗ with an initial value I3. (b) The stability of E∗ with an initial value I4.
(a)(b)


## 6. Discussion

In this paper, we formulated a novel R-D heroin epidemic model, which incorporates the relapse compartment and permanent immunization compartment. With the help of the next-generation matrix method, we obtained the basic reproduction numberR0. We obtained the global dynamics of the model by constructing some suitable Lyapunov functions. It is shown when R0≤1, the drug-free equilibrium is globally asymptotically stable; that is, the drug abuse will be eradicated; when R0>1, the endemic equilibrium is globally asymptotically stable, which means that drug abuse will be permanent.For the ODE system (3), k1,β2,k2,α1, and α2 represent the detoxification success rate, detoxification failure rate, relapse rate from Q to abusers, permanent withdrawal rates from S to R, and permanent withdrawal rates from Q to R, respectively. Figure 5 shows the relationship between R0 and k1, k2, β2, and α1, with other parameter values as given in Table 1. As shown in Figures 5(a)–5(c), R0 grows with k2 and β2 but decreases with k1, and it means that if we want to control heroin addiction, we should reduce k2 and β2 and increase k1, that is, increase the detoxification success rate and pay attention to people with a history of drug abuse to reduce their mental dependence on heroin. Moreover, as shown in Figure 5(d), R0 decreases with α2, and it means that if we want to control the heroin addiction, we should increase publicity to let people understand the harmful effect of heroin and take the initiative to stay away from drugs.Figure 5
The relationship betweenR0 and some parameter values. (a) The relationship between R0 and k1. (b) The relationship between R0 and k2. (c) The relationship between R0 and β2. (d) The relationship between R0 and α1.
(a)(b)(c)(d)The effect ofα1 and α2 on U1 is shown in Figure 6. It indicates that the values of α1 and α2 have a significant effect on the number of drug-addiction equilibrium. Let us observe Figure 6(a) first, the larger the α1 is, the fewer the people who use drugs in equilibrium is, and this tells us that we should not only pay attention to drug abuse but should pay more attention to those who do not use drugs. The government should strengthen publicity to raise people’s awareness of drug prevention. Comparing Figure 6(b) with Figure 6(a), we can get that although the effect of α2 on U1 is less than the effect of α1 on U1, it can still obtain that the larger the α2 is, the fewer the people who use drugs in equilibrium is; therefore, we also should pay attention to the people who have already quit drug abuse, so that they can stay away from drugs instead of reusing drugs.Figure 6
Effect of (a)α1 and (b) α2 on U1.
(a)(b)The diffusion phenomena make it more difficult for governments to control drugs. Fortunately, the R-D heroin epidemic model (2) still exists as a drug-free equilibrium, and it convinces us that the spread of drugs can be stopped. Additionally, the global dynamic behaviors of the R-D model (2) show that if R0>1 the endemic equilibrium Er∗ is global stability, and if R0≤1, the drug-free equilibrium Er is global stability, which indicates that the R-D model (2) may contain traveling waves connecting the steady states Er and Er∗. We shall conduct further research on this issue in the future.



---

*Source: 3781425-2020-06-08.xml*