MODELING AND CONTROL OF HEPATITIS B VIRUS TRANSMISSION DYNAMICS USING FRACTIONAL ORDER DIFFERENTIAL EQUATIONS

,


INTRODUCTION
The hepatitis B virus causes hepatitis B, a potentially fatal liver infection.It is a significant issue for world health.It can result in chronic liver disease, chronic infection, and high mortality rates from liver cancer and cirrhosis [1].Hepatitis B infections can only happen if the virus can get into the bloodstream and reach the liver.Once inside the liver, the virus multiplies and sends out a lot of fresh viruses into the bloodstream [2].
It has two stages of infection namely: Acute and chronic, According to the World Health Organization, the Hepatitis B virus is actively infecting more than one-third of the world's population.Additionally, more than 350 million of them have ongoing infections, and regrettably, 25 to 40 percent of them pass away from primary hepatocellular carcinoma (a form of liver cancer characterized by abnormal, dangerous growth(s) in the liver) or liver cirrhosis (scarring of the liver).[3].Hepatitis B is the tenth leading cause of death worldwide .Hepatocellular cancer (HCC) is the third most common cause of cancer death worldwide since it accounts for more than five hundred per year [4].
To improve knowledge of the pathophysiology (creation and progression) of hepatitis B infection, Long and Qi, in 2008 [5] suggested mathematical models.Based on Nowak's population dynamics model of immune response to persistent viruses, their work uses mathematical equations to describe the interaction between HBV and the immunological response to the virus.
Uninfected hepatocytes, infected hepatocytes, total host hepatocytes, free virus, and a CTL (Cytotoxic T lymphocyte utilized to kill virally infected or malignant cells) reaction are the five variables that make up the suggested model.Tahir khan et al in 2021 [6] looked into and evaluated the dynamics of hepatitis B, which has a number of different infection phases and transmission channels.The Caputo-Fabrizio operator and the idea of fractional calculus were used to fractionalize the model.The fixed point theory was used to discuss an extensive investigation of existence and uniqueness.To assist the analytical work with the aid of graphical representations, several numerical findings were made.Peijiang Liu et al in 2022 [7] formulated a five compartmental model of hepatitis B model using the fractional Calculus in the Caputo sense and detailed analysis of the model was carried out.The result obtained also showed that the fractional model is best suited to model the viral infection than the classical model.Elif Demirci in 2022 [7] presented a fractional order mathematical model to explain the spread of Hepatitis B Virus (HBV) in a non-constant population.The model included both vertical and horizontal transmission of the infection and also vaccination at birth and vaccination of the susceptible class.A frequency dependent transmission rate was used.Numerical simulations of the model are presented.The approach presented in this paper differ from those presented and references therein.We present a fractional order SV EI 1 RI 2 T (Susceptible-Exposed-Infected-Removed-Treated) model to discuss the dynamics of Hepatitis B and also show the impact of vaccination/Treatment on the population.This paper is organized as follows: A brief review of the fractional calculus is presented in section 2 with definitions.Section 3 discusses fractional order models while section 4 presents model analysis involving equilibrium points and stability.Section 5 is devoted to numerical simulations and discussion of results.Section 6 gives the concluding remarks.

FRACTIONAL ORDER CALCULUS
The concept of fractional order calculus is as old as the concept of integer order calculus.The complexity and lack of application delayed its advancement until a few decades ago.Although there has been a considerable amount of work done in simulating the dynamics of epidemic diseases, it has been limited to integer-order differential equations.Most dynamical systems based on integer order calculus have recently been changed into fractional order due to the flexibility that can be used to precisely fit the experimental data much better than integer order modeling.[8].Because the fractional derivative is a generalization of the integer-order derivative, fractional modeling is an effective approach that has been used to investigate the behavior of diseases.Most vaccination models are based on ordinary differential equations (ODEs), however we characterize the behavior of these systems using fractional order differential equations in this work.The fractional derivative is defined in various ways.Gruwald-Letnikov, Riemann-Liouville, and Caputo's fractional derivatives have been employed more frequently than others, but they are not always equal.[9].Comparing these three fractional derivatives, it is a fact that Caputo's derivative of a constant is equal to zero, which is not true for the Riemann-Liouville derivative.The main advantage of Caputo's approach is that the initial conditions for fractional differential equations with Caputo derivatives take on the same form as for integer-order differential equations.Having this in mind, we restrict our attention to the Caputo derivative of order α > 0, which is rather applicable to real world.For the purpose of this research work, we now gather some well-known definitions.[9].Definition 1.The Caputo Fractional derivative of order α of a function f :

Definition of terms
The formula for the Laplace transform of the Caputo derivative is given by Definition 3. The Fractional integral of order α of a function f : ℜ + → ℜ is given by The fractional integral of the Caputo Fractional derivative of order α of a function A two-parameter funvtion of the Mittag-Leffler type is defined by the series expansion

MODEL ASSUMPTIONS AND FORMULATION
The SV EI 1 RI 2 T model is based on the following assumptions: (1) The only way of entry into the population is through birth and the only way of exit is through death from natural causes or death from Hepatitis B-related causes.
(2) The the population mixes homogeneously.That is all individuals are equally likely to be infected by the infectious individuals in a case of contact except those who are Vaccinated or Removed.
(3) Any individual who recovers completely from the disease or who has been vaccinated receives a lifelong immunity from the disease.
(4) The proportion of people that moves from the susceptible to the removed class directly is assumed to have received the required three doses of Hepatitis B vaccine.
(5) That first dose of vaccine does not confer lifelong immunity.
(6) The treated class consists of people who are undergoing treatment to remain stable (7) The rate at which people die of the disease in the treated class is lesser than the rate at which people die of the disease in the chronically infectious class (that is δ 1 < δ 2 ).(8) The population in the treatment compartment will not recover from the illness.Taking into account the above descriptions and assumptions, the Fractional SV EI 1 RI 2 T model is described by By setting α = 1, the system of equation 3.6 can be reduced to integer order system.
With the non-negative initial condition: Invariant region and the positivity of the model solutions.
Lemma 3.1.(Generalized Mean Value Theorem [10]) and if Proof.The fractional derivative of the total human population, obtained by adding all the human equations of model 3.6, is given by  Taking the Laplace transform of (3.8) gives: Taking the inverse Laplace transform of (3.9), we have: where E α,β is the Mittag-Leffler function.But the fact that the Mittag-Leffler functions has an asymptotic behavior [9,11], it follows that: , α > 0 (3.12) Expanding (1.6), we have Expanding (1.7), we have Since Mittag-Leffler function has an asymptotic property, we have Taking limit as k−→ ∞, we have Then, it is clear that Ω is a positive invariant set.Therefore, all solutions of the model with initial conditions in Ω remain in Ω for all t > 0. Then , Ω = N(t) > 0 implies that it is feasible with respect to model (3.6).• If R 0 < 1, the the Disease free equilibrium is LAS (Locally asymptotically stable) and the disease cannot invade the population.

MODEL ANALYSIS
• If R ( 0) > 1, implies that the DFE (Diseases free equilibrium) is unstable and invasion is possible.
Diekmann et all and Van Driessche et al [12,13] provided a method for calculating R 0 , which is the formation of the next-generation matrix.It is comprised of two parts: F and V − 1, where The F i are the new infections, while the V i are transfer of infection from one compartment to another.x 0 is the disease free equilibrium point.R 0 is the spectral radius of the next generation matrix, which is the dominant Eigenvalue of the same matrix.To calculate this, we consider the infected compartments E(t), I 1 (t), I 2 (t), and T (t). Define: The Jacobian Matrices of F and V at DFE are given as Therefore, the dominant Eigenvalue of FV −1 is given as: 4.1.Equilibrium points and their stabilities.

4.1.2.
Local Stability of the Disease-free equilibrium point P 0 .In the previous section, we have seen that the basic reproduction number serves as a threshold parameter in determining the number of equilibria in system (3.6).We will show in this section that R 0 also determines the local stability of the equilibria.
Proof.We shall apply the method of linearization.The Jacobian matrix of system 3.6 at P 0 = ( Λ σ +µ , σ Λ µ(σ +µ) , 0, 0, 0, 0, 0) is given as: Next we find the characteristic equation which is given by|J 0 − λ I| = 0, where λ is the eigenvalue. Hence, and matrix reduces to: The characteristics cubic equation is given as: where Hence, a 3 is positive when R 0 < 1.By Routh-Hourwitz criterion [1], the Eigenvalues have negative real parts.Therefore, the disease-free equilibrium is locally asymptotically stable if 4.1.3.The Global Stability of the Disease-Free Equilibrium.The model equation is given as: The disease-free equilibrium of the model is , 0, 0, 0, 0, 0) Following the notation from theorem 2 [14], we have The matrix A can be written as The point P 0 = ( Λ σ +µ , σ Λ µ(σ +µ) , 0) is globally asymptotically stable for the system of uninfected individuals: satisfies the third equation.Also, satisfies the second equation.The solution of equation one is given as In addition, by Lemma 4 [14], E αα is nonnegative and so, by Theorem 2 [14], the disease-free equilibrium of the model (3.6) is globally asymptotically stable.
4.1.4.Endemic Equilibrium Point.The endemic equilibrium point 4.1.5.Local Stability of the Endemic Equilibrium when R 0 > 1.The Jacobian Matrix of the model 3.6 is given as The Jocobian Matrix of the Endemic Equilibrium point is Hence, λ 1,2 = −µ, λ 3 = −(δ 2 + µ) and the matrix reduces to Therefore, we have The characteristics equation becomes Further simplification yields If the coefficients are given as According to Routh Hurwitz's criterion [1], all the roots of the equation will be less than zero if the following conditions are met: • If all the coefficients and the constant term are greater than zero Then it follows that all the eigenvalues satisfy the condition |arg(λ 4.1.6.Global Stability of the Endemic Equilibrium Point when R 0 > 1.
Theorem 4.2.Suppose that R 0 > 1.Then the endemic equilibrium point P * is globally asymptotically stable in the interior of Ω.
Proof.To prove the global stability of P * , we use the Volterra type Lyapunov function approach [15].
The Lyapunov derivative of V along solutions of (3.6) is 2 and thus T = T * .This implies that V is negative definite with respect to P * .According to the LaSalle's invariance principle [15], if R 0 > 0, the endemic equilibrium P * is globally asymptotically stable.
4.1.7.Sensitivity Analysis.The effect of changing parameter values on the perceived usefulness of the reproduction number,R 0 , is demonstrated in this section.It is necessary to identify the important parameter, which may be a vital threshold to manage the illness.
The following are the mathematical representations of R 0 s sensitivity index towards the parameters β 1 ,β 2 ,Λ,γ,σ ,ε,ρ,µ,η,δ 1 : It may be deduced that some derivatives appear positive and that as any of the positive value parameters β 1 ,β 2 ,Λ,γ, described above is increased, the basic reproductive number,R 0 , increases.The proportionate reaction to the proportion stimulation is used to calculate the elasticity.
We have Consequently, we see that E β 1 ,E β 2 ,E Λ , and E γ are positive.This implies that raising the values of these parameters,β 1 ,β 2 ,Λ,γ, will increase the value of the basic reproduction number.R 0 .
The most frequently employed numerical method for solving fractional order initial value issues is the Adams-Bashforth-Moulton strategy.
Consider the fractional differential equation below: where G r j0 is the arbitrary real number, ν > 0 and the fractional differential operator D ν t is similar to the well-known Volterra integral equation in the Caputo sense.
In this study, we investigate the numerical solution of a fractional order SV EI 1 RI 2 T model with vaccination using the Adam's-Bashforth-Moulton predictor-corrector scheme.The algorithm is described below: Let h = T m , t n = nh, n = 0, 1, 2, ..., m.Corrector formulae: Predictor formulae: where

NUMERICAL SIMULATION AND DISCUSSION
In this section, we run rigorous numerical simulations to evaluate and validate our model system's analytical results 3.6.To achieve a numerical solution to the system 3.6, we used the mathematical software MATLAB (2018a version) and Adam's-Bashforth-Moulton predictorcorrector scheme.
We investigate numerical simulations of the model system 3.6 in the Caputo sense, using the parameters listed in Table 1.In the scenario, Table 1 is utilized for simulation.The following figures were produced to examine the behavior of the model 3.6 under various initial conditions.

CONCLUSION
In this paper we have discussed the fractional order SV EI 1 RI 2 T model with vaccination and treatment as control strtegies.Based on the data collected, we estimated the basic reproduction number.The fractional-order derivatives are typically more suitable in modeling because the option of derivative order allows one more degree of freedom, resulting in a better fit to real-time data with less inaccuracy than the integer-order model.The model shows that the Hepatitis B virus propagation is mostly determined by the population's contact rates with affected people.It has been observed that when the proportion of the population that is vaccinated increases, the spread of the virus is drastically reduced.Thus, if this is accomplished through widespread vaccination or making vaccine mandatory, the virus can be avoided.Treatment was also shown to be a control strategy for the spread of Hepaptitis B virus.Sensitivity analysis reveals that R 0 is directly proportional to the recruitment rate of susceptible individualsΛ, the rate of infection of susceptible individuals β 1 , β 2 and the rate of progression from exposed to infected individuals γ, all of which can be controlled through the effective implementation of vaccination drives.To obtain numerical solutions to the system, the Adam-Bashforth-Moulton predictor-corrector technique was applied.To validate the efficacy and influence of the control parameters, numerical simulations using MATLAB are given.

β 2 1 γ 1 ερ 1 ε( 1 1 δ 1 1 η 1 δ 2 1 µ
Interaction rate between the Susceptible and the Chronically Infected population day −Progression rate from exposed class to Acutely infected class day −Proportion of Acutely infected population day −− ρ) Proportion moving from the Acutely infected to the Chronically Infected day −Death rate as a result of the Infection in the Chronically infected class day −Progression rate from the chronically Infected to the Treatment Class day −Death rate as a result of the Infection in the treatment infected class day −Natural death rate day −1

4.0. 1 .
The basic reproduction number R 0 .The basic reproduction number R O is used in the study of disease transmission and control (epidemiology) to describe the average number of secondary infections caused by the introduction of one infectious person into a totally susceptible population.It has the following implications:

FIGURE 2 . 8 Figure 2 Figure 3
FIGURE 2. Comparison of dynamical behaviour of all individuals with respect to time for fractional order α = 0.82, σ = 0 and σ = 0.8

TABLE 1 .
Description of parameters and variables for model (3.6) β 1Interaction rate between the Susceptible and the Acutely Infected population day −1

TABLE 2 .
Estimated values of parameters and Varibles