Epidemic dynamics on a delayed multi-group heroin epidemic model ...
Available online at www.tjnsa.com J. Nonlinear Sci. Appl. 9 (2016), 2149–2160 Research Article
Epidemic dynamics on a delayed multi-group heroin epidemic model with nonlinear incidence rate Xianning Liua , Jinliang Wangb,∗ a
Key Laboratory of Eco-environments in Three Gorges Reservoir Region (Ministry of Education), School of Mathematics and Statistics, Southwest University, Chongqing 400715, China. b
School of Mathematical Science, Heilongjiang University, Harbin 150080, China. Communicated by Y. Je Cho
Abstract For a multi-group Heroin epidemic model with nonlinear incidence rate and distributed delays, we study some aspects of its global dynamics. By a rigorous analysis of the model, we establish that the model demonstrates a sharp threshold property, completely determined by the values of 0. In this paper, we study some aspects of the global dynamics for the following system of differential equations, n-group heroin epidemic model related to (1.1), to describe the epidemic dynamics of heroin users: n X dS(t) = Λi − βij G(Si (t))U1j (t) − µi Si (t), dt j=1 Z τ n X dU1i (t) = βij G(Si (t))U1j (t) − pi U1i (t) + fi (s)pi U1i (t − s)e−(µi +δ2i )s ds − (µi + δ1i )U1i , (2.1) dt 0 j=1 Z τ dU (t) 2i = pi U1i (t) − (µi + δ2i )U2i − fi (s)pi U1i (t − s)e−(µi +δ2i )s ds, i = 1, 2, · · · , n. dt 0 The parameters are described in Table (1). The initial conditions for model (2.1) are: (Si (θ), U1i (θ), U2i (θ)) = φi , i = 1, 2, · · · , n,
(2.2)
where φi = (φsi (θ), φ1i (θ), φ2i (θ)) ∈ C([−τ, 0], R3+ ). C([−τ, 0], R3+ ) stands for the space of continuous functions mapping [−τ, 0] into R3+ , equipped with the supremum norm, kφi k = maxθ∈[−τ,0] |φi (θ)|. For the continuity of the solutions of (2.1), we set Z τZ 0 φ2i (0) = fi (s)pi U1i (η)e(µi +δ2i )η dηds. (2.3) 0
−s
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160
2152
Table 1: The descriptions of the parameters in the model (2.1).
Parameter Λi µi pi δ1i δ2i βij
Description constant recruitment rate in group i; death rate of individuals in group i; the rate of drug users who enter treatment in group i; the sum of rates of drug-related deaths of heroin users and recovery in group i; the sum of rates of drug-related deaths of heroin users undergoing treatment and immunity to drug addition in group i; the rate of becoming infectious heroin user into group i through contact between susceptible individual in group i and drug user in group j.
It follows from the third equation of (2.1) and initial condition (2.3) that Z τZ t fi (s)pi U1i (η)e−(µi +δ2i )(t−η) dηds. U2i (t) = 0
t−s
In order to make the mathematics tractable, we make the following assumptions: Assumption 2.1. Consider system (2.1), we make the following assumptions on parameters: (i) Λi , µi , pi , δ1i , δ2i and βij are positive for all i = 1, 2, . . . , n. (ii) βij is nonnegative, and n-square matrix (βij )1≤i,j≤n is irreducible [1]. Denote by B = (βij )1≤i,j≤n the contact matrix, which encodes the patterns of contact and transmission among groups that are built into the model. Based on the fact that U2i does not appear in other equations of (2.1), it suffices only to work on two out of the three variables and therefore, the dynamics are governed by the reduced system n X dSi (t) = Λi − βij G(Si (t))U1j (t) − µi Si (t), i = 1, 2, · · · , n, dt j=1 (2.4) Z τ n X dU1i (t) −(µ +δ )s i 2i = βij G(Si (t))U1j (t) − pi U1i (t) + fi (s)pi U1i (t − s)e ds − (µi + δ1i )U1i . dt 0 j=1
The initial conditions for model (2.4) are: Si (θ) = φsi (θ), U1i (θ) = φ1i (θ), i = 1, 2, · · · , n,
(2.5)
where φi = (φsi (θ), φ1i (θ)) ∈ C([−τ, 0], R2+ ), the space of continuous functions mapping [−τ, 0] into R2+ . In the rest of this paper, we shall investigate the dynamics of system (2.4) with initial conditions (2.5). 3. Positivity and boundedness of solutions The standard fundamental theory of FDEs [7] implies that there is a unique solution (Si (t), U1i (t)), U1i (t)), i = 1, 2, · · · , n to the system with given initial conditions (φsi (θ), φ1i (θ), φ2i (θ)), i = 1, 2, . . . , n. The following results address the well-posedness of system (2.1) with (2.2). Theorem 3.1. Let (Si (t), U1i (t), U2i (t)), i = 1, 2, · · · , n be a solution of system (2.1) satisfying (2.2). Then it remains non-negative and bounded for all t ≥ 0.
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160
2153
Proof. First, for i = 1, 2, · · · , n, we prove that Si (t) > 0 for all t ≥ 0. If there exists a t0 > 0, such that min Si (t0 ) = 0, and Si (t) > 0 for all t ∈ (0, t0 ). Without the loss of generality, we may assume that Si (t0 ) = 0. Thus, we obtain Si0 (t0 ) = Λi > 0, which leads to a contradiction. Hence, for i = 1, 2, · · · , n, Si (t) > 0, for all t ≥ 0. Using the variation of constants formula to the second equation of system (2.1), we obtain U1i (t) = e−(µi +pi +δ1i )t φ1i (θ) X Z t Z n −(µi +pi +δ1i )(t−ξ) + e βij G(Si (ξ))U1j (ξ) + 0
τ
−(µi +δ2i )s
fi (s)pi U1i (ξ − s)e
ds dξ.
0
j=1
Hence, U1i (t) ≥ 0 if φ1i (θ) > 0. Recall that Z τZ t U2i (t) = fi (s)pi U1i (η)e−(µi +δ2i )(t−η) dηds ≥ 0. 0
t−s
Next, we prove the boundedness of the solutions. From the first equation of (2.1), we have dSi (t) ≤ Λi − µi Si (t), dt it follows that for i = 1, 2, · · · , n, we have lim supt→∞ Si (t) ≤ Si0 = Λi /µi . Furthermore, adding all equations of (2.1), we have dSi (t) dU1i (t) dU2i (t) + + = Λi − µi (Si + U1i (t) + U2i (t)) − δ1i U1i (t) − δ2i U2i (t) dt dt dt ≤ Λi − µi (Si (t) + U1i + U2i (t)), i = 1, 2, · · · , n. Hence we have lim supt→∞ (Si (t) + U1i (t) + U2i (t)) ≤ Λi /µi . Thus, the compact feasible region 3n 0 Γ := (Si , U1i , U2i ) ∈ R : 0 ≤ Si (t) ≤ Si , 0 ≤ Si (t) + U1i (t) + U2i (t) ≤ Λi /µi , i = 1, 2, · · · , n is positively invariant for system (2.1). From the biological significance, we only need to consider (2.4) in the following region 2n 0 Γ0 := (Si , U1i , ) ∈ R : 0 ≤ Si (t) ≤ Si , 0 ≤ Si (t) + U1i (t) ≤ Λi /µi , i = 1, 2, · · · , n . That is, the well-posedness of the system (2.4) directly follows. This completes the proof.
4. Basic reproduction number and equilibria It is easy to see system (2.4) always admits an equilibrium which is labeled as P0 = (S10 , 0, · · · , Sn0 , 0), i where Si0 = Λ µi , i = 1, 2, · · · , n, and it is called the drug-free equilibrium (DFE). It can be verified that solutions of (2.4) with initial condition (S1 (0), U11 (0), · · · , Sn (0), U1n (0)) ∈ R2n + remain non-negative. Therefore, in what follows, we consider model (2.4) in R2n + . Let Z τ Qi := fi (s)e−(µi +δ2i ) ds. 0
(4.1)
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160 It can be verified that Qi ∈ (0, 1) for all i. For system (2.4), let β11 G(S10 ) β12 G(S10 ) · · · .. .. .. F = . . . 0 0 βn1 G(Sn ) βn2 G(Sn ) · · · then the next generation matrix is β
0 11 G(S1 )
µ1 +δ11 +p1 −p1 Q1
FV −1 =
.. .
0) βn1 G(Sn µ1 +δ11 +p1 −p1 Q1
β1n G(S10 ) .. . 0 βnn G(Sn )
2154
and V = diag(µi + δ1i + pi − pi Qi ),
β12 G(S10 ) µ2 +δ12 +p2 −p2 Q2
.. .
··· .. .
βn2 G(S10 ) µ2 +δ12 +p2 −p2 Q2
···
β1n G(S10 ) µn +δ1n +pn −pn Qn
.. .
0) βnn G(Sn µn +δ1n +pn −pn Qn
.
Thus the basic reproduction number of model (2.4) is defined by the spectral radius of the next generation matrix (see e.g. [4, 22]): 0 for i = 1, 2, · · · , n. know that system Bv Set n X L= vi LEE , i=1
then n ∗ U (t) X U (t) U U (t) dLEE 1j 1j 1i 1i vi ≤ vi β¯ij − − ln ∗ ∗ ∗ dt (2.4) U1j U1i U1i (t)U1j i,j=1 i=1 X n n ∗ U (t) X U1j (t) U1i (t) U1i 1j ¯ ¯ − − vi βij ln = vi βij ∗ ∗ ∗ U1j U1i U1i (t)U1j
n X
i,j=1
i,j=1
=: H1 − H2 , where
n X U1j (t) U1i (t) ¯ H = vi βij − ; ∗ ∗ 1 U1j U1i i,j=1 n ∗ U (t) X U1i 1j ¯ vi βij ln ∗ . H2 = U1i (t)U1j i,j=1
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160
2159
¯ = 0, we have Next, we prove that H1 ≡ 0 for all U11 , U12 , · · · , U1n > 0. From Bv n X
β¯ji vj =
j=1
n X
β¯ik vi .
k=1
∗ into above yields Putting β¯ji = βji G(Sj∗ )U1i n X
∗ βji G(Sj∗ )U1i vj
j=1
=
n X
∗ βik G(Si∗ )U1k vi , i = 1, 2, · · · , n.
k=1
This implies that n X
vi βij G(Si∗ )U1j (t) =
i,j=1
n X i=1
U1i
n X j=1
βji G(Sj∗ )vj =
n n X U1i X
U∗ i=1 1i k=1
∗ βik G(Si∗ )U1k vi =
n X i,j=1
∗ vi βij G(Si∗ )U1j
U1i ∗ , U1i
and thus H1 ≡ 0 for all U11 , U12 , · · · , U1n > 0. Similar to the arguments in Section 5 of [19], just replacing k = 1, Ek = U1i and Ej = U1j in the equations (5.9) and (5.10) of [19], we can easily obtain H2 ≡ 0 for all U11 , U12 , · · · , U1n > 0, from the P classic graph-theoretic method based on Kirchhoff’s matrix tree theorem. Therefore, the function L = ni=1 vi LEE as defined in the Theorem 3.1 of [12] is a Lyapunov function for system (2.4), namely, dL dt ≤ 0 for all (S1 , U11 , S2 , U12 , · · · , Sn , U1n ) ∈ Γ. One can only show that the largest ∗ invariant subset where { dL dt = 0} is the singleton {P } using the same argument as that in [6, 12, 19]. By ∗ LaSalle’s invariance principle, P is globally asymptotically stable in Γ. That is, limt→∞ (Si (t), U1i (t)) = ∗ , · · · , S ∗ , U ∗ ). This completes the proof of Theorem 6.2. P ∗ = (S1∗ , U11 n 1n 7. Conclusion In this paper, we considered a delayed multi-group heroin epidemic model with relapse phenomenon and nonlinear incidence rate. The main contributions of the paper are the proofs of global stability of equilibria. The distributed delays are introduced by the time needed to return to an untreated drug user, which are not constants but vary according to drug users’ different temporal, social, and physical contexts. Although including the nonlinear incidence rate combined with the relapse distributed delays into the multigroup model leads to the analysis of the resulting system becoming very complex, we are able to make a rigorous analysis of the model and establish a sharp threshold property. By using the method of constructing Lyapunov functionals based on graph-theoretical approach for coupled systems, sufficient conditions for the global stability of equilibria are given. Acknowledgements The authors would like to thank the editors and the referees for their helpful comments. J. Wang was supported by National Natural Science Foundation of China (Nos. 11401182, 11471089), Science and Technology Innovation Team in Higher Education Institutions of Heilongjiang Province (No. 2014TD005), X. Liu was supported by the National Natural Science Foundation of China (No. 11271303). References [1] A. Berman, R. J. Plemmons, Nonnegative matrices in the mathematical sciences, Academic Press, New York, (1979). 2.1 [2] N. P. Bhatia, G. P. Szeg¨ o, Dynamical Systems: Stability Theory and Applications, in: Lecture Notes in Mathematics, vol. 35, Springer, Berlin, (1967). 6 [3] C. Chin, Control of Communicable Diseases Manual, American Public Health Association, Washington, (1999).1 [4] 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 (1990), 365–382. 4
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160
2160
[5] H. Guo, M. Y. Li, Z. Shuai, Global stability of the endemic equilibrium of multigroup SIR epidemic models, Can. Appl. Math. Q., 14 (2006), 259–284. 1, 2 [6] H. Guo, M. Y. Li, Z. Shuai, A graph-theoretic approach to the method of global Lyapunov functions, Proc. Amer. Math. Soc., 136 (2008), 2793–2802. 1, 2, 6 [7] J. K. Hale, Theory of Functional Differential Equations, Springer, New York, (1997). 3 [8] G. Huang, A. Liu, A note on global stability for heroin epidemic model with distributed delay, Appl. Math. Lett., 26 (2013), 687–691. 1 [9] G. Huang, J. Wang, J. Zu, Global dynamics of multi-group dengue disease model with latency distributions, Math. Meth. Appl. Sci., 38 (2015), 2703–2718. 1 [10] A. Korobeinikov, Global properties of infectious disease models with nonlinear incidence, Bull. Math. Biol., 69 (2007), 1871–1886. 2 [11] J. P. Lasalle, The Stability of Dynamical Systems, in: Regional Conference Series in Applied Mathematics, Philadelphia, SIAM, (1976). 5 [12] M. Y. Li, Z. Shuai, Global-stability problem for coupled systems of differential equations on networks, J. Differential Equations, 248 (2010), 1–20. 1, 2, 6 [13] M. Y. Li, Z. Shuai, C. Wang, Global stability of multi-group epidemic models with distributed delays, J. Math. Anal. Appl., 361 (2010), 38–47. 1 [14] J. Liu, T. Zhang, Global behaviour of a heroin epidemic model with distributed delay, Appl. Math. Lett., 24 (2011), 1685–1692. 1, 1 [15] S. W. Martin, Livestock Disease Eradication: Evaluation of the Cooperative State Federal Bovine Tuberculosis Eradication Program, National Academy Press, Washington, (1994). 1 [16] C. C. McCluskey, Global stability for an SIR epidemic model with delay and nonlinear incidence, Nonlinear Anal. Real World Appl., 11 (2010), 3106–3109. 2 [17] C. C. McCluskey, Global stability of an SIR epidemic model with delay and general nonlinear incidence, Math. Biosci. Eng., 7 (2010), 837–850. 2 [18] G. Mulone, B. Straughan, A note on heroin epidemics, Math. Biosci., 218 (2009), 138–141. 1 [19] H. Shu, D. Fan, J. Wei, Global stability of multi-group SEIR epidemic models with distributed delays and nonlinear transmission, Nonlinear Anal. Real World Appl., 13 (2012), 1581–1592. 6 [20] H. Smith, P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition, Cambridge University Press, Cambridge, (1995). 6 [21] R. Sun, Global stability of the endemic equilibrium of multigroup SIR models with nonlinear incidence, Comput. Math. Appl., 60 (2010), 2286–2291. 1 [22] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. 4, 4 [23] K. E. VanLandingham, H. B. Marsteller, G. W. Ross, F. G. Hayden, Relapse of herpes simplex encephalitis after conventional acyclovir therapy, JAMA, 259 (1988), 1051–1053. 1 [24] J. Wang, X. Liu, J. Pang, D. Hou, Global dynamics of a multi-group epidemic model with general exposed distribution and relapse, Osaka. J. Math., 52 (2015), 117–138. 1 [25] J. Wang, J. Pang, X. Liu, Modeling diseases with relapse and nonlinear incidence of infection: a multi-group epidemic model, J. Biol. Dynam., 8 (2014), 99–116. 1 [26] X. Wang, J. Yang, X. Li, Dynamics of a heroin epidemic model with very population. Appl. Math., 2 (2011), 732–738. 1 [27] J. Wang, J. Zu, X. Liu, G. Huang, J. Zhang, Global dynamics of a multi-group epidemic model with general relapse distribution and nonlinear incidence rate, J. Biol. Sys., 20 (2012), 235–258. 1, 6 [28] E. White, C. Comiskey, Heroin epidemics, treatment and ODE modelling, Math. Biosci., 208 (2007), 312–324. 1
Epidemic dynamics on a delayed multi-group heroin epidemic model with nonlinear incidence rate Xianning Liua , Jinliang Wangb,∗ a
Key Laboratory of Eco-environments in Three Gorges Reservoir Region (Ministry of Education), School of Mathematics and Statistics, Southwest University, Chongqing 400715, China. b
School of Mathematical Science, Heilongjiang University, Harbin 150080, China. Communicated by Y. Je Cho
Abstract For a multi-group Heroin epidemic model with nonlinear incidence rate and distributed delays, we study some aspects of its global dynamics. By a rigorous analysis of the model, we establish that the model demonstrates a sharp threshold property, completely determined by the values of 0. In this paper, we study some aspects of the global dynamics for the following system of differential equations, n-group heroin epidemic model related to (1.1), to describe the epidemic dynamics of heroin users: n X dS(t) = Λi − βij G(Si (t))U1j (t) − µi Si (t), dt j=1 Z τ n X dU1i (t) = βij G(Si (t))U1j (t) − pi U1i (t) + fi (s)pi U1i (t − s)e−(µi +δ2i )s ds − (µi + δ1i )U1i , (2.1) dt 0 j=1 Z τ dU (t) 2i = pi U1i (t) − (µi + δ2i )U2i − fi (s)pi U1i (t − s)e−(µi +δ2i )s ds, i = 1, 2, · · · , n. dt 0 The parameters are described in Table (1). The initial conditions for model (2.1) are: (Si (θ), U1i (θ), U2i (θ)) = φi , i = 1, 2, · · · , n,
(2.2)
where φi = (φsi (θ), φ1i (θ), φ2i (θ)) ∈ C([−τ, 0], R3+ ). C([−τ, 0], R3+ ) stands for the space of continuous functions mapping [−τ, 0] into R3+ , equipped with the supremum norm, kφi k = maxθ∈[−τ,0] |φi (θ)|. For the continuity of the solutions of (2.1), we set Z τZ 0 φ2i (0) = fi (s)pi U1i (η)e(µi +δ2i )η dηds. (2.3) 0
−s
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160
2152
Table 1: The descriptions of the parameters in the model (2.1).
Parameter Λi µi pi δ1i δ2i βij
Description constant recruitment rate in group i; death rate of individuals in group i; the rate of drug users who enter treatment in group i; the sum of rates of drug-related deaths of heroin users and recovery in group i; the sum of rates of drug-related deaths of heroin users undergoing treatment and immunity to drug addition in group i; the rate of becoming infectious heroin user into group i through contact between susceptible individual in group i and drug user in group j.
It follows from the third equation of (2.1) and initial condition (2.3) that Z τZ t fi (s)pi U1i (η)e−(µi +δ2i )(t−η) dηds. U2i (t) = 0
t−s
In order to make the mathematics tractable, we make the following assumptions: Assumption 2.1. Consider system (2.1), we make the following assumptions on parameters: (i) Λi , µi , pi , δ1i , δ2i and βij are positive for all i = 1, 2, . . . , n. (ii) βij is nonnegative, and n-square matrix (βij )1≤i,j≤n is irreducible [1]. Denote by B = (βij )1≤i,j≤n the contact matrix, which encodes the patterns of contact and transmission among groups that are built into the model. Based on the fact that U2i does not appear in other equations of (2.1), it suffices only to work on two out of the three variables and therefore, the dynamics are governed by the reduced system n X dSi (t) = Λi − βij G(Si (t))U1j (t) − µi Si (t), i = 1, 2, · · · , n, dt j=1 (2.4) Z τ n X dU1i (t) −(µ +δ )s i 2i = βij G(Si (t))U1j (t) − pi U1i (t) + fi (s)pi U1i (t − s)e ds − (µi + δ1i )U1i . dt 0 j=1
The initial conditions for model (2.4) are: Si (θ) = φsi (θ), U1i (θ) = φ1i (θ), i = 1, 2, · · · , n,
(2.5)
where φi = (φsi (θ), φ1i (θ)) ∈ C([−τ, 0], R2+ ), the space of continuous functions mapping [−τ, 0] into R2+ . In the rest of this paper, we shall investigate the dynamics of system (2.4) with initial conditions (2.5). 3. Positivity and boundedness of solutions The standard fundamental theory of FDEs [7] implies that there is a unique solution (Si (t), U1i (t)), U1i (t)), i = 1, 2, · · · , n to the system with given initial conditions (φsi (θ), φ1i (θ), φ2i (θ)), i = 1, 2, . . . , n. The following results address the well-posedness of system (2.1) with (2.2). Theorem 3.1. Let (Si (t), U1i (t), U2i (t)), i = 1, 2, · · · , n be a solution of system (2.1) satisfying (2.2). Then it remains non-negative and bounded for all t ≥ 0.
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160
2153
Proof. First, for i = 1, 2, · · · , n, we prove that Si (t) > 0 for all t ≥ 0. If there exists a t0 > 0, such that min Si (t0 ) = 0, and Si (t) > 0 for all t ∈ (0, t0 ). Without the loss of generality, we may assume that Si (t0 ) = 0. Thus, we obtain Si0 (t0 ) = Λi > 0, which leads to a contradiction. Hence, for i = 1, 2, · · · , n, Si (t) > 0, for all t ≥ 0. Using the variation of constants formula to the second equation of system (2.1), we obtain U1i (t) = e−(µi +pi +δ1i )t φ1i (θ) X Z t Z n −(µi +pi +δ1i )(t−ξ) + e βij G(Si (ξ))U1j (ξ) + 0
τ
−(µi +δ2i )s
fi (s)pi U1i (ξ − s)e
ds dξ.
0
j=1
Hence, U1i (t) ≥ 0 if φ1i (θ) > 0. Recall that Z τZ t U2i (t) = fi (s)pi U1i (η)e−(µi +δ2i )(t−η) dηds ≥ 0. 0
t−s
Next, we prove the boundedness of the solutions. From the first equation of (2.1), we have dSi (t) ≤ Λi − µi Si (t), dt it follows that for i = 1, 2, · · · , n, we have lim supt→∞ Si (t) ≤ Si0 = Λi /µi . Furthermore, adding all equations of (2.1), we have dSi (t) dU1i (t) dU2i (t) + + = Λi − µi (Si + U1i (t) + U2i (t)) − δ1i U1i (t) − δ2i U2i (t) dt dt dt ≤ Λi − µi (Si (t) + U1i + U2i (t)), i = 1, 2, · · · , n. Hence we have lim supt→∞ (Si (t) + U1i (t) + U2i (t)) ≤ Λi /µi . Thus, the compact feasible region 3n 0 Γ := (Si , U1i , U2i ) ∈ R : 0 ≤ Si (t) ≤ Si , 0 ≤ Si (t) + U1i (t) + U2i (t) ≤ Λi /µi , i = 1, 2, · · · , n is positively invariant for system (2.1). From the biological significance, we only need to consider (2.4) in the following region 2n 0 Γ0 := (Si , U1i , ) ∈ R : 0 ≤ Si (t) ≤ Si , 0 ≤ Si (t) + U1i (t) ≤ Λi /µi , i = 1, 2, · · · , n . That is, the well-posedness of the system (2.4) directly follows. This completes the proof.
4. Basic reproduction number and equilibria It is easy to see system (2.4) always admits an equilibrium which is labeled as P0 = (S10 , 0, · · · , Sn0 , 0), i where Si0 = Λ µi , i = 1, 2, · · · , n, and it is called the drug-free equilibrium (DFE). It can be verified that solutions of (2.4) with initial condition (S1 (0), U11 (0), · · · , Sn (0), U1n (0)) ∈ R2n + remain non-negative. Therefore, in what follows, we consider model (2.4) in R2n + . Let Z τ Qi := fi (s)e−(µi +δ2i ) ds. 0
(4.1)
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160 It can be verified that Qi ∈ (0, 1) for all i. For system (2.4), let β11 G(S10 ) β12 G(S10 ) · · · .. .. .. F = . . . 0 0 βn1 G(Sn ) βn2 G(Sn ) · · · then the next generation matrix is β
0 11 G(S1 )
µ1 +δ11 +p1 −p1 Q1
FV −1 =
.. .
0) βn1 G(Sn µ1 +δ11 +p1 −p1 Q1
β1n G(S10 ) .. . 0 βnn G(Sn )
2154
and V = diag(µi + δ1i + pi − pi Qi ),
β12 G(S10 ) µ2 +δ12 +p2 −p2 Q2
.. .
··· .. .
βn2 G(S10 ) µ2 +δ12 +p2 −p2 Q2
···
β1n G(S10 ) µn +δ1n +pn −pn Qn
.. .
0) βnn G(Sn µn +δ1n +pn −pn Qn
.
Thus the basic reproduction number of model (2.4) is defined by the spectral radius of the next generation matrix (see e.g. [4, 22]): 0 for i = 1, 2, · · · , n. know that system Bv Set n X L= vi LEE , i=1
then n ∗ U (t) X U (t) U U (t) dLEE 1j 1j 1i 1i vi ≤ vi β¯ij − − ln ∗ ∗ ∗ dt (2.4) U1j U1i U1i (t)U1j i,j=1 i=1 X n n ∗ U (t) X U1j (t) U1i (t) U1i 1j ¯ ¯ − − vi βij ln = vi βij ∗ ∗ ∗ U1j U1i U1i (t)U1j
n X
i,j=1
i,j=1
=: H1 − H2 , where
n X U1j (t) U1i (t) ¯ H = vi βij − ; ∗ ∗ 1 U1j U1i i,j=1 n ∗ U (t) X U1i 1j ¯ vi βij ln ∗ . H2 = U1i (t)U1j i,j=1
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160
2159
¯ = 0, we have Next, we prove that H1 ≡ 0 for all U11 , U12 , · · · , U1n > 0. From Bv n X
β¯ji vj =
j=1
n X
β¯ik vi .
k=1
∗ into above yields Putting β¯ji = βji G(Sj∗ )U1i n X
∗ βji G(Sj∗ )U1i vj
j=1
=
n X
∗ βik G(Si∗ )U1k vi , i = 1, 2, · · · , n.
k=1
This implies that n X
vi βij G(Si∗ )U1j (t) =
i,j=1
n X i=1
U1i
n X j=1
βji G(Sj∗ )vj =
n n X U1i X
U∗ i=1 1i k=1
∗ βik G(Si∗ )U1k vi =
n X i,j=1
∗ vi βij G(Si∗ )U1j
U1i ∗ , U1i
and thus H1 ≡ 0 for all U11 , U12 , · · · , U1n > 0. Similar to the arguments in Section 5 of [19], just replacing k = 1, Ek = U1i and Ej = U1j in the equations (5.9) and (5.10) of [19], we can easily obtain H2 ≡ 0 for all U11 , U12 , · · · , U1n > 0, from the P classic graph-theoretic method based on Kirchhoff’s matrix tree theorem. Therefore, the function L = ni=1 vi LEE as defined in the Theorem 3.1 of [12] is a Lyapunov function for system (2.4), namely, dL dt ≤ 0 for all (S1 , U11 , S2 , U12 , · · · , Sn , U1n ) ∈ Γ. One can only show that the largest ∗ invariant subset where { dL dt = 0} is the singleton {P } using the same argument as that in [6, 12, 19]. By ∗ LaSalle’s invariance principle, P is globally asymptotically stable in Γ. That is, limt→∞ (Si (t), U1i (t)) = ∗ , · · · , S ∗ , U ∗ ). This completes the proof of Theorem 6.2. P ∗ = (S1∗ , U11 n 1n 7. Conclusion In this paper, we considered a delayed multi-group heroin epidemic model with relapse phenomenon and nonlinear incidence rate. The main contributions of the paper are the proofs of global stability of equilibria. The distributed delays are introduced by the time needed to return to an untreated drug user, which are not constants but vary according to drug users’ different temporal, social, and physical contexts. Although including the nonlinear incidence rate combined with the relapse distributed delays into the multigroup model leads to the analysis of the resulting system becoming very complex, we are able to make a rigorous analysis of the model and establish a sharp threshold property. By using the method of constructing Lyapunov functionals based on graph-theoretical approach for coupled systems, sufficient conditions for the global stability of equilibria are given. Acknowledgements The authors would like to thank the editors and the referees for their helpful comments. J. Wang was supported by National Natural Science Foundation of China (Nos. 11401182, 11471089), Science and Technology Innovation Team in Higher Education Institutions of Heilongjiang Province (No. 2014TD005), X. Liu was supported by the National Natural Science Foundation of China (No. 11271303). References [1] A. Berman, R. J. Plemmons, Nonnegative matrices in the mathematical sciences, Academic Press, New York, (1979). 2.1 [2] N. P. Bhatia, G. P. Szeg¨ o, Dynamical Systems: Stability Theory and Applications, in: Lecture Notes in Mathematics, vol. 35, Springer, Berlin, (1967). 6 [3] C. Chin, Control of Communicable Diseases Manual, American Public Health Association, Washington, (1999).1 [4] 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 (1990), 365–382. 4
X. Liu, J. Wang, J. Nonlinear Sci. Appl. 9 (2016), 2149–2160
2160
[5] H. Guo, M. Y. Li, Z. Shuai, Global stability of the endemic equilibrium of multigroup SIR epidemic models, Can. Appl. Math. Q., 14 (2006), 259–284. 1, 2 [6] H. Guo, M. Y. Li, Z. Shuai, A graph-theoretic approach to the method of global Lyapunov functions, Proc. Amer. Math. Soc., 136 (2008), 2793–2802. 1, 2, 6 [7] J. K. Hale, Theory of Functional Differential Equations, Springer, New York, (1997). 3 [8] G. Huang, A. Liu, A note on global stability for heroin epidemic model with distributed delay, Appl. Math. Lett., 26 (2013), 687–691. 1 [9] G. Huang, J. Wang, J. Zu, Global dynamics of multi-group dengue disease model with latency distributions, Math. Meth. Appl. Sci., 38 (2015), 2703–2718. 1 [10] A. Korobeinikov, Global properties of infectious disease models with nonlinear incidence, Bull. Math. Biol., 69 (2007), 1871–1886. 2 [11] J. P. Lasalle, The Stability of Dynamical Systems, in: Regional Conference Series in Applied Mathematics, Philadelphia, SIAM, (1976). 5 [12] M. Y. Li, Z. Shuai, Global-stability problem for coupled systems of differential equations on networks, J. Differential Equations, 248 (2010), 1–20. 1, 2, 6 [13] M. Y. Li, Z. Shuai, C. Wang, Global stability of multi-group epidemic models with distributed delays, J. Math. Anal. Appl., 361 (2010), 38–47. 1 [14] J. Liu, T. Zhang, Global behaviour of a heroin epidemic model with distributed delay, Appl. Math. Lett., 24 (2011), 1685–1692. 1, 1 [15] S. W. Martin, Livestock Disease Eradication: Evaluation of the Cooperative State Federal Bovine Tuberculosis Eradication Program, National Academy Press, Washington, (1994). 1 [16] C. C. McCluskey, Global stability for an SIR epidemic model with delay and nonlinear incidence, Nonlinear Anal. Real World Appl., 11 (2010), 3106–3109. 2 [17] C. C. McCluskey, Global stability of an SIR epidemic model with delay and general nonlinear incidence, Math. Biosci. Eng., 7 (2010), 837–850. 2 [18] G. Mulone, B. Straughan, A note on heroin epidemics, Math. Biosci., 218 (2009), 138–141. 1 [19] H. Shu, D. Fan, J. Wei, Global stability of multi-group SEIR epidemic models with distributed delays and nonlinear transmission, Nonlinear Anal. Real World Appl., 13 (2012), 1581–1592. 6 [20] H. Smith, P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition, Cambridge University Press, Cambridge, (1995). 6 [21] R. Sun, Global stability of the endemic equilibrium of multigroup SIR models with nonlinear incidence, Comput. Math. Appl., 60 (2010), 2286–2291. 1 [22] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. 4, 4 [23] K. E. VanLandingham, H. B. Marsteller, G. W. Ross, F. G. Hayden, Relapse of herpes simplex encephalitis after conventional acyclovir therapy, JAMA, 259 (1988), 1051–1053. 1 [24] J. Wang, X. Liu, J. Pang, D. Hou, Global dynamics of a multi-group epidemic model with general exposed distribution and relapse, Osaka. J. Math., 52 (2015), 117–138. 1 [25] J. Wang, J. Pang, X. Liu, Modeling diseases with relapse and nonlinear incidence of infection: a multi-group epidemic model, J. Biol. Dynam., 8 (2014), 99–116. 1 [26] X. Wang, J. Yang, X. Li, Dynamics of a heroin epidemic model with very population. Appl. Math., 2 (2011), 732–738. 1 [27] J. Wang, J. Zu, X. Liu, G. Huang, J. Zhang, Global dynamics of a multi-group epidemic model with general relapse distribution and nonlinear incidence rate, J. Biol. Sys., 20 (2012), 235–258. 1, 6 [28] E. White, C. Comiskey, Heroin epidemics, treatment and ODE modelling, Math. Biosci., 208 (2007), 312–324. 1
