Qualitative analysis on a diffusive age-structured heroin transmission model

Nonlinear Analysis: Real World Applications 54 (2020) 103105

Contents lists available at ScienceDirect

âĂĺ

Nonlinear Analysis: Real World Applications www.elsevier.com/locate/nonrwa

Qualitative analysis on a diffusive age-structured heroin transmission model Xi-Chao Duan a ,∗, Xue-Zhi Li b , Maia Martcheva c a

College of Science, University of Shanghai for Science and Technology, Shanghai 200093, China College of Mathematics and Information Science, Henan Normal University, Xinxiang 453007, China Department of Mathematics, University of Florida, 358 Little Hall, PO Box 118105, Gainesville, FL 32611–8105, United States b c

article

info

Article history: Received 25 August 2019 Received in revised form 14 January 2020 Accepted 17 January 2020 Available online xxxx Keywords: Heroin transmission model Basic reproduction number Heterogeneous environment Principle eigenvalue Threshold dynamics

abstract In this paper, to understand the impact of spatial heterogeneity of treatment and movement of individuals on the persistence and extinction of heroin spread, we propose a new diffusive heroin transmission model with treatment-dependent age-structure. The basic reproduction number in heterogenous environment R0 of the system is defined, which is consistent with the one deduced from the next generation operator approach R(x). The threshold dynamics in terms of the basic reproduction number is established: if R0 ≤ 1, the drug-free steady state is globally asymptotically stable, if R0 > 1, heroin transmission is uniformly persistent if it is present initially. In particular, when the environment is homogeneous and R0 > 1, our system has a unique space-independent drug spread steady state and it is globally asymptotically stable. Finally, some numerical simulations are carried out to illustrate the main results. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction In the past two decades, China has faced a dramatic increase in illicit drug abuse accompanying rapid economic growth and development [1]. In 2000, heroin still was the first choice among drug users (rising from 83.4% in 1993 to 95.9% in 2000). Heroin is an opiate drug that is synthesized from morphine, a naturally occurring substance extracted from the seed pod of the Asian opium poppy plant. Heroin usually appears as a white or brown powder or as a black sticky substance, known as “black tar heroin” [2]. Heroin’s most frequent routes of delivery are intravenous injection (25%) and inhalation [3]. It crosses the blood brain barrier within 15–20 s, rapidly achieving a high level syndrome in the brain and the central nervous system which causes both the ‘rush’ experienced by users and the toxicity [4]. Heroin users are at high risk for addiction and it is estimated that about 23% of individuals who use heroin become dependent on it. This brings tremendous pressures and damages to society and public health. ∗ Corresponding author. E-mail addresses: [email protected] (X.-C. Duan), [email protected] (X.-Z. Li), [email protected] (M. Martcheva).

https://doi.org/10.1016/j.nonrwa.2020.103105 1468-1218/© 2020 Elsevier Ltd. All rights reserved.

2

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

In fact, the spread of heroin habituation and addiction can be well modeled by epidemic type models as “transmission” occurs in the form of peer pressure where established users recruit susceptible individuals into trying and using the drug [5–7]. As a result, mathematical modeling can provide a general insight into the dynamics of classes of drug users, and as such, could become a useful device to aid specialist teams in devising treatment strategies. Modeling heroin addiction and spread in epidemic fashion is not new [8]. Most of heroin models are ODE models and the population is divided into three classes, namely the number of susceptibles, S(t), the number of drug users not in treatment, U1 (t) and the number of drug users in treatment, U2 (t), respectively. More recently, to reveal the role of treatment-age in heroin spread dynamics, Fang et al. [9] proposed an age-structured heroin transmission model and established its global dynamics. Reaction–diffusion epidemic models with spatially heterogeneous environment can also play important role in modeling disease spread. Allen et al. have proposed a diffusive SIS epidemic model and studied the roles of spatial heterogeneity on the persistence and extinction of the disease in [10]. Recently, the authors in [11] proved global stability results of the steady states for a diffusive SIS epidemic model with varying total population. For more information about such models, one can refer to references [12–14] and related literatures therein. At present, there is little work on diffusive heroin transmission modeling. Note that the drug users under treatment cannot move freely and the treatment rate can vary according to location. It is necessary to formulate a diffusive heroin transmission model and study its dynamic behaviors. To account for the time-since-treatment of the treated individuals U2 (t), and the movement of the individuals, in this paper, we present a diffusive heroin transmission model with a treatment age. In this system, we assume the individuals in S(t) and U1 (t) can move freely and the individuals in treatment cannot move, the heroin relapse rate depends on how long an individual has been in treatment and a treatment rate depends on the spatial location. As a result, in this paper we consider the following system: ⎧ ∂S ⎪ = δ∆S + Λ − µS(x, t) − βS(x, t)U1 (x, t), ⎪ ⎪ ⎪ ∂t ⎪ ∫ ∞ ⎨ ∂U 1 = δ∆U1 + βS(x, t)U1 (x, t) − (µ + ν1 + γ(x))U1 (x, t) + k(θ)U2 (x, θ, t)dθ, (1.1) ⎪ ∂t 0 ⎪ ⎪ ⎪ ⎪ ⎩ ∂U2 (x, θ, t) + ∂U2 (x, θ, t) = −(µ + ν + k(θ))U (x, θ, t), 2 2 ∂θ ∂t for t > 0, x ∈ Ω , with the initial and boundary conditions: { ¯, U2 (x, 0, t) = γ(x)U1 (x, t), x∈Ω ¯, S(x, 0) = S0 (x), U1 (x, 0) = U10 (x), U2 (x, θ, 0) = U20 (x, θ), x ∈ Ω

(1.2)

and homogeneous Neumann boundary conditions ∂S ∂U1 = 0, = 0, t > 0, x ∈ ∂Ω . ∂n ∂n

(1.3)

In system (1.1)–(1.2), θ is the treatment age, that is the time that has elapsed since the drug users entered treatment; U2 (x, θ, t) is the density of drug users in treatment with age θ at location x and at time t; S(x, t) is the density of the susceptibles at location x and at time t; U1 (x, t) is the density of drug users not in treatment, never treated and relapsed drug users. Ω is a bounded domain in Rn with smooth boundary ∂Ω , ∂/∂n denotes the outward normal derivative on ∂Ω . The boundary conditions in (1.3) imply that the susceptible individuals and drug users do not move across the boundary ∂Ω . The parameter δ is the diffusion coefficient. The positive constant Λ is the recruitment of susceptibles, µ is the natural death rate of the general population, β is the addiction coefficient per contact per unit time. The parameter ν1 is a removal

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

3

rate that includes drug-related deaths of users not in treatment and a spontaneous recovery rate (individuals not in treatment who stop using drugs but are no longer susceptible), ν2 is a removal rate that includes the drug-related deaths of users in treatment and a rate of successful “cure” that corresponds to recovery to a drug-free life and immunity to drug addiction for the duration of the modeling time period. The function k(θ) is the probability of a drug user in treatment with treatment age θ relapsing to the untreated drug user; the function γ(x) is the rate of drug users who enter treatment at location x. Throughout the paper, we make the following assumptions: (A1) It is assumed that S0 and U10 are nonnegative continuous functions on Ω ; (A2) The initial condition U20 (x, θ) is uniformly bounded for θ ∈ (0, +∞) for each fixed x; (A3) The map θ → k(θ) is almost everywhere bounded and belongs to L∞ older continuous function on Ω and + ((0, +∞), R) \ {0L∞ }; (A4) The function γ(x) is a positive H¨ minx∈Ω¯ γ(x) = γ∗ > 0. To the best of our knowledge, there are few results on the global stability for age-structured reaction– diffusion epidemic models in the literature, which is what we address here. The paper is organized as follows. In the next section, we present some preliminary results for system (1.1)–(1.2). In Section 3, we discuss the threshold dynamics in terms of the basic reproduction number. In Section 4, we prove the uniform persistence of our system and the existence of a drug spread steady state. In Section 5, we prove global stability of the unique drug spread steady state when the parameter γ is a constant and the basic reproduction number is larger than one. Finally, in Section 6, a brief discussion and some numerical examples are presented. 2. Preliminary results In this section, we present some preliminary results of our system (1.1)–(1.2). In the following, we will prove the well-posedness, the uniform boundedness of the solutions of system (1.1)–(1.2) and its solution semiflow has a global attractor. Let C := U C(Ω , R) be the set of all bounded and uniformly continuous functions from Ω to R. C+ := U C(Ω , R+ ) be the set of all bounded and uniformly continuous functions from Ω to R+ , C 2 (Ω , R) denotes the set of all functions in U C whose derivatives up to order 2 are continuous, C 2 (Ω , R+ ) denotes the set of all functions in U C+ whose derivatives up to order 2 are continuous. Define X := C 2 (Ω , R) × C 2 (Ω , R) × L1 ((0, ∞), C), X+ := C 2 (Ω , R+ ) × C 2 (Ω , R+ ) × L1+ ((0, ∞), C+ ). Then we will show that the solutions of system (1.1)–(1.2) are uniformly bounded. Lemma 2.1. There exists a positive constant C1 independent on initial data such that the solution (S, U1 , U2 ) of (1.1) satisfies ∫ ∞     ∥S(·, t)∥L∞ (Ω) + ∥U1 (·, t)∥L∞ (Ω) +  U (·, θ, t)dθ ≤ C1 . (2.1) 2   ∞ 0

L

Proof . Let ∫ N (x, t) := S(x, t) + U1 (x, t) +

(Ω)

∞

U2 (x, θ, t)dθ 0

be the total number of individuals at location x ∈ Ω (= (0, L)) at time t, and ) ∫ ( ∫ ∞ ˆ (t) := N S(x, t) + U1 (x, t) + U2 (x, θ, t)dθ dx Ω

0

4

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

be the total number of individuals in Ω (= (0, L)) at time t. By adding the equations in (1.1) and integrating on Ω , we have that ) ∫ ∞ ˆ (t) ∫ ∂ ( dN = S(x, t) + U1 (x, t) + U2 (x, θ, t)dθ dx dt ∂t 0 ) ∫Ω ( ∫ ∞ ∂ ∂U2 ∂ = S + U1 + dθ dx ∂t ∂t ∂t 0 )) ( ∫Ω ∫ ( ∫ ∞ U2 (x, θ, t)dθ dx = (δ∆S + δ∆U1 )dx + Λ − µ S(x, t) + U1 (x, t) + 0 Ω Ω ∫ ( ∫ ( ∫ ∞ ⏐θ=∞ ) ) ⏐ − dx − γ(x)U1 (x, t) + U2 (x, θ, t)⏐ ν1 U1 (x, t) + ν2 U2 (x, θ, t)dθ dx θ=0 Ω 0 ∫ Ω ∫ ( ) ≤ ∆(δS + δU1 )dx + Λ − µN (x, t) dx Ω ∫Ω ∂ ˆ (t) = (δS + δU1 )dx + Λ|Ω | − µN ∂n ∂Ω ˆ (t). =Λ|Ω | − µN It then follows that

∫ ˆ (t) ≤ N (x, t)dx = lim sup N

lim sup t→+∞

t→+∞

Ω

Λ |Ω |, µ

ˆ (t) ≤ Λ |Ω | is satisfied for some t = t0 ∈ R, then which is independent of the initial data. Furthermore, if N µ it is satisfied for all t ≥ t0 . Therefore, system (1.1) is point dissipative (see [12, Definition 2.1.]). □ Lemma 2.2.

For the given continuous function S0 (x) in (1.2), we denote S0∗ ≜ max S0 (x).

(2.2)

x∈Ω

Then the solution (S, U1 , U2 ) of (1.1) satisfies S(x, t) ≤ Λ+ ,

∀x ∈ Ω ,

t > 0.

(2.3)

Proof . It follows from the first equation of (1.1) that we have ∂S − δ∆S = Λ − µS(x, t) − βS(x, t)U1 (x, t) ≤ Λ − µS(x, t). ∂t Observe that the solution S(t) of system (1.1) and max{Λ/µ, S0∗ } are a pair of sub- and super solutions to the initial–boundary value problem: ⎧ ⎪ ⎪ ∂v = δ∆v + Λ − µv(x, t), x ∈ Ω , t > 0, ⎪ ⎪ ⎨ ∂t ∂v (2.4) = 0, t > 0, x ∈ ∂Ω , t > 0, ⎪ ⎪ ∂n ⎪ ⎪ ⎩ v(x, 0) = S (x), x ∈ Ω . 0 { } Λ ∗ Since Λ is an equilibrium solution of (2.4). Then the maximum principle implies that v(x, t) ≤ max , S 0 . µ µ Thus, the comparison principle gives { } Λ ∗ S(x, t) ≤ v(x, t) ≤ max ,S ≜ Λ+ . µ 0 This completes the proof.

□

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

We define the set { ∫ ( ∫ D := (S, U1 , U2 ) ∈ X : S + U1 + Ω

0

∞

} Λ + . U2 dθ dx ≤ |Ω |, and S ≤ Λ µ

5

)

(2.5)

It follows from Lemmas 2.1 and 2.2 that we have the following Proposition 2.1. Assume that the set D, defined in (2.5), is positively invariant for the solution semiflow Ψ (t) of system (1.1)–(1.2), in the sense that Ψ (t)ψ ∈ D, ∀ t ≥ 0, ψ ∈ D. In what follows, we reformulate system (1.1)–(1.2) as a nonhomogeneous Cauchy problem in order to apply ˆ := C 2 (Ω , R) × C 2 (Ω , R) × C × L1 ((0, ∞), C) integrated semigroup theory. To this end, we pair the space X with the usual supremum norm, i.e., ∫ ∞ ∥χ∥X |U2 (x, θ)|dθ, ˆ = ∥S∥C 2 + ∥U1 ∥C 2 + sup x

0

ˆ where for χ ∈ X, ∥S∥C 2 = sup |S| + sup |Sx | + sup |Sxx |, ∥U1 ∥C 2 = sup |U1 | + sup |U1x | + sup |U1xx |. x

x

x

x

x

x

Set ˆ0 := C 2 (Ω , R) × C 2 (Ω , R) × {0C } × L1 ((0, ∞), C), X ˆ+ := C 2 (Ω , R+ ) × C 2 (Ω , R+ ) × C+ × L1+ ((0, ∞), C+ ), X ˆ0+ := X ˆ0 ∩ X ˆ+ . X ( To simplify matters, we define a vector U =

( S(·), U1 (·),

0 U2 (·, ·)

))T ˆ0 . The linear operator ∈ X

ˆ is defined as follows: A : Dom(A) → X ⎛

⎞ δ∆S − µS ⎜ ⎟ δ∆U1 − (µ + ν1 )U1 ( ) ⎟, AU = ⎜ ⎝ ⎠ −U2 (0) ∂U2 − ∂θ − (µ + ν2 + k(θ))U2

¯ )×C 2 (Ω ¯ )×{0}×W 1,1 ((0, +∞), U C) ⊂ X, ˆ where W 1,1 is a Sobolev space of functions with Dom(A) = C 2 (Ω whose first derivative in θ is integrable. To apply Theorem 6.2 in [15], it is easy to see that the solution of dU dt = AU, U(0) = 0, is U(t) ≡ 0 and Appendix A implies that ∥U(t)∥X ˆ ≤ c(t)∥U(0)∥X ˆ for U(0) = ψ ∈ X+ and t ≥ 0. From Theorem 6.2 in Thieme [15], it follows that the linear operator A generates a non-degenerate ˆ We define a map F : X ˆ0 → X ˆ by integrated semigroup T (t) on X. ⎛ ⎞ Λ − βS(x)U1 (x) ∫ +∞ ⎜ ⎟ ⎜ βS(x)U1 (x) + k(θ)U2 (x, θ)dθ − γ(x)U1 (x) ⎟ ⎟. F (U)(x) = ⎜ (2.6) ⎜ ⎟ 0( ) ⎝ ⎠ γ(x)U1 (x) 0L1 Then we can rewrite system (1.1)–(1.2) as the following nonhomogeneous abstract Cauchy problem dU(t) = AU(t) + F (U(t)), t > 0, dt

(2.7)

ˆ0+ . with U(0) = U 0 ∈ X From Proposition 4.16 in [16], we need the following lemma (see, for instance [17]) to obtain the well-posedness of system (2.7).

6

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

ˆ0 , we have Lemma 2.3. Assume that F is defined in (2.6) for ϕ ∈ X ⎞ ⎛ −βψ2 (x)ϕ1 (x) − βψ1 (x)ϕ2 (x) ∫ F ′ [ψ](ϕ) = ⎝ βψ2 (x)ϕ1 (x) + βψ1 (x)ϕ2 (x) + 0+∞ k(θ)ϕ3 (x, θ)dθ − γ(x)ϕ2 (x) ⎠ , γ(x)ϕ2 (x) ˆ0 . and then F is continuously Fr´echet differentiable on X Proposition 2.2. For any U0 in Dom(A), system (2.7) has a unique continuously differentiable solution ˆ0 such that U : [0, T0 ) → X ∫ t U(t) = etA U0 + et−s F (U(s))ds 0

and either T0 → +∞ or if T0 < ∞, the solution blows up. Lemma 2.4. Ω × [0, T0 ).

Let U(t, x; U0 ) be a solution of system (2.7). If U0 is positive, then U(t, x; U0 ) is positive on

Based on the above results and Proposition 2.1, we can conclude that the solution Ψ (t)(U0 ) = U(t, x; U0 ) is a global solution for all t ≥ 0 when U0 ∈ D. In the following, for the sake of convenience, we let ∫ ∞ ∫ ∞ ∫θ − (µ+ν2 +k(τ ))dτ Π (θ) = e 0 , Π (θ)dθ = K1 , and k(θ)Π (θ)dθ = K . (2.8) 0

0

Using integration, U2 (x, θ, t) satisfies the following Volterra formulation: ⎧ if t ≥ θ, ⎨ γ(x)U1 (x, t − θ)Π (θ), Π (θ) U2 (x, θ, t) = , if θ ≥ t. ⎩ U2 (x, θ − t, 0) Π (θ − t)

(2.9)

Lemma 2.5. The solution semiflow Ψ (t) is K-contracting in the sense that limt→+∞ K(Ψ (t)(B)) = 0 for any bounded set B ⊂ X+ . Proof . Recall the Kuratowski measure of noncompactness (see [18]), which is defined by K(B) := inf{r : B has a finite cover of diameter < r}

(2.10)

for any bounded set B. Let K(B) = +∞ whenever B is unbounded. It is easy to see that B is precompact ¯ is compact) if and only if K(B) = 0. (i.e., B Let B be a bounded subset in X+ , we need only to show that Ψ (t) is asymptotically compact on B in the sense that for any sequences ϕn ∈ B and tn → ∞, there exist subsequences ϕnk ∈ B and tnk → ∞ such that Ψ (tnk )(ϕnk ) converges in X+ as k → ∞ (see, Lemma 4.1 in [19]). Let ¯ . Since Sn (x, tn ) and U1n (x, tn ) (Sn (x, tn ), U1n (x, tn ), U2n (x, ·, tn )) = Ψ (t)(ϕn )(x), ∀ϕn ∈ X+ , t ≥ 0, x ∈ Ω are twice continuously differentiable with respect to x and (Sn , U1n , U2n ) ∈ B, we have that {Sn (x, tn )}n≥1 ¯ . It follows from (2.9) and similar arguments to those in Lemma and {U1n (x, tn )}n≥1 are equicontinuous on Ω ¯ . In fact, from (2.9), we have that 2.3 of [20] that we have {U2n (x, θ, tn )}n≥1 equicontinuous on Ω ∫ +∞ |U2n (x, θ, tn ) − U2n (y, θ, tn )|dθ ∫ 0tn (2.11) ≤ γ(x)Π (θ)|U1n (x, tn − θ) − U1n (y, tn − θ)|dθ 0 ∫ +∞ Π (θ) + |U2n (x, θ − tn , 0) − U2n (y, θ − tn , 0)| dθ. Π (θ − tn ) tn

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

7

Since U2n (x, θ, 0) is uniformly bounded, the last term clearly converges to 0 as t → +∞. It follows from assumption (A4) that γ(x) is H¨ older continuous. It further follows that {U1n (x, tn )}n≥1 is equicontinuous and we have that the last two terms of (2.11) converge to 0 as t → +∞. Therefore, Ψ (t) is asymptotically compact on B. Then from the proof of Lemma 4.1 in [19], we can deduce that the solution semiflow Ψ (t) is K-contracting with respect to the Kuratowski measure defined in (2.10), i.e., lim K(Ψ (t)(B)) = 0, t→+∞

for any bounded set B ⊂ X+ . □ Proposition 2.1 indicates that the semiflow Ψ (t) is point dissipative. It follows from Proposition 2.1, Lemma 2.5, and by Theorem 2.6 in [21] that we have the following: Theorem 2.1. The solution semiflow Ψ (t) admits a global attractor A, which attracts all bound subsets of X+ . 2.1. The basic reproduction number To study the dynamics of our system, we need to introduce the basic reproduction number R0 . According to the definition of the basic reproduction number in the existing literature [10,22,23], we define the basic reproduction number R0 as R0 =

βΛ µ

sup 0̸=ϕ∈H 1 (Ω)

(

(µ + ν1 ) + (1 − K )

∫ Ω

γ(x)ϕ2 (x)dx +

2

∫ Ω

δ|∇ϕ(x)| dx

).

(2.12)

Indeed, one can use the renewal process to obtain existence of the basic reproduction number [17] (see Appendix B). On the other hand, it is well known that the elliptic problem − δ∆x S = Λ − µS, x ∈ Ω ,

∂S = 0, x ∈ ∂Ω ∂n

(2.13)

0 0 admits a unique positive solution S 0 = Λ µ . Then E = (S , 0, 0) is a spatially homogeneous steady state of (1.1), which we call the drug-free steady state (DFS). In this paper, the stability of E 0 is completely determined by R0 .

Lemma 2.6. (i) R0 is a monotone decreasing function of δ with R0 → max x∈Ω

and R0 →

βΛ µ µ + ν1 + (1 − K )γ(x)

βΛ µ µ + ν1 + (1 − K )γ

as δ → 0

as δ → ∞ with γ =

1 |Ω |

∫ γ(x)dx. Ω

∗ (ii) If β Λ µ < µ + ν1 + (1 − K )γ (high-treatment domain), then there exists a threshold value δ > 0 such that R0 < 1 for δ > δ ∗ and R0 > 1 for δ < δ ∗ . (iii) If β Λ µ > µ + ν1 + (1 − K )γ (low-treatment domain), then R0 > 1 for δ > 0.

8

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

3. Stability of the drug-free steady state To study the stability of the DFS, we will consider first an eigenvalue problem which is associated with (1.1). The long term behavior of system (1.1) is determined by considering the total trajectories on the global compact attractor of system (1.1) [see 24, Smith and Thieme, p. 243]. Then we linearize (1.1) around the drug free steady state E 0 to obtain ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

Λ ∂η = δ∆η − µη(x, t) − β y(x, t), x ∈ Ω , ∂t µ ∫ ∞ ∂y Λ = δ∆y + β y(x, t) − (µ + ν1 + γ(x))y(x, t) + k(θ)z(x, θ, t)dθ, x ∈ Ω , ∂t µ 0 ∂z ∂z ¯ + = −(µ + ν2 + k(θ))z(x, θ, t), x ∈ Ω ∂t ∂θ ¯ z(x, 0, t) = γ(x)y(x, t), x∈Ω ∂η ∂y = = 0, x ∈ ∂Ω . ∂n ∂n

(3.1)

λt where η(x, t) = S(x, t) − Λ µ , y(x, t) = U1 (x, t) and z(x, θ, t) = U2 (x, θ, t). Suppose that η(x, t) = e ϕ(x), λt λt y(x, t) = e φ(x) and z(x, θ, t) = e ψ(x, θ) with λ ∈ R. We substitute this solution into the linearized equations and divide by eλt to get the following eigenvalue problem ⎧ Λ ⎪ ⎪ λϕ(x) = δ∆ϕ − µϕ(x) − β φ(x), x ∈ Ω , ⎪ ⎪ µ ⎪ ⎪ ∫ ∞ ⎪ ⎪ Λ ⎪ ⎪ ⎪ λφ(x) = δ∆φ + β φ(x) − (µ + ν1 + γ(x))φ(x) + k(θ)ψ(x, θ)dθ, x ∈ Ω , ⎪ ⎨ µ 0 (3.2) ∂ψ ⎪ = −(µ + ν2 + k(θ))ψ(x, θ) λψ(x, θ) + ⎪ ⎪ ⎪ ∂θ ⎪ ⎪ ψ(x, 0) = γ(x)φ(x), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ϕ = ∂φ = 0, x ∈ ∂Ω . ∂n ∂n

Note that the first equation of (3.2) is decoupled from the rest equations of (3.2). We consider the following system for the drug spread related compartments ∫ ∞ ⎧ ∂y Λ ⎪ ⎪ = δ∆y + β y(x, t) − (µ + ν + γ(x))y(x, t) + k(θ)z(x, θ, t)dθ, x ∈ Ω , 1 ⎪ ⎪ ∂t µ ⎪ 0 ⎪ ⎪ ⎪ ∂z ⎨ ∂z ¯ + = −(µ + ν2 + k(θ))z(x, θ, t), x ∈ Ω ∂t ∂θ ⎪ ⎪ ¯ ⎪ z(x, 0, t) = γ(x)y(x, t), x∈Ω ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂y = 0, x ∈ ∂Ω . ∂n

(3.3)

Then the generalized version of the system (3.3) is obtained as follows: ∫ ∞ ⎧ ∂y ⎪ ⎪ = δ∆y + βT (x)y(x, t) − (µ + ν1 + γ(x))y(x, t) + k(θ)z(x, θ, t)dθ, x ∈ Ω , ⎪ ⎪ ∂t ⎪ 0 ⎪ ⎪ ⎪ ∂z ⎨ ∂z ¯ + = −(µ + ν2 + k(θ))z(x, θ, t), x ∈ Ω ∂t ∂θ ⎪ ⎪ ¯ ⎪ z(x, 0, t) = γ(x)y(x, t), x∈Ω ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂y = 0, x ∈ ∂Ω , ∂n

(3.4)

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

9

¯ . It is easy see that (3.4) is a cooperative system. It follows form (3.4) that the where T (x) > 0, for ∀ x ∈ Ω associated eigenvalue problem: ∫ ∞ ⎧ ⎪ ⎪ λφ(x) = δ∆φ + βT (x)φ(x) − (µ + ν1 + γ(x))φ(x) + k(θ)ψ(x, θ)dθ, x ∈ Ω , ⎪ ⎪ ⎪ 0 ⎪ ⎪ ⎪ ∂ψ ⎨ = −(µ + ν2 + k(θ))ψ(x, θ) λψ(x, θ) + (3.5) ∂θ ⎪ ⎪ ⎪ ψ(x, 0) = γ(x)φ(x), ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂φ = 0, x ∈ ∂Ω . ∂n The following lemma presents the existence of the principle eigenvalue of (3.5). ¯ , the eigenvalue problem (3.5) has a principle eigenvalue, denoted by Lemma 3.1. For T (x) > 0, ∀ x ∈ Ω λ(T ) which is associated with a strongly positive eigenfunction. Proof . Denote Y := C 2 (Ω , R)×L1 ((0, ∞), C). For every initial value function ϕ = (ϕ1 , ϕ2 ) ∈ Y, the solution map Φt : Y → Y associated with the linear system (3.4) is defined by Φt (ϕ) = (y(·, t; ϕ), z(·, ·, t : ϕ)), ∀ ϕ ∈ Y, t ≥ 0. Following the same line of reasoning as in the proof of Lemma 2.5, we can obtain that for each t > 0, Φt is an K-contraction on Y in the sense that lim K(Φt (B)) = 0

t→∞

for any bounded set B in Y and K, defined in (2.10). Based on the above discussions, it is easy to see that the solution map Φt generated by (3.4) is K-condensing in the sense that K(Φt (B)) < K(B), for any bounded set B in Y with K(B) > 0, t > 0. By the generalized Krein–Rutman Theorem (see, e.g., [25, Lemma 2.2] and [26, Chapter II.14]), Eq. (3.5) has a principal eigenvalue, denoted by λ(T ), with an associated eigenvector ϑ∗ = (φ∗ , ψ ∗ ) ≫ 0. □ Choosing T (x) =

λ µ,

by Lemma 3.1, we have that the principal eigenvalue of (3.5) can be denoted as λ∗ := λ(T )|T = λ .

(3.6)

µ

After the above preparation, we now turn to the eigenvalue problem (3.2) to prove the local stability of E . From the last two equations of (3.2), we have 0

∗

ψ(x, θ) = γ(x)e−λ θ Π (θ)φ(x). Substituting (3.7) into the second equation of (3.2), it follows that ∫ ∞ ∗ Λ δ∆φ∗ + β φ∗ (x) − (µ + ν1 + γ(x))φ∗ (x) + γ(x) k(θ)e−λ θ Π (θ)dθφ∗ (x) = λ∗ φ∗ , µ 0

(3.7)

(3.8)

∗

for x ∈ Ω and ∂φ ∂n = 0 for x ∈ ∂Ω . From the expression of the basic reproduction number R0 in (2.12), we have that there exists a positive function Φ(x) ∈ C 2 (Ω ) such that − δ∆Φ(x) + (µ + ν1 + (1 − K )γ(x))Φ(x) −

βΛ/µ Φ(x) = 0, x ∈ Ω R0

(3.9)

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

10

= 0, x ∈ ∂Ω (similarly as in Lemma 2.3 in [10]). Recall that φ∗ and Φ are both positive on with ∂Φ(x) ∂n ∗ Ω and that ∂φ∂n(x) = ∂Φ(x) = 0 on ∂Ω . We multiply (3.8) by Φ and (3.9)) by φ∗ , add them together and ∂n integrate by parts on Ω to obtain: ( )∫ ∫ ∫ 1 Λ 1− φ∗ (x)Φ(x)dx + pˆ γ(x)φ∗ (x)Φ(x)dx = λ∗ φ∗ (x)Φ(x)dx β µ R0 Ω Ω Ω where

∫ pˆ =

∗

k(θ)e−λ θ Π (θ)dθ − K .

0

Since Ω γ(x)φ (x)Φ(x)dx and Ω φ (x)Φ(x)dx are both positive, we conclude that R0 − 1 and λ∗ have the same sign. In fact, if we assume that λ∗ > 0 when R0 − 1 < 0, then pˆ < 0 and therefore λ∗ < 0, which leads to a contradiction. If we assume that λ∗ < 0 when R0 − 1 > 0, then pˆ > 0 and therefore λ∗ > 0, which also leads to a contradiction. If we further assume that λ∗ ̸= 0 when R0 − 1 = 0, we also obtain a contradiction. In conclusion, we have the following result. ∫

∗

∞

Lemma 3.2.

∫

∗

R0 − 1 and λ∗ have the same sign.

We now show that the stability of the DFS E 0 is determined entirely by R0 . Lemma 3.3.

The DFS E 0 is locally stable if R0 < 1.

Proof . Suppose first that R0 < 1. We will show that the DFS is linearly stable. That is, if (λ, ϕ, φ) is any solution of (3.2), with at least one of ϕ or φ not identically zero on Ω , then Re(λ) must be negative. (i) If φ ≡ 0 and ϕ ̸= 0, it then follows from the first equation of (3.2) that we have: λϕ(x) = δ∆ϕ − µϕ(x), x ∈ Ω , and

∂ϕ = 0, x ∈ ∂Ω . ∂n

(3.10)

According to [27, pp. 147–148], the eigenvalue problem (3.10) has a principal eigenvalue, denoted by λ0 , with an associated eigenvector ϕ0 ≫ 0 and λ0 < 0. (ii) If ϕ ≡ 0 and φ ̸= 0, then (3.2) has no solution. (iii) If ϕ ̸= 0 and φ ̸= 0, it also follows from (3.8) that Reλ ≤ λ∗ < 0. The eigenfunction ϕ is a solution of the following equation Λ ∂ϕ λ∗ ϕ(x) = δ∆ϕ − µϕ(x) − β φ∗ (x), x ∈ Ω , with = 0, x ∈ ∂Ω . µ ∂n

(3.11)

If λ∗ = λ0 , Eq. (3.11) has no solution, then the principal eigenvalue of (3.10) λ0 is also the principal eigenvalue of (3.2). If λ∗ ̸= λ0 , Eq. (3.11) can be solved and the principle eigenvalue of (3.2) is max{λ∗ , λ0 } < 0. We conclude that the DFS E 0 is locally asymptotically stable if R0 < 1. □ The main purpose of this section is to establish the threshold dynamics of (1.1), which states that heroin use will vanish from the population if R0 < 1. Theorem 3.1.

Let (S, U1 , U2 ) be the unique solution of (1.1). If R0 < 1, then ∫ ∞ ( ) 0 lim S(x, t) − S = 0, lim U1 (x, t) = 0, lim U2 (x, θ, t)dθ = 0 t→∞

uniformly for x ∈ Ω .

t→∞

t→∞

0

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

11

Proof . Assume that R0 < 1. We will use the proof of [10, Lemma 2.5] to conclude that U1 (·, t) → 0, uniformly on Ω , as t → ∞. To begin we focus on the last two equations of system (1.1). We have that ∫ ∞ ⎧ ∂U 1 ⎪ k(θ)U2 (x, θ, t)dθ, = δ∆U1 + βS(x, t)U1 (x, t) − (µ + ν1 + γ(x))U1 (x, t) + ⎪ ⎪ ⎪ ∂t 0 ⎪ ⎪ ⎪ ⎪ ∂U2 (x, θ, t) ∂U2 (x, θ, t) ⎪ ⎪ + = −(µ + ν2 + k(θ))U2 (x, θ, t), x ∈ Ω ⎨ ∂θ ∂t U2 (x, 0, t) = γ(x)U1 (x, t), x ∈ Ω ⎪ ⎪ ⎪ ⎪ ⎪ ∂U1 ⎪ ⎪ = 0, t > 0, x ∈ ∂Ω ⎪ ⎪ ⎪ ∂n ⎩ U1 (x, 0) = U10 (x), U2 (x, θ, 0) = U20 (x, θ). It follows from S(x, t) ≤

where

Λ µ

(3.12)

in (2.3) that U1 (x, t) ≤ I1 (x, t), and U2 (x, θ, t) ≤ I2 (x, θ, t),

(3.13)

∫ ∞ ⎧ ∂I Λ 1 ⎪ = δ∆I + β I (x, t) − (µ + ν + γ(x))I (x, t) + k(θ)I2 (x, θ, t)dθ, ⎪ 1 1 1 1 ⎪ ⎪ ∂t µ 0 ⎪ ⎪ ⎪ ⎪ ∂I2 (x, θ, t) ∂I2 (x, θ, t) ⎪ ⎪ + = −(µ + ν2 + k(θ))I2 (x, θ, t), x ∈ Ω ⎨ ∂θ ∂t ⎪ I2 (x, 0, t) = γ(x)I1 (x, t), x ∈ Ω ⎪ ⎪ ⎪ ⎪ ∂I1 ⎪ ⎪ = 0, t > 0, x ∈ ∂Ω ⎪ ⎪ ⎪ ⎩ ∂n I1 (x, 0) = I10 (x) ≥ U10 (x), I2 (x, θ, 0) = U20 (x, θ).

(3.14)

Let I1 (x, t) = eλt φ(x), I2 (x, θ, t) = eλt ψ(x, θ) with λ ∈ R. We substitute this solution into (3.14) and divide by eλt to get the following linear eigenvalue problem ⎧ ∫ ∞ Λ ⎪ ⎪ λφ(x) = δ∆φ + β φ(x) − (µ + ν + γ(x))φ(x) + k(θ)ψ(x, θ)dθ, x ∈ Ω , ⎪ 1 ⎪ ⎨ µ 0 ∂ψ(x, θ) (3.15) = −(λ + µ + ν2 + k(θ))ψ(x, θ), x ∈ Ω , ⎪ ⎪ ∂θ ⎪ ⎪ ⎩ ψ(x, 0) = γ(x)φ(x, t), x ∈ Ω . In view of (1.3), we must have that ∂φ = 0, x ∈ ∂Ω . ∂n Solving the second and the third equations of (3.15), we have that ψ(x, θ) = γ(x)φ(x)e−λθ Π (θ). Substituting (3.16) into the second equation of (3.15), we have the following eigenvalue problem ⎧ Λ ⎪ ⎪ λφ(x) = δ∆φ + β φ(x) − (µ + ν1 + γ(x))φ(x) ⎪ ⎪ µ ∫ ∞ ⎪ ⎨ + γ(x) k(θ)e−λθ Π (θ)dθφ(x), x ∈ Ω , ⎪ ⎪ 0 ⎪ ⎪ ∂φ ⎪ ⎩ = 0, x ∈ ∂Ω , ∂n ∗

(3.16)

(3.17)

which agrees with (3.8). Then we can define I1 (x, t) = M eλ t φ∗ (x), where λ∗ < 0 by Lemma 3.2, φ∗ > 0 on Ω , and M is chosen so large that U1 (x, 0) ≤ I1 (x, 0) for every x ∈ Ω . Since λ∗ < 0 we have

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

12

∗

I1 (x, t) = M eλ t φ∗ (x) → 0 as t → ∞ for every x ∈ Ω . From (3.13), we have U1 (x, t) → 0 as t → ∞ for every x ∈ Ω . Next, we claim that S(x, t) → Λ/µ as t → ∞ for every x ∈ Ω . Given any small constant ϵ > 0, there exists a T > 0 and large enough, such that 0 ≤ U1 (x, t) ≤ ϵ for all x ∈ Ω , t > T . From the first equation in (1.1) and Lemma 2.2, it is easy to observe that U1 is a supersolution to ⎧ Λ ∂v ⎪ ⎪ = δ∆v + Λ − µv(x, t) − β ϵ, x ∈ Ω , t > T, ⎪ ⎪ ∂t µ ⎨ ∂v (3.18) = 0, t > 0, x ∈ ∂Ω , t > T, ⎪ ⎪ ⎪ ∂n ⎪ ⎩ v(x, T ) = S(x, T ), x ∈ Ω. and a subsolution to

⎧ ∂v ⎪ ⎪ = δ∆v + Λ − µv(x, t), x ∈ Ω , t > T, ⎪ ⎪ ⎨ ∂t ∂v = 0, t > 0, x ∈ ∂Ω , t > T, ⎪ ⎪ ∂n ⎪ ⎪ ⎩ v(x, T ) = S(x, T ), x ∈ Ω.

(3.19)

Let v1 (ϵ, x, t) and v2 (x, t) be the solutions of problems (3.18) and (3.19), respectively. It is easily checked that Λ−βΛ Λ µϵ v1 (ϵ, x, t) → , v2 (x, t) → , as t → ∞. µ µ Since v1 (ϵ, x, t) ≤ S(x, t) ≤ v2 (x, t), we have that Λ−βΛ µϵ µ

≤ S(x, t) ≤

Λ µ

for t large enough and any ϵ. Hence, we can obtain that S(x, t) → proof. □

Λ µ

= S 0 as t → ∞. This completes the

4. Existence of the drug spread steady state In this section, we will prove the existence of the drug spread steady state based on the persistence results. Lemma 4.1.

Suppose U(x, t; ϕ) is the solution of system (1.1) with U(·, 0; ϕ) = ϕ ∈ X+ .

(i) If there exists some t0 ≥ 0, such that S(·, t0 ; ϕ) ̸≡ 0, U1 (·, t0 ; ϕ) ̸≡ 0, then S(·, t; ϕ) > 0, U1 (x, t0 ; ϕ) > 0, ¯ t > t0 ; ∀x∈Ω ¯ t > t0 ; (ii) If there exists some t0 ≥ 0, such that U1 (·, t0 ; ϕ) ̸≡ 0, then U2 (·, ·, t; ϕ) > 0, ∀ x ∈ Ω ¯ (iii) For any ϕ ∈ X+ , we always have S(·, t; ϕ) > 0, ∀ x ∈ Ω t > t0 and lim inf S(·, t; ϕ) ≥ t→∞

Λ , µ + 2βS 0

where S 0 = Λ/µ. Proof . It is easy to establish that S and U1 satisfy the following { ∂S(x,t) ≥ δ∆S(x, t) − µS(x, t), x ∈ Ω , t > 0, ∂t ∂S ∂n

and

{

∂U1 (x,t) ≥ ∂t ∂U1 ∂n = 0,

= 0, x ∈ ∂Ω , t > 0 δ∆U1 (x, t) − (µ + ν1 + γ∗ )U1 (x, t), x ∈ Ω , t > 0, x ∈ ∂Ω , t > 0.

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

13

By the similar arguments as in [28, Lemma 2.1] and [28, Proposition 3.1], it follows from the strong maximum principle (see, e.g., [29, p. 172, Theorem 4]) and the Hopf boundary lemma (see, e.g., [29, p. 170, Theorem 3]) that part (i) is valid. From (2.9), we have U2 (x, θ, t) = γ(x)U1 (x, t − θ)Π (θ) + U2 (x, θ − t, 0)

Π (θ) . Π (θ − t)

This indicates that part (ii) is valid. From Proposition 2.1, it follows that there exists a t1 > 0 such that U1 (x, t) ≤ 2S 0 , ∀ t > t1 . From the first equation of (1.1), we get { ∂S(x,t) ≥ δ∆S(x, t) + Λ − µS(x, t) − 2βS 0 S(x, t), x ∈ Ω , t > 0, ∂t ∂S ∂n = 0, x ∈ ∂Ω , t > 0 Thus, we have lim inf S(·, t; ϕ) ≥ t→∞

This completes the proof.

Λ . µ + 2βS 0

□

Now we prove uniform persistence of our system (1.1). Define two sets ∂D0 = {U ∈ D|U10 (·) ≡ 0} and D0 = {U ∈ D|U10 (·) ̸≡ 0} Lemma 2.4 shows that Ψ (t)D0 ⊂ D0 , ∀ t ≥ 0 and that D0 = D\∂D0 . The semiflow Ψ (t) is called to be persist with respect to (D0 , ∂D0 ) if there exists a constant ϵ > 0 such that limt→∞ dist(Ψ (t)U, ∂D0 ) > ϵ for all U ∈ D0 . Lemma 4.2. If R0 > 1, then system (1.1) is uniformly weekly persistent in the sense that there exists a constant ϵ > 0 such that lim sup S(·, t) > ϵ, lim sup U1 (·, t) > ϵ, t→∞

t→∞

lim sup U2 (·, θ, t) > ϵΠ (θ)γ(·), t→∞

¯. for all x ∈ Ω ¯ , θ, t ∈ R+ . System (1.1) can be rewritten as: Proof . Suppose U0 ∈ ∂D0 , for x ∈ Ω ⎧ ∂S ⎪ ⎪ ⎪ ∂t = δ∆S(x, t) + Λ − µS(x, t), ⎪ ⎪ ⎪ ⎨ U (x, t) = 0, 1 ∂U ∂U2 (x, θ, t) ⎪ 2 (x, θ, t) ⎪ ⎪ + = −(µ + ν2 + k(θ))U2 (x, θ, t), ⎪ ⎪ ∂θ ∂t ⎪ ⎩ U2 (x, 0, t) = 0.

(4.1)

¯ , U1 (x, t) and U2 (x, θ, t) converge It follows that S(x, t) converges uniformly to S 0 as t → ∞ for each x ∈ Ω uniformly to 0L1 . Thus the omega limit set of the orbits {Ψ (t)U : t ≥ 0} is the singleton set {E 0 }. Suppose U0 ∈ D0 . Then, we can prove that   lim sup Ψ (t)U0 (·) − E 0 ∞ ≥ ϵ t→∞

14

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

for any ϵ > 0 and W S (E 0 ) ∩ D0 = ∅ where W S (E 0 ) be the manifold of E 0 . Assume by contradiction that (n) there exist some U0 ∈ D0 such that   1   (n) lim sup Ψ (t)U0 (·) − E 0  < . n ∞ t→∞ Then there exists a time t0 such that S0 −

1 1 ≤ S (n) (x, t) ≤ S 0 + , n n

(n)

U1 (x, t) ≤

1 n

¯ , we have for t > t0 . From the first equation of system (1.1), for t > t0 and x ∈ Ω ∂t S (n) (x, t) ≥ δ∆S (n) (x, t) + Λ − ε0 − µS (n) (x, t) with ε0 = β n1 S 0 + (

(4.2)

) 1

. Integrating Eq. (4.2), we have ∫ t ∫ ) Λ − ε0 ( ) Λ − ε0 ( 1 − e−µt ≥ 1 − e−µt0 S (n) (x, t) ≥ e−µθ Γ1 (x, y, θ)(Λ − ε0 )dθdy = µ µ 0 Ω n

where Γ1 is the Green function. From the second and the third equation of system (1.1), it follows that ⎧ (n) ) (n) ⎪ Λ − ε0 ( ∂U1 (n) (n) ⎪ ⎪ ≥ δ∆U1 + β 1 − e−µt0 U1 (x, t) − (µ + ν1 + γ(x))U1 (x, t) ⎪ ⎪ ∂t µ ⎪ ⎪ ∫ ∞ ⎪ ⎪ ⎪ (n) ⎪ + k(θ)U2 (x, θ, t)dθ, ⎪ ⎪ ⎪ ⎪ 0 ⎪ ⎪ ⎨ ∂U (n) (x, θ, t) ∂U (n) (x, θ, t) (n) 2 2 + = −(µ + ν2 + k(θ))U2 (x, θ, t), x ∈ Ω (4.3) ⎪ ∂θ ∂t ⎪ ⎪ (n) (n) ⎪ ⎪ U2 (x, 0, t) = γ(x)U1 (x, t), x ∈ Ω ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ∂U1(n) ⎪ ⎪ = 0, t > 0, x ∈ ∂Ω ⎪ ⎪ ∂n ⎪ ⎪ ⎩ (n) (n) U1 (x, 0) = U10 (x), U2 (x, θ, 0) = U20 (x, θ). (n)

(n)

Choose small enough m > 0 such that U1 (·, t0 ) ≥ mφ(x), U2 (x, θ, t) ≥ mψ(x, θ), where φ(x) and ψ(x, θ) are the solutions of the following system ∫ ∞ ⎧ ) Λ − ε0 ( −µt0 ⎪ ⎪ 1 − e φ(x) − (µ + ν + γ(x))φ(x) + k(θ)ψ(x, θ)dθ, x ∈ Ω , λφ(x) = δ∆φ + β 1 ⎪ ⎪ µ ⎪ 0 ⎪ ⎪ ⎪ ⎨ ∂ψ(x, θ) = −(λ + µ + ν2 + k(θ))ψ(x, θ), x ∈ Ω , (4.4) ∂θ ⎪ ⎪ ψ(x, 0) = γ(x)φ(x, t), x ∈ Ω ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂φ(x) = 0, t > 0, x ∈ ∂Ω . ∂n If R0 > 1, then there exists a λm > 0 corresponding to a positive eigenvector (φ(m) (x), ψ (m) (x, θ)). In fact, meλm t (φ(m) (x), ψ (m) (x, θ)) is a subsolution of system (4.3). Applying comparison principle, one obtains (n)

U1 (x, t) ≥ meλm t φ(m) (x), (n)

(n) ¯ , a, t ∈ R+ . U2 (x, θ, t) ≥ meλm t ψ (m) (x, θ), x ∈ Ω

(n)

Then, U1 (x, t) and U2 (x, θ, t) tend to infinity as t → ∞, which is a contradiction with the boundedness of the solution. Thus, there exists ϵ > 0 such that ¯. lim sup S(·, t) > ϵ, lim sup U1 (·, t) > ϵ, for all x ∈ Ω t→∞

t→∞

¯ . This completes the proof. □ It follows from (2.9) that we have U2 (·, θ, t) > ϵΠ (θ)γ(·), for all x ∈ Ω

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

15

From Lemma 2.5 and Theorem 3.1, we have the following theorem. Theorem 4.1. If R0 > 1, system (1.1) has at least one positive steady state E ∗ (x)= (S ∗ (x), U1∗ (x), U2∗ (x, θ)) . Further, system (1.1) is uniformly strongly persistent in D0 , in the sense that, there exists ϵ > 0, such that lim inf S(·, t) > ϵ, lim inf U1 (·, t) > ϵ, t→∞

t→∞

lim inf U2 (·, θ, t) > ϵΠ (θ)γ(·). t→∞

In addition, system (1.1) has a global attractor A0 ⊂ A. Proof . Define a continuous function ρ : X+ → [0, +∞) ρ(U0 ) := min U1 (x, 0), U0 ∈ X+ . ¯ x∈Ω

It follows from Lemma 4.1 that ρ−1 (0, +∞) ⊆ D0 and ρ has the property that if ρ(U0 ) > 0 or U0 ∈ D0 with ρ(U0 ) = 0, then ρ(Ψt (U0 )) > 0, ∀ t > 0. This means that ρ is a distance function for the semiflow Ψt (see, for instance [30]). Theorem 3.1 and Lemma 4.2 imply that E 0 is isolated and W S (E 0 ) ∩ D0 = ∅. Consequently, there is no cycle from E 0 to itself. Theorem 3 in [30] shows that there exists a positive constant ϵu such that for all U0 ∈ D0 min ρ(ϕ) > ϵu , ∀ U0 ∈ D0 , ϕ∈ω(U0 )

where ω(U0 ) is the omega limit set of the orbit {Ψ (t)U : t ≥ 0}. Therefore, lim inf t→∞ U1 (·, x) > ϵu , ∀ U0 ∈ D0 . One obtains from Lemma 4.2 that there exists a positive constant 0 < ϵ < ϵu such that for all x∈Ω lim inf S(·, t) > ϵ, lim inf U1 (·, t) > ϵ. t→∞

t→∞

Eq. (2.9) leads to lim inf U2 (·, θ, t) > ϵΠ (θ)γ(·). t→∞

Proposition 2.1 implies that the semiflow Ψ (t) is point dissipative. Theorem 3.7 in [29], shows that there exists a global attractor attracting all the points from D0 . Since the set D0 is convex and the solution maps Ψ (t) are K-contracting (see Lemma 2.5), it then follows from [21, Theorem 4.7] that Ψ (t) has a steady state E ∗ (·) ∈ D0 . Furthermore, Lemma 4.1 implies that E ∗ (·) is a positive steady state of system (1.1). This completes the proof. □ 5. Global stability of the drug spread steady state This section is devoted to the study of the global stability of the drug spread steady state (DSS) of (1.1) when the treatment rate is assumed to be a positive constant γ, i.e., γ(x) ≡ γ. Theorem 3.1 implies that system (1.1) admits no DSS if R0 < 1. Thus, it is easy to see that (S 0 , 0, 0) is the unique DFS of (1.1) and E ∗ = (S ∗ , U1∗ , U2∗ (θ)) is the unique DSS if and only if βΛ/µ > (µ + ν1 + γ(1 − K )), i.e., R0 = which satisfies

βΛ/µ > 1, µ + ν1 + γ(1 − K )

⎧ 0 = Λ − µS ∗ − βS ∗ U1∗ , ⎪ ⎪ ∫ ⎪ ⎪ ⎪ ⎪ ⎨ 0 = +βS ∗ U1∗ − (µ + ν1 + γ)U1∗ + ∂U2∗ (θ) ⎪ ⎪ = −(µ + ν2 + k(θ))U2∗ (θ), ⎪ ⎪ ∂θ ⎪ ⎪ ⎩ ∗ U2 (x, 0) = γU1∗ ,

0

(5.1)

∞

k(θ)U2∗ (θ)dθ, (5.2)

16

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

It then follows that

µ 1 Λ , U1∗ = (R0 − 1), U2∗ (θ) = γU1∗ Π (θ). R0 µ β In what follows, we shall construct suitable Lyapunov functionals to investigate the global stability of the unique DSS. S∗ =

Theorem 5.1. Suppose R0 is defined in (5.1). If R0 > 1, the unique space-independent drug spread steady state E ∗ is globally asymptotically stable. Proof . It follows from the third and the fourth equations of (5.2) that we have: U2∗ (θ) = γU1∗ Π (θ).

(5.3)

To prove the global stability of E ∗ , we need the following notations. Let ∫ ∞ ∫ ∞ ∫σ Π (σ) − (µ+ν2 +k(τ ))dτ θ dσ. α(θ) = k(σ)e dσ = k(σ) Π (θ) θ θ Note that α(θ) > 0 for all θ ∈ (0, ∞) and α(0) = K . Then we define ( ) S(x, t) VS (x, t) = S ∗ g , S∗ ) ( U1 (x, t) ∗ , VU1 (x, t) = U1 g U1∗ ) ( ∫ ∞ U2 (x, θ, t) ∗ dθ. VU2 (x, t) = α(θ)U2 (θ)g U2∗ (θ) 0 The function g(x) = x − 1 − ln x is a Volterra-type function and has the property that g(x) > 0 and g(1) = 0 for x > 0. Then we have that dVU1 S∗ U1∗ dVS =1− and =1− . dS S(x, t) dU1 U1 (x, t) For simplicity, set f1 (x, t) = Λ − µS(x, t) − βS(x, t)U1 (x, t), ∫ f2 (x, t) = βS(x, t)U1 (x, t) − (µ + ν1 + γ)U1 (x, t) +

∞

k(θ)U2 (x, θ, t)dθ. 0

Using the fact Λ = βS ∗ U1∗ + µS ∗ , we have ( ) dVS S∗ f1 = 1 − (Λ − βS(x, t)U1 (x, t) − µS(x, t)) dS S(x, t) ( ) [ ( ) ( )] S∗ = 1− µ S ∗ − S(x, t) + βS ∗ U1∗ − βS(x, t)U1 (x, t) S(x, t) ( )( ) µ(S(x, t) − S ∗ )2 S∗ S(x, t) U1 (x, t) =− + βS ∗ U1∗ 1 − 1− . S(x, t) S(x, t) S∗ U1∗

(5.4)

Next, we have dVU1 f2 = dU1

(

U1∗ 1− U1 (x, t)

){ ∫ ( ) βS(x, t)U1 (x, t) − µ + ν1 + γ U1 (x, t) + 0

Note that from the second equation of (5.2), µ + ν1 + γ =

1 U1∗

( ∫ βS ∗ U1∗ + 0

∞

) k(θ)U2∗ (θ)dθ .

∞

} k(θ)U2 (x, θ, t)dθ .

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

17

So we have dVU1 f2 dU1

){ U1∗ βS(x, t)U1 (x, t) U1 (x, t) } ) ∫ ∞ ( ∫ ∞ U1 (x, t) ∗ ∗ ∗ k(θ)U (x, θ, t)dθ k(θ)U (θ)dθ + βS U + − 2 2 1 U1∗ 0 0 ( ) ) { ( U1∗ S(x, t) U1 (x, t) U1 (x, t) = 1− − βS ∗ U1∗ U1 (x, t) S∗ U1∗ U1∗ ( ) } ∫ ∞ U2 (x, θ, t) U1 (x, t) + k(θ)U2∗ (θ) − dθ U2∗ (θ) U1∗ 0 ( = 1−

(5.5)

= βS(x, t)U1 (x, t) − βS ∗ U1 (x, t) − βS(x, t)U1∗ + βS ∗ U1∗ ) ( ∫ ∞ U1∗ U2 (x, θ, t) U2 (x, θ, t) U1 (x, t) − − + 1 dθ. + k(θ)U2∗ (θ) U2∗ (θ) U1∗ U1 (x, t) U2∗ (θ) 0 Now we turn to the derivative of VU2 (x, t), ∂VU2 ⏐⏐ ⏐ ∂t (1.1) ( ) ∫ ∞ ∂ U2 (x, θ, t) = α(θ)U2∗ (θ) · g dθ ∂t U2∗ (θ) 0 ) ( ∫ ∞ ∂U2 (x, θ, t) U2∗ (θ) dθ = α(θ) 1 − U2 (x, θ, t) ∂t 0 } ( ){ ∫ ∞ ) U2∗ (θ) ∂U2 (x, θ, t) ( = − µ + ν2 + k(θ) U2 (x, θ, t) dθ α(θ) 1 − − U2 (x, θ, t) ∂θ 0 ( } ) { ∫ ∞ ∗ ) U2 (θ) ∂θ U2 (x, θ, t) ( =− α(θ) 1 − U2 (x, θ, t) + µ + ν2 + k(θ) dθ U2 (x, θ, t) U2 (x, θ, t) 0 ( ){ } ∫ ∞ ) U2 (x, θ, t) ∂θ U2 (x, θ, t) ( =− α(θ)U2∗ (θ) · − 1 + µ + ν + k(θ) dθ, 2 U2∗ (θ) U2 (x, θ, t) 0

where ∂θ U2 (x, θ, t) denotes Note that ∂ g ∂θ

(

(5.6)

∂U2 (x, θ, t) . ∂θ

U2 (x, θ, t) U2∗ (θ)

)

( =

U2∗ (θ) 1− U2 (x, θ, t)

)

∂ · ∂θ

(

U2 (x, θ, t) U2∗ (θ)

)

∗ ) ∂U2 (x, θ, t) U ∗ (θ) − U2 (x, θ, t) ∂U2 (θ) 2 ∂θ ∂θ = 1− · ( ∗ )2 U2 (x, θ, t) U2 (θ) ( ) U2∗ (θ) = 1− U2 (x, θ, t) [ ( ) ] ∂θ U2 (x, θ, t)U2∗ (θ) − U2 (x, θ, t) − µ + ν2 + k(θ) U2∗ (θ) · ( ∗ )2 U2 (θ) ( ) { } ) U2 (x, θ, t) ∂θ U2 (x, θ, t) ( = − 1 · + µ + ν + k(θ) . 2 U2∗ (θ) U2 (x, θ, t)

(

U2∗ (θ)

(5.7)

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

18

Substituting (5.7) into (5.6) and using integration by parts, we obtain ) ( ∫ ∞ ∂VU2 ⏐⏐ ∂ U2 (x, θ, t) ∗ α(θ)U (θ) · = − dθ g ⏐ 2 ∂t (1.1) ∂θ U2∗ (θ) 0 ( ) ⏐θ=∞ U2 (x, θ, t) ⏐⏐ ∗ = −α(θ)U2 (θ)g ⏐ U2∗ (θ) ) ( θ=0 ∫ ∞ ( ) U2 (x, θ, t) ∂ + g α(θ)U2∗ (θ) dθ ∗ U2 (θ) ∂θ 0 ( )⏐ ( ) U2 (x, 0, t) U2 (x, θ, t) ⏐⏐ ∗ + α(0)U = −α(θ)U2∗ (θ)g (x, 0)g 2 ⏐ U2∗ (θ) U2∗ (x, 0) ) ( θ=∞ ∫ ∞ ( ) U2 (x, θ, t) ∂ g α(θ)U2∗ (θ) dθ. + U2∗ (θ) ∂θ 0 Note that

) ∂ ( d ∂ ∗ α(θ) · U2∗ (θ) + α(θ) · U2 (θ) α(θ)U2∗ (θ) = ∂θ dθ ∂θ [ ( ) ] = α(θ) µ + ν2 + k(θ) − k(θ) · U2∗ (θ) [ ( ) ] + α(θ) · − µ + ν2 + k(θ) U2∗ (θ)

(5.8)

(5.9)

= −k(θ)U2∗ (θ), and

∫ α(0) =

∞

k(θ)Π (θ)dθ,

U2∗ (x, 0) = γU1∗ ,

U2 (x, 0, t) = γU1 (x, t).

(5.10)

0

Substituting (5.9) and (5.10) into (5.8), we also have ( )⏐ ∂VU2 ⏐⏐ U2 (x, θ, t) ⏐⏐ ∗ = −α(θ)U2 (θ) · g ⏐ ⏐ ∂t (1.1) U2∗ (θ) ( θ=∞ ) ∫ ∞ U1 (x, t) + γU1∗ k(θ)Π (θ)dθ · g U1∗ ) ∫ ∞ 0( U2 (x, θ, t) − g k(θ)U2∗ (θ)dθ. ∗ (θ) U 0 2 Meanwhile, note that ( ) dVS ∂S dVS S∗ dVS ∂VS ⏐⏐ = = (δ∆S + f1 ) = 1 − δ∆S + f1 , ⏐ ∂t (1.1) dS ∂t dS S(x, t) dS ( ) dVU1 ∂U1 dVU1 dVU1 ∂VU1 ⏐⏐ U1∗ = = (δ∆U1 + f2 ) = 1 − δ∆U1 + f2 . ⏐ ∂t (1.1) dU1 ∂t dU1 U1 (x, t) dU1 Adding these terms together, we have ( ) ∂VU1 ∂VU2 ⏐⏐ ∂VS + + ⏐ ∂t ∂t ∂t (1.1) ( ) ( ) ∗ S U1∗ = 1− δ∆S(x, t) + 1 − δ∆U1 (x, t) S(x, t) U1 (x, t) ( ) µ(S(x, t) − S ∗ )2 S∗ S(x, t) − − βS ∗ U1∗ −2+ S(x, t) S(x, t) S∗ ( ) ∫ ∞ ∗ U1 U2 (x, θ, t) − k(θ)U2∗ (θ)g dθ U U ∗ (θ) 1 (x, t) 0 ( )⏐ 2 U2 (x, θ, t) ⏐ −α(θ)U2∗ (θ) · g ⏐ U2∗ (θ) θ=∞ [( ) ( ) ] ∗ S U1∗ ≤δ 1− ∆S(x, t) + 1 − ∆U1 (x, t) . S(x, t) U1 (x, t)

(5.11)

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

Define V (t) =

∫ (

) VS (x, t) + VU1 (x, t) + VU2 (x, t) dx.

19

(5.12)

Ω

Then, ) ∫ ( ∂VU1 ∂VU2 ∂VS dV (t) ⏐⏐ = + + dx ⏐ dt (1.1) ∂t ∂t ∂t Ω ) ( ) ] [( ∫ U1∗ S∗ ≤ ∆S(x, t) + 1 − ∆U1 (x, t) dx. δ 1− S(x, t) U1 (x, t) Ω Applying the Green’s first identity with a Neumann boundary condition, we have ⎡ ⎤ )2 )2 ∫ n ( n ( ∗ ∗ ∑ ∑ S ∂S(x, t) U ∂U (x, t) dV (t) ⏐⏐ 1 1 ⎣ ⎦ dx ≤ 0. ≤ −δ dx + ⏐ dt (1.1) (S(x, t))2 j=1 ∂xj (U1 (x, t))2 j=1 ∂xj Ω Therefore,

dV (t) dt

(5.13)

≤ 0, and from (5.11)–(5.13), the equality holds if and only if S(x, t) = S ∗ ,

U1 (x, t) = U1∗

and U2 (x, θ, t) = U2∗ (θ).

It then follows that the largest positively invariant subset of {(S, U1 , U2 ) ∈ A0 : dV (t)/dt = 0} is the singleton {E ∗ }. By LaSalle’s invariance principle [Theorem 4.2, 31], we obtain that E ∗ attracts all solutions of (1.1). The steady state E ∗ is globally asymptotically stable. □ 6. Simulations and discussions In this paper, we formulate a diffusive heroin transmission model with treatment age. The model allows the susceptibles and the drug users not in treatment to move freely but the drug users in treatment have a fixed spatial location. Besides that, the treatment rate of drug users depends on the location x, the heroin relapse rate of the individuals in treatment depends on the treatment age. From the existing literature [22,23], we define the basic reproduction number R0 for our heroin transmission model, which is consistent with the one obtained from the use of the renewal equation. By using measure theory, we obtained the smoothness of the solutions for our model and the existence of the principle eigenvalue for the local stability of the drug-free steady state. Then the threshold dynamics of our model are proved, i.e., if R0 ≤ 1, then the drug-free steady state is globally asymptotically stable, if R0 > 1, our model is uniformly persistent and there exists at least one drug spread steady state. When the treatment rate γ(x) is a constant, we proved the space-independent drug spread steady state is globally stable by the use of Lyapunov functionals. From the expression of the basic reproduction number R0 and Lemma 2.6, we find that the diffusion coefficient δ is a key factor in the reproduction number. Under the above threshold results, we see that the diffusion coefficient and heterogeneous treatment rate γ(·) play important roles in the dynamics of our model. Specifically, larger value of δ is more benefit for the control of heroin transmission. In what follows, to verify the theoretical results obtained in this paper, we present some numerical simulations. For simplicity, we assume that Ω = [0, L] and take a day as the unit time. Note that the function θ → k(θ) is almost everywhere bounded and belongs to L∞ + ((0, +∞), R) \ {0L∞ }. In this section we assume that the rate ⎧ ⎨ 0.1667(θ − 10)2 e−0.3(θ−10) , 10 < θ ≤ 40; k(θ) = (6.1) 0.0185, 40 < θ < 50; ⎩ 0, otherwise. Let the treatment rate of drug users at location x be γ(x) = γ0 + 0.1 sin(2πx),

(6.2)

20

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

Fig. 1. If R0 = 0.8969 < 1, the drug free steady state E 0 of system (1.1) is globally asymptotically stable with γ0 = 0.8 in (6.2) and initial conditions S0 (0) = 300, S0 (x) = 0 for x ∈ (0, 1], U10 ≡ 50 for x ∈ [0, 1].

for x ∈ [0, 1] and γ0 will be chosen later with γ0 ≥ 0.1. The other parameters in system (1.1) are as follows Λ = 102 , β = 7 × 10−4 , µ = 0.1, ν1 = 0.05, ν2 = 0.02, δ = 0.02.

(6.3)

We assume that there are no individuals in treatment in the initial time, i.e., the initial condition U20 (x, θ) ≡ 0, for x ∈ [0, 1] and θ ∈ [0, ∞). Besides that, with appropriate initial conditions, we will perform some simulations of the solutions of system (1.1) by Matlab. Choosing γ0 = 0.8, and the parameters adopted in (6.1), (6.2) and (6.3), we calculated that R0 = 0.8969 < 1. From Theorem 3.1-(i), we know that the drug-free steady state E 0 is globally asymptotically stable (see Fig. 1) and the heroin-transmission goes extinct in the end. When the treatment rate is homogeneous, there is a unique space-independent drug spread steady state E ∗ . Choosing γ(x) ≡ 0.5, we calculate that R0 = 1.2867 > 1. From Theorem 5.1, we have that the solutions

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

21

Fig. 2. If R0 = 1.2867 > 1, the drug spread steady state E ∗ of system (1.1) is globally asymptotically stable with γ(x) ≡ 0.5 for all x ∈ [0, 1] and initial conditions S0 (0) = 300, S0 (x) = 0, U10 (0) = 20, U10 (x) = 0 for x ∈ (0, 1].

of system (1.1) converge to space-independent drug spread steady state E ∗ (see Fig. 2), and the heroin transmission persists. To illustrate the role of spatially heterogeneous environment when R0 > 1, we have preformed some simulations (shown in Fig. 3), which demonstrate that the solutions of system (1.1) are uniformly bounded and seem to approach fixed profiles. Our theoretical results and numerical simulations support some useful insight into the heroin drug spread. However, there are some limitations, such as the only spatially explicit parameter is the treatment rate γ(x) in our system; the individuals in treatment cannot move freely; the movement of the susceptibles and the drug users follows a Fickian law. Other parameters in our model can be assumed as spatially explicit according to the necessities of the real world, and nonlocal diffusion can also be introduced in the heroin transmission model. These aspects can be worked as a valuable research direction and we leave them for the future.

22

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

Fig. 3. If R0 = 1.2867 > 1, the solutions of system (1.1) seem to approach fixed profiles with γ(x) adopted as (6.2) (γ0 = 0.5) and the initial conditions S0 (0) = 300, S0 (x) = 0 for x ∈ (0, 1], U10 ≡ 50 for x ∈ [0, 1].

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The authors would like to thank the handling editor and the referees for their valuable comments that have greatly improved the presentation and content of the paper. This work is supported partially by the

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

23

National Natural Science Foundation of China (11771017); M. Martcheva is supported partially through National Science Foundation grant DMS-1515661. Appendix A In this Appendix, we will prove that the solution of dU(t) = AU(t), U(0) = ψ ∈ X+ dt with Neumann boundary conditions (1.3) satisfies ∥U(t)∥X ˆ ≤ c(t)∥U(0)∥X ˆ with c(t) > 0 not depending on U(0). In fact, we need only to prove that ∥S∥C 2 ≤ c1 (t)∥S0 ∥C 2 where c1 (t) is a function of t that will be determined later. Now, we consider the following initial–boundary value problem ⎧ ∂S ⎪ = δ∆S − µS, 0 < x < L, ⎪ ⎪ ⎨ ∂t (A.1) S(x, 0) = S0 , 0 ≤ x ≤ L, ⎪ ⎪ ⎪ ⎩ Sx (0, t) = 0, Sx (L, t) = 0. ˆ t). Then S(x, ˆ t) should satisfy the following initial–boundary value problem Set S(x, t) = e−µt S(x, ⎧ ∂ Sˆ ⎪ ˆ 0 < x < L, ⎪ = δ∆S, ⎪ ⎪ ⎨ ∂t ˆ 0) = S0 , 0 ≤ x ≤ L, S(x, ⎪ ⎪ ⎪ ⎪ ⎩ Sˆx (0, t) = 0, Sˆx (L, t) = 0. Looking for separable solutions, we can solve problem (A.2) and have that ˆ t) = a0 + S(x,

∞ ∑

an e−(

nπ 2 δt L

)

cos

( nπx )

n=1

where a0 =

1 L

∫

L

S0 dx, an = 0

2 L

∫

L

S0 cos 0

( nπx ) L

L

,

dx, n = 1, 2, . . . .

It then follows that the solution of problem (A.1) S(x, t) satisfies that [ ] ∞ ( nπx ) ∑ 2 −( nπ δt −µt ) S(x, t) =e a0 + an e L cos , L n=1 [ ∞ ] ( nπx ) nπ ∑ 2 −( nπ δt −µt ) Sx (x, t) =e − an e L sin , L L n=1 [ ∞ ] ( nπx ) ( nπ )2 ∑ 2 −( nπ δt −µt Sxx (x, t) =e − an e L ) cos . L L n=1 By some simple calculations and |an | ≤ 2 supx S0 , we have that L |S(x, t)| ≤e−µt √ sup S0 , πδt x L |Sx (x, t)| ≤e−µt sup S0 , πδt x

(A.2)

24

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

L 1 |Sxx (x, t)| ≤e−µt √ sup S0 . πδt 2δt x It then follows that |S(x, t)| + |Sx (x, t)| + |Sxx (x, t)| ≤ c1 (t) sup S0 , x

and therefore ∥S(x, t)∥C 2 (Ω) ¯ ≤ c1 (t)∥S0 ∥C 2 (Ω) ¯ , where c1 (t) = e

−µt

L √ πδt

(

1 1 1+ √ + πδt 2δt

) > 0.

Similarly, the inequality ∥U1 ∥C 2 ≤ c2 (t)∥U10 ∥C 2 can be proved. It then follows from (2.9) that ∫

∫ ∞ Π (θ) |U20 (x, θ − t)| γ(x)|U1 (x, t − θ)|Π (θ)dθ + dθ Π (θ − t) t 0 ∫ t ∫ ∞ ≤γ ∗ c2 (t)∥U10 ∥C 2 Π (θ)dθ + |U20 (x, θ − t)|dθ 0 t ∫ ∞ ≤γ ∗ K1 c2 (t)∥U10 ∥C 2 + |U20 (x, θ)|dθ,

∞

∫

t

|U2 (x, θ, t)|dθ ≤ 0

0 ∗

where γ = maxx γ(x) and therefore ∫ ∞ ∫ sup |U2 (x, θ, t)|dθ ≤ γ ∗ K1 c2 (t)∥U10 ∥C 2 + sup x

x

0

∞

|U20 (x, θ)|dθ.

0

ˆ and some direct calculations, we can obtain that ∥U(t)∥ ≤ By use of the usual supremum norm of X ˆ X c(t)∥U(0)∥X for U(0) = ψ ∈ X and t ≥ 0. + ˆ Appendix B To get the renewal equation, we focus our attention to Eq. (3.3), i.e., ∫ ∞ ⎧ Λ ∂y ⎪ ⎪ = δ∆y + β y(x, t) − (µ + ν1 + γ(x))y(x, t) + k(θ)z(x, θ, t)dθ, x ∈ Ω , ⎪ ⎪ ∂t µ ⎪ 0 ⎪ ⎪ ⎪ ∂z ⎨ ∂z ¯ + = −(µ + ν2 + k(θ))z(x, θ, t), x ∈ Ω ∂t ∂θ ⎪ ⎪ ¯, ⎪ z(x, 0, t) = γ(x)y(x, t), x ∈ Ω ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂y(x, t) = 0, x ∈ ∂Ω . ∂n Solving the second and third equations, we have ⎧ ⎨ γ(x)y(x, t − θ)Π (θ), Π (θ) z(x, θ, t) = , ⎩ z(x, θ − t, 0) Π (θ − t)

(B.1)

if t ≥ θ, (B.2)

if θ ≥ t.

Substituting (B.2) into the first equation of (B.1), we have ∂y = δ∆y − c(x)y(x, t) + ∂t with the notation

∫

t

∫

∞

k(θ)γ(x)y(x, t − θ)Π (θ)dθ + 0

k(θ)z(x, θ − t, 0) t

Λ c(x) := µ + ν1 + γ(x) − β . µ

Π (θ) dθ Π (θ − t)

(B.3)

(B.4)

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

25

Solving (B.3), we have s

∫ t∫

k(θ)γ(x)y(x, s − θ)Π (θ)dθ × e(δ∆−c(x))(t−s) ds

y(x, t) = f (x, t) + 0

0

where f (x, t) = y(x, 0)e(δ∆−c(x))t +

∫ t∫

∞

k(θ)z(x, θ − s, 0) 0

s

Π (θ) dθds. Π (θ − s)

By exchanging the order of integrals, we have ∫ t ∫ t y(x, t) = f (x, t) + γ(x) k(θ)Π (θ) × y(x, s − θ)e(δ∆−c(x))(t−s) dsdθ 0 θ ∫ t ∫ t−θ = f (x, t) + γ(x) k(θ)Π (θ) × y(x, s)e(δ∆−c(x))(t−s−θ) dsdθ 0 0 ∫ t ∫ t = f (x, t) + γ(x) k(θ)Π (θ) × y(x, t − s)e(δ∆−c(x))(s−θ) dsdθ θ ∫0 t ∫ s = f (x, t) + γ(x) k(θ)Π (θ)e(δ∆−c(x))(s−θ) dθ × y(x, t − s)ds. 0

0

Then the next generation operator is ∞

∫

∫

s

k(θ)Π (θ)e(δ∆−c(x))(s−θ) dθds ∫ ∞ ∫ ∞ = γ(x) k(θ)Π (θ) × e(δ∆−c(x))(s−θ) dsdθ 0 θ ∫ ∞ ∫ ∞ = γ(x) k(θ)Π (θ)dθ × e(δ∆−c(x))s ds

R(x) = γ(x)

0

0

0

0

= K γ(x)(c(x) − δ∆)−1 )−1 ( Λ = K γ(x) µ + ν1 + γ(x) − β − δ∆ µ := P L where P =

K γ(x) µ+ν1 +γ(x)−β Λ µ

, L =

(

µ + ν1 + γ(x) − β Λ µ

)( )−1 µ + ν1 + γ(x) − β Λ − δ∆ . From Lemma 3.1 µ

in [17], R(x) is strictly positive and compact. ˜ 0 = r(P L ) is the By Krein–Rutman theorem [18], the compactness and positivity of R(x) shows that R ˜ only eigenvalue of R(x) with a positive eigenfunction. As shown in [Theorem 2.4, 32], R0 has a variational characterization form ⎫ ⎧ ∫ ⎬ ⎨ 1 ˜ 0 = sup [ ( ) ] R K γ(x)ϕ2 dx ∫ 2 ⎭ 0̸=ϕ∈H 1 Ω ⎩ δ|∇ϕ| + µ + ν + 1 + γ(x) − β Λ ϕ2 dx Ω µ

Ω

Moreover, from the basic reproduction number R0 defined in (2.12), we easily have the following equivalent ˜0: relationship between R0 and R ˜ 0 > 1 ⇔ R0 > 1, R

˜ 0 = 1 ⇔ R0 = 1, R

and

˜ 0 < 1 ⇔ R0 < 1. R

References [1] X. Li, Y. Zhou, B. Stanton, Illicit drug initiation among institutionalized drug users in China, Addiction 97 (2002) 575–582. [2] NIDA InfoFacts: Heroin, http://www.nida.nih.gov/infofacts/heroin.html.

26

X.-C. Duan, X.-Z. Li and M. Martcheva / Nonlinear Analysis: Real World Applications 54 (2020) 103105

[3] W. Hao, Z. Su, S. Xiao, et al., Longitudinal surveys of prevalence rates and use patterns of illicit drugs at selected high-prevalence areas in China from 1993 to 2000, Addiction 99 (2004) 1176–1180. [4] K.A. Sporer, Acute heroin overdose, Ann. Intern. Med. 130 (1999) 584–590. [5] C. Comiskey, National Prevalence of Problematic Opiate Use in Ireland, EMCDDA Tech. Report, 1999. [6] A. Kelly, M. Carvalho, C. Teljeur, Prevalence of Opiate Use in Ireland 2000-2001. A 3-Source Capture Recapture Study, A Report to the National Advisory Committee on Drugs, Subcommittee on Prevalence, Small Area Health Research Unit, Department of Public. [7] European monitoring centre for drugs and drug addiction (EMCDDA): Annual report, 2005, http://annualreport.emcd da.eu.int/en/homeen.html. [8] D.R. Mackintosh, G.T. Stewart, A mathematical model of a heroin epidemic: implications for control policies, J. Epidemiol. Community Health 33 (4) (1979) 299–304. [9] B. Fang, X.-Z. Li, M. Martcheva, L.-M. Cai, Global asymptotic properties of a heroin epidemic model with treat-age, Appl. Math. Comput. 263 (2015) 315–331. [10] L.J.S. Allen, B.M. Bolker, Y. Lou, A.L. Nevai, Asymptotic profiles of the steady states for an SIS epidemic reaction-diffusion model, Discrete Contin. Dyn. Syst. 21 (2008) 1–20. [11] H. Li, R. Peng, F.-B. Wang, Varying total population enhances disease persistence: Qualitative analysis on a diffusive SIS epidemicmodel, J. Differential Equations 262 (2017) 885–913. [12] R. Cui, K.-Y. Lam, Y. Lou, Dynamics and asymptotic profiles of steady states of an epidemic model in advective environments, J. Differential Equations 263 (2017) 2343–2373. [13] R. Cui, Y. Lou, A spatial SIS model in advective heterogeneous environments, J. Differential Equations 261 (2016) 3305–3343. [14] R. Peng, X. Zhao, A reaction-diffusion SIS epidemic model in a time-periodic environment, Nonlinearity 25 (2012) 1451–1471. [15] H.R. Thieme, Integrated semigroups and integrated solutions to abstract Cauchy problems, J. Math. Anal. Appl. 152 (1990) 416–447. [16] G.F. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, CRC Press, 1985. [17] J. Yang, X. Wang, Dynamics and asymptotical profiles of an age-structured viral infection model with spatial diffusion, Appl. Math. Comput. 360 (2019) 236–254. [18] K. Deimling, Nonlinear Functional Analysis, Springer-Verlag, Berlin, 1988. [19] S.B. Hsu, F.B. Wang, X.-Q. Zhao, Dynamics of a periodically pulsed bio-reactor model with a hydraulic storage zone, J. Dynam. Differential Equations 23 (2011) 817–842. [20] X. Duan, S. Yuan, K. Wang, Dynamics of a diffusive age-structured HBV model with saturating incidence, Math. Biosci. Eng. 13 (2016) 935–968. [21] P. Magal, X.-Q. Zhao, Global attractors and steady states for uniformly persistent dynamical systems, SIAM J. Math. Anal. 37 (2005) 251–275. [22] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, J. Math. Biol. 28 (4) (1990) 365–382. [23] W.-D. Wang, X.-Q. Zhao, Basic reproduction numbers for reactio-diffusion epidemic models, SIAM J. Appl. Dyn. Syst. 11 (4) (2012) 1652–1673. [24] H. Smith, H.R. Thieme, Dynamical System and Population Persistence, AMS, Providence, RI, 2010. [25] J. Jiang, X. Liang, X.-Q. Zhao, Saddle point behavior for monotone semiflows and reaction-diffusion models, J. Differential Equations 203 (2004) 313–330. [26] P. Hess, Periodic-parabolic boundary value problem and positivity, in: Pitman Res. Notes Math., vol. 247, Longman Scientific & Technical, Harlow, 1991, copublished in the United States withJohn Wiley & Sons, Inc., New York. [27] H.L. Smith, Monotone dynamical systems: An introduction to the theory of competitive and cooperative systems, in: Math. Surveys Monogr, vol. 41, American Mathematical Society Providence, RI, 1995. [28] S.-B. Hsu, J. Jiang, F.-B. Wang, On a system of reaction-diffusion equations arising from competition with internal storage in an unstirred chemostat, J. Differential Equations 248 (2010) 2470–2496. [29] M.H. Protter, H.F. Weinberger, Maximum Principles in Differential Equations, Corrected Reprint of the 1967 Original, Springer-Verlag, New York, 1984. [30] H.L. Smith, X.Q. Zhao, Robust persistence for semidynamical systems, Nonlinear Anal. 47 (2001) 6169–6179. [31] J.A. Walker, Dynamical Systems and Evolution Equations: Theory and Applications, Plenum Press, New York, 1980. [32] R.S. Cantrell, C. Cosner, Spatial Ecology via Reaction-Diffusion Equations, Wiley, USA, 2005.