# Analysis of a Diffusive Heroin Epidemic Model in a Heterogeneous Environment

Authors: Jinliang Wang; Hongquan Sun
Journal: Complexity
 (2020)

---

## Abstract

This paper is concerned with a reaction-diffusion heroin model in a bound domain. The objective of this paper is to explore the threshold dynamics based on threshold parameter and basic reproduction number (BRN)ℜ0, and it is proved that if ℜ0<1, heroin spread will be extinct, while if ℜ0>1, heroin spread is uniformly persistent and there exists a positive heroin-spread steady state. We also obtain that the explicit formula of ℜ0 and global attractiveness of constant positive steady state (PSS) when all parameters are positive constants. Our simulation results reveal that compared to the homogeneous setting, the spatial heterogeneity has essential impacts on increasing the risk of heroin spread.


---

## Body

## 1. Introduction

Nowadays, the prevalence of the use of heroin and other opiate drugs has become an increasing social problem. Every year, as reported in [1], approximately millions of drug users lose their ability to work and some of them even die due to drug addiction. It has been reported that in 2013, 2.7 per 100,000 persons of drug-poisoning deaths are related to heroin use in the USA which is a significant increase from 2010, when it was just only 0.7 per 100,000 person death.The prevalence of addiction to drugs should not be neglected and needs to attract attentions from the public health agency. From the view of mathematical models, heroin spread in a population can be regarded as a disease transmission in a population. In fact, we can use the way of the infection between susceptible individuals and infectious individuals to formulate the heroin spread mechanism. Suppose that a heroin drug user is introduced into a population, besides susceptible individuals, we divide the heroin drug user population into the following classes: without treatment class and with treatment class, respectively. Once susceptible individuals are contacted with heroin drug users, they are more or less to be influenced and become new heroin drug users. In a general setting, heroin drug users will receive treatment, and we label this subpopulation as heroin drug users with treatment. Due to the different social and physical contexts, heroin drug users with treatment may cease treatment and relapse to heroin user class, which brings in difficulties in controlling prevalence of addiction to drugs.Up to now, a large number of heroin epidemic models are often characterized by differential equation formulation. There have been quite a few publications along this line. To address heroin spread in ordinary differential equations (ODEs), White and Comiskey [2] used the standard incidence rate to model the interactions between susceptible individuals and heroin drug users. The stability of the positive equilibrium is analyzed, and the existence of a backward bifurcation is also addressed. Subsequently, the stability of equilibrium was further investigated by using the Poincaré–Bendixson theory in [3]. The model in [1] is in the form of ODEs with a nonlinear infection rate, and the main concern there is the investigation of the existence of various bifurcations.In order to describe the dynamics of heroin spread with delay differential equations (DDEs), time delay (which is interpreted as the time needed for relapse from the treatment class to the without treatment class) is introduced to ODE models (see, for example, [4]). The authors in [4] gave the stability analysis of the disease-free equilibrium by using BRN, but the stability analysis of the positive equilibrium is left as an open problem. Subsequently, Huang and Liu [5] solved the open problem that was left in [4] by using a suitable Lyapunov functional. By incorporating another time delay to describe the time needed for interactions between susceptibles and without treatment class, Fang et al. [6] carried out qualitative analysis of the model with two types of time delays.In the formulation of an age-structured model for heroin spread, besides timet, a continuous variable denoted by a is introduced to describe the treat age. Fang et al. in [7] investigated an age-structured model governed by a bilinear incidence function. In a recent work, Liu and Liu [8] formulated and investigated the infection age heroin model with susceptibility age and treat age. The threshold dynamics and control strategies in the light of BRN are achieved. We refer the readers to [9–11] for more literature on other types of heroin epidemic models.The aforementioned heroin epidemic models are only involved for a spatially homogeneous environment. However, spatial heterogeneity of the environment is ubiquitous due to the variance in environmental conditions such as temperature, humidity, and availability of resources. Mobility of host population and spatial heterogeneity are meaningful and important factors when considering a host population which lives in a spatially bounded domain (see, for example, [12–24]). This work intends to formulate and analyze a class of heroin epidemic model incorporating spatial heterogeneity and diffusion, and thus can be regarded as a continuation of the work [2, 5, 7, 11]. The topology heterogeneity is another important aspect in determining the spreading dynamics. Our work can also be extended to the information spreading dynamics on social networks (see, for example, [25–28]).In this paper, we useΩ⊂ℝn to denote a bounded domain, and the boundary of Ω is smooth and denoted by ∂Ω. Further, we suppose that the fluxes for each class are zero. At location x∈Ω and time t, the density of susceptible and heroin drug users without and with treatment classes is labeled by Ux,t, Vx,t, and Wx,t, respectively. We use the following initial conditions at location x∈Ω and time t=0:(1)Ux,0,Vx,0,Wx,0=U0x,V0x,W0x.To avoid overlapping writing, in the sequel, we always use the boundary condition at locationx∈∂Ω and time t>0:(2)dx∇U⋅n=dx∇V⋅n=dx∇W⋅n=0.Let us pay attention to the following diffusive heroin model at locationx∈Ω and time t>0:(3)Ut−∇⋅dx∇U=λx−μxU−βxUV,Vt−∇⋅dx∇V=βxUV−μx+pxV+kxW,Wt−∇⋅dx∇W=pxV−μx+kxW,where the symbol ∇ stands for the gradient operator; on the boundary of Ω, the label n stands for the outward normal vector; ∇⋅dx∇H stands for the divergence of dx∇H, where H=U,V,andW, respectively; we assume all model parameters are dependent of the location x allowing for spatial heterogeneity: dx is the diffusion rate, λx and μx are the recruitment rate and death rate, respectively, and the other parameters βx,px, and kx are the transmission rate, treatment rate, and relapse rate, respectively; all the above mentioned location-dependent parameters are subject to the following assumptions on Ω¯: positive, continuous, and uniformly bounded.We arrange the next section to give the preliminary results, including the uniform boundedness of solutions and the existence of a global solution. In Sections3 and 4, according to the theory of the next generation operator, we briefly determine BRN ℜ0 and further investigate its threshold role: if ℜ0<1, heroin spread will be extinct, while if ℜ0>1, heroin spread is uniformly persistent. Further, when ℜ0>1, the existence of the heroin-spread steady state is also confirmed. In Section 5, we also obtain that the explicit formula of ℜ0 and global attractiveness of constant heroin spread equilibrium with all model parameters are supposed to be positive constants (homogeneous case). The numerical examples are provided to support and validate theoretic results and explore the effect of the spatial heterogeneity in Section 6. This paper ends with some conclusions and discussion.


## 2. Well Posedness of Solutions of (3) with (2)

We place our problem on state spaceX≔CΩ¯,ℝ3, which equips with the supremum norm. The cone of X is denoted by X+≔CΩ¯,ℝ+3. Let ℤ≔CΩ¯,ℝ, and its cone is denoted by ℤ+≔CΩ¯,ℝ+. For convenience of notations, we denote θ1⋅,θ2⋅,θ3⋅=μ⋅,μ⋅+p⋅,μ⋅+k⋅. At location x∈Ω and time t>0, let Git,x,yi=1,2,3 denote the green function of the following equation Ht=∇⋅d⋅∇H−θi⋅H subject to the boundary condition (2). Hence, for i=1,2,3, we can define the solution semigroup Bit:ℤ⟶ℤ,t>0, as follows:(4)Bitψx=∫ΩGit,x,yψydy,ψ∈ℤ.According to [29] (Corollary 7.2.3), Bit,t>0 is strongly positive and compact.By the standard procedures, we further define the nonlinear operatorQ=Q1,Q2,Q3T:X+⟶X as(5)Q1ψ⋅=λ⋅−β⋅ψ1⋅ψ2⋅,Q2ψ⋅=β⋅ψ1⋅ψ2⋅+k⋅ψ3⋅,Q3ψ⋅=p⋅ψ2⋅,for any ψ∈X+. Under these settings, let v=U,V,WT, and we rewrite system (3) as(6)vt=ℬtψ+∫0tℬt−sQvsds,where ℬt=diagB1t,B2t,B3t.We then have the following preliminary theorem for (3) with (2).Theorem 1.
For anyψ∈X+, system (3) with (2) has a unique global solution v⋅,t,ψ on the interval 0,∞, with v⋅,0,ψ=ψ. Further, (3) generates a solution semiflow denoted by Φtψ≔v⋅,t,ψ, which possesses a global compact attractor.Proof.
Recall thatBit is defined in (4). Corresponding to which, ℬ is the linear part of system (3). As usual, we set the definition domain of ℬ as follows:(7)Dℬ=ϕ:dx∇ϕ⋅n=0on∂Ω,ℬϕ∈CΩ¯,ℝ3.
Hence,ℬ generates a C0-semigroup on X. Further, as standard procedures, we can check that(8)limh⟶0+1hdistψ+hQψ,X+=0,∀ψ∈X+.
Following the arguments stated in [30] (Corollary 4), we immediately obtain that on the interval0,τ∞, system (3) with a Neumann boundary condition has a local solution, where 0<τ∞≤∞. We next aim to give the proof on the uniform boundedness of the solution in Ω×0,τ∞, which in turn implies that the local solution is indeed a global solution. For this purpose, let ϒ=U+V+W, and we then add all equations of (3) to get(9)Πt−∇⋅dx∇Π=λx−μxΠ,x∈Ω,t>0,with boundary condition (2) and initial condition Πx,0=U0x+V0x+W0x≥0,x∈Ω. According to a standard argument as in [31], we conclude that (9) has a unique, positive, and global asymptotic stable steady state Π0, where Π0 in CΩ¯,ℝ satisfies(10)−∇⋅dx∇Π=λx−μxΠ,x∈Ω,dx∇Π⋅n=0,x∈∂Ω.
Hence, from the comparison principle, the boundedness of solutions of (3) is confirmed on the interval 0,τ∞. Moreover, the existence of a global solution of system (3) implies that system (3) generates a semiflow, denoted by Φt:X+⟶X+,t≥0, i.e., for ψ∈X,(11)Φtψ=v⋅,t;ψ=U⋅,t;ψ,V⋅,t;ψ,W⋅,t;ψT,∀x∈Ω¯,which is point dissipative. Thus, Φt has a global compact attractor which directly follows from the arguments as in [32] (Theorem 2.2.6).


## 3. Heroin-Free Dynamics

Heroin-free dynamics of (3) with (2) shall be explored in this section. From the classical procedures in [33, 34], let us first establish the basic reproduction number ℜ0. Letting V⋅,t=W⋅,t=0 in (3), from the first equation of (3), we know that U⋅,t is governed by (9). It follows that lim supt⟶∞U⋅,t−U0⋅=0. In such setting, the heroin-free steady state HFSS is denoted by U0⋅,0,0.Linearizing (3) around HFSS and further substituting the following solution V⋅,t,W⋅,t=eκtφ2⋅,eκtφ3⋅ into it yield the following linear system for x∈Ω,t>0:(12)κφ2−∇⋅d⋅∇φ2=β⋅U0φ2+k⋅φ3−μ⋅+p⋅φ2,κφ3−∇⋅d⋅∇φ3=p⋅φ2−μ⋅+k⋅φ3,subject to boundary condition (2) for x∈∂Ω,t>0. Since (12) is a cooperation system, we know that the linear eigenvalue problem problem (12) owns a principal eigenvalue κ0U0 with a positive eigenfunction φ20⋅,φ30⋅, which follows from the Krein–Rutman theorem.For convenience, we denoteY≔CΩ¯,ℝ2, and its cone is denoted by Y+≔CΩ¯,ℝ+2. Let us define Ψtψ≔A2tψ2,A3tψ3, ∀ψ2,ψ3∈Y, t≥0, where A2t and A3t are defined in (4). Hence, Ψt forms a positive C0-semigroup on Y. Precisely, for x∈Ω,t>0, Ψt is governed by(13)Vt−∇⋅d⋅∇V=−μ⋅+p⋅V,Wt−∇⋅d⋅∇W=−μ⋅+k⋅W,subject to boundary condition (2) for x∈∂Ω,t>0.Further, we define(14)F⋅=β⋅U0⋅k⋅p⋅0.Let us assume that both heroin drug users (without treatment and with treatment) are near the HFSS and introduce them by the distributionα=ψ2⋅,ψ3⋅ at time t=0. As time evolves, at time t, we can calculate the distribution of heroin drug users as(15)∑tα≔A2tψ2⋅,A3tψ3⋅.Hence, by the standard procedures as those in [34], the next-generation operator is regarded as the distribution of the total heroin drug users from the initial heroin drug users distribution α, which can be expressed as(16)ℒα⋅=F⋅∫0∞∑tαdt.This allows us to define the BRN by using the spectral radius of operatorℒ:(17)ℜ0=rℒ.The following result comes from [34], and we omit the proof here.Lemma 1.
(i)
ℜ0−1 has the same sign as κ0U0;(ii)
Ifℜ0>1, the HFSS becomes unstable;(iii)
Ifℜ0<1, the HFSS is locally stable.Denoting byλ∗≔minx∈Ω¯λx, μ∗≔maxx∈Ω¯μx, and β∗≔maxx∈Ω¯βx, we have the following result, which will be used later. Since the proof is similar to that in [35] (Lemma 2.3) with minor modification, we omit it here.Lemma 2.
(see [35], Lemma 2.3). For any ψ∈X+, we have the following statements on the solution of model (3):(i)
For anyt≥t0, if I⋅,t0,ψ,T⋅,t0,ψ≢0, we directly have U⋅,t,ψ,V⋅,t,ψ,W⋅,t,ψ>0,∀t≥t0,x∈Ω¯.(ii)
Further letδ1≔λ∗/μ∗+θβ∗ and θ=maxx∈Ω¯V⋅,t,ψ, then for x∈Ω¯, the following assertion holds:(18)liminft⟶∞U⋅,t,ψ≥δ1.
Now, let us consider heroin extinction with the thresholdℜ0.Theorem 2.
Letv⋅,t,ψ be the solution of system (3) with ψ∈X+, then we have the following statements:(i)
The unique HFSS,U0⋅,0,0, is globally asymptotically stable if ℜ0<1.(ii)
Ifℜ0>1, we have(19)limsupt⟶∞v⋅,t;ψ−U0⋅,0,0>ϵ0,for someϵ0>0.Proof.
We first prove (i). It follows fromℜ0<1 and Lemma 1 that κ0U0<0. From the continuity of κ0U0, we know that for ε>0, κ0U0+ϵ<0 holds. Corresponding to κ0U0+ϵ, there is a positive eigenvector φ˜ϵ⋅,ψ˜ϵ⋅. From the U-equation of (3), there is t0>0 such that for x∈Ω¯ and t>t0, U⋅,t<U0⋅+ϵ always hold. Hence, for x∈Ω,t>t0, V⋅,t and W⋅,t equations of (3) can be governed by(20)Vt−∇⋅d⋅∇V≤β⋅U0+ϵI+k⋅W−μ⋅+p⋅V,Wt−∇⋅d⋅∇W=p⋅V−μ⋅+k⋅W,with boundary condition (2). Choose a large enough number M such that(21)Mφ˜ϵ⋅,ψ˜ϵ⋅≥V⋅,t0,W⋅,t0.
Then, forx∈Ω and t>t0, V⋅,t,W⋅,t is a subsolution of the following linear auxiliary system:(22)V^t−∇⋅d⋅∇V^=β⋅U0+ϵV^+k⋅W^−μx+p⋅V^,W^t−∇⋅d⋅∇W^=p⋅V^−μ⋅+k⋅W^,V^⋅,t0=Mφ˜ϵ⋅,W^⋅,t0=Mψ˜ϵ⋅,with boundary condition (2). From the comparison principle (see, for example, [29]), we obtain that(23)V⋅,t,W⋅,t≤V^⋅,t,W^⋅,t=Meκ0U0+ϵt−t0φ˜ϵ⋅,ψ˜ϵ⋅.
Becauseκ0U0+ϵ<0, we directly have that as t⟶∞, V⋅,t,W⋅,t⟶0,0 uniformly for x∈Ω¯. By (9), we know that as t⟶∞, U⋅,t⟶U0⋅ uniformly for x∈Ω¯.
Next, we prove (ii). Fromℜ0>1, Lemma 1, and the continuity of κ0U0, we know that for ϵ0>0, κ0U0−ϵ0>0 holds. Corresponding to κ0U0−ϵ0, there is a positive eigenvector φ˜ϵ0x and ψ˜ϵ0⋅. Suppose, for the contrary, that there exists a positive solution of (3) such that(24)limsupt⟶∞v⋅,t;ϕ−U0,0,0<ϵ0.
It follows that there existst1>0 such that(25)U⋅,t>U0−ϵ0,0<V⋅,t<ϵ0,and0<W⋅,t<ϵ0,fort>t1.
Hence, forx∈Ω and t>t1, V⋅,t and W⋅,t equations of (3) can be governed by(26)Vt−∇⋅d⋅∇V≥β⋅U0−ϵ0V+k⋅W−μx+p⋅V,Wt−∇⋅d⋅∇W=p⋅V−μ⋅+k⋅W.
Choosing small enough as that with boundary condition (2). Without loss of generality, choosing η>0 small enough that(27)V⋅,t1,W⋅,t1≥ηφ˜ϵ0⋅,ψ˜ϵ0⋅.
Note that forx∈Ω and t>t1, ηeκ0U0−ϵ0t−t1φ˜ϵ0⋅,ψ˜ϵ0⋅ satisfies the following system:(28)V¯t−∇⋅d⋅∇V¯=β⋅U0−ϵ0V¯+k⋅W¯−μ⋅+p⋅V¯,W¯t−∇⋅d⋅∇W¯=p⋅V¯−μ⋅+k⋅W¯,V¯x,t1=ηφ˜ϵ0,W¯x,t1=ηψ˜ϵ0,with boundary condition (2). According to comparison principle (see, for example, [29]), we yield that for t≥t1,(29)V⋅,t,W⋅,t≥V¯⋅,t,W¯⋅,t=ηeκ0U0−ϵ0t−t1φ˜ϵ0⋅,ψ˜ϵ0⋅.
Consequently, we havelimt⟶∞V⋅,t⟶∞ and limt⟶∞W⋅,t⟶∞ as t⟶∞, which results in a contradiction. This proves (19).


## 4. Persistence of the Heroin Disease

Theorem 3.
For anyψ∈X+, we have the following statements:(i)
There exists a small enough numberδ>0 such that the solution of model (3) for ℜ0>1 is(30)liminft⟶∞w⋅,t,ψ≥δ,uniformly forx∈Ω¯,where w=U,V,andW, respectively.(ii)
System (3) has at least one PSS in Ω¯.Proof.
To proceed further, we divide the state spaceX+ into X0 and ∂X0 as follows:(31)X0≔ψ∈X+:ψ2⋅≠0andψ3⋅≠0,(32)∂X0≔X+X0=ψ∈X+:ψ2⋅≡0orψ3⋅≡0.
For these settings, we directly haveX=X0∪∂X0. From Lemma 2 (i), ΦtX0⊆X0 directly follows.
For the convenience of forthcoming discussions, we define the following set:(33)M∂=ψ∈∂X0:Φtψ∈∂X0,∀t≥0.
In what follows, we denote byωψ the omega limit set of ψ. We shall prove that ∪ψ∈M∂ωψ=HFSS. In fact, let ψ∈M∂ be given, and it follows that Φtψ∈∂X0, which indicates that V⋅,t,ψ≡0 or W⋅,t,ψ≡0.
In the case whereV⋅,t,ψ≡0 for all t≥0, we see from the U-equation of system (3), limt⟶∞U⋅,t=U0⋅ uniformly for x∈Ω¯. In view of the W-equation of system (3), we immediately yield that limt⟶∞W⋅,t,ψ=0 uniformly for x∈Ω¯.
In the case whereV⋅,t2,ψ≡0 for some t2≥0, Lemma 2 ensures that V⋅,t2,ψ>0 for all x∈Ω¯, ∀t≥t2. Thus, we have W⋅,t,ψ≡0, ∀t≥t2. In view of the W-equation of system (3), we immediately yield that limt⟶∞V⋅,t,ψ=0 uniformly for x∈Ω¯. By the U-equation of system (3), it then follows that limt⟶∞U⋅,t=U0⋅ uniformly for x∈Ω¯.
As a consequence,ωψ=HFSS, ∀ψ∈M∂.
Further from (ii) of Theorem2, it follows that if ℜ0>1, HFSS forms a uniform weak repeller for X0, i.e.,(34)limsupt⟶∞v⋅,t;ψ−U0x,0,0>ϵ0,∀ψ∈X+.
LetD⋅:X+⟶0,∞ be the continuous function defined by(35)Dψ=minminx∈Ω¯ψ2x,minx∈Ω¯ψ3x,∀ψ∈X+.
From the above definition, we directly have thatD−10,∞⊆X0, and D⋅ enjoys the following claim: if either Dψ=0 and ψ∈X0 or Dψ>0, then DΦtψ>0. Hence, D⋅ defines a distance function for Φt (see, for example, [36]).
Denote byWsHFSS the stable subset of HFSS. From the definition of M∂, HFSS is isolated in X+ in the sense that WsHFSS∩X0=∅. Moreover, there are no sets of HFSS from a cycle in ∂X0. Following the general results in [36], we can conclude that δ2>0 such that minDψ:ψ∈ωψ>δ2, that is,(36)liminft⟶∞C⋅,t,ψ≥δ2,C=V,W,∀ψ∈X0.
Combined with Lemma2 yields liminft⟶∞U⋅,t,ψ≥δ1, ∀ψ∈X0, where δ1 is a small enough number. Hence, we can choose δ=minδ1,δ2, and the uniform persistence result holds. According to [37] (Theorem 4.7), (3) with (2) has at least one PSS in X0. The positivity of the steady state directly follows from Lemma 2.


## 5. Special Case

In the case of all model parameters are positive constants, (3) admits a positive constant heroin-spread equilibrium. In the following, we aim to investigate the global behaviors of the heroin-spread equilibrium.Forx∈Ω,t>0, we shall study the following system:(37)Ut−dΔU=λ−βUV−μU,Vt−dΔV=βUV+kW−μ+pV,Wt−dΔW=pV−μ+kW,with the boundary condition(38)∂U∂n=∂V∂n=∂W∂n=0,x∈∂Ω,t>0,and initial condition(39)U⋅,0=U0⋅≥0,V⋅,0=V0⋅≥0,W⋅,0=W0⋅≥0.Clearly, (37) has a heroin-free equilibrium E0=λ/μ,0,0. The BRN of (37) takes the following form:(40)ℜ0=βλμ+kμμ+kμ+p−kμp,which is the same as the ODE system.We denote byE∗=U∗,V∗,W∗ the positive endemic equilibrium. Then, it satisfies(41)λ−βU∗V∗−μU∗=0,βU∗V∗+kW∗−μ+pV∗=0,pV∗−μ+kW∗=0.Direct calculation yields(42)U∗=μ2+kμ+pμβμ+k,V∗=μβℜ0−1,W∗=pμβμ+kℜ0−1.Clearly, ifℜ0>1, then (37) has a unique positive constant equilibrium E∗.In what follows, we give our main result in terms ofℜ0.Theorem 4.
Ifℜ0>1, then E∗ is globally attractive.Proof.
Define(43)LE∗t=∫ΩU∗gUU∗+V∗gVV∗+kW∗pV∗gWW∗dx,where gh=h−1−lnh, h>0.
Direct calculation of the derivative ofLE∗t yields(44)dLE∗tdt=∫Ω1−U∗UdΔU+1−V∗VdΔV+kW∗pV∗1−W∗WdΔWdx+∫Ω1−U∗UβU∗V∗+μU∗−βUV−μU+1−V∗VβUV+kW−βU∗V∗+kW∗VV∗+kW∗1−W∗WVV∗−WW∗dx.
Further, we have(45)dLE∗tdt=−d∫ΩU∗∇U2U2+V∗∇V2I2+kW∗2pV∗∇W2W2dx−∫ΩβU∗V∗+μU∗gU∗U+gUU∗+kW∗gWV∗W∗V+gW∗VWV∗dx≤0.
Here, we use the equality that2−A/B−B/A=−gA/B−gB/A, for A,B>0. It then follows from the above inequality that dLE∗t/dt=0 if and only if U,V,W=E∗. Hence, LE∗t defines a Lyapunov function. We can easily check that(46)U⋅,t,V⋅,t,W⋅,t⟶E∗inL2Ω3,ast⟶∞.
On the other hand,ϒ⋅,t=U⋅,t+V⋅,t+W⋅,t is bounded. It then follows from [38] (Theorem A2]) that for some positive constant C0, the following inequality holds:(47)U⋅,tC2Ω¯+V⋅,tC2Ω¯+W⋅,tC2Ω¯≤C0,∀t≥1.
Further, with the aid of the Sobolev embedding theorem, we get the conclusion that(48)U⋅,t,V⋅,t,W⋅,t⟶E∗inL∞Ω3,ast⟶∞,namely, E∗ attracts all solutions of (37). Thus, the global attractiveness of E∗ directly follows. □


## 6. Numerical Examples

This section is devoted to provide numerical examples to support our results. Further, we shall study the effect of the spatial heterogeneous environment. Throughout this section, we set the domain asΩ=−2,2⊂ℝ. We use labels S,I,T instead of U,V,W in the previous section to stand for the density of susceptible and heroin drug users without and with treatment classes.Example 1.
Threshold dynamics.
We will check the threshold dynamics obtained in Section5. To this end, we set parameters in model (37) as follows:(49)λ=10,μ=0.01,β=0.0001,d=0.0125,k=0.01,p=0.1.
With this setting, we can obtain thatℜ0=1.666667. From Theorem 4, the solutions of model (37) approach to E∗, that is, the heroin spread will be established, see Figure 1.
If we letp=0.2 in (49), we can obtain ℜ0=0.909091. From Theorem 2 (i), 1000,0,0 is globally asymptotically stable, that is, the heroin epidemic will be extinct, see Figure 2.Figure 1
(a)–(d) Dynamical behaviors of the solutions withℜ0=1.666667. λ=10,μ=0.01,β=0.0001,d=0.0125,k=0.01,andp=0.1. The initial value is S0,I0,T0=60,2,0.
(a)(b)(c)(d)Figure 2
(a)–(d) Dynamical behaviors of the solutions withℜ0=0.909091. λ=10,μ=0.01,β=0.0001,d=0.0125,k=0.01,andp=0.2. S0,I0,T0=60,2,0.
(a)(b)(c)(d)Example 2.
Influence of the diffusion ratedx.
We set(50)dx=1.25×10−2,1.25×10−3,and1.25×10−4,respectively. Under these settings, we will investigate the situations of spatial distribution of heroin drug user Ix,t prevalence with different diffusion rates. Figure 3 demonstrates the distribution of heroin drug users Ix,t. We find that the larger diffusion rate plays an increasingly important role in the higher prevalence.Figure 3
Ix,t prevalence with dx=1.25×10−2,1.25×10−3,and1.25×10−4, respectively. (a) ℜ0>1; (b) ℜ0<1.
(a)(b)Example 3.
Influence of the transmission rateβx.
We set(51)λ=5.5,μ=0.01,d=0.0125,k=0.01,p=0.1,(52)βx=0.0001⋅1+cβ⋅cosπx.
Here,cβ∈0,1 measures the magnitude of spatial heterogeneity. Obviously, cβ=0 means the spatial homogeneous case, and cβ∈0,1 means the higher spatial heterogeneity.
In Figure4(a), we find that ℜ0 increases in terms of cβ in βx, and ℜ0>1 when cβ exceeds the value cβ∗=0.747247. From Theorem 2, the spatial heterogeneous transmission rate βx can strongly enhance the risk of heroin spread.
In Figure4(b), we further set(53)dx=0.0125⋅1+cd⋅cosπx.
Here,cd∈0,1 measures the magnitude of spatial heterogeneity of the diffusion rate dx. The numerical results demonstrate that the effect of the diffusion rate dx is weaker than the transmission rate βx in the aspect of increasing ℜ0 more than one. It is found that the transmission rate βx strongly enhances the risk of heroin spread, while ℜ0 is less sensitive with respect to the diffusion rate dx.Figure 4
(a)ℜ0 with the transmission rate βx. λ=5.5,μ=0.01,d=0.0125,k=0.01,p=0.1, and βx=0.0001⋅1+cβ⋅cosπx; (b) the surface of ℜ0 with the transmission rate βx and diffusion rate dx. λ=5.5,μ=0.01,k=0.01,p=0.1, βx=0.0001⋅1+cβ⋅cosπx, and dx=0.0125⋅1+cd⋅cosπx.
(a)(b)Example 4.
Influence of the spatial heterogeneity inkx.
Recall thatkx denotes the relapse rate from heroin drug users with to without treatment class. We set(54)λ=6.02,μ=0.01,β=0.0001,d=0.0125,p=0.1,(55)kx=0.01⋅1+ck⋅cosπx,to check how ℜ0 will change with kx. Here, ck∈0,1 measures the magnitude of spatial heterogeneity of kx.
In Figure5(a), we can see that ℜ0 decreases with respect to ck. ℜ0<1 when ck exceeds the value ck∗=0.717217. From Theorem 2, we arrive at the conclusion that the spatially heterogeneous relapse rate kx can reduce the risk of heroin spread.
In Figure5(b), we further set(56)λ=6.01,dx=0.0125⋅1+cd⋅cosπx.
The numerical results demonstrate that the effect of the diffusion rate on increasingℜ0 is more sensitive than the relapse rate kx when ℜ0<1. It follows that the diffusion rate dx strongly enhances the risk of hroin spread, while kx will not do so.Figure 5
(a)ℜ0 with the relapse rate kx. λ=6.02,μ=0.01,β=0.0001,d=0.0125,p=0.1, and kx=0.01⋅1+ck⋅cosπx; (b) the surface of ℜ0 with a relapse rate kx and a diffusion rate dx. λ=6.01,μ=0.01,β=0.0001,p=0.1, kx=0.01⋅1+ck⋅cosπx, and dx=0.0125⋅1+cd⋅cosπx.
(a)(b)Example 5.
Influence of the spatial heterogeneity inpx. Recall that px denotes the treatment rate from without treatment class to the treatment class. We set(57)λ=5.5,μ=0.01,β=0.0001,d=0.0125,k=0.01,(58)px=0.1⋅1+cp⋅cosπx.
We shall check howℜ0 will change with respect to px. Here, cp∈0,1 represents the magnitude of spatial heterogeneity of px.
In Figure6(a), we find that ℜ0 increases in terms of cp. ℜ0>1 when cp exceeds the value cp∗=0.60101. From Theorem 2, we arrive at the conclusion that the heterogeneous treatment rate px can strengthen the risk of heroin spread.
In Figure6(b), we further set(59)λ=6.01,dx=0.0125⋅1+cd⋅cosπx.
The numerical results demonstrate that the effect of diffusion ratedx and treatment rate px on increasing the BRN is highly sensitive for both. It follows that the spatial heterogeneity of diffusion and treatment rates strongly enhance the risk of heroin spread.Figure 6
(a)ℜ0 with a treatment rate px. λ=5.5,μ=0.01,β=0.0001,d=0.0125,k=0.01, and px=0.1⋅1+cp⋅cosπx; (b) the surface of ℜ0 with a treatment rate px and a diffusion rate dx. λ=6.01,μ=0.01,β=0.0001,k=0.01, px=0.1⋅1+ck⋅cosπx, and dx=0.0125⋅1+cd⋅cosπx.
(a)(b)


## 7. Conclusions and Discussion

Our goal of this paper is to study a diffusive heroin epidemic model incorporating the spatial heterogeneity and spatial mobility of individuals. Unlike in ODEs and DDEs models, in this paper, we allow details of spatial heterogeneity in the sense that model parameters depends on the spatial variablex. This paper concerns with threshold dynamics of the model governed by BRN ℜ0, which determines whether heroin spread will occur in a bounded domain. We define ℜ0 and prove it plays a threshold role. The explicit formula of ℜ0 and global attractiveness of positive constant equilibrium in the homogeneous case are also addressed, see Figures 1 and 2.We have performed numerical simulation to verify our analytical results. Figure3 demonstrates the distribution of heroin drug users Ix,t. We find that the more larger the diffusion rate is used, the higher the prevalence becomes. Figures 4 and 6 indicate that transmission rate βx and treatment rate px strongly enhance the risk of heroin spread in the sense that ℜ0 will be larger than one when cβ and cp exceed the critical value.On the other hand, the influence of the diffusion rate on increasingℜ0 becomes remarkably stronger when compared with the relapse rate kx and the treatment rate px, see Figures 5(b) and 6(b). The implication is that the spatial heterogeneous environment and the mobility of individuals play an important role in increasing the heroin spread risk, compared to the homogeneous setting. In this regard, comparisons of reaction-diffusion models from ODE and DDE models (underestimate the heroin spread risk) may be instructive.However, since heroin spread in a population is a complex heterogeneous process, uniform diffusion rates do not allow the conclusion that diffusion rates alone are responsible for the heroin spread risk. Our model can also be extended to be a more realistic case. Practically, heroin drug users with treatment individuals shall disperse at different rates compared with susceptible and heroin drug users without treatment individuals due to the quarantine measures or control. We notice that, however, it will bring some difficulties and new challenges in the proof of the boundedness of solutions. We left this problem for future consideration.



---

*Source: 8268950-2020-05-20.xml*