EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6670 ISSN 1307-5543 – ejpam.com Published by New York Business Global Numerical Assessment of the Hepatitis B Virus Transmission Model Using a Non-standard Finite Difference Scheme and a Feedforward Neural Network Tahir Khan1, II Hyo Jung1,∗, Gul Zaman2, Ilyas Khan3 1 Department of Mathematics/Institute of Mathematical Science, Pusan National University, Busan, Republic of Korea 2 Department of Mathematics, University of Malakand Chakdara, Dir (L), Pakhtunkhwa, Pakistan 3 Faculty of Mathematics and Statistics, Ton Duc Thang University, Ho Chi Minh City, Vietnam Abstract. Hepatitis B virus (HBV) infection is one of the leading causes of death and is a conta- gious disease that produces chronic liver infection. The infection of hepatitis B has a complex nature involving multiple and long infectious periods, behavioral changes, and healthcare limitations. The multiple infection phases, immune system, behavioral changes, and healthcare saturation have a significant impact on the dynamics of hepatitis B virus transmission. Since HBV can persist in the body for years, the number of infectious individuals accumulates, and people take precautionary measures as awareness increases. Considering the multiple stages of the disease and the saturation level, we propose an innovative hybrid approach combining a mathematical model with a saturated incidence rate and a forward neural network to represent the transmission dynamics of the hepatitis B virus. First, we prove the biological and mathematical feasibility to show that the model under consideration is well-posed. We also investigate the dynamics of the model using linear stability analysis to derive the stability conditions. In addition, we perform the numerical assessment of the model using a novel hybrid approach of NSFD–FFNN framework to show the accuracy and large-scale numerical simulations. 2020 Mathematics Subject Classifications: 92D30, 65L12, 65M06, 68T07 Key Words and Phrases: Stability analysis, non-standard finite difference scheme, feed-forward neural network, numerical simulations 1. Introduction Living things, particularly humans, are constantly at risk from infectious diseases. Every person may experience an infectious disease once or multiple times during their ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6670 Email addresses: tahirmaths200014@gmail.com (T. Khan), ilhjung@pusan.ac.kr (I. H. Jung), gzaman@uom.edu.pk (G. Zaman), ilyaskhanqau@yahoo.com (I. Khan) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 2 of 22 lifetime. The infection could be acute or chronic, i.e., short-term or long-term, but cause mortality. The second leading cause of death and disability in human beings and other living organisms is infectious diseases around the globe, according to the WHO. In the current century, millions of people suffer due to the consequences of infectious diseases. Among these are tuberculosis, diarrheal diseases, malaria, hepatitis, HIV/AIDS, etc. The liver is one of the important organs in the human body, located in the upper right-hand part of the abdomen. Hepatitis B is one of the severe types of hepatitis that produces inflammation of the liver and leads to cirrhosis [1]. Hepatitis B is transmitted from one individual to another by direct contact or indirect contact [2]. Hepatitis B virus (HBV) infection progresses through multiple phases, most notably the acute and chronic stages, which are critical in understanding its clinical and epidemiological impact. During the acute phase, individuals often remain asymptomatic and may recover spontaneously with- out the need for antiviral treatment. However, a subset of infected individuals fail to clear the virus, leading to the development of a chronic infection, which may persist for life [3]. Globally, an estimated 300 million people are chronic carriers of HBV, facing signifi- cantly increased risks of HBV-related mortality. Known as a silent killer, HBV frequently progresses unnoticed, particularly in its chronic form. The acute phase typically lasts up to six months, during which the immune system may successfully eliminate the virus from the host. In contrast, the chronic phase is more severe and symptomatic, with the virus persisting in the body. This long-term infection can lead to serious complications, including liver cirrhosis and liver failure [4, 5]. While antiviral treatment is generally unnecessary during the acute phase, where rest, hydration, and proper nutrition are rec- ommended—immediate medical intervention becomes essential if the infection progresses to the chronic stage. Mathematical modeling refers to the process of representing real-world phenomena us- ing mathematical expressions, and it is widely applied across various disciplines, including physics, computer science, and biology [6–9]. In the context of epidemiology, mathemat- ical models play a crucial role in understanding the transmission dynamics of infectious diseases, evaluating control strategies, and predicting both short and long-term outcomes [10, 11]. These models consist of differential equations, which describe a framework to show the temporal evolution of disease spread (see, e.g., [12, 13]). However, analyzing the mathematical properties of such models and finding their exact or approximate so- lutions is challenging, but one must also look for alternative approaches while dealing with their solutions. To address this, alternative computational techniques, particularly neural network-based methods, have been increasingly employed to obtain accurate and efficient solutions. Given that many real-world systems are inherently governed by differ- ential equations, solving them using neural networks is an emerging and rapidly growing field of research [14]. In recent years, significant progress has been made, highlighting the promising interplay between neural networks and differential equation solvers [15]. Neural networks are computational models inspired by the architecture and function- ing of the human brain. They consist of interconnected artificial neurons organized into multiple layers, typically including an input layer, one or more hidden layers, and an output layer [16]. Input data passes through these layers, where mathematical operations T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 3 of 22 produce an output. Each artificial neuron processes incoming signals by applying weighted sums to the inputs, followed by the application of an activation (or merit) function, which introduces nonlinearity and determines whether the signal passes to the next layer. This layered structure enables the network to learn and model complex, nonlinear relationships in data. Among neural network architectures, the feed-forward neural network (FFNN) is particularly significant. In FFNNs, data flows in one direction from input to output, with- out feedback loops. They are widely used in the analysis of epidemiological models, where they approximate nonlinear dynamics involved in the transmission of infectious diseases. By learning from data, FFNNs offer a powerful, data-driven approach for forecasting dis- ease trends and enhancing the understanding of epidemic processes with high predictive accuracy. With the aid of nonlinear activation functions and multiple layers of neurons, FFNNs can approximate the disease transmission mechanisms to help estimate the key epi- demic parameters (infection rate, recovery rate, and transmission probability). Moreover, FFNNs are highly adaptable and can integrate real-world datasets, including time series infection data, demographic information, and vaccination, making them more effective in forecasting outbreaks in uncertain environments. When applied to epidemiological mod- els, FFNNs improve predictive accuracy, decision-making, intervention strategies, policy implementation, assisting public health authorities in resource allocation, and early warn- ing for outbreaks. Various neural network procedures are utilized to analyze the models as reported in [17–19]. However, the proposed epidemiological problem has never been considered or solved with the help of the proposed approaches. The main objective of the proposed work is to introduce a model representing the dynamics of the hepatitis B virus with reasonable assumptions and to analyze it with the aid of dynamical system theory, non-standard finite difference (NSFD) scheme, and feed- forward neural network (FFNN) framework. Precisely, we have the following contributions. 1. Since the multiple phases of hepatitis B have a significant impact on the transmission dynamics of the disease, we will consider the different infectious phases (latent, acute, and chronic) of the disease while formulating the proposed model. 2. Incorporating the saturated incidence rate is also an important consideration of the proposed work because the bilinear incidence only assumes that the incidence rate is proportional to the amount of susceptible and infectious population; however, this is not the case in hepatitis B virus transmission. Besides this one, the saturated incidence assumes the fact that transmission may plateau as the ratio of the infectious population reaches a certain threshold. Because the virus persists in individuals for a long time, the ratio of infected individuals may not be directly proportional to the number of susceptible individuals. 3. The framework of NSFD and FFNN will be utilized to present the accurate dynam- ics and to show the close match between predicted and true values, as well as to investigate the robustness of the proposed framework while modeling the dynamics of HBV. The paper is organized as follows. The detailed formulation of the model is provided T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 4 of 22 in section 2. We also discuss the biological and mathematical feasibility of the model in subsection 2.1. We linearize the proposed model and find the reproductive number in subsections 2.2 and 2.3 respectively. The dynamical properties of the model take place in subsections 2.4 and 2.5. We provide the framework of the NSFD and FFNN in section 3. We conclude our work in section 4. 2. Mathematical model formulation In this section, we will present an epidemic model for the transmission of HBV. Keeping the complex nature and multiple phases of HBV according to the disease characteristics, we assume three infected population groups of latently, acutely and chronically infected populations to acknowledge the biological importance of the latent, acute and chronic stage in HBV infection. Moreover, the infection has long infectious periods and can persist in the body for years, and as awareness increases, the people take precautionary measures, so the number of infectious individuals accumulates, therefore, we use saturated inci- dence rate βSI/(1 + γI), to better reflects the transmission scenarios. The susceptible individuals may transition to latent, acute and then chronically infected group of popu- lation after successful interaction with the infected individual. Since, there is no need of proper treatment in acute stage and those who recover naturally will transit to recover compartment, however, those who produce complication will leads to chronically infected population group. These infected individuals leave their compartments due death or full recovery. Parentally infected individuals will directly transit to chronically infected pop- ulation. Vaccine of hepatitis B are also very effective and provide immunity, so those who vaccinated successfully will transit to vaccinated compartment. Thus, the total population groups are distributed in six epidemiological subclasses i.e. s, l, a, c, r and v respectively, represent the susceptible, latent, acute infected, chronic carrier infectious, recovered and vaccinated individuals, whose the schematic process of the disease propagation among dif- ferent epidemiological groups of population is illustrated by a flowchart as given in Figure 1. Eventually, the dynamics of these compartments are governed by the following system of nonlinear differential equations: ds(t) dt = gη(1− φc(t)) + πv(t)− αs(t)a(t) 1 + βc(t) − λαs(t)c(t) 1 + βc(t) − (µ0 + λ3)s(t), dl(t) dt = αs(t)a(t) 1 + βc(t) + λαs(t)c(t) 1 + βc(t) − (µ0 + θ)l(t), da(t) dt = θl(t)− (µ0 + λ1)a(t), dc(t) dt = gηφc(t) + qλ1a(t)− (µ0 + µ1 + λ2)c(t), dr(t) dt = λ2c(t) + (1− q)λ1a(t)− µ0r(t), dv(t) dt = g(1− η) + λ3s(t)− (µ0 + π)v(t), (1) T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 5 of 22 Figure 1: The plot represent the schematic process of the HBV propagation among different epidemiological groups of population. along with non-negative initial population sizes s(0), v(0) > 0, l(0), a(0), c(0), r(0) ≥ 0. (2) In the model (1), g represents the birth rate while the rate of birth without successful vaccination is represented by η. The quantity of maternally distorted individuals is de- noted by φ, and the rate of waning vaccine-induced immunity is π. We symbolize the transmission rate from susceptible to infected, denoted by the parameters α and λ, as the reduced transmission rate. µ0 represents rate of natural death. The rate of vaccination is λ3, and the amount at which the latent individuals lead to acute class is represented by θ. In contrast, the portion of individuals who transition from the acute population to the chronic carrier is represented by λ1. We also denote the rate at which the individuals move from chronic carrier to the immune class by λ2, and the parameter µ1 is the death that occurs from hepatitis B. Moreover, we assume that q is the average rate of those individuals who become unsuccessful in recovering from hepatitis B in the primary stage of hepatitis B (acute) and move to the chronic stage. Having presented the detailed model formulation, we proceed to analyze the proposed epidemic model by establishing its well-posedness, ensuring that it is both biologically and mathematically meaningful. 2.1. Well-posedness We discuss the well-posedness of the proposed system to show that the epidemic model is feasible in both senses, biologically and mathematically. For this, we assume that the total population is represented by n(t) with n = v+l+a+s+r+c. The model (1), together T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 6 of 22 with initial conditions (2) satisfies that n(t) ≥ 0 implies that n(t) (total population) is bounded and positive for t > 0, thus the time derivative of n(t) gives dn(t) dt = g − µ0n(t)− µ1c(t), (3) which shows that as t → ∞, n(t) ≤ g µ0 , hence Ω = { (s(t), l(t), a(t), c(t), r(t), v(t)) ∈ R6 +, n(t) ≤ g µ0 } . (4) is the feasible region. To show the existence analysis, we can express the proposed epidemic problem as initial value problem as follows: H′ (t) = F(y), y(0) = y0, (5) where y = (s, l, a, c, r, v)t and F = (f1, f2, . . . , f6) t. To proceed further, we recall the definition of Lipschitz continuity as follows. Definition 1. [20] The function F(y) is said to be Lipschitz continuous if there exists a constant M > 0 such that ∥F(y1)−F(y2)∥ ≤ M∥y1 − y2∥. Thus, regarding existence of solution, we prove the following proposition. Proposition 2.1. There exists a unique solution of the system (5) analogous to the proposed model (1) and (2). Proof. To show that the model posses a unique solution, it is sufficient to prove that the function F(y) satisfies the Lipschitz condition provided in the above definition. Keeping in mind the set Ω and n(t) ≤ g/µ0, which implies that all the state variables for i = 1, 2, si, li, ai, ci, ri and vi are bounded. One can proceed as follows: ∥f1(y1)− f1(y2)∥ = ∥∥∥∥gηφ(c2 − c1) + π(v1 − v2)− α ( s1a1 1 + βc1 − s2a2 1 + βc2 ) −λα ( s1c1 1 + βc1 − s2c2 1 + βc2 ) − (µ0 + λ3)(s1 − s2) ∥∥∥∥ , which implies that ∥f1(y1)− f1(y2)∥ ≤ gηφ∥c1 − c2∥+ π∥v1 − v2∥+ (µ0 + λ3)∥s1 − s2∥ + α ∥∥∥∥ s1a1 1 + βc1 − s2a2 1 + βc2 ∥∥∥∥+ λα ∥∥∥∥ s1c1 1 + βc1 − s2c2 1 + βc2 ∥∥∥∥ . (6) To simplify the above inequality, we use some algebraic manipulation as follows s1a1 1 + βc1 − s2a2 1 + βc2 = s1a1 − s2a2 1 + βc1 − s2a2 ( 1 1 + βc1 − 1 1 + βc2 ) , s1c1 1 + βc1 − s2c2 1 + βc2 = s1c1 − s2c2 1 + βc1 − s2c2 ( 1 1 + βc1 − 1 1 + βc2 ) . T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 7 of 22 Also ∥s1a1 − s2a2∥ ≤ ∥s1∥∥a1 − a2∥+ ∥a2∥∥s1 − s2∥, ∥s1c1 − s2c2∥ ≤ ∥s1∥∥c1 − c2∥+ ∥c2∥∥s1 − s2∥, and ∥∥∥∥ 1 1 + βc1 − 1 1 + βc2 ∥∥∥∥ ≤ β∥c1 − c2∥ ∥(1 + βc1)(1 + βc2)∥ ≤ β∥c1 − c2∥. Using the above algebraic manipulations, the equation (6) leads to the following inequality ∥f1(y1)− f1(y2)∥ ≤ { µ0 + λ3 + α∥a2∥+ λα∥c2∥ } ∥s1 − s2∥+ α∥s1∥∥a1 − a2∥ + { gηφ+ αβ∥s2∥∥a2∥+ αβλ∥s2∥∥c2∥+ αλ∥s1∥ } ∥c1 − c2∥+ π∥v1 − v2∥. Let us assume that the M1, M2 and M3 are any positive constants, then the above equation takes the form ∥f1(y1)− f1(y2)∥ ≤ M1∥s1 − s2∥+M2∥a1 − a2∥+M3∥c1 − c2∥+ π∥v1 − v2∥. Let M = max{M1,M2,M3, π}, the last inequality can be re-written as ∥f1(y1)− f1(y2)∥ ≤ M ( ∥s1− s2∥+ ∥a1− a2∥+ ∥c1− c2∥+ ∥v1− v2∥ ) ≤ M∥y1− y2∥. (7) The second equation of the model (1), implies that ∥f2(y1)− f2(y2)∥ ≤ { α∥a2∥+ λα∥c2∥ } ∥s1 − s2∥+ (µ0 + θ)∥l1 − l2∥+ α∥s1∥∥a1 − a2∥ + { αβ∥s2∥∥a2∥+ αβλ∥s2∥∥c2∥+ αλ∥s1∥ } ∥c1 − c2∥. If N1, N2, N3 and N4 are positive constants, and N = max{N1,N2,N3,N4}, then the above inequality gives ∥f2(y1)− f2(y2)∥ = N { ∥s1 − s2∥+ ∥l1 − l2∥+ ∥∥a1 − a2∥+ ∥c1 − c2∥ ≤ N∥y1 − y2∥ } . (8) From the third equation of the model (1), we obtain ∥f3(y1)− f3(y2)∥ ≤ θ∥l1 − l2∥+ (µ0 + λ1)∥a1 − a2∥ ≤ O∥y1 − y2∥, (9) where O = {θ, (µ0 + λ1)}. The fourth equation of the model (1) implies that ∥f4(y1)− f4(y2)∥ ≤ qλ1∥a1 − a2∥+ (gηφ+ µ0 + µ1 + λ2) ∥c1 − c2∥ ≤ P∥y1 − y2∥, (10) where P = {qλ1, (gηφ+ µ0 + µ1 + λ2)}. In a similar way, the last two equation of the proposed model gives ∥f5(y1)− f5(y2)∥ ≤ (1− q)λ1∥a1 − a2∥+ λ2∥c1 − c2∥+ µ0∥r1 − r2∥ ≤ Q∥y1 − y2∥, ∥f6(y1)− f6(y2)∥ ≤ λ3∥s1 − s2∥+ (µ0 + π)∥v1 − v2∥ ≤ S∥y1 − y2∥, (11) T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 8 of 22 where Q = {(1− q)λ1, λ2, µ0} and S = {λ3, (µ0 + π)}. Combining equation (7) to (11), we obtain ∥F(y1)−F(y2)∥ = H∥y1 − y2∥, where H = max{M,N ,O,P,Q,S}, ensure that the function F(y) is Lipschitz continuous, and hence the model posses a unique solution. Proposition 2.2. Let us consider that (s(t), l(t), a(t), c(t), r(t), v(t)) is the solution of the model with non-negative initial conditions, and the closed set Ω is the attracting under the flow described by equation (1) that is positive invariant. Proof. Let us assume a function L(t) given by L(t) = s(t) + l(t) + a(t) + c(t) + r(t) + v(t), (12) then the time differentiation with the implementation of equation (1) leads to dL(t) dt = g − µ0l(t)− µ1c(t), (13) which implies that dL(t) dt ≤ g − µ0l(t) ≤ 0. (14) It is clear from the equation (14), that dL(t) dt ≤ 0, which indicates that the set φ is positively invariant. Also, 0 ≤ l(t) ≤ l(0)e−µ0t + g µ0 (1− e−µ0t). Thus, as t → ∞, 0 ≤ L(t) ≤ g µ0 , which is sufficient to prove that Ω is an attracting set. The above proposition confirms that the proposed model is well-posed. We now proceed with the linearization of the model, computation of the threshold quantity, identification of the model equilibria, and its stability analysis, as outlined below. 2.2. Linearization of the model Using the methods of dynamical systems to linearize the model that is under consid- eration. To make our calculation simpler, we take the Jacobian matrix of the reduced system (without the 5th equation), which looks like J =  −x11 0 −x13 x14 π x21 −x22 x23 −x24 0 0 θ −x33 0 0 0 0 qλ1 −x44 0 λ3 0 0 0 −x55  , where x11 = αa(t) 1 + βc(t) + λαc(t) 1 + βc(t) + µ0 + λ3, x13 = αs(t) 1 + βc(t) = x23, x14 = βαs(t)a(t) 1 + βc(t)2 − λαs(t) 1 + βc(t)2 = −x24, x21 = x11 − (µ0 + λ3), x22 = µ0 + θ, x33 = µ0 + λ1, x44 = gηφ− µ0 − µ1 − λ2, x55 = µ0 + π. T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 9 of 22 2.3. Basic reproduction number In epidemiological models, the threshold quantity, called the disease’s basic repro- duction number, is an essential quantity which describes the predicted average rate of infections induced by a single infective after introducing into an entirely susceptible com- munity. To calculate this quantity for the proposed problem, we follow the methodology introduced by Watmough and Driessche [21], therefore assume that y = (l(t), a(t), c(t)), then the model (1) yields dy dt = F − V. (15) In equation (15), the matrices are define as: F =  βs(t)a(t) 1+βc(t) + λαs(t)c(t) 1+βc(t) 0 0  , V =  (θ + µ0)l(t) (µ0 + λ1)a(t)− θl(t) (µ0 + µ1 + λ2)c(t)− qλ1a(t)− gηφc(t)  . Let F and V are the Jacobian of F and V at F0, then it becomes F = Jacobian of F =  0 αs0 λαs0 0 0 0 0 0 0  , V = Jacobian of V =  x11 0 0 −θ x22 0 0 −qλ1 x33  , where x11 = θ + µ0, x22 = µ0 + λ1 and x33 = µ0 + µ1 + λ2 − gηφ. ˜Thus, the basic reproductive number r0 is the spectral radius of K = FV−1 (r0 = ρ(FV−1)). So the basic reproductive quantity r0 of our suggested model (1) looks like r0 = r1 + r2, (16) where r1 = θαs0 (µ0 + θ)(µ0 + λ1) , r2 = θαs0λλ1q (µ0 + θ)(µ0 + λ1)(µ0 + µ1 + λ2 − gηφ) . 2.4. Equilibrium analysis To study the dynamics of the model that is under consideration (1), we first find the model equilibria. For this, we assume that the infection-free state of system (1) is T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 10 of 22 symbolized by F0 and F0 = (s0, 0, 0, 0, 0, v0), where s0 = g{π + ηµ0} µ0{π + µ0 + λ3} , v0 = g{µ0 − µ0η + λ3} µ0{π + λ3 + µ0} . (17) Similarly, to find the endemic states of the model, we assume P1 = (µ0+λ3), P2 = (µ0+θ), P3 = (µ0 + λ1), P4 = (µ0 + µ1 + λ2 − gηφ) and P5 = (µ0 + π) for the shake of simplicity, then the endemic state becomes F1 = (s∗, l∗, a∗, c∗, r∗, v∗), where s∗ = {P5gη + πg(1− η)} {P1P5 − πλ3} , l∗ = P2P3P4µ0{µ0 + π + λ3}(r0 − 1) + P2P3P4πλ3 P2θ{P4P5α+ P5λαqλ1 + P1P5βqλ1 − βqλ1πλ3} , a∗ = P2P4µ0{µ0 + π + λ3}(r0 − 1) + P2P3P4πλ3 P2{αP4P5 + P5λαqλ1 + P1P5βqλ1 − βqλ1πλ3} , r∗ = 1 µ0 {λ2c+ (1− q)λ1a}, c∗ = P2qλ1µ0{µ0 + π + λ3}(r0 − 1) + P2P3P4πqλ1λ3 P2{αP4P5 + P5 + P5λαqλ1 + P1P5βqλ1 − βqλ1πλ3} , v∗ = 1 P5 {g(1− η) + λ3s}. (18) Now we describe the asymptotic stabilities of the model that is under consideration in the upcoming sections. 2.5. Dynamical properties of the model We discuss the local analysis of the model to investigate the local asymptotic stability of the proposed problem at F0 and F1. To do this, we follow linearization and Routh- Hurwitz criteria. Thus, we state the following results. Theorem 1. The proposed model at F0 is locally asymptotic stable whenever r0 < 1 and gηφ > µ0 + µ1 + λ2. Proof. The characteristic equation of the matrix J at F0 takes the following form: (λ+ x44)(λ 4 + x1λ 3 + x2λ 2 + x3λ+ x4) = 0, (19) where x1 = x11 + x22 + x33 + x44 + x55, x2 = x11x22 + x11x33 + x11x55 + x22x55 + x33x55 + x22x33(1− r1), x3 = µ0(µ0 + π + λ3)(x22 + x33) + x22x33(x11 + x55)(1− r1), x4 = µ0(µ0 + π + λ3)x22x33(1− r1). Clearly, one eigenvalue λ1 has negative real part, if gηφ > µ0+µ1+λ2. For the remaining, if r0 < 1, we have 0 < rj < 1, j = 1, 2. So xi > 0, for i = 1, 2, 3, 4, and x1x2x3 > x23+x21x4 holds, which implies that the Routh-Hurwitz criterion for stability of the model holds. Hence, the proposed model is locally asymptotically stable at the F0, if r0 < 1. Following the same procedure, we can discuss the local dynamics of the model at the endemic state by the subsequent theorem. T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 11 of 22 Theorem 2. For r0 > 1, the model is stable locally asymptotically if the following condi- tions hold: 1. βαs∗ < (1 + βc∗)(µ0 + θ)(µ0 + λ1). 2. βαs∗ < (µ0 + θ){αa∗ + λαc∗ + (1 + βc∗)(µ0 + λ3)}. Proof. The characteristic equation of the matrix J at endemic equilibrium F1 looks like λ5 + d1λ 4 + d2λ 3 + d3λ 2 + d4λ+ d5 = 0, (20) where d1 =x11 + x22 + x33 + x44 + x55, d2 ={µ0 + π}x21 + µ0{µ0 + λ3 + π}+ x11x22 + x11x33 + x11x44 + x22x33 + x22x44 + x22x55 + x33x44 + x33x55 + x44x55 − θx22, d3 =θqλ1x24 + x11{x22x33 − θx23}+ x44{x11x22 − θx23}+ x11x33x44 + θx21x13 + {x22 + x33 + x44}{(µ0 + π) + µ0(µ0 + π + λ3)}+ x55{x22x33 − θx23} + x22x44x55 + x33x44x55 + x22x33x44 + x21, d4 =θx21x31x55 + {(µ0 + π)x21 + µ0(µ0 + π + λ3)}{x22x33 + x22x44 + x33x44} + θqλ1(x55 − x14) + x44x55(x22x33 − θx23) + x11x44{x22x33 − θx23}+ θx21x13x44 + θπλ3x23 + θqλ1x11x24 − θx11x23x55, d5 ={(µ0 + π)x21 + µ0(µ0 + π + λ3)}x22x33x44 + θλ3πx23x44θλ3πqλ1x24 + θqλ1x55{x11x24 − x14}+ θx13x24x44x55 − θx11x23x44x55. All di > 0, for i = 1, 2, 3, 4, 5 and H0 : { r0 > 1, d1d2d3 > d23 + d21d4, (d1d4 − d5) ( d1d2d3 − d23 − d21d4 ) > d5 ( d1d2 − d23 )2 + x1x 2 5 } , holds implies that H0 is satisfied if conditions 1 and 2 hold. Thus, it ensures that the Routh-Hurwitz stability criteria are satisfied, and we conclude that the endemic state F1 is locally asymptotically stable. Upon establishing the analytical results, we now proceed to numerically simulate the model in order to validate our theoretical findings, as detailed in the following section. 3. Numerical simulation In this section, we carry out numerical experiments to validate the analytical results and demonstrate the robustness of the proposed hybrid approach. We begin by discretizing the model using the nonstandard finite difference (NSFD) scheme, followed by a summary of the architecture of the proposed feed-forward neural network (FFNN). T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 12 of 22 3.1. Non-standard finite difference scheme and the architecture of the neural network In this section, we present the model discretization with the aid of NSFD and provide the architecture of the neural network that will be used for the learning process. It is obvious that the proposed model has six state variables i.e. s(t), l(t), a(t), c(t), r(t), and v(t), and g, η, π, α, β, λ, λ1, λ2, λ3, µ0, θ, q and µ1 are the parameters. To approximate the model, we use a non-standard finite difference scheme given by yn+1 = yn +∆t.h(yn), (21) for every y ∈ {s, l, a, c, r, v}, and h(yn) represents the right side of proposed model. Let ∆t is the discrete time step represent the approximation at t = n∆t, then sn ≈ s(n∆t), ln ≈ l(n∆t), an ≈ a(n∆t), cn ≈ c(n∆t), rn ≈ r(n∆t), vn ≈ v(n∆t). Using NSFD discretization, the proposed model takes the following form sn+1 = sn +∆t { gη(1− φcn + πvn)− αsnan 1 + βcn − λαsncn 1 + βcn − (µ0 + λ3)sn } , ln+1 = ln +∆t { αsnan 1 + βcn + λαsncn 1 + βcn − (µ0 + θ)ln } , an+1 = an +∆t {θln − (µ0 + λ1)an} , cn+1 = cn +∆t {gηφcn + qλ1an − (µ0 + µ1 + λ2)cn} , rn+1 = rn +∆t {λ2cn + (1− q)λ1an − µ0rn} , vn+1 = vn +∆t {g(1− η) + λ3sn − (µ0 + π)vn} . (22) Using the hypothetical values of the model parameters as: g = 0.7, η = 0.85, φ = 0.05, π = 0.08, α = 0.4, β = 0.15, λ = 0.35, λ1 = 0.12, λ2 = 0.06, λ3 = 0.03, µ0 = 0.02, θ = 0.15, q = 0.75, µ1 = 0.025. To approximate the proposed system numerically, we then use a standard feed-forward neural network (FFNN) with one hidden layer consisting of input layer t ∈ R, number of neurons nh = 10, and activation function (ReLU) σ(.). Let W (1) ∈ Rnh×1 represent the weights from input to the hidden layer and b(1) ∈ Rnh is the biases in the hidden layer, then the hidden layer out is given by h = σ ( W 1t+ b(1) ) ∈ Rnh . (23) Likewise, the output layer is a six-dimensional vector (s, l, a, c, r, v). Let W (2) ∈ R6×nh represent the weights from hidden to output layer, and b2 ∈ R6 represent the biases in output layer. Then the predicted output becomes ŷ(t) = W (2)h+ b(2) ∈ R6. (24) Combining all, we get the following expression ŷ(t) = W (2)hσ ( W (1)t + b(1) ) + b(2), (25) T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 13 of 22 where ŷ(t) = [ ŝ(t), l̂(t), â(t), ĉ(t), r̂(t), v̂(t) ]t and σ(.) is the activation function. Further, to minimize the Mean Squared Error (MSE), we train the network to minimize the error between the true and predicted values over the training set {(ti, yi)}Ni=1, which is given by MSE = 1 6N N∑ i=1 ∥yi − ŷ(ti)∥2, (26) where yi = [s(ti), l(ti), a(ti), c(ti), r(ti), v(ti)] t and ŷ(ti) is the network output. 3.2. Numerical experiment To perform the numerical experiment, we use the discretized system (22) to generate the data for the learning process, and then use a feed-forward neural network (FFNN) to forecast the dynamics of the proposed model in the long run. For each state (s, l, a, c, r, v) of the model, we present the dynamics as reported in Figures 2 to 7. Each figure represents the dynamics of the compartmental population of the proposed model with absolute error. The blue line demonstrates the true values of the model states over time using NSFD, and the red dashed line shows the predicted values of the model states over time using FFNN. Clearly, the results show how the true and approximated values are closer because the blue line follows the red-dashed line, which guarantees a better prediction. We also provided the associated absolute error graphs as shown by a magenta line for each state of the model, which shows the magnitude of error, describing how far off the prediction is (see Figures 2 to 7). To quantify the overall error prediction for each model state variable (s, l, a, c, r, and v), a Mean Square Error (MSE) bar graph is presented as shown by Figure 8. This provides the overall performance of the proposed FFNN on each state variable. On the x-axis, the indices of the variables s, l, a, c, r, and v are taken, while on the y-axis, the Mean Square Error (MSE) values are plotted, which clearly shows that the low MSE guarantees better FFNN performance. In addition, the validation checks, damping factor (µ), and Gradient interpretation are provided as shown in Figure 9. More precisely, this demonstrates the performance of the validation set against epochs. Every time the check occurs, if the performance fails to improve for a specific number of epochs, however, the training stops, if the performance on these validations set does not improve over these checks. Also, to show the training efficiency, the damping factor µ is visualized as shown in Figure 9. Similarly, to show the performance of the gradient with respect to the biases and weights, that is, how steep the slope is at the current point in the error surface, we interpret the gradient. In the end, the regression plots are provided for the FFNN training to present the correlation between the actual and predicted values for the data set and R-value as illustrated in Figure 9, which shows the performance of the trained FFNN. In every subplot of Figure 10, representing perfect agreement. The regression lines fitted to the data lie almost exactly on the identity line, with negligible intercepts, suggesting minimal prediction bias. Moreover, the close alignment verifies that the network has learned well the underlying mapping with high precision. The results highlight the accuracy and robustness of the T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 14 of 22 hybrid approach of NSFD–FFNN, demonstrating its effectiveness in capturing the complex nonlinear dynamics of the modeled system. 0 20 40 60 80 100 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 t s True Predicted 0 20 40 60 80 100 0 0.005 0.01 0.015 0.02 0.025 t A b s o lu te E rr o r Figure 2: The graphs represents the temporal dynamics of the state s and its associated absolute error. Overall, we observe that the feed-forward neural network (FFNN) learned from the data simulated by NSFD exhibits an excellent learning and generalization capacity. The validation checks describe minimal overfitting, and the trajectory of the damping factor (µ) shows adaptive learning and stable behavior, as well as effective convergence, as seen by the gradient descent pattern. The regression analysis shows that the clustering of points around the diagonal line indicates that the proposed network has efficiently captured the nonlinear dynamics of the proposed epidemiological model simulated by NSFD, which confirms the robustness of the underlying hybrid NSFD–FFNN framework. 0 20 40 60 80 100 0 0.5 1 1.5 2 2.5 3 3.5 4 t l True Predicted 0 20 40 60 80 100 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 t A b s o lu te E rr o r Figure 3: The time dynamics of the state l with its absolute error. T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 15 of 22 0 20 40 60 80 100 0 0.5 1 1.5 2 2.5 3 3.5 4 t a True Predicted 0 20 40 60 80 100 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 t A b s o lu te E rr o r Figure 4: This demonstrate the dynamics of state a and its absolute error. 0 20 40 60 80 100 0 0.5 1 1.5 2 2.5 3 3.5 4 t c True Predicted 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 x 10 −3 t A b s o lu te E rr o r Figure 5: The graphs visualizes the dynamics of the model state c with its associated absolute error. T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 16 of 22 0 20 40 60 80 100 0 2 4 6 8 10 12 14 t r True Predicted 0 20 40 60 80 100 0 0.5 1 1.5 2 2.5 x 10 −3 t A b s o lu te E rr o r Figure 6: Dynamics of the model state r with absolute error. 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 t v True Predicted 0 20 40 60 80 100 0 1 2 3 4 5 6 7 x 10 −3 t A b s o lu te E rr o r Figure 7: The dynamics of the model state v with its associated absolute error. T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 17 of 22 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Error Histogram with 20 Bins In st a n ce s Errors = Targets − Outputs − 0 .0 2 1 1 − 0 .0 1 9 3 2 − 0 .0 1 7 5 4 − 0 .0 1 5 7 6 − 0 .0 1 3 9 8 − 0 .0 1 2 2 − 0 .0 1 0 4 2 − 0 .0 0 8 6 4 − 0 .0 0 6 8 6 − 0 .0 0 5 0 8 − 0 .0 0 3 3 − 0 .0 0 1 5 2 0 .0 0 0 2 5 8 0 .0 0 2 0 3 8 0 .0 0 3 8 1 8 0 .0 0 5 5 9 8 0 .0 0 7 3 7 8 0 .0 0 9 1 5 9 0 .0 1 0 9 4 0 .0 1 2 7 2 Training Validation Test Zero Error Figure 8: The graphs show time dynamics of the compartmental population of the model that is under consid- eration at the disease free equilibrium state. T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 18 of 22 10 −10 10 0 10 10 g ra d ie n t Gradient = 8.8342e−06, at epoch 1000 10 −10 10 −5 10 0 m u Mu = 1e−06, at epoch 1000 0 200 400 600 800 1000 0 0.5 1 v a l fa il 1000 Epochs Validation Checks = 0, at epoch 1000 Figure 9: The graphs show time dynamics of the compartmental population of the model that is under consid- eration at the disease free equilibrium state. T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 19 of 22 2 4 6 8 10 12 2 4 6 8 10 12 Target O ut pu t ~ = 1* Ta rg et + 5 .9 e− 07 Training: R=1 Data Fit Y = T 2 4 6 8 10 12 2 4 6 8 10 12 Target O ut pu t ~ = 1* Ta rg et + 9 .5 e− 05 Validation: R=1 Data Fit Y = T 2 4 6 8 10 12 2 4 6 8 10 12 Target O ut pu t ~ = 1* Ta rg et + 3 e− 06 Test: R=1 Data Fit Y = T 2 4 6 8 10 12 2 4 6 8 10 12 Target O ut pu t ~ = 1* Ta rg et + 1 .5 e− 05 All: R=1 Data Fit Y = T Figure 10: The graphs show time dynamics of the compartmental population of the model that is under consideration at the disease free equilibrium state. 4. Conclusion In this paper, we investigated the transmission dynamics of hepatitis B virus (HBV) using the integration of an epidemic model FFNN. The model has been carefully con- structed to reflect the key biological features of HBV, including its various infection stages (latent, acute, and chronic), as well as a saturated incidence rate to capture the behavioral and healthcare constraints associated with HBV transmission. We established the well- posedness of the model to ensure biological and mathematical feasibility of the problem. We also derived the basic reproduction number, and studied a comprehensive stability analysis to assess the conditions under which the disease-free and endemic states are sta- ble. To validate the analytical findings and explore the model’s behavior, we employed a novel hybrid numerical framework combining the nonstandard finite difference (NSFD) scheme with a feed-forward neural network (FFNN). The numerical results demonstrated a high degree of accuracy, evidenced by low mean squared error and a strong agreement between the true trajectories produced by NSFD and predicted trajectories provided by T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 20 of 22 FFNN. This confirms the robustness and reliability of both the proposed model and the computational approach. The study not only provides valuable insights into the complex dynamics of HBV transmission but also provides evidence of the effectiveness of a hy- brid approach capturing the progression of infectious diseases. Moreover, the proposed approach is computationally stable, sound, and adaptable, which offers a promising tool for other epidemiological models. In future, we will extend the model to incorporate the stochastic perturbation to capturing disease uncertainty in heterogeneous environments, and the integrate it with machine learning with the help of feed-forward neural network. The FFNNs will be used to learn the underlying nonlinear dynamics governed by the associated stochastic differential equations (SDEs) and predict for the long run. Acknowledgements This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean Government (MSIT) (No. RS-2024-00342113 and 2022R1A5A1033624). References [1] Scott D Holmberg, Anil Suryaprasad, and John W Ward. Updated cdc recommen- dations for the management of hepatitis b virus infected health care providers and students. 2012. [2] Sonia Altizer, Andrew Dobson, Parviez Hosseini, Peter Hudson, Mercedes Pascual, and Pejman Rohani. Seasonality and the dynamics of infectious diseases. Ecology letters, 9(4):467–484, 2006. [3] Brian J McMahon. Epidemiology and natural history of hepatitis b. In Seminars in liver disease, volume 25, pages 3–8. Published in 2005 by Thieme Medical Publishers, Inc., 333 Seventh Avenue . . . , 2005. [4] Yukihiko Nakata and Toshikazu Kuniya. Global dynamics of a class of seirs epidemic models in a periodic environment. Journal of Mathematical Analysis and Applications, 363(1):230–237, 2010. [5] Marc Ringehan, Jane A McKeating, and Ulrike Protzer. Viral hepatitis and liver cancer. Philosophical Transactions of the Royal Society B: Biological Sciences, 372(1732):20160274, 2017. [6] Abubakar Mwasa and Jean M Tchuenche. Mathematical analysis of a cholera model with public health interventions. Biosystems, 105(3):190–200, 2011. [7] Tailei Zhang, Kai Wang, and Xueliang Zhang. Modeling and analyzing the transmis- sion dynamics of hbv epidemic in xinjiang, china. PloS one, 10(9):e0138765, 2015. [8] Kamil Shah, Jamal Shah, Ebenezer Bonyah, Tmader Alballa, Hamiden Abd El- Wahed Khalifa, Usman Khan, and Hameed Khan. Optimal control of covid-19 through strategic mathematical modeling: Incorporating harmonic mean incident rate and vaccination. AIP Advances, 14(9), 2024. T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 21 of 22 [9] Muhammad Farhan, Fahad Aljuaydi, Zahir Shah, Ebraheem Alzahrani, Ebenezer Bonyah, and Saeed Islam. A fractional modeling approach to a new hepatitis b model in light of asymptomatic carriers, vaccination and treatment. Scientific African, 24:e02127, 2024. [10] Lan Zou, Weinian Zhang, and Shigui Ruan. Modeling the transmission dynamics and control of hepatitis b virus in china. Journal of theoretical biology, 262(2):330–338, 2010. [11] Junjie Zhu, Misbah Ullah, Saif Ullah, Muhammad Bilal Riaz, Abdul Baseer Saqib, Atif M Alamri, and Salman A AlQahtani. A novel numerical solution of nonlinear stochastic model for the propagation of malicious codes in wireless sensor networks using a high order spectral collocation technique. Scientific Reports, 15(1):228, 2025. [12] Anwarud Din and Yongjin Li. Optimizing hiv/aids dynamics: stochastic control strategies with education and treatment. The European Physical Journal Plus, 139(9):812, 2024. [13] Wan Yang, Alicia Karspeck, and Jeffrey Shaman. Comparison of filtering methods for the modeling and retrospective forecasting of influenza epidemics. PLoS compu- tational biology, 10(4):e1003583, 2014. [14] Suthep Suantai, Zulqurnain Sabir, Muhammad Asif Zahoor Raja, and Watchara- porn Cholamjiak. A stochastic bayesian neural network for the mosquito dispersal mathematical system. Fractal and Fractional, 6(10):604, 2022. [15] Atifa Asghar, Mohsan Hassan, Zulqurnain Sabir, Shahid Ahmad Bhat, and Sharifah E Alhazmi. A design of computational stochastic framework for the mathematical severe acute respiratory syndrome coronavirus model. Biomedical Signal Processing and Control, 100:107049, 2025. [16] Muhammad Farhan, Saif Ullah, Waseem, Muath Suliman, Abdul Baseer Saqib, and Mohammed Qeshta. Deep learning-driven insights into the transmission dynamics of hepatitis b virus with treatment. Scientific Reports, 15(1):23741, 2025. [17] Zulqurnain Sabir, Thongchai Botmart, Muhammad Asif Zahoor Raja, R Sadat, Mo- hamed R Ali, Abdulaziz A Alsulami, Abdullah Alghamdi, et al. Artificial neural network scheme to solve the nonlinear influenza disease model. Biomedical Signal Processing and Control, 75:103594, 2022. [18] Kanit Mukdasai, Zulqurnain Sabir, Muhammad Asif Zahoor Raja, Peerapongpat Singkibud, R Sadat, and Mohamed R Ali. A computational supervised neural network procedure for the fractional siq mathematical model. The European Physical Journal Special Topics, 232(5):535–546, 2023. [19] Muhammad Farhan, Zhi Ling, Saif Ullah, Almetwally M Mostafa, Salman A AlQah- tani, et al. A novel intelligent framework for assessing within-host transmission dy- namics of chikungunya virus using an unsupervised stochastic neural network ap- proach. Computational Biology and Chemistry, 117:108380, 2025. [20] Jamiu Adeyemi Ademosu, Samson Olaniyi, Sulaimon Femi Abimbade, Fu- raha Michael Chuma, Richard Chinedu Ogbonna, Ramoshweu Solomon Lebelo, and Kazeem Oare Okosun. Modelling population dynamics of substance abuse in the pres- ence of addicted immigrant with real data of rehabilitation cases. Journal of Applied T. Khan et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6670 22 of 22 Mathematics, 2025(1):2241780, 2025. [21] 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.