Threshold dynamics of a delayed multi-group heroin epidemic model in heterogeneous populations

The aim of this paper is to investigate the threshold dynamics of a heroin epidemic in heterogeneous populations. 
The model is described by a delayed multi-group model, which allows us to model interactions both within-group and inter-group separately. Here we are able to prove the existence of heroin-spread equilibrium and the uniform persistence of the model. The proofs of main results come from suitable applications of graph-theoretic approach to the method of Lyapunov functionals and Krichhoff's matrix tree theorem. Numerical simulations are performed to support the results of the model for the case where $n=2$.

existence, uniqueness, and the stability of equilibria, provide the theoretical basis for the specialist teams to design suitable and reasonable strategies for the control of heroin.
Most results on heroin epidemic model focus on the existence and long time behavior of solutions of single group. There are few results on the heroin epidemic in heterogeneous populations. The purpose of this paper is to establish the qualitative properties of a delayed multi-group heroin epidemic model. To this end, we would like to mention a recent work of Fang et al. in [4]. Denote by S(t), U 1 (t) and U 2 (t) the numbers of susceptible individuals, heroin users not in treatment, and heroin users undergoing treatment, at time t, respectively. The delayed heroin epidemic model formulated in [4] takes the following form where the parameters Λ, β, µ, p, δ 1 , δ 2 represent the entering flux of susceptible individuals, the rate of becoming heroin users through contacts between susceptible individuals and heroin users, the natural death rate of populations, the rate of heroin users becoming treatment individuals, the removal rate including recovery and heroin-related death of heroin users not in treatment, the removal rate including immunity to heroin addiction during treatment period and heroin-related death of heroin users in treatment. The time that a susceptible individual becomes a heroin user is described by τ1 0f (θ)U 1 (t − θ)e −(µ+δ1+p)θ dθ, where τ 1 is the maximum of the delay, the kernel functionf (θ) denotes the distribution of the infectivity of heroin users at time θ, and e −(µ+δ1+p)θ accounts the probability surviving the progression stage to be a heroin user. The time that a heroin user who is under treatment relapsing to untreated compartment is described by where τ 2 is the maximum of this delay, the kernel functionĝ(θ) denotes the distribution function over the interval [0, τ 2 ], and e −(µ+δ2)θ denotes the probability surviving untreated until relapsing to untreated compartment at time θ.
Multi-group epidemic models have been developed to explore transmission dynamics of infectious diseases in heterogeneous populations. This multi-group feature or phenomenon can be particularly shown for some diseases, such as heroin users, measles, HIV/AIDS or malaria (vector borne disease). The heterogeneity of population allows us to divide population into several homogeneous groups according to epidemiological or geographical such as communities, cities, and countries. Thus, we can model interactions both within-group and inter-group separately. Li et al. [11] investigated a class of multi-group epidemic models with distributed delays to describe the disease spread in a heterogeneous host population with general age-structure and varying infectivity. Feng et al. [5] established a general class of multi-group epidemic models with latency and relapse in heterogeneous populations. The global dynamics of equilibria of multi-group models involving basic reproduction number and a graph-theoretic approach to the method of global Lyapunov functionals can be found in the literatures. We refer to [6,10,19,20,3,16,17,25,9,24,26,23] and references therein.
Motivated by the above works, it comes natural questions that: How does finite time delay affect the heroin spread in heterogeneous populations? and whether sustained oscillations can occur by introducing heterogeneity? To this end, we will decompose the population into n homogeneous groups. Within k-th group, S k , U 1k and U 2k denote the numbers of susceptible, heroin users not in treatment, and heroin users undergoing treatment at time t, respectively.
Thus, for k = 1, · · · , n, we formulate the following delayed multi-group epidemic model (based on (1)) to describe the spread of heroin, where all parameters are assumed to be nonnegative, and their epidemic meanings are described as follows: • Λ k : the influx of individuals in k-th group; • β kj : the rate of becoming heroin users into k-th group through contacts between susceptible individuals in k-th group and heroin users in j-th group; • µ k : the natural death rate in k-th group; • p k : the rate of heroin users entering into the treatment compartment in k-th group; • δ 1k : the removal rate including recovery and heroin-related death of heroin users not in treatment in k-th group; • δ 2k : the removal rate including immunity to heroin addiction during treatment period and heroin-related death of heroin users in treatment in k-th group.
We assume that β kj is nonnegative and n-square matrix (β kj ) n×n is irreducible and provides the patterns of contact and transmission among groups, which means that for any two distinct groups k and j, individuals in U 1k can infect those in S j directly or indirectly. The term n j=1 β kj S k τ1 0f j (θ)U 1j (t − θ)e −(µj +δ1j +pj )θ dθ denotes the infection incidence in k-th group and τ2 0ĝ k (θ)U 1k (t − θ)e −(µ k +δ 2k )θ dθ denotes the rate of heroin users in treatment relapsing to heroin users in k-th group. Furthermore, for all j, k = 1, · · · , n, we always assume that τ1 0f j (θ)dθ = 1 and τ2 0ĝ k (θ)dθ = 1 hold true.
The initial conditions for model (2) take the form According to the standard theory of functional differential equations [7], model (2) has a unique solution (S 1 , U 11 , U 21 , · · · , S n , U 1n , U 2n ) with the initial conditions (3). Note that the case n = 1 is exactly the model investigated in Fang et al. [4]. The organization of this paper is as follows. In Section 2, some preliminary results are presented for model (2). The main results concerning on the global attractivity are shown in Section 3. Finally, a brief conclusion is made in Section 4.

2.
Preliminaries. This section is devoted to the positivity and boundedness of solutions, the existence of equilibria and the explicit expression of the basic reproduction number.
It follows from the first equation of system (2), we have which implies that lim sup t→∞ S k ≤ Λ k /µ k . Furthermore, adding all equations of system (2), we have For easy of presentations, we denote that and Clearly, σ k , δ k ∈ (0, 1] for all 1 ≤ k ≤ n. Note that the variable U 2k do not appear in the first and the second equations of system (2). Therefore, in the rest of the paper, we can consider the following subsystem of system (2) with the initial conditions We consider system (7) in the phase space . It is easy to check that all solutions of system (7) in X with initial conditions (8) remain nonnegative for all t ≥ 0. Similarly, it follows from the first equation of system (7) that one has that lim sup Therefore, the following set is positively invariant set of system (7). Here, we define • Γ as the interior of Γ.
Thus, summarizing the above discussions yield the following result.

Remark 1.
Here, 0 is the spectral radius of the matrix β kj σ k S 0 k D k n×n and has a clearly biological meaning. Note that 1/λ k is the average time in the untreated compartment on the first pass, and p k /λ k denotes the probability of entering the treatment compartment, and δ k is the probability of relapsing into the untreated compartment. Thus, the total average time D k in the untreated compartment at the k-th group on multiple passes Multiplying (13) by the adequate contact rate β kj σ k and the susceptible S 0 k at the k-th group produces each entry of M 0 , and thus gives 0 . It denotes the average number of new heroin users produced by one drug user not in treatment compartment at the j-th group introduced into the susceptible compartment at the k-th group.

Remark 2.
Consider the case where n = 1, according to (11), 0 reduces to Actually, 01 is exactly the same as that in [4].

Remark 3.
Consider the case where n = 2. The next generation matrix (M 0 2 ) has the following form Direct calculations show that 3. Main results. In this section, we devote to the global attractiveness of equilibria by constructing proper Lyapunov functionals and using Krichhoff's matrix tree theorem. Now, we are in the position to state our main results.
Theorem 3.1. Assume that B = (β kj ) n×n is irreducible. If 0 ≤ 1, then the HFE is unique and globally attracting in Γ and if 0 > 1, it is unstable.
Proof. We will solve this problem by using the Volterra type function which takes the following form Obviously, H(x) ≥ 0 for x > 0 and H (x) = 1 − 1/x. Thus, H(x) has its unique global minimum at x = 1 with H(1) = 0. Constructing the following Lyapunov functional Calculating the time derivative of L HF E along with the solutions of system (7), we have and similarly, Then, combining Λ k = µ k S 0 k and (14)-(15) into L HF E | (7) yield Since B is irreducible, we know that matrix M 0 is also irreducible, thus M 0 has a positive left eigenvector (ω 1 , ω 2 , ..., ω n ) corresponding to the spectral radius ρ(M 0 ) ≤ 1. Let Define
By Theorem 3.2 and using similar arguments that in [14,Lemma 4.1], the following corollary holds. Corollary 1. Assume that 0 > 1. Let (S i (t), ϕ it ) be a solution of system (7), then there exists a positive constant η > 0 such that S(t), U 1 (t) > η for all t > 0. Uniform persistence of system (7) from Theorem 3.1 and uniform boundedness of solutions in • Γ imply that there exists at least a heroin-spread equilibrium (positive equilibrium) of system (7) in • Γ (see Theorem 2.8.6 in [1]). Thus, we have the following corollary.
Finally, we have the following result for system (7). Proof. Let us mention that the global attractiveness of the HSE implies that the HSE is unique in the • Γ whenever it exists. Constructing a Volterra-type Lyapunov functional as follows Calculating the time derivative of L HSE along with the solutions of system (7) yields and similarly, we have Thus, together with (9) and (23)-(23), we have Note thatB is the Laplacian matrix of the matrix B (see [10] for more details). Since B is irreducible, matrixB also is irreducible. Let C kj denote the cofactor of the (k, j) entry ofB. We know that systemBv = 0 has a positive solution v = (v 1 , · · · , v n ), where v k = C kk > 0 for k = 1, 2, · · · , n. Define Using (24), the time derivative of V HSE along with the solutions of system (7) is Here, we denote that In the following, we will show that Φ 1 ≡ 0 and Φ 2 ≡ 0.
We firstly claim that Φ 1 ≡ 0. It follows fromBv = 0 that Therefore, it follows from (25) that which implies that Φ 1 ≡ 0 for all U 1k , k = 1, 2, · · · , n. Next we will use Krichhoff's matrix tree theorem (see [10,11]) to show that Φ 2 ≡ 0. The matrix (β kj ) n×n can be viewed as a directed graph denoted by G with vertices 1, 2, ..., n and a directed arc (k, j) from k to j iff β kj = 0. Furthermore, the set of all directed arcs of G is denoted by E(G). The directed spanning subtree of G that is rooted at vertex k is denoted by T , and Krichhoff's matrix tree theorem implies that the sum of weights of all T can be expressed as v kβkj . According to adding a directed arc (k, j) from the root k to vertex j in a subtree T , one can obtain such a unicyclic subgraph Q of G with the weight w(Q) = v kβkj . There exists a unique cycle CQ of Q and the arc (k, j) is part of the unique cycle CQ, thus, the fact that each arc of CQ is added to a corresponding rooted subtree T can form the same unicyclic subgraph Q. Therefore, the double sum in Φ 2 can be reorganized as a sum over all unicyclic subgraphs Q containing vertices 1, 2, ..., n, that is, .
Since E(CQ) denotes the set of all arcs of CQ, we have which implies that Φ Q = 0 for each Q, and thus Φ 2 ≡ 0 for each U 1k , k = 1, · · · , n.
By the Lyapunov-LaSalle invariance principle, the HSE is globally attracting in Remark 5. Biologically, Theorem 3.3 implies that the heroin epidemic will remain spread and persist at the level of unique HSE in all groups, independent of the initial values if 0 > 1.
4. Numerical simulation. This section is devote to numerical simulation to support our analysis. We focus on the effects of two delays on the level of heroin spread. We consider the case where n = 2. Note that the expression of 02 (described in Remark 3) takes the following form
Taking τ 2 = 1,ĝ(s) = 1. According to the effects of τ 1 and τ 2 on σ and δ, we simulate that 02 increases as the values of σ and δ increase (see Fig. 1). Furthermore, Fig. 2 shows that τ 1 has the positive effects on the values of σ and 02 . Taking τ 1 = 1 andf (s) = 1. The other parameter values are the same as in Fig.  1. Simulation result shows that τ 2 has the positive effects on the values of δ and 02 (see Fig. 3). Thus, the increasing delays will intense the spread heroin and raise the ultimately spread level. 5. Conclusions and discussions. In this paper, we have analyzed a delayed multi-group heroin epidemic model to describe the spread of heroin in heterogeneous populations. We provide a rigorous analysis of the model and establish its global dynamics. The threshold parameter, which corresponds to the well-known basic reproduction number 0 , completely determines the global dynamics of the model. That is, if 0 ≤ 1, then the HFE is globally attracting; while if 0 > 1, then the HSE is globally attracting by combining a proper Lyapunov functional and Krichhoff's matrix tree theorem.  Figure 2. Influence of τ 1 on σ ( Fig. 2(a)) and 02 (Fig. 2(b)).  Figure 3. Influence of τ 1 on δ (Fig. 3(a)) and 02 (Fig. 3(b)). Theorem 3.1 and Theorem 3.3 show that the model has the sharp threshold property determined by the basic reproduction number 0 . Thus, the heroin either is extinct or remains spread completely depending on the value of the basic reproduction number. Biologically, Theorem 3.2 implies that the heroin will persist in all groups of the population and will eventually settle at a constant level in each group. Since the heterogeneity of populations plays an important role on 0 , it affects the global dynamics of the model. However, it does not cause any sustained oscillations.