Analytic solution of a fractional-order hepatitis model using Laplace Adomian decomposition method and optimal control analysis

Nnaemeka S. Aguegboh; Phineas Roy Kiogora; Mutua Felix; Walter Okongo; Boubacar Diallo

doi:10.1515/cmb-2023-0114

Article Open Access

Analytic solution of a fractional-order hepatitis model using Laplace Adomian decomposition method and optimal control analysis

Nnaemeka S. Aguegboh , Phineas Roy Kiogora , Mutua Felix , Walter Okongo and Boubacar Diallo

Published/Copyright: July 29, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 12 Issue 1

Abstract

Infectious illnesses like hepatitis place a heavy cost on global health, and precise mathematical models must be created in order to understand and manage them. The Adomian decomposition method (ADM) and an optimal control strategy are utilized to solve a fractional-order hepatitis model in this research. By adding fractional derivatives to account for memory effects and non-integer order dynamics, the fractional-order model expands the conventional compartmental model to take into account the complexity of hepatitis dynamics. The fractional-order hepatitis model is resolved using the ADM, a powerful and effective analytical approach. This approach offers a series solution that converges quickly, enabling the model’s precise analytical solution to be derived. To identify crucial criteria and enhance control mechanisms for the management of hepatitis, an optimal solution strategy is also introduced. The optimization procedure tries to lessen the disease’s spread and its negative effects on public health. We can find the best interventions, immunization schedules, and treatment regimens to effectively reduce the hepatitis pandemic by integrating the ADM solution with an optimization framework. The findings of this study show that the suggested method may be used to solve the fractional-order hepatitis model and optimize control measures. The analytical solution produced by ADM offers important insights into the underlying dynamics of hepatitis transmission, and the optimization process produces suggestions that public health professionals and politicians may put into practice. In the end, this research presents a promising direction for improving disease control efforts in a fractional-order context and contributes to a deeper understanding of hepatitis epidemiology. The importance of this method is that it gives solutions that coincide with that obtained using the numerical approach.

Keywords: hepatitis B; fractional-order differential equations; Laplace Adomian decomposition; optimal control; numerical simulations; predictor-corrector method

MSC 2010: 37A50; 60H10; 60J65; 60J70

1 Introduction

Over the past few years, mathematical modeling of infectious diseases using differential equations of integer order has received a lot of interest [1,2]. However, fractional derivatives and integrals have been successfully used to develop and analyze epidemiological models as well as other models in science and engineering [3,4]. One of these models that has attracted researchers the most in recent years is the hepatitis B model. More than one-third of the world’s population is reportedly actively infected with the hepatitis B virus (HBV), according to the World Health Organization. More than 350 million of them also have persistent infections, and sadly, 25–40% of them die from primary hepatocellular carcinoma (a type of liver cancer characterized by harmful growth(s) in the liver) or liver cirrhosis (liver scarring) [5]. The tenth leading cause of death in the world is hepatitis B. More than 500 people worldwide die from hepatocellular cancer each year, making it the third most frequent type of cancer [6]. The main mechanism of HBV transmission is a crucial variable in determining the incidence of chronic HBV in a certain location. For instance, intravenous drug use and unprotected sexual activity are the main means of transmission in regions with low prevalence, such as the United States and Western Europe. The disease is primarily distributed among youngsters in regions with moderate incidence, such as Japan and Russia, where 2–7% of the population has a chronic infection. Transmission occurs most frequently during childbirth in regions with high prevalence, such as Africa and China. In locations with high endemicity, at least 8% of the populace is chronically infected [7].

Mathematical models can be used to detect the spread of the HBV and forecast its effects on the community in order to aid in its prevention. Bernoulli initially proposed the concept of mathematical modeling for the spread of disease in 1766, which formed the groundwork for contemporary epidemiology. Ross and Hamer also introduced their infectious illness models at the turn of the twentieth century, and they made use of the law of mass action to describe how epidemic models behaved. The Reed-Frost epidemic model, which describes the connection between susceptible, infected, and immune individuals in a society, was later developed by Reed and Frost. A fractional-order model for an HBV infection with cured infected cells was taken into consideration by Salman and Yousef in 2017 [8]. The equilibria’s local asymptotic stability was examined. The peak of the infection is reduced when a fractional-order derivative is added to the model, but the disease will take longer to completely vanish, according to their numerical findings. Ain et al. [9] used an optimal variational iteration method, which was applied to solving parametric boundary vary problem, and a convergent solution was obtained. Ain and Chu [10] presented a fractional-order hepatitis model where two opertors were compared and the effect of vaccination was studied. Based on the traits of HBV transmission, Zhao et al. [11], created a five-compartmental mathematical model. A collection of first-order partial differential equations were used to represent the model. Susceptible, latently infected, temporary HBV carriers, chronic HBV carriers, and immune are the compartments. The infectious classes were located in the third and fourth compartments. Their model made the assumption that the immunological status lasts a lifetime and that the birth rate is stable. Since the likelihood of intrauterine HBV infection (infection that occurs in the womb) is extremely low, in their work, their model was simplified by disregarding it. The model presented in [11] was then modified in [12] with the incorporation of the treatment class. The modified model is given as

(1.1) D α S ( t ) = Λ − β 1 S I 1 − β 2 S I 2 − σ S − μ S , D α V ( t ) = σ S − μ V , D α E ( t ) = β 1 S I 1 + β 2 S I 2 − γ E − μ E , D α I 1 ( t ) = γ E − ε ρ I 1 − ε ( 1 − ρ ) I 1 − μ I 1 , D α R ( t ) = ε ρ I 1 − μ R , D α I 2 ( t ) = ε ( 1 − ρ ) I 1 − η I 2 − δ 1 I 2 − μ I 2 , D α T ( t ) = η I 2 − δ 2 T − μ T ,

With the non-negative initial condition:

S ( 0 ) = S 0 = n 1 , V ( 0 ) = V 0 = n 2 , E ( 0 ) = E 0 = n 3 , I 1 ( 0 ) = I 1 0 = n 4 , I 2 ( 0 ) = I 2 0 = n 5 , R ( 0 ) = R 0 = n 6 , and T ( 0 ) = T 0 = n 7 .

Table 1 shows the model variables and parameters descriptions.

Table 1

Description of parameters and variables for Model (1.1)

Variables	Description	Unit
S	Susceptible human population	People
V	Vaccinated human population	People
E	Exposed human population	People
I 1	Acutely infected human population	People
R	Recovered human population	People
I 2	Chronically infected human population	People
T	Human population under treatment	People

Parameters	Description	Unit
Λ	Recruitment rate	day ‒ 1
σ	Vaccination rate	day ‒ 1
β 1	Interaction rate between the susceptible and the acutely infected population	day ‒ 1
β 2	Interaction rate between the susceptible and the chronically infected population	day ‒ 1
γ	Progression rate from exposed class to acutely infected class	day ‒ 1
ε ρ	Proportion of acutely infected population	day ‒ 1
ε ( 1 ‒ ρ )	Proportion moving from the acutely infected to the chronically infected	day ‒ 1
δ 1	Death rate as a result of the infection in the chronically infected class	day ‒ 1
η	Progression rate from the chronically infected to the treatment class	day ‒ 1
δ 2	Death rate as a result of the infection in the treatment infected class	day ‒ 1
μ	Natural death rate	day ‒ 1

We present an analytical solution to a model utilizing the Laplace Adomian decomposition approach in this work, which was inspired by the aforementioned literature. We also do an optimal control analysis.

Adomian created the ADM in 1980, which is a useful technique for identifying a numerical and explicit solution to a set of differential equations that describe physical issues. The Laplace transform method is a potent tool in applied mathematics and engineering. The Laplace Adomain decomposition method (LADM), which is a potent technique, is produced by the combination of ADM with the Laplace transform. We convert differential equations to algebraic equations with the use of the Laplace transform, and the nonlinear terms are then broken down into Adomain polynomials. This numerical technique effectively solves a set of stochastic differential equations as well as deterministic differential equations. More specifically, it can be applied to a system of classical as well as nonclassical linear and nonlinear ordinary and partial differential equations. This approach does not call for any disruption or liberalization. Additionally, unlike fourth-order Runge Kutta method, it does not require pre-defined step sizes. Additionally, unlike the homotopy perturbation technique and homotopy analysis method, this method is not dependent on a particular parameter. However, the solutions acquired using this method are identical to those obtained using ADM [12].

This article is organized as follows: important definitions as related to the work are presented in Section 2. ADM is presented in Section 3. Section 4 is devoted to optimal control analysis of the model. Graphical plots of the obtained series solution are presented in Section 5 and the convergence analysis in Section 6. Section 7 gives the concluding remarks.

2 Preliminaries

The definitions of certain basic fractional calculus operations are covered in this section. For a more in-depth analysis, please see

Definition 1

The Caputo fractional derivative of order α of a function f : ℜ + → ℜ is given by

(2.1) D t α f ( t ) = 1 Γ ( α − n ) ∫ α t f ( n ) ( τ ) d τ ( t − τ ) α + 1 − n ( n − 1 < α ≤ n ) .

Definition 2

The formula for the Laplace transform of the Caputo derivative is given by

(2.2) ∫ 0 ∞ e − p t { D t α f ( t ) } d t = p α F ( p ) − ∑ k = 0 n − 1 p α − k − 1 f ( k ) ( 0 ) , ( n − 1 < α ≤ n ) .

Definition 3

The fractional integral of order α of a function f : ℜ + → ℜ is given by

(2.3) J α ( f ( x ) ) = 1 Γ ( α ) ∫ 0 x ( x − t ) α − 1 f ( t ) d t , α > 0 , x > 0 .

Definition 4

The fractional integral of the Caputo fractional derivative of order α of a function f : ℜ + → ℜ is given by

(2.4) J α { D α f ( t ) } = f ( t ) − ∑ k = 0 n − 1 f ( k ) ( 0 ) t k k ! , t > 0 .

Definition 5

A two-parameter function of the Mittag-Leffler type is defined by the series expansion

(2.5) E α , β ( z ) = ∑ k = 0 ∞ z k Γ ( α k + β ) , ( α , β > 0 ) .

3 LADM

The overall operation of Model (1.1) with specified initial conditions is the focus of this section. When we apply the Laplace transform to the model’s two sides (1.1), we have

(3.1) ℒ { D α 1 S ( t ) } = ℒ { Λ − β 1 S I 1 − β 2 S I 2 − σ S − μ S } , ℒ { D α 2 V ( t ) } = ℒ { σ S − μ V } , ℒ { D α 3 E ( t ) } = ℒ { β 1 S I 1 + β 2 S I 2 − γ E − μ E } , ℒ { D α 4 I 1 ( t ) } = ℒ { γ E − ε ρ I 1 − ε ( 1 − ρ ) I 1 − μ I 1 } , ℒ { D α 5 R ( t ) } = ℒ { ε ρ I 1 − μ R } , ℒ { D α 6 I 2 ( t ) } = ℒ { ε ( 1 − ρ ) I 1 − η I 2 − δ 1 I 2 − μ I 2 } , ℒ { D α 7 T ( t ) } = ℒ { η I 2 − δ 2 T − μ T } ,

which implies that

(3.2) s α 1 ℒ { D α 1 S ( t ) } − s α 1 − 1 S ( 0 ) = ℒ { Λ − β 1 S I 1 − β 2 S I 2 − σ S − μ S } , s α 1 ℒ { D α 2 V ( t ) } − s α 1 − 1 V ( 0 ) = ℒ { σ S − μ V } , s α 1 ℒ { D α 3 E ( t ) } − s α 1 − 1 E ( 0 ) = ℒ { β 1 S I 1 + β 2 S I 2 − γ E − μ E } , s α 1 ℒ { D α 4 I 1 ( t ) } − s α 1 − 1 I 1 ( 0 ) = ℒ { γ E − ε ρ I 1 − ε ( 1 − ρ ) I 1 − μ I 1 } , s α 1 ℒ { D α 5 R ( t ) } − s α 1 − 1 R ( 0 ) = ℒ { ε ρ I 1 − μ R } , s α 1 ℒ { D α 6 I 2 ( t ) } − s α 1 − 1 I 2 ( 0 ) = ℒ { ε ( 1 − ρ ) I 1 − η I 2 − δ 1 I 2 − μ I 2 } , s α 1 ℒ { D α 7 T ( t ) } − s α 1 − 1 T ( 0 ) = ℒ { η I 2 − δ 2 T − μ T } .

Now, using the initial conditions and taking inverse Laplace transform of System (3.2), we obtain

(3.3) S ( t ) = S 0 + ℒ − 1 1 s α 1 ℒ { Λ − β 1 S I 1 − β 2 S I 2 − σ S − μ S } , V ( t ) = V 0 + ℒ − 1 1 s α 2 ℒ { σ S − μ V } , E ( t ) = E 0 + ℒ − 1 1 s α 3 ℒ { β 1 S I 1 + β 2 S I 2 − γ E − μ E } , I 1 ( t ) = I 1 0 + ℒ − 1 1 s α 4 ℒ { γ E − ε ρ I 1 − ε ( 1 − ρ ) I 1 − μ I 1 } , R ( t ) = R 0 + ℒ − 1 1 s α 5 ℒ { ε ρ I 1 − μ R } , I 2 ( t ) = I 2 0 + ℒ − 1 1 s α 6 ℒ { ε ( 1 − ρ ) I 1 − η I 2 − δ 1 I 2 − μ I 2 } , T ( t ) = T 0 + ℒ − 1 1 s α 7 ℒ { η I 2 − δ 2 T − μ T } .

Using the values of initial condition in (1.1), we obtain

(3.4) S ( t ) = n 1 + ℒ − 1 1 s α 1 ℒ { Λ − β 1 S I 1 − β 2 S I 2 − σ S − μ S } , V ( t ) = n 2 + ℒ − 1 1 s α 2 ℒ { σ S − μ V } , E ( t ) = n 3 + ℒ − 1 1 s α 3 ℒ { β 1 S I 1 + β 2 S I 2 − γ E − μ E } , I 1 ( t ) = n 4 + ℒ − 1 1 s α 4 ℒ { γ E − ε ρ I 1 − ε ( 1 − ρ ) I 1 − μ I 1 } , R ( t ) = n 5 + ℒ − 1 1 s α 5 ℒ { ε ρ I 1 − μ R } , I 2 ( t ) = n 6 + ℒ − 1 1 s α 6 ℒ { ε ( 1 − ρ ) I 1 − η I 2 − δ 1 I 2 − μ I 2 } , T ( t ) = n 7 + ℒ − 1 1 s α 7 ℒ { η I 2 − δ 2 T − μ T } .

Assuming that the solutions, S ( t ) , V ( t ) , E ( t ) , I 1 ( t ) , R ( t ) , I 2 ( t ) , and T ( t ) , in the form of infinite series given by

(3.5) S ( t ) = ∑ i = 0 ∞ S n , V ( t ) = ∑ i = 0 ∞ V n , E ( t ) = ∑ i = 0 ∞ E n , I 1 ( t ) = ∑ i = 0 ∞ I 1 n , R ( t ) = ∑ i = 0 ∞ R n , I 2 ( t ) = ∑ i = 0 ∞ I 2 n , and T ( t ) = ∑ i = 0 ∞ T n .

And the non-linear terms in the model S ( t ) I 1 ( t ) and S ( t ) I 2 ( t ) are decomposed by Adomian polynomial as:

(3.6) S ( t ) I 1 ( t ) = ∑ i = 0 ∞ A n and S ( t ) I 2 ( t ) = ∑ i = 0 ∞ B n ,

where A n ad B n are Adomian polynomials given by

(3.7) A n = 1 Γ ( n + 1 ) d n d t n ∑ k = 0 n λ k S k , ∑ k = 0 n λ I 1 k λ = 0 ,

(3.8) B n = 1 Γ ( n + 1 ) d n d t n ∑ k = 0 n λ k S k , ∑ k = 0 n λ I 2 k λ = 0 .

Substituting (3.5) and (3.6) into Model (3.3), we obtain

(3.9) ℒ ( S 0 ) = n 1 s , ℒ ( V 0 ) = n 2 s , ℒ ( E 0 ) = n 3 s , ℒ ( I 1 0 ) = n 4 s , ℒ ( R 0 ) = n 5 s , ℒ ( I 2 0 ) = n 6 s , ℒ ( T 0 ) = n 7 s ℒ ( S 1 ) = ℒ { Λ − β 1 A 0 − β 2 B 0 − σ S 1 − μ S 0 } 1 s α 1 + 1 , ℒ ( V 1 ) = ℒ { σ S 0 − μ V 0 } 1 s α 2 + 1 , ℒ ( E 1 ) = ℒ { β 1 A 0 + β 2 B 0 − γ E 0 − μ E 0 } 1 s α 3 + 1 , ℒ ( I 1 1 ) = ℒ { γ E 0 − ε ρ I 1 0 − ε ( 1 − ρ ) I 1 0 − μ I 1 0 } 1 s α 4 + 1 , ℒ ( R 1 ) = ℒ { ε ρ I 1 0 − μ R 0 } 1 s α 5 + 1 , ℒ ( I 2 1 ) = ℒ { ε ( 1 − ρ ) I 1 0 − η I 2 0 − δ 1 I 2 0 − μ I 2 0 } 1 s α 6 + 1 , ℒ ( T 1 ) = ℒ { η I 2 0 − δ 2 T 0 − μ T 0 } 1 s α 7 + 1 , ℒ ( S 2 ) = ℒ { Λ − β 1 A 1 − β 2 B 1 − σ S 1 − μ S 1 } 1 s α 1 + 1 , ℒ ( V 2 ) = ℒ { σ S 1 − μ V 1 } 1 s α 2 + 1 , ℒ ( E 2 ) = ℒ { β 1 A 1 + β 2 B 1 − γ E 1 − μ E 1 } 1 s α 3 + 1 , ℒ ( I 1 2 ) = ℒ { γ E 1 − ε ρ I 1 1 − ε ( 1 − ρ ) I 1 1 − μ I 1 1 } 1 s α 4 + 1 , ℒ ( R 2 ) = ℒ { ε ρ I 1 1 − μ R 1 } 1 s α 5 + 1 , ℒ ( I 2 2 ) = ℒ { ε ( 1 − ρ ) I 1 1 − η I 2 1 − δ 1 I 2 1 − μ I 2 1 } 1 s α 6 + 1 , ℒ ( T 2 ) = ℒ { η I 2 1 − δ 2 T 1 − μ T 1 } 1 s α 7 + 1 , ⋮ ℒ ( S n + 1 ) = ℒ { Λ − β 1 A n − β 2 B n − σ S n − μ S n } 1 s α 1 + 1 , ℒ ( V n + 1 ) = ℒ { σ S n − μ V n } 1 s α 2 + 1 , ℒ ( E n + 1 ) = ℒ { β 1 A n + β 2 B n − γ E n − μ E n } 1 s α 3 + 1 , ℒ ( I 1 n + 1 ) = ℒ { γ E n − ε ρ I 1 n − ε ( 1 − ρ ) I 1 n − μ I 1 n } 1 s α 4 + 1 , ℒ ( R n + 1 ) = ℒ { ε ρ I 1 n − μ R n } 1 s α 5 + 1 , ℒ ( I 2 n + 1 ) = ℒ { ε ( 1 − ρ ) I 1 n − η I 2 n − δ 1 I 2 n − μ I 2 n } 1 s α 6 + 1 , ℒ ( T n + 1 ) = ℒ { η I 2 n − δ 2 T n − μ T n } 1 s α 7 + 1 .

Taking the inverse Laplace of (3.15), we obtain

S 0 = n 1 , V 0 = n 2 , E 0 = n 3 , I 1 0 = n 4 , R 0 = n 5 , I 2 0 = n 6 , T 0 = n 7 S 1 = ( Λ − β 1 n 1 n 4 − β 2 n 1 n 6 − σ n 1 − μ n 1 ) t α 1 Γ ( α 1 + 1 ) , V 1 = ( σ n 1 − μ n 2 ) t α 2 Γ ( α 2 + 1 ) , E 1 = ( β 1 n 1 n 4 + β 2 n 1 n 6 − γ n 3 − μ n 3 ) t α 3 Γ ( α 3 + 1 ) , I 1 1 = ( γ n 3 − ε ρ n 4 − ε ( 1 − ρ ) n 4 − μ n 4 ) t α 4 Γ ( α 4 + 1 ) , R 1 = ( ε ρ n 4 − μ n 5 ) t α 5 Γ ( α 5 + 1 ) , I 2 1 = ( ε ( 1 − ρ ) n 4 − η n 6 − δ 1 n 6 − μ n 6 ) t α 6 Γ ( α 6 + 1 ) , T 1 = ( η n 6 − δ 2 n 7 − μ n 7 ) t α 7 Γ ( α 7 + 1 ) ,

S 2 = Λ − β 1 n 1 ( γ n 3 − ε ρ n 4 − ε ( 1 − ρ ) n 4 − μ n 4 ) t α 1 + α 4 Γ ( α 1 + α 4 + 1 ) + n 4 ( Λ − β 1 n 1 n 4 − β 2 n 1 n 6 − σ n 1 − μ n 1 ) t 2 α 1 Γ ( 2 α 1 + 1 ) − β 2 n 1 ( ε ( 1 − ρ ) n 4 − η n 6 − δ 1 n 6 − μ n 6 ) t α 1 + α 6 Γ ( α 1 + α 6 + 1 ) + n 6 ( Λ − β 1 n 1 n 4 − β 2 n 1 n 6 − σ n 1 − μ n 1 ) t 2 α 1 Γ ( 2 α 1 + 1 ) − σ ( Λ − β 1 n 1 n 4 − β 2 n 1 n 6 − σ n 1 − μ n 1 ) t 2 α 1 Γ ( 2 α 1 + 1 ) − μ ( Λ − β 1 n 1 n 4 − β 2 n 1 n 6 − σ n 1 − μ n 1 ) t 2 α 1 Γ ( 2 α 1 + 1 )

V 2 = σ ( Λ − β 1 n 1 n 4 − β 2 n 1 n 6 − σ n 1 − μ n 1 ) t α 1 + α 2 Γ ( α 1 + α 2 + 1 ) − μ ( σ n 1 − μ n 2 ) t 2 α 2 Γ ( 2 α 2 + 1 ) E 2 = β 1 n 1 ( γ n 3 − ε ρ n 4 − ε ( 1 − ρ ) n 4 − μ n 4 ) t α 1 + α 4 Γ ( α 1 + α 4 + 1 ) + n 4 ( Λ − β 1 n 1 n 4 − β 2 n 1 n 6 − σ n 1 − μ n 1 ) t 2 α 1 Γ ( 2 α 1 + 1 ) + β 2 n 1 ( ε ( 1 − ρ ) n 4 − η n 6 − δ 1 n 6 − μ n 6 ) t α 1 + α 6 Γ ( α 1 + α 6 + 1 ) + n 6 ( Λ − β 1 n 1 n 4 − β 2 n 1 n 6 − σ n 1 − μ n 1 ) t 2 α 1 Γ ( 2 α 1 + 1 ) − γ ( β 1 n 1 n 4 + β 2 n 1 n 6 − γ n 3 − μ n 3 ) t 2 α 3 Γ ( 2 α 3 + 1 ) − μ ( β 1 n 1 n 4 + β 2 n 1 n 6 − γ n 3 − μ n 3 ) t 2 α 3 Γ ( 2 α 3 + 1 ) I 1 2 = γ ( β 1 n 1 n 4 + β 2 n 1 n 6 − γ n 3 − μ n 3 ) t α 3 + α 4 Γ ( α 3 + α 4 + 1 ) − ε ρ ( γ n 3 − ε ρ n 4 − ε ( 1 − ρ ) n 4 − μ n 4 ) t 2 α 4 Γ ( 2 α 4 + 1 ) − ε ( 1 − ρ ) ( γ n 3 − ε ρ n 4 − ε ( 1 − ρ ) n 4 − μ n 4 ) t 2 α 4 Γ ( 2 α 4 + 1 ) − μ ( γ n 3 − ε ρ n 4 − ε ( 1 − ρ ) n 4 − μ n 4 ) t 2 α 4 Γ ( 2 α 4 + 1 )

(3.10) R 2 = ε ρ ( γ n 3 − ε ρ n 4 − ε ( 1 − ρ ) n 4 − μ n 4 ) t α 4 + α 5 Γ ( α 4 + α 5 + 1 ) − μ ( ε ρ n 4 − μ n 5 ) t 2 α 5 Γ ( 2 α 5 + 1 ) I 2 2 = ε ( 1 − ρ ) ( γ n 3 − ε ρ n 4 − ε ( 1 − ρ ) n 4 − μ n 4 ) t α 4 + α 6 Γ ( α 4 + α 6 + 1 ) − η ( η n 6 − δ 2 n 7 − μ n 7 ) t 2 α 6 Γ ( 2 α 6 + 1 ) − δ 1 ( η n 6 − δ 2 n 7 − μ n 7 ) t 2 α 6 Γ ( 2 α 6 + 1 ) − μ ( η n 6 − δ 2 n 7 − μ n 7 ) t 2 α 6 Γ ( 2 α 6 + 1 ) T 2 = η ( ε ( 1 − ρ ) n 4 − η n 6 − δ 1 n 6 − μ n 6 ) t α 6 + α 7 Γ ( α 6 + α 7 + 1 ) − δ 2 ( η n 6 − δ 2 n 7 − μ n 7 ) t 2 α 7 Γ ( 2 α 7 + 1 ) , − μ ( η n 6 − δ 2 n 7 − μ n 7 ) t 2 α 7 Γ ( 2 α 7 + 1 )

We may calculate the remaining terms by repeating the process in the same way. In this case, the needed approximation solution of the suggested Model (1.1) only required the computation of three terms.

4 Optimal control analysis

Treatment and vaccination are important tools to fight against infectious illnesses. Hepatitis B vaccine and treatment have recently been shown to be efficient ways to control the spread of the disease. In this section, our aim is to minimize the number of hepatitis B-infected individuals and, concurrently, reduce the associated cost. This is accomplished by (i) the use of vaccinated as an effective time-dependent measure control u 1 ( t ) , and (ii) the use of treatment on the acutely infected individuals, a control measure, u 2 ( t ) . Therefore, we consider the following fractional optimal control problem:

(4.1) min ( J ) = ∫ 0 t f ( k 1 I 1 + k 2 I 2 + k 3 u 1 2 + k 4 u 2 2 ) d t

subject to

(4.2) D α S ( t ) = Λ − β 1 S I 1 − β 2 S I 2 − σ S − μ S , S ( 0 ) = S 0 > 0 , D α V ( t ) = σ S − μ V , V ( 0 ) = V 0 > 0 , D α E ( t ) = β 1 S I 1 + β 2 S I 2 − γ E − μ E , E ( 0 ) = E 0 > 0 , D α I 1 ( t ) = γ E − ε ρ I 1 − ε ( 1 − ρ ) I 1 − μ I 1 , I 1 ( 0 ) = I 1 0 > 0 , D α R ( t ) = ε ρ I 1 − μ R , R ( 0 ) = R 0 > 0 , D α I 2 ( t ) = ε ( 1 − ρ ) I 1 − η I 2 − δ 1 I 2 − μ I 2 , I 2 ( 0 ) = I 2 0 > 0 , D α T ( t ) = η I 2 − δ 2 T − μ T , T ( 0 ) = T 0 > 0 ,

where 0 ≤ u 1 ≤ 1 and 0 ≤ u 2 ≤ 1 .

The problem is resolved using the Pontryagin maximum principle for fractional optimum control [13].

Theorem 4.1

[13,14] Let u 1 ( t ) and u 2 ( t ) be the measurable control functions on [ 0 , t g ] , with u 1 ( t ) and u 2 ( t ) having a value in [ 0 , 1 ] . Then, optimal controls u 1 * and u 2 * minimizing the objective function ( J ) of (4.1) with

u 1 * = max { min { u 1 ¯ , 1 } , 0 } a n d u 2 * = max { min { u 2 ¯ , 1 } , 0 } , u 1 ¯ = ( ξ 1 − ξ 2 ) S 2 k 3 a n d u 2 ¯ = ( ξ 6 − ξ 7 ) S 2 k 4 ,

where ( S , V , E , I 1 , R , I 2 , T ) is the corresponding solution of System (4.16).

Proof

The Hamiltonian of the resulting optimal control problem is given as

(4.3) ℋ = k 1 I 1 + k 2 I 2 + k 3 u 1 2 + k 4 u 2 2 + ξ 1 ( Λ − β 1 S I 1 − β 2 S I 2 − u 1 S − μ S ) + ξ 2 ( u 1 S − μ V ) + ξ 3 ( β 1 S I 1 + β 2 S I 2 − γ E − μ E ) + ξ 4 ( γ E − ε ρ I 1 − ε ( 1 − ρ ) I 1 − μ I 1 ) + ξ 5 ( ε ρ I 1 − μ R ) + ξ 6 ( ε ( 1 − ρ ) I 1 − u 2 I 2 − δ 1 I 2 − μ I 2 ) + ξ 7 ( u 2 I 2 − δ 2 T − μ T ) ,

with ξ i ( t ) , with i = 1 , 2 , 3 , 4 , 5 , 6 , 7 being the adjoint variables with ξ i ( t g ) = 0 , expressed in terms of canonical equations:

(4.4) D α ξ 1 = − ∂ ℋ ∂ S = ξ 1 ( β 1 I 1 + β 2 I 2 + u 1 + μ ) − ξ 2 ( u 1 ) − ξ ( β 1 I 1 + β 2 I 2 ) , D α ξ 2 = − ∂ ℋ ∂ V = ξ 2 ( μ ) , D α ξ 3 = − ∂ ℋ ∂ E = ξ 3 ( γ + μ ) − ξ 4 ( γ ) , D α ξ 4 = − ∂ ℋ ∂ I 1 = ξ 4 ( ε ρ + ε ( 1 − ρ ) + μ ) − ξ 5 ( ε ρ ) − ξ 6 ( ε ( 1 − ρ ) ) , D α ξ 5 = − ∂ ℋ ∂ R = ξ 5 ( μ ) , D α ξ 6 = − ∂ ℋ ∂ I 2 = ξ 6 ( u 2 − δ 1 − μ ) − ξ 7 ( u 2 ) , D α ξ 6 = − ∂ ℋ ∂ T = ξ 7 ( δ 2 + μ ) .

We then establish the following optimum conditions using the Pontryagin principle:

(4.5) ∂ ℋ ∂ u 1 = 2 k 3 u 1 − ( ξ 1 − ξ 2 ) S = 0 ,

which may be solved using state and adjoint variables to yield

(4.6) u 1 ¯ = ( ξ 2 − ξ 2 ) S 2 k 3 .

Similarly,

(4.7) u 2 ¯ = ( ξ 6 − ξ 7 ) I 2 2 k 4 .

For the best controls u 1 * and u 2 * , we take into account the control constraints as well as the signs of ∂ ℋ ∂ u 1 and ∂ ℋ ∂ u 2 .

Therefore, we have

(4.8) u 1 * = 0 , if ∂ ℋ ∂ u 1 < 0 , u 1 ¯ , if ∂ ℋ ∂ u 1 = 0 , 1 , if ∂ ℋ ∂ u 1 > 0 ,

and u 1 * = max { min { u 1 ¯ , 1 } , 0 } , where u 1 ¯ = ( ξ 1 − ξ 2 ) S 2 k 3 .

Similarly,

(4.9) u 2 * = 0 , if ∂ ℋ ∂ u 2 < 0 , u 2 ¯ , if ∂ ℋ ∂ u 2 = 0 , 1 , if ∂ ℋ ∂ u 2 > 0 ,

and u 2 * = max { min { u 2 ¯ , 1 } , 0 } , where u 2 ¯ = ( ξ 6 − ξ 7 ) I 2 2 k 4 .

By substituting u 1 * and u 2 * into the equation, the optimal condition may be determined for System (4.16).□

5 Numerical simulations

After simplification for three terms using the following values [11],

n 1 = 100 , n 2 = 15 , n 3 = 5 , n 4 = 6 , n 5 = 10 , n 6 = 3 , n 7 = 4 , Λ = 0.0260 , σ = 0.8 , β 1 = 0.095 , β 2 = 0.25 , γ = 0.03 , ε δ = 3.6 , ε ( 1 − ρ ) = 0.4 , δ 1 = 0.0063 , η = 0.025 , δ 2 = 0.051 , and μ = 0.01890 ,

System (3.16) can be written as:

S 0 = 100 , V 0 = 15 , E 0 = 5 , I 1 0 = 6 , R 0 = 10 , I 2 0 = 3 , T 0 = 4 , S 1 = − 213.864 t α Γ ( α + 1 ) , V 1 = 79.7165 t α Γ ( α + 1 ) , E 1 = 137.7555 t α Γ ( α + 1 ) , I 1 1 = − 23.9634 t α Γ ( α + 1 ) , R 1 = 21.411 t α Γ ( α + 1 ) , I 2 1 = 2.2494 t α Γ ( α + 1 ) , T 1 = − 0.2046 t α Γ ( α + 1 ) S 2 = 5166.58 t 2 α Γ ( 2 α + 1 ) , V 2 = − 172.598 t 2 α Γ ( 2 α + 1 ) , E 2 = − 337.269307 t 2 α Γ ( 2 α + 1 ) , I 1 2 = 100.439 t 2 α Γ ( 2 α + 1 ) , R 2 = − 86.67 t 2 α Γ ( 2 α + 1 ) , I 2 2 = − 9.698 t 2 α Γ ( 2 α + 1 ) , and T 2 = 0.07053654 t α Γ ( α + 1 )

Thus, solution after three terms is given as follows:

S ( t ) = 100 − 213.864 t α Γ ( α + 1 ) + 5166.58 t 2 α Γ ( 2 α + 1 ) , V ( t ) = 15 + 79.7165 t α Γ ( α + 1 ) − 172.598 t 2 α Γ ( 2 α + 1 ) , E ( t ) = 5 + 137.7555 t α Γ ( α + 1 ) − 337.269307 t 2 α Γ ( 2 α + 1 ) , I 1 ( t ) = 6 − 23.9634 t α Γ ( α + 1 ) + 100.439 t 2 α Γ ( 2 α + 1 ) , R ( t ) = 10 + 21.411 t α Γ ( α + 1 ) − 86.67 t 2 α Γ ( 2 α + 1 ) , I 2 ( t ) = 3 + 2.2494 t α Γ ( α + 1 ) − 9.698 t 2 α Γ ( 2 α + 1 ) , T ( t ) = 4 − 0.2046 t α Γ ( α + 1 ) + 0.07053654 t α Γ ( α + 1 ) .

In particular, for taking α = 1 , the solution of proposed Model (1.1) after three terms is given by Figures 1, 2, 3, 4, 5, 6, 7.

(5.1) S ( t ) = 100 − 213.86 t + 2583.29 t 2 , V ( t ) = 15 + 79.72 t − 86.30 t 2 , E ( t ) = 5 + 137.76 t − 188.64 t 2 , I 1 ( t ) = 6 − 23.97 t + 50.23 t 2 , R ( t ) = 10 + 21.41 t − 43.34 t 2 , I 2 ( t ) = 3 + 2.25 t − 4.85 t 2 , T ( t ) = 4 − 0.20 t + 0.04 t 2 .

$Figure 1 Evolutions of the susceptible population with time for different values of α \alpha .$

Figure 1

Evolutions of the susceptible population with time for different values of α .

$Figure 2 Density of the vaccinated population with time for different values of α \alpha .$

Figure 2

Density of the vaccinated population with time for different values of α .

$Figure 3 Density of the exposed population with time for different values of α \alpha .$

Figure 3

Density of the exposed population with time for different values of α .

$Figure 4 Dynamics of the acutely infected population with time for different values of α \alpha .$

Figure 4

Dynamics of the acutely infected population with time for different values of α .

$Figure 5 Evolution of the recovered population with time for different values of α \alpha .$

Figure 5

Evolution of the recovered population with time for different values of α .

$Figure 6 Dynamics of the chronically infected population with time for different values of α \alpha .$

Figure 6

Dynamics of the chronically infected population with time for different values of α .

$Figure 7 Dynamics of the acutely population under treatment with time for different values of α \alpha .$

Figure 7

Dynamics of the acutely population under treatment with time for different values of α .

6 Convergence analysis

The series solution obtained earlier converges rapidly, and it converges uniformly to the exact solution. To verify the convergence of the series, we use the following theorems by using the approach in [15].

Theorem 6.1

Let X be a Banach space and F : X → X be a contractive nonlinear operator such that for all x , x ′ ε X , ∣ ∣ F ( x ) − F ( x ′ ) ∣ ∣ ≤ k ∣ ∣ x − x ′ ∣ ∣ , 0 < k < 1 . Then, F has a unique point x such that F x = x , where x = ( S , V , E , I 1 , R , I 2 , T ) . The series obtained in (3.16) can be written by applying ADM as

x n = F x n − 1 , x n − 1 = ∑ i = 1 n − 1 x i , n = 1 , 2 , 3 , … ,

and assume that x 0 ∈ B r ( x ) , where B r ( x ) = { x ′ ∈ X : ∣ ∣ x ′ − x ∣ ∣ < r } , then, we have

x n ∈ B r ( x ) ;
lim x → ∞ x n = x .

Proof

For (1), using mathematical induction for n = 1 , we obtain

∣ ∣ x 0 − x ∣ ∣ = ∣ ∣ F ( x 0 ) − F ( x ) ∣ ∣ ≤ k ∣ ∣ x 0 − x ∣ ∣ .

Let the result be true for m − 1 , then

∣ ∣ x 0 − x ∣ ∣ ≤ k m − 1 ∣ ∣ x 0 − x ∣ ∣ .

We have

∣ ∣ x m − x ∣ ∣ = ∣ ∣ F ( x m − 1 ) − F ( x ) ∣ ∣ ≤ k ∣ ∣ x m − 1 − x ∣ ∣ ≤ k m ∣ ∣ x 0 − x ∣ ∣ ,

which can also be written as

∣ ∣ x n − x ∣ ∣ ≤ k n ∣ ∣ x 0 − x ∣ ∣ ≤ k n r < r ,

which implies that x n ∈ B r ( x ) .

For (2), since ∣ ∣ x n − x ∣ ∣ ≤ k n ∣ ∣ x 0 − x ∣ ∣ and as lim x → ∞ k n = 0 , hence, we have lim x → ∞ ∣ ∣ x n − x ∣ ∣ = 0 , which implies that lim x → ∞ x n = x .□

7 Conclusion

In this article, we have considered a fractional SV EI 1 RI 2 T , describing the spread of hepatitis B. LADM was used in obtaining the analytical solution of the fractional model. The optimal solution analysis was carefully carried out with the inclusion of the cost-effectiveness parameter. The comparison for some varying values of α has been obtained. A convergence analysis is also offered to show how effective the method is.

Funding information: This research received no specific grant from any funding agency, commercial, or nonprofit sectors.
Author contributions: Nnaemeka Stanley Aguegboh conceived the idea, performed the analysis, and prepared draft of the manuscript. Walter Oknongo and Boubacar Diallo carried out the simulations. Phineas Roy Kiogora and Mutua Felix verified all the results and edited the draft.
Conflict of interest: The authors have no conflicts of interest to disclose.
Ethical approval: This research did not require ethical approval.

References

[1] Megal, P., & Ruan, S. (2014). Susceptible-infectious-recovered models revisited: From the individual level to the population level. Mathematical Biosciences, 250, 26–40. 10.1016/j.mbs.2014.02.001Search in Google Scholar PubMed

[2] Ngwa, A. G., & Shu, W. S. (2000). A mathematical model for endemic malaria with variable human and mosquito populations. Mathematical and Computer Modelling, 32, 747–763. 10.1016/S0895-7177(00)00169-2Search in Google Scholar

[3] Eric Okyere, F. T. (2016). Fractional order SIR model with constant population. British Journal of Mathematics and Computer Science, 14(2), 1–12. 10.9734/BJMCS/2016/23017Search in Google Scholar

[4] Nnaemeka, A. A., & Amanso, O. R. (2021). Analysis of a model on the transmission dynamics (with prevention and control) of Hepatitis B. Journal of Fractional Calculus and Applications, 12(1), 76–89. Search in Google Scholar

[5] Medicins Sans Frontiers- Middle East. (2016). Retrieved May 9, 2018, from https://www.msf-me.org/issues/hepatitis. Search in Google Scholar

[6] Parkin, D. M., Bray, F., Ferlay, J., & Pisani, P. (2005). Global cancer statistics, 2002. CA: A Cancer Journal for Clinicians, 55(2), 74–108. 10.3322/canjclin.55.2.74Search in Google Scholar PubMed

[7] Kafui, L. A. (2014). Mathematical of Hepatitis B in the North Tongu District of the Volta region of Ghana. Msc Thesis. Kwame Nkrumah university of science and technology. Search in Google Scholar

[8] Salam, S. M., & Yousef, A. M. (2017). On a fractional-order model for HBV infection with cure of infected cells. Journal of the Egyptian Mathematical Society, 25(4), 447–451. 10.1016/j.joems.2017.06.003Search in Google Scholar

[9] Ain, Q. T., Nadeem, M., Karim, S., Akgül, A., & Jarad, F., (2022). Optimal variational iteration method for parametric boundary value problem. AIMS Mathematics, 7(9), 16649–16656. 10.3934/math.2022912Search in Google Scholar

[10] Ain, Q. T., & Chu, Y. (2023). On fractal fractional Hepatitis B epidemic model with modified vaccination effects. Fractals, 31(10), 2340006. 10.1142/S0218348X23400066Search in Google Scholar

[11] Zhao, S., Xu, Z., & Lu, Y. (2000). A mathematical model of hepatitis b virus transmission and its application for vaccination strategy in China. International Journal of Epidemiology 29(4), 744–752. 10.1093/ije/29.4.744Search in Google Scholar PubMed

[12] Nnaemeka, S. A., Kiogora, R. P., Mutua, F., & Boubacar, D. (2023). Modeling and control of Hepatitis B Virus transmission dynamics using fractional-order differential equations. Communications in Mathematical Biology and Neuroscience, 2023, 107. Search in Google Scholar

[13] Fazal, H., Kamal, S., Ghaus, R., & Muhammad, S. (2018). Numerical solution of fractional-order smoking model via Laplace Adomian decomposition method. Alexandria Engineering Journal, 57(2), 1061–1069. 10.1016/j.aej.2017.02.015Search in Google Scholar

[14] Bergounioux, M., & Bourdin, L. (2020). Pontryagin maximum principle for general Caputo fractional optimal control problems with Bolza cost and terminal constraints, ESAIM - Control, Optimisation and Calculus of Variations, 26, 35. 10.1051/cocv/2019021Search in Google Scholar

[15] Abdelrazec, A., & Pelinovsky, D. (2011). Convergence of the Adomian decomposition method for initial-value problems. Numerical Methods for Partial Differential Equations, 27(4), 749–766. 10.1002/num.20549Search in Google Scholar

Received: 2023-10-10

Revised: 2023-12-07

Accepted: 2024-06-21

Published Online: 2024-07-29

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmb-2023-0114

Keywords for this article

hepatitis B; fractional-order differential equations; Laplace Adomian decomposition; optimal control; numerical simulations; predictor-corrector method

Creative Commons

BY 4.0