Analysis of a Diffusive Heroin Epidemic Model in a Heterogeneous Environment

+is paper is concerned with a reaction-diffusion heroin model in a bound domain. +e objective of this paper is to explore the threshold dynamics based on threshold parameter and basic reproduction number (BRN) R0, and it is proved that if R0 < 1, heroin spread will be extinct, while ifR0 > 1, heroin spread is uniformly persistent and there exists a positive heroin-spread steady state. We also obtain that the explicit formula of R0 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.


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. e 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.
ere 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. e 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]. e 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]). e 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 time t, 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. e threshold dynamics and control strategies in the light of BRN are achieved. We refer the readers to [9][10][11] for more literature on other types of heroin epidemic models. e 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][13][14][15][16][17][18][19][20][21][22][23][24]). is 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]. e 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][26][27][28]).
In this paper, we use Ω ⊂ R n to denote a bounded domain, and the boundary of Ω is smooth and denoted by zΩ. 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 U(x, t), V(x, t), and W(x, t), respectively. We use the following initial conditions at location x ∈ Ω and time t � 0: To avoid overlapping writing, in the sequel, we always use the boundary condition at location x ∈zΩ and time t > 0: Let us pay attention to the following diffusive heroin model at location x ∈ Ω and time t > 0: where the symbol ∇ stands for the gradient operator; on the boundary of Ω, the label n stands for the outward normal vector; ∇ · (d(x)∇H) stands for the divergence of d(x)∇H, where H � U, V, and W, respectively; we assume all model parameters are dependent of the location x allowing for spatial heterogeneity: d(x) is the diffusion rate, λ(x) and μ(x) are the recruitment rate and death rate, respectively, and the other parameters β(x), p(x), and k(x) 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 Sections 3 and 4, according to the theory of the next generation operator, we briefly determine BRN R 0 and further investigate its threshold role: if R 0 < 1, heroin spread will be extinct, while if R 0 > 1, heroin spread is uniformly persistent. Further, when R 0 > 1, the existence of the heroin-spread steady state is also confirmed. In Section 5, we also obtain that the explicit formula of R 0 and global attractiveness of constant heroin spread equilibrium with all model parameters are supposed to be positive constants (homogeneous case). e numerical examples are provided to support and validate theoretic results and explore the effect of the spatial heterogeneity in Section 6. is paper ends with some conclusions and discussion.
By the standard procedures, we further define the nonlinear operator Q � (Q 1 , Q 2 , Q 3 ) T : X + ⟶ X as 2 Complexity for any ψ ∈ X + . Under these settings, let v � (U, V, W) T , and we rewrite system (3) as where We then have the following preliminary theorem for (3) with (2). (4). Corresponding to which, B is the linear part of system (3). As usual, we set the definition domain of B as follows: Hence, B generates a C 0 -semigroup on X. Further, as standard procedures, we can check that Following the arguments stated in [30] (Corollary 4), we immediately obtain that on the interval [0, τ ∞ ), 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 with boundary condition (2) and initial condition According to a standard argument as in [31], we conclude that (9) has a unique, positive, and global asymptotic stable steady state 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, which is point dissipative. us, Φ(t) has a global compact attractor which directly follows from the arguments as in [32] ( eorem 2.2.6).
For convenience, we denote Y ≔ C(Ω, R 2 ), and its cone is denoted by subject to boundary condition (2) for x ∈ zΩ , t > 0. Further, we define 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 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 is allows us to define the BRN by using the spectral radius of operator L: e following result comes from [34], and we omit the proof here.
Lemma 2 (see [35], Lemma 2.3). For any ψ ∈ X + , we have the following statements on the solution of model (3): , then for x ∈ Ω, the following assertion holds: Now, let us consider heroin extinction with the threshold R 0 .

Persistence of the Heroin Disease
Theorem 3. For any ψ ∈ X + , we have the following statements: (i) ere exists a small enough number δ > 0 such that the solution of model (3) where w � U, V, and W, respectively.

(ii) System (3) has at least one PSS in Ω.
Proof. To proceed further, we divide the state space X + into X 0 and zX 0 as follows: For these settings, we directly have X � X 0 ∪ zX 0 . From Lemma 2 (i), Φ(t)X 0 ⊆ X 0 directly follows.

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.
For x ∈ Ω, t > 0, we shall study the following system: with the boundary condition zU zn and initial condition Clearly, (37) has a heroin-free equilibrium E 0 � (λ/μ, 0, 0). e BRN of (37) takes the following form: which is the same as the ODE system. We denote by E * � (U * , V * , W * ) the positive endemic equilibrium. en, it satisfies Direct calculation yields Clearly, if [R 0 ] > 1, then (37) has a unique positive constant equilibrium E * .
In what follows, we give our main result in terms of [R 0 ]. (44) Here, we use the equality that 2 − A/B − B/A � − g(A/B) − g(B/A), for A, B > 0. It then follows from the above inequality that dL E * (t)/dt � 0 if and only if (U, V, W) � E * . Hence, L E * (t) defines a Lyapunov function. We can easily check that On the other hand, Υ(·, t) � U(·, t) + V(·, t) + W(·, t) is bounded. It then follows from [38] ( eorem A2]) that for some positive constant C 0 , the following inequality holds: Further, with the aid of the Sobolev embedding theorem, we get the conclusion that namely, E * attracts all solutions of (37). us, the global attractiveness of E * directly follows. □ □

Numerical Examples
is section is devoted to provide numerical examples to support our results. Further, we shall study the effect of the spatial heterogeneous environment. roughout this section, we set the domain as Ω � [− 2, 2] ⊂ R. 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.

(49)
With this setting, we can obtain that R 0 � 1.666667. From eorem 4, the solutions of model (37) approach to E * , that is, the heroin spread will be established, see Figure 1.
If we let p � 0.2 in (49), we can obtain R 0 � 0.909091. From eorem 2 (i), (1000, 0, 0) is globally asymptotically stable, that is, the heroin epidemic will be extinct, see Figure 2. respectively. Under these settings, we will investigate the situations of spatial distribution of heroin drug user I(x, t) prevalence with different diffusion rates. Figure 3 demonstrates the distribution of heroin drug users I(x, t). We find that the larger diffusion rate plays an increasingly important role in the higher prevalence.

Complexity
Example 3. Influence of the transmission rate β(x). We set 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 Figure 4(a), we find that R 0 increases in terms of c β in β(x), and R 0 > 1 when c β exceeds the value c * β � 0.747247.
From eorem 2, the spatial heterogeneous transmission rate β(x) can strongly enhance the risk of heroin spread.
In Figure 4(b), we further set Here, c d ∈ [0, 1] measures the magnitude of spatial heterogeneity of the diffusion rate d(x). e numerical results demonstrate that the effect of the diffusion rate d(x) is weaker than the transmission rate β(x) in the aspect of increasing R 0 more than one. It is found that the transmission rate β(x) strongly enhances the risk of heroin spread, while R 0 is less sensitive with respect to the diffusion rate d(x).

Example 4. Influence of the spatial heterogeneity in k(x).
Recall that k(x) denotes the relapse rate from heroin drug users with to without treatment class. We set to check how R 0 will change with k(x). Here, c k ∈ [0, 1] measures the magnitude of spatial heterogeneity of k(x).
In Figure 5(a), we can see that R 0 decreases with respect to c k . R 0 < 1 when c k exceeds the value c * k � 0.717217. From eorem 2, we arrive at the conclusion that the spatially heterogeneous relapse rate k(x) can reduce the risk of heroin spread.
In Figure 5(b), we further set λ � 6.01, e numerical results demonstrate that the effect of the diffusion rate on increasing R 0 is more sensitive than the relapse rate k(x) when R 0 < 1. It follows that the diffusion rate d(x) strongly enhances the risk of hroin spread, while k(x) will not do so.
We shall check how R 0 will change with respect to p(x). Here, c p ∈ [0, 1] represents the magnitude of spatial heterogeneity of p(x).
In Figure 6(a), we find that R 0 increases in terms of c p . R 0 > 1 when c p exceeds the value c * p � 0.60101. From eorem 2, we arrive at the conclusion that the heterogeneous treatment rate p(x) can strengthen the risk of heroin spread.
In Figure 6(b), we further set λ � 6.01, e numerical results demonstrate that the effect of diffusion rate d(x) and treatment rate p(x) 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.

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 variable x. is paper concerns with threshold dynamics of the model governed by BRN R 0 , which determines whether heroin spread will occur in a bounded domain. We define R 0 and prove it plays a threshold role. e explicit formula of R 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. Figure 3 demonstrates the distribution of heroin drug users I(x, 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 p(x) strongly enhance the risk of heroin spread in the sense that R 0 will be larger than one when c β and c p exceed the critical value.
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.

Data Availability
All data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.