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) , and it is proved that if , heroin spread will be extinct, while if , heroin spread is uniformly persistent and there exists a positive heroin-spread steady state. We also obtain that the explicit formula of 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.

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 time , a continuous variable denoted by 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 [911] 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, [1224]). 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, [2528]).

In this paper, we use 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 and time , the density of susceptible and heroin drug users without and with treatment classes is labeled by , , and , respectively. We use the following initial conditions at location and time :

To avoid overlapping writing, in the sequel, we always use the boundary condition at location and time :

Let us pay attention to the following diffusive heroin model at location and time :where the symbol stands for the gradient operator; on the boundary of , the label stands for the outward normal vector; stands for the divergence of , where , respectively; we assume all model parameters are dependent of the location allowing for spatial heterogeneity: is the diffusion rate, and are the recruitment rate and death rate, respectively, and the other parameters , and 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 and further investigate its threshold role: if , heroin spread will be extinct, while if , heroin spread is uniformly persistent. Further, when , the existence of the heroin-spread steady state is also confirmed. In Section 5, we also obtain that the explicit formula of 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 space , which equips with the supremum norm. The cone of is denoted by . Let , and its cone is denoted by . For convenience of notations, we denote . At location and time , let denote the green function of the following equation subject to the boundary condition (2). Hence, for , we can define the solution semigroup , as follows:

According to [29] (Corollary 7.2.3), is strongly positive and compact.

By the standard procedures, we further define the nonlinear operator asfor any . Under these settings, let , and we rewrite system (3) aswhere .

We then have the following preliminary theorem for (3) with (2).

Theorem 1. For any , system (3) with (2) has a unique global solution on the interval , with . Further, (3) generates a solution semiflow denoted by , which possesses a global compact attractor.

Proof. Recall that is defined in (4). Corresponding to which, is the linear part of system (3). As usual, we set the definition domain of as follows:Hence, generates a -semigroup on . Further, as standard procedures, we can check thatFollowing the arguments stated in [30] (Corollary 4), we immediately obtain that on the interval , system (3) with a Neumann boundary condition has a local solution, where . We next aim to give the proof on the uniform boundedness of the solution in , which in turn implies that the local solution is indeed a global solution. For this purpose, let , and we then add all equations of (3) to getwith 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 , where in satisfiesHence, from the comparison principle, the boundedness of solutions of (3) is confirmed on the interval . Moreover, the existence of a global solution of system (3) implies that system (3) generates a semiflow, denoted by , i.e., for ,which is point dissipative. Thus, 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 . Letting in (3), from the first equation of (3), we know that is governed by (9). It follows that . In such setting, the heroin-free steady state HFSS is denoted by .

Linearizing (3) around HFSS and further substituting the following solution into it yield the following linear system for :subject to boundary condition (2) for . Since (12) is a cooperation system, we know that the linear eigenvalue problem problem (12) owns a principal eigenvalue with a positive eigenfunction , which follows from the Krein–Rutman theorem.

For convenience, we denote , and its cone is denoted by . Let us define , , , where and are defined in (4). Hence, forms a positive -semigroup on . Precisely, for , is governed bysubject to boundary condition (2) for .

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 at time . As time evolves, at time , 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

This allows us to define the BRN by using the spectral radius of operator :

The following result comes from [34], and we omit the proof here.

Lemma 1. (i) has the same sign as ;(ii)If , the HFSS becomes unstable;(iii)If , the HFSS is locally stable.

Denoting by , , and , 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 , we have the following statements on the solution of model (3):(i)For any , if , we directly have .(ii)Further let and , then for , the following assertion holds:Now, let us consider heroin extinction with the threshold .

Theorem 2. Let be the solution of system (3) with , then we have the following statements:(i)The unique HFSS, , is globally asymptotically stable if .(ii)If , we have

Proof. We first prove (i). It follows from and Lemma 1 that . From the continuity of , we know that for , holds. Corresponding to , there is a positive eigenvector . From the -equation of (3), there is such that for and , always hold. Hence, for , and equations of (3) can be governed bywith boundary condition (2). Choose a large enough number such thatThen, for and , is a subsolution of the following linear auxiliary system:with boundary condition (2). From the comparison principle (see, for example, [29]), we obtain thatBecause , we directly have that as , uniformly for . By (9), we know that as , uniformly for .
Next, we prove (ii). From , Lemma 1, and the continuity of , we know that for , holds. Corresponding to , there is a positive eigenvector and . Suppose, for the contrary, that there exists a positive solution of (3) such thatIt follows that there exists such thatHence, for and , and equations of (3) can be governed byChoosing small enough as that with boundary condition (2). Without loss of generality, choosing small enough thatNote that for and , satisfies the following system:with boundary condition (2). According to comparison principle (see, for example, [29]), we yield that for ,Consequently, we have and as , which results in a contradiction. This proves (19).

4. Persistence of the Heroin Disease

Theorem 3. For any , we have the following statements:(i)There exists a small enough number such that the solution of model (3) for iswhere , respectively.(ii)System (3) has at least one PSS in .

Proof. To proceed further, we divide the state space into and as follows:For these settings, we directly have . From Lemma 2 (i), directly follows.
For the convenience of forthcoming discussions, we define the following set:In what follows, we denote by the omega limit set of . We shall prove that . In fact, let be given, and it follows that , which indicates that or .
In the case where for all , we see from the -equation of system (3), uniformly for . In view of the -equation of system (3), we immediately yield that uniformly for .
In the case where for some , Lemma 2 ensures that for all , . Thus, we have , . In view of the -equation of system (3), we immediately yield that uniformly for . By the -equation of system (3), it then follows that uniformly for .
As a consequence, , .
Further from (ii) of Theorem 2, it follows that if , forms a uniform weak repeller for , i.e.,Let be the continuous function defined byFrom the above definition, we directly have that , and enjoys the following claim: if either and or , then . Hence, defines a distance function for (see, for example, [36]).
Denote by the stable subset of . From the definition of , HFSS is isolated in in the sense that . Moreover, there are no sets of from a cycle in . Following the general results in [36], we can conclude that such that , that is,Combined with Lemma 2 yields , , where is a small enough number. Hence, we can choose , and the uniform persistence result holds. According to [37] (Theorem 4.7), (3) with (2) has at least one PSS in . 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.

For , we shall study the following system:with the boundary conditionand initial condition

Clearly, (37) has a heroin-free equilibrium . The BRN of (37) takes the following form:which is the same as the ODE system.

We denote by the positive endemic equilibrium. Then, it satisfies

Direct calculation yields

Clearly, if , then (37) has a unique positive constant equilibrium .

In what follows, we give our main result in terms of .

Theorem 4. If , then is globally attractive.

Proof. Definewhere , .
Direct calculation of the derivative of yieldsFurther, we haveHere, we use the equality that , for . It then follows from the above inequality that if and only if . Hence, defines a Lyapunov function. We can easily check thatOn the other hand, is bounded. It then follows from [38] (Theorem A2]) that for some positive constant , the following inequality holds:Further, with the aid of the Sobolev embedding theorem, we get the conclusion thatnamely, attracts all solutions of (37). Thus, the global attractiveness of 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 . We use labels instead of 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 Section 5. To this end, we set parameters in model (37) as follows:With this setting, we can obtain that . From Theorem 4, the solutions of model (37) approach to , that is, the heroin spread will be established, see Figure 1.
If we let in (49), we can obtain . From Theorem 2 (i), is globally asymptotically stable, that is, the heroin epidemic will be extinct, see Figure 2.

Example 2. Influence of the diffusion rate .
We setrespectively. Under these settings, we will investigate the situations of spatial distribution of heroin drug user prevalence with different diffusion rates. Figure 3 demonstrates the distribution of heroin drug users . We find that the larger diffusion rate plays an increasingly important role in the higher prevalence.

Example 3. Influence of the transmission rate .
We setHere, measures the magnitude of spatial heterogeneity. Obviously, means the spatial homogeneous case, and means the higher spatial heterogeneity.
In Figure 4(a), we find that increases in terms of in , and when exceeds the value . From Theorem 2, the spatial heterogeneous transmission rate can strongly enhance the risk of heroin spread.
In Figure 4(b), we further setHere, measures the magnitude of spatial heterogeneity of the diffusion rate . The numerical results demonstrate that the effect of the diffusion rate is weaker than the transmission rate in the aspect of increasing more than one. It is found that the transmission rate strongly enhances the risk of heroin spread, while is less sensitive with respect to the diffusion rate .

Example 4. Influence of the spatial heterogeneity in .
Recall that denotes the relapse rate from heroin drug users with to without treatment class. We setto check how will change with . Here, measures the magnitude of spatial heterogeneity of .
In Figure 5(a), we can see that decreases with respect to . when exceeds the value . From Theorem 2, we arrive at the conclusion that the spatially heterogeneous relapse rate can reduce the risk of heroin spread.
In Figure 5(b), we further setThe numerical results demonstrate that the effect of the diffusion rate on increasing is more sensitive than the relapse rate when . It follows that the diffusion rate strongly enhances the risk of hroin spread, while will not do so.

Example 5. Influence of the spatial heterogeneity in . Recall that denotes the treatment rate from without treatment class to the treatment class. We setWe shall check how will change with respect to . Here, represents the magnitude of spatial heterogeneity of .
In Figure 6(a), we find that increases in terms of . when exceeds the value . From Theorem 2, we arrive at the conclusion that the heterogeneous treatment rate can strengthen the risk of heroin spread.
In Figure 6(b), we further setThe numerical results demonstrate that the effect of diffusion rate and treatment rate 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.

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 variable . This paper concerns with threshold dynamics of the model governed by BRN , which determines whether heroin spread will occur in a bounded domain. We define and prove it plays a threshold role. The explicit formula of 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 . We find that the more larger the diffusion rate is used, the higher the prevalence becomes. Figures 4 and 6 indicate that transmission rate and treatment rate strongly enhance the risk of heroin spread in the sense that will be larger than one when and exceed the critical value.

On the other hand, the influence of the diffusion rate on increasing becomes remarkably stronger when compared with the relapse rate and the treatment rate , 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.

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

J. Wang was supported by the National Natural Science Foundation of China (no. 11871179), Natural Science Foundation of Heilongjiang Province (nos. LC2018002 and LH2019A021), and Heilongjiang Provincial Key Laboratory of the Theory and Computation of Complex Systems. H. Sun was supported by the Science and Technology Project of Jiangxi Provincial Department of Education (no. GJJ190923).