Mathematical analysis of an age structured heroin-cocaine epidemic model

This paper is devoted to studying the dynamics of a certain age structured heroin-cocaine epidemic model. More precisely, this model takes into account the following unknown variables: susceptible individuals, heroin users, cocaine users and recovered individuals. Each one of these classes can change or interact with others. In this paper, firstly, we give some results on the existence, uniqueness and positivity of solutions. Next, we obtain a threshold value \begin{document}$ r(\Psi'[0]) $\end{document} such that an endemic equilibrium exists if \begin{document}$ r(\Psi'[0]) > 1 $\end{document} . We then show that if \begin{document}$ r(\Psi'[0]) , then the disease-free equilibrium is globally asymptotically stable, whereas if \begin{document}$ r(\Psi'[0]) > 1 $\end{document} , then the system is uniformly persistent. Moreover, for \begin{document}$ r(\Psi'[0]) > 1 $\end{document} , we show that the endemic equilibrium is globally asymptotically stable under an additional assumption that epidemic parameters for heroin users and cocaine users are same. Finally, some numerical simulations are presented to illustrate our theoretical results.

1. Introduction. The mathematical and statistical techniques applied to the modeling of infectious diseases had an extensive development in last decades. In 1927, Kermack and McKendrick [15] applied the mathematical modeling to study the dynamics of human infectious disease transmission. Their model, as well as several modifications of it, has been widely investigated in [1,5,6,14,27] and references therein by introducing the age dependence and in [2, 8, 12, 17, 29-31, 33, 35] by considering the space and diffusion effects.
Recently, several mathematical models have been developed to describe the drug epidemics. In [4], Burattini et al. constructed an ODEs (ordinary differential equations) system to study the impact of the crack-cocaine use on the final prevalence of HIV/AIDS in a community of drug users. In [32], White and Comiskey constructed an ODEs system to consider the heroin epidemics. They obtained the basic reproduction number R 0 for their model and performed the sensitivity analysis of it. Moreover, they showed that the backward bifurcation of the endemic equilibrium can occur at R 0 = 1. In [24], Mulone and Straughan showed the stability of the equilibria in the White-Comiskey model. In [25], Samanta modified the White-Comiskey model into a nonautonomous model with distributed time delay and obtained sufficient conditions for the permanence and the extinction of the epidemic. In [19], Liu and Zhang studied the global behavior of a heroin epidemic model with distributed time delay: the drug-free equilibrium is globally asymptotically stable if the basic reproduction number R 0 is less than 1, whereas a disease endemic equilibrium is locally asymptotically stable and the system is uniformly persistent if R 0 > 1. For other studies on the drug epidemic models with time delay, see, e.g., [9,13,20,21].
The above models are systems of ODEs or RDEs (retarded differential equations). Recently, some authors have focused on the effect of the age that implies the time elapsed since the entering into a particular class. In [10] and [11], Fang et al. studied the global stability of heroin epidemic models with age-dependent susceptibility and treat-age, respectively. In [34], Yang et al. studied a heroin epidemic model with treat-age and nonlinear incidence rate. In [7], Djilali et al. considered a heroin epidemic model with treat-age and general nonlinear incidence rate. In [18], Liu and Liu studied the global behavior of a heroin epidemic model with age-dependent susceptibility and treat-age with relapse.
Most of the above studies focused on the epidemics of just one kind of drugs such as heroin. In this paper, we formulate a new mathematical model describing the dynamics of heroine and cocaine epidemics. The model enables us to consider the effect of the interaction between the epidemics of these two drugs. In our model, a susceptible to the use of drug is defined as an individual who is exposed to both of heroin and cocaine. The affected vector is defined as an individual using either one of these two drugs. We consider the dynamics in each class given by a system of age structured partial differential equations where the age structure means here the time spent by consuming heroin and cocaine. We obtain a threshold value allowing to indicate the global stability of each equilibrium: the disease-free equilibrium and the endemic equilibrium. More precisely, each situation correspond either for the eradication of the consummation of drug from the population or the persistence and spreading of it in the population. It is therefore of interest to explore the problem mathematically.
In this paper, we denote by S the density of individuals in the susceptible class, i 1 the density of individuals using heroin, i 2 the density of individuals using cocaine and R the density of recovered individuals. We assume that all newborns, given and noted by A per unit time, are susceptible. Let µ be the per capita natural death rate per unit time. We model the contact or the rate of becoming a drug user by the law of mass action given by a transmission rate β. Let θ 1 and θ 2 be the per capita age-specific recovery rates from the consumption of heroin and cocaine, respectively. A part of them becomes new user of the other drug by a rate k 1 for the conversion from heroin to cocaine users and by a rate k 2 for the converse situation. The part (1 − k 1 ) of heroin and 1 − k 2 of cocaine users becomes totally cured. We design also by the rates δ 1 and δ 2 to describe the individuals who have relapsed into one of the two infected classes (see Table 1 and Figure 1). The model is given by,

Symbol
Description Density of susceptible individuals at time t i 1 (t, ξ 1 ) Density of heroin users at time t and age ξ 1 i 2 (t, ξ 2 ) Density of cocaine users at time t and age ξ 2 R(t) Density of recovered individuals at time t A Number of all newborns per unit time µ Natural death rate per capita and unit time β(ξ i ) (i = 1, 2) Transmission rate for drug users with age ξ i (i = 1, 2) Recovery rate from the consumption of heroin at age ξ 1 θ 2 (ξ 2 ) Recovery rate from the consumption of cocaine at age ξ 2 k 1 Rate at which an individual recovered from the consumption of heroin becomes a cocaine user k 2 Rate at which an individual recovered from the consumption of cocaine becomes a heroin user δ 1 Rate at which a recovered individual becomes a heroin user δ 2 Rate at which a recovered individual becomes a cocaine user Table 1. Description of each symbol in model (1). Figure 1. Transfer diagram for model (1).
where ξ 1 is the time spent consuming heroin and ξ 2 is the time spent consuming cocaine. The boundary and initial conditions are given by, for t > 0, and For the above system, we make the following assumptions, 4. β(·), θ 1 (·) and θ 2 (·) are Lipschitz continuous on R + , that is, there exist positive Lipschitz constants L β , L θ1 and L θ2 such that hold for all t ≥ 0 and h > 0.

There exist nonempty intervals
In this paper, we will be interested in analyzing the dynamics of model. In the next section 2, we formulate the main model as an abstract equation and show the existence and uniqueness of solutions. Section 3 is devoted to the study of equilibria. Our main result in this section shows that there exists a threshold value that guarantees the existence of an endemic equilibrium. In Section 4, we prove the existence of a compact attractor. In Section 5, we prove the global stability of the disease-free equilibrium when the threshold value is less than 1. In Section 6, we prove the uniform persistence of the system when the threshold value is greater than 1. In Section 7, in a special case that θ 1 (·) = θ 2 (·) and k 1 = k 2 , we prove the global stability of the endemic equilibrium when the threshold value is greater than 1. Finally, in Section 8, some numerical simulations are presented to illustrate our theoretical results.
2. Existence and uniqueness of the solution. To prove the existence and uniqueness of the solution of system (1), we use the approach in [22,28]. Let X := L 1 (R + , R) and let X + := L 1 (R + , R + ) be the positive cone of X. Let Y := R × X × X × R and let Y + := R + × X + × X + × R + be the positive cone of Y .
Lemma 2.1. A is a Hille-Yosida operator, that is, there exist real constants M ≥ 1 and ω ∈ R such that (ω, +∞) ⊂ ρ(A), and where ρ(·) denotes the resolvent set of an operator, and I d denotes the identity operator.
Using Lemmas 2.1 and 2.2, we prove the following proposition on the existence and uniqueness of a local integral solution of system (1). Proposition 1. For any C > 0, there exists a T C > 0 such that system (1) has the unique integral solution u ∈ C([0, Proof. By Lemmas 2.1 and 2.2, we can apply [22, Lemma 3.1 and Proposition 3.2] to obtain that system (1) has the unique integral solution u ∈ C([0, T C ], Y 0 ). To show its positivity, it is sufficient from [22,Proposition 3.6] to show the following two properties. (5). We can easily see that (ii) holds by choosing γ C := µ + 2β ∞ C. This completes the proof.
To prove the existence and uniqueness of the global classical solution of system (1), we need the following lemma.
Proof. For any y := (ϕ, 0, φ 1 (·), 0, φ 2 (·), ψ) T ∈ Y 0 , let us define the linear operator, We then have It is easy to see that F [y] is a bounded operator for any fixed y. Thus, F is continuously differentiable. This completes the proof.
Using Proposition 1 and Lemma 2.3, we establish the following theorem on the existence and uniqueness of the global classical solution of system (1).
Proof. LetC := u 0 W . Then, by Proposition 1 and Lemma 2.3, we can apply [22,Theorem 4.4] (see also [28,Theorem 3.7]) to obtain that there exists a TC > 0 such that system (1) has the unique local classical solution This completes the proof.

Equilibria.
It is easy to see that system (1) always has the disease-free equi- the endemic equilibrium of system (1). Then, the following equations hold.
which implies that Φ(·) is point dissipative. Let C ⊂ Y + 0 be an arbitrary bounded set such that for some K > 0, y W ≤ K holds for any y ∈ C. As in the proof of Theorem 2.4, we have, for any u 0 ∈ C, which implies that Φ(·) is eventually bounded on C. This completes the proof.
To show the existence of a compact attractor, we have to prove the asymptotic smoothness of semiflow Φ(·). Let J 1 (t) := i 1 (t, 0) and J 2 (t) := i 2 (t, 0) for all t ≥ 0. Now we have the following lemma.
Lemma 4.2. S(·), J 1 (·), J 2 (·) and R(·) are Lipschitz continuous on R + , that is, there exist positive Lipschitz constants L S , L J1 , L J2 and L R such that  (1), we have For any t ≥ 0 and h > 0, we have Note that we have, for j = 1, 2, and hence, Hence, we have In a similar way, we have Consequently, from (12), we have which implies that J 1 (·) is Lipschitz continuous on R + . In a similar way, we can show that J 2 (·) is Lipschitz continuous on R + . This completes the proof.
From Lemma 4.1 and Proposition 2, we can apply [26, Theorem 2.33] to establish the following proposition on the existence of a compact attractor. We end this section by describing the total trajectories of the system (1)-(2)-(3). Let U : R → Y + 0 such that U(t) = (S(t), 0, i 1 (t, .), 0, i 2 (t, .), R(t)) T be a total trajectory of semiflow Φ(·), that is, U(·) is a function satisfying U(r + t) = Φ(t)U(r) for all r ∈ R and t ≥ 0. For U(·), by a simple computation see also [26] and for all t ∈ R, we have and 5. Global stability of the disease-free equilibrium. From a well-known fact and Proposition 3, the compact attractor A 0 is consisting of total trajectories (see Proposition 2-34 [26]), that is, for any U 0 ∈ A 0 , there exists a total trajectory such that U(0) = U 0 and U(t) ∈ A 0 for all t ∈ R. Let J j := sup t∈R J j (t), j = 1, 2, where J j (t) = i j (t, 0) for all t ∈ R and j = 1, 2. We now prove the following lemma.
Proof. The nonnegativity of S(t) and R(t) for all t ∈ R is obvious since A 0 ⊂ Y + 0 . From the first equation of (1), we have dS(t)/dt ≤ A − µS(t), t ∈ R. This implies that S(t) ≤ A/µ, t ∈ R. Similarly, from the fourth equation of (1), we have This completes the proof.
By using Lemma 5.1, we prove the global asymptotic stability of the disease-free equilibrium for r(Ψ [0]) < 1.
Proof. First recall that the compact attractor of bounded set is the union of bounded total orbits. Let U(t) = (S(t), 0, i 1 (t, .), 0, i 2 (t, .), R(t)) T be a total trajectory of semiflow Φ(·). By (15)-(16) and Lemma 5.1, we have, for t ∈ R, Taking the supremum of the right-hand sides of these inequalities, we have by [26,Theorem 2.39]. This implies that the E 0 is globally asymptotically stable. This completes the proof.
6. Uniform persistence of the system. Let : Y + 0 → R be a function defined by We first prove the following proposition.
By adding the equations in (2), we have, for all t ∈ R, where, for all t ∈ R, Note thatÎ(0) = (U(0)). We next prove the following lemma.
Lemma 6.2. For total trajectory U(·) ∈ A 0 , we have, for all t ∈ R, Proof. By (11) in the proof of Lemma 4.1, we see that U(t) W ≤ A/µ for all t ∈ R. This implies that +∞ 0 i j (t, ξ j )dξ j ≤ A/µ for all t ∈ R and j = 1, 2. We then have from the first equation in (1) that, for all t ∈ R, For any t ∈ R and r < t, we then have Letting r → −∞, we obtain S(t) ≥ S L . This completes the proof.
By using Lemmas 6.1 and 6.2, we prove the following proposition.
Proposition 5. Assume that β(·) is not zero almost everywhere. Either (U(·)) is identically zero or strictly positive on R.
Proof. By (24) and Lemma 6.2, we have, for all t ∈ R, By Lemma 6.1 and appropriate shifts, for any r ∈ R, we see that if (U(t)) = 0 for all t ≤ r, then U(·) is identically zero. Suppose that there exists a decreasing sequence {r j } +∞ j=1 such that r 1 < r and r j → −∞ as j → +∞ and (U(r j )) > 0 for all j ∈ N. For any t ∈ R and j ∈ N, let U j (t) := U(t + r j ) andÎ j (t) :=Î(t + r j ). We then have from (25) that, for all t ∈ R, Note thatÎ j (·) is continuous at 0 andÎ j (0) =Î(r j ) = (U(r j )) > 0. By [26,Corollary B.6], we see that there exists a b > 0 such that (U j (t)) > 0 for all t > b. This implies that (U(t)) > 0 for all t > b + r j . Since r j → −∞ as j → +∞, we obtain that (U)(t) > 0 for all t ∈ R. This implies that (U(·)) is strictly positive on R. This completes the proof.
Proof. By Propositions 4 and 5, we can apply [26,Theorem 5.2] to conclude that the uniform weak -persistence implies the uniform strong -persistence. This completes the proof.
From [26,Theorem 5.7], we have the following result. Theorem 6.4. There exists a compact attractor A 1 that attracts all solutions with initial condition belonging to Ω 0 . Moreover there exists some ε > 0 such that, for any total trajectory U(·) in A 1 , A 1 is called the persistence attractor [23,Section 8]. Note that (26) states that the sum of i 1 (t, 0) and i 2 (t, 0) is bounded below by a positive constant. To show that both of i 1 (t, 0) and i 2 (t, 0) are bounded below by a positive constant is important for constructing a Lyapunov function. The main purpose of the following proposition is to prove that our solution is bounded away from 0 in the persistence attractor A 1 .
Proposition 6. For any total trajectory U(·) ∈ A 1 , the following estimates hold for all t ∈ R: Proof. The first inequality follows directly from Lemma 6.2. By (26) we have i 1 (t, 0) + i 2 (t, 0) ≥ ε for all t ∈ R. Without loss of generality, we suppose that i 1 (t, 0) ≥ ε for all t ∈ R. Hence from the equation of i 1 in (15), we get Since i * 1 (ξ 1 ) = Θ 1 (ξ 1 )i * 1 (0) for all ξ 1 ∈ R + , we have Similarly, in view of the equation of i 2 (t, 0) in (16) and combining with the above result, From the equation of R in (15) and combining with the above results, we have so by a simple computation we conclude that This completes the proof.
Lemma 7.1. The following equality holds, where, for simplicity, we write S = S(t) and R = R(t).
Let g(x) = x − 1 − ln x, x > 0. We now prove the following theorem.
Let V 3 (R) := R * g(R/R * ). Using the fourth equation in (28), the derivative of V 3 (R) is calculated as follows.
Let V (U) := V 1 (S) + V 2 (i) + CV 3 (R), where C > 0 is a constant defined below. By (32)- (34), the derivative of V (U) is calculated as follows. Let We then have from (35) that We then obtain from (37) that Thus, we obtain V (U) ≤ 0. This implies that V (U(t)) is nonincreasing for all t ∈ R.
Finally, from the third equality in (39), we conclude that M consists of only the endemic equilibrium. Since V attains its minimum at the endemic equilibrium and it is nonincreasing, U(t) = (S * , 0, i * (·), R * ) holds for all t ∈ R. Therefore, the compact attractor A 1 is reduced to the endemic equilibrium. In addition by [26,Theorem 2.39] the endemic equilibrium is also locally stable. Uniqueness of this one is a direct consequence of the fact that V (U(t)) = 0 holds only on the line S = S * . This completes the proof. For these parameters r(Ψ [0]) = 0.9414 < 1. According to the result stated, the disease-free equilibrium is globally asymptotically stable (Figures 3 and 4).  2(a − 12) 2 e −0.12(a−12) × 10 −3 , if a > 12. In this case the recovery rate θ 1 is more important than the recovery rate θ 2 (see Figure 2), then r(Ψ [0]) = 1.4163 > 1. From Figures 5 and 6, it is seen that the solution converge to the endemic equilibrium and the density of individuals i 1 is significantly lower than the density of individuals i 2 at equilibrium.