A New Modeling of Fractional-Order and Sensitivity Analysis for Hepatitis-B Disease with Real Data

Abstract

In this study, we propose new illustrative and effective modeling to point out the behaviors of the Hepatitis-B virus (Hepatitis-B). Not only do we consider the mathematical modeling, equilibria, stabilities, and existence–uniqueness analysis of the model, but also, we make numerical simulations by using the Adams–Bashforth numerical scheme. However, we apply the parameter estimation method to determine our model parameters and find the curve that best fits the model. Additionally, in this study, the stability analysis of the aforementioned model is considered, and also the sensitivity analysis of R₀ is examined. The results point out that the order of the fractional derivative has an essential effect on the dynamical process of the constructed model for Hepatitis-B.

1. Introduction

Recently, studies of infectious disease research have increased in the literature. It is increasingly difficult to control the spread of infectious diseases among humans. A notable disease caused by liver infection is an infection caused by the Hepatitis-B virus, which is widespread worldwide. In some people, Hepatitis-B can survive in the body after infection without causing illness. In some people, although the disease is cured, the virus cannot be eliminated from the body. Although the virus does not cause any symptoms in these people, the person remains a carrier, and the risk of transmitting the infection continues. Hepatitis-B disease is divided into acute and chronic. Acute Hepatitis-B is a short-term disease that occurs within the first 6 months after exposure to the Hepatitis-B virus. Acute Hepatitis-B disease causes symptoms such as fever, fatigue, loss of appetite, nausea or vomiting, jaundice (yellow skin or eye color, dark urine, and clay-colored stools), and muscle, joint, and stomach pain. Chronic Hepatitis-B, on the other hand, is the form of Hepatitis-B virus that cannot enter the body for more than 6 months and can cause bad consequences, such as liver damage (cirrhosis), liver cancer, and death [1]. There is no known cure for acute Hepatitis-B. However, chronic Hepatitis-B can be controlled using medication.

Hepatitis-B infection is transmitted from person to person both horizontally and vertically. Being contagious from one person to another, it can be re-transmitted through blood or water via sexual contact or the reuse of unsafe syringes or needles. This type of transmission is known as horizontal transmission. A mother infected with Hepatitis-B also transmits the virus to her newborn baby, and this is called vertical transmission. In the acute stage, most people cannot show any symptoms. Common ways of transmission of Hepatitis-B include percutaneous (parenteral) transmission, early childhood infection, sexual transmission, contaminated water or food, blood transfusion, unsafe infection, occupational exposure of healthcare workers, and other transmission routes and risk groups [2].

When the history of medicine is examined, it is seen that jaundice diseases are mentioned for the first time in the documents belonging to ancient Greece and the Romans. Viral hepatitis was first described by Hippocrates as epidemic (infectious) jaundice, described in the 5th century BC. While probably most of these are due to the Hepatitis-A virus, the epidemic transmission of Hepatitis-B has begun to be observed in places where the use of blood and blood products is common [3]. The hepatitis form transmitted directly with blood and blood products was first described by Lurman in 1883. The importance of Hepatitis-B infection to society was first revealed during a vaccination campaign against smallpox in Bremen, Germany [4,5]. In 1947, the terms Hepatitis-A and -B were introduced by Mccallum and Bauer to distinguish between infectious (epidemic) and serum hepatitis [6]. These terms were adopted in 1973 by the scientific group of the World Health Organization (WHO) working on viral hepatitis. At the National Institutes of Health (USE), Blumberg et al. showed that the serum of an Australian person has a precipitating antigen on agar gel with the serum of a patient with multiple blood transfusions; this protein is now known as “Hepatitis-B surface antigen HBsAg”, and they called the protein “Australian antigen-Au antigen”. Over time, it became clear that this protein is associated with type B hepatitis. In 1968, researchers led by Prince, Okochi, and Murakam showed that Au antigen (hepatitis B surface antigen) was found only in the serum of patients infected with type B hepatitis [7]. Viral hepatitis is a serious public health infection that is common all over the world and is very closely related to the economies of the country. July 28 is designated as World Hepatitis Day to draw attention to hepatitis, as the disease throughout the world shows symptoms in the late period and the majority of those infected are not aware of their diseases. According to World Health Organization (WHO) data, more than 250 million people globally live with chronic Hepatitis-B infection, and cost the lives of approximately 900,000 people each year. Therefore, the Hepatitis-B epidemic continues to threaten public health globally [8]. The WHO Western Pacific Region and WHO African Region have the highest Hepatitis-B prevalence rates, with 6.2% and 6.1% of the adult population afflicted, respectively. According to estimates, 3.3%, 2.0%, and 1.63% of the general population in the WHO Eastern Mediterranean Region, WHO Southeast Asia Region, and WHO European Region, respectively, are infected. In the WHO Americas Region, 0.7% of the population is infected [8].

WHO divided countries into low ($<2$%), medium (2–8%) and high ($>8$%) endemic regions in terms of Hepatitis-B infection carrier rates [9]. The Hepatitis-B virus can last for at least 7 days outside the body. During this time, the virus can still cause infection if it enters the body of a person who has not received the vaccine. Hepatitis-B virus incubation lasts, on average, 75 days, but it can last anywhere between 30 and 180 days. The virus can persist and cause chronic Hepatitis-B, and it can be found within 30 to 60 days of infection [8]. The Hepatitis-B virus (Hepatitis-B), which has a circular genome that is partially composed of double-stranded DNA, reproduces through an RNA intermediate form by reverse copying [9,10,11,12]. In other words, although Hepatitis-B is a DNA virus, it encodes the “reverse transcriptase” enzyme and, thanks to this enzyme, replicates through the RNA mediator. It is found in the nucleus of an infected cell as a mini chromosome. Replication and all transcriptions occur through a DNA chain that has the feature of a mediator molecule, called covalently linked circular DNA [13]. For chronic Hepatitis-B carriers, there is no readily accessible, efficient treatment. The most crucial preventative step is Hepatitis-B vaccination [14]. Based on the utilization of a viral envelope protein, several vaccines have been created to prevent Hepatitis-B infection (HBsAg). From infants to adults, the vaccination offers protection to 85–90% of those who receive it [14]. Base vaccines consist of the three-dose HB vaccine and the HB birth dose (that is, within 24 h after birth). WHO recommended in 1991 that Hepatitis-B vaccination be a part of the national immunization program in all nations with HBsAg carrier prevalence rates of 8% or above by 1995 and in all nations by 1997. By 2002, 154 nations had routine HB immunization programs in place [15].

Epidemiological study plays an important role in understanding the impact of infectious diseases on society. Understanding the mechanics of disease can be facilitated by mathematical models. Mathematical modeling is widely utilized in the social sciences (including economics, psychology, sociology, and political science), in the natural sciences (including physics, biology, earth science, meteorology), and engineering fields (including computer science and artificial intelligence). Additionally, to explain and analyze nonlinear processes, mathematicians are frequently used by physicists, engineers, statisticians, operations research analysts, and economists [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]. The model can help explain a system, examine the effects of different components, and make predictions about the behavior of a real-life problem. In mathematical modeling, we investigate models by creating models, estimating parameters, checking the precision of models with variable parameters, and calculating numerical simulations. This type of research helps to understand the spread of disease in the population and to control its parameters. These types of disease patterns are often referred to as infectious diseases (disease passed from one person to another) [9,10]. Many researchers have studied the directions of the propagation of Hepatitis-B in various regions and the immune response during infection with mathematical models. Anderson and May illustrated the consequences of the transmission of Hepatitis-B in carriers using a straightforward mathematical model [16,17]. A variable combination of sexual activity and age is included in Anderson’s and Williams’ descriptions of models of sexual transmission of Hepatitis-B [18]. The link between the age of Hepatitis-B infection and the emergence of the carrier state was demonstrated by Edmunds et al. [19]. A feedback mechanism model developed by Medley et al. connects the likelihood of carrier-class creation following infection to the rate of transmission, the average age of infection, and the prevalence of infection [20]. To forecast chronic hepatitis-B infection in New Zealand, Thornley et al. used the Medley model [21] and Din and Abidin [22] modeled a vaccinated Hepatitis-B epidemic with the Mittag–Leffler kernel.

The spread of Hepatitis-B differs from that in industrialized countries in developing nations. Edmunds et al. modeled the transmission rate of childhood Hepatitis-B as the main determinant of epidemic level and sexual contact rates in developing countries [23]. Little is known about sexual contact rates in developing countries. McLean, Blumberg, and Edmunds studied Hepatitis-B transmission patterns in developing countries, and Williams described a Hepatitis-B model in the UK [17,24,25]. To study acute Hepatitis-B infection, the timing of effector cell activation and progression, and the function of pre-existing or vaccine-induced antibodies in preventing Hepatitis-B infection, Coupe et al. developed a typical immune response model [26,27,28]. While Gourley et al. produced a time-delayed extension of this model, Min et al. employed a conventional model function rather than a mass action to explain vulnerability to Hepatitis-B infection [29,30]. For the Hepatitis-B population, Hews et al. applied a logistical growth and standard model to improve the model’s representation of the available data and achieve more accurate baseline reproduction number estimates [34].

It has been proved that mathematical models created by using ordinary differential equations of integer order are not sufficient in understanding the propagation of biological systems. Using fractional-order mathematical models for biological systems gives closer results because it gives more comprehensive results than integer-order mathematical models. The use of fractional derivative and partial differential equations, especially fractional-order derivatives, which are included in important applications of mathematics that are intertwined with real life, has become widespread. Recently, many scientists have turned their studies to real-life problems [35,36,37,38,39]. This trend ranges from new model structures to real-life problems and the application of fractional-order differential equations. The orders of differential equations specifying physical events determine the rate of change of the physical event in question. Additionally, fractional-order differential equations play an important role in filling some deficiencies in explaining some nonlinear phenomena and in understanding physical phenomena [40,41,42,44,45,46]. The concepts of the fractional derivative and fractional integral have many applications in applied sciences, engineering, finance, geology, thermal sciences, seismology, fluid fluids, elastic theory, thermodynamics, and hydrodynamics [47,48,49,50,51,52]. On the other hand, the earliest systematic studies in this area were performed by Liouville, Riemann, and Holmgren in the 19th century [40]. Apart from Leibniz, many famous mathematicians, such as Riemann, Euler, Liouville, Laplace, and Fourier, have also worked on the fractional analysis. Different fractional derivatives, such as Riemann–Liouville, Caputo, and Grünwald–Letnikov, are used according to the properties of these processes and systems in obtaining the fractional mathematical models of physical processes and systems because each fractional derivative has different properties and advantages over the others [53]. The integer-order derivative provides local modeling, and the fractional-order derivative provides global modeling. Therefore, different studies have been performed by using differential equation models that include fractional derivatives [54,55,56]. The commonly used derivatives related to fractional derivatives in literature are Caputo and Riemann–Liouville fractional derivatives. These derivative definitions have been effectively applied to many real-life problems, such as substance transport, fluid mechanics, population models, control systems, and financial changes. The motivation for this study came from the aforementioned conversation. Our objective is to analyze a fractional-order mathematical model and support it by demonstrating that it fits the real data more closely than the integer-order model. In this context, we propose a fractional-order mathematical model that describes the dynamics of Hepatitis-B in Section 3. Although fractional-order derivative definitions are widely used in the literature, there are some shortcomings arising from definitions. Caputo–Fabrizio (2015) and Atangana–Baleanu (2016) derivative definitions were introduced to eliminate the inadequate aspects of fractional-order derivatives [57,58,59,60,61,62,63,64,65]. These derivative definitions were created by modifying the core function.

The sections in the paper are arranged as follows: In Section 2, we present the most important definitions of fractional calculus. The fractional model is given in Section 3. The existence and uniqueness of the model solution, positivity and boundedness of the solution and the stability of the equilibrium points of the proposed model are given in Section 4. Additionally, the basic reproduction number is calculated. The sensitivity analysis is investigated according to the parameters to study the extent of its effect on the reproduction number in Section 5. The parameter estimation method is studied in Section 6. The numerical method of the proposed model is given in Section 7. In Section 8, the numerical solutions for our model are presented by using the Adams–Bashforth Moulton method for the fitted parameter values in Table 1. Finally, a summary of the present work is presented in Section 9.

2. Preliminary Results

This section will focus on some and basic concepts that will be used throughout the study.

Definition 1.

The one-parameter Mittag–Leffler function $E_{\alpha }\left(y\right)$ is defined by [58]

$$E_{\alpha }\left(y\right)=\sum _{k=0}^{\infty }\frac{y^k}{\Gamma (\alpha k+1)},(y\in \mathbb{C}).$$

(1)

Definition 2.

The two-parameter Mittag–Leffler function $E_{\alpha ,\beta }\left(y\right)$ is defined by [58]

$$E_{\alpha ,\beta }\left(y\right)=\sum _{k=0}^{\infty }\frac{y^k}{\Gamma (\alpha k+\beta )},(y\in \mathbb{C}).$$

(2)

where $\alpha ,\beta >0$.

Definition 3.

The Riemann–Liouville form of fractional integral operator of order $\alpha >0$ of a function $f:(0,\infty )\to \mathcal{R}$ is defined by [41]

$${{}_0^{RL}D}_t^{-\alpha }f\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\tau )}^{\alpha -1}f\left(\tau \right)d\tau ,t>0,$$

(3)

or

$${{}_0^{RL}\mathbb{I}}_t^{\vartheta }f\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\tau )}^{\alpha -1}f\left(\tau \right)d\tau ,t>0,$$

(4)

where $\alpha >0$ and $\Gamma \left(.\right)$ is a Gamma function.

Definition 4.

The Riemann–Liouville fractional derivative of order α for a given $f:(0,\infty )\to \Re ,$ is stated as [41]

$${}_0^{RL}{D}_t^{\alpha }f\left(t\right)=\left\{\begin{array}{c}\frac{1}{\Gamma (n-\alpha )}{\left(\frac{d}{dt}\right)}^n{\int }_0^t\frac{f\left(\tau \right)}{{(t-\tau )}^{\alpha -n+1}}d\tau ,0\le n-1<\alpha <n,n=\left[\alpha \right],n\in N,\hfill \\ {\left(\frac{d}{dt}\right)}^{n}f\left(t\right),\alpha =n,n\in N.\hfill \end{array}\right.$$

(5)

Definition 5.

The Caputo fractional derivative of order α for a given $f:(0,\infty )\to \Re ,$ is stated as [41]

$${}_0^C{D}_t^{\alpha }f\left(t\right)=\left\{\begin{array}{c}\frac{1}{\Gamma (n-\alpha )}{\int }_0^t\frac{{\left(\frac{d}{d\tau }\right)}^{n}f\left(\tau \right)}{{(t-\tau )}^{\alpha -n+1}}d\tau ,0\le n-1<\alpha <n,n=\left[\alpha \right],\hfill \\ {\left(\frac{d}{dt}\right)}^{n}f\left(t\right),\alpha =n,n\in N.\hfill \end{array}\right.$$

(6)

Remark 1.

The Laplace transform (LT) of the Caputo operator of $f\left(t\right)$ order $\alpha >0$ is given as [63]

$$L\left[{}_0^C{D}_t^{\alpha }f\left(t\right)\right]=s^{\alpha }F\left(s\right)-\sum _{k=0}^{n-1}{f}^k\left(0\right)s^{\alpha -k-1}.$$

(7)

Theorem 1.

The Laplace transform of the function $t^{{\alpha }_1-1}{E}_{\alpha ,{\alpha }_1}(\pm \lambda t^{\alpha })$ is defined as [63]

$$L\left[t^{{\alpha }_1-1}{E}_{\alpha ,{\alpha }_1}(\pm \lambda t^{\alpha })\right]=\frac{s^{\alpha -{\alpha }_1}}{s^{\alpha }\mp \lambda },$$

(8)

where $E_{\alpha ,{\alpha }_1}$ is the two-parameter Mittag–Leffler function with $\alpha ,{\alpha }_1>0.$

3. Model Formulation

In this study, the host population has been divided into six epidemiological groups: the first compartment is susceptible to infection, the second one is latently infected, and the third compartment is individuals with acute infections. Individuals in this population are infected individuals who recover in a short time after exposure to Hepatitis-B disease or can infect susceptible individuals. Individuals in this group may also be exposed to death from Hepatitis-B. There are carrier individuals in the fourth compartment (carriers). In the fifth compartment (Recovered) are individuals who have recovered from Hepatitis-B. In the sixth compartment (vaccination), there are immunized individuals after vaccination. Very few of these individuals can be exposed to the disease again. The model is then presented by the following ODEs [66]:

$$\begin{array}{c}\hfill {\displaystyle \frac{dS}{dt}}=\lambda \rho \left(1-\kappa A\right)+\eta V-\left(\tau +\beta A+\omega \beta C+\gamma \right)S,\\ \hfill {\displaystyle \frac{dL}{dt}}=\left(\beta A+\omega \beta C\right)S-\left(\tau +\psi \right)L,\\ \hfill {\displaystyle \frac{dA}{dt}}=\psi L+(\lambda \rho \kappa -\tau -{\tau }_H-q-r)A,\\ \hfill {\displaystyle \frac{dC}{dt}}=qA-\phi C-\tau C,\\ \hfill {\displaystyle \frac{dR}{dt}}=rA-\tau R+\phi C,\\ \hfill {\displaystyle \frac{dV}{dt}}=\lambda \left(1-\rho \right)+\gamma S-\tau V-\eta V,\end{array}$$

(9)

where the parameters are given in Table 1. In this context, the total population can be considered as $N\left(t\right)=S\left(t\right)+L\left(t\right)+A\left(t\right)+C\left(t\right)+R\left(t\right)+V\left(t\right)$. We assume that in the Hepatitis-B acute group, the newborn virus population is less than the sum of the deaths of individuals with acute disease and the sum of the population from acute to immune state [62,66].

Otherwise, as long as there is an infection, individuals with acute disease will continue to increase rapidly. So, it becomes $dA/dt>0$ for $C\ne 0$, $A\ne 0$ and $t\le 0$.

Most of the studies modeling the biological systems are restricted to integer-order ordinary differential equations. Mathematical models, using ordinary differential equations with integer order, have been proved valuable in understanding the dynamics of diseases. However, models with fractional-order differential equations provide more advantages than integer-order mathematical models. The fractional-order models include a memory effect unlike the integer-order models. This special property will be convenient because of the fact that the essential features of the immune system also involve memory. Moreover, in the fractional systems, dimensional consistency is a very important tool, in which the units of measurement from the left- and right-hand sides of the equations are coherent. This consistency can be provided by modifying the parameters involved in the right-hand side of the equations, e.g., raising them to power $\alpha$. In this context, we extended Model (9) to the fractional order, which is presented in the following system:

$$\begin{array}{c}\hfill {}_0^C{D}_t^{\alpha }S\left(t\right)={\lambda }^{\alpha }{\rho }^{\alpha }\left(1-{\kappa }^{\alpha }A\right)+{\eta }^{\alpha }V-\left({\tau }^{\alpha }+{\beta }^{\alpha }A+{\omega }^{\alpha }{\beta }^{\alpha }C+{\gamma }^{\alpha }\right)S,\\ \hfill {}_0^C{D}_t^{\alpha }L\left(t\right)=\left({\beta }^{\alpha }A+{\omega }^{\alpha }{\beta }^{\alpha }C\right)S-\left({\tau }^{\alpha }+{\psi }^{\alpha }\right)L,\\ \hfill {}_0^C{D}_t^{\alpha }A\left(t\right)={\psi }^{\alpha }L+({\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }-{\tau }^{\alpha }-{\tau }_H^{\alpha })A-q^{\alpha }A-r^{\alpha }A,\\ \hfill {}_0^C{D}_t^{\alpha }C\left(t\right)=q^{\alpha }A-{\phi }^{\alpha }C-{\tau }^{\alpha }C,\\ \hfill {}_0^C{D}_t^{\alpha }R\left(t\right)=r^{\alpha }A-{\tau }^{\alpha }R-{\phi }^{\alpha }C,\\ \hfill {}_0^C{D}_t^{\alpha }V\left(t\right)={\lambda }^{\alpha }\left(1-{\rho }^{\alpha }\right)+{\gamma }^{\alpha }S-{\tau }^{\alpha }V-{\eta }^{\alpha }V,\end{array}$$

(10)

subject to the initial conditions,

$$S\left(0\right)=S_0,L\left(0\right)=L_0,A\left(0\right)=A_0,C\left(0\right)=C_0,R\left(0\right)={\mathcal{R}}_0,V\left(0\right)=V_0,$$

(11)

where $0<\alpha \le 1$, $(\left(S\left(t\right),L\left(t\right),A\left(t\right),C\left(t\right),R\left(t\right),V\left(t\right)\right)\in {\mathcal{R}}_{+}^6$ and if $\alpha =1$, then System (10) turns into an integer-order one (9). In the Equation (10) system, $S\left(t\right)$, $L\left(t\right)$, $A\left(t\right)$, $C\left(t\right)$, $R\left(t\right)$, $V\left(t\right)$ functions and Caputo fractional derivatives are assumed to be continuous when $t\ge 0$

$$\lambda \rho \kappa <\tau +{\tau }_H+r.$$

(12)

Table 1. Estimated values of the parameters.

Parameter	Parameter Meaning	Value	Reference	Unit
$\lambda$	Birth Rate	0.0121	[66]	year${}^{-1}$
$\tau$	Natural mortality rate	0.000034857	[8]	year${}^{-1}$
${\tau }_H$	Hepatitis-B related mortality rate	0.1019	Fitted	year${}^{-1}$
$\beta$	Transmission coefficient of the disease	0.00014334	Fitted	year${}^{-1}$
$\psi$	Transition rate from Latent population to Acute population	0.1989	Fitted	year${}^{-1}$
q	Transition rate of individuals with Acute infection to carrier-class	0.3387	Fitted	year${}^{-1}$
$\phi$	Recovery rate of individuals in the carrier class	0.0741	Fitted	year${}^{-1}$
$\gamma$	Vaccination rate	0.8569	Fitted	year${}^{-1}$
$\rho$	Rate of births without successful vaccination	0.00043102	Fitted	year${}^{-1}$
$\kappa$	Infected rate of mothers with HB Acute virus	0.0137	Fitted	year${}^{-1}$
$\eta$	The rate of decrease in immunity with the effect of vaccine	0.9472	Fitted	year${}^{-1}$
$\omega$	Reduced transmission rate compared to Acute	0.7534	Fitted	year${}^{-1}$
r	Recovery rate of individuals with Acute infection	0.0277	Fitted	year${}^{-1}$

Table 1. Estimated values of the parameters.

Parameter	Parameter Meaning	Value	Reference	Unit
$\lambda$	Birth Rate	0.0121	[66]	year${}^{-1}$
$\tau$	Natural mortality rate	0.000034857	[8]	year${}^{-1}$
${\tau }_H$	Hepatitis-B related mortality rate	0.1019	Fitted	year${}^{-1}$
$\beta$	Transmission coefficient of the disease	0.00014334	Fitted	year${}^{-1}$
$\psi$	Transition rate from Latent population to Acute population	0.1989	Fitted	year${}^{-1}$
q	Transition rate of individuals with Acute infection to carrier-class	0.3387	Fitted	year${}^{-1}$
$\phi$	Recovery rate of individuals in the carrier class	0.0741	Fitted	year${}^{-1}$
$\gamma$	Vaccination rate	0.8569	Fitted	year${}^{-1}$
$\rho$	Rate of births without successful vaccination	0.00043102	Fitted	year${}^{-1}$
$\kappa$	Infected rate of mothers with HB Acute virus	0.0137	Fitted	year${}^{-1}$
$\eta$	The rate of decrease in immunity with the effect of vaccine	0.9472	Fitted	year${}^{-1}$
$\omega$	Reduced transmission rate compared to Acute	0.7534	Fitted	year${}^{-1}$
r	Recovery rate of individuals with Acute infection	0.0277	Fitted	year${}^{-1}$

4. Analysis of the Model

4.1. Existence, Uniqueness, Positivity and Boundedness

Let ${\mathcal{R}}_{+}^6=\left\{\chi \left(t\right)\in {\mathcal{R}}^6\right.:\chi \left(t\right)\ge \left.0\right\}$ and $\chi \left(t\right)={\left[S\left(t\right),L\left(t\right),A\left(t\right),C\left(t\right),R\left(t\right),V\left(t\right)\right]}^T$. We will review the following lemma (generalized mean value theorem [60]) in order to demonstrate the proof of the basic theorem regarding the non-negativity of the solutions for Model (10).

Lemma 1.

Suppose that $f\left(t\right)\in C[a,b]$ and Caputo fractional derivative ${{}_0^{C}D}_t^{\alpha }f\left(t\right)\in C(a,b]$ for $0<\alpha \le 1,$ then

$$f\left(t\right)=f\left(\omega \right)+\frac{1}{\Gamma \left(\alpha \right)}{{}_0^{C}D}_t^{\alpha }f\left(\tau \right){(t-\omega )}^{\alpha },$$

with $0\le \tau \le t,\forall t\in (a,b].$

Remark 2.

If $f\left(t\right)\in C[0,b]$ and ${{}_0^{C}D}_t^{\alpha }f\left(t\right)\in (0,b]$ for $0<\alpha \le 1$. It is clear from Lemma 1 that if ${{}_0^{C}D}_t^{\alpha }f\left(t\right)\ge 0,\forall t\in (0,b],$ then $f\left(t\right)$ is non-decreasing and if ${{}_0^{C}D}_t^{\alpha }f\left(t\right)\le 0,\forall t\in (0,b],$ then $f\left(t\right)$ is non-increasing for all $t\in [0,b]$ [61].

Theorem 2.

The solution of Model (10) with the initial Conditions (11) is unique and in ${\mathcal{R}}_{+}^6.$

Proof. 

Lemma 1 and Remark 2 can be used to demonstrate the existence and uniqueness of the solution of Systems (10) and (11) on the time interval $\left(0,\infty \right)$. Following the method of [60], we now explain that the non-negative region ${\mathcal{R}}_{+}^6,$ is a positively invariant zone. From Model (10), we find

$$\begin{array}{c}\hfill {}_0^C{D}_t^{\alpha }S\left(t\right){|}_{S=0}={\lambda }^{\alpha }{\rho }^{\alpha }\left(1-{\kappa }^{\alpha }A\right)+{\eta }^{\alpha }V\ge 0,\\ \hfill {}_0^C{D}_t^{\alpha }L\left(t\right){|}_{L=0}=\left({\beta }^{\alpha }A+{\omega }^{\alpha }{\beta }^{\alpha }C\right)S\ge 0,\\ \hfill {}_0^C{D}_t^{\alpha }A\left(t\right){|}_{A=0}={\psi }^{\alpha }L\ge 0,\\ \hfill {}_0^C{D}_t^{\alpha }C\left(t\right){|}_{C=0}=q^{\alpha }A\ge 0,\\ \hfill {}_0^C{D}_t^{\alpha }R\left(t\right){|}_{R=0}=r^{\alpha }A+{\phi }^{\alpha }C\ge 0,\\ \hfill {}_0^C{D}_t^{\alpha }V\left(t\right){|}_{V=0}={\lambda }^{\alpha }-{\lambda }^{\alpha }{\rho }^{\alpha }+{\gamma }^{\alpha }S\ge 0.\end{array}$$

(13)

For further reading, one can see [62]. Since $(1-{\rho }^{\alpha })$ indicates the successful vaccine birth rate, the expression ${\lambda }^{\alpha }-{\lambda }^{\alpha }{\rho }^{\alpha }$ is always positive or zero. □

If $\left(S\left(0\right),L\left(0\right),A\left(0\right),C\left(0\right),R\left(0\right),V\left(0\right)\right)\in {\mathcal{R}}_{+}^6,$ then according to Equation (13) and Remark 2, the solution $\left(S\left(t\right),L\left(t\right),A\left(t\right),C\left(t\right),R\left(t\right),V\left(t\right)\right)$ cannot escape from the hyperplanes $S\ge 0,L\ge 0,\ge 0,A\ge 0,C\ge 0,R\ge 0,V\ge 0.$ Additionally, on each hyperplane enclosing the non-negative orthant, the vector field points into ${\mathcal{R}}_{+}^6$, meaning that the domain ${\mathcal{R}}_{+}^6$ is a positively invariant set.

Theorem 3.

Taking into account the positive invariant set ${\mathcal{R}}_{+}^6$, the solution of System (10)

$$\begin{array}{c}\hfill \mathcal{A}=\left\{\right(S\left(t\right),L\left(t\right),A\left(t\right),C\left(t\right),R\left(t\right),V\left(t\right))\in {\mathcal{R}}_{+}^6|0<S\left(t\right)+L\left(t\right)+A\left(t\right)\\ \hfill +C\left(t\right)+R\left(t\right)+V\left(t\right)\le \frac{{\lambda }^{\alpha }}{{\tau }^{\alpha }}\}\end{array}$$

is a positive invariant set for System (10).

Proof. 

We then search for situations where trajectories are confined to a compact set and do not “escape to infinity”. This boundedness, which in biology denotes that no population may develop indefinitely, allows the model to accurately represent the dynamics. Take a look at the new function N = S + L + A + C + R + V, which stands for the total population of the 6D system. Therefore, the Caputo fractional derivative of the total population (N) is as follows:

$${}_0^C{D}_t^{\alpha }N\left(t\right)={\lambda }^{\alpha }-{\tau }^{\alpha }N-{\tau }_H^{\alpha }A.$$

From here, we have a total population (N) as follows:

$${}_0^C{D}_t^{\alpha }N\left(t\right)\le {\lambda }^{\alpha }-{\tau }^{\alpha }N.$$

(14)

Applying the Laplace transform to Equation (14), we obtain

$$s^{\alpha }\tilde{N}\left(s\right)-s^{\alpha -1}N\left(0\right)\le \frac{{\lambda }^{\alpha }}{s}-{\tau }^{\alpha }\tilde{N}\left(s\right).$$

It also has the following form:

$$\tilde{N}\left(s\right)\le \frac{s^{-1}}{s^{\alpha }+{\tau }^{\alpha }}{\lambda }^{\alpha }+\frac{s^{\alpha -1}}{s^{\alpha }+{\tau }^{\alpha }}N\left(0\right).$$

In the last equation, if we take the inverse Laplace transform, then we have the following:

$$\begin{array}{c}\hfill N\left(t\right)\le {\lambda }^{\alpha }{t}^{\alpha }{E}_{\alpha ,\alpha +1}(-{\tau }^{\alpha }{t}^{\alpha })+E_{\alpha ,1}(-{\tau }^{\alpha }{t}^{\alpha })N\left(0\right)\\ \hfill \le \frac{{\lambda }^{\alpha }}{{\tau }^{\alpha }}\left({\tau }^{\alpha }{t}^{\alpha }{E}_{\alpha ,\alpha +1}(-{\tau }^{\alpha }{t}^{\alpha })\right)+E_{\alpha ,1}(-{\tau }^{\alpha }{t}^{\alpha })\le \frac{{\lambda }^{\alpha }}{{\tau }^{\alpha }}\frac{1}{\Gamma \left(1\right)}\le \frac{{\lambda }^{\alpha }}{{\tau }^{\alpha }}.\end{array}$$

□

4.2. Equilibria and Stability

By equating the right side of System (10) to zero, the equilibrium points are obtained:

$$\begin{array}{c}\hfill {\lambda }^{\alpha }{\rho }^{\alpha }\left(1-{\kappa }^{\alpha }A\right)+{\eta }^{\alpha }V-\left({\tau }^{\alpha }+{\beta }^{\alpha }A+{\omega }^{\alpha }{\beta }^{\alpha }C+{\gamma }^{\alpha }\right)S=0,\\ \hfill \left({\beta }^{\alpha }A+{\omega }^{\alpha }{\beta }^{\alpha }C\right)S-\left({\tau }^{\alpha }+{\psi }^{\alpha }\right)L=0,\\ \hfill {\psi }^{\alpha }L+{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }A-{\tau }^{\alpha }A-{\tau }_H^{\alpha }A-q^{\alpha }A-r^{\alpha }A=0,\\ \hfill q^{\alpha }A-{\phi }^{\alpha }C-{\tau }^{\alpha }C=0,\\ \hfill r^{\alpha }A-{\tau }^{\alpha }R+{\phi }^{\alpha }C=0,\\ \hfill {\lambda }^{\alpha }\left(1-{\rho }^{\alpha }\right)+{\gamma }^{\alpha }S-{\tau }^{\alpha }V-{\eta }^{\alpha }V=0.\end{array}$$

(15)

After simplification, the disease-free equilibrium (DFE), namely $E_0=(S_0,0,0,0,0,V_0)$, where

$$\begin{array}{c}\hfill S_0={\displaystyle \frac{{\lambda }^{\alpha }({\rho }^{\alpha }{\tau }^{\alpha }+{\eta }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\gamma }^{\alpha }+{\eta }^{\alpha })}},V_0={\displaystyle \frac{{\lambda }^{\alpha }({\tau }^{\alpha }+{\gamma }^{\alpha }-{\tau }^{\alpha }{\rho }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\gamma }^{\alpha }+{\eta }^{\alpha })}},\end{array}$$

(16)

The endemic equilibrium (EE) $E^{*}=(S^{*},L^{*},A^{*},C^{*},R^{*},V^{*})$, where

$$\begin{array}{c}\hfill S^{*}={\displaystyle \frac{(-{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+{\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha })({\tau }^{\alpha }+{\psi }^{\alpha })({\phi }^{\alpha }+{\tau }^{\alpha })}{{\psi }^{\alpha }+{\beta }^{\alpha }({\phi }^{\alpha }+{\tau }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })}}={\displaystyle \frac{S_0}{{\mathcal{R}}_0}},\\ \hfill L^{*}={\displaystyle \frac{(-{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+{\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha })({\phi }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\eta }^{\alpha })S_0({\lambda }^{\alpha }{\rho }^{\alpha }-S^{*}({\tau }^{\alpha }+{\gamma }^{\alpha }))}{{\psi }^{\alpha }({\tau }^{\alpha }+{\eta }^{\alpha })S_0(S^{*}{\beta }^{\alpha }(({\phi }^{\alpha }+{\tau }^{\alpha })+{\omega }^{\alpha }{q}^{\alpha })+({\phi }^{\alpha }+{\tau }^{\alpha }){\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha })}}\\ \hfill +{\displaystyle \frac{{\eta }^{\alpha }({\phi }^{\alpha }+{\tau }^{\alpha })({\lambda }^{\alpha }-S_0({\lambda }^{\alpha }{\rho }^{\alpha }+{\gamma }^{\alpha }))}{{\psi }^{\alpha }({\tau }^{\alpha }+{\eta }^{\alpha })S_0(S^{*}{\beta }^{\alpha }(({\phi }^{\alpha }+{\tau }^{\alpha })+{\omega }^{\alpha }{q}^{\alpha })+({\phi }^{\alpha }+{\tau }^{\alpha }){\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha })}},\\ \hfill A^{*}={\displaystyle \frac{({\phi }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\eta }^{\alpha })S_0({\lambda }^{\alpha }{\rho }^{\alpha }-S^{*}({\tau }^{\alpha }+{\gamma }^{\alpha }))+{\eta }^{\alpha }({\phi }^{\alpha }+{\tau }^{\alpha })({\lambda }^{\alpha }-S_0({\lambda }^{\alpha }{\rho }^{\alpha }+{\gamma }^{\alpha }))}{({\tau }^{\alpha }+{\eta }^{\alpha })S_0(S^{*}{\beta }^{\alpha }(({\phi }^{\alpha }+{\tau }^{\alpha })+{\omega }^{\alpha }{q}^{\alpha })+({\phi }^{\alpha }+{\tau }^{\alpha }){\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha })}},\\ \hfill C^{*}={\displaystyle \frac{q^{\alpha }({\tau }^{\alpha }+{\eta }^{\alpha })S_0({\lambda }^{\alpha }{\rho }^{\alpha }-S^{*}({\tau }^{\alpha }+{\gamma }^{\alpha }))+{\eta }^{\alpha }({\phi }^{\alpha }+{\tau }^{\alpha })({\lambda }^{\alpha }{\rho }^{\alpha }+{\gamma }^{\alpha }))}{({\tau }^{\alpha }+{\eta }^{\alpha })S_0(S^{*}{\beta }^{\alpha }({\phi }^{\alpha }+{\tau }^{\alpha })+{\omega }^{\alpha }{q}^{\alpha })+({\phi }^{\alpha }+{\tau }^{\alpha }){\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha })}},\\ \hfill R^{*}={\displaystyle \frac{r({\phi }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\eta }^{\alpha })S_0({\lambda }^{\alpha }{\rho }^{\alpha }-S^{*}({\tau }^{\alpha }+{\gamma }^{\alpha }))+{\eta }^{\alpha }({\phi }^{\alpha }+{\tau }^{\alpha })({\lambda }^{\alpha }-S_0({\lambda }^{\alpha }{\rho }^{\alpha }+{\gamma }^{\alpha }))}{({\tau }^{\alpha }+{\eta }^{\alpha })S_0(S^{*}{\beta }^{\alpha }(({\phi }^{\alpha }+{\tau }^{\alpha })+{\omega }^{\alpha }{q}^{\alpha })+({\phi }^{\alpha }+{\tau }^{\alpha }){\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha })}}\\ \hfill +{\displaystyle \frac{{\phi }^{\alpha }{q}^{\alpha }({\tau }^{\alpha }+{\eta }^{\alpha })S_0({\lambda }^{\alpha }{\rho }^{\alpha }-S^{*}({\tau }^{\alpha }+{\gamma }^{\alpha }))+{\eta }^{\alpha }({\phi }^{\alpha }+{\tau }^{\alpha })({\lambda }^{\alpha }{\rho }^{\alpha }+{\gamma }^{\alpha }))}{{\tau }^{\alpha }({\tau }^{\alpha }+{\eta }^{\alpha })S_0(S^{*}{\beta }^{\alpha }({\phi }^{\alpha }+{\tau }^{\alpha })+{\omega }^{\alpha }{q}^{\alpha })+({\phi }^{\alpha }+{\tau }^{\alpha }){\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha })}},\\ \hfill V^{*}={\displaystyle \frac{{\lambda }^{\alpha }-{\lambda }^{\alpha }{\rho }^{\alpha }+{\gamma }^{\alpha }{S}^{*}}{{\tau }^{\alpha }+{\eta }^{\alpha }}}.\end{array}$$

(17)

4.3. Basic Reproduction Number

For the local stability of the disease-free equilibrium, we first compute the basic reproduction number by using the next-generation matrix method [65]. It is defined as the number of cases occurring in a population that is fully susceptible to any infectious individual. In biological models, if ${\mathcal{R}}_0<1$, infection disappears, if ${\mathcal{R}}_0>1$, there is infection and the disease continues. To determine ${\mathcal{R}}_0$, which is considered the spectral radius of the next-generation matrix $F{V}^{-1}$, we assemble the compartments which are infected from System (10) and decompose the right-hand side as $\mathcal{F}-\mathcal{V}$, where $\mathcal{F}$ is the transmission part, expressing the production of a new infection, and $\mathcal{V}$ is the transition part which describes the change in the state. Therefore,

$$\mathcal{F}\left(\chi \right)=\left(\begin{array}{c}0\\ ({\beta }^{\alpha }A+{\omega }^{\alpha }{\beta }^{\alpha }C)S\\ 0\\ 0\\ 0\\ 0\end{array}\right),\mathcal{V}\left(\chi \right)=\left(\begin{array}{c}-{\lambda }^{\alpha }{\rho }^{\alpha }\left(1-{\kappa }^{\alpha }A\right)-{\eta }^{\alpha }V+P_{1}S\\ \left({\tau }^{\alpha }+{\psi }^{\alpha }\right)L\\ -{\psi }^{\alpha }L-{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }A+P_{2}A\\ -q^{\alpha }A+{\phi }^{\alpha }C+{\tau }^{\alpha }C\\ {\tau }^{\alpha }R-r^{\alpha }A-{\phi }^{\alpha }C\\ V({\eta }^{\alpha }+{\tau }^{\alpha })-{\lambda }^{\alpha }(1-{\rho }^{\alpha })-{\gamma }^{\alpha }S\end{array}\right),$$

where $P_1={\tau }^{\alpha }+{\beta }^{\alpha }A+{\omega }^{\alpha }{\beta }^{\alpha }C+{\gamma }^{\alpha }$ and $P_2={\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }.$ By the next-generation matrix method [65], the matrices F and V at the disease-free equilibrium point $E_0$ are obtained by $F=\left[\frac{\partial {\mathcal{F}}_i\left(E_0\right)}{\partial {\chi }_j}\right]$ and $V=\left[\frac{\partial {\mathcal{V}}_i\left(E_0\right)}{\partial {\chi }_j}\right]$, $1\le i,j\le 3$, where $\chi$ represents the number of individuals in each compartment. This implies

$$\mathbf{F}=\left(\begin{array}{ccc}0& {\beta }^{\alpha }{S}_0& {\omega }^{\alpha }{\beta }^{\alpha }{S}_0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),\mathbf{V}=\left(\begin{array}{ccc}{\tau }^{\alpha }+{\psi }^{\alpha }& 0& 0\\ -{\psi }^{\alpha }& {\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }& 0\\ 0& -q^{\alpha }& {\phi }^{\alpha }+{\tau }^{\alpha }\end{array}\right).$$

The expression for ${\mathcal{R}}_0$ is the spectral radius of the matrix $F{V}^{-1}$ and is written as follows:

$$\begin{array}{c}\hfill {\mathcal{R}}_0=\frac{{\lambda }^{\alpha }{\beta }^{\alpha }{\psi }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}.\end{array}$$

(18)

4.4. Stability of Equilibria

We present the local stability findings of equilibrium points as theorems with justifications in this subsection.

Theorem 4.

The suggested fractional-order Hepatitis-B epidemic model’s disease-free equilibrium $E_0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Proof. 

The Jacobian matrix at $E_0$ is

$$\mathcal{J}\left(E_0\right)=\left(\begin{array}{cccccc}-({\tau }^{\alpha }+{\gamma }^{\alpha })& 0& -{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }& -{\omega }^{\alpha }{\beta }^{\alpha }{S}_0& 0& {\eta }^{\alpha }\\ 0& -({\tau }^{\alpha }+{\psi }^{\alpha })& {\beta }^{\alpha }{S}_0& {\omega }^{\alpha }{\beta }^{\alpha }{S}_0& 0& 0\\ 0& {\psi }^{\alpha }& {\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }-({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha })& 0& 0& 0\\ 0& 0& q^{\alpha }& -({\phi }^{\alpha }+{\tau }^{\alpha })& 0& 0\\ 0& 0& r^{\alpha }& {\phi }^{\alpha }& -{\tau }^{\alpha }& 0\\ {\gamma }^{\alpha }& 0& 0& 0& 0& -({\tau }^{\alpha }+{\eta }^{\alpha })\end{array}\right).$$

The characteristic equation is

$$P\left({\lambda }^{\alpha }\right)={({\lambda }^{\alpha }+{\tau }^{\alpha })}^2({\lambda }^{\alpha }+{\gamma }^{\alpha }+{\eta }^{\alpha }+{\tau }^{\alpha })({\left({\lambda }^{\alpha }\right)}^3+A{\left({\lambda }^{\alpha }\right)}^2+B{\lambda }^{\alpha }+C)=0,$$

where

$$\begin{array}{c}\hfill A=2{\tau }^{\alpha }+{\phi }^{\alpha }+{\psi }^{\alpha }-{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+P_2,\\ \hfill B={\left({\tau }^{\alpha }\right)}^2-({\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+P_2)(2{\tau }^{\alpha }{\phi }^{\alpha }+{\psi }^{\alpha })+{\tau }^{\alpha }{\phi }^{\alpha }+{\tau }^{\alpha }{\psi }^{\alpha }+{\phi }^{\alpha }{\psi }^{\alpha }-\frac{{\beta }^{\alpha }{\lambda }^{\alpha }({\tau }^{\alpha }{\rho }^{\alpha }+{\eta }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\gamma }^{\alpha }+{\eta }^{\alpha })}{\psi }^{\alpha },\\ \hfill C=-\frac{{\beta }^{\alpha }{\lambda }^{\alpha }({\tau }^{\alpha }{\rho }^{\alpha }+{\eta }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\gamma }^{\alpha }+{\eta }^{\alpha })}{\eta }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })-({\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+P_2)({\left({\tau }^{\alpha }\right)}^2+{\tau }^{\alpha }+{\phi }^{\alpha }+{\tau }^{\alpha }{\psi }^{\alpha }+{\phi }^{\alpha }{\psi }^{\alpha }),\end{array}$$

in which

$$\begin{array}{c}\hfill A=2{\tau }^{\alpha }+{\phi }^{\alpha }+{\psi }^{\alpha }-{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+{\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }>0,\\ \hfill C=(1-{\mathcal{R}}_0)({\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }-({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }))({\left({\tau }^{\alpha }\right)}^2+{\tau }^{\alpha }{\phi }^{\alpha }+{\tau }^{\alpha }{\psi }^{\alpha }+{\phi }^{\alpha }{\psi }^{\alpha })>0,\\ \hfill AB-C=({\left({\tau }^{\alpha }\right)}^2-({\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+P_2)(2{\tau }^{\alpha }{\phi }^{\alpha }+{\psi }^{\alpha })+{\tau }^{\alpha }{\phi }^{\alpha }+{\tau }^{\alpha }{\psi }^{\alpha }+{\phi }^{\alpha }{\psi }^{\alpha }-\frac{{\beta }^{\alpha }{\lambda }^{\alpha }({\tau }^{\alpha }{\rho }^{\alpha }+{\eta }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\gamma }^{\alpha }+{\eta }^{\alpha })}{\psi }^{\alpha })\\ \hfill \times (2{\tau }^{\alpha }+{\phi }^{\alpha }+{\psi }^{\alpha }-{\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+P_2)+(\frac{{\beta }^{\alpha }{\lambda }^{\alpha }({\tau }^{\alpha }{\rho }^{\alpha }+{\psi }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\gamma }^{\alpha }+{\eta }^{\alpha })}{\eta }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })\\ \hfill +({\lambda }^{\alpha }{\rho }^{\alpha }{\kappa }^{\alpha }+P_2)({\left({\tau }^{\alpha }\right)}^2+{\tau }^{\alpha }+{\phi }^{\alpha }+{\tau }^{\alpha }{\psi }^{\alpha }+{\phi }^{\alpha }{\psi }^{\alpha }))>0.\end{array}$$

Therefore, by Routh–Hurwitz criteria, all roots of $P\left({\lambda }^{\alpha }\right)$ have negative real parts, and $E_0$ is stable. Furthermore, if ${\mathcal{R}}_0>1$, we have $A<0$ and $E_0$ is unstable. □

Remark 3.

The endemic equilibrium $E^{*}$ of the suggested Hepatitis-B model of fractional order is locally asymptotically stable if ${\mathcal{R}}_0>1$ and is unstable if ${\mathcal{R}}_0<1$.

5. Sensitivity Analysis

In this section, we examine the sensitivity of ${\mathcal{R}}_0$ according to parameters that impact the reproduction number. In these analyses, only the effects of ${\beta }^{\alpha }$, ${\lambda }^{\alpha }$, ${\omega }^{\alpha }$, ${\kappa }^{\alpha }$, ${\psi }^{\alpha }$, ${\phi }^{\alpha }$, ${\rho }^{\alpha }$, ${\gamma }^{\alpha }$, ${\eta }^{\alpha }$ and $q^{\alpha }$, $r^{\alpha }$ parameters on the relevant reproduction number values are taken into account. We apply the same method as in [67] and obtain the following:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\beta }^{\alpha }}}=\frac{{\lambda }^{\alpha }{\psi }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }q)({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}>0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\lambda }^{\alpha }}}=\frac{({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha }){\beta }^{\alpha }{\psi }^{\alpha }}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}\\ \hfill +\frac{{\lambda }^{\alpha }{\beta }^{\alpha }{\psi }^{\alpha }{\kappa }^{\alpha }{\rho }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })}{({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}>0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\omega }^{\alpha }}}=\frac{{\lambda }^{\alpha }{q}^{\alpha }{\beta }^{\alpha }{\psi }^{\alpha }({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\eta }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}>0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\kappa }^{\alpha }}}=\frac{{{\lambda }^{\alpha }}^2{\psi }^{\alpha }{\rho }^{\alpha }{\beta }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha }){({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}^2}>0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\psi }^{\alpha }}}=\frac{{\lambda }^{\alpha }{\beta }^{\alpha }{\psi }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })\left(1-({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })\right)}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha }){({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })}^2({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}>0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\phi }^{\alpha }}}=-\frac{{\lambda }^{\alpha }{\beta }^{\alpha }{\psi }^{\alpha }{\omega }^{\alpha }{q}^{\alpha }({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })}{{\left(({\tau }^{\alpha }+{\phi }^{\alpha })\right)}^2({\tau }^{\alpha }+{\psi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}<0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\eta }^{\alpha }}}=\frac{{\lambda }^{\alpha }{\beta }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })}{({\tau }^{\alpha }+{\phi }^{\alpha }){\left(({\tau }^{\alpha }+{\psi }^{\alpha })\right)}^2({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}>0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\rho }^{\alpha }}}=\frac{{\lambda }^{\alpha }{\beta }^{\alpha }{\eta }^{\alpha }{\tau }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}>0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial {\gamma }^{\alpha }}}=-\frac{{\lambda }^{\alpha }{\beta }^{\alpha }{\psi }^{\alpha }({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha }){({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })}^2({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}<0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial q^{\alpha }}}=\frac{{\lambda }^{\alpha }{\beta }^{\alpha }{\psi }^{\alpha }({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })({\omega }^{\alpha }-{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha })({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}>0,\\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial r^{\alpha }}}=-\frac{{\lambda }^{\alpha }{\beta }^{\alpha }{\psi }^{\alpha }({\eta }^{\alpha }+{\rho }^{\alpha }{\tau }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha }+{\omega }^{\alpha }{q}^{\alpha })}{{\tau }^{\alpha }({\tau }^{\alpha }+{\psi }^{\alpha })({\tau }^{\alpha }+{\phi }^{\alpha })({\eta }^{\alpha }+{\gamma }^{\alpha }+{\tau }^{\alpha }){({\tau }^{\alpha }+{\tau }_H^{\alpha }+q^{\alpha }+r^{\alpha }-{\kappa }^{\alpha }{\lambda }^{\alpha }{\rho }^{\alpha })}^2}<0.\end{array}$$

(19)

In the sensitivity analysis, one can obtain that the values of ${\mathcal{R}}_0$ increase and decrease in proportion to the growth of ${\beta }^{\alpha }$, ${\lambda }^{\alpha }$, ${\omega }^{\alpha }$, ${\kappa }^{\alpha }$, ${\psi }^{\alpha }$, ${\eta }^{\alpha }$, ${\rho }^{\alpha }$, $q^{\alpha }$ and ${\phi }^{\alpha }$, ${\gamma }^{\alpha }$, $r^{\alpha }$ values, respectively. In addition, Figure 1a shows the changes of the reproduction number according to the rate of decrease in immunity with the vaccination strategy (${\eta }^{\alpha }$) versus the transition rate of individuals with acute infection to carrier class ($q^{\alpha }$). Moreover, Figure 1b depicts the changes of the reproduction number according to the vaccination rate (${\gamma }^{\alpha }$) versus the infected rate of mothers with the HB acute virus (${\kappa }^{\alpha }$).

The variables are taken into account in the analysis of the pertinent reproduction number within the established bounds. After considering the analyses and visuals, it is determined that reasonable actions should be taken to stop the disease from spreading by reducing the conditions that make the derivative of reproduction numbers positive and increasing the ones that make it negative.

6. Parameter Estimation

An essential component of an epidemiological model’s validation is the fitting of its parameters. This improves the model’s accuracy to better understand the transmissions of the epidemic and predict the future directions of diseases. Therefore, in this section, we aim to explain the determining of parameters with the least squares curve-fitting technique. A total of 13 different parameters are stated in the proposed model for Hepatitis-B infectious disease. We have taken two of these parameters ($\tau =1/(78.6\times 365)$ and $\lambda =0.0121$) from the literature, and the rest of them ensure the best estimation based on the real Hepatitis-B cases in Türkiye. The initial conditions are determined by dividing the individuals in each compartment by $N\left(t\right)$ (to the total population). In this context, the initial conditions are $S\left(0\right)=0.8565$, $L\left(0\right)=0.0122363$, $A\left(0\right)=0.1067495$, $C\left(0\right)=0.0048945$, $R\left(0\right)=0.0073417$ and $V\left(0\right)=0.0122363$. With the help of the least square curve-fitting method, there are 13 biological parameters predicted, and as shown in Figure 2, they ensure that the solution of the proposed HB model best fits the real pandemic cases. The ideal values of the pertinent parameters are reached to minimize the mean absolute relative error between Hepatitis-B cases and the model’s solution for the infectious class. Real HB cases are shown with solid red circles, while the model’s best-fit curve is shown with a blue line. The biological parameters that are taken into account for the model are presented in Table 1 along with the best estimates that could be made using the least squares method.

7. Numerical Technique

To approximate the Caputo fractional derivative of order $\alpha$, the fractional variant of the Adams–Bashforth approach is discussed in this section. The two-step Lagrange polynomial and the fundamental theorem of fractional calculus are combined to create this numerical approach [68]. Take into account the Caputo fractional-order system as

$${}_0^C{D}_t^{\alpha }f\left(t\right)=G(t,f\left(t\right)),f\left(0\right)=f_0,0<t<T<\infty ,$$

(20)

By using the fundamental theorem of fractional calculus, Equation (20) may be transformed into an integral type to obtain

$$f\left(t\right)-f\left(0\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t}g\left(\tau ,f\left(\tau \right)\right){\left(t-\tau \right)}^{\alpha -1}d\tau .$$

(21)

At $t=t_{n+1},n=0,1,2,\dots ,$ we have

$$f\left(t_{n+1}\right)-f\left(0\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_{n+1}}g(\tau ,f\left(\tau \right)){(t_{n+1}-\tau )}^{\alpha -1}d\tau .$$

(22)

Similarly, at $t=t_n,n=0,1,2,\dots ,$ we obtain

$$f\left(t_n\right)-f\left(0\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}g(\tau ,f\left(\tau \right)){(t_n-\tau )}^{\alpha -1}d\tau .$$

(23)

Further simplification is achieved by applying the Lagrange interpolation on $g(\tau ,f(\tau \left)\right)$ and removing Equation (23) from Equation (22). This results in

$$f\left(t_{n+1}\right)=f\left(t_n\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_{n+1}}g(\tau ,f\left(\tau \right)){(t_{n+1}-\tau )}^{\alpha -1}d\tau -\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^{t_n}g(\tau ,f\left(\tau \right)){(t_n-\tau )}^{\alpha -1}d\tau .$$

(24)

As a result, the final approximation scheme is as follows:

$$\begin{array}{c}\hfill f\left(t_{n+1}\right)=f\left(t_n\right)+\frac{g(t_n,f_n)}{h\Gamma \left(\alpha \right)}\left\{\frac{2h{t}_{n+1}^{\alpha }}{\alpha }-\frac{t_1^{\alpha +1}n+1}{\alpha +1}+\frac{h{t}_n^{\alpha }}{\alpha }-\frac{t_n^{\alpha +1}}{\alpha }\right\}\\ \hfill +\frac{g(t_{n-1},f_{n-1})}{h\Gamma \left(\alpha \right)}\left\{\frac{h{t}_{n+1}^{\alpha }}{\alpha }-\frac{t_1^{\alpha +1}n+1}{\alpha +1}+\frac{t_n^{\alpha }}{\alpha +1}\right\}.\end{array}$$

(25)

For the numerical approximation of the Caputo derivative, this method is referred to as the two-step fractional Adams–Bashforth method [68,69,70,71,72]. The suggested model Equation (10) can be solved using the fractional Adams–Bashforth approach as explained above.

8. Numerical Analysis and Results

The results of such biological studies should be presented numerically. We will discuss the biological significance of the stable equilibrium points of the system given by Equation (10) and examine the behavior of the results obtained in this section. The numerical parameter values of the $\alpha =(0,1]$ fractional-order system given by Equation (10) are given in Table 1.

Additionally, the initial value conditions of System (10) are S(0) = 0.8565, L(0) = 0.01223, A(0) = 0.10675, C(0) = 0.00489, R(0) = 0.00734 and V(0) = 0.01224. In Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8, the behavior of System (10) for order values $\alpha =0.85$, $\alpha =0.90$, $\alpha =0.95$ and $\alpha =1$ is obtained.

In Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 for fractional order $\alpha$ values, the time-dependent variation of the population of the disease compartments in susceptible, latent, acute, carrier, recovered and vaccinated individuals, respectively, is discussed, and the future course is tried to be estimated. In Figure 3, it is seen that, as the order of $\alpha$ decreases in the long run, there are more susceptible individuals in the case of fractional derivatives than in the case of integer derivatives. In Figure 4 and Figure 5, it is observed that the number of latent individuals and acute individuals decreases over time, and there are more individuals in the case of fractional derivatives in the long run. In addition, it is seen that fractional cases and integer-order cases approach each other and become stable over time. According to Figure 6, the peak seen in carrier individuals occurred between about 4 and 6 years for the $\alpha =1$ value, while it occurred in earlier years for the fractional values of $\alpha$. For example, for $\alpha =0.85$, carrier individuals seem to have peaked before the 4th year. It can be said that this is because fractional derivatives have a memory effect. In Figure 7 and Figure 8, changes over time in recovered and vaccinated individuals are examined for the fractional case and the integer-order case, respectively. Up to the 10th year, it is observed that the number of recovered individuals decreases as the $\alpha$ level increases, whereas in the long term, as the $\alpha$ increases, the number of recovered individuals increases. In Figure 8, there is a noticeable increase in the number of vaccinated individuals in the long term. The estimated number of vaccinated individuals for the fractional case appears to be higher than the number of vaccinated individuals for the integer-order case.

In Figure 9, the variation of the population of susceptible individuals over time have been investigated for the values of $\beta$ = 0.00014334, $\beta$ = 0.3, $\beta$ = 0.4 and $\beta$ = 0.5. Here, the most appropriate $\beta$ value for the real data is calculated as $\beta$ = 0.00014334 by the parameter estimation method. According to this value, it is observed that the number of susceptible individuals increased over time compared to other values of $\beta$ in the simulation. Similarly, in Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14, for different $\beta$ =0.00014334, $\beta$ = 0.3, $\beta$ = 0.4 and $\beta$ = 0.5 values of $\beta$, which is the transmission coefficient of the disease, the time-dependent variation of the population in susceptible, latent, acute, carrier, recovered and vaccinated individuals have been examined, respectively. The estimated number of latent individuals for the $\beta =0.00014334$ value obtained by the parameter estimation method is less than the estimated latent individuals for the $\beta$ = 0.3, $\beta$ =0.4, $\beta$ = 0.5, values. Additionally, for the value $\beta$ = 0.00014334, the number of latent individuals appears to be more stable than the others. It is observed that the number of acute individuals is also significantly less for the value of $\beta =0.00014334$ compared to other values and approaches to zero at the end of the simulation period.

In Figure 12 and Figure 13, it is seen that the estimated number of carrier and recovered individuals for $\beta =0.00014334$ is less than the estimated number of carrier and recovered individuals for other $\beta$ values. According to Figure 14, it is seen that the number of vaccinated individuals increases as the $\beta$ value increases.

In Figure 15, for different values of $\gamma$ ($\gamma$ = 0.8569, $\gamma$ = 0.6, $\gamma$ = 0.5 and $\gamma$ = 0.4) which is the vaccination rate in the (10) system, the variation over time in susceptible individuals is considered. Here, it is observed that the number of susceptible individuals in the population for different $\gamma$ values at the beginning is close to each other, and as the $\gamma$ value increases in the following years, the number of susceptible individuals in the population decreases. This situation can be interpreted as individuals who have never been infected with the disease (susceptible) getting full protection against the HB virus (effective vaccine) as a result of vaccination. Similarly, the variation of latent individuals in the population over time is investigated for different $\gamma$ values in Figure 16. According to Figure 16, it is observed that $\gamma$ has no significant effect on the latent population. This is interpreted as a result of the vaccination not being applied to latent individuals. In Figure 17, the time-dependent variation of vaccinated individuals in the population for different $\gamma$ values is examined. Here, it is seen that as the value of $\gamma$ increases, the number of vaccinated individuals in the population increases. This can be interpreted as the vaccine being effective against the HB virus and providing full protection.

In Figure 18, the time variation of susceptible individuals for different $\psi$ = 0.1989, $\psi$ = 0.3, $\psi$ = 0.4 and $\psi$ = 0.5 values of $\psi$ is examined. Here, $\psi$ = 0.1989 is the value obtained after applying the parameter estimation method. It is seen that as the value of $\psi$ increases, the number of latent individuals in the population decreases.

In Figure 19 and Figure 20 for different values of $\eta$ which is the failed vaccination rate, $\eta$ = 0.9472, $\eta$ = 0.2, $\eta$ =0.3, and $\eta$ =0.4, the population time in susceptible individuals’ related variation is studied. Here, the most appropriate $\eta$ value for the real data is calculated as $\eta$ = 0.9472 by parameter estimation method. According to Figure 19, it is observed that as the $\eta$ value increases, the number of susceptible individuals in the population also increases. Similarly, in Figure 20, the time variation of the number of vaccinated individuals for different $\eta$ values is examined. According to the simulation results, as $\eta$ values increase over time, the number of vaccinated individuals also decreases. Here, the opposite of the situation observed in Figure 19 is expected to occur. In other words, an increase in $\eta$ values (increasing number of failed vaccines) will cause a decrease in the number of vaccinated individuals and increase the number of susceptible individuals.

9. Concluding Remarks

In this study, mathematical analysis of a developed new Hepatitis-B mathematical model is made using Caputo fractional derivative, and the parameter values of this model are estimated with real data from Türkiye. Numerical simulations are carried out using the estimated parameter values for Equation (10) system and the future processes of the Hepatitis-B epidemic are tried to be predicted. In this part of the paper, the results and recommendations obtained throughout the study are presented.

• A new mathematical structure modeling Hepatitis-B disease is developed by using appropriate parameters, and the efficiency and accuracy of this model are examined. The structured model consists of susceptible (S), carrier (C) and recovered (R) individuals, which are considered to be the most basic components of Hepatitis-B disease, as well as latent (L), includes acute (A) and vaccinated (V) individuals. In this sense, it can be said that the model created is a very effective and productive model for Hepatitis-B. It is seen that the model created as a result of the examination intuitively models the processes of the Hepatitis-B disease and provides predictions about its future course.
• In order to make a more detailed analysis of the developed model and to take into account the memory effect, the fractional derivative is used to account for the memory effect. Considering the results obtained with the help of graphical methods, it is seen that the fractional derivative gives more meaningful results than the classical (integer) order. Thus, the connection between the fractional and integer orders in terms of the future course of Hepatitis-B disease is revealed.
• The non-negative solution region and the limitations of the model’s compartments are discussed to show the biological significance of the system forming the model. In addition, the existence and uniqueness of the solution of the relevant system are examined. Thus, the necessary conditions are obtained in the system created for the solution to exist and be unique.
• The model’s equilibrium points for diseased (endemic) and disease-free states are computed, and an investigation of their stability is performed. Thus, it is established under which circumstances the system’s disease-free equilibrium points are stable.
• The fundamental reproduction number, also referred to as the secondary infection rate in epidemics, is calculated to be ${\mathcal{R}}_0=0.000000035447$ and provides crucial information about how the disease will develop in the future.
• The parameters of the Hepatitis-B model are estimated by the "least squares curve fitting" method. Numerical simulations are run in accordance with these estimated values using actual data from Türkiye. Numerical simulations are used to forecast how the Hepatitis-B disease would progress in the future and to determine how the parameters affect each compartment.

The following suggestions can be offered to guide researchers in future similar studies:

• Prediction course of Hepatitis-B can be considered by using different fractional derivative operators, such as Riemann–Liouville, Caputo–Fabrizio or Atangana–Baleanu instead of the Caputo fractional derivative operator, which is already used in the study. Thus, the relationship between different derivative operators can be revealed.
• More data sets can be used for parameter estimation on different types of biological models. In addition to the least squares curve-fitting method, other methods, such as maximum likelihood, can be used.