EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 1, Article Number 5509 ISSN 1307-5543 – ejpam.com Published by New York Business Global Numerical Investigation of the Dynamical Behavior of Hepatitis B Virus via Caputo-Fabrizio Fractional Derivative Imtiaz Ahmad1, Rashid Jan2,3, Normy Norfiza Abdul Razak2, Aziz Khan4, Thabet Abdeljawad5,4,6,7 1 Institute of Informatics and Computing in Energy (IICE), Universiti Tenaga Nasional, Kajang, Selangor, Malaysia 2 Institute of Energy Infrastructure (IEI), Department of Civil Engineering, College of Engineering, Universiti Tenaga Nasional (UNITEN), Putrajaya Campus, Jalan IKRAM-UNITEN, 43000 Kajang, Selangor, Malaysia 3 Mathematics Research Center, Near East University TRNC, Mersin 10, Nicosia 99138, Turkey 4 Department of Mathematics and Sciences, Prince Sultan University, P.O. Box 66833, 11586 Riyadh, Saudi Arabia 5 Department of Mathematics, Saveetha School of Engineering, Saveetha Institute of Medical and Technical Sciences, Saveetha University, Chennai 602105, Tamil Nadu, India 6 Center for Applied Mathematics and Bioinformatics (CAMB), Gulf University for Science and Technology, Hawally, 32093, Kuwait 7 Department of Mathematics and Applied Mathematics, Sefako Makgatho Health Sciences University, Garankuwa, Medusa 0204, South Africa Abstract. Hepatitis B is a viral infection that primarily targets the liver, potentially leading to acute or chronic liver diseases with severe complications, such as cirrhosis and liver cancer. Its persistent prevalence underscores its status as a significant global health issue. This research constructs a mathematical model for the progression of Hepatitis B using fractional derivatives, accounting for a two-dose vaccine regimen. The basic concepts of the Caputo-Fabrizio (CF) deriva- tive are presented for the model’s analysis. The infection-free equilibrium is investigated, and the endemic indicator of the system, R0, is determined using the next-generation matrix method. The model exhibits local asymptotic stability at the infection-free equilibrium when R0 < 1, and insta- bility otherwise. Conditions ensuring the existence and uniqueness of solutions for the proposed fractional dynamics model are established. A new numerical method for analyzing the time series of the system is also presented. The study elucidates the impact of input variables on the system’s dynamic behavior and identifies critical factors within the model, highlighting key parameters that can be targeted for the control and management of Hepatitis B infection. 2020 Mathematics Subject Classifications: 92D30, 34A08, 34D20, 37N25, 65L05 DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5509 Email addresses: imtiazkakakhil@gmail.com (I. Ahmad), rashid.jan@uniten.edu.my (R. Jan), tabdeljawad@psu.edu.sa (T. Abdeljawad) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 2 of 27 Key Words and Phrases: Hepatitis B, fractional dynamics, Caputo-Fabrizio derivative, model analysis, dynamical behavior 1. Introduction Several globally significant infectious diseases, including hepatitis B (HB) [22], pose a threat to human life. Hepatitis B results from an infection instigated by the hepati- tis B virus (HBV), leading to hepatic impairment and the potential progression of both acute and chronic infections. Despite the asymptomatic nature in the early phases of HB infection for many individuals, certain cases may subsequently manifest abrupt clin- ical symptoms, encompassing nausea, jaundice, fatigue, darkened urine, and abdominal discomfort [40]. Roughly one-third of the world’s population is estimated to experience infection at some point in their lives, with 240-350 million individuals afflicted by chronic infections [40]. In comparison to the extensive global effort directed towards HIV by pub- lic health planners, it is noteworthy that HBV and HCV, while equally transmissible, are significantly more infectious-100 times more for HBV and 10 times more for HCV [30]. In areas with a significant prevalence of HBV, the primary transmission routes involve the transfer of the virus from mother to child, exposure to contaminated blood or bodily fluids through percutaneous means, and sexual contact [30]. In the year 2015, there were ap- proximately 900,000 fatalities attributed to HBV infection, predominantly resulting from the development of hepatocellular carcinoma (HCC) and cirrhosis [29]. The implemen- tation of HB birth immunization is pivotal in preventing the transmission of HBV from mother to child and between individuals. It is imperative for healthcare systems to enhance their efforts in the treatment and prevention of viral hepatitis [36]. Additionally, current investigations are exploring the application of antiviral treatments, including tenofovir, telbivudine, and lamivudine, as prospective approaches to alleviate perinatal transmission of HBV in expectant women exhibiting elevated HBV DNA levels. Screening for HBV is crucial due to the highly contagious nature of the disease and its severe public health implications. Emphasizing vaccination as a preventive measure against HBV infection and its associated complications is imperative [36, 40]. The primary objectives of HBV screening encompass the detection of individuals with chronic infections, who could benefit from treatment or receive education on lifestyle modifications to reduce transmission risks. Furthermore, it strives to pinpoint close contacts of infected individuals for vaccination and offer care for chronic HB cases, along with ongoing monitoring of disease activity and surveillance for HCC [32]. Focused screening of at-risk populations, such as migrants, could enhance the identification of cases. A cost-effectiveness evaluation of screening for viral hepatitis among immigrants originating from regions with moderate to high prevalence emphasizes the significance of screening as a secondary prevention measure [33]. Numerous scholars have explored the patterns of HB transmission in different geographical areas, as well as the immune system’s response to infection, through the application of mathematical models. Anderson and May, in their research, elucidated the outcomes of HB transmission among carriers by employing a basic mathematical model I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 3 of 27 [11]. Anderson and Williams, in their modeling of sexual transmission of HB, integrated a range of variables encompassing sexual behavior and age [41]. The age of HB infection and the development of the carrier state were found to be related by Edmunds et al. [21]. Medley et al. introduced a feedback mechanism model that links the probability of carrier- class formation after infection to factors such as transmission rate, average age of infection, and infection prevalence [31], whereas Din and Abidin [18] used the Mittag-Leffler kernel to model an epidemic of HB vaccination. Transmission dynamics of hepatitis B differ between industrialized nations and de- veloping countries. In their study, Edmunds et al. [19] primarily considered childhood HB transmission rates as a key determinant of epidemic levels in developing countries, complemented by the assessment of sexual contact rates. It is noteworthy that data on sexual contact rates in developing nations are limited. While Williams et al. [20] provided a mathematical model for assessing the transmission dynamics of HB in the UK. Coupe et al. [17] constructed a standard model of immune response to investigate the temporal aspects of effector cell activation, the course of the disease, the induction of antibodies by the vaccine, and the role of preexisting conditions in preventing HB. Subsequently, Gourley et al. later expanded on this model with a time-delayed variation, while Min et al. used a traditional model function instead of mass action to elucidate susceptibility to HB infection [23]. In order to improve the precision of the baseline reproduction number estimations for the HB population, Hews et al. employed a logistic growth model and a conventional methodology, thereby enhancing the model’s concordance with the available data [25]. In this study, we will formulate a mathematical model to analyze the trans- mission dynamics of Hepatitis B (HB), incorporating double-dose vaccination and the presence of asymptomatic carriers. The model will account for the role of asymptomatic individuals in disease spread and evaluate the effectiveness of a two-dose vaccination strat- egy in controlling the infection. By integrating these factors, the proposed model aims to provide a comprehensive understanding of HB transmission and inform public health strategies for disease prevention and management. The application of fractional derivatives extends beyond theoretical mathematics and finds relevance in various real-world problems across diverse fields [4–7, 10, 14, 26, 34]. These applications highlight the versatility of fractional calculus in addressing real-world problems across disciplines, offering a powerful tool for modeling complex phenomena and improving the understanding of intricate dynamic systems [1–3, 8, 9, 15, 37, 39]. The uti- lization of fractional derivatives in epidemic models represents a sophisticated approach to capturing non-integer order dynamics inherent in certain biological and epidemiological systems [27, 35]. Unlike traditional integer-order derivatives, fractional derivatives involve non-integer orders, allowing for a more nuanced representation of complex behaviors and long-range dependencies observed in epidemic spread. In the context of epidemic mod- eling, the fractional derivative is employed to describe the rate of change of infected or susceptible populations with respect to time, considering fractional orders that may re- flect memory effects or non-local interactions. This mathematical framework offers a more accurate portrayal of epidemic dynamics, particularly when conventional models fall short in capturing the intricate patterns exhibited by real-world epidemics. Therefore, we have I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 4 of 27 chosen to model the dynamics of Hepatitis B (HB) using a fractional-order framework to achieve more accurate and precise results. The fractional approach allows for better rep- resentation of the memory effects and complex dynamics inherent in disease transmission, providing enhanced insights into the behavior of the infection and the impact of control measures. The present study is organized as follows: Section 2 elucidates definitions and con- cepts pertinent to fractional theory in the context of this research. Section 3 introduces the developed HB model. Section 4 employs Caputo-Fabrizio derivatives to portray the dynamics, computes the threshold parameter R0, and investigates the equilibria. Sec- tion 5 delves into the uniqueness and existence of solutions, employing fixed point theory. Finally, Section 6 provides a synthesis of the findings and concludes the research. 2. Theory and concepts Here, we introduce the fundamental concept of the Caputo-Fabrizio (CF) fractional operator for the examination of our proposed HB model. The theory and principles un- derlying the fractional CF derivative are outlined as follows: Definition 1. Assume a function k ∈ H1(a, b), subject to the condition that a is less than b. In this case, the CF derivative is given by Dς t (k(t)) = U(ς) 1− ς ∫ t a k′(x) exp ( − ς t− x 1− ς ) dx, (1) with the condition that ς ∈ [0, 1] and U(τ) indicates the normality holding U(0) = U(1) = 1 [16]. Otherwise, if k /∈ H1(a, b), then Dς t (k(t)) = ςU(ς) 1− ς ∫ t a (k(t)− k(x)) exp ( − ς t− x 1− ς ) dx. (2) Remark 1. If α = 1−ς ς ∈ [0,∞) and ς = 1 1+α ∈ [0, 1], then the above (2) suggest Dς t (k(t)) = M(α) α ∫ t a k′(x)e[− t−x α ]dx, M(0) = M(∞) = 1. (3) Moreover, we have lim α−→0 1 α exp ( − t− x α ) = δ(x− t). (4) Definition 2. [28]. The fractional integral of a given function is expressed as follows. Iςt (k(t)) = 2(1− ς) (2− ς)U(ς) k(t) + 2ς (2− ς)U(ς) ∫ t 0 k(u)du, t ≥ 0. (5) Next, we move to the following Remark: I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 5 of 27 Remark 2. The Definition 2 provided above yields the following result. 2(1− ς) (2− ς)U(ς) + 2ς (2− ς)U(ς) = 1, (6) with the condition that U(ς) = 2 2−ς . In [28], new operator has been developed by utilizing (6) as follows Dς t (k(t)) = 1 1− ς ∫ t 0 k′(x) exp ( ς t− x 1− ς ) dx. (7) In this context, the fractional order adheres to the requirement that 0 < ς < 1. In the literature, various authors have introduced several fractional operators. Some novel fractional operators are defined below: Definition 3. The Atangana-Beleanu fractional operator in the Caputo form is given by: ABC p Dς t g(t) = B(ς) 1− ς ∫ t p g′(κ)Eς [ − ς (t− κ)ς 1− ς ] dκ. in which g ∈ H1(p, q), q > p, and ∈ [0, 1]. Definition 4. [24]. The He’s fractional operator is defined as follows: Dς t g(t) = 1 Γ(1− ς) dn dtn ∫ t t0 (s− t)n−ς−1[g0(s)− g(s)]ds. where ς is the fractional order and g0 is a known function. 3. Formulation of the model We constructed a compartmental mathematical model to simulate the dynamics of hepatitis B transmission based on various transmission routes. The total population, denoted as N(t), is partitioned into six distinct classes: S(t), V1(t), V2(t), I(t), C(t), and R(t). Consequently, N(t) can be described as follows: N(t) = S(t) + V1(t) + V2(t) + I(t) + C(t) +R(t). (8) Additionally, S(t) stands for individuals who have not encountered the infection yet and remain susceptible to infection. V1(t) denotes those who have been administered the initial vaccine dose, whereas V2(t) represents individuals who have successfully received both doses of the vaccine, usually with a one-month interval between them. I(t) denotes those who are in the infectious acute stage of HB, which typically spans the first six months following infection. C(t) represents individuals in the chronic infectious stage, and R(t) designates those who have successfully recovered from hepatitis B. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 6 of 27 The underlying assumptions of this model include the transmission of HB via both vertical and horizontal methods. Additionally, the model disregards infection from acutely infected mothers, as the acute infection stage is relatively short compared to the duration of pregnancy. The model accounts for disease-induced mortality within the chronic class and incorporates a natural mortality rate for each compartment. Furthermore, it postulates a constant generation of new susceptible individuals, represented by the symbol Π. All parameters involved in the model are nonnegative. Individuals at risk can contract the infection through horizontal transmission, driven by the force of infection represented as λ = ξ(κA + C). Here, the effective contact rate is denoted by ξ, while the relative infectiousness of acute infections versus chronic ones is represented by κ. When κ > 1, it indicates that acute infections have a higher likelihood of transmission to susceptible individuals than chronic infections. Both acute and chronic infections have an identical chance of spreading to individuals who are susceptible when κ = 1. In contrast, if κ < 1, it means that chronic infections are more probable than acute infections to spread to individuals who are susceptible. In the scenario of vertical transmission, we make the assumption that newborn infants receive successful vaccination at a rate represented by Πω, or unsuccessful vaccination at a rate represented by Π(1− ω). Here, the birth rate is denoted by Π, and the proportion of effectively immunized individuals is denoted by ω, where 0 < ω < 1. Consequently, Πω individuals transition to the vaccinated first dose class, while Π(1 − ω)τC move into the chronic class. The remaining fraction, Π(1− ω)(1− τC), enters the susceptible class. In this context, the symbol τ denotes the fraction of infants born to carrier mothers who remain unvaccinated and subsequently develop a chronic condition. Furthermore, we have taken into account that the percentage of unvaccinated newborns moving into a chronic state is lower than the sum of the mortality rates from hepatitis B and the percentage of patients who recover from the chronic category. Furthermore, the model encompasses a range of parameters, each with specific inter- pretations: Π represents the birth rate, µ signifies the natural death rate, σ denotes the natural recovery rate from the acute class, ϑ characterizes the rate of recovery from the chronic class due to treatment, δ represents the speed at which individuals move from the acute class to the chronic class once they exhibit symptoms, α represents the rate at which individuals in the initially vaccinated first dose class receive their second vaccine dose whereas the rate of recovery from the second dose class of vaccination is represented by ϵ. ψ represents the disease-related death rate in the chronic class, φ1 designates the rate at which members of the susceptible population receive their first dosage of the vaccine, and φ2 signifies the rate at which vaccine immunity diminishes in the vaccinated first dose I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 7 of 27 class. The system of ODEs, under the specified assumptions, is expressed as follows [13]: dS dt = Π(1− ω)(1− τC) + φ2V1 − (λ+ µ+ φ1)S, dV1 dt = φ1S +Πω − (α+ µ+ φ2)V1, dV2 dt = αV1 − (ϵ+ µ)V2, dA dt = λS − (δ + σ + µ)A, dC dt = δA+Π(1− ω)τC − (ϑ+ µ+ ψ)C, dR dt = σA+ ϑC + ϵV2 − µR, (9) with the following values 0 ≤ S(0), 0 ≤ V1(0), 0 ≤ V2(0) 0 ≤ A(0), 0 ≤ C(0), 0 ≤ R(0). This incorporation facilitates a comprehensive examination of epidemic dissemination in complex scenarios featuring memory effects and long-range interactions, which enhances our understanding of the underlying mechanisms. Consequently, we have applied a frac- tional framework to describe the dynamics of HB disease. In our proposed model, we have employed the CF fractional operator in the following manner. CF 0 Dς tS = Π(1− ω)(1− τC) + φ2V1 − (λ+ µ+ φ1)S, CF 0 Dς tV1 = φ1S +Πω − (α+ µ+ φ2)V1, CF 0 Dς tV2 = αV1 − (ϵ+ µ)V2, CF 0 Dς tA = λS − (δ + σ + µ)A, CF 0 Dς tC = δA+Π(1− ω)τC − (ϑ+ µ+ ψ)C, CF 0 Dς tR = σA+ ϑC + ϵV2 − µR, (10) with positive initial conditions. 3.1. Analysis of Disease-Free Equilibrium The system (10) reaches a disease-free equilibrium (DFE) in the absence of infections. To find the DFE, one must solve the system by setting the left-hand side of equation (10) to zero. At the DFE, both A and C are equal to zero, causing system (10) to simplify to 0 = Π(1− ω) + φ2V1 − (λ+ µ+ φ1)S, 0 = φ1S +Πω − (α+ µ+ φ2)V1, 0 = αV1 − (ϵ+ µ)V2, 0 = ϵV2 − µR, (11) therefore, based on (11), the DFE for model (10) can be expressed as E0 = ( S0, V 0 1 , V 0 2 , A 0, C0, R0 ) , = ( Π(l−ωv) µl+φ1v , Π(φ1+ωµ) µl+φ1v , Πα(φ1+ωµ) (ϵ+µ)(µl+φ1v) , 0, 0, Πεα(φ1+ωµ) (ε+µ)(µl+φ1v) ) , (12) where l = α+ µ+ φ2, v = α+ µ, and α+ µ < ω(α+ µ). I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 8 of 27 The basic reproduction number, R0, for the system described in (10) is calculated using the approach outlined in [38]. This parameter is defined as R0 = ρ(FV −1), where ρ(FV −1) represents the spectral radius of the matrix FV −1. In our analysis, we focus only on the infectious classes within the differential equation system (10), which can be expressed as: dA dt = λS − (δ + σ + µ)A, dC dt = δA+Π(1− ω)τC − (ϑ+ µ+ ψ)C. (13) Consider f to represent the rate at which new infections enter the system, and let v quantify the disparity between the immigration rate into the compartments and the rate of inter-compartmental transfer. f = [ λS 0 ] , and v = [ (δ + σ + µ)A (ϑ+ µ+ ψ)C −Π(1− ω)τC − δA ] , where λ = ξ(κA + C). Next, calculate the Jacobian matrix for both f and v concerning the variables A and C at the DFE denoted as E0. F = [ κ 1 0 0 ] , and V = [ δ + σ + µ 0 −δ (ϑ+ µ+ ψ)−Π(1− ω)τ ] . Next, by calculating the matrix V inverse, we obtain the following result. V −1 = 1 (δ + σ + µ)((ϑ+ µ+ ψ)−Π(1− ω)τ) [ (ϑ+ µ+ ψ)−Π(1− ω)τ 0 δ δ + σ + µ ] , then FV −1 is given by FV −1 = [ κξS0 δ+σ+µ + ξδS0 (δ+σ)((ϑ+µ+ψ)−Π(1−ω)τ) ξS0 (ϑ+µ+ψ)−Π(1−ω)τ 0 0 ] . The endemic indicator R0 corresponds to the highest eigenvalue of the matrix mentioned above. R0 = ξS0 ( κ δ+σ+µ + δ (δ+σ+µ)((ϑ+µ+ψ)−Π(1−ω)τ) ) . (14) Theorem 1. The state of DFE exhibits local asymptotic stability when R0 < 1, and it becomes unstable when R0 exceeds 1. Proof. The Jacobian matrix of the system described in (10) at the equilibrium point denoted as E0 is given by J(E0) =  −(µ+ φ1) φ2 0 −ξκS0 −Π(1− ω)τ − ξS0 0 φ1 −(α+ µ+ φ2) 0 0 0 0 0 α −(µ+ ϵ) 0 0 0 0 0 0 ξκS0 − (δ + σ + µ) ξS0 0 0 0 0 δ Π(1− ω)τ − (ϑ+ µ+ ψ) 0 0 0 ϵ σ ϑ −µ  , I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 9 of 27 then the characteristic equation of the above is |J(E0)− λ| = 0. From which, we derived λ1 = −µ < 0, λ2 = −(µ+ ϵ) < 0. Furthermore, employing the Routh-Hurwitz criteria, we established that both λ3 and λ4 are negative. As for the remaining two eigenvalues, they constitute the eigenvalues of the given system. J1 = [ ξκS0 − (δ + σ + µ) ξS0 δ Π(1− ω)τ − (ϑ+ µ+ ψ) ] . We can now employ the standard criteria to guarantee that the eigenvalues of J1 exhibit a negative real component. Specifically, we seek a condition where both the trace of J1 < 0, and the determinant of J1 > 0. From the second set of inequalities, the following condition has been derived: (ξκS0 − (δ + σ + µ))(Π(1− ω)τ − (ϑ+ µ+ ψ))− δξS0 > 0, −(δ + σ + µ)(Π(1− ω)τ − (ϑ+ µ+ ψ)) > (ϑ+ µ+ ψ −Π(1− ω)τ)ξκS0 + δξS0, 1 > ξS0 ( κ δ+σ+µ + δ (δ+σ+µ)(ϑ+µ+ψ−Π(1−ω)τ) ) = R0. (15) Observing the condition R0 < 1, we can determine that the trace of J1 is negative (J1 < 0), and the determinant of J1 is positive (J1 > 0). These indications imply that the remaining two eigenvalues have negative real parts. Consequently, when R0 < 1, the DFE point exhibits local asymptotic stability, while for R0 > 1, the DFE point becomes unstable. 3.2. Existence theory Our focus is on analyzing the solutions of the proposed fractional HBV system. The existence of the solution (10) is investigated using fixed point theory. The procedure is as follows: S(t)− S(0) = CF 0 Iςt { Π(1− ω)(1− τC) + φ2V1 − (λ+ µ+ φ1)S } , V1(t)− V1(0) = CF 0 Iςt { φ1S +Πω − (α+ µ+ φ2)V1 } , V2(t)− V2(0) = CF 0 Iςt { αV1 − (ϵ+ µ)V2 } , A(t)−A(0) = CF 0 Iςt { λS − (δ + σ + µ)A } , C(t)− C(0) = CF 0 Iςt { δA+Π(1− ω)τC − (ϑ+ µ+ ψ)C } , R(t)−R(0) = CF 0 Iςt { σA+ ϑC + ϵV2 − µR } . (16) After that, we have S(t)− S(0) = 2(1− ς) (2− ς)U(ς) { Π(1− ω)(1− τC) + φ2V1 − (λ+ µ+ φ1)S } I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 10 of 27 + 2ς (2− ς)U(ς) ∫ t 0 { Π(1− ω)(1− τC) + φ2V1 − (λ+ µ+ φ1)S } dy, V1(t)− V1(0) = 2(1− ς) (2− ς)U(ς) { φ1S +Πω − (α+ µ+ φ2)V1 } + 2ς (2− ς)U(ς) ∫ t 0 { φ1S +Πω − (α+ µ+ φ2)V1 } dy, V2(t)− V2(0) = 2(1− ς) (2− ς)U(ς) { αV1 − (ϵ+ µ)V2 } + 2ς (2− ς)U(ς) ∫ t 0 { αV1 − (ϵ+ µ)V2 } dy, A(t)−A(0) = 2(1− ς) (2− ς)U(ς) { λS − (δ + σ + µ)A } + 2ς (2− ς)U(ς) ∫ t 0 { λS − (δ + σ + µ)A } dy, C(t)− C(0) = 2(1− ς) (2− ς)U(ς) { δA+Π(1− ω)τC − (ϑ+ µ+ ψ)C } + 2ς (2− ς)U(ς) ∫ t 0 { δA+Π(1− ω)τC − (ϑ+ µ+ ψ)C } dy, R(t)−R(0) = 2(1− ς) (2− ς)U(ς) { σA+ ϑC + ϵV2 − µR } + 2ς (2− ς)U(ς) ∫ t 0 { σA+ ϑC + ϵV2 − µR } . (17) Moreover,  L1(t, S) = Π(1− ω)(1− τC) + φ2V1 − (λ+ µ+ φ1)S, L2(t, V1) = φ1S +Πω − (α+ µ+ φ2)V1, L3(t, V2) = αV1 − (ϵ+ µ)V2, L4(t, A) = βC + γI − (d+ τ)R, L5(t, C) = δA+Π(1− ω)τC − (ϑ+ µ+ ψ)C, L6(t, R) = σA+ ϑC + ϵV2 − µR. (18) Theorem 2. The kernels L1,L2,L3,L4,L5, and L6 satisfy the conditions of Lipschitz and contraction if the following criteria are met. 0 ≤ (ξ(κM +M) + µ+ φ1) < 1. Proof. To achieve the specified outcomes, we consider the variables denoted as S and S1, initiating the process from the initial state L1 is as follows: L1(t, S)− L1(t, S1) = −λ{S(t)− S(t1)} − µ{S(t)− S(t1)} I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 11 of 27 −φ1{S(t)− S(t1)}. (19) Upon simplifying equation (19), we attain the below ∥L1(t, S)− L1(t, S1)∥ ≤ ∥λ{S(t)− S(t1)}||+ µ∥{S(t)− S(t1)}∥ +φ1∥{S(t)− S(t1)}∥ ≤ ∥ξ(κA+ C)∥∥{S(t)− S(t1)}||+ µ∥{S(t)− S(t1)}∥ +φ1∥{S(t)− S(t1)}∥ ≤ ξ(κM +M)∥{S(t)− S(t1)}||+ µ∥{S(t)− S(t1)}∥ +φ1∥{S(t)− S(t1)}∥ ≤ (ξ(κM +M) + µ+ φ1)∥{S(t)− S(t1)}∥. Taking ϖ1 = (ξ(κM +M)+µ+φ1), where ∥A∥ ≤M and ∥C∥ ≤M due to boundedness, we have ||L1(t, S)− L1(t, S1)|| ≤ ϖ1||S(t)− S(t1)||. (20) Therefore, we have established the Lipschitz condition for L1, and additionally, contraction is derived from the condition 0 ≤ (M + d+ ν) < 1. Similarly, the Lipschitz conditions can be determined as follows ||L2(t, V1)− L2(t, V11)|| ≤ ϖ2||V1(t)− V1(t1)||, ||L3(t, V2)− L3(t, V21)|| ≤ ϖ3||V2(t)− V2(t1)||, ||L4(t, A)− L4(t, A1)|| ≤ ϖ4||A(t)−A(t1)||, ||L5(t, C)− L5(t, C1)|| ≤ ϖ5||C(t)− C(t1)||, ||L6(t, R)− L6(t, R1)|| ≤ ϖ6||R(t)−R(t1)||. (21) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 12 of 27 Utilizing (17), the following can be derived: Sm(t) = 2 (1−ς) (2−ς)U(ς)L1(t, S(m−1)) + 2 ς (2−ς)U(ς) ∫ t 0(L1(z, S(m−1)))dz, V1m(t) = 2 (1−ς) (2−ς)U(ς)L2(t, V1(m−1)) + 2 ς (2−ς)U(ς) ∫ t 0(L2(z, V1(m−1)))dz, V2m(t) = 2 (1−ς) (2−ς)U(ς)L3(t, V2(m−1)) + 2 ς (2−ς)U(ς) ∫ t 0(L3(z, V2(m−1)))dz, Am(t) = 2 (1−ς) (2−ς)U(ς)L4(t, A(m−1)) + 2 ς (2−ς)U(ς) ∫ t 0(L4(z,A(m−1)))dz, Cm(t) = 2 (1−ς) (2−ς)U(ς)L5(t, C(m−1)) + 2 ς (2−ς)U(ς) ∫ t 0(L5(z, C(m−1)))dz, Rm(t) = 2 (1−ς) (2−ς)U(ς)L6(t, C(m−1)) + 2 ς (2−ς)U(ς) ∫ t 0(L6(z,R(m−1)))dz, (22) with the below initial values S0(t) = S(0), V 0 1 (t) = V1(0), V 0 2 (t) = V2(0), A 0(t) = A(0), C0(t) = C(0), R0(t) = R(0). The difference expressions are derived in the following manner. κ1m(t) = Sm(t)− S(m−1)(t) = 2(1− ς) (2− ς)U(ς) (L1(t, S(m−1))− L1(t, S(m−2))) +2 ς (2− ς)U(ς) ∫ t 0 (L1(z, S(m−1))− L1(z, S(m−2)))dz, κ2m(t) = V1m(t)− V1(m−1)(t) = 2(1− ς) (2− ς)U(ς) (L2(t, V1(m−1))− L2(t, V1(m−2))) +2 ς (2− ς)U(ς) ∫ t 0 (L2(z, V1(m−1))− L2(z, V1(m−2)))dz, κ3m(t) = V2m(t)− V2(m−1)(t) = 2(1− ς) (2− ς)U(ς) (L3(t, V2(m−1))− L3(t, V2(m−2))) +2 ς (2− ς)U(ς) ∫ t 0 (L3(z, V2(m−1))− L3(z, V2(m−2)))dz, κ4m(t) = Am(t)−A(m−1)(t) = 2(1− ς) (2− ς)U(ς) (L4(t, A(m−1))− L4(t, A(m−2))) +2 ς (2− ς)U(ς) ∫ t 0 (L4(z,A(m−1))− L4(z,A(m−2)))dz, κ5m(t) = Cm(t)− C(m−1)(t) = 2(1− ς) (2− ς)U(ς) (L5(t, C(m−1))− L5(t, C(m−2))) +2 ς (2− ς)U(ς) ∫ t 0 (L5(z, C(m−1))− L5(z, C(m−2)))dz, κ6m(t) = Rm(t)−R(m−1)(t) = 2(1− ς) (2− ς)U(ς) (L6(t, R(m−1))− L6(t, R(m−2))) +2 ς (2− ς)U(ς) ∫ t 0 (L6(z,R(m−1))− L6(z,R(m−2)))dz. (23) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 13 of 27 Noticing Sm(t) = ∑m j=1 κ1j(t), V1n(t) = ∑m j=1 κ2j(t), V2m(t) = ∑m j=1 κ3J(t), Am(t) = ∑m j=1 κ4j(t), Cm(t) = ∑m j=1 κ5j(t), Rm(t) = ∑m j=1 κ6j(t). (24) Similarly, ||κ1m(t)|| = ||Sm(t)− S(m−1)(t)|| = ∥∥∥∥2 (1− ς) (2− ς)U(ς) (L1(t, S(m−1))− L1(t, S(m−2))) +2 ς (2− ς)U(ς) ∫ t 0 (L1(z, S(m−1))− L1(z, S(m−2)))dz ∥∥∥∥. (25) Equation (25) suggests that ∥Sm(t)− S(m−1)(t)∥ ≤ 2 (1− ς) (2− ς)U(ς) ∥(L1(t, S(m−1))− L1(t, S(m−2)))∥ +2 ς (2− ς)U(ς) ∥∥∥∥∫ t 0 (L1(z, S(m−1))− L1(z, S(m−2)))dz ∥∥∥∥.(26) The aforementioned results in ∥Sm(t)− Sm−1(t)∥ ≤ 2 (1− ς) (2− ς)U(ς) ϖ1∥S(m−1) − S(m−2)∥+ 2 ς (2− ς)U(ς) ϖ1 × ∫ t 0 ∥S(m−1) − S(m−2)∥dz. (27) Furthermore ∥κ1m(t)∥ ≤ 2 (1− ς) (2− ς)U(ς) ϖ1∥κ1(m−1)(t)∥+ 2 ς (2− ς)U(ς) ϖ1 ∫ t 0 ∥κ1(m−1)(z)∥dz.(28) Similarly ∥κ2m(t)∥ ≤ 2 (1− ς) (2− ς)U(ς) ϖ2∥κ2(m−1)(t)∥+ 2 ς (2− ς)U(ς) ϖ2 ∫ t 0 ∥κ2(m−1)(z)∥dz, ∥κ3m(t)∥ ≤ 2 (1− ς) (2− ς)U(ς) ϖ3∥κ3(m−1)(t)∥+ 2 ς (2− ς)U(ς) ϖ1 ∫ t 0 ∥κ3(m−1)(z)∥dz, ∥κ4m(t)∥ ≤ 2(1− ς) (2− ς)U(ς) ϖ4∥κ4(m−1)(t)∥+ 2 ς (2− ς)U(ς) ϖ4 ∫ t 0 ∥κ4(m−1)(z)∥dz, ∥κ5m(t)∥ ≤ 2 (1− ς) (2− ς)U(ς) ϖ5∥κ5(m−1)(t)∥+ 2 ς (2− ς)U(ς) ϖ5 ∫ t 0 ∥κ5(m−1)(z)∥dz, ∥κ6m(t)∥ ≤ 2 (1− ς) (2− ς)U(ς) ϖ6∥κ6(m−1)(t)∥+ 2 ς (2− ς)U(ς) ϖ6 ∫ t 0 ∥κ6(m−1)(z)∥dz.(29) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 14 of 27 Theorem 3. If there exists a t0 such that the following condition satisfies 2 (1− ς) (2− ς)U(ς) ϖ1 + 2 ς (2− ς)U(ς) ϖ1t0 < 1, then, exact coupled solutions for the formulated fractional system (10) are obtained. Proof. Given the satisfaction of the Lipschitz condition and the boundedness of S(t), V1(t), V2(t), A(t), C(t), and R(t), we deduce the following from (28) and (29). ∥κ1m(t)∥ ≤ ||Sm(0)|| [( 2 (1− ς) (2− ς)U(ς) ϖ1 ) + ( 2 ς (2− ς)U(ς) ϖ1t )]m , ∥κ2m(t)∥ ≤ ||V1m(0)|| [( 2 (1− ς) (2− ς)U(ς) ϖ2 ) + ( 2 ς (2− ς)U(ς) ϖ2t )]m , ∥κ3m(t)∥ ≤ ||V2m(0)|| [( 2 (1− ς) (2− ς)U(ς) ϖ3 ) + ( 2 ς (2− ς)U(ς) ϖ3t )]m , ∥κ4m(t)∥ ≤ ||Am(0)|| [( 2 (1− ς) (2− ς)U(ς) ϖ4 ) + ( 2 ς (2− ς)U(ς) ϖ4t )]m , ∥κ5m(t)∥ ≤ ||Cm(0)|| [( 2 (1− ς) (2− ς)U(ς) ϖ5 ) + ( 2 ς (2− ς)U(ς) ϖ5t )]m , ∥κ6m(t)∥ ≤ ||Rm(0)|| [( 2 (1− ς) (2− ς)U(ς) ϖ6 ) + ( 2 ς (2− ς)U(ς) ϖ6t )]m . (30) Consequently, the continuity and existence of solutions are established. Additionally, it is imperative to demonstrate that the aforementioned is a solution of the model (10), and we proceed with the following steps. S(t)− S(0) = Sm(t)−W1m(t), V1(t)− V1(0) = V1m(t)−W2m(t), V2(t)− V2(0) = V2m(t)−W3m(t), A(t)−A(0) = Am(t)−W4m(t), C(t)− C(0) = Cm(t)−W5m(t), R(t)−R(0) = Rm(t)−W6m(t). (31) Subsequently, we proceed to the following stage. ∥Bm(t)∥ = ∣∣∣∣∣∣ 2(1− ς) (2− ς)U(ς) (L1(t, Sm)− L1(t, S(m−1))) + 2ς (2− ς)U(ς) ×∫ t 0 (L1(z, Sm)− L1(z, S(m−1)))dz ∣∣∣∣∣∣, I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 15 of 27 ≤ 2(1− ς) (2− ς)U(ς) ∥(L1(t, Sm)− (L1(t, S(m−1)))∥+ 2ς (2− ς)U(ς) ×∫ t 0 ||(L1(z, S)− L1(z, S(m−1)))||dz, ≤ 2(1− ς) (2− ς)U(ς) ϖ1∥S − S(m−1)∥+ 2ς (2− ς)U(ς) ϖ1∥S − S(m−1)∥t. (32) Moreover, ∥W1m(t)∥ ≤ ( 2(1− ς) (2− ς)U(ς) + 2ς (2− ς)U(ς) t )nm+1 µm+1 1 a. (33) At the temporal point t0, we have ∥W1m(t)∥ ≤ ( 2(1− ς) (2− ς)U(ς) + 2ς (2− ς)U(ς) t0 )m+1 µm+1 1 a. (34) Employing identical procedures and utilizing (34), we have ∥W1m(t)∥ −→ 0, n→ ∞. Similarly, we deduce that the functions W2m(t),W3m(t),W4m(t),W5m(t),W6m(t) con- verge to zero as the parameter m approaches ∞. In order to establish the uniqueness of the solution for the system (10), we assume an alternative solution denoted as (S1(t), V11(t), V21(t), A1(t), C1(t) R1(t)). S(t)− S1(t) = 2(1− ς) (2− ς)U(ς) (L1(t, S)− L1(t, S1)) + 2ς (2− ς)U(ς) ×∫ t 0 (L1(z, S)− L1(z, S1)) dz. (35) Applying a norm to (35), we obtain ∥S(t)− S1(t)∥ ≤ 2(1− ς) (2− ς)U(ς) ∥L1(t, S)− L1(t, S1)∥+ 2ℓ (2− ℓ)U(ℓ) ×∫ t 0 ∥L1(z, Sh)− L1(z, S1h)∥dz. (36) Here, we establish the following via the Lipschitz condition. ∥S(t)− S1(t)∥ ≤ 2(1− ς) (2− ς)U(ς) ϖ1∥Sh(t)− S1(t)∥+ 2ς (2− ς)U(ς) ×∫ t 0 ϖ1t∥S(t)− S1(t)∥dz. (37) This implies that ∥S(t)− S1(t)∥ ( 1− 2(1− ς) (2− ς)U(ς) ϖ1 − 2ς (2− ℓ)U(ς) ϖ1t ) ≤ 0. (38) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 16 of 27 Theorem 4. If the condition expressed by the following equation is satisfied( 1− 2(1− ς) (2− ς)U(ς) ϖ1 − 2ς (2− ς)U(ς) ϖ1t ) > 0. (39) then the system (10) possesses a unique solution. Proof. Assuming that the expression in (39) holds, the subsequent (38) yields the following expression ∥S(t)− S1(t)∥ = 0. (40) Consequently, we obtain S(t) = S1(t). (41) Similarly, we achieve the following V1(t) = V11(t), V2(t) = V21(t), A(t) = A1(t), C(t) = C1(t), R(t) = R1(t). 4. Computational scheme for the model Here we provide the numerical solution for the proposed fractional model (10) utilizing the approach outlined in [12]. The subsequent steps adhere to the procedure outlined for solving (10): S(t)− S(0) = (1− ς) B(ς) F1(t, S) + ς B(ς) ∫ t 0 F1(ς, S)dς. (42) At each time step tm+1, where m = 0, 1, 2, ...,, we acquire S(tm+1)− S0 = 1− ς B(ς) F1(tm, Sm) + ς B(ς) ∫ tm+1 0 F1(t, S)dt. (43) The difference between successive terms is expressed as follows: Sm+1 − Sm = 1− ς B(ς) {F1(tm, Sm)−F1(tm−1, Sm−1)}+ ς B(ς) ∫ tm+1 tn F1(t, S)dt. (44) Within the close interval [tk, t(k+1)], one can estimate the function F1(t, S) by employing an interpolation polynomial. Pk(t) ∼= f(tk, zk) h (t− tk−1)− f(tk−1, zk−1) h (t− tk), (45) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 17 of 27 where h = tm − tm−1. By applying the previously mentioned polynomial approximation, we can compute the integral in equation (44).∫ tm+1 tm F1(t, S)dt = ∫ tm+1 tm F1(tm, Sm) h (t− tm−1)− F1(tm−1, Sm−1) h (t− tm)dt = 3h 2 F1(tm, Sm)− h 2 F1(tm−1, Sm−1). (46) Substituting (46) into (44) and subsequent simplification yielded the following result. Sm+1 = Sm + ( 1− ς B(ς) + 3h 2B(ς) ) F1(tm, Sm)− ( 1− ς B(ς) + ςh 2B(ς) ) F1(tm−1, Sm−1).(47) Similarly, for the remaining equations within the system (10), we derived the recursive formula as presented below. V1(m+1) = V10 + ( 1− ς B(ς) + 3h 2B(ς) ) F2(tm, V1m)− ( 1− ς B(ς) + ςh 2B(ς) ) F2(tm−1, V1(m−1)), V2(m+1) = V20 + ( 1− ς B(ς) + 3h 2B(ς) ) F3(tm, V2m)− ( 1− ς B(ς) + ςh 2B(ς) ) F3(tm−1, V2(m−1)), Am+1 = A0 + ( 1− ς B(ς) + 3h 2B(ς) ) F4(tm, Am)− ( 1− ς B(ς) + ςh 2B(ς) ) F4(tm−1, Am−1), Cm+1 = C0 + ( 1− α B(ς) + 3h 2B(ς) ) F5(tm, Cm)− ( 1− ς B(ς) + ςh 2B(ς) ) F5(tm−1, Cm−1), Rm+1 = R0 + ( 1− ς B(ς) + 3h 2B(ς) ) F6(tm, Rm)− ( 1− ς B(ς) + ςh 2B(ς) ) F6(tm−1, Rm−1).(48) This numerical scheme will be utilized to examine the dynamic behavior of the proposed model under different input parameters. This investigation aims to predict the spread of infectious diseases within populations, providing critical insights into their propagation patterns. Understanding these dynamics is essential for developing strategies to control or mitigate the disease’s impact. Such findings equip public health officials with the knowledge needed to implement effective measures to reduce infection rates and alleviate the associated public health burden. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 18 of 27 Time in days 0 10 20 30 40 50 60 O n e d o s e v a c c in a te d 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 ζ= 0.85 ζ=0.90 ζ= 0.95 ζ = 1.00 (a) Time in days 0 10 20 30 40 50 60 T w o d o s e v a c c in a te d 0.25 0.3 0.35 0.4 0.45 0.5 0.55 ζ= 0.85 ζ=0.90 ζ= 0.95 ζ = 1.00 (b) Time in days 0 10 20 30 40 50 60 A c u te i n fe c ti o u s 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ζ= 0.85 ζ=0.90 ζ= 0.95 ζ = 1.00 (c) Time in days 0 10 20 30 40 50 60 C h ro n ic i n fe c ti o u s 0.2 0.3 0.4 0.5 0.6 0.7 0.8 ζ= 0.85 ζ=0.90 ζ= 0.95 ζ = 1.00 (d) Time in days 0 10 20 30 40 50 60 R e c o v e re d i n d iv id u a ls 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 ζ= 0.85 ζ=0.90 ζ= 0.95 ζ = 1.00 (e) Figure 1: Visualization of the solution trajectories of vaccinated, acutely infected, chronically infected, and recovered individuals within the proposed infection system under varying values of the parameter ς, specifically, ς = 1.00, 0.95, 0.90, 0.85. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 19 of 27 Time in days 0 10 20 30 40 50 60 O n e d o s e v a c c in a te d 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 ζ= 0.5 ζ=0.6 ζ= 0.7 ζ = 0.8 (a) Time in days 0 10 20 30 40 50 60 T w o d o s e v a c c in a te d 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 ζ= 0.5 ζ=0.6 ζ= 0.7 ζ = 0.8 (b) Time in days 0 10 20 30 40 50 60 A c u te i n fe c ti o u s 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ζ= 0.5 ζ=0.6 ζ= 0.7 ζ = 0.8 (c) Time in days 0 10 20 30 40 50 60 C h ro n ic i n fe c ti o u s 0.2 0.3 0.4 0.5 0.6 0.7 0.8 ζ= 0.5 ζ=0.6 ζ= 0.7 ζ = 0.8 (d) Time in days 0 10 20 30 40 50 60 R e c o v e re d i n d iv id u a ls 0.2 0.25 0.3 0.35 0.4 0.45 0.5 ζ= 0.5 ζ=0.6 ζ= 0.7 ζ = 0.8 (e) Figure 2: Visualization of the solution trajectories of vaccinated, acutely infected, chronically infected, and recovered individuals within the proposed infection system under varying values of the parameter ς, specifically, ς = 0.5, 0.6, 0.7, 0.8. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 20 of 27 Time in days 0 10 20 30 40 50 60 O n e d o s e v a c c in a te d 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ϕ1= 0.25 ϕ1=0.35 ϕ1= 0.45 ϕ1 = 0.55 (a) Time in days 0 10 20 30 40 50 60 T w o d o s e v a c c in a te d 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 ϕ1= 0.25 ϕ1=0.35 ϕ1= 0.45 ϕ1 = 0.55 (b) Time in days 0 10 20 30 40 50 60 A c u te i n fe c ti o u s 0 0.2 0.4 0.6 0.8 1 1.2 ϕ1= 0.25 ϕ1=0.35 ϕ1= 0.45 ϕ1 = 0.55 (c) Time in days 0 10 20 30 40 50 60 C h ro n ic i n fe c ti o u s 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 ϕ1= 0.25 ϕ1=0.35 ϕ1= 0.45 ϕ1 = 0.55 (d) Time in days 0 10 20 30 40 50 60 R e c o v e re d i n d iv id u a ls 0.43 0.44 0.45 0.46 0.47 0.48 0.49 ϕ1= 0.25 ϕ1=0.35 ϕ1= 0.45 ϕ1 = 0.55 (e) Figure 3: Graphical analysis of vaccinated, acutely infected, chronically infected, and recovered individu- als within the proposed infection model, considering varying values of the parameter φ1, specifically φ1 = 0.25, 0.35, 0.45, 0.55. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 21 of 27 Time in days 0 10 20 30 40 50 60 O n e d o s e v a c c in a te d 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 ξ= 0.45 ξ=0.55 ξ= 0.65 ξ = 0.75 (a) Time in days 0 10 20 30 40 50 60 T w o d o s e v a c c in a te d 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 ξ= 0.45 ξ=0.55 ξ= 0.65 ξ = 0.75 (b) Time in days 0 10 20 30 40 50 60 A c u te i n fe c ti o u s 0.2 0.4 0.6 0.8 1 1.2 1.4 ξ= 0.45 ξ=0.55 ξ= 0.65 ξ = 0.75 (c) Time in days 0 10 20 30 40 50 60 C h ro n ic i n fe c ti o u s 0.2 0.3 0.4 0.5 0.6 0.7 0.8 ξ= 0.45 ξ=0.55 ξ= 0.65 ξ = 0.75 (d) Time in days 0 10 20 30 40 50 60 R e c o v e re d i n d iv id u a ls 0.455 0.46 0.465 0.47 0.475 0.48 0.485 0.49 ξ= 0.45 ξ=0.55 ξ= 0.65 ξ = 0.75 (e) Figure 4: Graphical analysis of vaccinated, acutely infected, chronically infected, and recovered individu- als within the proposed infection model, considering varying values of the parameter ξ, specifically ξ = 0.45, 0.55, 0.65, 0.75. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 22 of 27 Time in days 0 10 20 30 40 50 60 A c u te i n fe c ti o u s 0 0.2 0.4 0.6 0.8 1 1.2 δ= 0.126 δ=0.226 δ= 0.326 δ = 0.426 (a) Time in days 0 10 20 30 40 50 60 C h ro n ic i n fe c ti o u s 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 δ= 0.126 δ=0.226 δ= 0.326 δ = 0.426 (b) Time in days 0 10 20 30 40 50 60 R e c o v e re d i n d iv id u a ls 0.43 0.44 0.45 0.46 0.47 0.48 0.49 δ= 0.126 δ=0.226 δ= 0.326 δ = 0.426 (c) Figure 5: Plotting the dynamic response of the solution trajectories of the proposed system as a function of the system’s input parameter δ i.e., δ = 0.126, 0.226, 0.326, 0.426. 5. Numerical findings In this section, we will examine the dynamical behaviour of the recommended system to show the impact of input factors on the output of the system. These analysis will provide a roadmap for responding to and managing outbreaks effectively, benefiting public health on both local and global scales. In a progressively interconnected global landscape, the capacity to model and forecast disease dissemination is imperative for global preparedness. Epidemic models play a crucial role in guiding early detection and containment strategies, thereby mitigating the swift worldwide transmission of infectious diseases.In this context, specific values for input parameters and state variables were assumed to facilitate the sim- ulations. All computations were performed on an HP PC laptop equipped with an Intel(R) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 23 of 27 Core(TM) i5-7200U CPU @ 2.50 GHz (2.71 GHz) and 8 GB of RAM. The simulations and numerical analysis were carried out using MATLAB (R2012a) software. In the first simulation presented in Figure 1 and Figure 2, we have show the effect of fractional order on the vaccinated, acute, chronic and recovered individuals of the system. In Figure 1, the value of zeta is considered to be ς = 1.00, 0.95, 0.90 and 0.85 while in Figure 2, we assume the value of fractional order ς = 0.8, 0.7, 0.6 and 0.5. In these Figures, we noticed that decreasing the value of fractional parameter can decrease the infection level in the society. This implies that this input factor is fruitful and can control the infection level in the society. In the second simulation presented in Figure 3, we have shown the impact of vaccination on the dynamics of the infection. We assumed the values of φ1 to be 0.25, 0.35, 0.45 and 0.55. It can be seen that vaccination is also an important factor and can significantly contribute to reduce the burden of the infection in the society. In Figure 4, we illustrated the solution pathways of vaccinated, acute, chronic and recovered individuals with the variation of the ξ. In this simulation, the values of ξ are assumed to be 0.45, 0.55, 0.65 and 0.75. We observed that the increasing value of this parameter increases the infection level of the hepatitis B which means that ξ is critical parameter and contribute to the spread of HB. In the last simulation illustrated in Figure 5, the impact of δ has been conceptualized on the vaccinated, acute, chronic and recovered individuals of the recommend fractional system. In our increasingly interconnected world, accurately modeling and predicting disease transmission is essential for global preparedness. Epidemic models enable early detection and guide containment strategies, playing a crucial role in preventing the rapid spread of infectious diseases. Understanding epidemic dynamics not only improves the communica- tion of the importance of preventive measures to the public but also promotes adherence to recommendations such as vaccination, treatment, quarantine, and other preventive ac- tions, ultimately helping to mitigate the impact of outbreaks. 6. Conclusion Hepatitis B infection constitutes a severe and potentially fatal viral disease, necessi- tating a thorough investigation into its transmission dynamics. This research has inves- tigated the dynamics of hepatitis B by introducing a two-shot vaccination series and has employed the fractional Caputo-Fabrizio operator to model the system. The fractional Caputo-Fabrizio model has been derived, and key findings related to this fractional model have been presented. Notably, the basic characteristics of the hepatitis B disease model have been determined, with the threshold parameter R0 having been determined using the next-generation matrix method. Stability analysis of infection-free equilibrium has been examined, and the uniqueness and existence of the solution for the hepatitis B dis- ease system have been investigated. A novel numerical approach has been introduced to explore the system of hepatitis B disease. The use of different orders of derivatives has re- vealed that the hypothesized model with the Caputo-Fabrizio derivative has yielded more accurate and efficient graphical results and analyses compared to integer-order deriva- tives. The graphical representations have also highlighted the impact and significance of I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 24 of 27 various input parameters and have visualized the most important factors of the system. The findings of this study are expected to provide critical insights for policymakers and healthcare authorities working to control the spread of hepatitis B within affected pop- ulations. Additionally, future extensions of this model may incorporate time-dependent control strategies, grounded in optimal control theory, to further enhance the effectiveness of disease management and intervention efforts. Declarations Competing interests: Authors declare no conflict of interest. Availability of data and materials: Data will be provided on request to the corresponding author. Acknowledgments and Funding Aziz Khan and Thabet Abdeljawad would like to thank Prince Sultan University for paying the APC and support through TAS research lab. References [1] Sneha Agarwal and Lakshmi Narayan Mishra. Attributes of residual neural networks for modeling fractional differential equations. Heliyon, 10(19), 2024. [2] Imtiaz Ahmad, Ihteram Ali, Rashid Jan, Sahar Ahmed Idris, and Mohamed Mousa. Solutions of a three-dimensional multi-term fractional anomalous solute transport model for contamination in groundwater. Plos one, 18(12):e0294348, 2023. [3] Imtiaz Ahmad, Asmidar Abu Bakar, Rashid Jan, and Salman Yussof. Dynamic behaviors of a modified computer virus model: Insights into parameters and network attributes. Alexandria Engineering Journal, 103:266–277, 2024. [4] Imtiaz Ahmad, Siraj-ul-Islam, Mehnaz, and Sakhi Zaman. Local meshless differential quadrature collocation method for time-fractional PDEs. Discrete & Continuous Dynamical Systems-Series S, 13(10), 2020. [5] Mustafa Ahmed and Shireen Jawad. The role of antibiotics and probiotics supple- ments on the stability of gut flora bacteria interactions. Commun. Math. Biol. Neu- rosci., 2023:Article–ID, 2023. [6] Saim Ahmed, Ahmad Taher Azar, Mahmoud Abdel-Aty, Hasib Khan, and Jehad Alzabut. A nonlinear system of hybrid fractional differential equations with applica- tion to fixed time sliding mode control for leukemia therapy. Ain Shams Engineering Journal, 15(4):102566, 2024. [7] Ahmed Ali, Shireen Jawad, Ali Hasan Ali, and Matthias Winter. Stability analysis for the phytoplankton-zooplankton model with depletion of dissolved oxygen and strong Allee effects. Results in Engineering, 22:102190, 2024. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 25 of 27 [8] Najat Almutairi and Sayed Saber. On chaos control of nonlinear fractional Newton-Leipnik system via fractional Caputo-Fabrizio derivatives. Scientific Reports, 13(1):22726, 2023. [9] Najat Almutairi and Sayed Saber. Existence of chaos and the approximate solution of the lorenz–lü–Chen system with the Caputo fractional operator. AIP Advances, 14(1), 2024. [10] Abdulrahman Obaid Alshammari, Imtiaz Ahmad, Rashid Jan, and Sahar Ahmed Idris. Fractional-calculus analysis of the dynamics of CD4+ T cells and human im- munodeficiency viruses. The European Physical Journal Special Topics, pages 1–13, 2024. [11] RM Anderson, GF Medley, and DJ Nokes. Preliminary analyses of the predicted impacts of various vaccination strategies on the transmission of hepatitis B virus. The control of Hepatitis B: The role of prevention in Adolescence, 95130, 1992. [12] Abdon Atangana and Kolade M. Owolabi. New numerical approach for fractional differential equations. Mathematical Modelling of Natural Phenomena, 13:3, 2018. [13] M.A. Belay, O.J. Abonyo, and D.M. Theuri. Mathematical model of hepatitis b dis- ease with optimal control and cost-effectiveness analysis. Computational and Mathe- matical Methods in Medicine, 2023(1):5215494, 2023. [14] Imtiyaz Ahmad Bhat and Lakshmi Narayan Mishra. A comparative study of dis- cretization techniques for augmented Urysohn type nonlinear functional volterra in- tegral equations and their convergence analysis. Applied Mathematics and Computa- tion, 470:128555, 2024. [15] Imtiyaz Ahmad Bhat, Lakshmi Narayan Mishra, Vishnu Narayan Mishra, Mahmoud Abdel-Aty, and Montasir Qasymeh. A comprehensive analysis for weakly singu- lar nonlinear functional volterra integral equations using discretization techniques. Alexandria Engineering Journal, 104:564–575, 2024. [16] Michele Caputo and Mauro Fabrizio. A new definition of fractional derivative without singular kernel. Progress in Fractional Differentiation & Applications, 1(2):73–85, 2015. [17] Stanca M Ciupe, Ruy M Ribeiro, and Alan S Perelson. Antibody responses during hepatitis B viral infection. PLoS computational biology, 10(7):e1003730, 2014. [18] Anwarud Din and Muhammad Zainul Abidin. Analysis of fractional-order vaccinated hepatitis-B epidemic model with Mittag-Leffler kernels. Mathematical Modelling and Numerical Simulation with Applications, 2(2):59–72, 2022. [19] W John Edmunds, GF Medley, D James Nokes, CJ O’callaghan, HC Whittle, and Andrew James Hall. Epidemiological patterns of hepatitis B virus (HBV) in highly endemic areasr. Epidemiology & Infection, 117(2):313–325, 1996. [20] WJ Edmunds, GF Medley, and DJ Nokes. The transmission dynamics and control of hepatitis B virus in the Gambia. Statistics in medicine, 15(20):2215–2233, 1996. [21] WJ Edmunds, GF Medley, DJ Nokes, AJ Hall, and HC Whittle. The influence of age on the development of the hepatitis B carrier state. Proceedings of the Royal Society of London. Series B: Biological Sciences, 253(1337):197–201, 1993. [22] Engida Endriyas Endashaw and Temesgen Tibebu Mekonnen. Modeling the effect I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 26 of 27 of vaccination and treatment on the transmission dynamics of hepatitis B virus and HIV/AIDS coinfection. Journal of Applied Mathematics, 2022(1):5246762, 2022. [23] Stephen A Gourley, Yang Kuang, and John D Nagy. Dynamics of a delay differen- tial equation model of hepatitis B virus infection. Journal of Biological Dynamics, 2(2):140–153, 2008. [24] J.H. He. A tutorial review on fractal spacetime and fractional calculus. International Journal of Theoretical Physics, 53:3698–3718, 2014. [25] Sarah Hews, Steffen Eikenberry, John D Nagy, and Yang Kuang. Rich dynamics of a hepatitis B viral infection model with logistic hepatocyte growth. Journal of Mathematical Biology, 60:573–590, 2010. [26] Hasib Khan, Jehad Alzabut, JF Gómez-Aguilar, and Abdulwasea Alkhazan. Essential criteria for existence of solution of a modified-ABC fractional order smoking model. Ain Shams Engineering Journal, 15(5):102646, 2024. [27] Hasib Khan, Altaf Hussain Rajpar, Jehad Alzabut, Muhammad Aslam, Sina Etemad, and Shahram Rezapour. On a fractal–fractional-based modeling for influenza and its analytical results. Qualitative Theory of Dynamical Systems, 23(2):70, 2024. [28] Jorge Losada and Juan J. Nieto. A new definition of fractional derivative without singular kernel. Progr. Fract. Differ. Appl, 1(2):87–92, 2015. [29] Manoochehr Makvandi. Update on occult hepatitis B virus infection. World journal of gastroenterology, 22(39):8720, 2016. [30] Graham F Medley, Nathan A Lindop, W John Edmunds, and D James Nokes. Hepatitis-B virus endemicity: heterogeneity, catastrophic dynamics and control. Na- ture medicine, 7(5):619–624, 2001. [31] Graham F Medley, Nathan A Lindop, W John Edmunds, and D James Nokes. Hepatitis-B virus endemicity: heterogeneity, catastrophic dynamics and control. Na- ture medicine, 7(5):619–624, 2001. [32] Douglas K Owens, Karina W Davidson, Alex H Krist, Michael J Barry, Michael Cabana, Aaron B Caughey, Chyke A Doubeni, John W Epling, Alex R Kemper, Martha Kubik, et al. Screening for hepatitis B virus infection in pregnant women: US Preventive Services Task Force reaffirmation recommendation statement. Jama, 322(4):349–354, 2019. [33] Monica Preboth. PHS guidelines for management of occupational exposure to HBV, HCV and HIV: management of occupational blood exposures. American Family Physician, 64(12):2012–2014, 2001. [34] Sayed Saber. Control of Chaos in the Burke-Shaw system of fractal-fractional order in the sense of Caputo-Fabrizio. Journal of Applied Mathematics and Computational Mechanics, 23(1):83–96, 2024. [35] Rafel Ibrahim Salih, Shireen Jawad, Kaushik Dehingia, and Anusmita Das. The effect of a psychological scare on the dynamics of the tumor-immune interaction with optimal control strategy. An International Journal of Optimization and Control: Theories & Applications (IJOCTA), 14(3):276–293, 2024. [36] Sarah Schillie, Trudy V Murphy, Mark Sawyer, Kathleen Ly, Elizabeth Hughes, Ruth Jiles, Marie A de Perio, Meredith Reilly, Kathy Byrd, John W Ward, et al. CDC I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5509 27 of 27 guidance for evaluating health-care personnel for hepatitis B virus protection and for administering postexposure management. MMWR Recomm Rep, 62(10):1–19, 2013. [37] Siraj-ul-Islam and Imtiaz Ahmad. Local meshless method for PDEs arising from models of wound healing. Applied Mathematical Modelling, 48:688–710, 2017. [38] Pauline Van den Driessche and James Watmough. Reproduction numbers and sub- threshold endemic equilibria for compartmental models of disease transmission. Math- ematical Biosciences, 180(1-2):29–48, 2002. [39] F. Wang, J. Zhang, I. Ahmad, A. Farooq, and H. Ahmad. A novel meshfree strategy for a viscous wave equation with variable coefficients. Frontiers in Physics, 9:701512, 2021. [40] Steven Wiersma. Hepatitis B virus: preventing liver disease with the first vaccine against cancer. Protect. Against Cancer-Causing Infect, 6, 2010. [41] JR Williams, DJ Nokes, GF Medley, and RM Anderson. The transmission dynamics of hepatitis B in the UK: a mathematical model for evaluating costs and effectiveness of immunization programmes. Epidemiology & Infection, 116(1):71–89, 1996.