Qualitative and Quantitative Analysis of Fractional Dynamics of Infectious Diseases with Control Measures

Abstract

Infectious diseases can have a significant economic impact, both in terms of healthcare costs and lost productivity. This can be particularly significant in developing countries, where infectious diseases are more prevalent, and healthcare systems may be less equipped to handle them. It is recognized that the hepatitis B virus (HBV) infection remains a critical global public health issue. In this study, we develop a comprehensive model for HBV infection that includes vaccination and hospitalization through a fractional framework. It has been shown that the solutions of the recommended system of HBV infection are positive and bounded. We examine the steady states of the model and determine the basic reproduction number; denoted by ${\mathcal{R}}_0$. The qualitative and quantitative behavior of the model is demonstrated using mathematical skills and numerical techniques. It has been proved that the infection-free steady state of the system is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable otherwise. Furthermore, the Ulam–Hyers stability (UHS) of the recommended fractional models is investigated and the significant conditions are provided. We present an iterative technique to visualize the dynamical behavior of the system. We perform different simulations to illustrate the effect of different input factors on the solution pathways of the system of HBV infection to conceptualize the role of parameters in the control and prevention of the infection.

1. Introduction

Infectious diseases are infections brought on by organisms, often microscopic in size, such as bacteria, viruses, fungi, or parasites that are transferred from host to host either indirectly or directly. Hepatitis B is a dangerous liver disease among infectious diseases which results in cancer, scarring of the organ, and liver failure. This disease affects the entire world and is one of the dominant causes of mortality. A vaccination for this illness has been available since 1982 and is 95% effective [1]. About 80% of primary liver cancer cases are caused by this illness. About a quarter of all paediatric illnesses progress to chronic infections. In China, approximately 93 million individuals are affected by the infection of HBV [2]. This infection transmits from one person to another when they use infectious blood, infectious body fluids of infected individuals, and through contaminated tools including syringes and in operations, etc. [1,3]. A vertical sort of transmission occurs when a newborn infant becomes infected from their infected mother. Sometimes a person infected by the HBV has no symptoms, which leads them to a severe unstable status or even death. However, after two to five months, the HBV infection reveals its symptoms, moreover, liver cancer and cirrhosis are the main causes of HBV infection.

Mathematical models are very important tools for effective investigation of the transmission and management of any infectious disease. Many mathematicians and biologists have created different epidemic models to represent the transmission dynamics of HBV. Medley et al. [4] explored catastrophic behavior in the dynamics of HBV for the first time. The authors [5] proposed an important technique for controlling the infection. Pang et al. investigated a model with immunized class and some other controls to eradicate the HBV disease [3]. The occurrence of the HBV infection was modelled and investigated by researchers in China [6]. In [7], the authors studied and explored the CTL immunological responses using a mathematical model. In [8,9], the researchers constructed and investigated diffusion models for HBV infection with a time delay. Zhao et al. established a model for calculating the vaccination’s long-term effectiveness [10], while Khan et al. considered the consequences of migrants on the dynamics of the infection [11]. In [12], the authors formulated a model to investigate the dynamics of the infection in Xinjiang. In addition to this, several studies have been presented in the literature to provide effective control interventions for this infectious disease [13,14]. The authors in [15] introduced a model of HBV with the effect of heroin and provided some control measures through control theory. This viral infection is the leading cause of liver cancer and is responsible for over 780,000 deaths each year. In addition to this, HBV infection also imposes a significant economic burden, with the costs associated with treating HBV-related liver diseases and lost productivity estimated to be in the billions of dollars each year. Therefore, we aim to develop an epidemic model for the intricate transmission phenomena of the HBV infection that includes vaccination and medication administered during hospitalization, to explore how these strategies affect the transmission and prevalence of HBV, and how they can be optimized to achieve the greatest reduction in the disease burden.

Fractional calculus involves the use of non-integer derivatives and integrals. The application of fractional derivatives is becoming increasingly important in many fields, including physics, engineering, finance, and biology [16,17]. In physics and engineering, fractional derivatives are used to model and analyze complex systems that exhibit non-linear behavior, such as viscoelastic materials, signal processing, and control systems [18]. Fractional calculus provides a more accurate description of these systems compared to traditional integer-order calculus. It has been presented in the literature [19,20] that non-integer derivatives can accurately depict biological processes because of their inherent characteristics and non-local behavior. In this study, our objective is to construct an epidemic model that accurately represents the intricate transmission dynamics of HBV infection with the effect of vaccination and medication administered during hospitalization through a fractional framework. Furthermore, we aim to demonstrate how the dynamical behavior of the system is affected by various input factors and memory index. The purpose of this work is to highlight the most effective scenario for the prevention of HBV infection and to illustrate whether the fractional parameter may be utilized as a control parameter or not.

The remaining part of this paper is structured as follows: In Section 2, we provide the rudimentary theory and concepts of the Caputo fractional derivative. We construct an epidemic model of HBV infection with vaccination and medication through hospitalization in Section 3. The suggested dynamics of HBV infection is represented in a fractional framework for accurate outcomes. We investigate the epidemic model for its basic properties and determine the reproduction parameter ${\mathcal{R}}_0$. The local and global stability results are established in Section 4 of this work. We investigate the solution of the system and demonstrated the Ulam–Hyers stability in Section 5. The tracking path behavior of the model is illustrated in Section 6 through an iterative method. We also conceptualize the impact of different input factors on the dynamics of HBV infection. Finally, the ending remarks of the overall work are given in Section 7.

2. Theory of Fractional Calculus

Here, we will introduce the rudimentary ideas and terminology of fractional theory, which will be utilized for the analysis of the system. Additionally, it has several beneficial uses across a variety of scientific areas. The essential concepts of the fractional derivative (FD) of Liouville and Caputo, introduced in [21,22,23], are as follows:

Definition 1. 

Let us take a function $g\left(r\right)$ such that $g\left(r\right)\in L^1([c,d],R)$, then, the fractional integral of Liouville–Caputo (LC) is denoted by I and is given by

$${}^{\xi }{I}_{c+}^c\left(g\left(r\right)\right)=\frac{1}{\mathsf{\Gamma }\left(\xi \right)}{\int }_0^r{(r-s)}^{\xi -1}g\left(s\right)ds,$$

(1)

where $0<\xi \le 1$ is the order of the fractional derivative and r is assumed for time.

Definition 2. 

Consider $g\left(r\right)$ in a manner that $g\left(r\right)\in C^n[c,d]$; then, the fractional derivative of $g\left(r\right)$ in the LC form is denoted by ${}^{LC}{D}_{0^{+}}^{\xi }g\left(r\right)$ and is

$${}^{LC}{D}_{0^{+}}^{\xi }g\left(r\right)=\frac{1}{\mathsf{\Gamma }(n-\xi )}{\int }_0^r{(r-s)}^{n-\xi -1}{g}^n\left(s\right)ds,$$

(2)

subject to the constraint that $n-1<\xi \le n$.

Lemma 1. 

Let us take the setup described below

$$\left\{\begin{array}{cc}\hfill {}^{LC}{D}_{0^{+}}^{\xi }g\left(r\right)=& w\left(r\right),r\in [0,\tau ],n-1<\xi <n,\\ \hfill g\left(0\right)=w_0\end{array}\right.$$

(3)

the above system has a solution given by

$$g\left(r\right)=\sum _{j=0}^{n-1}{d}_j{r}^j,where{d}_j\in R,j=0,1,\dots ,n-1.$$

(4)

Theorem 1 

([23]). Select X to be a Banach space and $T:X\to X$ to be compact and continuous. Then, T has a fixed point if

$$\mathrm{G}=\{L\in X:L=\eta TL,\eta \in (0,1\left)\right\},$$

is bounded.

3. Evaluation of the Model

To construct the dynamics of HBV, we indicate the total human population by $\mathrm{N}\left(t\right)$, which is further distributed in seven different classes: susceptible, vaccinated, exposed, acutely infected, chronic HBV case, hospitalized, and recovered individuals, expressed by $\mathrm{S}\left(t\right)$, $\mathrm{V}\left(t\right)$, $\mathrm{E}\left(t\right)$, $\mathrm{A}\left(t\right)$, $\mathrm{C}\left(t\right)$, $\mathrm{H}\left(t\right)$, and $\mathrm{R}\left(t\right)$, respectively.

The birth rate of newborn babies ($\mathsf{\Pi }$) generates the susceptible population of humans, and successfully immunized newborns ($\mathsf{\Pi }v$, where $0<v<1)$ transfer to $\mathrm{V}\left(t\right)$ while the un-immunized newborn babies are expressed, with a factor q, as carriers of the HBV infection. In addition to this, the fraction of newborns $\left(\mathsf{\Pi }\right(q-qv\left)\mathrm{C}\right)$ joins the chronic class, members of which are not fully immunized. This population is decreased due to the effective contact of the susceptible individuals with the population that are infected with HBV at a rate $\frac{(f_{1}A+f_{2}C+f_3\mathrm{H})}{\mathrm{N}}$, where $f_1,f_2$, and $f_3$ are the different rates of transmission due to $\mathrm{A}\left(t\right)$, $\mathrm{C}\left(t\right)$, and $\mathrm{H}\left(t\right)$; moreover, the said portion moves to $\left(\mathrm{E}\right(t\left)\right)$ to generate the exposed class. The vaccinated people leave this class and join the vaccinated class at the rate $g_1$. In addition, those vaccinated individuals who lost their immunity join the $\mathrm{S}\left(t\right)$ class at the rate $g_2$. The rate of death which occurs naturally, is considered to be d in all epidemiological classes. The number of individuals in $\mathrm{E}\left(t\right)$ decreases by the transfer of exposed people to the acutely infected class at the rate $ϵ$ and also the rate d. Similarly, the acutely infected $\mathrm{A}\left(t\right)$ population is decreased at the rate d and treatment rate ${\delta }_1$. Furthermore, this class is also decreased by $\alpha h_1$ to $\mathrm{C}\left(t\right)$ and a portion $(1-\alpha )h_1$, where $0<\alpha <1$, to the recovered class. The $\mathrm{C}\left(t\right)$ class is also joined by the hospitalized at a rate $\sigma$, which expresses medication failure in the hospitalized class. It is also reduced by the rate d and the rate of recovery $h_2$, as well as by medication at ${\delta }_2$. The hospitalized class $\mathrm{H}\left(t\right)$ is reduced by the rate d and recovery of individuals at the rate $\varphi$, as well as the incomplete rate of medication $\sigma$. Finally, the remaining two classes, of vaccination $\mathrm{V}\left(t\right)$ and recovered $\mathrm{R}\left(t\right)$, are decreased by the death rate d. By utilizing all of this above-mentioned information about the dynamics of a life threatening infection, we have the following system of differential equations for the vaccinated and hospitalized classes. Consequently, all of the equations may be expressed as the nonlinear ODE system:

$$\left\{\begin{array}{ccc}\hfill \frac{d\mathrm{S}}{dt}& =& \mathsf{\Pi }(1-v)(1-q\mathrm{C})-(f_1\mathrm{A}+f_2\mathrm{C}+f_3\mathrm{H})\frac{\mathrm{S}}{\mathrm{N}}-(g_1+d)\mathrm{S}+g_2\mathrm{V},\hfill \\ \hfill \frac{d\mathrm{V}}{dt}& =& \mathsf{\Pi }v+g_1\mathrm{S}-(d+g_2)\mathrm{V},\hfill \\ \hfill \frac{d\mathrm{E}}{dt}& =& (f_1\mathrm{A}+f_2\mathrm{C}+f_3\mathrm{H})\frac{\mathrm{S}}{\mathrm{N}}-(d+ϵ)\mathrm{E},\hfill \\ \hfill \frac{d\mathrm{A}}{dt}& =& ϵ\mathrm{E}-(d+h_1+{\delta }_1)\mathrm{A},\hfill \\ \hfill \frac{d\mathrm{C}}{dt}& =& h_1\alpha \mathrm{A}+\mathsf{\Pi }(1-v)q\mathrm{C}-(d+h_2+{\delta }_2)\mathrm{C}+\sigma \mathrm{H},\hfill \\ \hfill \frac{d\mathrm{H}}{dt}& =& {\delta }_1\mathrm{A}+{\delta }_2\mathrm{C}-(d+\varphi +\sigma )\mathrm{H},\hfill \\ \hfill \frac{d\mathrm{R}}{dt}& =& h_1(1-\alpha )\mathrm{A}+h_2\mathrm{C}+\varphi \mathrm{H}-d\mathrm{R},\hfill \end{array}\right.$$

(5)

with

$$\mathrm{S}\left(0\right)={\mathrm{S}}_0\ge 0,\mathrm{V}\left(0\right)={\mathrm{V}}_0\ge 0,\mathrm{E}\left(0\right)={\mathrm{E}}_0\ge 0,\mathrm{A}\left(0\right)={\mathrm{A}}_0\ge 0,$$

and

$$\mathrm{C}\left(0\right)={\mathrm{C}}_0\ge 0,\mathrm{H}\left(0\right)={\mathrm{H}}_0\ge 0,\mathrm{R}\left(0\right)={\mathrm{R}}_0\ge 0.$$

A fractional derivative is a generalization of an ordinary derivative that can capture the memory effects and long-range correlations present in complex biological phenomena and epidemics. One of the most important aspects of epidemic modeling is to accurately capture the dynamics of the spread of the disease over time. Fractional derivatives allow us to incorporate memory and non-local effects that are not accounted for in classical models. For example, the fractional derivative can model the non-exponential decay of infectiousness of a disease over time, which is a common characteristic of many epidemics.

Moreover, fractional derivatives have the capability to account for heterogeneity within populations, encompassing variations in susceptibility and contact rates among individuals. This becomes especially crucial when modeling the spread of diseases in populations characterized by diverse demographics. When it comes to data fitting and analysis, the fractional framework offers an additional parameter that enhances the flexibility of the epidemic system. For our study on HBV infection dynamics, we have opted to employ fractional derivatives to depict the intricacies of the system. Our objective is to investigate the impact of different input factors of HBV on the system and to examine whether the fractional derivative order can be utilized as a viable control parameter or not. Therefore, we represent our model through Caputo FD as

$$\left\{\begin{array}{ccc}\hfill {}_0^{LC}{D}_t^{\beta }\mathrm{S}& =& \mathsf{\Pi }(1-v)(1-q\mathrm{C})-(f_1\mathrm{A}+f_2\mathrm{C}+f_3\mathrm{H})\frac{\mathrm{S}}{\mathrm{N}}-(g_1+d)\mathrm{S}+g_2\mathrm{V},\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{V}& =& \mathsf{\Pi }v+g_1\mathrm{S}-(d+g_2)\mathrm{V},\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{E}& =& (f_1\mathrm{A}+f_2\mathrm{C}+f_3\mathrm{H})\frac{\mathrm{S}}{\mathrm{N}}-(d+ϵ)\mathrm{E},\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{A}& =& ϵ\mathrm{E}-(d+h_1+{\delta }_1)\mathrm{A},\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{C}& =& h_1\alpha \mathrm{A}+\mathsf{\Pi }(1-v)q\mathrm{C}-(d+h_2+{\delta }_2)\mathrm{C}+\sigma \mathrm{H},\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{H}& =& {\delta }_1\mathrm{A}+{\delta }_2\mathrm{C}-(d+\varphi +\sigma )\mathrm{H},\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{R}& =& h_1(1-\alpha )\mathrm{A}+h_2\mathrm{C}+\varphi \mathrm{H}-d\mathrm{R},\hfill \end{array}\right.$$

(6)

where $\beta$ is the order of the LC derivative and $0<\beta \le 1$, moreover, system (6) becomes system (5) for $\beta =1$. The theory of fractional calculus offers a wide range of practical applications and provides a more precise representation of the dynamics of biological phenomena compared to integer-order derivatives. In the recent literature, several new fractional operators have been introduced. Our future work aims to explore the dynamics of HBV infection with the help of these novel operators, and we will further investigate and compare the outcomes produced by these fractional operators. The results presented below are on the non-negativity of the solution of the above-mentioned HBV infection model (5).

Theorem 2. 

The solutions $(\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})$ of the recommended system (6) of HBV infection are non-negative and bounded for positive initial values of state variables.

Proof of Theorem 2. 

To prove the non-negativity and boundedness of the solutions of the recommended fractional system (6) of HBV infection, we proceed as follows:

$$\left\{\begin{array}{ccc}\hfill {}_0^{LC}{D}_t^{\beta }\mathrm{S}{\mid }_{\mathrm{S}=0}& =& \mathsf{\Pi }(1-v)(1-q\mathrm{C})+g_2\mathrm{V}\ge 0,\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{V}{\mid }_{\mathrm{V}=0}& =& \mathsf{\Pi }v+g_1\mathrm{S}\ge 0,\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{E}{\mid }_{\mathrm{E}=0}& =& (f_1\mathrm{A}+f_2\mathrm{C}+f_3\mathrm{H})\frac{\mathrm{S}}{\mathrm{N}}\ge 0,\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{A}{\mid }_{\mathrm{A}=0}& =& ϵ\mathrm{E}\ge 0,\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{C}{\mid }_{\mathrm{C}=0}& =& h_1\alpha \mathrm{A}+\sigma \mathrm{H}\ge 0,\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{H}{\mid }_{\mathrm{H}=0}& =& {\delta }_1\mathrm{A}+{\delta }_2\mathrm{C}\ge 0,\hfill \\ \hfill {}_0^{LC}{D}_t^{\beta }\mathrm{R}{\mid }_{\mathrm{R}=0}& =& h_1(1-\alpha )\mathrm{A}+h_2\mathrm{C}+\varphi \mathrm{H}\ge 0,\hfill \end{array}\right.$$

(7)

thus, the solutions of our fractional model of HBV infection are non-negative. To determine the boundedness of the solutions, we add all the equations of the system and obtain the following

$${}_0^{LC}{D}_t^{\beta }(\mathrm{S}+\mathrm{V}+\mathrm{E}+\mathrm{A}+\mathrm{C}+\mathrm{H}+\mathrm{R})\le \mathsf{\Pi }-d(\mathrm{S}+\mathrm{V}+\mathrm{E}+\mathrm{A}+\mathrm{C}+\mathrm{H}+\mathrm{R}).$$

From the above, we have

$$\begin{gathered} \left(\mathrm{S}\left(t\right)+\mathrm{V}\left(t\right)+\mathrm{E}\left(t\right)+\mathrm{A}\left(t\right)+\mathrm{C}\left(t\right)+\mathrm{H}\left(t\right)+\mathrm{R}\left(t\right)\right)\le (\mathrm{S}\left(0\right)+\mathrm{V}\left(0\right)+\mathrm{E}\left(0\right)+\mathrm{A}\left(0\right)+\mathrm{C}\left(0\right)+\mathrm{H}\left(0\right)+\mathrm{R}\left(0\right)-\frac{\mathsf{\Pi }}{d}) \\ \times E_{\beta }(-d{t}^{\beta })+\frac{\mathsf{\Pi }}{d}, \end{gathered}$$

Further, we obtain the following through the properties of the Mittag–Leffler function [21]:

$$\left(\mathrm{S}\left(t\right)+\mathrm{V}\left(t\right)+\mathrm{E}\left(t\right)+\mathrm{A}\left(t\right)+\mathrm{C}\left(t\right)+\mathrm{H}\left(t\right)+\mathrm{R}\left(t\right)\right)\le \frac{\mathsf{\Pi }}{d}\doteq \mathrm{M},$$

this implies that $\left(\mathrm{S}\left(t\right)+\mathrm{V}\left(t\right)+\mathrm{E}\left(t\right)+\mathrm{A}\left(t\right)+\mathrm{C}\left(t\right)+\mathrm{H}\left(t\right)+\mathrm{R}\left(t\right)\right)\le \mathrm{M}.$ Thus, the solutions of the recommended fractional model (6) of HBV infection are non-negative and bounded.

□

Model Analysis

The steady states, reproduction number, and stability of the steady states of our model will be investigated in this subsection of the paper. Epidemic models have two meaningful steady states, namely, disease-free equilibrium (DFE) and endemic equilibrium (EE). The DFE is denoted by ${\mathcal{E}}_0$ and is given by

$${\mathcal{E}}_0(\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})=\left(\frac{\mathsf{\Pi }(d(1-v)+g_2)}{d(d+g_1+g_2)},\frac{\mathsf{\Pi }(dv+g_1)}{d(d+g_1+g_2)},0,0,0,0,0\right).$$

Here, we determine the expressions for the basic reproduction number, symbolized by ${\mathcal{R}}_0$, with the help of the next generation method [24]. In our system of HBV infection, we have four infected compartments $\mathrm{E},$$\mathrm{A},$$\mathrm{C}$, and $\mathrm{H}$, thus $m=4$, then we have

$$\mathcal{F}=\left[\begin{array}{cccc}0& \frac{{\mathrm{S}}^0{f}_1}{{\mathrm{S}}^0+{\mathrm{V}}^0}& \frac{{\mathrm{S}}^0{f}_2}{{\mathrm{S}}^0+{\mathrm{V}}^0}& \frac{{\mathrm{S}}^0{f}_3}{{\mathrm{S}}^0+{\mathrm{V}}^0}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]and\mathcal{V}=\left[\begin{array}{cccc}{q}_1& 0& 0& 0\\ -ϵ& q_2& 0& 0\\ 0& -\alpha h_1& q_3& -\sigma \\ 0& -{\delta }_1& -{\delta }_2& q_4\end{array}\right].$$

The spectral radius of $\mathcal{F}{\mathcal{V}}^{-1}$ is the reproduction number, given by

$${\mathcal{R}}_0=\frac{ϵ\{d(1-v)+g_2\}{f}_1}{q_1{q}_2(d+g_1+g_2)}+\frac{ϵ\{d(1-v)+g_2\}(\alpha h_1{q}_4+\sigma {\delta }_1)f_2}{q_1{q}_2(d+g_1+g_2)(q_3{q}_4-\sigma {\delta }_2)}+\frac{ϵ\{d(1-v)+g_2\}(\alpha h_1{\delta }_2+q_3{\delta }_1)f_3}{q_1{q}_2(d+g_1+g_2)(q_3{q}_4-\sigma {\delta }_2)}.$$

in which $q_1=d+ϵ,q_2=d+h_1+{\delta }_1,q_3=d+h_2+{\delta }_2-\mathsf{\Pi }(1-v)q$, and $q_4=d+\varphi +\sigma .$ We represent

$${\mathcal{R}}_{01}=\frac{ϵ\{d(1-v)+g_2\}{f}_1}{q_1{q}_2(d+g_1+g_2)},$$

$${\mathcal{R}}_{02}=\frac{ϵ\{d(1-v)+g_2\}(\alpha h_1{q}_4+\sigma {\delta }_1)f_2}{q_1{q}_2(d+g_1+g_2)(q_3{q}_4-\sigma {\delta }_2)},$$

and

$${\mathcal{R}}_{03}=\frac{ϵ\{d(1-v)+g_2\}(\alpha h_1{\delta }_2+q_3{\delta }_1)f_3}{q_1{q}_2(d+g_1+g_2)(q_3{q}_4-\sigma {\delta }_2)}.$$

The EE of (5) is indicated by

$${\mathcal{D}}^{★}=({\mathrm{S}}^{★}\left(t\right),{\mathrm{V}}^{★}\left(t\right),{\mathrm{E}}^{★}\left(t\right),{\mathrm{A}}^{★}\left(t\right),{\mathrm{C}}^{★}\left(t\right),{\mathrm{H}}^{★}\left(t\right),{\mathrm{R}}^{★}\left(t\right)),$$

where ${\mathrm{S}}^{★}\left(t\right),{\mathrm{V}}^{★}\left(t\right),{\mathrm{E}}^{★}\left(t\right),{\mathrm{A}}^{★}\left(t\right),{\mathrm{C}}^{★}\left(t\right),{\mathrm{H}}^{★}\left(t\right)$, and ${\mathrm{R}}^{★}\left(t\right)$ represent the endemic states of the state variables and are calculated as

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\mathrm{S}}^{★}=\hfill & \frac{(d(1-v)+g_2){\mathrm{N}}^{★}}{{\mathcal{R}}_0},\hfill \\ {\mathrm{E}}^{★}=\hfill & \frac{q_2}{ϵ}{\mathrm{A}}^{★},\hfill \\ {\mathrm{V}}^{★}=\hfill & \frac{(\mathsf{\Pi }v+g_1)(d(1-v)+g_2)N^{★}}{(d+g_2){\mathrm{R}}_0},\hfill \\ {\mathrm{C}}^{★}=\hfill & \frac{(\alpha h_1{q}_4+\sigma {\delta }_1){\mathrm{A}}^{★}}{q_3{q}_4-\sigma {\delta }_2},\hfill \\ {\mathrm{H}}^{★}=\hfill & \frac{(\alpha h_1{\delta }_2+q_3{\delta }_1){\mathrm{A}}^{★}}{q_3{q}_4-\sigma {\delta }_2},\hfill \\ {\mathrm{R}}^{★}=\hfill & \left\{\frac{(1-\alpha )v_1}{d}+\frac{(\alpha h_1{q}_4+\sigma {\delta }_1)h_2}{d(q_3{q}_4-\sigma {\delta }_2)}+\frac{(\alpha h_1{\delta }_2+q_3{\delta }_1)\varphi }{d(q_3{q}_4-\sigma {\delta }_2}\right\}{\mathrm{A}}^{★}.\hfill \end{array}\right.\end{array}$$

From system (6), we calculate the below

$${\mathrm{A}}^{★}=\frac{\mathsf{\Pi }{(q_3{q}_4-\sigma {\delta }_2)}^2{q}_1{q}_2(d+g_1+g_2)({\mathcal{R}}_0-1)ϵ}{({\mathrm{L}}_1+{\mathrm{L}}_2+{\mathrm{L}}_3+{\mathrm{L}}_4+{\mathrm{L}}_5+{\mathrm{L}}_6+{\mathrm{L}}_7)},$$

in which

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\mathrm{L}}_1=\hfill & {ϵ}^2\mathsf{\Pi }(1-v)q(d+g_2)\{(q_3{q}_4-\sigma {\delta }_2)\sigma {\delta }_1{h}_1+{\alpha }^2{h}_1^2{q}_4(q_4{f}_2+{\delta }_2{f}_3)+\sigma {\delta }_1^2(\sigma f_2+q_3{f}_3)\},\hfill \\ {\mathrm{L}}_2=\hfill & \alpha h_1{ϵ}^2\mathsf{\Pi }(1-v)q(d+g_2)\{q_4((q_3{q}_4-\sigma {\delta }_2)f_1+2\sigma {\delta }_1{f}_2)+{\delta }_1(q_3{q}_4-\sigma {\delta }_2)f_2\},\hfill \\ {\mathrm{L}}_3=\hfill & h_1ϵq_1{q}_2{q}_4(q_3{q}_4-\sigma {\delta }_2)\{(\alpha (d+h_2-\mathsf{\Pi }(1-v)q)+(1-\alpha )q_3)(d+{\gamma }_1+{\gamma }_2)+\alpha f_2(d+{\gamma }_2)\},\hfill \\ {\mathrm{L}}_4=\hfill & q_1{q}_2(q_3{q}_4-\sigma {\delta }_2){\delta }_2\{\alpha f_3(d+g_2)+(d+g_1+g_2)(\alpha q_4-\sigma )\},\hfill \\ {\mathrm{L}}_5=\hfill & q_1{q}_2{(q_3{q}_4-\sigma {\delta }_2)}^2\{ϵ(f_1(d+g_2)+(d+g_1+g_2)(d(ϵ+q_2)-q_1{q}_2)\},\hfill \\ {\mathrm{L}}_6=\hfill & q_1{q}_2(q_3{q}_4-\sigma {\delta }_2)ϵ\sigma {\delta }_1\{(d+h_2-\mathsf{\Pi }(1-v)q)(d+g_1+g_2)+(d+g_2)f_2\},\hfill \\ {\mathrm{L}}_7=\hfill & q_1{q}_2(q_3{q}_4-\sigma {\delta }_2)ϵq_3{\delta }_1\{(d+\varphi )g_1+(d+g_2)(d+\varphi +f_3)\}.\hfill \end{array}\right.\end{array}$$

In the next step, we determine the local stability of the DFE of our system of HBV infection.

Theorem 3. 

If ${\mathcal{R}}_0<1$, then ${\mathcal{E}}_0$ of the recommended HBV model presented in (5) is locally asymptotically stable $\left(LAS\right)$.

Proof of Theorem 3. 

The Jacobian matrix associated with the model, denoted by $J_{{\mathcal{E}}_0}$, evaluated around the FDE of system (6) is as follows:

$$\begin{array}{c}\hfill J_{{\mathcal{D}}_0}=\left(\begin{array}{ccccccc}-d-{\gamma }_1& 0& -\frac{(d-vd+{\gamma }_2)f_1}{d_1}& -\frac{(d-vd+g_2)f_2}{d_1}-\mathsf{\Pi }(q-vq)& -\frac{(d-vd+g_2)f_3}{d_1}& g_2& 0\\ 0& -q_1& \frac{(d(1-v)+g_2)f_1}{d+g_1+g_2}& \frac{(d(1-v)+g_2)f_2}{d+g_1+g_2}& \frac{(d(1-v)+g_2)f_3}{d+g_1+g_2}& 0& 0\\ 0& ϵ& -q_2& 0& 0& 0& 0\\ 0& 0& \alpha h_1& -q_3& \sigma & 0& 0\\ 0& 0& {\delta }_1& {\delta }_2& -q_4& 0& 0\\ g_1& 0& 0& 0& 0& -d-g_2& 0\\ 0& 0& (1-\alpha )h_1& h_2& \varphi & 0& -d\end{array}\right),\end{array}$$

where $d_1=d+g_1+g_2$ and the associated characteristic equation of $J_{{\mathcal{E}}_0}$ is constructed as

$$(\lambda +d)({\lambda }^6+{\mathrm{C}}_1{\lambda }^5+{\mathrm{C}}_2{\lambda }^4+{\mathrm{C}}_3{\lambda }^3+{\mathrm{C}}_4{\lambda }^2+{\mathrm{C}}_5\lambda +{\mathrm{C}}_6)=0.$$

(8)

Obviously, one of the eigenvalues is $-d$, which is a negative real number and obviously has a negative real part. The coefficients involved in the above Equation (8) are calculated as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\mathrm{C}}_1=\hfill & q_1+q_2+q_3+q_4+2d+g_1+g_2,\hfill \\ {\mathrm{C}}_2=\hfill & d^2+2d(q_3+q_4)+(d+q_3+q_4)(g_1+g_2)+(2d+q_3+q_4+g_1+g_2)(q_1+q_2)\hfill \\ & +(q_3{q}_4-\sigma {\delta }_2)+q_1{q}_2(1-{\mathcal{R}}_{01}),\hfill \\ {\mathrm{C}}_3=\hfill & d(d+g_1+g_2)(q_1+q_2+q_3+q_4)+(2d+g_1+g_2)(q_1{q}_3+q_1{q}_4+q_2{q}_3+q_2{q}_4)\hfill \\ & +(2d+q_1+q_2+g_1+g_2)(q_3{q}_4-\sigma {\delta }_2)+q_1{q}_2(2d+q_3+q_4+g_1+g_2)(1-{\mathcal{R}}_{01})\hfill \\ & -\frac{ϵ(d(1-v)+g_2)(\alpha h_1{f}_2+{\delta }_1{f}_3)}{d+g_1+g_2},\hfill \\ {\mathrm{C}}_4=\hfill & q_1{q}_2{q}_3{g}_1+d(d+g_1+g_2)(q_1{q}_3+q_1{q}_4+q_2{q}_3+q_2{q}_4)+(d^2+2d(q_1+q_2)\hfill \\ & +(d+q_1+q_2)(g_1+g_2))(q_3{q}_4-\sigma {\delta }_2)+\{(d^2+2d(q_3+q_4)+(d+q_4)g_1+(d+q_3+q_4)g_2)\hfill \\ & +q_1{q}_2(q_3{q}_4-\sigma {\delta }_2)\}(1-{\mathcal{R}}_{01})-\frac{ϵ(d(1-v)+g_2)(\alpha h_1{f}_2+{\delta }_1{f}_3)(d+d_1)}{d_1},\hfill \\ {\mathrm{C}}_5=\hfill & d{q}_1{q}_2\{(q_3+q_4)d_1+(q_3{q}_4-\sigma {\delta }_2)\}(1-{\mathcal{R}}_{01})+d{q}_1{q}_2(q_3{q}_4-\sigma {\delta }_2)(1-{\mathcal{R}}_0)\hfill \\ & +d\{(q_1+q_2)(d+g_1+g_2)+q_1{q}_2(g_1+g_2)\}(q_3{q}_4-\sigma {\delta }_2)+\frac{ϵ(d(1-v)+g_2)(\sigma {\delta }_2(g_1+g_2)-q_3{q}_4)f_1}{d+g_1+g_2}\hfill \\ & -ϵ(d(1-v)+g_2)\{(\alpha h_1(q_4+1)+\sigma {\delta }_1)f_2+(\alpha h_1{\delta }_2+(q_3+d){\delta }_1)f_3\},\hfill \\ {\mathrm{C}}_6=\hfill & d{q}_1{q}_2(d+g_1+g_2)(q_3{q}_4-\sigma {\delta }_2)(1-{\mathcal{R}}_0).\hfill \end{array}\right.\end{array}$$

Clearly, if ${\mathcal{R}}_0<1$, then ${\mathrm{C}}_i\in {\mathbb{R}}_{+}^{*},i=1,2,\dots ,6$. Furthermore, the other will have a negative real part if Equation (8) follows the associated Routh–Hurtwiz criterion [25]. Hence, it can be concluded from the above calculations that the DFE is $LAS$ for ${\mathcal{R}}_0<1$ if the associated Routh–Hurtwiz criterion holds.

□

4. Existence Theory

Here, the qualitative analysis of the recommended system (6) of HBV infection will be presented. We assumed the following:

$$\left\{\begin{array}{ccc}\hfill {\mathcal{Y}}_1(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})& =& \mathsf{\Pi }(1-v)(1-q\mathrm{C})-(f_1\mathrm{A}+f_2\mathrm{C}+f_3\mathrm{H})\frac{\mathrm{S}}{\mathrm{N}}-(g_1+d)\mathrm{S}+g_2\mathrm{V},\hfill \\ \hfill {\mathcal{Y}}_2(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})& =& \mathsf{\Pi }v+g_1\mathrm{S}-(d+g_2)\mathrm{V},\hfill \\ \hfill {\mathcal{Y}}_3(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})& =& (f_1\mathrm{A}+f_2\mathrm{C}+f_3\mathrm{H})\frac{\mathrm{S}}{\mathrm{N}}-(d+ϵ)\mathrm{E},\hfill \\ \hfill {\mathcal{Y}}_4(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})& =& ϵ\mathrm{E}-(d+h_1+{\delta }_1)\mathrm{A},\hfill \\ \hfill {\mathcal{Y}}_5(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})& =& h_1\alpha \mathrm{A}+\mathsf{\Pi }(1-v)q\mathrm{C}-(d+h_2+{\delta }_2)\mathrm{C}+\sigma \mathrm{H},\hfill \\ \hfill {\mathcal{Y}}_6(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})& =& {\delta }_1\mathrm{A}+{\delta }_2\mathrm{C}-(d+\varphi +\sigma )\mathrm{H},\hfill \\ \hfill {\mathcal{Y}}_7(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R})& =& h_1(1-\alpha )\mathrm{A}+h_2\mathrm{C}+\varphi \mathrm{H}-d\mathrm{R},\hfill \end{array}\right.$$

(9)

then, system (6) of HBV infection, in view of the above, can be expressed as:

$$\left\{\begin{array}{cc}\hfill {}^{LC}{D}_{0^{+}}^{\beta }\mathcal{Y}\left(r\right)=& \mathcal{P}(r,\mathcal{Y}(r\left)\right),r\in [0,\tau ],0<\beta \le 1,\\ \hfill \mathcal{Y}\left(0\right)={\mathcal{Y}}_0,\end{array}\right.$$

(10)

in which $[0,\tau ]$ is the bounded time interval and

$$\left\{\begin{array}{cc}\hfill \mathcal{Y}\left(r\right)=& \mathrm{S},\\ & \mathrm{V},\\ & \mathrm{E},\\ & \mathrm{A},\\ & \mathrm{C},\\ & \mathrm{H},\\ & \mathrm{R}.\end{array}\right.\left\{\begin{array}{cc}\hfill ,{\mathcal{Y}}_0\left(r\right)=& {\mathrm{S}}_0,\\ & {\mathrm{V}}_0,\\ & {\mathrm{E}}_0,\\ & {\mathrm{A}}_0,\\ & {\mathrm{C}}_0,\\ & {\mathrm{H}}_0,\\ & {\mathrm{R}}_0.\end{array}\right.\left\{\begin{array}{c}\hfill \mathcal{P}(r,\mathcal{Y}\left(r\right))={\mathcal{Y}}_1(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R}),\\ \hfill {\mathcal{Y}}_2(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R}),\\ \hfill {\mathcal{Y}}_3(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R}),\\ \hfill {\mathcal{Y}}_4(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R}),\\ \hfill {\mathcal{Y}}_5(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R}),\\ \hfill {\mathcal{Y}}_6(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R}),\\ \hfill {\mathcal{Y}}_7(r,\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{A},\mathrm{C},\mathrm{H},\mathrm{R}),\end{array}\right.$$

(11)

the aforementioned (10) comparable integral converts into:

$$\mathcal{Y}\left(r\right)={\mathcal{Y}}_0\left(r\right)+\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-x)}^{\beta -1}\mathcal{P}(x,\mathcal{Y}\left(x\right))dx.$$

(12)

The key steps for analyzing our suggested system are the following hypotheses:

• (A1) Constants $V_{\mathcal{P}},{\mathcal{N}}_{\mathcal{P}}$, and $q\in [0,1)$ are considered such that

  $$\left|\mathcal{P}(r,\mathcal{Y}\left(r\right))\right|\le V_{\mathcal{P}}{\left|\mathcal{Y}\right|}^q+{\mathcal{N}}_{\mathcal{P}}.$$

  (13)

  (A2) Constants $L_{\mathcal{P}}>0$, and $\mathcal{Y}$, $\stackrel{-}{\mathcal{Y}}\in X$ are assumed in a way that

  $$|\mathcal{P}(r,\mathcal{Y})-\mathcal{P}(r,\stackrel{-}{\mathcal{Y}})|\le L_{\mathcal{P}}|\mathcal{Y}-\stackrel{-}{\mathcal{Y}}|.$$

  (14)

  In this case, we define a map T on X as follows:

  $$T\mathcal{Y}\left(r\right)={\mathcal{Y}}_0\left(r\right)+\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-x)}^{\beta -1}\mathcal{P}(x,\mathcal{Y}\left(x\right))dx.$$

  (15)

If the above A1 and A2 hold true, then there is at least one solution of (10). In the next step, we will examine the solution of our recommended system.

Theorem 4. 

If A1 and A2 are satisfied, then the recommended system (6) of HBV infection has at least one solution.

Proof of Theorem 4. 

We will apply Schaefer’s fixed-point theorem to demonstrate the result. The following steps are taken to prove the theorem.

• S1: We shall demonstrate the continuity of T in the first step. For this, we assume that ${\mathcal{Y}}_j$ is continuous for $j=1,2,\dots ,9$, which also means that $\mathcal{P}(r,\mathcal{Y}(r\left)\right)$ is continuous. Here, take ${\mathcal{Y}}_i$, $\mathcal{Y}\in X$ in a way that $T{\mathcal{Y}}_i\to T\mathcal{Y}$. Further, take the following:

  $$\begin{array}{cc}\hfill ∥T{\mathcal{Y}}_i-T\mathcal{Y}∥& =\underset{r\in [0,\tau ]}{max}\left|\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-x)}^{\beta -1}{\mathcal{P}}_i\left(x,{\mathcal{Y}}_i\left(x\right)\right)dx-\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-x)}^{\beta -1}\mathcal{P}(x,\mathcal{Y}\left(x\right))dx\right|\hfill \\ \hfill & \le \underset{r\in [0,\tau ]}{max}{\int }_0^r\left|\frac{{(r-x)}^{\beta -1}}{\mathsf{\Gamma }\left(\beta \right)}\right|\left|{\mathcal{P}}_i\left(x,{\mathcal{Y}}_i\left(x\right)\right)-\mathcal{P}(x,\mathcal{Y}\left(x\right))\right|dx\hfill \\ \hfill & \le \frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}∥{\mathcal{Y}}_i-\mathcal{Y}∥\to 0\mathrm{as}i\to \infty .\hfill \end{array}$$

  (16)
• This implies that the operator T is continuous, because $T{\mathcal{Y}}_i\to T\mathcal{Y}$ is true as long as $\mathcal{P}$ is continuous.
• S2: The boundedness of the operator T will be investigated in the second step. If we select any $\mathcal{Y}\in X$, then, we have

  $$\begin{array}{cc}\hfill \parallel T\mathcal{Y}\parallel & =\underset{r\in [0,\tau ]}{max}\left|{\mathcal{Y}}_0\left(r\right)+\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-x)}^{\beta -1}\mathcal{P}(x,\mathcal{Y}\left(x\right))dx\right|\hfill \\ \hfill & \le \left|{\mathcal{Y}}_0\right|\underset{r\in [0,\tau ]}{max}\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r\left|{(r-x)}^{\beta -1}\parallel \mathcal{P}(x,\mathcal{Y}\left(x\right))\right|dx\hfill \\ \hfill & \le \left|{\mathcal{Y}}_0\right|+\frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}\left[L_{\mathcal{P}}{\parallel \mathcal{Y}\parallel }^q+N_{\mathcal{P}}\right].\hfill \end{array}$$

  (17)

Here, we will prove the boundedness of $T\left(S\right)$, where S is a bounded subset of X. Then, for any $\mathcal{Y}\in S$, we have $V\ge 0$ such that

$$\left|\right|\mathcal{Y}\left|\right|\le V,\forall \mathcal{Y}\in S.$$

(18)

Consequently, if we use the aforementioned condition to any $\mathcal{Y}\in S$, we obtain the following:

$$\left|\right|T\mathcal{Y}\left|\right|\le |{\mathcal{Y}}_0|+\frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}[V_{\mathcal{P}}{\left|\right|\mathcal{Y}\left|\right|}^q+N_{\mathcal{P}}]\le |{\mathcal{Y}}_0|+\frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}[V_{\mathcal{P}}{V}^q+N_{\mathcal{P}}].$$

(19)

This implies that $T\left(S\right)$ is bounded.

• S3: In the third stage, we select that $r_1,r_2\in [0,\tau ]$ to demonstrate the equi-continuity, and that $r_1\ge r_2$, otherwise, the following:

  $$\begin{array}{cc}\hfill \left|T\mathcal{Y}\left(r_1\right)-T\mathcal{Y}\left(r_2\right)\right|& =\left|\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{r_1}\right|{\left(r_1-x\right)}^{\beta -1}\left|\right|\mathcal{P}(x,\mathcal{Y}\left(x\right))\mid dx\hfill \\ \hfill & -\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{r_2}\left|{\left(r_2-x\right)}^{\beta -1}\parallel \mathcal{P}(x,\mathcal{Y}\left(x\right))\right|dx\mid \hfill \\ \hfill & \le \left|\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{r_1}\right|{\left(r_1-x\right)}^{\beta -1}\left|-\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{r_2}\right|{\left(r_2-x\right)}^{\beta -1}||\left|\mathcal{P}(x,\mathcal{Y}\left(x\right))\right|dx\hfill \\ \hfill & \le \frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}\left[V_{\mathcal{P}}{\parallel \mathcal{Y}\parallel }^q+N_{\mathcal{P}}\right]\left[r_1^{\beta }-r_2^{\beta }\right]\to 0\mathrm{as}{r}_1\to r_2.\hfill \end{array}$$

  (20)

Therefore, the relative compactness of $T\left(S\right)$ is proved by the Arzela–Ascoli theorem.

• S4: In the final step of the theorem, we take the following:

  $$\mathcal{D}=\{\mathcal{Y}\in X:\mathcal{Y}=\lambda T\mathcal{Y},\lambda \in (0,1\left)\right\}.$$

Here, we take $\mathcal{Y}\in \mathcal{I}$ to demonstrate the boundedness of $\mathcal{D}$. Then, for each r in the range of $r\in [0,\tau ]$, the following holds true:

$$\left|\right|\mathcal{Y}\left|\right|=\lambda \left|\right|T\mathcal{Y}\left|\right|\le \lambda \left[|{\mathcal{Y}}_0|\frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}[V_{\mathcal{P}}{\left|\right|\mathcal{Y}\left|\right|}^q+N_{\mathcal{P}}]\right].$$

(21)

This suggests that the set $\mathcal{D}$ is bounded. Schaefer’s theorem shows that the operator T has a fixed point as a result, and consequently, our recommended model (6) of HBV infection has at least one solution.

□

Remark 1. 

In case ($A1$) holds for $q=1$, the result of Theorem (4) is fulfilled for $\frac{{\tau }^{\beta }{V}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}<1$.

Theorem 5. 

The recommended fractional system (10) of HBV infection has a unique solution if $\frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}<1$ is satisfied.

Proof of Theorem 5. 

Applying Banach’s contraction theorem, and assuming $\mathcal{Y},\stackrel{-}{\mathcal{Y}}$ to obtain the desired result, we have

$$\begin{array}{ccc}\hfill \left|\right|T\mathcal{Y}-T\stackrel{-}{\mathcal{Y}}\left|\right|& \le & \underset{r\in [0,\tau ]}{max}\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{|(r-x)}^{\beta -1}\left|\right|\mathcal{P}(x,\mathcal{Y}\left(x\right))-\mathcal{P}(x,\stackrel{-}{\mathcal{Y}}\left(x\right))|dx,\hfill \end{array}$$

(22)

also

$$\underset{r\in [0,\tau ]}{max}\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{|(r-x)}^{\beta -1}\left|\right|\mathcal{P}(x,\mathcal{Y}\left(x\right))-\mathcal{P}(x,\stackrel{-}{\mathcal{Y}}\left(x\right))|dx\le \frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}\left|\right|\mathcal{Y}-\stackrel{-}{\mathcal{Y}}\left|\right|,$$

this implies that

$$\left|\right|T\mathcal{Y}-T\stackrel{-}{\mathcal{Y}}\left|\right|\le \frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}\left|\right|\mathcal{Y}-\stackrel{-}{\mathcal{Y}}\left|\right|.$$

Thus, T has a fixed point. This implies that our system (10) of HBV has a unique solution.

□

5. Ulam–Hyers Stability

Here, the Ulam–Hyers stability for the recommended fractional system of HBV infection will be examined. The idea of this stability was initially presented by Ulam in 1940; Hyers developed it [26,27]. Many researchers have embraced the Ulam–Hyers stability hypothesis in a variety of disciplines of study [28,29]. Its basic ideas and definitions include the following.

Assume that the operator $\mathcal{H}:X\to X$ behaves in a way that

$$\mathcal{HY}=\mathcal{Y}for\mathcal{Y}\in X.$$

(23)

Definition 3. 

The above (23) is Ulam–Hyers stable (UHS) if for $ϵ>0$, and takes any solution $\mathcal{Y}\in X$ of

$$\left|\right|\mathcal{Y}-\mathcal{HY}\left|\right|\le ϵ,forr\in [0,\tau ].$$

(24)

Then, there exists a unique solution $\stackrel{-}{\mathcal{Y}}$ of (23), with $C_q>0$ fulfilling the following

$$\left|\right|\stackrel{-}{\mathcal{Y}}-\mathcal{Y}\left|\right|\le C_qϵ,r\in [0,\tau ].$$

(25)

Definition 4. 

The aforementioned (23) is a generalized Ulam–Hyers stable (GUHS), if every solution $\mathcal{Y}$ of (24) and any other $\stackrel{-}{\mathcal{Y}}$ of (23) satisfies

$$\left|\right|\stackrel{-}{\mathcal{Y}}-\mathcal{Y}\left|\right|\le \mathcal{U}\left(ϵ\right),$$

(26)

where zero is the image of zero, i.e., $\mathcal{U}\left(0\right)=0$ and $\mathcal{U}\in C(R,R)$.

Remark 2. 

The solution $\stackrel{-}{\mathcal{Y}}\in X$ fulfills (24) if

• (1) $\left|{\rm Y}\right(r\left)\right|\le ϵ,\forall r\in [0,\tau ]$ and ${\rm Y}\in C\left(\right[0,\tau ];R)$,
• (2) ${}^{LC}{D}_{0^{+}}^{\beta }\stackrel{-}{\mathcal{Y}}\left(r\right)=\mathcal{P}(r,\stackrel{-}{\mathcal{Y}}\left(r\right))+{\rm Y}\left(r\right).$

After a minor perturbation, system (10) can be expressed as follows:

$$\left\{\begin{array}{cc}\hfill {}^{LC}{D}_{0^{+}}^{\beta }\mathcal{Y}\left(r\right)=& \mathcal{P}(r,\mathcal{Y}(r\left)\right)+{\rm Y}\left(r\right),\\ \hfill \mathcal{Y}\left(0\right)={\mathcal{Y}}_0.\end{array}\right.$$

(27)

Lemma 2. 

The above system (27) satisfies the condition

$$|\mathcal{Y}\left(r\right)-T\mathcal{Y}\left(r\right)|\le c\mathsf{\Omega }\left(t\right)ϵ,wherec=\frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}.$$

(28)

Proof. 

From the second part of Remark 2 for $\mathcal{Y}\in X$, we have

$${}^{LC}{D}_{0^{+}}^{\beta }\mathcal{Y}\left(r\right)=\mathcal{P}(r,\mathcal{Y}\left(r\right))+{\rm Y}\left(r\right),$$

applying fractional theory [21], we have

$$\begin{array}{ccc}\hfill \mathcal{Y}\left(r\right)& =& {\mathcal{Y}}_0+\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-\tau )}^{\beta -1}\mathcal{P}(\tau ,\mathcal{Y}\left(\tau \right))d\tau \hfill \\ \\ & +& \frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-\tau )}^{\beta -1}{\rm Y}\left(\tau \right)d\tau ,\hfill \end{array}$$

simplifying, we have

$$\begin{array}{ccc}\hfill \left|\mathcal{Y}\right(r)& -& {\mathcal{Y}}_0-\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-\tau )}^{\beta -1}\mathcal{P}(\tau ,\mathcal{Y}\left(\tau \right))d\tau |\hfill \\ \\ & \le & \frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^r{(r-\tau )}^{\beta -1}\left|{\rm Y}\left(\tau \right)\right|d\tau .\hfill \end{array}$$

Using the definition of $T\mathcal{Y}\left(r\right)$ and the first part of the above Remark 2, we have

$$|\mathcal{Y}\left(r\right)-T\mathcal{Y}\left(r\right)|\le c\mathsf{\Omega }\left(t\right)ϵ,wherec=\frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}.$$

Thus, the required result is proved. □

Theorem 6. 

The solution of system (10) of HBV infection is UHS and GUHS on Lemma 2 if $\frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}<1$ holds true.

Proof of Theorem 6. 

To prove that the solution of system (10) of HBV infection is UHS and GUHS, we assume that $\mathcal{Y}\in X$ is any solution and $\stackrel{-}{\mathcal{Y}}\in X$ is a unique solution of system (10), then

$$\begin{array}{ccc}\hfill |\mathcal{Y}\left(r\right)-\stackrel{-}{\mathcal{Y}}\left(r\right)|& =& |\mathcal{Y}\left(r\right)-T\stackrel{-}{\mathcal{Y}}\left(r\right)|\hfill \\ \\ & \le & |\mathcal{Y}\left(r\right)-T\mathcal{Y}\left(r\right)|+|T\mathcal{Y}\left(r\right)-T\stackrel{-}{\mathcal{Y}}\left(r\right)|.\hfill \end{array}$$

Applying Lemma 2, we have

$$\begin{array}{ccc}\hfill |\mathcal{Y}\left(r\right)-\stackrel{-}{\mathcal{Y}}\left(r\right)|& \le & cϵ+\frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}|\mathcal{Y}\left(r\right)-\stackrel{-}{\mathcal{Y}}\left(r\right)|\hfill \\ & \le & \frac{cϵ}{1-\frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}}.\hfill \end{array}$$

(29)

This implies that the solution of system (10) of HBV infection is both UHS and GUHS.

□

Definition 5. 

Assume that $\mathsf{\Omega }\in C\left(\right[0,\tau ],R)$, then (23) is Ulam–Hyers–Rassias stable (UHRS) if and only if any solution $\mathcal{Y}\in X$ satisfies the below

$$\left|\right|\mathcal{Y}-\mathcal{HY}\left|\right|\le \mathsf{\Omega }\left(r\right)ϵ,forr\in [0,\tau ]andϵ>0.$$

(30)

Then, system (23) has a unique solution $\stackrel{-}{\mathcal{Y}}$ with $C_q>0$ and fulfilling

$$\left|\right|\stackrel{-}{\mathcal{Y}}-\mathcal{Y}\left|\right|\le C_q\mathsf{\Omega }\left(r\right)ϵ,\forall r\in [0,\tau ].$$

(31)

Definition 6. 

Let $\mathcal{Y}$ be any solution and $\stackrel{-}{\mathcal{Y}}$ be a unique solution of (23), moreover, if there exists $C_{q,\mathsf{\Omega }}$ for $\mathsf{\Omega }\in C\left(\right[0,\tau ],R)$ with $ϵ>0$ in a manner that

$$\left|\right|\stackrel{-}{\mathcal{Y}}-\mathcal{Y}\left|\right|\le C_{q,\mathsf{\Omega }}\mathsf{\Omega }\left(r\right)ϵ,\forall r\in [0,\tau ].$$

(32)

Then, the above (23) is generalized Ulam–Hyers–Rassias stable (GUHRS).

Remark 3. 

The solution $\stackrel{-}{\mathcal{Y}}\in X$ satisfies (24) and the below

• (a) $\left|{\rm Y}\right(r\left)\right|\le ϵ\mathsf{\Omega }\left(r\right),\forall r\in [0,\tau ]$, in which ${\rm Y}\left(r\right)\in C\left(\right[0,\tau ];R)$,
• (b) $\mathcal{H}\stackrel{-}{\mathcal{Y}}\left(r\right)=\stackrel{-}{\mathcal{Y}}+{\rm Y}\left(r\right),\forall r\in [0,\tau ].$

Lemma 3. 

The perturbed system (10) satisfies the inequality below

$$|\mathcal{Y}\left(r\right)-T\mathcal{Y}\left(r\right)|\le c\mathsf{\Omega }\left(r\right)ϵ,wherec=\frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}.$$

(33)

The use of Lemma (1) and Remark (3) makes it simple to demonstrate the above result.

Theorem 7. 

The system solution (10) is UHR stable and GUHR stable on Lemma (3) if $\frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}<1$.

Proof of Theorem 7. 

To prove the result, let $\stackrel{-}{\mathcal{Y}}\in X$ be a unique solution and $\mathcal{Y}\in X$ be any other solution of system (10), then we have

$$\begin{array}{ccc}\hfill |\mathcal{Y}\left(r\right)-\stackrel{-}{\mathcal{Y}}\left(r\right)|& =& |\mathcal{Y}\left(r\right)-T\stackrel{-}{\mathcal{Y}}\left(r\right)|\hfill \\ & \le & |\mathcal{Y}\left(r\right)-T\mathcal{Y}\left(r\right)|+|T\mathcal{Y}\left(r\right)-T\stackrel{-}{\mathcal{Y}}\left(r\right)|\hfill \\ & \le & a\mathsf{\Omega }\left(r\right)ϵ+\frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}|\mathcal{Y}\left(r\right)-\stackrel{-}{\mathcal{Y}}\left(r\right)|\hfill \\ & \le & \frac{c\mathsf{\Omega }\left(r\right)ϵ}{1-\frac{{\tau }^{\beta }{L}_{\mathcal{P}}}{\mathsf{\Gamma }(\beta +1)}},\hfill \end{array}$$

(34)

where a is a positive constant. Thus, the solution of system (10) is therefore both GUHRS and UHRS. □

6. Numerical Results and Discussions

In this section of the paper, we will find numerical results to illustrate the dynamical behavior of the recommended fractional model (6) of HBV infection. In the literature, several numerical methods have been introduced for the analysis of an LC fractional derivative. We will use the recently developed numerical scheme presented below

$$\begin{array}{c}\hfill {}_0^{LC}{D}_t^{\beta }\mathrm{x}\left(t\right)=\hslash (t,\mathrm{x}\left(t\right)).\end{array}$$

(35)

The above implies that

$$\begin{array}{c}\hfill \mathrm{x}\left(t\right)-\mathrm{x}\left(0\right)=\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^t\hslash (\rho ,\mathrm{x}\left(\rho \right)){(t-\rho )}^{\beta -1}d\rho ,\end{array}$$

(36)

for $t=t_{n+1}$ with $n=0,1,\dots$, we have

$$\begin{array}{c}\hfill \mathrm{x}\left(t_{n+1}\right)-\mathrm{x}\left(0\right)=\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_{n+1}}{(t_{n+1}-t)}^{\beta -1}\hslash (t,\mathrm{x}\left(t\right))dt,\end{array}$$

(37)

and

$$\begin{array}{c}\hfill \mathrm{x}\left(t_n\right)-\mathrm{x}\left(0\right)=\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_n}{(t_n-t)}^{\beta -1}\hslash (t,\mathrm{x}\left(t\right))dt.\end{array}$$

(38)

The above (38) and (37) implies that

$$\begin{array}{ccc}\hfill \mathrm{x}\left(t_{n+1}\right)& =& \mathrm{x}\left(t_n\right)+\underset{{\mathcal{B}}_{\beta ,1}}{\underset{.}{\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_{n+1}}{(t_{n+1}-t)}^{\beta -1}\hslash (t,\mathrm{x}\left(t\right))dt}}\hfill \\ & & -\underset{{\mathcal{B}}_{\beta ,2}}{\underset{.}{\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_n}{(t_n-t)}^{\beta -1}\hslash (t,\mathrm{x}\left(t\right))dt}}.\hfill \end{array}$$

(39)

in which

$$\begin{array}{c}\hfill {\mathcal{B}}_{\beta ,1}=\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_{n+1}}{(t_{n+1}-t)}^{\beta -1}\hslash (t,\mathrm{x}\left(t\right))dt,\end{array}$$

(40)

and

$$\begin{array}{c}\hfill {\mathcal{B}}_{\beta ,2}=\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_n}{(t_n-t)}^{\beta -1}\hslash (t,\mathrm{x}\left(t\right))dt.\end{array}$$

(41)

Here, we apply the well-known Lagrange approximation for $\hslash (t,\mathrm{x}(t\left)\right)$ as follows:

$$\begin{array}{ccc}\hfill \mathcal{L}\left(t\right)& \simeq & \hslash (t_n,{\mathrm{x}}_n)\frac{t-t_{n-1}}{t_n-t_{n-1}}+\hslash (t_{n-1},{\mathrm{x}}_{n-1})\frac{t-t_n}{t_{n-1}-t_n}\hfill \\ & =& \frac{\hslash (t_n,{\mathrm{x}}_n)}{h}(t-t_{n-1})-\frac{\hslash (t_{n-1},{\mathrm{x}}_{n-1})}{h}(t-t_n).\hfill \end{array}$$

(42)

The following is obtained with the help of above

$$\begin{array}{ccc}\hfill {\mathcal{B}}_{\beta ,1}& =& \frac{\hslash (t_n,{\mathrm{x}}_n)}{h\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_{n+1}}{(t_{n+1}-t)}^{\beta -1}(t-t_{n-1})dt\hfill \\ & -& \frac{\hslash (t_{n-1},{\mathrm{x}}_{n-1})}{h\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_{n+1}}{(t_{n+1}-t)}^{\beta -1}(t-t_n)dt,\hfill \end{array}$$

(43)

further, we have

$$\begin{array}{ccc}\hfill {\mathcal{B}}_{\beta ,1}& =& \frac{\hslash (t_n,{\mathrm{x}}_n)}{h\mathsf{\Gamma }\left(\beta \right)}\left[\frac{2h}{\beta }{t}_{n+1}^{\beta }-\frac{t_{n+1}^{\beta +1}}{\beta +1}\right]\hfill \\ & -& \frac{\hslash (t_{n-1},{\mathrm{x}}_{n-1})}{h\mathsf{\Gamma }\left(\beta \right)}\left[\frac{h}{\beta }{t}_{n+1}^{\beta }-\frac{1}{\beta +1}{t}_{n+1}^{\beta +1}\right].\hfill \end{array}$$

(44)

Following the same steps, we obtain

$$\begin{array}{ccc}\hfill {\mathcal{B}}_{\beta ,2}& =& \frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^{t_n}{(t_n-t)}^{\beta -1}[\frac{\hslash (t_n,{\mathrm{x}}_n)}{h}(t-t_{n-1})\hfill \\ & -& \frac{\hslash (t_{n-1},{\mathrm{x}}_{n-1})}{h}(t-t_n)]dt.\hfill \end{array}$$

(45)

The following is obtained through simplification

$$\begin{array}{ccc}\hfill {\mathcal{B}}_{\beta ,2}& =& \frac{\hslash (t_n,{\mathrm{x}}_n)}{h\mathsf{\Gamma }\left(\beta \right)}\left[\frac{h}{\beta }{t}_n^{\beta }-\frac{t_n^{\beta +1}}{\beta +1}\right]\hfill \\ & +& \frac{\hslash (t_{n-1},{\mathrm{x}}_{n-1})}{h\mathsf{\Gamma }\left(\beta \right)}\left[\frac{1}{\beta +1}{t}_n^{\beta +1}\right].\hfill \end{array}$$

(46)

By substituting Equations (45) and (46) into (39), we have obtained the ultimate approximate solution for the fractional system

$$\begin{array}{ccc}\hfill \mathrm{x}\left(t_{n+1}\right)& =& \mathrm{x}\left(t_n\right)+\frac{\hslash (t_n,{\mathrm{x}}_n)}{h\mathsf{\Gamma }\left(\beta \right)}\left[\frac{2h{t}_{n+1}^{\beta }}{\beta }-\frac{t_{n+1}^{\beta +1}}{\beta +1}+\frac{h}{\beta }{t}_n^{\beta }-\frac{t_{n+1}^{\beta +1}}{\beta +1}\right]\hfill \\ & +& \frac{\hslash (t_{n-1},{\mathrm{x}}_{n-1})}{h\mathsf{\Gamma }\left(\beta \right)}\left[-\frac{h}{\beta }{t}_{n+1}^{\beta }+\frac{t_{n+1}^{\beta +1}}{\beta +1}+\frac{t_n^{\beta +1}}{\beta +1}\right].\hfill \end{array}$$

(47)

This method has been recently introduced in the literature [30] and has been effectively utilized in different research work. In this work, our main concern is to represent the dynamical behavior of the recommended system of HBV infection with this numerical scheme. However, other aspects of the method will be discussed in future work. Understanding the dynamical behavior of an epidemic model is crucial for predicting the course of an outbreak, evaluating the effectiveness of interventions, and informing public health policies. By analyzing the model, conducting simulations, and fitting the model to real-world data, researchers can gain insights into the key factors driving the epidemic’s dynamics and make informed decisions to mitigate its impact. For simulation purposes, the state variables are assumed to be $\mathrm{S}\left(0\right)=0.8,\mathrm{V}\left(0\right)=0.3,\mathrm{E}\left(0\right)=0.3,\mathrm{A}\left(0\right)=0.3,\mathrm{C}\left(0\right)=0.3,\mathrm{H}\left(0\right)=0.3$, and $\mathrm{R}\left(0\right)=0.3$.

In Figure 1, we illustrate the impact of the fractional parameter $\beta$ on the solution pathways of the system of HBV infection to check whether this parameter of the system can be used as a preventive measure or not. The values of $\beta$ are chosen to be $0.7,$$0.8,$$0.9$, and $1.0$ in the first simulation. It can be seen that the infection level can be decreased by decreasing the value of the input parameter $\beta$. We noticed that this input factor contributes significantly and can be used as a control parameter, therefore, suggested this to the policymakers. In Figure 1, a comparative analysis of integer and non-integer derivatives are also illustrated. The curve for $\beta =1$ represents the integer case, while curves with other values of $\beta$ demonstrate non-integer cases. The effect of vaccination on the dynamics of HBV infection is shown in Figure 2. In this simulation, we assumed different values of the vaccination factor v, i.e., $v=0.45,$$0.55,$$0.65$, and $0.75$. We observed that an increase in v decreases the exposed and chronically infected individuals of the system. We recommended increasing the efficiency of vaccination for better control of the disease. In Figure 3, the effect of the loss rate of immunity due to vaccination is illustrated on the vaccinated and exposed individuals of the system. It can be seen that this factor is dangerous and can increase the risk of infection in society.

To visualize the impact of $h_1$ and $h_2$ on the system, we performed two different simulations, shown in Figure 4 and Figure 5. In Figure 4, we take $h_1=0.019,$$0.039,$$0.059$, and $0.079$ while the values of $h_2$ are assumed to be $0.025,$$0.050,$$0.075$, and $0.100$ in FigureFigure 5. It is clear from the figures how the infected individuals behave with the variation of these parameters. In the final simulation, presented in Figure 6, the solution pathways of the system are highlighted with the variation in transmission rate $f_1$. In Figure 6, we assumed the values of $f_1$ to be $0.033$, $0.063$, $0.093$, and $0.123$ to identify the role of $f_1$ on the dynamics of HBV infection. We noticed that this parameter is critical and highly increases the endemic level of the infection in society. This implies that the transmission rate $f_1$ is sensitive and needs to be controlled.

It is evident that HBV is a dangerous infectious disease, therefore, effective interventions are needed to eradicate and eliminate this viral infection from society. The control of the transmission rate and fractional parameter will reduce the intensity of the infection. Therefore, the control of these factors is recommended to officials. We also suggest increasing the efficiency of vaccination and medication through hospitalization for better prevention of the infection.

7. Conclusions

In this work, we constructed an epidemic model for HBV infection with vaccination and medication through hospitalization. The model is structured in a fractional framework for more accurate outcomes. We examined the steady states and determined the threshold parameter, indicated by ${\mathcal{R}}_0$. The stability result of the infection-free steady state was established. The existence and uniqueness of the hypothesized system’s solution are investigated using the fixed-point theorem in the context of Banach’s and Schaefer’s theorems. In our HBV system, we established the essential conditions for Ulam–Hyers stability. A time-series analysis of the system is presented through a numerical method. We have shown the impact of different parameters on the solution pathways of the system and recommended the most critical scenario of the system for the prevention and control of the infection. This model can assist in saving lives and prevent the spread of the infection by offering valuable information on the transmission dynamics of HBV and the efficacy of diverse intervention strategies. HBV infection is often associated with other comorbidities, such as HIV, hepatitis C, and liver cancer. These comorbidities can affect the course of HBV infection and complicate the design of effective control strategies. In future work, we will explore how comorbidities interact with HBV infection and how they can be incorporated into the model to better understand the overall disease burden and design more effective interventions.