EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 1, 2024, 59-92 ISSN 1307-5543 – ejpam.com Published by New York Business Global Mathematical modelling and Numerical Simulation of Hepatitis B Viral Infection: the case of Burkina Faso Adama Kiemtoré1, Wendabbo Olivier Sawadogo2,∗, Fatima Aqel3, Hamza Alaa4, Kounpielime Sosthène Somda5 1 LANIBIO Laboratory, University Joseph Ki-Zerbo, Ouagadougou, Burkina Faso 2 LARED Laboratory, University of Ouahigouya, Burkina Faso 3 Computer, Networks, Mobility and Modeling laboratory (IR2M), Faculty of Sciences and Technics, Hassan First University, Settat 26000, Morocco 4 LAMAI Laboratory, Faculty of Sciences and Technology, University of Cadi Ayyad, Abdelkrim, Khattabi, Marrakech 40000, Morocco 5 Department of Biochemistry and Microbiology, Pietro Annigonni Biomolecular Research Center, Ouagadougou, Burkina Faso Abstract. Hepatitis B is a viral infection that can cause inflammation of the liver and lead to severe liver damage and even death. The study of hepatitis in Burkina Faso is crucial for several reasons. Indeed, understanding the epidemiology of hepatitis in Burkina Faso can help develop effective prevention and control strategies. Its study can also contribute to a better understanding of the global burden of the disease and the development of effective interventions in other parts of the world. To this aim, a new differential susceptibility and infectivity mathematical model of Hepatitis B transmission was developed in order to simulate the potential spread of the Hepati- tis B virus in the population of Burkina Faso. Once the mathematical model is presented, the existence and uniqueness of non-negative solutions are proved. The model has a globally asymp- totically stable disease-free equilibrium when the basic reproduction number R0 < 1 and a globally asymptotically stable endemic equilibrium when R0 > 1. The global asymptotic stability of the disease-free equilibrium has been studied using the Castillo Chavez method. The Lyapunov func- tion and the LaSalle invariance principle are used to prove the global asymptotic stability of the endemic equilibrium. To simulate the proposed model, a Matlab numerical code has been devel- oped. Numerical simulations are performed using data of Burkina Faso. The obtained numerical results confirm analytical results as well as the evolution of hepatitis B in Burkina Faso. 2020 Mathematics Subject Classifications: 35A01, 92B05, 34C11, 47H10, 00A71. Key Words and Phrases: Hepatitis B, susceptibility, infectivity, stability, simulation ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v17i1.4987 Email addresses: kiemtoreadama@gmail.com (Adama Kiemtoré), wenddabo81@gmail.com (Wendabbo Olivier Sawadogo),fatima.aqel@uhp.ac.ma (Fatima Aqel), hamza.alaa@ced.uca.ma (Hamza Alaa), sosthenesomda@gmail.com (Kounpielime Sosthène Somda) https://www.ejpam.com 59 © 2024 EJPAM All rights reserved. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 60 1. Introduction Hepatitis B virus (HBV) is a global health threat, and its elimination by 2030 is a World Health Organization priority. Hepatitis B is not easily curable. Most people who are infected will not have any symptoms of infection [11], [12], [10]. They are called asymp- tomatic carriers and play a very important role in the transmission of this disease. These are people who are affected by the virus but do not show any clinical signs. [33] Serological testing is therefore essential. Thus, many people infected with HBV may not know it and may transmit it to other individuals [13], [9], [33], [41]. According to the World Health Organization (WHO), it has been estimated that 2.5 billion people are infected or have been infected with HBV, representing 30% million of the global population. Among them, there are 257 million people who suffer from chronic hepatitis B (global report 2017), of which 60 million are in Africa with 136,000 deaths in 2015 [33]. The HBV infection is the second known human carcinogenic after tobacco. The hepatitis B virus (HBV) primarily targets the liver, and infection with this virus can cause damage to the liver through either acute or chronic infection. The major damage from HBV is due to the long-term consequences of the chronic infection, mainly cirrhosis and hepatocellular carcinoma. According to the World Health Organization (WHO), about 25% of people with chronic hepatitis B infection develop cirrhosis or liver cancer, which can lead to a premature death if left untreated [19, 34]. Indeed, the risk of becoming chronic, and of complications such as liver cirrhosis and/or hepatocellular carcinoma (liver cancer), makes it a serious pathology since it causes more than 5000 deaths per year in Burkina Faso [32]. This latter is classified as one of the high-endemicity countries (seroprevalence > 8 %) [34]. Furthermore, in 2011’s Demo- graphic Health Survey, HBV seroprevalence was estimated to be 9.1% [31]. Also, a study conducted in Ouagadougou ((the capital city of Burkina Faso)) revealed a prevalence of 14.5 % in the general population [39]. Thus, nearly 2 million people in Burkina Faso are infected with HBV. Nevertheless, the actual number of infected people in Burkina Faso may be higher, as many people with chronic HBV infection may not know that they are infected due to the absence of symptoms or awareness of the disease. The distribu- tion of the prevalence of HBV infection in the regions of Burkina faso is shown in Figure 1. The main mode of hepatitis B virus (HBV) transmission in highly endemic regions such as Burkina Faso is through mother-to-child transmission during childbirth, or through close contact with infected family members in early childhood. This mode is known as vertical transmission of HBV. In the endemic regions, HBV is often transmitted from infected mothers to their babies during delivery, more particularly when the mother is tested positive for hepatitis B surface antigen (HBsAg) and hepatitis B e antigen (HBeAg) [33]. Transmission of HBV during perinatal period or the first years of childhood has a very high probability of being transmitted to a chronic disease. A pregnant woman with HBV can infect her child at delivery or during the first months of the infant’s life, either through breastfeeding or through extremely close contact (perinatal transmission). The W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 61 risk of transmission in the mother’s womb is related to the maternal viremia. The risk of infection with HBV may also be higher for individuals who engage in high-risk behaviours, such as injecting drugs, having unprotected sex with infected partners, or getting tattoos or piercings with unsterilised equipments. In addition, transmission between individuals living in the same household could be linked to the sharing of toiletries that cause small wounds, as well as razors, toothbrushes and nail clippers. Promiscuity and unsanitary conditions may be risk factors for the transmission of HBV, although they are not the only ones. Figure 1: Burkina Faso distribution of HBV prevalence [31] The hepatitis B virus is an enveloped virus belonging to the Hepadnaviridae family and the Orthohepadnavirus genus. It is a very contagious virus. The reservoir of HBV is strictly human, meaning that the virus only infects humans and has no other natural host. It can survive in the outdoor environment for at least 7 days and remain infectious on surfaces, such as needles or other medical equipments, for even longer [7]. Hepatitis B is considered as a DNA virus (deoxyribonuclease acid) with 3 antigenic systems, HBsAg (hepatitis B surface antigen) for the outer capsule, HBcAg (hepatitis B core antigen) and HBeAg (hepatitis B ”e” antigen) for the internal nucleocapsid. HBeAg is associated with high viral replication, but due to certain mutations, HBeAg is not detected in the serum although viral replication is presented. In clinical practices in Burkina Faso, HBeAg negative patients represent about 90%. Thus, the best marker of viral multiplication is the W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 62 detection of viral DNA (deoxyribonuclease acid) in the serum [33]. The high transmission of hepatitis B virus is explained by: • the length of the incubation period (6 weeks to 4 months) [28], • the very high infectious blood titer (0.0001 mL of plasma can transmit HBV) [16], • the large number of asymptomatic carriers [16], • HBV is 50 to 100 times more infectious than HIV [16], • HBV has the ability to retain its infectivity for more than 7 days at room temperature [16]. Figure 2: structure of an HBV virus particle (James A. Prekins, medical and scientific illustrations, 2002) [29] In this context, Therefore, which contribution could mathematics make to improve the situation in Burkina Faso? The answer to this question is to propose a mathematical model to better understand the dynamics of HBV transmission and predict the infection’s evolution. The remainder of this paper is organized as follows: section two is dedicated to the presentation of the proposed transmission network model and to the calculus of the re- production rate. Section three is devoted to the proof of global existence and uniqueness of solutions for the proposed system, as well as to the calculus of the reproduction rate. The global stability analysis of the model problem for the disease free equilibrium point is presented in the fourth section. In the fifth section, the global stability analysis for W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 63 the endemic equilibrium point is proved. In the sixth section, a numerical algorithm is developed to numerically simulate the proposed model for different involved parameters. The obtained numerical results are presented in the sixth section. 2. The proposed model of HBV transmission in Burkina Faso A SEIR (Susceptible S, Exposed E, Infected I, Recovered R) model was formulated which takes into account the epidemiology of HBV infection. Due to the fact that vertical transmission is the most frequent mode of transmission in Burkina Faso, a model with vertical transmission was developed. The total population births is denoted by Λ. The introduction of vertical transmission decreases the number of births by an amount of λII + λCC which no longer become susceptible, as babies from these births are acutely and chronically infected by vertical transmission. Thus, both λII and λCC appear in acute EI and chronic EC compartment of acutely infected exposed [16]. HBV is transmitted vertically by infected individuals (in the acute phase) and chronic carriers of the virus. The description of the variables is given in Table 1 and that of the parameters in Table 2. 2.1. Assumptions Throughout this paper, it is assumed that : (H1) Susceptibles individuals were subdivided into three age classes Si, i = 1, 2, 3. (H2) The population in each compartment is likely to fall into two groups, with a pro- portion αi, i = 1, 2, 3 of this population developing an acute infection i.e. the acute exposed EI who will evolve to an acute infected state I and a proportion (1 − αi), i = 1, 2, 3 of this population that will become chronically infected (the chronic exposed EC that will evolve to a chronic infected state C ). (H3) The high prevalence (> 8%) in Burkina Faso is mainly due to the vertical transmis- sion, hence the choice of the vertical transmission model [16, 31]. (H4) The general recruitment of the population is Λ− λII − λCC. (H5) A Part of the births from the acute and chronic HBV infections participate with an amount λI + λC , which will move to the exposed compartments of both acute and chronic infections. (H6) Vertical transmission occurs in both acute and chronically infected individuals. (H7) The adequate contact coefficient for a susceptible person to be infected with HBV, depends on the number of infected individuals [2, 16], i.e. the EI , EC , I, C. One has the following (i) β1 = β1,1EI +β1,2EC +β1,3I+β1,4C the adequate coefficient of contact so that a susceptible S1 becomes contaminated by the HBV. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 64 (ii) β2 = β2,1EI + β2,2EC + β2,3I + β2,4C the appropriate contact coefficient for a susceptible S2 to become contaminated with HBV. (iii) β3 = β3,1EI + β3,2EC + β3,3I + β3,4C the appropriate contact coefficient for a susceptible S3 to become contaminated with HBV. (H8) Acute infection includes fulminant hepatitis, which is very rare [28]. (H9) All individuals who do not die from HVB, i.e. individuals in compartments S1, S2, S3, EI , EC , R have the same mortality rate µ. (H10) Λ, λI , λC , µ, µ1, µ2, µ3, µ4, δ1, δ2, γi, γc are assumed to be strictly positive. (H11) The prognosis for the evolution of HBV infection to chronicity in Burkina Faso is generally high, as it is in many other sub-Saharan African countries. The prognosis depends on the age at which the individual was infected. Thus, susceptible infants aged 0-1 year become chronic with a probability of 0.9; children aged 1-5 years with a probability of 0.3; and those aged over 5 years with a probability of 0.05 [17]. (H12) The total population N is not constant. 2.2. The interactions In the following sub-section, we start by presenting the different parameters, and the different variables involved in our mathematical model. Next, detailed explanations of the different compartments that compose the proposed model as well as their interactions, are given. Table 1 contains the different variables with their meanings, while Table 2 contains the different parameters with their significance. Table 1: Variables used in the model Variables description S1(t) babies aged 0-1 year susceptible to disease S2(t) infants 1 to 5 years old susceptible to disease S3(t) children over 5 years of age, adolescents and adults at risk of disease EI(t) Exposed infected individuals who will progress to an acute infection state EC(t) Exposed chronic individuals who will progress to a chronically infected state I(t) infected individuals acute C(t) chronic carriers of the virus R(t) recovered persons W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 65 Table 2: Parameters used in the model Parameters description Λ births of total population λI proportion of births from the infected that will become infected acute λC proportion of births from the chronic that will become chronic p1 Proportion of infants aged 0 − 1year that grow healthy to enter the S2 compartment without being contaminated p2 Proportion of infants aged 1− 5 years who grow up healthy to enter the compartment S3 without being contaminated β1,i; i = 1, .., 4 Adequate contact coefficient for people likely S1 to be contaminated. β2,i; i = 1, .., 4 Adequate contact coefficient for individuals likely S2 to be contaminated. β3,i; i = 1, .., 4 Adequate contact coefficient for individuals likely S3 to be contaminated. α1 Probability of S1 that will be contaminated, and become exposed infected EI 1− α1 Probability of S1 becoming contaminated, and exposed chronic EC α2 Probability of S2 that will be contaminated, and become exposed infected EI 1− α2 Probability of S2 that will be contaminated and become exposed chronic EC α3 Probability of S3 becoming exposed to infection EI 1− α3 Probability of S3 that will be contaminated and become chronically exposed EC δ1 Proportion of exposed infected who become infected (acute) I. δ2 Proportion of exposed infected that become chronic C. γI Proportion of acute patients that are cured recovered from the infection process. γC Proportion of chronic carriers cured and recovered from the infection process. µ Natural mortality rate, i.e., mortality not caused by HBV. µI Death rate caused by acute HBV infection µC Mortality rate caused by chronic HBV infection The model is composed of eight compartments. In Burkina Faso, the age at which an individual is infected influences his or her fate in terms of the disease’s chronicity. It is therefore necessary to divide susceptible individuals into three age classes, i.e. three compartments of susceptible individuals. Dividing the susceptible population into 3 com- partments provides a good approximation of the reality of hepatitis B epidemiology in Burkina Faso. • The compartment S1 represents the class of infants aged 0 to 1 year. In this com- partment, the total population births is Λ− λII − λCC. It emerges the infants who grow up with a proportion p1 to enter the compartment S2, the natural mortality µS1, and also the babies who are infected by vertical transmission with a coefficient β1,i, i = 1, ..., 4. The latter become exposed infected EI with probability α1, and chronic exposed EC with probability (1− α1). • The compartment S2 represents children aged 1 to 5 years. In this compartment, W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 66 infants from compartment S1 with a proportion p1 enter, and babies who grow up with a proportion of p2 move to compartment S3, natural mortality µS2 and those infected with a coefficient β1,i for i = 1, ..., 4. The latter ones become exposed infected EI with probability (1 − α2), and chronic exposed EC with probability (1− α2). • The S3 compartment is composed of children over 5 years old, adolescents and sus- ceptible adults. In this compartment, the input stream is p2S2, which represents children in compartment S2 who grow up with a proportion p2. The output stream is the natural mortality of susceptible S3, and infected individuals who become in- fected with the coefficient β3,i, i = 1, 2, 3. The latter ones become exposed infected EI with probability α3, and chronically exposed EC with probability (1− α3). There are two compartments of exposed people: • The EI compartment represents exposed infected individuals, i.e. exposed that will evolve towards an infectious state. The entry of this compartment is characterized by the susceptible that have been infected according to the coefficients αiβi,j , i = 1, 2, 3, j = 1, 2, 3, and acute contaminated people λII who are infected by vertical trans- mission. The output is composed of natural mortality EI , as well as the evolution towards an infectious state with a proportion δ1. • The EC compartment represents the chronically exposed, i.e. exposed that will evolve towards a chronic state. The input of this compartment is characterized by susceptible individuals who have been infected according to the coefficients (1 − αi)βi,j for (i = 1, 2, 3 and j = 1, 2, 3). Also, the input of the EC compartment is characterized by individuals who are chronically infected through vertical transmis- sion λCC. The exit is constituted on one hand by the natural mortality EC , and on the other hand by the evolution towards an infectious state with a proportion δ2. Three further compartments remain: • The acute infected compartment I, characterized by an input δ1EI , and an output composed of mortality (µ+µI)I (natural mortality and mortality due to the disease), and recovery γII. • The compartment represents the chronically infected C. The entry of this compart- ment is EC , and the output is composed of mortality (µ+ µC)C (natural mortality and mortality due to the disease), and recovery γCC. • The compartment R (Recovered) represents individuals who are cured of the disease. There enters the cured (γII and γCC) and emerges the mortality µR. The described interactions are illustrated in Figure 3. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 67 Figure 3: HBV model with vertical transmission 3. Mathematical Analysis of the model 3.1. Mathematical model By doing a mass balance through the different considered compartments, one obtains the following system: dS1 dt = Λ− λII − λCC − µS1 − p1S1 − β1,1EIS1 − β1,2ECS1 − β1,3IS1 − β1,4CS1 dS2 dt = p1S1 − µS2 − β2,1EIS2 − β2,2ECS2 − β2,3IS2 − β2,4CS2 − p2S2 dS3 dt = p2S2 − µS3 − β3,1EIS3 − β3,2EIS3 − β3,3IS3 − β3,4CS3 dEI dt = α1(β1,1EIS1 + β1,2ECS1 + β1,3IS1 + β1,4CS1) + α2(β2,1EIS2 + β2,2ECS2 + β2,3IS2 + β2,4CS2) +α3(β3,1EIS3 + β3,2ECS3 + β3,3IS3 + β3,4CS3)− µEI − δ1EI + λII dEC dt = (1− α1)(β1,1EIS1 + β1,2ECS1 + β1,3IS1 + β1,4CS1) + (1− α2)(β2,1EIS2 + β2,2ECS2 + β2,3IS2 +β2,4CS2) + (1− α3)(β3,1EIS3 + β3,2ECS3 + β3,3IS3 + β3,4CS3) + λCC − µEC − δ2EC dI dt = δ1EI − µI − µII − γII dC dt = δ2EC − µC − µCC − γCC dR dt = γII + γCC − µR (1) W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 68 The total number N(t) of the population studied is not constant and: N(t) = S1(t) + S2(t) + S3(t) + EI(t) + EC(t) + I(t) + C(t) +R(t) (2) Since people and viruses have different dimensions, they therefore have to carry out the following normalization: s1 = S1 N0 , s2 = S2 N0 , s3 = S3 N0 , ei = EI N0 , ec = EC N0 , i = I N0 , c = C N0 , r = R N0 , λ = Λ N0 , bi,j = βi,jN0, i = 1, ..., 3; j = 1, ..., 4 Where N0 = N(0) is the total number of the initial population. Finally the normalized model becomes:  ds1 dt = λ− iλi − cλc − (µ+ p1)s1 − b1,1eis1 − b1,2ecs1 − b1,3is1 − b1,4cs1 ds2 dt = p1s1 − (µ+ p2)s2 − b2,1eis2 − b2,2ecs2 − b2,3is2 − b2,4cs2 ds3 dt = p2s2 − µs3 − b3,1eis3 − b3,2ecs3 − b3,3is3 − b3,4cs3 dei dt = α1s1(b1,1ei + b1,2ec + b1,3i+ b1,4c) + α2s2(b2,1ei + b2,2ec + b2,3i+ b2,4c) +α3s3(b3,1ei + b3,2ec + b3,3i+ b3,4c)− µei − δ1ei + iλi dec dt = (1− α1)s1(b1,1ei + b1,2ec + b1,3i+ b1,4c) + (1− α2)s2(b2,1ei + b2,2ec + b2,3i+ b2,4c) +(1− α3)s3(b3,1ei + b3,2ec + b3,3i+ b3,4c) + cλc − µec − δ2ec di dt = δ1ei − µi− γii− µii dc dt = δ2ec − µc− γcc− µcc dr dt = γii+ γcc− µr (3) In the following, the positivity and existence of solutions for the normalized system (3) is given. 3.2. Positivity, boundedness and global existence of solutions Theorem 1. Let the initial value (s1,0, s2,0, s3,0, ei,0, ec,0, i0, c0, r0) ∈ R8 +, such that s1,0 + s2,0 + s3,0 + ei,0 + ec,0 + i0 + c0 + r0 = 1, and 0 ≤ λ, λi, λc, µ, p1, p2, µi, µc, δ1, δ2, γi, γc ≤ 1 0 ≤ b1,i, b2,i, b3,i ≤ 1 for all 1 ≤ i ≤ 4, then there exists a unique, nonnegative, bounded global solution to system (3). Moreover for all t ≥ 0 W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 69 i) 0 ≤ n(t) ≤ λ µ + 1 (4) where n(t) = s1(t) + s2(t) + s3(t) + ei(t) + ec(t) + i(t) + c(t) + r(t) ii) Moreover, if s1,0 ≤ λ µ+ p1 , s2,0 ≤ p1λ (µ+ p1)(µ+ p2) and s3,0 ≤ p1p2λ µ(µ+ p1)(µ+ p2) then s1(t) ≤ λ µ+ p1 , s2(t) ≤ p1λ (µ+ p1)(µ+ p2) , s3(t) ≤ p1p2λ µ(µ+ p1)(µ+ p2) . Proof. i) For the local existence, all the functions of system (3) are locally Lipschitz contin- uous. Thus, there exists a unique local solution on t ∈ [0, Tmax), where Tmax is the explosion time. The analysis of this kind of system is based on elementary methods of ordinary differential equations. The existence of unique solutions is guaranteed by various fixed point theorems on a maximal interval [0, Tmax). By proving that the components of the solution vector (s1(t), s2(t), s3(t), ei(t), ec(t), i(t), c(t), r(t)) are uniformly bounded on any bounded interval [0, Tmax), one ensures that Tmax = ∞. We remark that the components of the vector F (s1, s2, s3, ei, ec, i, c, r) =  F1(s1, s2, s3, ei, ec, i, c, r) F2(s1, s2, s3, ei, ec, i, c, r) F3(s1, s2, s3, ei, ec, i, c, r) F4(s1, s2, s3, ei, ec, i, c, r) F5(s1, s2, s3, ei, ec, i, c, r) F6(s1, s2, s3, ei, ec, i, c, r) F7(s1, s2, s3, ei, ec, i, c, r) F8(s1, s2, s3, ei, ec, i, c, r)  (5) Where F1(s1, s2, s3, ei, ec, i, c, r) = λ− iλi − cλc − (µ+ p1)s1 − b1,1eis1 − b1,2ecs1 − b1,3is1 − b1,4cs1 F2(s1, s2, s3, ei, ec, i, c, r) = p1s1 − (µ+ p2)s2 − b2,1eis2 − b2,2ecs2 − b2,3is2 − b2,4cs2 F3(s1, s2, s3, ei, ec, i, c, r) = p2s2 − µs3 − b3,1eis3 − b3,2ecs3 − b3,3is3 − b3,4cs3 F4(s1, s2, s3, ei, ec, i, c, r) = α1s1(b1,1ei + b1,2ec + b1,3i+ b1,4c) + α2s2(b2,1ei + b2,2ec + b2,3i+ b2,4c)+ + α3s3(b3,1ei + b3,2ec + b3,3i+ b3,4c) + λii− µei − δ1ei W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 70 F5(s1, s2, s3, ei, ec, i, c, r) = (1− α1)s1(b1,1ei + b1,2ec + b1,3i+ b1,4c) + (1− α2)s2(b2,1ei + b2,2ec + b2,3i+ b2,4c) + +(1− α3)s3(b3,1ei + b3,2ec + b3,3i+ b3,4c)+ + λcc− µec − δ2ec F6(s1, s2, s3, ei, ec, i, c, r) = δ1ei − µi− γii− µii F7(s1, s2, s3, ei, ec, i, c, r) = δ2ec − µc− γcc− µcc F8(s1, s2, s3, ei, ec, i, c, r) = γii+ γcc− µr are quasi-positive. Consequently, since the initial conditions are nonnegative, this implies that the solution components are nonnegative for all t ∈ [0, Tmax). Now, let the function n be defined as n(t) = s1(t) + s2(t) + s3(t) + ei(t) + ec(t) + i(t) + c(t) + r(t) By taking the sum of the first eight equations in (S), we observe dn dt ≤ λ− µS(t) n(0) = 1. (6) Integrating equation (6) over (0, t) for all 0 < t < T , one can get the following n(t)eµt − 1 ≤ λ µ (eµt − 1), which implies that n(t) ≤ e−µt + λ µ (1− e−µt). Therefore n(t) ≤ (1− λ µ )e−µt + λ µ . Here, two different cases are distinguished. If λ µ < 1, then the following inequality is satisfied n(t) ≤ 1− λ µ + λ µ ≤ 1, otherwise, if λ µ ≥ 1 then n(t) ≤ λ µ . W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 71 Finally, one can get the following n(t) ≤ λ µ + 1. (7) Hence, Tmax = ∞ and the existence of unique, non-negative and bounded global solution are proved. ii) We remark that s1 satisfies the following ds1 dt ≤ λ− (µ+ p1)s1(t) s1(0) = s1,0. (8) By integrating (8) over (0, t) for all 0 < t < T , we obtain s1(t)e (µ+p1)t ≤ λ (µ+ p1) (e(µ+p1)t − 1) + s1,0. Since s1,0 ≤ λ (µ+ p1) , we obtain s1(t) ≤ λ (µ+ p1) . The same reasoning is applied for s2 and s3 which concludes the proof of this theorem. 3.3. Reproduction rate and its interpretations Infectious diseases are generally modeled mathematically by compartmental models. Individuals can move from one compartment to another, and the population is represented by compartments with labels. The first models date back to the beginning of the 20th century by [36] in 1916, [37, 38] in 1917, [22] in 1927 and [21] in 1956. Ordinary differential equations (deterministic, stochastic or fractional) are used to an- alyze these compartmental models mathematically. These mathematical models can be used to predict the spread of a disease, the total number of people infected or the dura- tion of an epidemic, as well as to estimate various epidemiological parameters such as the reproduction rate [1], [3]. This section is devoted to the calculation of the reproduction number R0 of our proposed model. For this we will use the generation matrix method which was developed by [8] and then adopted by [40] for finite dimensional systems. Theorem 2. Consider the system (3) with the given parameters 0 ≤ λ, λi, λc, µ, p1, p2, µi, µc, δ1, δ2, γi, γc ≤ 1 0 ≤ b1,i, b2,i, b3,i ≤ 1 for all 1 ≤ i ≤ 4, Then, W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 72 (i) the disease free equilibrium point (DFE) is X0 = ( λ µ+ p1 , p1λ (µ+ p1)(µ+ p2) , p1p2λ µ(µ+ p1)(µ+ p2) , 0, 0, 0, 0, 0 ) , (9) (ii) the reproduction rate is R0 = η1,1 + η2,2 + ( (η1,1 − η2,2) 2 + 4η1,2η2,1 ) 1 2 2 (10) Where η1,1, η2,2, η1,2, η2,1 are given by (17) Proof. Here, we consider the proposed mathematical model (3) with eight homogeneous compartments. This model can be written as d dt (s1, s2, s3, ei, ec, i, c, r) = F (s1, s2, s3, ei, ec, i, c, r). Where F is defined by (5). The point X0 defined by (9) satisfies F (X0) = 0. Infected and chronic exposed are infected but not infectious, so they do not participate in the transmission of the virus. Thus bk,1 = bk,2 = 0; k = 1, 2, 3. To compute the repro- duction rate R0, one must consider only the infected and infectious compartments that satisfy the following fourth-order system : d dt  ei ec i c  =  α1s1(b1,3i+ b1,4c) + α2s2(b2,3i+ b2,4c) + α3s3(b3,3i+ b3,4c)− µei − δ1ei + iλi (1− α1)s1(b1,3i+ b1,4c) + (1− α2)s2(b2,3i+ b2,4c) + (1− α3)s3(b3,3i+ b3,4c) + cλc − µec − δ2ec δ1ei − µi− γii− µii δ2ec − µc− γcc− µcc  . The rate of occurrence of new infection in the four compartments (ei, ec, i, c) is repre- sented by the vector U as follows U =  α1s1(b1,3i+ b1,4c) + α2s2(b2,3i+ b2,4c) + α3s3(b3,3i+ b3,4c) + λii (1− α1)s1(b1,3i+ b1,4c) + (1− α2)s2(b2,3i+ b2,4c) + (1− α3)s3(b3,3i+ b3,4c) + λcc 0 0  (11) W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 73 and the transfer rate of individuals into and out of the infected compartments is given by the vector T T =  µei + δ1ei µec + δ2ec −δ1ei + µi+ γii+ µii −δ2ec + µc+ γcc+ µcc  (12) Thus, all the epidemiological events leading to new infections are incorporated into the model via the matrix Θ, while all the other events are included in the matrix Γ. Progression to either death or immunity ensures that Γ is invertible. By substituting the variables (s1, s2, s3) with the free equilibrium point X0, we obtain the two matrices Θ and Γ which are expressed as follows : Θ =  0 0 Θ1,3 Θ1,4 0 0 Θ2,3 Θ2,4 0 0 0 0 0 0 0 0  (13) Where Θ1,3 = [ α1b1,3 + α2b2,3 p1 µ+ p2 + α3b3,3 p1p2 µ(µ+ p2) ] λ µ+ p1 + λi Θ1,3 = [ α1b1,4 + α2b2,4 p1 µ+ p2 + α3b3,4 p1p2 µ(µ+ p2) ] λ µ+ p1 Θ2,3 = [ (1− α1)β1,3 + (1− α2)b2,3 p1 µ+ p2 + (1− α3)b3,3 p1p2 µ(µ+ p2) ] λ µ+ p1 Θ2,4 = [ (1− α1)b1,4 + (1− α2)b2,4 p1 µ+ p2 + (1− α3)b3,4 p1p2 µ(µ+ p2) ] λ µ+ p1 + λc Γ =  µ+ δ1 0 0 0 0 µ+ δ2 0 0 −δ1 0 µ+ γi + µi 0 0 −δ2 0 µ+ γc + µc  (14) By calculating Γ−1, we obtain : W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 74 Γ−1 =  1 µ+ δ1 0 0 0 0 1 µ+ δ2 0 0 δ1 (µ+ γi + µi)(µ+ δ1) 0 1 µ+ γi + µi 0 0 δ2 (µ+ γc + µc)(µ+ δ2) 0 1 µ+ γc + µc  (15) We calculate the matrix of the new generation, we obtain : ΘΓ−1 =  η1,1 η1,2 η1,3 η1,4 η2,1 η2,2 η2,3 η2,4 0 0 0 0 0 0 0 0  (16) Where η1,1 = [( α1β1,3 + α2b2,3 p1 µ+ p2 + α3b3,3 p1p2 µ(µ+ p2) ) λ µ+ p1 + λi ] δ1 (µ+ γi + µi)(µ+ δ1) η1,2 = [ α1b1,4 + α2b2,4 p1 µ+ p2 + α3b3,4 p1p2 µ(µ+ p2) ] λ (µ+ p1) δ2 (µ+ γc + µc)(µ+ δ2) η1,3 = [( α1b1,3 + α2b2,3 p1 µ+ p2 + α3b3,3 p1p2 µ(µ+ p2) ) λ µ+ p1 + λi ] 1 µ+ γi + µi η1,4 = [ α1b1,4 + α2b2,4 p1 µ+ p2 + α3b3,4 p1p2 µ(µ+ p2) ] λ (µ+ p1)(µ+ γc + µc) η2,1 = [ (1− α1)b1,3 + (1− α2)b2,3 p1 µ+ p2 + (1− α3)b3,3 p1p2 µ(µ+ p2) ] λ (µ+ p1) δ1 (µ+ γi + µi)(µ+ δ1) η2,2 = [( (1− α1)b1,4 + (1− α2)b2,4 p1 µ+ p2 + (1− α3)b3,4 p1p2 µ(µ+ p2) ) λ µ+ p1 + λc ] δ2 (µ+ γc + µc)(µ+ δ2) η2,3 = [ (1− α1)b1,3 + (1− α2)b2,3 p1 µ+ p2 + (1− α3)b3,3 p1p2 µ(µ+ p2) ] λ (µ+ p1)(µ+ γi + µi) η2,4 = [( (1− α1)b1,4 + (1− α2)b2,4 p1 µ+ p2 + (1− α3)b3,4 p1p2 µ(µ+ p1) ) λ (µ+ p1) + λc ] 1 µ+ γc + µc (17) By the next generation matrix approach, R0 for the model. By calculating the spectral radius of the next generation matrix ΘΓ−1, we get : W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 75 R0 = ρ ( ΘΓ−1 ) = max{|ωi|, i = 1, 2, 3} = max  ∣∣∣∣∣∣∣ η1,1 + η2,2 + ( (η1,1 − η2,2) 2 + 4η1,2η2,1 ) 1 2 2 ∣∣∣∣∣∣∣ , ∣∣∣∣∣∣∣ η1,1 + η2,2 − ( (η1,1 − η2,2) 2 + 4η1,2η2,1 ) 1 2 2 ∣∣∣∣∣∣∣ , 0  (18) Then the reproduction rate R0 of the normalized system is given as : R0 = η1,1 + η2,2 + ( (η1,1 − η2,2) 2 + 4η1,2η2,1 ) 1 2 2 (19) Where η1,1, η2,2, η1,2, η2,1 are given by (17). Remark 1. The reproduction number R0 has a crucial role in the control of Hepatitis B viral infection. When this number is less than 1, an infected individual infects less than one other individual on average, therefore the disease is under control. Conversely, when this number exceeds 1, the disease spreads through the population and becomes epidemic [40]. 4. Global stability analysis of the model problem for the disease free equilibrium point In the section, the global stability analysis of the model for the disease free equilibrium is shown. In the following, we aim to provide a brief investigation of the Castillo Chavez technique [5, 6] to prove the stability of system (3) in the global sense at the disease free equilibrium point. Therefore, by applying the Castillo Chavez technique [6], the given problem (3) is converted into the following sub-models: dX1 dt = F (X1, X2), dX2 dt = G(X1, X2), G(X1, 0) = 0. (20) Where X1 and X2 designate the population of uninfected individuals, and infected indi- viduals, respectively. The following conditions (A1) and (A2) must be satisfied to guarantee the local asymptotic stability. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 76 (A1) If dX1 dt = F (X1, 0) then X0 is globally asymptotically stable. (A2) G(X1, X2) = AX2 − Ĝ(X1, X2). Where Ĝ(X1, X2) ≥ 0 for (X1, X2) ∈ Ω, the matrix A is a M-matrix whose off diagonal elements are nonnegative. In the proposed system (3), X1 = (s1, s2, s3, r) ∈ R4 and X1 = (ei, ec, i, c) ∈ R4. According to the results obtained in Section 3, the free disease equilibrium point (DFE) was denoted by X0, it is defined by the following sense X0 = ( λ µ+ p1 , p1λ (µ+ p1)(µ+ p2) , p1p2λ µ(µ+ p1)(µ+ p2) , 0, 0, 0, 0, 0 ) which can be written as X0 = (s0, 0). In order to ensure the globally asymptotically stability of DFE point, the results given above [6] were applied. Theorem 3. If R0 < 1 and s1,0 ≤ λ µ+ p1 , s2,0 ≤ p1λ (µ+ p1)(µ+ p2) , s3,0 ≤ p1p2λ µ(µ+ p1)(µ+ p2) , then the DFE point X0 of the model (3) is globally asymptotically stable. Proof. In System 20, since infected and chronically exposed individuals are infected but not infectious. Therefore, they do not participate in the transmission of the virus. Hence bk,1 = 0 and bk,2 = 0 for k = 1, 2, 3. Here, one can set the following F (X1, X2) =  λ− iλ1 − cλ2 − (µ+ p1)s1 − b1,3is1 − b1,4cs1 p1s1 − (µ+ p2)s2 − b2,3is2 − b2,4cs2 p2s2 − µs3 − b3,3is3 − b3,4cs3 γ1i+ γ2c− µr  and G(X1, X2) =  G1(X1, X2) G2(X1, X2) G3(X1, X2) G4(X1, X2)  Where each component is defined by the following sense G1(X1, X2) = α1s1(b1,3i+ b1,4c) + α2s2(b2,3i+ b2,4c) + α3s3(b3,3i+ b3,4c)− µei − δ1ei + iλi G2(X1, X2) = (1− α1)s1(b1,3i+ b1,4c) + (1− α2)s2(b2,3i+ b2,4c) + (1− α3)s3(b3,3i+ b3,4c) + cλc − µec − δ2ec G3(X1, X2) = δ1ei − µi− γii− µii G4(X1, X2) = δ2ec − µc− γcc− µcc. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 77 At the DFE point, it is clear that G(X1, 0) = 0. Now, we need to prove that (A1) is satified. To this aim, we calculate the eignvalues of the Jacobian matrix associated to F at the DFE point. DX1F (X0) =  −(µ+ p1) 0 0 0 p1 −(µ+ p2) 0 0 0 p2 −µ 0 0 0 0 −µ  (21) Since the eigenvalues ζ1 = −(µ + p1), ζ2 = −(µ + p2), ζ3 = −µ, ζ4 = −µ of the jacobian matrix are negatives, then the DFE point X0 is globally asymptotically stable. To derive the condition (A2), we first calculate the matrix A = DX2G(X0) at the DFE point A =  −(µ+ δ1) 0 λi +Θ1 Θ2 0 −(µ+ δ2) Θ3 λc +Θ4 δ1 0 −(µ+ γi + µi) 0 0 δ2 0 −(µ+ γc + µc)  (22) Where Θ1 = α1b1,3 λ µ+ p1 + α2b2,3 p1λ (µ+ p1)(µ+ p2) + α3b3,3 p1p2λ µ(µ+ p1)(µ+ p2) Θ2 = α1b1,4 λ µ+ p1 + α2b2,4 p1λ (µ+ p1)(µ+ p2) + α3b3,4 p1p2λ µ(µ+ p1)(µ+ p2) Θ3 = (1− α1)b1,3 λ µ+ p1 + (1− α2)b2,3 p1λ (µ+ p1)(µ+ p2) + (1− α3)b3,3 p1p2λ µ(µ+ p1)(µ+ p2) Θ4 = (1− α1)b1,4 λ µ+ p1 + (1− α2)b2,4 p1λ (µ+ p1)(µ+ p2) + (1− α3)b3,4 p1p2λ µ(µ+ p1)(µ+ p2) (23) It is clear that the matrix A given above is an M-matrix. Now, one can calculate the following function Ĝ(X1, X2). We have Ĝ(X1, X2) =  Ĝ1(X1, X2) Ĝ2(X1, X2) 0 0  (24) W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 78 Ĝ1(X1, X2) = −α1s1(b1,3i+ b1,4c)− α2s2(b2,3i+ b2,4c)− α3s3(b3,3i+ b3,4c) + iΘ1 + cΘ2 Ĝ2(X1, X2) = −(1− α1)s1(b1,3i+ b1,4c)− (1− α2)s2(b2,3i+ b2,4c)− (1− α3)s3(b3,3i+ b3,4c)+ + iΘ3 + cΘ4. By replacing Θ1, Θ2, Θ3 and Θ4 by their expressions given in (23), one obtains the following result Ĝ1(X1, X2) = α1b1,3i( λ µ+ p1 − s1) + α2b2,3i( p1λ (µ+ p1)(µ+ p2) − s2) + α3b3,3i( p1p2λ µ(µ+ p1)(µ+ p2) − s3)+ + α1b1,4c( λ µ+ p1 − s1) + α2b2,4c( p1λ (µ+ p1)(µ+ p2) − s2) + α3b3,4c( p1p2λ µ(µ+ p1)(µ+ p2) − s3) Ĝ2(X1, X2) = (1− α1)b1,3i( λ µ+ p1 − s1) + (1− α2)b2,3i( p1λ (µ+ p1)(µ+ p2) − s2)+ + (1− α3)b3,3i( p1p2λ µ(µ+ p1)(µ+ p2) − s3) + (1− α1)b1,4c( λ µ+ p1 − s1)+ + (1− α2)b2,4c( p1λ (µ+ p1)(µ+ p2) − s2) + (1− α3)b3,4c( p1p2λ µ(µ+ p1)(µ+ p2) − s3) Thanks To Theorem 1, ii), one gets the following result λ µ+ p1 − s1 ≥ 0, p1λ (µ+ p1)(µ+ p2) − s2 ≥ 0, p1p2λ µ(µ+ p1)(µ+ p2) − s3 ≥ 0. Since i, c, b1,3, b1,3, b2,3, b3,3, b1,4, b2,4, b3,3 are nonnegatives, and 0 ≤ αk ≤ 1 for k = 1, 2, 3, the following result is obtained Ĝ1(X1, X2) ≥ 0, Ĝ2(X1, X2) ≥ 0. Consequently, the hypothesis (A1) and (A2) are satisfied. Moreover, one uses Castillo Chavez [6] technique to conclude that if R0 < 1 then the DFE point is globally asymp- totically stable. 5. Endemic Equilibrium Point E∗ In this section, the global stability analysis of the proposed model problem for the endemic equilibrium point is studied. To prove this result, Lyapunov functions are used. The concept of global asymptotic stability refers to the property that if a system starts near the endemic point, it will eventually converge to that point. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 79 5.1. Exixtence of E∗ Theorem 4. System (1) admits a unique positive emdemic equilibrium E∗ = (S∗ 1 , S ∗ 2 , S ∗ 3 , E ∗ I , E ∗ C , I ∗, C∗, R∗) whenever R0 > 1 Proof. By seting the right hand side of system (1) equal to zero, and keeping each state variable non-zero (S1 ̸= 0, S2 ̸= 0, S3 ̸= 0, EI ̸= 0, EC ̸= 0, I ̸= 0, C ̸= 0, R ̸= 0), then one obtains S∗ 1 = Λ− ΛII ∗ − ΛCC ∗ µ+ p1 + β1,1E∗ I + β1,2E∗ C + β1,3I∗ + β1,4C∗ , S∗ 2 = p1S ∗ 1 µ+ p2 + β2,1E∗ I + β2,2E∗ C + β2,3I∗ + β2,4C∗ , S∗ 3 = p2S ∗ 2 µβ3,1E∗ I + β3,2E∗ C + β3,3I∗ + β3,4C∗ , E∗ I = ΛII ∗ + ΛCC ∗α1 (β1,2E ∗ CS ∗ 1 − β1,3I ∗S∗ 1 − β1,4C ∗S∗ 1) µ+ δ1 − α1β1,1S∗ 1 − α2β2,1S∗ 2 − α3β3,1S∗ 3 + α2(β2,2E ∗ CS ∗ 2 − β2,3I ∗S∗ 2 − β2,4C ∗S∗ 2) + α3(β3,2E ∗ CS ∗ 3 − β3,3I ∗S∗ 3 − β3,4C ∗S∗ 3 + β2,5T ∗S∗ 3) µ+ δ1 − α1β1,1S∗ 1 − α2β2,1S∗ 2 − α3β3,1S∗ 3 , E∗ C = ΛII ∗ + λCC ∗ + (1− α1)(β1,1E ∗ IS ∗ 1 − β1,3I ∗S∗ 1 − β1,4C ∗S∗ 1) µ+ δ2 − α1β1,1S∗ 1 − α2β2,1S∗ 2 − α3β3,1S∗ 3 + (1− α2)(β2,1E ∗ IS ∗ 2 − β2,3I ∗S∗ 2 − β2,4C ∗S∗ 2) + (1− α3)(β3,1E ∗ IS ∗ 3 − β3,3I ∗S∗ 3 − β3,4C ∗S∗ 3) µ+ δ2 − α1β1,1S∗ 1 − α2β2,1S∗ 2 − α3β3,1S∗ 3 , I∗ = δ1E ∗ I µ+ γI + µI , C∗ = δ2E ∗ C µ+ γC + µC , R∗ = γII ∗ + γCC ∗ µ 5.2. Global stability of E∗ In this subsubsection, by constructing suitable Lyapunov function, we will prove the global asymptotic stability of the endemic equilibrium point E∗. Theorem 5. If R0 > 1, the global endemic equilibrium point E∗ of the system (1) is globally asymptotically stable.. Proof. When R0 > 1, one can define the following Lyapunov function as in [4, 23, 25, 26, 30] : W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 80 V = (S1 − S∗ 1) + (S2 − S∗ 2) + (S3 − S∗ 3) + (EI − E∗ I ) + (EC − E∗ C) + (I − I∗)+ + (C − C∗) + (R−R∗)− (S∗ 1 + S∗ 2 + S∗ 3 + E∗ I + E∗ C + I∗ + C∗ ++R∗)× ln ( S1 + S2 + S3 + EI + EC + I + C +R S∗ 1 + S∗ 2 + S∗ 3 + E∗ I + E∗ C + I∗ + C∗ +R∗ ) (25) As N = S1 + S2 + S3 + EI + EC + I + C +R, one can set N∗ = S∗ 1 + S∗ 2 + S∗ 3 + E∗ I + E∗ C + I∗ + C∗ +R∗. Then, the Lyapunov function can also be rewriten as follows: V = N −N∗ −N∗ln N N∗ V = N∗ ( N N∗ − 1− ln N N∗ ) (26) We shall use the family of Volterra-type Lyapunov function defined by g(x) = x−1−ln(x), x ∈ R+ which admits a global mimnimun at x = 1, and satisfies g(1) = 0. Since S1(t) > 0, S2(t) > 0, S3(t) > 0, EI(t) > 0, EC(t) > 0, I(t) > 0, C(t) > 0, R(t) > 0 , then one can obtain the following V = N −N∗ −N∗ln N N∗ ≥ 0 (27) Therefore, the Lyapunov function V derivative is given by the following sense dV dt = ( 1− N∗ N ) dN dt (28) Note that according to system (1) dN dt = Λ− µII − µCC − µN (29) As at the endemic equilibrium point dN dt = 0, then one obtains Λ = µII ∗ + µCC ∗ + µN∗ (30) From (28), (29), (30), and by assuming that N −N∗ ≥ 0, I − I∗ ≥ 0, C − C∗ ≥ 0, one has dV dt = ( 1− N∗ N ) (µII ∗ + µCC ∗ + µN∗ − µII − µCC − µN) dV dt = − ( N −N∗ N ) [µa (I − I∗) + µC (C − C∗) + µ (N −N∗)] dV dt ≤ 0 (31) W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 81 From (31) and by using the fact that dV dt = 0 if and only if S1 = S∗ 1 , S2 = S∗ 2 , S3 = S∗ 3 , EI = E∗ I , EC = E∗ C , I = I∗, C = C∗, R = R∗, then dV dt converges when t −→ ∞ . Thanks to LaSalle’s invariance principle theorem [26], the endemic equilibrium point E∗ is said to be globally asymptotically stable when R0 > 1 [4, 14, 18, 27] W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 82 6. Numerical results and discussion For the numerical simulation of the proposed system, the Rung Kutta 4 method was used. In table 3 the values of the different parameters used in the model problem (1) are based on data from Burkina Faso and the literature. With these parameters, the basic reproduction number R0 = 3.6625. We set N0 = 20818036 : total number of Burkina Faso’s population in 2020 [20]. Initial conditions were estimated using prevalence in 2020 and biological parameters. This gives S1(0) = 1191200, S2(0) = 2384100, S3(0) = 11938000, EI(0) = 1551400, EC(0) = 172604, I(0) = 1665482, C(0) = 416354, R(0) = 1498896. Table 3: Parameters used in the model Parameters Values References N0 20818036 [20] λ 0.0444 [20] λi 0.091 Assumed λc 0.0051 Assumed p1 0.9 Assumed p2 0.9 Assumed b1,j , j = 1, 2 0 [16] b2,j , j = 1, 2 0 [16] b3,j , j = 1, 2 0 [16] b1,j , j = 3, 4 0.159 [[24],[15]] b2,j , i = 3, 4 0.144 [[24],[15]] b3,j , j = 3, 4 0.5 Assumed α1 0.1 [[19],[17] ] α2 0.7 [[19],[17] ] α3 0.95 [[19],[17] ] δ1 0.6 [[24],[15]] δ2 0.6 [[24],[15]] γi 0.8 [[24],[15]] γc 0.023 [[24],[15]] µ 0.009 [20] µi 0.00461 [[24],[15]] µc 0.01 [[24],[15]] We carried out simulations to see the dynamics of the hepatitis B virus in the popu- lation of Burkina Faso over a period of 50 years. Figure 4 represents the evolution of the three susceptible compartments, figure 5 repre- sents the evolution of infected and chronic exposed compartments, while figure 6 is the one of infected and chronic individuals. The last figure 7 is that of recovered individuals. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 83 0 5 10 15 20 25 30 35 40 45 50 Time[y] 0 5 10 15 S u s c e p ti b le 10 6 Population of susceptibles S1 S2 S3 Figure 4: Evolution of susceptible individuals In Figure 4, the number of all susceptible individuals is decreasing, confirming the high prevalence of HBV infection in Burkina Faso. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 84 0 5 10 15 20 25 30 35 40 45 50 Time[y] 0 2 4 6 8 10 12 14 16 P o p u la ti o n o f E x p o s e d 105 Population of Exposed Exposed Acutes Exposed Chronicles Figure 5: Evolution of acutes exposed together with chronic exposed individuals According to the results presented in this figure, the number of cases in the acute and chronic infected exposed compartments (Figure 5) increases as they come from susceptible individuals, before decreasing to become acute and chronically infected. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 85 0 5 10 15 20 25 30 35 40 45 50 Time[y] 0 0.5 1 1.5 2 2.5 3 3.5 P o p u la ti o n o f A c u te s a n d C h ro n ic le s 10 6 Population of Acutes and Chronicles Acutes Chronicles Figure 6: Evolution of acutes and chronics individuals As shown in this figure, the number of people in the acutes and chronics compartments (6) increases, since we are in a zone of high endemicity and this confirms the high preva- lence of heaptitis B in Burkina Faso [31]. The number of chronic carriers of the virus is rising sharply, while the number of acutely infected individuals is falling, as they either recover or become chronic. W.-O Sawadogo et al. / Eur. J. Pure Appl. Math, 17 (1) (2024), 59-92 86 0 5 10 15 20 25 30 35 40 45 50 Time[y] 0 0.5 1 1.5 2 2.5 3 3.5 re c o v e re d 10 7 Population of recovered Recovered Figure 7: Evolution of recovered individuals Here, the number of cured or recovered patients increases significantly (7), as 90−95% of acutely infected patients recover, and 30% of those (5 − 10%) who become chroni- cally infected are healthy carriers [35], and therefore recovered from the hepatitis B virus transmission process. 7. Conclusion In this paper, a new differential mathematical model of susceptibility and infectivity of hepatitis B transmission was developed to simulate the potential spread of hepatitis B virus in the population of Burkina Faso. The global existence and uniqueness of the model’s solutions have been proven. Mathematical analysis of the model shows that dis- ease progression is governed by the basic reproduction number R0, which is a key concept in epidemiology. When R0 < 1, the disease disappears from the population, resulting in a disease-free equilibrium that has been proven to be globally asymptotically stable. On the other hand, when R0 > 1 , the disease persists, leading to an endemic state that is also globally asymptotically stable. The numerical simulation of the model used to predict the spread of infection in the population showed that : REFERENCES 87 • R0 = 3.6625 > 1, which proves that the disease is persisting in Burkina Faso; • the number of susceptible people falls sharply as they become infected due to the persistence of the disease; • the number of chronic hepatitis B virus carriers increases gradually over the fifty-year simulation period. These numerical results confirm the high endemicity of hepatitis B viral infection in Bur- kina Faso, which is why this pathology remains a public health problem for the country. In perspective, we intend to improve this model by taking into account vaccination, treatment and awareness, in order to assess the impact of the strategies implemented in Burkina Faso to combat the spread of the virus. Acknowledgements The authors thank World Bank for its financial support. References [1] H. Alaa, N.-E. Alaa, and F. Aqel. Development and simulation of a mathematical model to simulate the phase transmissibility of covid19 in morocco. Annals of the University of Craiova-Mathematics and Computer Science Series, 49(1):75–83, 2022. [2] R.-M. Anderson and R.-M. May. Infectious diseases of humans: dynamics and control. 1992. [3] F. Aqel, H. Alaa, and N.-E. Alaa. Mathematical model of covid-19 transmissibility during the vaccination period. Eurasian Journal of Mathematical and Computer Applications, 11(1):4–28, 2023. [4] C. Bounkaicha, K. Allali, Y. Tabit, and Danane J. Global dynamic of spatio-temporal fractional order seir model. MATHEMATICAL MODELING AND COMPUTING, 10(2):299–310, 2023. [5] C. Castillo-Chavez, S. Blower, P. van den Driessche, D. Kirschner, and A.-A. Yakubu. Mathematical approaches for emerging and reemerging infectious diseases: models, methods, and theory, volume 126. Springer Science & Business Media, 2002. [6] C.-C. Chavez, Z. Feng, and W. Huang. On the computation of r0 and its role on global stability. Mathematical approaches for emerging and re-emerging infection diseases: an introduction, 125:31–65, 2002. [7] D. Debray. Hepatitis b and c in children: natural history and treatment, pediatric hepatology, necker hospital. January 2015. REFERENCES 88 [8] O. Diekmann and J.-A.-P. Heesterbeek. Mathematical epidemiology of infectious diseases: model building, analysis and interpretation. John Wiley & Sons, 5:262, 2000. [9] A. Din. The stochastic bifurcation analysis and stochastic delayed optimal control for epidemic model with general incidence function. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(12):123101, 2021. [10] A. Din and Y. Li. Stationary distribution extinction and optimal control for the stochastic hepatitis b epidemic model with partial immunity. Physica Scripta, 96(7):074005, 2021. [11] A. Din and Y. Li. Stochastic optimal control for norovirus transmission dynamics by contaminated food and water. Chinese Physics B, 31(2):020202, 2022. [12] A. Din, Y. Li, F. Muhammad Khan, Z.-U. Khan, and P. Liu. On analysis of frac- tional order mathematical model of hepatitis b using atangana–baleanu caputo (abc) derivative. Fractals, 30(01):2240017, 2022. [13] A. Din, Y. Li, A. Yusuf, and A.-I. Ali. Caputo type fractional operator applied to hepatitis b system. Fractals, 30(01):2240023, 2022. [14] X. Duan, S. Yuan, and X. Li. Global stability of an svir model with age of vaccination. Applied Mathematics and Computation, 226:528–540, 2014. [15] W.-J. Edmunds, G.-F. Medley, D.-J. Nokes, C.-J. Ocallaghan, H.-C. Whittle, and A.- J. Hall. Epidemiological patterns of hepatitis b virus (hbv) in highly endemic areas. Epidemiol Infect., 117(2):313–325, 1996. [16] A. A. Fall. Studies of some epidemiological models: application to the transmission of the hepatitis b virus in sub-saharan africa (case of senegal), 2013. [17] J. Feld, H. L. Janssen, Z. Abbas, A. Elewaut, P. Ferenci, V. Isakov, ..., and A. LeMair. World gastroenterology organisation global guideline hepatitis b: September 2015. Journal of Clinical Gastroenterology, 50(9):691–703, 2016. [18] H. Guo and M.-Y. Li. Global dynamics of a staged-progression model with amelio- ration for infectious diseases. Journal of Biological Dynamics, 2(2):154–168, 2008. [19] Lancet Gastroenterol Hepatol. Global, regional, and national burden of hepatitis b, 1990–2019: a systematic analysis for the global burden of disease study 2019. https://doi.org/10.1016/S2468-1253(22)00124-8, pages 796–829, 2022. [20] Burkina Faso Institut National de Statistique et de Démographie. Fifth general census of population and housing of burkina faso. synthesis of final results. https://insd.bf/, 2019. REFERENCES 89 [21] D.-G. Kendall. Deterministic and stochastic epidemics in closed populations. In Pro- ceedings of the third Berkeley symposium on mathematical statistics and probability, volume 4, pages 149–165. University of California Press Berkeley, 1956. [22] W.-O. Kermack and A.-G. McKendrick. A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character, 115(772):700–721, 1927. [23] A. Korobeinikov. Lyapunov functions and global stability for sir and sirs epidemio- logical models with non-linear transmission. Bull. Math. Biol., 68(3):615–626, 2006. [24] J.-P.-I. Kouenkam, J. Mbang, and Y. Emvudu. Global dynamics of a model of hep- atitis b virus infection in a sub-saharan african rural area. International Journal of Biomathematics, 13(6):313–325, 2020. [25] O. Koutou and B. Sangaré. Mathematical analysis of the impact of the media coverage in mitigating the outbreak of covid-19. Mathematics and Computers in Simulation, 205:600–618, 2023. [26] J.-P. La Salle. The stability of dynamical systems. SIAM, Philadelphia, 1976. [27] V. Lakshmikantham, S. Leela, and A.-A. Martynyuk. Stability analysis of nonlinear systems. Marcel Dekker, New York, 1989. [28] M. Magane. Incidental discovery of hbs antigen at the y.o. university hospital: socio- demographic, clinical and paraclinical aspects, 2018. [29] J. Martinet. Cellules dendritiques plasmocytöıdes et infection par le virus de l’hépatite b, rôle physiopathologique et potentiel vaccinal, 2012. [30] C.-C. McCluskey. Global stability of an sir epidemic model with delay and general nonlinear incidence. Math. Biosci. Eng.l, 7(4):837–850, 2010. [31] N. Meda, E. Tuaillon, D. Kania, A. Tiendrebeogo, A. Pisoni, S. Zida, ..., and P. Dujols. Hepatitis b and c virus seroprevalence, burkina faso: a cross-sectional study. Bulletin of the World Health Organization, 96(11):750, 2018. [32] Ministry of health. Burkina: hepatitis b, cirrhosis and primary liver cancer kill nearly 5,000 people a year by editorial board. https://fasozine.com, 2017. [33] Ministry of Health. normes et protocoles de prise en charge des hepatites virales au burkina faso. 2019. [34] J.-J. Ott, G.-A. Stevens, J. Groeger, and S.-T. Wiersma. Global epidemiology of hepatitis b virus infection: new estimates of age-specific hbsag seroprevalence and endemicity. Vaccine, 30(12):2212–2219, 2012. REFERENCES 90 [35] Stanislas POl. Natural history of infection hepatitis b virus infection. Presse Med, online www.masson.fr/revues/pm, 35:308–316, 2006. [36] R. Ross. An application of the theory of probabilities to the study of a priori path- ometry.—part i. Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character, 92(638):204–230, 1916. [37] R. Ross and H.-P. Hudson. An application of the theory of probabilities to the study of a priori pathometry : Part iii. Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character, 89(621):225–240, 1917. [38] R. Ross and H.-P. Hudson. An application of the theory of probabilities to the study of a priori pathometry: Part ii. Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character, 93(650):212–225, 1917. [39] I. Tao, T. Compaoré, B. Diarra, F. Djigma, T.-M. Zohoncon, M. Assih, D. Ouermi, V. Pietra, S.-D. Karou, and J. Simpore. Seroepidemiology of hepatitis b and c viruses in the general population of burkina faso. Hepatitis research and treatment, 2014, 2014. [40] P. Van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences, 180(1-2):29–48, 2002. [41] I. Zada, M. Naeem Jan, N. Ali, D. Alrowail, K. Sooppy Nisar, and G. Zaman. Math- ematical analysis of hepatitis b epidemic model with optimal control. Advances in Difference Equations, 2021(1):1–29, 2021. A. Appendix % hepatis.m function dtdy = hepatis(~,u) global N l l1 l2 r mu mu1 mu2 a1 a2 a3 p1 p2 global b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 global d1 d2 g2 g3 dtdy = zeros (8,1); dS1 = l -l1*u(6)-l2*u(7) -(mu+p1)*u(1) -b11*u(1)*u(4) -b12*u(5)*u(1)-b13*u(6)*u(1)-b14*u(7)*u(1); dS2 = p1*u(1) -(mu+p2)*u(2) -b21*u(4)*u(2)- b22*u(5)* u(2)-b23*u(6)*u(2) -b24*u(7)*u(2); dS3 = p2*u(2) -b31*u(4)*u(3) -b32*u(3)*u(5) -b33*u (3)*u(6)-b34*u(3)*u(7) -(mu)*u(3); REFERENCES 91 dEI = a1*b11*u(1)*u(4) +a1*b12*u(1)*u(5) +a1*b13*u(1) *u(6)+a1*b14*u(1)*u(7) ... + a2*b21*u(2)*u(4)+a2*b22*u(2)*u(5)+a2*b23*u(2)*u(6)+ a2*b24*u(2)*u(7) ... + a3*b31*u(3)*u(4)+a3*b32*u(3)*u(5)+a3*b33*u(3)*u(6)+ a3*b34*u(3)*u(7) -(mu+d1)*u(4) + r*(l1*u(6)+l2*u(7) ); dEC = (1-a1)*b11*u(1)*u(4) +(1-a1)*b12*u(1)*u(5) +(1- a1)*b13*u(1)*u(6)+(1-a1)*b14*u(1)*u(7) ... + (1-a2)*b21*u(2)*u(4)+(1-a2)*b22*u(2)*u(5)+(1-a2)* b23*u(2)*u(6)+(1-a2)*b24*u(2)*u(7) ... + (1-a3)*b31*u(3)*u(4)+(1-a3)*b32*u(3)*u(5)+(1-a3)* b33*u(3)*u(6)+(1-a3)*b34*u(3)*u(7) +(1-r)*(l1*u(6) + l2*u(7))- (mu+d2)*u(5); dI =d1*u(4) - (mu+g2+mu1)*u(6); dC = d2*u(5) -(mu+g3+mu2)*u(7); dR =g2*u(6) +g3*u(7)- (mu)*u(8); dtdy=[dS1;dS2;dS3;dEI;dEC;dI;dC;dR]; % main.m global N l l1 l2 r mu mu1 mu2 mu3 a1 a2 a3 p1 p2 global b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 global d1 d2 g2 g3 N=20818036; l=0.0444; l1= 0.091; l2 =0.0051; r=0.1; [mu , mu1 , mu2 , mu3]= deal (0.009 , 0.00461 , 0.01, 0.005); [a1 ,a2 ,a3]= deal (0.1, 0.7, 0.95); [p1 , p2]= deal (0.9 , 0.9); [b11 , b12 , b13 , b14]=deal (0.0, 0.0 , 0.159 , 0.159) ; [b21 , b22 , b23 , b24]=deal (0.0, 0.0 , 0.144 , 0.144) ; [b31 , b32 , b33 , b34]= deal (0.0, 0.0 , 0.5, 0.5); [d1 , d2]=deal (0.6, 0.6); [g2 , g3] =deal (0.8, 0.023); y0 =[1191200 ,2384100 ,11938000 ,1551400 ,172604 , 1665482 ,416354 ,1498896]/N; REFERENCES 92 M=200; tspan=linspace (0,50,M); [t,y1]= ode45(@hepatis ,tspan ,y0); figure (1) plot(t,N*y1(:,1),'r',t,N*y1(:,2),'g',t,N*y1(:,3),'b', 'Linewidth ', 2.5);legend('S1','S2','S3');title('\ bf Population of susceptibles '); xlabel('\bf{Time [y]}');ylabel('\bf{Susceptible}'); figure (2) plot(t,N*y1(:,4),'g',t,N*y1(:,5),'r','Linewidth ', 2.5);legend('Exposed Acutes ','Exposed Chronicles ') ;title('\bf Population of Exposed ');xlabel('\bf{ Time[y]}');ylabel('\bf{Population of Exposed}'); figure (3) plot(t,N*y1(:,6),'g',t,N*y1(:,7),'r','Linewidth ', 2.5);legend('Acutes ','Chronicles ');title('\bf Population of Acutes and Chronicles ');xlabel('\bf{ Time[y]}');ylabel('\bf{Population of Acutes and Chronicles}'); figure (4) plot(t,N*y1(:,8),'Linewidth ', 2.5);legend('Recovered ' );title('\bf Population of recovered ');xlabel('\bf {Time[y]}');ylabel('\bf{recovered}');