EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 16, No. 1, 2023, 207-232 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Mathematical Framework for Modeling the Spread of HIV-disease within Two Different Age Classes Mahamadou Alassane1,∗, Abdoulaye Samaké1, Amadou Mahamane1, Ouaténi Diallo1 1 Université des Sciences, des Techniques et des Technologies de Bamako, BPE 423, Bamako, Mali Abstract. In this paper, we propose a mathematical model for the spread of HIV disease within two different age classes. We define a basic reproduction number R0 that depends on the charac- teristics of the two age classes. We prove that if R0 < 1, then the disease is extinct in both age classes. In contrast, we prove that if R0 > 1, then the disease is endemic in both age classes. 2020 Mathematics Subject Classifications: 34A12,37N25,65L05,92B05 Key Words and Phrases: HIV-AIDS, epidemic model, basic reproduction number, asymptotic stability, global stability, numerical simulations 1. Introduction Acquired Immunodeficiency Syndrome (AIDS) is one the most deadly disease caused by a Human Immunodeficiency Virus (HIV). The virus destroys all the immune system, in particular the CD4 T-lymphocytes, and leaves individuals susceptible to any other infections. It multiplies within those lymphocytes and eventually destroys them. Once the lymphocytes are depleted, then the immune system stops functioning properly. As a result, the individual can catch any kind of disease that might kill him easily because of the failure of its immune system. However, there are drugs that can slow down the progression of the virus. HIV-AIDS is usually transmitted in three different ways, namely the sexual contact, blood transfusion, mother-to-child exchanges during pregnancy, childbirth and breastfeeding. Many mathematical models are used to study the impact of preventive control strategies on the spread of HIV-AIDS in given populations [1, 3, 5, 6, 9–14]. Some of these models have shown that a change in risky behavior is necessary to prevent the spread of HIV-AIDS, even in the presence of a treatment [2, 8]. In this paper, we study the spread of HIV-AIDS in the age structured populations. In fact, we consider two different age classes: A first class that corresponds to individuals aged ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v16i1.4550 Email addresses: alassanemaiga@yahoo.fr (M. Alassane), abdoulaye.samake@usttb.edu.ml (A. Samaké), moulaye.ahmad@gmail.com (A. Mahamane), ouateni@yahoo.fr (O. Diallo) https://www.ejpam.com 207 © 2023 EJPAM All rights reserved. M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 208 25 years or less and a second class that corresponds to individuals aged over 25 years. We suppose that each age class is composed of susceptible and infected individuals. According to the Center of Disease Control (CDC), there are three main classes of HIV-AIDS infected individuals based on CD4 T-lymphocyte counts, see Table 1. Table 1: Main classes of HIV-AIDS infected individuals Stages of infection CD4 T-lymphocytes/mm3 Stage 1 > 500 Stage 2 200 < CD4 T < 500 Stage 3 < 200 The first stage of infection occurs between two and six weeks after HIV infection. The infected individual begins to produce antibodies that are detectable by HIV tests. The individual is then called HIV-positive. The second stage of infection is characterized by a reduction in the number of viral particles in the blood, marking the beginning of the clinical latency phase of the infection. Finally, the third stage of infection is characterized by the presence of major infections. We aim to capture the spread of HIV-AIDS in a population divided into two different age classes by a system of ordinary differential equations. The paper is organized as follows. In Section 2, we formulate a mathematical model for HIV-AIDS. The basic properties of the model are given in Section 3. In Section 4, the disease-free equilibrium point (DFE) and the basic reproduction number R0 are calculated. In section 5, we prove the local extinction of the infected populations in both age classes when R0 < 1. In Section 6, the global extinction of the disease in both age classes is studied and followed by some concluding results. Sections 7 and 8 deal with the persistence of the disease in both classes when R0 > 1. In Section 9, some numerical results are presented. The main conclusions are recapped in Section 10. 2. Mathematical Model for HIV-AIDS In this section, we formulate a mathematical model for HIV-AIDS. We divide the total population N into two age classes. The first age class is denoted by C1 and the second by C2. In each age class, there is one compartment of susceptible individuals and three compartments of infected individuals. In the class C1, S1 represents the compartment of susceptible individuals and I11, I 2 1 and I31 are the compartments of infected individuals at stage 1, 2 and 3 of infection, respectively. Similarly, in the class C2, S2 is the compart- ment of susceptible individuals and I12, I 2 2 and I32 represent the compartments of infected individuals at stage 1, 2 and 3 of infection, respectively. The total population N can be expressed as the following sum: N = 2∑ j=1 Sj + 3∑ j=1 Ij1 + 3∑ j=1 Ij2. M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 209 The Spread of HIV-AIDS within two different age classes in the population is illustrated in Figure 1. S2 m I12 m I22 m I32 m S1 m I11 m I21 m I31 m β2 a3 a4 γ β1 a1 a2 c1 c2 c3 bN Figure 1: Flow diagram of HIV-AIDS transmission dynamics Using the above representation, we formulate the corresponding dynamical model as follows: M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 210  d S1 dt = bN− S1 N β1 3∑ j=1 Ij1 + β2 3∑ j=1 Ij2 − (m + γ)S1 d I11 dt = S1 N β1 3∑ j=1 Ij1 + β2 3∑ j=1 Ij2 − (m + a1 + c1)I 1 1 d I21 dt = a1I 1 1 − (m + a2 + c2)I 2 1 d I31 dt = a2I 2 1 − (m + c3)I 3 1 d S2 dt = γS1 − S2 N β1 3∑ j=1 Ij1 + β2 3∑ j=1 Ij2 −mS2 d I12 dt = S2 N β1 3∑ j=1 Ij1 + β2 3∑ j=1 Ij2 − (m + a3)I 1 2 + c1I 1 1 d I22 dt = a3I 1 2 − (m + a4)I 2 2 + c2I 2 1 d I32 dt = a4I 2 2 −mI32 + c3I 3 1 (1) The system (1) is completed with the following initial conditions: S1 ≥ 0, I11 ≥ 0, I21 ≥ 0, I31 ≥ 0, S2 ≥ 0, I12 ≥ 0, I22 ≥ 0, I32 ≥ 0. (2) Summing the equations of system (1), we obtain: d N dt = (b−m)N. The parameters of the model are reported in Table 2. They are all positive. The parameter b represents the rate at which young begin sexual activity. The pa- rameter m is the death rate. The parameters β1 and β2 are infection rates in C1 and C2, respectively. The rate at which susceptible individuals in C1 get older and reach the age of becoming susceptible individuals in C2 is given by γ. The parameters a1 and a2 are the rates at which infected individuals in the class C1 move from stage 1 to stage 2 and from stage 2 to stage 3, respectively, within this class. Similarly, the parameters a3 and M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 211 Table 2: Description of the model parameters Parameters Description b, m Recruitment rate, natural mortality rate γ, c1, c2, c3 Transfer rates from C1 to C2 a1, a2 Transfer rates within C1 a3, a4 Transfer rates within C2 β1, β2 Transmission rates a4 describe the rates at which infected individuals in the class C2 move from stage 1 to stage 2 and from stage 2 to stage 3, respectively, within this class. The parameters ck, k = 1, 2, 3, denote the rates at which infected individuals in the stage k of C1 get older and reach the age of becoming infected individuals in the stage k of C2. For the mathematical analysis of the model, we introduce the following scalings: s1 = S1 N , i11 = I11 N , i21 = I21 N , i31 = I11 N , s2 = S2 N , i12 = I12 N , i22 = I22 N , i32 = I32 N . (3) From the equation 3, it holds: 2∑ j=1 sj + 3∑ j=1 ij1 + 3∑ j=1 ij2 = 1. (4) M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 212 According to the new variables, the system (1) can be rewritten as follows: d s1 dt = b - s1 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ1s1 d i11 dt = s1 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ2i 1 1 d i21 dt = a1i 1 1 − κ3i 2 1 d i31 dt = a2i 2 1 − κ4i 3 1 d s2 dt = γs1 − s2 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − bs2 d i12 dt = s2 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ5i 1 2 + c1i 1 1 d i22 dt = a3i 1 2 − κ6i 2 2 + c2i 2 1 d i32 dt = a4i 2 2 − b i32 + c3i 3 1 (5) where κ1 = b+γ, κ2 = b+a1+c1, κ3 = b+a2+c2, κ4 = b+c3, κ5 = b+a3, κ6 = b+a4. After normalization of the initial data, we obtain: 2∑ j=1 sj(0) + 3∑ j=1 ij1(0) + 3∑ j=1 ij2(0) = 1 (6) and s1(0) ≥ 0 , i11(0) ≥ 0 , i21(0) ≥ 0 , i31(0) ≥ 0 , s2(0) ≥ 0 , i12(0) ≥ 0 , i22(0) ≥ 0 , i32(0) ≥ 0. (7) The variables of the model (5) are reported in Table 3. M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 213 Table 3: Variables for the re-scaled HIV-AIDS model Variables Description s1 Proportion of susceptible individuals in class C1 s2 Proportion of susceptible individuals in class C2 i11 Proportion of individuals at stage 1 of infection in C1 i21 Proportion of individuals at stage 2 of infection in C1 i31 Proportion of individuals at stage 3 of infection in C1 i12 Proportion of individuals at stage 1 of infection in C2 i22 Proportion of individuals at stage 2 of infection in C2 i32 Proportion of individuals at stage 3 of infection in C2 3. Basic Properties Theorem 1. The feasible region Γ defined by Γ = e = ( s1, i11, i21, i31, s2, i12, i22, i32 ) ∈ R8 + : 0 ≤ 2∑ j=1 sj + 3∑ j=1 ij1 + 3∑ j=1 ij2 ≤ 1  with the initial conditions s1(0) ≥ 0 , i11(0) ≥ 0 , i21(0) ≥ 0 , i31(0) ≥ 0 , s2(0) ≥ 0 , i12(0) ≥ 0 , i22(0) ≥ 0 , i32(0) ≥ 0 is a positively invariant set for the system (5). Proof. 1. Positivity of Solutions We show by absurd that for all t ≥ 0, e(t) ≥ 0. Suppose that for a time t′ > 0, we have e(t’) < 0. The function e being continuous, from the intermediate value theorem, there exists a time t1 ∈]0, t’[ such that e(t1) = 0. Consider the equations of system (5) and let: ξ1(t) = exp ∫ t 0 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 + κ1 dτ  , ξ3(t) = exp (∫ t 0 (κ2 − β1s1) dτ ) , ξ2(t) = exp (κ3t) , ξ4(t) = exp (κ4t) , ξ5(t) = exp (κ6t) , ξ7(t) = exp (bt) , ξ6(t) = exp (∫ t 0 (κ5 − β2s2) dτ ) , ξ8(t) = exp ∫ t 0 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 + b  dτ  . M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 214 By differentiating each of the expressions s1ξ1, i 1 1ξ3, i 2 1ξ2, i 3 1ξ4, s2ξ8, i 1 2ξ5, i 2 2ξ6 and i32ξ7 with respect to time t, we obtain: ds1ξ1 dt = bξ1 (8a) di11ξ3 dt = s1ξ3 β1 3∑ j=2 ij1 + β2 3∑ j=1 ij2  (8b) di21ξ2 dt = a1i 1 1ξ2 (8c) di31ξ4 dt = a2i 2 1ξ4 (8d) ds2ξ8 dt = γs1ξ8 (8e) di12ξ6 dt = ξ6 β1s2 3∑ j=1 ij1 + β2s2 3∑ j=2 ij2 + c1i 1 1  (8f) di22ξ5 dt = ( a3i 1 2 + c2i 2 1 ) ξ5 (8g) di32ξ7 dt = ( a4i 2 2 + c3i 3 1 ) ξ7 (8h) By integrating the equations (8a)-(8h) between 0 and t1, it holds: s1(t1) = 1 ξ1(t1) [ s1(0) + ∫ t1 0 bξ1dt ] > 0 (9a) i11(t1) = 1 ξ3(t1) i11(0) + ∫ t1 0 s1ξ3 β1 3∑ j=2 ij1 + β2 3∑ j=1 ij2 dt  > 0 (9b) i21(t1) = 1 ξ2(t1) [ i21(0) + ∫ t1 0 a1ξ2i 1 1dt ] > 0 (9c) i31(t1) = 1 ξ4(t1) [ i31(0) + ∫ t1 0 a2ξ4i 2 1dt ] > 0 (9d) s2(t1) = 1 ξ8(t1) [ s2(0) + ∫ t1 0 γs1ξ8dt ] > 0 (9e) i12(t1) = 1 ξ6(t1) i12(0) + ∫ t1 0 ξ6 β1s2 3∑ j=1 ij1 + β2s2 3∑ j=2 ij2 + c1i 1 1 dt  > 0 (9f) i22(t1) = 1 ξ5(t1) [ i22(0) + ∫ t1 0 ( a3i 1 2 + c2i 2 1 ) ξ5dt ] > 0 (9g) i32(t1) = 1 ξ7(t1) [ i32(0) + ∫ t1 0 ( a4i 2 2 + c3i 3 1 ) ξ7dt ] > 0 (9h) M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 215 From (9a)-(9h), it follows that e(t1) > 0. This is a contradiction according to the start- ing hypothesis. Then, ∀t ≥ 0, e(t) ≥ 0. Therefore, all solutions initiated in R8 + are positive. 2. Invariant Region Summing the equations of system (5), we obtain: d dt s1 + 3∑ j=1 ij1 + s2 + 3∑ j=1 ij2  = 0. (10) Integrating (10) using initial conditions, it holds: ∀t ≥ 0, s1(t) + 3∑ j=1 ij1(t) + s2t + 3∑ j=1 ij2(t) ≤ 1. This achieves the proof. Consequently, in Γ, the model (5) is epidemiologically and mathematically well-posed. Therefore, it is sufficient to study the dynamics of the model in Γ. M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 216 4. Disease-Free Equilibrium (DFE) E0 and Reproduction Number R0 The Disease-Free Equilibrium (DFE) E0 of the model (5) is determined by solving the following system:  d s1 dt = 0 d i11 dt = 0 d i21 dt = 0 d i31 dt = 0 d s2 dt = 0 d i12 dt = 0 d i22 dt = 0 d i32 dt = 0 (11) In the case of absence of disease, i.e. the population size is zero in all compartments except the susceptible compartment, the solution of (11) is given by: E0 ( b κ1 , 0, 0, 0, γ κ1 , 0, 0, 0 ) . We determine the basic reproduction number R0 using the next generation matrix method at Disease-Free Equilibrium [12]. According to this method, R0 is defined as the effective number of secondary infections caused by typical infected individual during his/her entire period of infectioussness [1, 4]. Let X be the vector of infected classes: X = ( i11, i21, i31, i12, i22, i32 )T . F = s1 ( β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 ) , 0, 0, s2 ( β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 ) , 0, 0 T M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 217 denotes the vector of terms corresponding to new infections. V = ( κ2i 1 1, κ3i 2 1 − a1i 1 1, κ4i 3 1 − a2i 2 1, κ5i 1 2 − c1i 1 1, κ6i 2 2 − c2i 2 1 − a3i 1 2, b i32 − c3i 3 1 − a4i 2 2 )T refers to the vector of terms corresponding to individuals entering a given compartment and individuals leaving. The partial derivatives of F and V with respect to i11, i 2 1, i 3 1, i 1 2, i 2 2 and i32 are given by the following matrices F and V: F =  β1 b κ1 β1 b κ1 β1 b κ1 β2 b κ1 β2 b κ1 β2 b κ1 0 0 0 0 0 0 0 0 0 0 0 0 β1 γ κ1 β1 γ κ1 β1 γ κ1 β2 γ κ1 β2 γ κ1 β2 γ κ1 0 0 0 0 0 0 0 0 0 0 0 0  , V =  κ2 0 0 0 0 0 −a1 κ3 0 0 0 0 0 −a2 κ4 0 0 0 −c1 0 0 κ5 0 0 0 −c2 0 −a3 κ6 0 0 0 −c3 0 −a4 b  The next-generation matrix is defined by: K = FV−1 = 1 M7  M1b M2b M3b M4b M5b M6b 0 0 0 0 0 0 0 0 0 0 0 0 M1γ M2γ M3γ M4γ M5γ M6γ 0 0 0 0 0 0 0 0 0 0 0 0  where M1 = β1bκ3κ4κ5κ6 + a1β1bκ4κ5κ6 + β1a1a2bκ5κ6 + β2c1bκ3κ4κ6 + β2c2a1bκ4κ5 +β2a3c1bκ3κ4 + β2c3a2a1κ5κ6 + β2a4a1c2κ4κ5 + β2a4a3c1κ3κ4 (12) M2 = β1bκ2κ4κ5κ6 + β1a2bκ2κ5κ6 + β2c2bκ2κ4κ5 + β2c3a2κ2κ5κ6 +β2a4c2κ2κ4κ5 (13) M3 = β1bκ2κ3κ5κ6 + β2c3κ2κ3κ5κ6 (14) M4 = β2bκ2κ3κ4κ6 + β2a3bκ2κ3κ4 + β2a3a4κ2κ3κ4 (15) M5 = β2bκ2κ3κ4κ5 + β2a4κ2κ3κ4κ5 (16) M6 = β2κ2κ3κ4κ5κ6 (17) M7 = bκ1κ2κ3κ4κ5κ6 (18) The basic reproduction number R0, computed from the spectral raduis of the next- generation matrix K, is given by: R0 = b M1 + γM4 M7 (19) where M1, M4 and M7 are explicitly defined by (12), (15) and (18). M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 218 5. Local Stability of Disease-Free Equilibrium (DFE) Theorem 2. If R0 < 1, then the disease-free equilibrium point E0 ( b κ1 , 0, 0, 0, γ κ1 , 0, 0, 0 ) is locally asymptotically stable. Proof. Consider the system (5). The jacobian of this system at E0 is denoted by J(E0) and its eigenvalues are the solutions of the following equations: r8X 8 + r7X 7 + r6X 6 + r5X 5 + r4X 4 + r3X 3 + r2X 2 + r1X+ r0 = 0. Let us define: f1 = −a1κ5a4c2κ4β2b− a1a2κ5κ6b 2β1 − γβ2κ6bκ4κ3κ2 − a3γa4κ4κ3κ2β2a1κ5κ6b 2κ4β1 −a3a4κ4κ3c1β2b f2 = −β2b 2κ6κ4κ3c1 − a1a2c3κ5κ6β2b− a3β2b 2κ4κ3c1 − κ5κ6b 2κ4κ3β1 − a3γβ2bκ4κ3κ2 −a1κ5b 2c2κ4β2 f3 = bκ1κ2κ3κ4κ5κ6 f4 = −a1κ5a4c2β2b− β2b 2κ6κ3c1 − a3β2b 2κ3c1 − κ5κ6b 2κ3β1 − a3γβ2bκ3κ2 − a1κ5b 2c2β2 f5 = −γβ2κ6bκ3κ2 − a3γa4κ3κ2β2 − a1κ5κ6b 2β1 − a3a4κ3c1β2b + κ1κ5κ6bκ3κ2 f6 = −a3γβ2bκ4κ3 − γβ2κ6bκ4κ3 − a3γa4κ4κ3β2 + κ1κ5κ6bκ4κ3 f7 = −β2b 2κ6κ4c1 − a3β2b 2κ4c1 − κ5κ6b 2κ4β1 − a3γβ2bκ4κ2 f8 = −γβ2κ6bκ4κ2 − a3γa4κ4κ2β2 − a3a4κ4c1β2b + κ1κ5κ6bκ4κ2 f9 = −a1a4c2κ4β2b− a1a2c3κ6β2b− κ6b 2κ4κ3β1 − a1b 2c2κ4β2 − a1a2κ6b 2β1 − a1κ6b 2κ4β1 +κ1κ6bκ4κ3κ2 f10 = −β2b 2κ4κ3c1 − a1a2c3κ5β2b− κ5b 2κ4κ3β1 − a1a2κ5b 2β1 − γβ2bκ4κ3κ2 − a1κ5b 2κ4β1 +κ1κ5bκ4κ3κ2 f11 = −a1a2κ5κ6β1b− a3β2bκ4κ3c1 − κ5κ6bκ4κ3β1 − a3γβ2κ4κ3κ2 − a1κ5c2κ4β2b f12 = −β2bκ6κ4κ3c1 − a1κ5κ6κ4β1b− γβ2κ6κ4κ3κ2 + κ1κ5κ6κ4κ3κ2 f13 = −β2b 2κ6c1 − a3β2b 2c1 − κ5κ6b 2β1 − a3γβ2bκ2 − γβ2κ6bκ2 − a3γa4κ2β2 − a3a4c1β2b +κ1κ5κ6bκ2 f14 = −a3γβ2bκ3 − γβ2κ6bκ3 − a3γa4κ3β2 + κ1κ5κ6bκ3 f15 = −β2bκ4κ3c1 − κ5bκ4κ3β1 − a1a2κ5β1b− γβ2κ4κ3κ2 − a1κ5κ4β1b + κ1κ5κ4κ3κ2 f16 = −β2bκ6κ4c1 − a3β2bκ4c1 − κ5κ6bκ4β1 − a3γβ2κ4κ2 − γβ2κ6κ4κ2 + κ1κ5κ6κ4κ2 f17 = −β2bκ6κ3c1 − a3β2bκ3c1 − κ5κ6bκ3β1 − a3γβ2κ3κ2 − a1κ5c2β2b− γβ2κ6κ3κ2 −a1κ5κ6β1b + κ1κ5κ6κ3κ2 M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 219 f18 = −a1a4c2β2b− κ6b 2κ3β1 − a1b 2c2β2 − a1κ6b 2β1 + κ1κ6bκ3κ2 f19 = −a3γβ2bκ4 − γβ2κ6bκ4 − a3γa4κ4β2 + κ1κ5κ6bκ4 f20 = −κ6b 2κ4β1 + κ1κ6bκ4κ2 f21 = −a3γβ2κ4κ3 − γβ2κ6κ4κ3 + κ1κ5κ6κ4κ3 f22 = −γβ2bκ4κ3 + κ1κ5bκ4κ3 f23 = −β2b 2κ3c1 − κ5b 2κ3β1 − γβ2bκ3κ2 − a1κ5b 2β1 + κ1κ6bκ4κ3 + κ1κ5bκ3κ2 f24 = −κ6bκ4κ3β1 − a1c2κ4β2b− a1a2κ6β1b− a1κ6κ4β1b + κ1κ6κ4κ3κ2 f25 = −b2κ4κ3β1 − a1a2b 2β1 − a1b 2κ4β1 − a1a2c3β2b + κ1bκ4κ3κ2 f26 = −κ5b 2κ4β1 − γβ2bκ4κ2 − β2b 2κ4c1 + κ1κ5bκ4κ2 f27 = −γβ2bκ4 + κ1κ5bκ4 − κ6b 2β1 + κ1κ6bκ2 − κ4κ3β1b− a1a2β1b− a1κ4β1b + κ1κ4κ3κ2 f28 = −b2κ3β1 − a1b 2β1 + κ1bκ3κ2 − κ6bκ3β1 − a1c2β2b− a1κ6β1b + κ1κ6κ3κ2 f29 = −κ6bκ4β1 + κ1κ6κ4κ2 − β2b 2c1 − κ5b 2β1 − γβ2bκ2 + κ1κ5bκ2 f30 = −β2bκ6c1 − a3β2bc1 − κ5κ6bβ1 − a3γβ2κ2 − γβ2κ6κ2 + κ1κ5κ6κ2 f31 = −a3γβ2κ3 − γβ2κ6κ3 + κ1κ5κ6κ3 − a3γβ2b− γβ2κ6b− a3γa4β2 + κ1κ5κ6b f32 = −γβ2κ4κ3 + κ1κ5κ4κ3 − β2bκ4c1 − κ5bκ4β1 − γβ2κ4κ2 + κ1κ5κ4κ2 f33 = −β2bκ3c1 − κ5bκ3β1 − γβ2κ3κ2 − a1κ5β1b + κ1κ5κ3κ2 − a3γβ2κ4 − γβ2κ6κ4 + κ1κ5κ6κ4 f34 = −γβ2bκ3 + κ1κ5bκ3 − b2κ4β1 + κ1bκ4κ2 f35 = −b2β1 + κ1bκ2 − κ4β1b + κ1κ4κ2 − γβ2b + κ1κ5b− a3γβ2 − γβ2κ6 + κ1κ5κ6 − γβ2κ4 +κ1κ5κ4 f36 = −κ6bβ1 + κ1κ6κ2 − κ3β1b− β1ba1 + κ1κ3κ2 − β2bc1 − κ5bβ1 − γβ2κ2 + κ1κ5κ2 − γβ2κ3 +κ1κ5κ3 f37 = κ1κ6 − β1b + κ1κ5 − γβ2 + κ1b + κ1κ4 + κ2κ1 + κ1κ3 and F1 = f1 + f2 + f3; F2 = f13 + f14 + f15 + f16 + f17 + f18 + f19 + f20 + f21 + f22 + f23 + f24 + f25 + f26 F3 = f4 + f5 + f6 + f7 + f8 + f9 + f10 + f11 + f12; F4 = f27 + f28 + f29 + f30 + f31 + f32 + f33 + f34 F5 = f35 + f36; F6 = f37. Then r8 = 1, r7 = 1 κ1 ( F6 + κ1b + κ21 ) , r6 = 1 κ1 [ F5 + F6 (κ1 + b) + κ21b ] , r5 = 1 κ1 [F4 + F5 (κ1 + b) + F6κ1b] , r4 = 1 κ1 [F2 + F4 (κ1 + b) + F5κ1b], r3 = 1 κ1 [F3 + F2 (κ1 + b) + F4κ1b] , r2 = 1 κ1 [F1 + F3 (κ1 + b) + F2κ1b] , r1 = 1 κ1 [F1 (κ1 + b) + F3κ1b] , r0 = F1b. M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 220 We now need to verify that all the coefficients r0, r1, r2, r3, r4, r5 and r6 are positive. For this, it is sufficient to prove that F1, F2, F3 , F4, F5 and F6 are positive. We can clearly see that all the Fj, j = 1 to 6, are positive when R0 < 1. Therefore, all the rj, j = 1 to 6, are positive when R0 < 1. Therefore, from Routh Hurwitz Criterion, E0 is locally asymptotically stable if R0 < 1. 6. Global Stability In Theorem 2, we have proved that the disease-free equilibrium point E0 is locally asymptotically stable if R0 < 1. We will now prove that, independently of the initial population size, if R0 < 1, then the disease will die out. Let us define: x1 = b κ1 − s1 and x2 = γ κ2 − s2. The system (5) becomes: M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 221  d x1 dt = ( b κ1 − x1 )β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ1x1 d i11 dt = ( b κ1 − x1 )β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ2i 1 1 d i21 dt = a1i 1 1 − κ3i 2 1 d i31 dt = a2i 2 1 − κ4i 3 1 d x2 dt = −γ ( b κ1 − x1 ) + ( γ κ1 − x2 )β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − bx2 d i12 dt = ( γ κ1 − x2 )β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ5i 1 2 + c1i 1 1 d i22 dt = a3i 1 2 − κ6i 2 2 + c2i 2 1 d i32 dt = a4i 2 2 − b i32 + c3i 3 1 (20) E0 ( b κ1 , 0, 0, 0, γ κ1 , 0, 0, 0 ) is globally asymptotically stable for system (5) if and only if E0 (0, 0, 0, 0, 0, 0, 0, 0) is globally asymptotically stable for system (20) . Theorem 3. If R0 < 1, then the disease-free equilibrium point E0 (0, 0, 0, 0, 0, 0, 0, 0) is globally asymp- totically stable for the system (20). Proof. Consider the following function V : Γ −→ R+ defined by: V = κ−1 2 (M7 − γM4) i 1 1 + b M2i 2 1 +M4i 1 2 +M5i 2 2 + (M3 +M6)  2∑ j=1 xj − 3∑ j=1 ij1 − 2∑ j=1 ij2  M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 222 If R0 < 1, then M7 − γM4, M2, M3, M4, M5 and M6 are positive. Consequently, the function V is positive and vanishes at the disease-free equilibrium. The derivative of this Lyapunov function V along the trajectories of the ordinary differential system is V̇ = κ−1 2 (M7 − κγM4) (bκ−1 1 − x1 )β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ2i 1 1 + bκ−1 2 M2 ( a1i 1 1 − κ3i 2 1 ) + bκ−1 2 M3 ( a2i 2 1 − κ4i 3 1 ) + bκ−1 2 M4 (γκ−1 1 − x2 )β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ5i 1 2 + c1i 1 1  + bκ−1 2 M5 ( a3ii 1 2 − κ6i 2 2 + c2i 2 1 ) + bκ−1 2 M6 ( a4i 2 2 − bi32 + c3i 3 1 ) . We can also write V̇ = κ−1 2 (M7 − γM4) b κ1 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − κ−1 2 (M7 − γM4) β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2  x1 − (M7 − γM4) i 1 1 + bκ−1 2 M2a1i 1 1 − bκ−1 2 M2κ3i 2 1 + bκ−1 2 M3a2i 2 1 − bκ−1 2 M3κ4i 3 1 − bκ−1 2 M4κ5i 1 2 +bκ−1 2 γM4κ1 −1 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 − bκ−1 2 M4 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2  x2 + bκ−1 2 M4c1i 1 1 +bκ−1 2 M5a3i 1 2 − bκ−1 2 M5κ6i 2 2 + bκ−1 2 M5c2i 2 1 + bκ−1 2 M6a4i 2 2 − bκ−1 2 M6bi 3 2 + bκ−1 2 M6γ3i 3 2. Following algebraic manipulations, it holds: V̇ = −κ−1 2 (M7 − γM4) β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2  x1 − bκ−1 2 M4 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2  x2 + ( −M7 + γM4 + bκ−1 2 M2a1 + bκ−1 2 M4c1 + β1bκ −1 1 κ−1 2 M7 ) i11 or V̇ = −κ−1 2 (M7 − γM4) β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2  x1−bκ−1 2 M4 (kh + kf ) x2−(M7 − γM4 − bM1) i 1 1. If R0 < 1, then M7 − κθM4 and M7 − κθM4 − bM1 are positive, consequently, V̇ is negative definite along the trajectories of the system (5). Therefore, the DFE E0 is globally asymptotically stable for the system (5) if R0 < 1. This ends the proof of Theorem 3. 7. Existence of Endemic Equilibrium In this section, we analyze the existence of non-trivial endemic equilibrium E∗ (s∗1, i11∗, i21∗, i31∗, s∗2, i12∗, i22∗, i32∗) of the system (5). M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 223 Theorem 4. If R0 > 1, then there exists an endemic equilibrium point E∗ for the system (5). Proof. Solving the equations of system (5) at equilibrium state, we obtain: s1 = κ2i 1 1 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 (21a) s2 = κ5i 1 2 − c1i 1 1 β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 (21b) i21 = a1 κ3 i11 (21c) i31 = a1a2 κ3κ4 i11 (21d) i22 = a3 κ6 i12 + a1c2 κ6κ3 i11 (21e) i32 = a3a4 bκ6 i12 + ( a4a1c2 bκ6κ3 + a1a2c3 bκ4κ3 ) i11 (21f) i12 = −i11 f38 + f39i 1 1 β2κ3κ4 ( −f40 + f41i 1 1 ) (21g) β1 3∑ j=1 ij1 + β2 3∑ j=1 ij2 = b− κ1s1 s1 (21h) where f38 = −β1κ6b 2κ4κ3 − β1κ6b 2a1κ4 − β1κ6b 2a1a2 − β2a1c2b 2κ4 − bβ2a1a4c2κ4 − bβ2a1a2c3κ6 +κ1κ2κ6κ3bκ4 f39 = κ2β1κ6bκ4κ3 + κ2β1κ6ba1κ4 + κ2β1κ6ba1a2 + κ2β2a1c2bκ4 + κ2β2a1a4c2κ4 + κ2β2a1a2c3κ6 f40 = b2κ6 + a3b 2 + ba3a4 f41 = κ2bκ6 + κ2a3b + κ2a3a4. From (21b)-(21g), let f44 ( i11 )2 + f43i 1 1 + f42 = 0 (22) where f42 = b2 (1− R0) bκ1κ2κ3κ4κ5κ6 f43 = κ2b ( 2b2κ6κ5β1a1κ4 + 2b2κ6κ5β1a1a2 + 2b2κ6κ5β1κ4κ3 + 2b2κ6β2κ3κ4c1 + 2b2κ5β2a1c2κ4 +2b2β2κ3κ4c1a3 + 2bβ2κ3κ4c1a3a4 + 2bκ6κ5β2a1a2c3 − bκ6c1β2κ3κ4κ1 − bc1β2κ3κ4κ1a3 M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 224 −bκ6κ5κ1β1κ4κ3 + 2bκ2β2κ3κ4γa3 − bκ6κ5κ1β1a1κ4 − bκ6κ5κ1β1a1a2 − bκ2κ6κ5κ1κ3κ4 +2bκ5β2a1a4c2κ4 − bκ5κ1β2a1c2κ4 + 2bκ2κ6β2κ3κ4γ − κ6κ5κ1β2a1a2c3 + 2κ2β2κ3κ4γa3a4 +κ2κ6κ5κ 2 1κ3κ4 − c1β2κ3κ4κ1a3a4 − κ5κ1β2a1a4c2κ4 ) f44 = κ22 ( − κ2β2κ3κ4γa3a4 − bκ2κ6β2κ3κ4γ − bκ2β2κ3κ4γa3 + bκ6c1β2κ3κ4κ1 + bc1β2κ3κ4κ1a3 +c1β2κ3κ4κ1a3a4 − b2κ6κ5β1κ4κ3 − b2κ6κ5β1a1a2 + bκ6κ5κ1β1a1κ4 + κ6κ5κ1β2a1a2c3 −b2κ6κ5β1a1κ4 + bκ6κ5κ1β1a1a2 + κ5κ1β2a1a4c2κ4 − b2β2κ3κ4c1a3 + bκ6κ5κ1β1κ4κ3 −bβ2κ3κ4c1a3a4 − b2κ6β2κ3κ4c1 − b2κ5β2a1c2κ4 − bκ6κ5β2a1a2c3 − bκ5β2a1a4c2κ4 +bκ5κ1β2a1c2κ4 ) . By replacing κ1, κ2, κ3, κ4, κ5 and κ6 by their expressions into f44, we get: f44 = (a1 + c1 + b)2 bγ (−β2 + β1) (b + a4) (b + a3) ( b2 + ba2 + bc2 +bc3 + a1b + a1c3 + a1a2 + c2c3 + a2c3 ) . If R0 > 1, we clearly see that f42 < 0. We also note that: (i) If β1 > β2, then f44 > 0. Consequently, the discriminant of (22) is positive and the product of the solutions is negative. So, there exists a positive solution i11 for (22). (ii) If β1 < β2, then f44 < 0 and f43 > 0. Consequently, the discriminant of (22) is positive, the sum and the product of the solutions of (22) are positive. Therefore, there exist two distinct solutions. (iii) If β1 = β2, then f44 = 0 and f43 > 0. So, (22) becomes: f43i 1 1 + f42 = 0. (23) Therefore, there exists a unique positive solution i11. In all the cases discussed above, there exists at least a positive solution i11 of (22). Let us denote by i11 ∗ this positive solution. We will now define and prove the positivity of i12 ∗ , i31 ∗ , i22 ∗ , i32 ∗ , s1 ∗ and s2 ∗ if R0 > 1. From equation (21g), we get: i12 ∗ = −i11 ∗ f38 + f39i 1 1 ∗ β2κ3κ4 ( −f40 + f41i 1 1 ∗) If R0 > 1, then f38 > 0. From (21a) and (21h), it follows that −f40 + f41i 1 1 ∗ < 0, hence i12 ∗ > 0. From equation (21c), it follows that i21 ∗ = a1 κ3 i11 ∗ , hence i21 ∗ > 0. M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 225 From equation (21d), it follows that, i31 ∗ = a1a2 κ3κ4 i11 ∗ , hence i31 ∗ > 0. From equation (21e), it follows that, i22 ∗ = a3 κ6 i12 ∗ + a1c2 κ6κ3 i11 ∗ , hence i22 ∗ > 0. From equation (21f), it follows that, i32 ∗ = a3a4 bκ6 i12 ∗ + ( a4a1c2 bκ6κ3 + a1a2c3 bκ4κ3 ) i11 ∗ , hence i32 ∗ > 0. From equation (21a), it follows that, s∗1 = κ2i 1 1 ∗ β1 3∑ j=1 ij1 ∗ + β2 3∑ j=1 ij2 ∗ , hence s∗1 > 0. From equation (21b), it follows that, s∗2 = κ5i 1 2 ∗ − c1i 1 1 ∗ β1 3∑ j=1 ij1 ∗ + β2 3∑ j=1 ij2 ∗ , hence s∗2 > 0. Therefore, if R0 > 1, there exists a positive solution E∗ (s∗1, i11∗, i21∗, i31∗, s∗2, i12∗, i22∗, i32∗) for the system (5). 8. Stability Analysis of Endemic Equilibrium Theorem 5. If R0 > 1, then the endemic equilibrium point E∗ ( s∗1, i 1 1 ∗ , i21 ∗ , i31 ∗ , s∗2, i 1 2 ∗ , i22 ∗ , i32 ∗ ) is locally asymptotically stable. Proof. Let us compute the jacobian of the system at the point E∗ (s∗1, i11∗, i21∗, i31∗, s∗2, i12∗, i22∗, i32∗). M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 226 J(E∗) =  −H∗ − κ1 −β1s ∗ 1 −β1s ∗ 1 −β1s ∗ 1 0 −β2s ∗ 1 −β2s ∗ 1 −β2s ∗ 1 H∗ β1s ∗ 1 − κ2 β1s ∗ 1 β1s ∗ 1 0 β2s ∗ 1 β2s ∗ 1 β2s ∗ 1 0 a1 −κ3 0 0 0 0 0 0 0 a2 −κ4 0 0 0 0 γ −β1s ∗ 2 −β1s ∗ 2 −β1s ∗ 2 −H∗ − b −β2s ∗ 2 −β2s ∗ 2 −β2s ∗ 2 0 β1s ∗ 2 + c1 β1s ∗ 2 H∗ β2s ∗ 2 − κ5 β2s ∗ 2 β2s ∗ 2 0 0 c2 0 0 a3 −κ6 0 0 0 0 c3 0 0 a4 −b  where H∗ = β1 3∑ j=1 ij1 ∗ + β2 3∑ j=1 ij2 ∗ . The characteristic equation of J(E∗) is given by: A8X 8 +A7X 7 +A6X 6 +A5X 5 +A4X 4 +A3X 2 +A2X 2 +A1X+A0 = 0 (24) where A0, A1, A2, A3, A4, A5, A6, A7, A8 are obtained as a result of a boring calculation. They are all positive. Applying the Routh-Hurwitz criterion [7], it follows that all eigenvalues of the characteristic equation (24) have negative real part if R0 > 1. Therefore, the endemic solution E∗ is locally asymptotically stable if R0 > 1. 9. Numerical Simulations and Discussion In this section, we perform numerical simulations to support the theoretical results from the mathematical analysis of model (5). In addition to the verification of the theoretical results, these numerical solutions are very important from a practical point of view. We first consider the case where R0 = 0.67893 < 1 using the parameter values reported in Table 4. Using different initial conditions, the dynamics of the susceptible and infected popu- lations of the model are plotted in Figures 2, 3, 4 and 5. In Figure 2, we can observe that the proportions of susceptible individuals in classes C1 and C2 are consistent, (s1 = 0.63343, s2 = 0.36657). In contrast, as shown in Figures 3, 4 and 5, the proportions of infected individuals in classes C1 and C2 decline to zero( i11 = 0, i21 = 0, i31 = 0, i12 = 0, i22 = 0, i32 = 0 ) , i.e. approach the disease-free equilibrium M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 227 Table 4: Parameter values Parameter Value Source Parameter Value Source b 0.0432 Estimated a2 0.01 Estimated m 0.0096 Estimated a3 0.01 Estimated γ 0.025 Estimated a4 0.01 Estimated β1 0.035 Estimated c1 0.02 Estimated β2 0.025 Estimated c2 0.02 Estimated a1 0.01 Estimated c3 0.02 Estimated 0 200 400 600 800 1000 1200 1400 1600 time (t) 0.2 0.3 0.4 0.5 0.6 0.7 0.8 P ro p o rt io n o f s u s c e p ti b le i n d iv id u a ls o f c la s s C 1 (a) Proportion of susceptible individuals in C1 0 200 400 600 800 1000 1200 1400 1600 time (t) 0.1 0.15 0.2 0.25 0.3 0.35 0.4 P ro p o rt io n o f s u s c e p ti b le i n d iv id u a ls o f c la s s C 2 (b) Proportion of susceptible individuals in C2 Figure 2: Time series plots of the proportions of susceptible individuals in classes C1 and C2 for R0 = 0.67893 < 1 using various initial conditions and parameter values reported in Table 4. 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.01 0.02 0.03 0.04 0.05 0.06 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 1 o f c la s s C 1 (a) Proportion of individuals infected at stage 1 in C1 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 1 o f c la s s C 2 (b) Proportion of individuals infected at stage 1 in C2 Figure 3: Time series plots of the proportions of infected individuals at stage 1 in classes C1 and C2 for R0 = 0.67893 < 1 using various initial conditions and parameter values reported in Table 4. (DFE). They show that DFE is locally asymptotically stable when R0 < 1. These numer- ical simulations support the result stated in Theorem 2 on the stability of DFE. Further using the parameter values given in Table 5, we consider the case when R0 = M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 228 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 2 o f c la s s C 1 (a) Proportion of individuals infected at stage 2 in C1 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.02 0.04 0.06 0.08 0.1 0.12 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 2 o f c la s s C 2 (b) Proportion of individuals infected at stage 2 in C2 Figure 4: Time series plots of the proportions of infected individuals at stage 2 in classes C1 and C2 for R0 = 0.67893 < 1 using various initial conditions and parameter values reported in Table 4. 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 3 o f c la s s C 1 (a) Proportion of individuals infected at stage 3 in C1 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 3 o f c la s s C 2 (b) Proportion of individuals infected at stage 3 in C2 Figure 5: Time series plots of the proportions of infected individuals at stage 3 in classes C1 and C2 for R0 = 0.67893 < 1 using various initial conditions and parameter values reported in Table 4. 1.408 > 1. Table 5: Parameter values Parameter Value Source Parameter Value Source b 0.0432 Estimated a2 0.01 Estimated m 0.0096 Estimated a3 0.01 Estimated γ 0.025 Estimated a4 0.01 Estimated β1 0.075 Estimated c1 0.02 Estimated β2 0.05 Estimated c2 0.02 Estimated a1 0.01 Estimated c3 0.02 Estimated Using different initial conditions, the dynamics of the susceptible and infected popu- M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 229 lations of the model are plotted in Figures 6, 7, 8 and 9. 0 200 400 600 800 1000 1200 1400 1600 time (t) 0.2 0.3 0.4 0.5 0.6 0.7 0.8 P ro p o rt io n o f s u s c e p ti b le i n d iv id u a ls o f c la s s C 1 (a) Proportion of susceptible individuals in C1 0 200 400 600 800 1000 1200 1400 1600 time (t) 0.1 0.15 0.2 0.25 0.3 0.35 0.4 P ro p o rt io n o f s u s c e p ti b le i n d iv id u a ls o f c la s s C 2 (b) Proportion of susceptible individuals in C2 Figure 6: Time series plot of the proportions of susceptible individuals in classes C1 and C2 for R0 = 1.408 > 1 using various initial conditions and parameter values reported in Table 5. 0 200 400 600 800 1000 1200 1400 1600 time (t) 0.02 0.04 0.06 0.08 0.1 0.12 0.14 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 1 o f c la s s C 1 (a) Proportion of individuals infected at stage 1 in C1 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 1 o f c la s s C 2 (b) Proportion of individuals infected at stage 1 in C2 Figure 7: Time series plot of the proportions of individuals infected at stage 1 in classes C1 and C2 for R0 = 1.408 > 1 using various initial conditions and parameter values reported in Table 5. As shown in Figures 6, 7, 8 and 9, the proportions of susceptible and infected in- dividuals in classes C1 and C2 are consistent, [( s1 ∗, i11 ∗ , i21 ∗ , i31 ∗ , s2 ∗, i12 ∗ , i22, i32 ∗ ) =( 0.49548, 0.12853, 0.017559, 0.0027783, 0.19918, 0.11941, 0.029047, 0.0080157 )] , i.e., the population tends to endemic equilibrium E∗ when R0 > 1. This indicates that, regard- less of initial conditions, the infected population eventually reaches endemic equilibrium over time and the disease-free equilibrium point becomes unstable when R0 > 1. These numerical simulations support our theoretical results. 10. Conclusions In this paper, we have developed a mathematical model for the spread of HIV disease within two different age classes. We proposed a basic reproduction number that depends M. Alassane et al. / Eur. J. Pure Appl. Math, 16 (1) (2023), 207-232 230 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 2 o f c la s s C 1 (a) Proportion of individuals infected at stage 2 in C1 0 200 400 600 800 1000 1200 1400 1600 time (t) 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 2 o f c la s s C 2 (b) Proportion of individuals infected at stage 2 in C2 Figure 8: Time series plot of the proportions of individuals infected at stage 2 in classes C1 and C2 for R0 = 1.408 > 1 using various initial conditions and parameter values reported in Table 5. 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 3 o f c la s s C 1 (a) Proportion of individuals infected at stage 3 in C1 0 200 400 600 800 1000 1200 1400 1600 time (t) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 P ro p o rt io n o f in d iv id u a ls i n fe c te d a t s ta g e 3 o f c la s s C 2 (b) Proportion of individuals infected at stage 3 in C2 Figure 9: Time series plot of the proportions of individuals infected at stage 3 in classes C1 and C2 for R0 = 1.408 > 1 using various initial conditions and parameter values reported in Table 5. on the characteristics of the two age classes. We have proved that if the Routh-Hurwitz cri- terion are satisfied, then the disease-free equilibrium (DFE) E0 is locally asymptotically stable. We constructed a Lyapunov function to prove that the disease-free equilibrium (DFE) E0 is globally stable when R0 < 1. For R0 > 1, we obtain from mathematical analysis a quadratic equation in i11. It has been proven that the existence of an endemic equilibrium depends on the existence of at least one real positive value for i11. The stability analysis of endemic equilibrium produces that if the Routh-Hurwitz criterion are satisfied, then the endemic equilibrium E∗ is locally asymptotically stable. The important math- ematical results in this paper were all corroborated by numerical simulations performed using MATLAB. Indeed, we verified through numerical experiments that the disease-free equilibrium E0 is stable when R0 < 1. On the other hand, we numerically verified that, if R0 > 1, then the endemic equilibrium E∗ becomes stable. REFERENCES 231 References [1] Neterindwa Ainea, Estomih S Massawe, and Oluwole Daniel Makinde. Modelling the effect of treatment and infected immigrants on the spread of hepatitis c virus disease with acute and chronic stages. American Journal of Computational and Applied Mathematics, 2(1):10–20, 2012. [2] Mahamadou Alassane, Ouaténi Diallo, and Jérôme Pousin. Individual behavior and epidemiological model. Journal of Applied Mathematics and Bioinformatics, 1(2):57, 2011. [3] RM Anderson, GF Medley, RM May, and AM Johnson. A preliminary study of the transmission dynamics of the human immunodeficiency virus (hiv), the causative agent of aids. Mathematical Medicine and Biology: a Journal of the IMA, 3(4):229– 263, 1986. [4] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz. On the definition and the computation of the basic reproduction ratio r0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology, 28(4):365–382, 1990. [5] AB Gumel, Connell C McCluskey, and Pauline van den Driessche. Mathematical study of a staged-progression hiv model with imperfect vaccine. Bulletin of Mathe- matical Biology, 68(8):2105–2128, 2006. [6] James M Hyman and E Ann Stanley. Using mathematical models to understand the aids epidemic. Mathematical Biosciences, 90(1-2):415–473, 1988. [7] Mark Kot. Elements of mathematical ecology. Cambridge University Press, 2001. [8] Christopher M Kribs-Zaleta and Jorge X Velasco-Hernández. A simple vaccination model with multiple endemic states. Mathematical biosciences, 164(2):183–201, 2000. [9] Rachid Lounès and Héctor De Arazoza. A two-type model for the cuban national programme on hiv/aids. Mathematical Medicine and Biology: A Journal of the IMA, 16(2):143–154, 1999. [10] C Connell McCluskey. A model of hiv/aids with staged progression and amelioration. Mathematical biosciences, 181(1):1–16, 2003. [11] JYT Mugisha. The mathematical dynamics of hiv/aids epidemic in age-structured population. In Proceedings of XI SAMSA Symposium on the Potential of Mathemat- ical Modeling of Problems from the SAMSA Region, pages 8–23, 1997. [12] Z Mukandavire, W Garira, and C Chiyaka. Asymptotic properties of an hiv/aids model with a time delay. Journal of Mathematical Analysis and Applications, 330(2):916–933, 2007. REFERENCES 232 [13] Arni SR Srinivasa Rao. Mathematical modelling of aids epidemic in india. Current Science, 84(9):1192–1197, 2003. [14] Robert J Smith, Jing Li, Richard Gordon, and Jane M Heffernan. Can we spend our way out of the aids epidemic? a world halting aids model. BMC Public Health, 9(1):1–17, 2009.