Asymptotic analysis of hepatitis B epidemic model using Caputo Fabrizio fractional operator

Ting Cui; Peijiang Liu; Anwarud Din

doi:10.1515/phys-2022-0033

Open Access Published by De Gruyter Open Access May 6, 2022

Asymptotic analysis of hepatitis B epidemic model using Caputo Fabrizio fractional operator

Ting Cui , Peijiang Liu and Anwarud Din

From the journal Open Physics

https://doi.org/10.1515/phys-2022-0033

Abstract

A mathematical model representing the temporal dynamics of hepatitis B virus (HBV) is discussed in this research work. This is based on the asymptomatic carriers and symptomatic individuals keeping in view the characteristics of the disease. We also incorporate the vaccination parameter to vaccinate susceptible individuals. Moreover, we use fractional calculus to extend the model to its associated fractional-order. For this, we particularly use the fractional operator of the Caputo-Fabrizio type to fractionalize the proposed model. First, the model formulation has been derived in classic order and then extended to its associated fractional-order version for generalization. The model equilibria was calculated, and the basic reproductive number was found. Then we will discuss the existence with properties of the uniqueness of the proposed fractional version of the model that is under consideration. The positivity with boundedness is shown to investigate that the considered model is feasible biologically as well as mathematically. Finally, we use the Mittag–Leffler approach to visualize the model of fractional-order and to support the results carried out in the theocratical part. We also demonstrate the solution curves for different values of the fractional parameter to differentiate between integer-order and fractional-order on the disease transmission.

Keywords: epidemic model; hepatitis B virus; CF-fractional derivative; asymptomatic and symptomatic; numerical simulation

1 Research background

Hepatitis B is a contagious disease among other infectious diseases. This causes inflammation of the liver. Around 350 million individuals globally have chronic HBV, with 25–40 percent of these cases leading to chronic liver disease, hepatocellular carcinoma, and cirrhosis [1,2]. As a result, HBV infection has become a huge public health issue all over the world. A patient with a chronic stage has a serum burden [3]. If we assume the average mass of the liver to be 1.5 kg, a human liver should have the same number of cells. Due to the large numbers, we must use the conventional incidence rate rather than the mass action incidence rate.

Epidemiological models yield some intriguing mathematical insights, particularly in the prevention and management of infectious diseases, see, for instance, [4,5,6, 7,8,9]. Varied models of infectious diseases with various properties and characteristics have been produced and studied in the literature. Different elements (such as temperature, precipitation, absolute humidity, and others) have a major impact on the force of infection in human viral infections. Various deterministic and stochastic models have been proposed by different authors and investigated the dynamics of infectious diseases with some real applications (see, for instance, [10,11, 12,13,14, 15,16]). However, models with fractional-order are proved to be useful for describing the complicated and cross behavior because the range of derivative order is optional in this case [17,18,19, 20,21,22]. The information about the spectrum between the integer orders has been provided via fractional-order calculus [23]. Fractional calculus is more realistic from integer order while producing more accurate outputs [24,25,26, 27,28,29]. Therefore, problems with real-world situations can be studied with derivative of fractional-order, with a greater degree of freedom. Instead of classical cases nowadays, researchers are very much interested in fractional-order derivatives while applying them to real-world problems. With this, many analyses has been made, i.e., the existence theory, uniqueness theory, approximating solutions, and stabilities of the fractional equation. Physical problems in various forms have been formulated using differential, functional, and integral equations with fractional-order, such as a logistic model, model of microorganism, dingy problem, and Lotka–Volterra equation. These equations are foundation models in the field of formulating models in case of biological situations [31,32,33, 34,35]. Further various concepts are examined for FODEs using analytical, semi-analytical, and approximation analysis. Many more methods like Euler, Taylor, Adams-Bashforth, predictor-corrector, etc., are the transforms with many more other well-known techniques [36,37]. Various fractional derivatives have been applied in different shapes of biological models to show the temporal dynamics of distinct infectious diseases (see, for more detail, [38,39,40, 41,42]).

In the present article, we analyze and formulate a novel model to represent the temporal dynamics of HBV using the CF fractional operator based on the asymptomatic, symptomatic, and vaccination assumptions from the findings of the previous studies. All of the models in the literature have distinct perspectives for studying HBV dynamics, but we use asymptomatic and symptomatic carriers with vaccination of susceptible individuals to examine the dynamics of the disease more thoroughly. HBV-infected people can live for 30–180 days without showing symptoms and can infect others, which can lead to mortality in humans [43]. As a result, the goal of the study is to forecast and develop a novel model with asymptomatic, symptomatic carries, and vaccination parameter for immunization. The model is established first with integer order and then extending to the generalized form of fractional derivatives keeping in view the CF derivative operator. We calculate the basic reproductive number and examine the equilibria of the classical model. Once the reported model fractionalizes, we will then study the existence analysis with unique properties. Finally, we perform some numerical simulations and show a graphical representation of the analytical results.

We organized the reported studies as the formulation of the novel HBV model, and primary definitions of the CF derivative are presented in Section 2. In Section 3, we discuss the analysis of existence with unique properties of a positive solution. In Section 4, some numerical visualization has been provided and also compared the results to show for different fractional values. We end our work with the conclusion in Section 5.

2 HBV transmission models

We consider a model to describe the time dynamics of hepatitis B virus transmission by following the work presented in ref. [44]. For this purpose, we will introduce some notations in the form of state variables and parameters. Particularly, we assume that W ( t ) , X ( t ) , Y ( t ) , and Z ( t ) , respectively, denote the vulnerable population, the fraction of population with HBV positive, the portion of the population with positive HBV chronic, and recovered individuals. Further, the vulnerable are those individuals who at any time may catch the disease. Likewise, we divide the total population as symbolized by N ( t ) in various compartments of W , X , Y , and Z . Obviously, the entire compartments are dependent on time, which can be written as N ( t ) = Z ( t ) + Y ( t ) + X ( t ) + W ( t ) , and so the model looks like

(1) d W ( t ) d t = Λ − β W ( t ) X ( t ) N − β d 2 W ( t ) Y ( t ) N − ( d 0 + ν ) W ( t ) − δ Λ Y ( t ) , d X ( t ) d t = β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) , d Y ( t ) d t = r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) Y ( t ) , d Z ( t ) d t = d 1 X ( t ) − d 0 Z ( t ) + ν W ( t ) + κ Y ( t ) ,

where the description of the parameters is as Λ is the rate of birth and ν is assumed to be the rate of vaccination. The rate of offspring prenatally infected portion is denoted by δ , while d 2 , d 0 , and d 1 are reduced rate of disease transmission, natural, and disease death rates respectively. Further, β and d 1 are the transmission co-efficient and recovery rate, but κ and r 1 symbolize the transfer rate of acute to chronic population.

To shorten our calculation, we assume that γ 1 = ν + d 0 , γ 1 = d 1 + r 1 + d 0 , and γ 1 = κ + d 1 + d − Λ η , then using the theory of the dynamical system, we always admits a disease-free point for our proposed model as symbolized by E 0 , which looks like

(2) E 0 = Λ γ 1 , 0 , 0 , ν Λ d 0 γ 1 .

We use the disease-free state and applying the next-generation approach to calculate the reproductive number. Let us assume that this quantity is denoted by R 0 D , then

(3) R 0 D = β Λ γ 1 γ 2 + d 2 γ 2 Λ β γ 1 γ 2 γ 3 .

Further, the dynamical properties of the reported model (1) has been described as follows:

The equilibrium point E 0 is table if R 0 D < 1 , but if R 0 D > 1 , then the model is unstable.
If R 0 D > 1 , the unique endemic state E ∗ = ( W ∗ , X ∗ , Y ∗ , Z ∗ ) is stable both locally and globally, implies that the infection of HBV persists.

We will now study the fractional-order version of the afore-reported model, so first we present some preliminary knowing of CF fractional derivatives to derive acquainted with customary notations [7,23].

Definition 1

Let ϕ ∈ H 1 ( 0 , T ) , if n − 1 < ϑ < n , ϑ > 0 and n ∈ N . The derivative in sense of Caputo as well as in Caputo-Fabrizio operators are defined by

(4) D 0 , t ϑ CF { φ ( t ) } = K ( ϑ ) ( 1 − ϑ ) ∫ 0 t φ ′ ( z ) exp ( z − t ) ϑ 1 − ϑ d z

and

(5) D 0 , t ϑ C { φ ( t ) } = 1 Γ ( − ϑ + n ) ∫ 0 t ( t − z ) − 1 + n − ϑ φ n ( z ) d z ,

where CF and C symbolizes the Caputo-Fabrizio and Caputo, respectively.

Definition 2

Let 0 < ϑ < 1 , then the integral

(6) J 0 , t ϑ RL { φ ( t ) } = 1 Γ ( ϑ ) ∫ 0 t ( t − z ) ϑ − 1 φ ( z ) d z

is the Riemann-Liouville and

(7) J 0 , t ϑ CF { φ ( t ) } = 2 ( 2 − ϑ ) K ( ϑ ) × ( 1 − ϑ ) φ ( t ) + ϑ ∫ 0 t φ ( z ) d z

is the Caputo-Fabrizio-Caputo (CF) integral operators.

We now incorporate the aforementioned fractional operator into model (1) to study the associated fraction version, which looks like

(8) D 0 , t ϑ CF W ( t ) = Λ − β W ( t ) X ( t ) N − β d 2 W ( t ) Y ( t ) N − ( d 0 + ν ) W ( t ) − δ Λ Y ( t ) , D 0 , t ϑ CF X ( t ) = β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) , D 0 , t ϑ CF Y ( t ) = r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) Y ( t ) , D 0 , t ϑ CF Z ( t ) = d 1 X ( t ) − d 0 Z ( t ) + ν W ( t ) + κ Y ( t ) .

3 Existence and uniqueness

The proposed fractional-order system will be transformed to its associated integral equation system to show the existence analysis. We will apply the fixed-point theorem and show existence uniqueness. Thus taking the transformation of the aforementioned system leads to the integral system as given by

W ( t ) = J 0 , t ϑ CF Λ − β W ( t ) X ( t ) N − β d 2 W ( t ) Y ( t ) N − ( d 0 + ν ) W ( t ) − δ Λ Y ( t ) } + W ( 0 ) , X ( t ) = X ( 0 ) + J 0 , t ϑ CFC β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) } , Y ( t ) = J 0 , t ϑ CFC { r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) Y ( t ) } + Y ( 0 ) , Z ( t ) = Z ( 0 ) + J 0 , t ϑ CFC { d 1 X ( t ) − d 0 Z ( t ) + ν W ( t ) + κ Y ( t ) } .

To simplify our calculations, we assume ϒ = K ( ϑ ) ( 2 − ϑ ) for simplicity, and by applying integration in CF sense, we obtain the following equivalent system:

W ( t ) = W ( 0 ) + 2 ( 1 − ϑ ) ϒ Λ − β W ( t ) X ( t ) N − β d 2 W ( t ) Y ( t ) N − ( d 0 + ν ) W ( t ) − δ Λ Y ( t ) } + 2 ϑ ∫ 0 t Λ − β W ( x ) X ( x ) N − β d 2 W ( x ) Y ( x ) N − ( d 0 + ν ) W ( x ) − δ Λ Y ( x ) } d x 1 ( 2 − ϑ ) K ( ϑ ) ,

X ( t ) = 2 ( 1 − ϑ ) K ( ϑ ) β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) } 1 ( 2 − ϑ ) + X ( 0 ) + 1 K ( ϑ ) ∫ 0 t β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) } d x 2 ϑ ( 2 − ϑ ) , Y ( t ) = ( 1 − ϑ ) K ( ϑ ) { r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) Y ( t ) } 2 ( 2 − ϑ ) + Y ( 0 ) + 2 ϑ ( 2 − ϑ ) K ( ϑ ) ∫ 0 t { r 1 X ( x ) − ( d + κ + d 0 − δ Λ ) Y ( x ) } d x , Z ( t ) = Z ( 0 ) + 2 ( 1 − ϑ ) ϒ { d 1 X ( t ) − d 0 Z ( t ) + ν W ( t ) + κ Y ( t ) } + 1 ( 2 − ϑ ) ∫ 0 t { d 1 X ( x ) − d 0 Z(x) + ν W ( x ) + κ Y ( x ) } d x 2 ϑ K ( ϑ ) .

Let l 1 , l 2 , l 3 , and l 4 represents the kernels, then

(9) l 1 ( W ( t ) , t ) = Λ − β W ( t ) X ( t ) N − β d 2 W ( t ) Y ( t ) N − ( d 0 + ν ) W ( t ) − δ Λ Y ( t ) } , l 2 ( X ( t ) , t ) = β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) } , l 3 ( Y ( t ) , t ) = { r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) Y ( t ) } , l 4 ( Z ( t ) , t ) = { d 1 X ( t ) − d 0 Z ( t ) + ν W ( t ) + κ Y ( t ) } .

Theorem 1

Lipschitz conditions hold for the aforementioned kernels l i , i − 1 , 2 , 3 , 4 .

Proof

We assume W and W 1 , X , and X 1 , Y and Y 1 , Z , and Z 1 are the functions of l i ; thus, we have

l 1 ( W ( t ) , t ) = l 1 ( W 1 ( t ) , t ) + Λ − β ( W ( t ) − W 1 ( t ) ) X ( t ) N − β d 2 ( W ( t ) − W 1 ( t ) ) Y ( t ) N − ( d 0 + ν ) ( W ( t ) − W 1 ( t ) ) − δ Λ Y ( t ) } , l 1 ( X ( t ) , t ) = l 1 ( X 1 ( t ) , t ) + β W ( t ) ( X ( t ) − X 1 ( t ) ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) ( X ( t ) − X 1 ( t ) ) , l 1 ( Y ( t ) , t ) = l 1 ( Y 1 ( t ) , t ) + { r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) ( Y ( t ) − Y 1 ( t ) ) } , l 1 ( Z ( t ) , t ) = l 1 ( Z 1 ( t ) , t ) + { d 1 X ( t ) − d 0 ( Z ( t ) − Z 1 ( t ) ) + ν W ( t ) + κ Y ( t ) } .

We apply the Cauchy’s inequality implies that

‖ l 1 ( W ( t ) , t ) − l 1 ( W 1 ( t ) , t ) ‖ ≤ Λ − β ( W ( t ) − W 1 ( t ) ) X ( t ) N − β d 2 ( W ( t ) − W 1 ( t ) ) Y ( t ) N − ( d 0 + ν ) ( W ( t ) − W 1 ( t ) ) − δ Λ Y ( t ) } ∥ , ‖ l 1 ( X ( t ) , t ) − l 1 ( X 1 ( t ) , t ) ‖ ≤ β W ( t ) ( X ( t ) − X 1 ( t ) ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) ( X ( t ) − X 1 ( t ) ) } ∥ , ‖ l 1 ( Y ( t ) , t ) − l 1 ( Y 1 ( t ) , t ) ‖ ≤ ‖ { r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) ( Y ( t ) − Y 1 ( t ) ) } ‖ , ‖ l 1 ( Z ( t ) , t ) − l 1 ( Z 1 ( t ) , t ) ‖ ≤ ‖ { d 1 X ( t ) − d 0 ( Z ( t ) − Z 1 ( t ) ) + ν W ( t ) + κ Y ( t ) } ‖ .

We then obtain the relations given by

(10) W ( t ) = 2 ( 1 − ϑ ) l 1 ( W n − 1 ( t ) , t ) ϒ + 2 ϑ ϒ ∫ 0 t l 1 ( W n − 1 ( x ) , x ) d x , X ( t ) = 2 ( 1 − ϑ ) l 2 ( X n − 1 ( t ) , t ) ϒ + 2 ϑ ϒ ∫ 0 t l 2 ( X n − 1 ( x ) , x ) d x , Y ( t ) = 2 ( 1 − ϑ ) l 3 ( Y n − 1 ( t ) , t ) ( 2 − ϑ ) K ( ϑ ) + 2 ϑ ( 2 − ϑ ) K ( ϑ ) ∫ 0 t l 3 ( Y n − 1 ( x ) , x ) d x , Z ( t ) = 2 ( 1 − ϑ ) l 4 ( Z n − 1 ( t ) , t ) ( 2 − ϑ ) K ( ϑ ) + 2 ϑ K ( ϑ ) ( 2 − ϑ ) ∫ 0 t l 4 ( Z n − 1 ( x ) , x ) d x .

Taking the difference between successive terms and applying norm as well as majorizing, we then lead to the assertions

(11) ‖ A n ( t ) ‖ = ‖ W n ( t ) − W 1 , n − 1 ( t ) ‖ ≤ 2 ( 1 − ϑ ) ϒ ‖ l 1 ( W n − 1 ( t ) , t ) − l 1 ( W 1 , n − 2 ( t ) , t ) ‖ + 2 ϑ K ( ϑ ) ( 2 − ϑ ) ∫ 0 t [ l 1 ( W n − 1 ( x ) , x ) − l 1 ( W 1 , n − 2 ( x ) , x ) ] d x ∥ , ‖ B n ( t ) ‖ = ‖ X n ( t ) − X 1 , n − 1 ( t ) ‖ ≤ ‖ l 2 ( X n − 1 ( t ) , t ) − l 2 ( X 1 , n − 2 ( t ) , t ) ‖ 2 ( 1 − ϑ ) K ( ϑ ) ( 2 − ϑ ) + 2 ϑ K ( ϑ ) ( 2 − ϑ ) ∫ 0 t [ l 2 ( X n − 1 ( x ) , x ) − l 2 ( X 1 , n − 2 ( x ) , x ) ] d x ∥ , ‖ C n ( t ) ‖ = ‖ Y n ( t ) − Y 1 , n − 1 ( t ) ‖ ≤ 2 ( 1 − ϑ ) K ( ϑ ) ( 3 − ϑ ) ‖ l 3 ( Y n − 1 ( t ) , t ) − l 1 ( Y 1 , n − 2 ( t ) , t ) ‖ + 2 ϑ K ( ϑ ) ( 2 − ϑ ) ∫ 0 t [ l 3 ( Y n − 1 ( x ) , x ) − l 3 ( Y 1 , n − 2 ( x ) , x ) ] d x ∥ , ‖ D n ( t ) ‖ = ‖ Z n ( t ) − Z 1 , n − 1 ( t ) ‖ ≤ ‖ l 4 ( Z n − 1 ( t ) , t ) − l 4 ( Z 1 , n − 2 ( t ) , t ) ‖ 2 ( 1 − ϑ ) K ( ϑ ) ( 2 − ϑ ) + 2 ϑ ϒ ∫ 0 t [ l 4 ( Z n − 1 ( x ) , x ) − l 4 ( Z 1 , n − 2 ( x ) , x ) ] d x ∥ .

In the aforementioned system

(12) ∑ i = 0 ∞ A i ( t ) = W n ( t ) , ∑ i = 0 ∞ B i ( t ) = X n ( t ) , ∑ i = 0 ∞ C i ( t ) = Y n ( t ) , ∑ i = 0 ∞ D i ( t ) = Z n ( t ) .

Since, the kernels l i satisfy the conditions of Lipschitz, so

(13) ‖ A n ( t ) ‖ = ‖ W n ( t ) − W 1 , n − 1 ( t ) ‖ ≤ 2 ( 1 − ϑ ) ϒ τ 1 ‖ W n − 1 ( t ) − W 1 , n − 2 ( t ) ‖ + 2 ϑ ( 2 − ϑ ) K ( ϑ ) τ 2 ∫ 0 t ∥ W n − 1 ( x ) − W 1 , n − 2 ( x ) ∥ d x , ‖ B n ( t ) ‖ = ‖ X n ( t ) − X 1 , n − 1 ( t ) ‖ ≤ τ 3 ‖ X n − 1 ( t ) − X 1 , n − 2 ( t ) ‖ 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ τ 4 ∫ 0 t ‖ X n − 1 ( x ) − X 1 , n − 2 ( x ) ‖ d x , ‖ C n ( t ) ‖ = ‖ Y n ( t ) − Y 1 , n − 1 ( t ) ‖ ≤ τ 5 ‖ Y n − 1 ( t ) − Y 1 , n − 2 ( t ) ‖ 2 ( 1 − ϑ ) ϒ + 2 ϑ ( 2 − ϑ ) K ( ϑ ) τ 6 ∫ 0 t ‖ Y n − 1 ( x ) − Y 1 , n − 2 ( x ) ‖ d x , ‖ D n ( t ) ‖ = ‖ Z n ( t ) − Z 1 , n − 1 ( t ) ‖ ≤ τ 7 ‖ Z n − 1 ( t ) − Z 1 , n − 2 ( t ) ‖ ( 1 − ϑ ) ( 2 ) ( 2 − ϑ ) K ( ϑ ) + τ 8 ∫ 0 t ‖ Z n − 1 ( x ) − Z 1 , n − 2 ( x ) ‖ d x 2 ϑ ( 2 − ϑ ) K ( ϑ ) .□

Theorem 2

Model (8) poses a solution.

Proof

Using Eq. (12) with the recursive scheme leads to the assertion as follows:

‖ A n ( t ) ‖ ≤ ‖ W ( 0 ) ‖ + 2 ( 1 − ϑ ) τ 1 ϒ n + 2 t τ 2 ϑ ϒ n , ‖ B n ( t ) ‖ ≤ ‖ X ( 0 ) ‖ + 2 τ 3 ( 1 − ϑ ) ϒ n + 2 τ 4 t ϑ ( 2 − ϑ ) K ( ϑ ) n ,

(14) ‖ C n ( t ) ‖ ≤ ‖ Y ( 0 ) ‖ + 2 τ 5 ( 1 − ϑ ) ϒ n + 2 τ 6 ϑ t K ( ϑ ) ( 2 − ϑ ) n , ‖ D n ( t ) ‖ ≤ ‖ Z ( 0 ) ‖ + 2 τ 7 ( 1 − ϑ ) ( 2 − ϑ ) K ( ϑ ) n + 2 τ 8 t ϑ K ( ϑ ) ( 2 − ϑ ) n .

Furthermore, let us assume that ϒ 1 , i ( t ) , i = 1 , … , 5 represent the remainder terms, and then to prove that Eq. (14) satisfies the model (8), we make the substitution as follows:

(15) W ( t ) = W n ( t ) + ϒ 1 , n ( t ) , X ( t ) = X n ( t ) + ϒ 2 , n ( t ) , Y ( t ) + ϒ 3 , n ( t ) = Y n ( t ) , Z ( t ) = + ϒ 4 , n ( t ) + Z n ( t ) .

Thus,

(16) W ( t ) − W n − 1 ( t ) = 2 ( 1 − ϑ ) l 1 ( W ( t ) − ϒ 1 , n ( t ) ) K ( ϑ ) ( 2 − ϑ ) + 2 ϑ K ( ϑ ) ( 2 − ϑ ) ∫ 0 t l 1 ( W ( x ) − ϒ 1 , n ( x ) ) d x , X ( t ) − X n − 1 ( t ) = 2 l 2 ( X ( t ) − ϒ 2 , n ( t ) ) ( 1 − ϑ ) K ( ϑ ) ( 2 − ϑ ) + 2 ϑ ( 2 − ϑ ) K ( ϑ ) ∫ 0 t l 2 ( X ( x ) − ϒ 2 , n ( x ) ) d x , Y ( t ) − Y n − 1 ( t ) = 2 l 3 ( Y ( t ) − ϒ 3 , n ( t ) ) ( 1 − ϑ ) K ( ϑ ) ( 2 − ϑ ) + 2 ϑ ( 2 − ϑ ) K ( ϑ ) ∫ 0 t l 3 ( Y ( x ) − ϒ 3 , n ( x ) ) d x , Z ( t ) − Z n − 1 ( t ) = 2 ( 1 − ϑ ) l 4 ( Z ( t ) − ϒ 4 , n ( t ) ) ( 2 − ϑ ) K ( ϑ ) + 2 ϑ ( 2 − ϑ ) K ( ϑ ) ∫ 0 t l 4 ( Z ( x ) − ϒ 4 , n ( x ) ) d x .

The application of norm with Lipschitz conditions may lead to

(17) W ( t ) − 2 l 1 ( W ( t ) , t ) ( 1 − ϑ ) K ( ϑ ) 1 ( 2 − ϑ ) − W ( 0 ) − 2 ϑ ϒ ∫ 0 t l 1 ( W ( x ) , x ) d x ≤ ‖ ϒ 1 , n ( t ) ‖ 1 + 2 ( 1 − ϑ ) τ 1 ϒ + 2 ϑ τ 2 t ϒ , X ( t ) − 2 l 2 ( X ( t ) , t ) ( 1 − ϑ ) ϒ − 2 ϑ ( 2 − ϑ ) K ( ϑ ) ∫ 0 t l 2 ( X ( x ) , x ) d x − X ( 0 ) ≤ ‖ ϒ 2 , n ( t ) ‖ 1 + 2 ( 1 − ϑ ) τ 3 ϒ + 2 ϑ τ 4 t ϒ , Y ( t ) − 2 l 3 ( Y ( t ) , t ) ( 1 − ϑ ) ϒ − Y ( 0 ) − 2 ϑ ( 2 − ϑ ) K ( ϑ ) ∫ 0 t l 3 ( Y ( x ) , x ) d x ≤ ‖ ϒ 3 , n ( t ) ‖ 1 + 2 τ 5 ( 1 − ϑ ) ϒ + 2 ϑ τ 6 t ϒ , Z ( t ) − 2 l 4 ( Z ( t ) , t ) ( 1 − ϑ ) ϒ − 2 ϑ ( 2 − ϑ ) K ( ϑ ) ∫ 0 t l 4 ( Z ( x ) , x ) d x − Z ( 0 ) ≤ ‖ ϒ 4 , n ( t ) ‖ 1 + 2 ( 1 − ϑ ) τ 7 ϒ + 2 ϑ τ 8 t ϒ .

Taking lim as t approaches ∞ , we reach to

(18) W ( t ) = 2 l 1 ( W ( t ) , t ) ( 1 − ϑ ) ϒ + ∫ 0 t l 1 ( W ( x ) , x ) d x 2 ϑ ϒ + W ( 0 ) , X ( t ) = 2 l 2 ( X ( t ) , t ) ( 1 − ϑ ) ϒ + ∫ 0 t l 2 ( X ( x ) , x ) d x 2 ϑ ϒ + X ( 0 ) , Y ( t ) = 2 l 3 ( Y ( t ) , t ) ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t l 3 ( Y ( x ) , x ) d x + Y ( 0 ) , Z ( t ) = 2 ( 1 − ϑ ) l 4 ( Z ( t ) , t ) ϒ + 2 ϑ ϒ ∫ 0 t l 4 ( Z ( x ) , x ) d x + Z ( 0 ) .

This completes the conclusion.□

Theorem 3

The proposed model as reported by system (8)exists a unique solution.

Proof

Let ( W + ( t ) , X + ( t ) , Y + ( t ) , Z + ( t ) ) is also the solution, then

(19) W ( t ) = W + ( t ) + { l 1 ( W ( t ) , t ) − l 1 ( W + ( t ) , t ) } 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t { l 1 ( W ( y ) , y ) − l 1 ( W + ( y ) , y ) } d x , X ( t ) = X + ( t ) + { l 2 ( X ( t ) , t ) − l 2 ( X + ( t ) , t ) } 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t { l 2 ( X ( y ) , y ) − l 2 ( X + ( x ) , x ) } d x , Y ( t ) = Y + ( t ) + { l 3 ( Y ( t ) , t ) − l 3 ( Y + ( t ) , t ) } 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t { l 3 ( Y ( x ) , x ) − l 3 ( Y + ( x ) , x ) } d x , Z ( t ) = Z + ( t ) + { l 4 ( Z ( t ) , t ) − l 4 ( Z + ( t ) , t ) } 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t { l 4 ( Z ( x ) , x ) − l 4 ( Z + ( x ) , x ) } d x .

The application of majorizing gives

(20) ‖ W ( t ) − W + ( t ) ‖ = ‖ l 1 ( W ( t ) , t ) − l 1 ( W + ( t ) , t ) ‖ 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t ‖ l 1 ( W ( x ) , x ) − l 1 ( W + ( x ) , x ) ‖ d x , ‖ X ( t ) − X + ( t ) ‖ = ‖ l 2 ( X ( t ) , t ) − l 2 ( X + ( t ) , t ) ‖ 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t ‖ l 2 ( X ( x ) , x ) − l 2 ( X + ( x ) , x ) ‖ d x , ‖ Y ( t ) − Y + ( t ) ‖ = ‖ l 3 ( Y ( t ) , t ) − l 3 ( Y + ( t ) , t ) ‖ 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t ‖ l 3 ( Y ( x ) , x ) − l 3 ( Y + ( x ) , x ) ‖ d x , ‖ Z ( t ) − Z + ( t ) ‖ = ‖ l 4 ( Z ( t ) , t ) − l 4 ( Z + ( t ) , t ) ‖ 2 ( 1 − ϑ ) ϒ + 2 ϑ ϒ ∫ 0 t ‖ l 4 ( Z ( x ) , x ) − l 4 ( Z + ( x ) , x ) ‖ d x .

By using Theorems 1 and 2, we may lead to the assertion as described by

(21) ‖ W ( t ) − W + ( t ) ‖ ≤ 2 τ 1 ψ 1 ( 1 − ϑ ) ϒ + 2 ϑ τ 2 t ϕ 2 ϒ n , ‖ X ( t ) − X + ( t ) ‖ ≤ 2 τ 3 ( 1 − ϑ ) ψ 3 ϒ + 2 τ 4 ϑ ϕ 4 t ϒ n , ‖ Y ( t ) − Y + ( t ) ‖ ≤ 2 τ 5 ψ 5 ( 1 − ϑ ) ϒ + 2 ϑ τ 6 ϕ 6 t ϒ n , ‖ Z ( t ) − Z + ( t ) ‖ ≤ 2 τ 7 ψ 7 ( 1 − ϑ ) ϒ + 2 ϑ τ 8 ϕ 8 t ϒ n .

This implies that the aforementioned formulas hold for every n , so we obtain

(22) W ( t ) = W + ( t ) , X ( t ) = X + ( t ) , Y ( t ) = Y + ( t ) , Z ( t ) = Z + ( t ) . □

Further, we will discuss that the model under consideration is feasible by Lemma 1.

Lemma 1

Since ( W ( t ) , X ( t ) , Y ( t ) , Z ( t ) ) is the solution of system (8). Let us assume that the initial population sizes are nonnegative, then ( W ( t ) , X ( t ) , Y ( t ) , Z ( t ) ) is also non-negative for every t ≥ 0 .

Proof

Since ϑ is the fractional parameter, then

(23) D 0 , t ϑ G ( W ( t ) ) = Λ − β W ( t ) X ( t ) N − β d 2 W ( t ) Y ( t ) N − ( d 0 + ν ) W ( t ) − δ Λ Y ( t ) , D 0 , t ϑ G ( X ( t ) ) = β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) , D 0 , t ϑ G ( Y ( t ) ) = r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) Y ( t ) , D 0 , t ϑ G ( Z ( t ) ) = d 1 X ( t ) − d 0 Z ( t ) + ν W ( t ) + κ Y ( t ) ,

where G is a general fractional operator. Thus, we have

(24) D 0 , t ϑ G ( W ( t ) ) ∣ κ ( W ( t ) ) = Λ > 0 , D 0 , t ϑ G ( X ( t ) ) ∣ κ ( X ( t ) ) = β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N ≥ 0 D 0 , t ϑ G ( Y ( t ) ) ∣ κ ( Y ( t ) ) = r 1 X ( t ) ≥ 0 , D 0 , t ϑ G ( Z ( t ) ) ∣ κ ( Z ( t ) ) = d 1 X ( t ) + ν W ( t ) + κ Y ( t ) ≥ 0 ,

where κ ( ξ ) = { ξ = 0 and W ( t ) , X ( t ) , Y ( t ) , Z ( t ) is in C ( R + × R + ) } and ξ ∈ { W ( t ) , X ( t ) , Y ( t ) , Z ( t ) } . We then reach to our conclusion by using the same methodology as adopted in [30], that the solutions of the considered model are nonnegative for all t ≥ 0 .□

4 Parameters of estimation and numerical simulation

In this section, we will present some numerical simulations to show the significance of our model. We will first estimate the model parameters and then use the estimated value to present the long-run numerical simulation.

4.1 Estimation of parameters

We use real data to estimate our model parameters that are important to show the validation of the reported model. We used a number of infected cases as reported in district Swat Khyber Pakhtunkhwa Pakistan from January 1, 2020, to December 28, 2020. To estimate the value of epidemic parameters, the method of parameter estimation technique under lsqcurvefit routine via MATLAB software has been used, while to measure the goodness of fit, the associated average relative error of formula 1 12 ∑ k = 1 12 x k real − x k approximate x k real ≈ 1.5685 × 10⁻¹ is used. This verifies the suitability of small relative error’s value (1.5685 × 10⁻¹). Once to perform the estimation analysis as presented in Figure 1, we estimate the epidemic parameters presented in Table 1.

Figure 1

The graph visualizes the fitting of model parameters to real data for the reported HBV model.

Table 1

The estimated value of epidemic parameters

Parameter	Description	Value	Source
Λ	Birth rate	260,479	Assumed
β	Transmission coefficient	14.24	Fitted
d 2	Reduced transmission coefficient	7.02529	Fitted
ν	Vaccination rate	5.8480	Fitted
r 1	Portion of acute leads to chronic	5.5650	Fitted
d 1	Recovery rate of acute	0.97	Fitted
d 0	Natural death	0.0005	Fitted
d 1	Death from disease	0.001253133	Fitted
δ	Offspring prenatally infected portion	0.0004	Fitted
κ	Rate of recovery of chronic	0.00001	Fitted

4.2 Numerical simulation

We will present some graphic descriptions of the analytical results to present the justification as well as the authenticity of the analytical work. We present the simulation of the reported fractional HBV epidemic model (8). For this, we assume that the units of time is 50, while the value of parameters is taken from Table 2. Moreover, we also use various values of fractional-order parameters to show the impact of fractional-order on disease transmission. We further assume that time step integration is h = 1 0 − 3 and [ 0 , t ] is the interval of time. Moreover, let n = T h , n ∈ N , and u = 0 , 1 , 2 , … , n , then the discretization of the proposed model looks like:

(25) W u + 1 CF = W ( 0 ) + ( 1 − ϑ ) Λ − β W ( t ) X ( t ) N − β d 2 W ( t ) Y ( t ) N − ( d 0 + ν ) W ( t ) − δ Λ Y ( t ) } + ϑ h ∑ k = 0 u Λ − β W ( t ) X ( t ) N − β d 2 W ( t ) Y ( t ) N − ( d 0 + ν ) W ( t ) − δ Λ Y ( t ) } , X u + 1 CF = X ( 0 ) + ( 1 − ϑ ) β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) } + ϑ h ∑ k = 0 u β W ( t ) X ( t ) N + β d 2 W ( t ) Y ( t ) N − ( d 1 + r 1 + d 0 ) X ( t ) } , Y u + 1 CF = Y ( 0 ) + ( 1 − ϑ ) { r 1 X ( t ) − ( d + κ + d 0 − δ Λ ) Y ( t ) }