Estimation of Parameters for the Mathematical Model of the Spread of Hepatitis B in Burkina Faso Using Grey Wolf Optimizer

. In this paper, we developed a mathematical model of di ﬀ erential susceptibility, taking into account vaccination and treatment, to simulate the transmission of the hepatitis B virus in the population of Burkina Faso. The existence and uniqueness of non-negative solutions are proved. The model has a globally asymptotically stable disease free equilibrium when the basic reproduction number R 0 < 1 and an endemic equilibrium when R 0 > 1. We estimated the parameters of the model based on hepatitis B cases from 1997 to 2020 by using a Grey Wolf Optimizer Algorithm (GWO). The results demonstrated the e ﬃ cacy of the GWO algorithm in estimating the model parameters. A sensitivity analysis was carried out to determine the determining factors in the spread of hepatitis B in Burkina Faso. The estimated parameters were used to simulate the spread of hepatitis B in Burkina Faso from 1997 to 2020.


Introduction
Viral hepatitis B remains a major public health problem in the world, especially in Africa and particularly in Burkina Faso.HBV is second only to tobacco as a known human carcinogen.
It damages the liver through acute and chronic infection.The chronic infection, through its complications such as cirrhosis and/or hepatocellular carcinoma (liver cancer), makes it a really serious pathology since it causes more than 5000 deaths per year in Burkina Faso [43].It is well noted that one over four people with chronic HBV infection is likely to die prematurely of cirrhosis or liver cancer.Hepatitis B endemicity is defined as the seroprevalence of HBsAg carrier statusantigen in the general population.Burkina Faso is classified as a highly endemic country ( [39], [6]).In addition, in the 2011 Demographic Health Survey, HBV seroprevalence was estimated at 9.1% [5].Similarly, a study conducted in Ouagadougou (capital city of Burkina Faso) revealed a prevalence of 14.5% in the general population [7].Thus, nearly 2 million people in Burkina Faso are infected with HBV, whether they know it or not [41].In these highly endemic areas, transmission of hepatitis B virus occurs mainly by vertical transmission, that is, from infected mothers to their infants during the perinatal period, and the rate of chronic carriage is inversely proportional to the age at which the individual becomes infected ( [6,39]).The three main axes of the fight against the spread of HBV in Burkina Faso are: vaccination, treatment, and awareness.
Mathematical models are valuable tools for understanding the dynamics of hepatitis B in the population.They can be used to assess the impact of measures taken to combat the disease.
Edmunds et al. [19] proposed a mathematical model to study the transmission dynamics and control of hepatitis B virus in the Gambia.Kamyad et al. [22] Proposed a mathematical model for hepatitis b virus taking acount vaccination and treatment.Khan et al. [23] presented a model for the hepatitis B virus dynamic by dividing the infectious class into three sub-classes: acute, chronic and carrier.Horizontal and vertical transmissions have also been taken into account.
Fall [9] proposed a differential susceptibility and infectivity mathematical model of Hepatitis B transmission to study the dynamic of virus in sub-Saharan Africa (case of Senegal).
Recently, a mathematical model for analysis of effective intervention strategies on transmission dynamics of hepatitis B virus has been proposed by Firaol et al. [20].
In this paper, a differential susceptibility mathematical model of Hepatitis B ( [3], [9]) transmission taking acount vaccination and treatement was developed in order to simulate the potential spread of the Hepatitis B virus in the population of Burkina Faso presented.This depend on parameters.To identify these parameters, a cost function to minimize or maximize is built by using Burkina Faso hepatitis B cases for the period 1997 to 2020.Generally these functions have a complex structure with respect to the parameters, and sometimes they suffer from a lack of regularity.Thus, the minimization of these functions using classical methods becomes impossible, which is precisely the main problem encountered here.Hence, the remedy is to use meta-heuristic methods.Meta-heuristic optimization techniques are adapted better to the problems of optimization in which the size of the space of research is important, where the parameters interact in a complex way and where very little information on the function to be optimized is available [21].
In this work we used Grey Wolf Optimizer (GWO) algorithm, a meta-heuristic proposed by S.
Mirjalili et al [21].This algorithm is inspired by grey wolves.He mimics the leadership hierarchy and hunting mechanism of grey wolves in nature.
The remainder of this paper is organized as follows: section 2 is dedicated to the presentation of the proposed model.In Section 3, the mathematical analysis of the model is done focusing on uniqueness, existence, positivity and boundedness of the solutions, global stability of equilibriums.
We also determined the number of basic reproduction.Section 4 is devoted to the parameters estimation.After constructing the function to be minimized, we presented the algorithm used.
The estimation result was presented.Data of hepatitis B cases for the period 2007 to 2020 of Burkina Faso are used.In Section 5, discussion of the obtained numerical results is presented.The conclusion is given in section 6.

Model Formulation
Here we formulate a deterministic mathematical model with differential susceptibility [3] SVLACTR, i.e.Susceptible-Vaccinated-Latent-Acute-Chronic-Treatment-Recovered, composed of nine ordinary differential equations (ODE) to illustrate the dynamics of HBV in Burkina Faso.The cumulative population at time t, represented by N(t), is classified into nine different compartments, the meaning of each of which is given in Table 1.We draw the compartmental diagram given in Figure 1, from which we obtain the model 2.1 of hepatitis B transmission dynamics with vaccination and treatment in the Burkina Faso population.In this work, we make these considerations: (i) Susceptible individuals were subdivided into three age classes S i , i = 1, 2, 3. (ii) The high prevalence (> 8%) in Burkina Faso is mainly due to the vertical transmission, hence the choice of the vertical transmission model ( [5]), [9]).
(iv) A Part of the births from the acute and chronic HBV infections participate with an amount λ 1 A + λ 2 C + λ 3 T, which will move to the exposed compartments of both acute and chronic infections.
(v) Vertical transmission occurs in acutely infected, chronic carriers and treated.
(vi) The adequate contact β coefficient for a susceptible person to be infected with HBV, depends of infectious ( [45], [3]) individuals i.e. the A, C, T r .One has: adequate coefficient of contact so that a susceptible S i becomes contaminated by the HBV.
(viii) All individuals who do not die from HBV, i.e. individuals in compartments S 1 , S 2 , S 3 , V, L, R have the same mortality rate µ.
(ix) The risk of hepatitis B virus reactivation in recovered individuals can be high in cases of immunosuppression ( [31]), [32]). Let: and 3) The normalized model of the model 2.1 is given by: 3. Mathematical analysis 3.1.Positivity, boundedness and global existence of solutions.
(i) For the local existence, all the functions of system 2.4 are locally Lipschitz continuous.
Thus, there exists a unique local solution on t ∈ [T 0 , T max ), where T max 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 [T 0 , T max ) [2].By proving that the components of the solution vector (s 1 (t), s 2 (t), s 3 (t), v(t), l(t), a(t), c(t), t r (t), r(t)) are uniformly bounded on any bounded interval [0, T max ), one ensures that T max = ∞.We remark that the components of the vector Where are quasi-positive.
Consequently, since the initial conditions are non-negative, this implies that the solution components are non-negative for all t ∈ [T 0 , T max ).Now, let the function Ψ be defined as By taking the sum of the first nine equations in 2.4, we observe Integrating equation 3.6 over (0, t) for all t 0 < t < T, one can get the following which implies that Here, two different cases are distinguished.If λ µ < 1, then the following inequality is satisfied Finally, one can get the following which Hence, T max = ∞ and the existence of unique, non-negative and bounded global solution are proved.
(ii) We remark that s 1 satisfies the following: By integrating (11) over [0, t] for all t > 0 , we obtain: Since The same reasoning is applied for s 2 (t), s 3 (t) and v(t) which concludes the proof of this theorem.
The basic reproduction number R 0 of model is: Where k 5 , k 6 , k 7 , k 8 are given by 2.2 and H 1 , H 2 , H 3 are given by: Proof.Here, we consider the proposed mathematical model 2.4 with nine homogeneous compartments.This model can be written as (s 1 , s 2 , s 3 , v, l, a, c, t r , r) T = F(s 1 , s 2 , s 3 , v, l, a, c, tr, r) T Where F is defined by 3.5.The point E 0 defined by 3.2 satisfies F(E 0 ) = 0.
To calculate the reproduction rate R 0 , let's use the next generation matrix method described in [1].
Thus let us consider only the infected and infectious compartments which satisfy the following order 4 system: The rate of occurrence of new infection in the four compartments (l, a, c, t r ) is represented by the vector F as follows and the transfer rate of individuals into and out of the infected compartments is given by the Thus, all the epidemiological events leading to new infections are incorporated into the model via the matrix F , while all the other events are included in the matrix G. Progression to either death or immunity ensures that G is invertible.By substituting the variables (s 1 , s 2 , s 3 ) with the free equilibrium point E 0 , we obtain the two matrices F and G which are expressed as follows: Where H 1 , H 2 , 3 are given in 3.14 By calculating G −1 , we obtain: We calculate the matrix of the new generation, we obtain: Where By the next generation matrix approach, R 0 for the model.By calculating the spectral radius of the next generation matrix F G −1 , we get: Where k 5 , k 6 , k 7 , k 8 are given by 2.2 and H 1 , H 2 , H 3 are given by 3.14.

Endemic equilibria.
In presence of infected individuals, model 2.1 is said to be exhibiting an endemic equilibrium point r , R * are given by: 3.4.Global stability of disease-free equilibrium point E 0 .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 ( [15], [14], [10]) to prove the stability of system 2.4 in the global sense at the disease free equilibrium point.
Therefore, by applying the Castillo Chavez technique ( [15], [10]), the given problem 2.4 is converted into the following sub-models: Where X 1 and X 2 designate the population of uninfected individuals, and infected individuals, respectively.
The following conditions (H1) and (H2) must be satisfied to guarantee the local asymptotic stability.
then the DFE point E 0 of the model 2.4 is globally asymptotically stable.
Proof.In System 2.4, one can set the following Where each component is defined by the following sense At the DFE point, it is clear that G(X 1 , 0) = 0. To derive the condition (A 2 ), we first calculate the matrix A = D X 2 G(X 0 ) at the DFE point.
It is clear that the matrix A given above is an M-matrix.Now, one can calculate the following function: By replacing H 1 , H 2 , H 3 by their expressions given in 3.14, one obtain the following result: Thanks to Theorem 1, (ii) , one gets the folllowing result: Since Theorem3.1, the following result is obtained Ḡ(X 1 , X 2 ) ≥ 0 .Consequently, the hypothesis (H 1 ) and (H 2 ) are satisfied.Moreover, one uses Castillo Chavez technique ( [10], [15], [13]) to conclude that If R 0 < 1, the Disease free-equilibrium point is globally asymptotically stable.
3.5.Global stability of endemic equilibrium point E * .In this subsection, by constructing suitable Lyapunov function, we will prove the global asymptotic stability of the endemic equilibrium point Theorem 3.4.If R 0 > 1, the endemic equilibrium point E * of system 2.1 is globally asymptotically stable.
Proof.When R 0 > 1, one defines the following Lyapunov function as in ( [17], [12] , [16], [11], [47]): Then, the Lyapunov function can also be rewriten as follows: 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 , then one can obtain the following Therefore, the Lyapunov function V derivative is given by the following sense Note that according to system (1) As at the endemic equilibrium point dN dt = 0, then one obtains From 2.1, 2.2, and by assuming that From (2.4) and by using the fact that dV dt = 0 if and only if Thanks to LaSalle's invariance principle theorem [47], the endemic equilibrium point ↑ * is said to be globally asymptotically stable when R 0 > 1 ( [17], [46]).

Parameters estimation
4.1.Problem to solve.Let thus M observations of values of infected I obs (t j ) at the moments t j , j = 1, . . ., M.
• the omega wolves always have to submit to all the other dominant wolves.
To model mathematically the hunting mechanism of grey wolves, three steps were considered: • Tracking, chassing and approaching the prey.
• Pursuing, encircling and harassing the prey until it stops moving.
• Attacking the prey.
To respect the hierarchy, the best solution is alpha, the second solution is beta and the third delta.
The optimum being the position of the prey.The grey wolves encircle prey during the hunt.The mathematical model is given: where t indicates the current iteration, A = 2a.r 1 , C = 2. r 1 ; a are linearly decreased from 2 to 0 over the course of iterations and r 1 , r 2 are random vectors in [0, 1].X p is the position vector of the prey, and X indicates the position vector of a grey wolf.
For better exploration of candidate solutions which tend to diverge when | A| > 1 and to converge when | A| < 1.
Grey wolves have the ability to recognize the location of prey and encircle them.Over the course of iterations, the first three fittest solutions we obtain so far are considered as α, β and δ respectively, which guide the optimization process (the hunting) and are assumed to take the position of the optimum (the prey).To model this process, we adapt the positions of population using the following formula: where: • C 1 , C 2 and C 3 are random vectors.
• X α , X α and X δ , the positions of alpha, beta and delta respectively.
• X the position of prey( current solution).
The next position of the best solution is given by: where A 1 , A 2 and A 3 are random vectors.
To accelerate convergence, we take µ natural mortality rate 0.0165 [42] µ a acute HBV related mortality rate 0.005 [42] µ c acute HBV related mortality rate 0.02 [42] µ tr treated patients related mortality rate 0.008 [42] τ rate of reactivation of hbv infection after recovery 0.0002 Assumed   The figure (3) shows the observed data and the estimated data.This demonstrates the robustness of the Gray Wolf Optimizer (GWO) algorithm.The numerical simulation based on the estimated parameters provides, on one hand, the figure (4), illustrating the increasing trend in the number of acute infections, chronic carriers, and individuals treated chronically in the population of Burkina Faso.On the other hand, the obtained value for the basic reproduction number is R 0 = 2.5088 > 1.These results confirm the current trend in the endemic dynamics of hepatitis B in Burkina Faso.

Sensitivity analysis.
For hepatitis B virus transmission, the impact of each parameter on the endemic threshold was described in the sensitivity analysis.Sensitivity analysis is an essential method for complex systems, and has been used to determine the strength of model parameters 2.4.The parameter with higher sensitivity index magnitude is more influential than that with smaller magnitude of sensitive index.The sign of the sensitivity indices of R 0 with respect to the parameters shows the positive or negative impact of the parameters.Here, we calculate sensitivity index for each parameter included in effective reproduction number using the parameters given on Table 3 above.The standard equation for the sensitivity index of a parameter Υ for R 0 is given in 5.2 [37,38]: Given the complexity of the expression for R 0 , we have used numerical differentiation.Thus the numerical values of the sensitivity indices are given in table (5).

Sensitivity index value
γ c = (∂R 0 /∂γ c ) × (γ c /R 0 ) -0.066392 χ R 0 γ tr = (∂R 0 /∂γ tr ) × (γ tr /R 0 ) -0.124634 Table 5. Sensitivity indices of the model parameters Which show on the one hand that infectious contacts with acutely infected people and carriers as well as infected births from acute and chronic mothers, on the other hand, the evolution of acutely infected people towards chronicity, strongly contribute to the increase in the endemic threshold and therefore the spread of the disease.The most negative sensitive parameter is χ R 0 γ a = −0.965469with basic reproduction number R 0 = 2.5088, which proves that all strategies that can help reduce the number of acute infections in the population can prevent the spread of the disease.In addition to vaccination, which will significantly reduce infectious contacts, we recommend mass screening of the population and early treatment of acute infections.

Conclusion
In this work, a mathematical model of hepatitis B virus transmission with differential susceptibility according to three age groups, including vertical transmission from acutely, chronically, and treated infected mothers, was developed and applied to the case of Burkina Faso.This model took into account the vaccination of the three susceptible groups and the treatment of chronically ill people.Our modeling process first began with the study of some biological aspects of the disease that enabled us to better approach the modeling problem.Mathematical analysis of the model shows that the course of the disease is governed by the basic reproduction number R 0 < 1.Indeed, when R 0 the disease disappears from the population giving a disease-free equilibrium which has been proved to be and globally asymptotically stable.Moreover, when R 0 exceeds 1, the disease persists and leads to an endemic state which is also and globally asymptotically stable.
Using data on acutely, chronically, and treated infected individuals in Burkina Faso from 1997 to 2020, we successfully estimated the optimal parameters of our model by solving a global optimization problem using the Grey Wolf Optimizer algorithm.The numerical simulation which showed a strong increase in infectious individuals in the population with a basic reproduction number R 0 = 2.5088 > 1, which confirms the high endemicity of hepatitis B in Burkina Faso.
Furthermore, the sensitivity analysis of R 0 highlighted the influence of each parameter on the dynamics of hepatitis B in Burkina Faso.Thus, the most positively sensitive parameters include, on the one hand, the appropriate contact coefficients for a susceptible to be infected by an acutely ill individual β 1 or by a chronically infected carrier β 2 , and, on the other hand, the proportions λ 1 and λ 2 of births infected through vertical transmission from acutely and chronically infected mothers, respectively.
To reverse the current trend of the endemic dynamics of hepatitis B virus infection in Burkina Faso, we recommend: • The rigorous implementation of systematic vaccination for children at birth, especially in this period of security crisis accompanied by a massive displacement of populations (majority of whom are women and children), in order to break the chain of transmission.
This would significantly reduce the proportions of infected births λ 1 and λ 2 .
• Mass screening of populations, so that those who test negative for hepatitis B virus infection can be vaccinated, and those who test positive can take measures to avoid being a reservoir of infection or undergo treatment to reduce their infectivity or achieve recovery.The goal here would be to reduce the appropriate contact coefficients β 1 and β 2 .

3. 2 .Theorem 3 . 2 .
Disease free-equilibrium and Basic reproduction number.Consider the model 2.4 with the given parameters in 3.1 Then, (i) The disease free equilibrium point (DFE) is

. 7 ) 4 . 3 . 5 . Results and discussion 5 . 1 .
Problem solving algorithm.In summary the resolution of our problem by GWO algorithm is given by the algorithm below: Algorithm 4.1 Algorithm Initialize the input parameters for GWO (N p , d, lb, ub, Maxiter) Initialize Alpha, Beta and Delta Position and Score.Initialize the random position of search agents.k ← 0 while k < Maxiter do for i=1 to N do solving direct problem (2.1) evaluate the score of each search agent using objective function (4.2) if f itness < AlphaScore then Update alpha end if if f itness > AlphaScore and f itness < BetaScore then Update beta end if if f itness > AlphaScore and f itness > BetaScore and f itness < DeltaScore then N p do Update the Position of search agents including omegas using equation (4.4-4.6)Update the position of prey using equation(4.7)end for k ← k + 1 end whileReturn the position of Alpha as the fittest optimum Results of parameter estimation.The data used were obtained using the evolution of hepatitis B prevalence and population trends in Burkina Faso from 1997 to 2020[25,42].The initial values used are S 1 (0) = 236764, S 2 (0) = 1304902, S 3 (0) = 5327470, V(0) = 18258, L(0) = 1355550, A(0) = 756067, C(0) = 457206, T(0) = 126802, R(0) = 886318 .The parameters values used are given in the table below:

Figure 3 .
Figure 3.The yearly Hepatitis B cases time series in Burkina Faso from 1997 to 2020 and the proposed model best fitted curve

Figure 4 .
Figure 4. Evolution of acute, chronic and treated individuals

Table 2
Parameters description λ birth rate of total population λ 1 proportion of births infected by vertical transmission, from acute mothers.λ 2 proportion of births infected by vertical transmission, from chronic mothers.λ 3 proportion of births infected by vertical transmission, from treated mothers.θ i , i = 1, 2, 3 transfer rate of S i , i = 1, 2, 3 vaccinated individuals to the vaccinated class V φ i , i = 1, 2, 3 transfer rate of vaccinated individuals and non-immunized to the susceptible S i p 1 proportion of S 1 individuals aged 0 − 1 who grow up healthy enough to enterS 2 p 2 proportion of individuals aged 1 − 5 years who grow up healthy to enter in S 3 β 1 adequate contact coefficient for a susceptible to be infected by a acute β 2 adequate contact coefficient for a susceptible to be infected by a chronic β 3 adequate contact coefficient for a susceptible to be infected by a treated

Table
HBV model with vertical transmission, vaccination and treatment

Table 3 .
Parameters Value used in the modelThe values of the other parameters were obtained by solving the problem 4.2.The following table shows the values obtained.