EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6379 ISSN 1307-5543 – ejpam.com Published by New York Business Global Stochastic Modelling of Seasonal Influenza Dynamics: Integrating Random Perturbations and Behavioural Factors Rania Saadeh1, Alshaikh A. Shokeralla2, Naseam Al-Kuleab3, Walid S. Hamad4, Mawada Ali5, Mohamed A. Abdoon6, Fathelrhman El Guma2,∗ 1 Department of Mathematics, Faculty of Science, Zarqa University, Zarqa 13110, Jordan 2 Department of Mathematics, Faculty of Science, Al-Baha University, Al-Aqiq 65931, Saudi Arabia 3 Department of Mathematics and Statistics, College of Science, King Faisal University, Al-Ahsa 31982, Saudi Arabia. 4 Department of Preparatory Year, Al-Ghad College for Applied Medical Sciences, Najran, Saudi Arabia 5 Department of Mathematical Modelling, Faculty of Mathematical Sciences and Statistics, Al-Neelain University, Khartoum, Sudan 6 Department of Basic Sciences, Common First Year Deanship, King Saud University, P.O. Box 1142, Riyadh 12373, Saudi Arabia Abstract. The study proposes a stochastic model to investigate the seasonality of influenza in Saudi Arabia. In contrast to the classical deterministic model, we incorporate internal stochastic ambient noise through white noise perturbations in order to present a more realistic portrayal of the oscillation of sick individuals. The model correctly reproduces the empirical seasonal peak in the number of influenza cases, which is most strongly expressed in epidemiologic week 30, also showing the seasonal outbreak periodicity. Sensitivity analysis demonstrates that the magnitudes of stochastic fluctuations play a crucial role in the prediction uncertainty and the outbreak variabil- ity. Transmission rates and the recovery parameter are the main determinants of the magnitude, timing, and impact on hospitalization of the epidemic wave. A greater transmission rate is consis- tently associated with more intense and prolonged breakouts, whereas a larger transmission rate facilitates a more rapid ascent and descent of the peak. These findings underscore the significance of stochastic factors and pinpoint essential parameters for enhancing public health interventions and resource allocation strategies. Our results provide practical recommendations for enhancing influenza preparation using data-driven stochastic modelling. 2020 Mathematics Subject Classifications: 60H10, 92D30, 34F05 Key Words and Phrases: Stochastic Modeling, Seasonal Influenza, Forecast Uncertainty, Trans- mission Dynamics, White Noise Perturbations, Parameter Sensitivity ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6379 Email addresses: rsaadeh@zu.edu.jo (Rania Saadeh), Ashokralla@bu.edu.sa (Alshaikh A. Shokeralla), nkuleab@ksu.edu.sa (Naseam Al-Kuleab), shamad@gc.edu.sa (Walid S. Hamad), mawadaali23@neelain.edu.sd (Mawada Ali), mabdoon.c@ksu.edu.sa (Mohamed A. Abdoon), fgumu@bau.edu.sa (Fathelrhman El Guma) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 2 of 34 1. Introduction Seasonal flu is an infectious viral respiratory disease that spreads during the au- tumn and winter seasons [1]. The illness comes with typical symptoms such as fever, headache,chills, cough, nasal congestion, body aches and pains. It has been demonstrated that in severe cases of seasonal influenza, mortality can occur [1]. Persons with multiple chronic medical conditions are more likely to become infected with seasonal influenza [2] and experience an increase in its severity with higher rates of emergency room visits and hospitalizations [3]. Seasonal influenza is responsible for killing between 290,000 and 650,000 people annually, according to the WHO in 2022 [4]. In 2020, World Health Ranking estimated that influenza and pneumonia, together, caused the death of 4.58% of the population of Saudi Arabia and its seasonal influenza death rate was 30.9 per 100 000, so that it is 82nd most-frequented country in the world [5]. A systematic review and meta-analysis of 18 studies in Hajj pilgrims reporting on 62,431 individuals estimated a prevalence of 5.9% for influenza A and 3.6% for influenza B [6] Among currently approved viral vaccines, vaccine for the influenza virus is distinc- tive. Vaccines that target influenza disease have to be reformulated annually to provide matching between the vaccine strain and the wild-type viruses that will be present in a season and this is owing to the influenza virus’s ongoing antigenic drift of its HA and NA surface glycoproteins. Licensed inactivated influenza vaccines (IIV) have been available in the United States (US) since 1945[7]. A live attenuated influenza vaccine has been available since it was approved in 2003 [8]. The objective of the WHO Global Vaccine Action Plan (GVAP) is to reduce transmission of this global epidemic by increasing the seasonal influenza vaccine coverage rate [9]. Annual vaccination against influenza is the most effective measure to prevent infec- tions and to reduce the adverse events associated with them, such as visit to the doctor, hospitalization, and deaths [10]. There is variation in VE from year to year depending on the degree of match between the vaccine viruses and those in circulation. Accordingly, vaccine effectiveness (VE) could be weakened by lack of compatible with the strains of influenza which are presently spreading and, hence, Et of the influenza-related outcomes increases [11–14]. However studies suggest that the uptake of seasonal influenza vaccine in Saudi Arabia is low. The references included six Saudi citizens’ studies published from October 19, 2017, to October 18, 2022. Seasonal influenza vaccine coverage was 12.7%-55.0% in these studies [15–20]. Taken together, these studies exemplify that seasonal flu vaccination cov- erage is below that recommended by vaccination guidelines, despite the different methods of measuring vaccine status. As such, monitoring the status of seasonal influenza vacci- nation should be regularly implemented among the residents of the Saudi population in order to guide the efforts of boosting vaccination coverage. Due to the increasing use of fractional calculus in modeling diffusion, control pro- cesses, and viscoelasticity, applied mathematics has gained notable popularity over the R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 3 of 34 past few decades, particularly in physics and engineering research, where fractional dif- ferential equations are widely employed [21–25]. This mathematical discovery has had considerable impact on the practice of epidemiological modelling, given that models play a key role in evaluating the transmission dynamics of infectious diseases [26–30]. The majority of current models are classical deterministic and fractal models, which may ig- nore the inherent natural variability of disease transmission and the background noise or ambient noise. Environmental randomicity is an important factor in the occurrence and progression of epidemics, and there is a greater attention to stochastic models containing randomicity. Stochastic systems and fractional-order mathematical models have been more widely used in the study of the transmission dynamics of infectious diseases in recent years. For example, fractional differential equations and stochastic systems have been used in modeling complex diseases such as visceral leishmaniasis, flu pandemics, and COVID- 19. These methodologies enable better understanding of disease persistence, extinction conditions, and evaluation of control strategies [31–37]. Ali, M. et al. [07] proposed a new stochastic model change point SEIHR to predict COVID-19 in South Africa and to evaluate the effectiveness of NPIs. The study developed solutions using stochastic Lyapunov function theory and analyzed data gathere between April and September 2021 to explore the dynamics of the pandemic. The model tests the efficacy of lockdowns in preventing the spread of viruses. Alnafisah and El-Shahed [08] proposed a stochastic Hantavirus infection model. were established for the solution in the space determined space. They also established conditions for a unique ergodic stationary distribution and conditions for extinction of Hantavirus infection. With the Milstein method, they highlighted the role of environmental noise in the model. Liu and Jiang [09] considered a SIR epidemic model with logistic birth rate and pos- sibilitic formulation, and established the global stability of the positive equilibrium by using the analysing method based on the Lyapunov functions and constructed a globally asymptotically stable solution. The stochastic basic reproduction number RS 0 is shown to play a role in determining the threshold dynamics of the system in given circumstances, and sufficient conditions are established for the extinction of the disease and the stably existence of the positive solutions. Kang et al. [38] studied the SIS model in randomly fluctuating environments. They used parameter perturbation to investigate the impact of uncertainties in parameters, such as infection rate, recovery rate. they demonstrated that Stratonovich’s SDE is not well defined once the fluctuation is Gaussian. The reason is that the Stratonovich SDE is more appropriate for the parameter changes in the epidemic model. In contrast, the influence of the Itô SDE scales with the variance. We take one step further in influenza transmission modeling by incorporating stochas- ticity in contrast to the classical deterministic modeling. The study investigates what the intensity of white noise does to variability, similar to what happens with social behavior and environmental change. Looking at different levels of noise, the model includes the randomness of the trends in the epidemics and finds the basic parameters ( R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 4 of 34 beta, alpha) that determine the timing, severity, and maximum hospitalization of the outbreaks. The results confirm that stochastic modeling is essential for a successful public health in- tervention. 2. Model formulation We consider an SEIHVR epidemic model, where the total human population at any time t is denoted by N (t). subdivides into six compartments Susceptible individuals (S) , exposed (E), infected (I), hospitalized (H), vaccinated (V ), and recovered (R). Thus, N(t) = S(t) + E(t) + I(t) +H(t) + V (t) +R(t). Each class incurs a fixed natural death at rate µ. The susceptible class is assumed to be increasing by recruitment process at rate λ and by waning immunity of the individuals in classes V and R at rate ρ. It is decreased by vaccination at rate α and by the force of infection f = β (E + I)(S + σV ) N , where σ is vaccination inefficacy, i.e., 1 − σ is vaccine efficacy, and β is effective contact rate of a susceptible or vaccinated individual with an infected or exposed individual . Individuals in the exposed class are recruited by the force of infection f , and move to class I (after spending the incubation period) at rate w and to class R at rate γ1. Vaccination does not confer a long lasting immunity nor full protection is guaranteed. Individuals in class V are recruited from class S at rate α, exposed to the virus at rate σ, and become susceptible again at rate ρ due to waning immunity. Infected individuals are increased at rate ω from class E, recover at rate γ, and ϵ fraction of them are hospitalized. The recovered population is increased from classes E, I, and H at rates γ1, γ2, and γ3, respectively; and become susceptible again at rate ρ. From the above description we derive system (1) of ordinary differential equations to mathematically describe the transmission dynamics of influenza disease. dS dt = λ+ ρ(V +R)− β (E + I)(S + σV ) N − (α+ µ)S dE dt = β (E + I)(S + σV ) N − (ω + γ1 + µ)E dI dt = ωE − (γ2 + ϵ+ µ)I (1) dH dt = ϵI − (γ3 + µ)H dV dt = αS − (ρ+ σ + µ)V dR dt = γ1E + γ2I + γ3H − (ρ+ µ)R R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 5 of 34 with initial conditions S(0) ≥ 0, E (0) ≥ 0, I(0) ≥ 0, H(0) ≥ 0, V (0) ≥ 0, R(0) ≥ 0. Table 1: Description Parameters in the model Parameter Description λ Recruitment rate of individuals into the population µ Natural death rate β Average effective contact rate α Vaccination rate σ Vaccine inefficacy ϵ Fraction of infected individuals who hospitalized γ1 Recovery rate for exposed individuals γ2 Recovery rate for infected individuals γ3 Recovery rate for hospitalized individuals 1 ω Average latent or incubation period ρ Rate at which individuals lose immunity The stochastic differential equation model (2) correspondence of the deterministic model (1) is given by dS(t) = ( λ+ ρ(V +R)− β (E + I)(S + σV ) N − (α+ µ)S ) dt +δ1 S(t) dB1(t) dE(t) = ( β (E + I)(S + σV ) N − (ω + γ1 + µ)E ) dt+ δ2E(t) dB2(t) dI(t) = (ωE − (γ2 + ϵ+ µ)I) dt+ δ3 I(t) dB3(t) (2) dH(t) = (ϵI − (γ3 + µ)H) dt+ δ4H(t) dB4(t) dV (t) = (αS − (ρ+ σ + µ)V ) dt+ δ5 V (t) dB5(t) dR(t) = (γ1E + γ2I + γ3H − (ρ+ µ)R) dt+ δ6R(t) dB6(t) Where B1(t), B2(t), B3(t), B4(t), B5(t), B6(t) as independent standard Brownian motions, and δ1, δ2, δ3, δ4, δ5, δ6 as the intensities of the standard Gaussian white noises,respectively. 3. Analysis of the Model 3.1. Existence and uniqueness of solution to the stochastic model This section discusses the existence and uniqueness of solution of the proposed stochas- tic model (2). Lemma 1. ([39]) For any initial condition (S(0), E(0), I(0), H(0), V (0), R(0)) ∈ R6 +, there exists a unique solution (S(t), E(t), I(t), H(t), V (t), R(t)) of stochastic model (2) on t ≥ 0, which remains in R6 + with probability one for all t ≥ 0. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 6 of 34 Proof. As for initial condition of the state variables (S(0), E(0), I(0), H(0) , V (0), R(0)) ∈ R6 +, the coefficients used in equations are continuous and local Lipschitz condition. Hence, there must exists a local unique solution (S(t), E(t), I(t), H(t), V (t), R(t)) of the stochastic model (2) over t ∈ [o, τe), τe is the explosion time. In order to prove that the solution is global, we only need to prove that τe = ∞ a.s. Assume that we suppose that there exists z0 such that all of the initial conditions on the state lie within { 1 z0 , z0} L̇et us define for each positive integer z ≥ z0, define the finishing time as follows τz = { t ∈ [o, τe) : min{S,E, I,H, V,R} ≤ 1 z or max{S,E, I,H, V,R} ≥ z } Let inf ϕ = ∞ where ϕ denotes the null set. Definition of τz and as z → ∞, we say that it is increasing. Assume τ∞ = limz→∞ τz, obviously τ∞ ≤ τe a.s . Upon showing τ∞ = ∞ a.s., we suppose that τe = ∞ and hence (S(t), E(t), I(t), H(t), V (t), R(t)) will lie in ∈ R6 + a.s. ∀t ≥ 0. Hence, it is sufficient to prove that τe = ∞ a.s. if not, there must exists tow positive constants ζ ∈ (0, 1) and T such that P{T ≥ τ∞} > ζ, we can fined an integer z1 ≥ z0, such that P{τz ≤ T} ≥ ζ, ∀z ≥ z1. Define a C2− function G : R6 + −→ R+ by G(S(t), E(t), I(t), H(t), V (t), R(t)) = S + E + I +H + V +R− lnS − 6 − lnE − ln I − lnH − lnV − lnR (3) We note that G is a non-negative function, such that 0 ≤ g − ln g − 1,∀g > 0. Assume that z0 ≤ z and T > 0 are arbitrary. Upon applying Itô’s formula to equation (3), So we get dG(S,E, I,H, V,R) = LG(S,E, I,H, V,R) + δ1(S − 1)dB1 + δ2(E − 1)dB2 +δ3(I − 1)dB3 + δ4(H − 1)dB4 + δ5(V − 1)dB5 +δ6(R− 1)dB6 (4) In equation (4), LG(S,E, I,H, V,R) : R6 + −→ R+is defined by LG = (1− 1 S ) ( λ+ ρ(V +R)− β (E + I)(S + σV ) N − (α+ µ)S ) +(1− 1 E ) ( β (E + I)(S + σV ) N − (ω + γ1 + µ)E ) +(1− 1 I ) (ωE − (γ2 + ϵ+ µ)I) +(1− 1 H ) (ϵI − (γ3 + µ)H) +(1− 1 V ) ( αS − (ρ+ σ + µ)V ) +(1− 1 R ) (γ1E + γ2I + γ3H − (ρ+ µ)R) R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 7 of 34 + δ21 + δ22 + δ23 + δ24 + δ25 + δ26 2 ≤ λ+ α+ 5µ+ ω + γ1 + γ2 + γ3 + ϵ+ σ + 2ρ+ β + βσ + δ21 + δ22 + δ23 + δ24 + δ25 + δ26 2 = K (5) Then dG(X) ≤ Kdt+ δ1(S − 1)dB1 + δ2(E − 1)dB2 + δ3(I − 1)dB3 +δ4(H − 1)dB4 + δ5(V − 1)dB5 + δ6(R− 1)dB6 Integrating both sides from 0 to τz ∧ T , and taking expectations, we obtain E[G(S(τz ∧ T )), E(τz ∧ T ), I(τz ∧ T ), H(τz ∧ T ), V (τz ∧ T ), R(τz ∧ T ))] ≤ G(S(0), E(0), I(0), H(0), V (0), R(0)) + E (∫ τz∧T 0 Kdt ) ≤ G(S(0), E(0), I(0), H(0), V (0), R(0)) +KT (6) Then E[G] ≤ G(S(0), E(0), I(0), H(0), V (0), R(0)) +KT for any positive z ≥ z1, we set Ωz = (τz < T ), this leads to P (Ωz) ≥ ζ. Note that for each ω ∈ Ωz there must exist one more than one S(τz, ω), E(τz, ω), I(τz, ω), H(τz, ω), V (τz, ω), R(τz, ω) which equals 1 z or z, Consequently, G(S(τz), E(τz), I(τz), H(τz), V (τz), R(τz)) is no less then 1 z − 1 + ln z or z − 1− ln z. Therefore, G(S(τz), E(τz), I(τz), H(τz), V (τz), R(τz)) ≥ ( 1 z − 1 + ln z) ∧ (z − 1− ln z) So we obtain G(S(0), E(0), I(0), H(0), V (0), R(0)) +KT ≥ E[IΩωG(S(τz), E(τz), I(τz), H(τz), V (τz), R(τz))] = P (ΩK)G(S(τz), E(τz), I(τz), H(τz), V (τz), R(τz)) > ζ [ (1z − 1 + ln z) ∧ (z − 1− ln z) ] IΩω is the indicator function of Ωω. Set z −→ ∞, we have ∞ > G(S(0), E(0), I(0), H(0), V (0), R(0)) +KT > ∞ showing that t∞ = ∞, a.s. Lemma 2. For any positive solution (S(t), E(t), I(t), H(t), V (t), R(t)) of stochastic model (2) with initial value (S(0), E(0), I(0), H(0), V (0), R(0)) ∈ R6 + we have max{limt→∞ supS(t), limt→∞ supE(t), limt→∞ sup I(t), limt→∞ supH(t), limt→∞ supV (t), limt→∞ supR(t)} ≤ λ µ , a.s. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 8 of 34 Proof. From stochastic model (2), we have lim t→∞ d(S(t) + E(t) + I(t) +H(t) + V (t) +R(t)) dt = lim t→∞ (Λ− µ(S(t) +E(t) + I(t) +H(t) +V (t) +R(t))− σV ) lim t→∞ (S(t) + E(t) + I(t) +H(t) + V (t) +R(t)) ≤ λ µ − λ µ e−µt ≤ λ µ Then obviously we obtain limt→∞ supS(t) ≤ λ µ , limt→∞ supE(t) ≤ λ µ , limt→∞ sup I(t) ≤ λ µ limt→∞ supH(t) ≤ λ µ , limt→∞ supV (t) ≤ λ µ , limt→∞ supR(t) ≤ λ µ , a.s. 3.2. Extinction In this section, we investigate the conditions for the extinction of the disease in the stochastic model (2) under the white noise stochastic disturbance. The following theorem gives conditions for I(t) and H(t) got extinction. Theorem 1. If (γ2 + ϵ + µ) + 1 2δ 2 3 > ωλ µ , then I(t) go to extinction almost surely. If (γ3 + µ) + 1 2δ 2 4 > ϵλ µ , then H(t) go to extinction almost surely. Proof. Let (S(t), E(t), I(t), H(t), V (t), R(t)) be a solution of stochastic model (2) with initial value (S(0), E(0), I(0), H(0), V (0), R(0)) ∈ R6 +. Applying Itô’s formula to the third equation of stochastic model (2), we get dI(t) = ( ωE I − (γ2 + ϵ+ µ) ) I(t)dt+ δ3I(t)dB3(t) d(ln I(t)) = ( f(t)− 1 2 g2(t) ) dt+ g(t)dB(t) = ( ωE I − (γ2 + ϵ+ µ)− 1 2 δ23 ) dt+ δ3 dB3(t) integration from 0 to t , then ln I(t)− ln I(0) = ∫ t 0 ( ωE I − (γ2 + ϵ+ µ)− 1 2 δ23 ) du+ ∫ t 0 δ3dB3 ln I(t) ≤ ln I(0) + ∫ t 0 ( ωλ µ − (γ2 + ϵ+ µ)− 1 2 δ23 ) du+ ∫ t 0 δ3dB3 ≤ ln I(0)− ( (γ2 + ϵ+ µ) + 1 2 δ23 − ωλ µ ) t+ δ3B3 R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 9 of 34 Dividing by t, we have ln I(t) t ≤ ln I(0) t − ( (γ2 + ϵ+ µ) + 1 2 δ23 − ωλ µ ) + δ3B3 t By using the strong law of large numbers [40], we obtain lim t→∞ δ3B3 t = 0 lim t→∞ ln I(t) t ≤ − ( (γ2 + ϵ+ µ) + 1 2 δ23 − ωλ µ ) (7) Since (γ2 + ϵ+ µ) + 1 2δ 2 3 > ωλ µ , taking the limit superior of both sides leads to lim t→∞ sup ln I(t) t ≤ − ( (γ2 + ϵ+ µ) + 1 2 δ23 − ωλ µ ) < 0, a.s. (8) which implies limt→∞ I(t) = 0 From the fourth equation of system dH(t) = ( ϵI H − (γ3 + µ) ) H(t)dt+ δ4H(t)dB4(t) d(lnH(t)) = ( ϵI H − (γ3 + µ)− 1 2 δ24 ) dt+ δ4 dB4(t) integration from 0 to t , then lnH(t)− lnH(0) = ∫ t 0 ( ϵI H − (γ3 + µ)− 1 2 δ24 ) du+ ∫ t 0 δ4dB4 lnH(t) ≤ lnH(0)− ( (γ3 + µ) + 1 2 δ24 − ϵλ µ ) t+ δ4B4 Dividing by t, we have lnH(t) t ≤ lnH(0) t − ( (γ3 + µ) + 1 2 δ24 − ϵλ µ ) + δ4B4 t lim t→∞ lnH(t) t ≤ − ( (γ3 + µ) + 1 2 δ24 − ϵλ µ ) (9) Since (γ3 + µ) + 1 2δ 2 4 > ϵλ µ , taking the limit superior of both sides leads to lim t→∞ sup lnH(t) t ≤ − ( (γ3 + µ) + 1 2 δ24 − ϵλ µ ) < 0, a.s. (10) which implies limt→∞H(t) = 0 R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 10 of 34 3.3. Persistence and Stationary Distributions As for as the stochastic model are concerned, they have no endemic equilibrium. Thus, the stability analysis cannot be used as a tool for studying the disease persistence. As a result, we shall establish sufficient conditions for the existence of a unique ergodic station- ary distribution.which in some sense, will work for persistence of the disease. Define a parameter R∗ 0 = σβωα (α+ µ+ 1 2δ 2 1)(ω + γ1 + µ+ 1 2δ 2 2)(γ2 + ϵ+ µ+ 1 2δ 2 3)(ρ+ σ + µ+ 1 2δ 2 5) Theorem 2. The solution (S(t), E(t), I(t), H(t), V (t), R(t)) of the stochastic model (2) is ergodic as well as there is a unique stationary distribution π(·) whenever R∗ 0 > 1. Proof. In view of Lemma (1) we have obtained that for any initial value (S(0), E(0), I(0), H(0), V (0), R(0)) ∈ R6 +, there is a unique solution (S(t), E(t), I(t), H(t), V (t), R(t)) ∈ R6 + . The diffusion ma- trix of stochastic model (2) is given by A =  δ21S 2 0 0 0 0 0 0 δ22E 2 0 0 0 0 0 0 δ23I 2 0 0 0 0 0 0 δ24H 2 0 0 0 0 0 0 δ25V 2 0 0 0 0 0 0 δ26R 2  (11) If we choose M = min (S,E,I,H,V,R)∈D̄⊂R6 + { δ21S 2, δ22E 2, δ23I 2, δ24H 2, δ25V 2, δ26R 2 } We obtain 6∑ i,j=0 aij(S,E, I,H, V,R)ξi ξj ≥ δ21S 2 ξ21 + δ22 E 2 ξ22 + δ23I 2 ξ23 +δ24 H 2 ξ24 + δ25 V 2 ξ25 + δ26 R 2 ξ26 ≥ M |ξ|2 where (S,E, I,H, V,R) ∈ D̄, ξ = {ξ1, ξ2, ξ3, ξ4, ξ5, ξ6} ∈ R6 + Then the first condition in Lemma 4.1 in [41] is satisfied. To prove the second condition we construct a C2− function V : R6 + −→ R+ in the following form V1(S,E, I,H, V,R) = S + E + I +H + V +R− w1 lnS −w2 lnE − w3 ln I − w4 lnV (12) R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 11 of 34 where w1, w2, w3, w4 are the positive constant . Using Itô’s formula, we have LV1 = w1 ( −λ S − ρ(V +R) S + β (E + I)(S + σV ) SN + (α+ µ) + 1 2 δ21 ) +λ− µN − σV + w2(−β (E + I)(S + σV ) EN + (ω + γ1 + µ) + 1 2 δ22) +w3(− ωE I + (γ2 + ϵ+ µ) + 1 2 δ23) + w4(− αS V + (ρ+ σ + µ) + 1 2 δ25) ≤ λ− 5 ( σV × w1λ S × w2βI E × w3ωE I × w4αS V ) 1 5 − w2β +w1(− ρ(V +R) S + β (E + I)(S + σV ) SN + (α+ µ) + 1 2 δ21) +w2((ω + γ1 + µ) + 1 2 δ22) + w3((γ2 + ϵ+ µ) + 1 2 δ23) +w4((ρ+ σ + µ) + 1 2 δ25) Using the inequality x+ y ≥ 2 √ xy, x, y > 0, leads to LV1 ≤ −5 (σw1λw2βw3ωw4α) 1 5 + λ+ w1β(E + I)(S + σV ) SN +w1(α+ µ+ 1 2 δ21) + w2(ω + γ1 + µ+ 1 2 δ22) +w3(γ2 + ϵ+ µ+ 1 2 δ23) + w4(ρ+ σ + µ+ 1 2 δ25) let w1 = λ (α+ µ+ 1 2δ 2 1) , w2 = λ (ω + γ1 + µ+ 1 2δ 2 2) w3 = λ (γ2 + ϵ+ µ+ 1 2δ 2 3) , w4 = λ (ρ+ σ + µ+ 1 2δ 2 5) Then LV1 ≤ w1β(E + I)(S + σV ) SN − 5λ (( σβ (α+ µ+ 1 2δ 2 1)(ω + γ1 + µ+ 1 2δ 2 2) × ωα (γ2 + ϵ+ µ+ 1 2δ 2 3)(ρ+ σ + µ+ 1 2δ 2 5) ) 1 5 − 1  let R∗ 0 = σβωα (α+ µ+ 1 2δ 2 1)(ω + γ1 + µ+ 1 2δ 2 2)(γ2 + ϵ+ µ+ 1 2δ 2 3)(ρ+ σ + µ+ 1 2δ 2 5) Then LV1 ≤ −5λ(R∗ 0 1 5 − 1) + w1β(E + I)(S + σV ) SN (13) R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 12 of 34 In addition, we obtain V2(S,E, I,H, V,R) = w5(S + E + I +H + V +R− w1 lnS − w2 lnE −w3 ln I − w4 lnV )− lnS − lnH − lnR+ S + E +I +H + V +R = (w5 + 1) (S + E + I +H + V +R) −(w1w5 + 1) lnS − w2w5 lnE − w3w5 ln I −w4w5 lnV − lnH − lnR Where w5 is positive constant. It is easy to check that Notice that, let z → ∞,then we have lim z→∞,(S,E,I,H,V,R)∈R6 +\Uz inf V2(S,E, I,H, V,R) = +∞ (14) where Uz = (1z , z)× (1z , z)× (1z , z)× (1z , z)× (1z , z)× (1z , z). The next step is to prove that V2(S,E, I,H, V,R) has one and only one minimum value V2(S(0), E(0), I(0), H(0), V (0), R(0)). The partial derivative of V2(S,E, I,H, V,R) with respect to S,E, I,H, V,R is as follow ∂V2(S,E, I,H, V,R) ∂S = 1 + w5 − (w5w1 + 1) S ∂V2(S,E, I,H, V,R) ∂E = 1 + w5 − w5w2 E ∂V2(S,E, I,H, V,R) ∂I = 1 + w5 − w5w3 I ∂V2(S,E, I,H, V,R) ∂H = 1 + w5 − 1 H ∂V2(S,E, I,H, V,R) ∂V = 1 + w5 − w5w4 V ∂V2(S,E, I,H, V,R) ∂R = 1 + w5 − 1 R We can easily show that V2 has unique stagnation point (S(0), E(0), I(0), H(0), V (0), R(0)) =( w5w1 + 1 1 + w5 , w5w2 1 + w5 , w5w3 1 + w5 , 1 1 + w5 , w5w4 1 + w5 , 1 1 + w5 ) (15) R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 13 of 34 Moreover, the Hessain matrix of V2(S,E, I,H, V,R) at (S(0), E(0), I(0), H(0), V (0), R(0)) is H =  w5w1 + 1 S2(0) 0 0 0 0 0 0 w5w2 E2(0) 0 0 0 0 0 0 w5w3 I2(0) 0 0 0 0 0 0 1 H2(0) 0 0 0 0 0 0 w5w4 V 2(0) 0 0 0 0 0 0 1 R2(0)  The Hessian matrix is positive definite. Thus, V2(S,E, I,H, V,R) has a minimum value V2(S(0), E(0), I(0), H(0), V (0), R(0)). From the continuity of V2(S,E, I,H, V,R) and according to Equation (14), one can say that V2(S,E, I,H, V,R) has one and only one minimum value V2(S(0), E(0), I(0), H(0), V (0), R(0)) in R6 +. Now , we will define a non-negative Lyapunov C2− function V : R6 + −→ R+ as follows V = V2(S,E, I,H, V,R)− V2(S(0), E(0), I(0), H(0), V (0), R(0)) Applying the Itô’s formula and using the stochastic model (2), we obtain LV ≤ −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ+ 3µ+ α+ γ3 + ρ −µ(S + E + I +H + V +R)− σV − λ S − ρ(V +R) S − ϵI H −γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 (16) Where w6 = 5λ(R∗ 0 1 5 − 1) > 0 We next define the bounded closed set D = { ε1 < S < 1 ε1 , ε2 < E < 1 ε2 , ε3 < I < 1 ε3 , ε4 < H < 1 ε4 , ε5 < V < 1 ε5 , ε6 < R < 1 ε6 } where εi > 0, (i = 1, 2, 3, 4, 5, 6), we divide the whole R6 +\D into the following domains D1 = { (S,E, I,H, V,R) ∈ R6 + : 0 < S ≤ ε1 } D2 = { (S,E, I,H, V,R) ∈ R6 + : , S ≥ ε1, V ≥ ε5, 0 < I ≤ ε3 } D3 = { (S,E, I,H, V,R) ∈ R6 + : , S ≥ ε1, V < ε5, 0 < I ≤ ε3 } D4 = { (S,E, I,H, V,R) ∈ R6 + : S ≥ ε1, I ≥ ε3, H ≥ ε4, 0 < E < ε2 } R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 14 of 34 D5 = { (S,E, I,H, V,R) ∈ R6 + : 0 < R < ε6, I ≥ ε3, H ≥ ε4 } D6 = { (S,E, I,H, V,R) ∈ R6 + : 0 < R < ε6, I ≥ ε3, H < ε4 } D7 = { (S,E, I,H, V,R) ∈ R6 + : S ≥ 1 ε1 } D8 = { (S,E, I,H, V,R) ∈ R6 + : E ≥ 1 ε2 } D9 = { (S,E, I,H, V,R) ∈ R6 + : I ≥ 1 ε3 } D10 = { (S,E, I,H, V,R) ∈ R6 + : C ≥ 1 ε4 } D11 = { (S,E, I,H, V,R) ∈ R6 + : V ≥ 1 ε5 } D12 = { (S,E, I,H, V,R) ∈ R6 + : R ≥ 1 ε6 } where εi > 0, is small enough. In the set R6 +\D , we can choose εi > 0 sufficiently small and satisfying − λ ε1 + C1 ≤ −j (17) −w5w6 + (w5w1 + 1) βε3 N + C2 ≤ −j (18) −w5w6 + (w5w1 + 1) βε3 N + (w5w1 + 1) βε3σε1 N + C3 ≤ −j (19) (w5w1 + 1) βε2 N + C4 ≤ −j (20) −γ2ε6 − γ3ε6 + C5 ≤ −j (21) −γ2ε4 − ϵε6 + C6 ≤ −j (22) − µ ε1 + C7 ≤ −j (23) − µ ε2 + C8 ≤ −j (24) − µ ε3 + C9 ≤ −j (25) − µ ε4 + C10 ≤ −j (26) −(µ+ σ) ε5 + C11 ≤ −j (27) − µ ε6 + C12 ≤ −j (28) Next, we will show that LV(S,E, I,H, V,R) ≤ −j on R6 +\D, which is equivalent to proving it on the above twelfth domains. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 15 of 34 Case 1. If (S,E, I,H, V,R) ∈ D1 ,then by equation (16) we get LV ≤ −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ+ 3µ+ α+ γ3 + ρ −µ(S + E + I +H + V +R)− σV − λ S − ρ(V +R) S − ϵI H −γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 ≤ −λ S + C1 ≤ − λ ε1 + C1 where C1 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ+ 3µ +α+ γ3 + ρ− µ(S + E + I +H + V +R)− σV − ρ(V +R) S −ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } According to (17), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D1. Case 2. If (S,E, I,H, V,R) ∈ D2 ,then by equation (16) we get LV ≤ −w5w6 + (w5w1 + 1) βI N + C2 ≤ −w5w6 + (w5w1 + 1) βε3 N + C2 where C2 = sup (S,E,I,H,V,R)∈R6 + { (w5w1 + 1) βE(S + σV ) SN + (w5w1 + 1) βIσV SN + λ +3µ+ α+ γ3 + ρ− µ(S + E + I +H + V +R)− σV − λ S −ρ(V +R) S − ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } According to (18), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D2. Case 3. If (S,E, I,H, V,R) ∈ D3 ,then LV ≤ −w5w6 + (w5w1 + 1) βI N + (w5w1 + 1) βIσV SN + C3 ≤ −w5w6 + (w5w1 + 1) βε3 N + (w5w1 + 1) βε3σε5 ε1N + C3 ≤ −w5w6 + (w5w1 + 1) βε3 N + (w5w1 + 1) βε3σε1 N + C3 Choosing ε5 = ε21, where C3 = sup (S,E,I,H,V,R)∈R6 + { (w5w1 + 1) βE(S + σV ) SN + λ+ 3µ+ α+ γ3 + ρ R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 16 of 34 −µ(S + E + I +H + V +R)− σV − λ S − ρ(V +R) S − ϵI H −γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } According to (19), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D3. Case 4. If (S,E, I,H, V,R) ∈ D4 , then LV ≤ (w5w1 + 1) βE N + C4 ≤ (w5w1 + 1) βε2 N + C4 where C4 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) βI(S + σV ) SN + λ+ 3µ+ α +γ3 + ρ+ (w5w1 + 1) βEσV SN − µ(S + E + I +H + V +R) −σV − λ S − ρ(V +R) S − ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } By (20), we can conclude that LV ≤ −j on D4 . Case 5. If (S,E, I,H, V,R) ∈ D5 , then LV ≤ −γ2I R − γ3H R + C5 ≤ −γ2ε3 ε6 − γ3ε4 ε6 + C5 ≤ −γ2ε6 − γ3ε6 + C5 Choosing ϵ3 = ϵ26, ϵ4 = ϵ26 where C5 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ+ 3µ +α+ γ3 + ρ− µ(S + E + I +H + V +R)− σV − λ S −ρ(V +R) S − ϵI H − γ1E R + δ21 + δ24 + δ26 2 } By (21), we can conclude that LV ≤ −j on D5 . Case 6. If (S,E, I,H, V,R) ∈ D6 , then LV ≤ −γ2I R − ϵI H + C6 R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 17 of 34 ≤ −γ2ε3 ε6 − ϵε3 ε4 + C6 ≤ −γ2ε4 − ϵε6 + C6 Choosing ϵ3 = ϵ6ϵ4 where C6 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ+ 3µ +α+ γ3 + ρ− µ(S + E + I +H + V +R)− σV − λ S −ρ(V +R) S − γ3H R − γ1E R + δ21 + δ24 + δ26 2 } According to (22), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D6. Case 7. If (S,E, I,H, V,R) ∈ D7 , then LV ≤ −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ+ 3µ+ α+ γ3 +ρ− µ(S + E + I +H + V +R)− σV − λ S − ρ(V +R) S −ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 ≤ −µS + C7 ≤ − µ ε1 + C7 where C7 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ +3µ+ α+ γ3 + ρ− µ(E + I +H + V +R)− σV − λ S −ρ(V +R) S − ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } According to (23), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D7. Case 8. If (S,E, I,H, V,R) ∈ D8 , then LV ≤ −µE + C8 ≤ − µ ε2 + C8 where C8 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ +3µ+ α+ γ3 + ρ− µ(S + I +H + V +R)− σV − λ S R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 18 of 34 −ρ(V +R) S − ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } According to (24), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D8. Case 9. If (S,E, I,H, V,R) ∈ D9 , then LV ≤ −µI + C9 ≤ − µ ε3 + C9 where C9 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ +3µ+ α+ γ3 + ρ− µ(S + E +H + V +R)− σV − λ S −ρ(V +R) S − ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } According to (25), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D9. Case 10. If (S,E, I,H, V,R) ∈ D10 , then LV ≤ −µH + C10 ≤ − µ ε4 + C10 where C10 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ +3µ+ α+ γ3 + ρ− µ(S + E + I + V +R)− σV − λ S −ρ(V +R) S − ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } According to (26), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D10. Case 11. If (S,E, I,H, V,R) ∈ D11 , then LV ≤ −µV − σV + C11 ≤ −(µ+ σ) ε5 + C11 where C11 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ +3µ+ α+ γ3 + ρ− µ(S + E + I +H +R)− ρ(V +R) S −λ S − ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 19 of 34 According to (27), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D11. Case 12. If (S,E, I,H, V,R) ∈ D12 , then LV ≤ −µR+ C12 ≤ − µ ε6 + C12 where C12 = sup (S,E,I,H,V,R)∈R6 + { −w5w6 + (w5w1 + 1) β(E + I)(S + σV ) SN + λ +3µ+ α+ γ3 + ρ− µ(S + E + I +H + V )− σV − λ S −ρ(V +R) S − ϵI H − γ1E R − γ2I R − γ3H R + δ21 + δ24 + δ26 2 } According to (28), we have LV ≤ −j for any (S,E, I,H, V,R) ∈ D12. Obviously, from equations (17)–(28), one can see for a sufficiently small εi that LV(S,E, I,H, V,R) ≤ −j for all (S,E, I,H, V,R) ∈ R6 +\D. Consequently, condition two in Lemma 4.1 in [41] is satisfied. This show that stochastic model (2) is ergodic and has a unique stationary distribution. 4. Numerical Simulations This section concerns the estimation of parameters for the deterministic model (1) through weekly data of reported seasonal influenza cases from Saudi Arabia in 2022. The data were obtained from the official WHO database for influenza surveillance[42]. Param- eters have been estimated via a Bayesian inference approach and the Monte Carlo Markov chain methodology [43–45]. This facilitates the assessment of parameter uncertainty while calibrating the model to the observed time series data. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 20 of 34 Table 2: Estimated values from the data Influenza parameters Estimated value Reference λ 100 Fitted µ 0.01246293 Fitted β 0.36983069 Fitted α 0.18302843 Fitted σ 0.89689644 Fitted ϵ 0.69978458 Fitted γ1 0.16348621 Fitted γ2 0.49999414 Fitted γ3 0.18814877 Fitted ω 0.10000064 Fitted ρ 0.00038610 Fitted 0 5 10 15 20 25 30 35 40 45 50 Week in year 2022 0 20 40 60 80 100 120 W ee kl y In fe ct ed c as es Observed I(t) Fitted I(t) Figure 1: Model fit to weekly influenza cases in Saudi Arabia, 2022. The fit of the deterministic model model(1) to the weekly reported influenza cases in Saudi Arabia for the year 2022 is shown in Figure1. The observed data spots are marked by red dots, and the model trend is drawn as the blue curve. The. But it does capture very well the epidemic nature. The seasonal peak (height during epidemiological week 30) appears well and is well described by the model. This alignment emphasizes the model’s potential for identifying the grand epidemiological profile, as well as for yielding guidance for public health interventions. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 21 of 34 0 5 10 15 20 25 30 35 40 45 50 Weeks 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 S( t) 1e7 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 200 400 600 800 1000 1200 E( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 20 40 60 80 100 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 50 100 150 200 250 300 350 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 1 2 3 4 5 6 V( t) 1e6 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 1000 2000 3000 4000 R( t) Figure 2: Stochastic simulations vs. deterministic model (1) (red line) using estimated parameters from Table 2 with δi = 0.1. Figure 2 shows how the determined and the stochastic model results are associated in six compartments of the influenza system in 52 weeks. The deterministic trajectories are represented as red solid lines, indicating a singular solution trajectory derived from the average behaviour. The stochastic encoding (for δ1 = δ2 = δ3 = δ4 = δ5 = δ6 = 0.1), on the other hand, provides a distribution of potential epidemic outcomes, applying randomness coming from fluctuations in the environment and the population behavior. The strong alignment between the outcomes of the stochastic framework and the deterministic solution indicates that the mean-field model is valid in the low-noise limit. The findings indicate R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 22 of 34 that the stochastic model effectively encapsulates the uncertainty and variability inherent in the epidemic’s progression, offering a more comprehensive understanding of the potential dynamics of the outbreak. 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 20 40 60 80 100 120 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 20 40 60 80 100 120 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 20 40 60 80 100 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 25 50 75 100 125 150 175 I(t ) Figure 3: Stochastic simulations of model (2) with δ1 = δ2 = δ5 = δ6 = 0.1 and varying δ3 = δ4 ∈ {0.05, 0.2, 0.25, 0.3}. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 23 of 34 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 50 100 150 200 250 300 350 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 50 100 150 200 250 300 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 50 100 150 200 250 300 350 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 100 200 300 400 500 H( t) Figure 4: Stochastic simulations of model (2) with δ1 = δ2 = δ5 = δ6 = 0.1 and varying δ3 = δ4: 0.05, 0.2 (top), 0.25, 0.3 (bottom). The effect of varying levels of white noise intensities (δ3, δ4) on the dynamics of in- fluenza transmission over 52 weeks is presented using stochastic simulations in Figures3 and4. Each subplot shows simulated trajectories of the epidemic for a range of noise levels in panel (a) and panel (i) with the x-axis the time in weeks and the y-axis the infected or hospitalized counts. We set δ1 = 0.1, δ2 = 0.1, δ5 = 0.1 and δ6 = 0.1 and for the sensitivity analysis we vary δ3 and δ4 systematically. With increasing noise intensity (0.05 ≤ δ3 = δ4 ≤ 0.3), the variation in trajectories also increases, indicating increasing prediction uncertainty for both infection and hospi- talization dynamics. This diffusion shows a method’s dependency on randomly occurring perturbations, highlighting the essential significance of randomness originating from the environment and behavior. In particular, larger levels of noise result in more deviation from the deterministic mean curve, showing that randomness may intensify or dampen the impact of outbreaks depending on the initial condition and the system behavior. These results underscore the need to consider uncertainty when developing epidemic forecasting models and devising resilient public health responses. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 24 of 34 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 250 500 750 1000 1250 1500 1750 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 50 100 150 200 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 2 4 6 8 10 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 1 2 3 4 I(t ) Figure 5: Stochastic simulations with varying β: Top—β = 0.5 and 0.4 yield higher peaks (801–1790 and 170–236 cases); Bottom—β = 0.3 shows minimal infections, while β = 0.2 leads to extinction. Figures5 & 6 show the effect of decreasing the transmission rate β on infection and hospitalization dynamics over 52 weeks of simulation. Figure 5 shows the trajectories of the infection: • With β = 0.5, we observed infections peak during weeks 32–35 with between 801 to 1790 cases. • When β = 0.4, the peak is shifted to between weeks 38–43, with a height between 170–236 cases. • Infectious peaks are much lower and estimated at 11 cases at β = 0.3. • When β = 0.2, the infection dies out: it goes extinct. Figure 6 shows the same statistics for hospitalizations: • The parameter values that provide β = 0.5 and β = 0.4 lead to the cases peaking from 2305 to 4659 in hospitals. • At β = 0.3, there are at most 32 new hospitalisations in a week. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 25 of 34 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 1000 2000 3000 4000 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 100 200 300 400 500 600 700 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 5 10 15 20 25 30 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 1 2 3 4 H( t) Figure 6: Stochastic simulations of model (2) for varying transmission rates β. Top: (left) β = 0.5, hospitaliza- tions peak between 2305 and 4659 cases; (right) β = 0.4, similar peak range observed. Bottom: (left) β = 0.3, peak hospital cases limited to 32; (right) β = 0.2, the epidemic fades out. • For β = 0.2 the disease vanishes completely. These simulations verify that large β causes early and strong outbreaks, while small β suppresses transmission. This sensitivity analysis highlights the importance of the trans- mission parameter to determine outbreak size and guide intervention efforts under the public health strategy. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 26 of 34 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 50 100 150 200 250 300 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 5 10 15 20 25 30 35 40 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 2 4 6 8 10 12 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 1 2 3 4 5 I(t ) Figure 7: Stochastic simulations of model (2) for varying α. As α increases from 0.1 to 0.7, peak infections occur earlier and total cases decrease, indicating faster disease extinction. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 27 of 34 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 100 200 300 400 500 600 700 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 20 40 60 80 100 120 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 10 20 30 40 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 H( t) Figure 8: Stochastic model (2) with varying α. Higher α leads to earlier, smaller peaks in hospital cases, decreasing from 678 (α = 0.1) to 21 (α = 0.7). The effect of the recovery rate α is depicted in Figures 7–10 on the infection and hospitalization dynamics for 52 weeks. Infection paths as α increases are visualized in Figure 7. For α = 0.1, infections are maintained with number of the cases ranging from 123 to 305. In this case, for α = 0.5, 0.7, the peak moves to an earlier time and the total infection decreases albeit naturally (the recovery occurs faster). Hospitalizations are presented in Fig 8. For α = 0.1, the peak hospitalizations range from 320 to 678. Larger α values give rise to lower and earlier peaks, while deaths in hospital have decreased considerably by α = 0.7. 9 dives into infection heterogeneities at α = 0.9, 0.7, 0.5, and 0.3. A larger α corre- sponds to sharper peaks and peaks that occur earlier in time, and also to fewer total cases, emphasizing the role of rapid recovery in outbreak suppression. Figure 10 verifies the above pattern for hospital cases. With α = 0.3 and 0.5, hospital burdens are small, but α = 0.9 leads to brief yet intense peaks. Collectively, higher α leads to a lower epidemic duration and severity, which emphasizes it as a key factor that affects both infection and hospitalization outcomes. These findings confirm the relevance of the stochastic model in representing the time and magnitude of epidemic fadeout at different recovery settings. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 28 of 34 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 20 40 60 80 100 120 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 5 10 15 20 25 30 35 40 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 2 4 6 8 10 12 14 16 I(t ) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 1 2 3 4 I(t ) Figure 9: Stochastic model (2) for different σ values. As σ decreases from 0.9 to 0.3, infection peaks occur later and at lower levels, with extinction observed at σ = 0.3. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 29 of 34 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 50 100 150 200 250 300 350 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 20 40 60 80 100 120 140 160 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 5 10 15 20 25 30 35 40 H( t) 0 5 10 15 20 25 30 35 40 45 50 Weeks 0 1 2 3 4 5 H( t) Figure 10: Stochastic model (2) with varying σ. As σ decreases from 0.9 to 0.3, hospital cases peak later and decline in magnitude, with extinction at σ = 0.3. 5. Conclusion This study uses a stochastic model to assess Saudi Arabia’s weekly seasonal influenza cases in 2022. The model closely aligns with the general trend of observed influenza cases, showing a peak in infections around week 30, which suggests a seasonal epidemic. When comparing deterministic and stochastic models, it becomes evident that the stochastic model excels at capturing a broader spectrum of potential outcomes due to its natural variability. This study delves into the substantial influence of white noise intensity on the variability of infection predictions within the stochastic model, highlighting the crucial role of stochastic factors like environmental or social influences in modeling endeavors. Moreover, the parameter β plays a crucial role in impacting infection and hospitaliza- tion rates, as higher β values are linked to more severe epidemic outcomes. Finally, the parameter α plays a vital role in forecasting the timing and intensity of infection and hospitalization peaks. Higher α values result in an earlier peak and a swift decrease in cases, which may be associated with recovery rates or other factors that help resolve the disease. These findings demonstrate the stochastic model’s sensitivity to critical factors needed to predict illness progression and help build effective public health treatments and epidemic control measures. in the future, we intend to solve some new models, such as in [46–48] and make comparisons with other methods [49–53]. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 30 of 34 Acknowledgements The authors thank the deanship of Zarqa University. This research is funded fully by Zarqa University. Conflicts of Interest: The authors declare that they have no conflict of interest. References [1] World Health Organization. Influenza (seasonal): Overview. influenza (seasonal), 2023. [2] William Schaffner, Janet McElhaney, Albert A Rizzo, Margot Savoy, Allen J Taylor, and Melissa Young. The dangers of influenza and benefits of vaccination in adults with chronic health conditions. Infectious Diseases in Clinical Practice, 26(6):313– 322, 2018. [3] Aimee M Near, Jenny Tse, Yinong Young-Xu, David K Hong, and Carolina M Reyes. Burden of influenza hospitalization among high-risk groups in the united states. BMC Health Services Research, 22(1):1209, 2022. [4] World Health Organization. Global influenza programme, 2022. [5] World Health Rankings. Global influenza programme, 2020. [6] Hamid Safarpour, Meysam Safi-Keykaleh, Iman Farahi-Ashtiani, Jafar Bazyar, Salman Daliri, and Ali Sahebi. Prevalence of influenza among hajj pilgrims: A sys- tematic review and meta-analysis. Disaster medicine and public health preparedness, 16(3):1221–1228, 2022. [7] Anthony E Fiore, Carolyn B Bridges, Jacqueline M Katz, and Nancy J Cox. Inac- tivated influenza vaccines. Vaccines. 6th ed. Philadelphia, PA: Elsevier Inc, pages 257–92, 2012. [8] Plotkin Sl. A short history of vaccination. Vaccines: Expert Consult, pages 5–16, 2008. [9] World Health Organization. Who global vaccine action plan 2011–2020. 2017, 2018. [10] Organisation mondiale de la Santé, World Health Organization, et al. Vaccines against influenza: Who position paper–may 2022–vaccins antigrippaux: note de synthèse de l’oms–mai 2022. Weekly Epidemiological Record= Relevé épidémiologique hebdo- madaire, 97(19):185–208, 2022. [11] Huong Q McLean, Mark G Thompson, Maria E Sundaram, Burney A Kieke, Manjusha Gaglani, Kempapura Murthy, Pedro A Piedra, Richard K Zimmerman, Mary Patricia Nowalk, Jonathan M Raviotta, et al. Influenza vaccine effectiveness in the united states during 2012–2013: variable protection by age and virus type. The Journal of infectious diseases, 211(10):1529–1540, 2015. [12] Manjusha Gaglani, Jessica Pruszynski, Kempapura Murthy, Lydia Clipper, Anne Robertson, Michael Reis, Jessie R Chung, Pedro A Piedra, Vasanthi Avadhanula, Mary Patricia Nowalk, et al. Influenza vaccine effectiveness against 2009 pandemic influenza a (h1n1) virus differed by vaccine type during 2013–2014 in the united states. The Journal of infectious diseases, 213(10):1546–1556, 2016. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 31 of 34 [13] Brendan Flannery, Rebecca J Garten Kondor, Jessie R Chung, Manjusha Gaglani, Michael Reis, Richard K Zimmerman, Mary Patricia Nowalk, Michael L Jackson, Lisa A Jackson, Arnold S Monto, et al. Spread of antigenically drifted influenza a (h3n2) viruses and vaccine effectiveness in the united states during the 2018–2019 season. The Journal of infectious diseases, 221(1):8–15, 2020. [14] Mark W Tenforde, Rebecca J Garten Kondor, Jessie R Chung, Richard K Zimmer- man, Mary Patricia Nowalk, Michael L Jackson, Lisa A Jackson, Arnold S Monto, Emily T Martin, Edward A Belongia, et al. Effect of antigenic drift on influenza vaccine effectiveness in the united states—2019–2020. Clinical Infectious Diseases, 73(11):e4244–e4250, 2021. [15] Ali M Alhazmi, Sulaiman A Alshammari, Hanan A Alenazi, Shaffi A Shaik, Hala M AlZaid, Nouf S Almahmoud, and Hotoon S Alshammari. Community’s compliance with measures for the prevention of respiratory infections in riyadh, saudi arabia. Journal of Family and Community Medicine, 26(3):173–180, 2019. [16] Mazin A Barry, Khalid I Aljammaz, and Abdulaziz A Alrashed. Knowledge, attitude, and barriers influencing seasonal influenza vaccination uptake. Canadian Journal of Infectious Diseases and Medical Microbiology, 2020:1–6, 2020. [17] Alaa A Aljamili. Knowledge and practice toward seasonal influenza vaccine and its barriers at the community level in riyadh, saudi arabia. Journal of Family Medicine and Primary Care, 9(3):1331–1339, 2020. [18] Amel Ahmed Fayed, Abeer Salem Al Shahrani, Leenah Tawfiq Almanea, Nardeen Ibrahim Alsweed, Layla Mohammed Almarzoug, Reham Ibrahim Almuwal- lad, and Waad Fahad Almugren. Willingness to receive the covid-19 and seasonal influenza vaccines among the saudi population and vaccine uptake during the ini- tial stage of the national vaccination campaign: a cross-sectional survey. Vaccines, 9(7):765, 2021. [19] Ibrahim A Sales, Wajid Syed, Majed F Almutairi, and Yazed Al Ruthia. Public knowledge, attitudes, and practices toward seasonal influenza vaccine in saudi arabia: a cross-sectional study. International journal of environmental research and public health, 18(2):479, 2021. [20] Faisal Minshawi, Mohammed Samannodi, Hassan Alwafi, Hamza M Assaggaf, Mo- hammed A Almatrafi, Emad Salawati, Radi Alsafi, Ruba A Alharbi, Raghad F Aldu- ais, Muruj Alrehaili, et al. The influence of covid-19 pandemic on influenza immuniza- tion in saudi arabia: cross-sectional study. Journal of Multidisciplinary Healthcare, pages 1841–1849, 2022. [21] Ahmad Qazza and Rania Saadeh. On the analytical solution of fractional sir epidemic model. Applied Computational Intelligence and Soft Computing, 2023(1):6973734, 2023. [22] Emad Salah, Rania Saadeh, Ahmad Qazza, and Raed Hatamleh. Direct power series approach for solving nonlinear initial value problems. Axioms, 12(2):111, 2023. [23] Rania Saadeh, Osama Ala’yed, and Ahmad Qazza. Analytical solution of coupled hirota–satsuma and kdv equations. Fractal and Fractional, 6(12):694, 2022. [24] Rania Saadeh, Mohammad Abu-Ghuwaleh, Ahmad Qazza, and Emad Kuffi. A fun- R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 32 of 34 damental criteria to establish general formulas of integrals. Journal of Applied Math- ematics, 2022(1):6049367, 2022. [25] Mohamed A Abdoon. First integral method: a general formula for nonlinear fractional klein-gordon equation using advanced computing language. American Journal of Computational Mathematics, 5(2):127–134, 2015. [26] Mohamed A Abdoon et al. Programming first integral method general formula for the solving linear and nonlinear equations. Applied Mathematics, 6(03):568, 2015. [27] Mohammad Abu-Ghuwaleh, Rania Saadeh, and Ahmad Qazza. General master the- orems of integrals with applications. Mathematics, 10(19):3547, 2022. [28] Sayed Saber and Abdullah Alahmari. Impact of fractal-fractional dynamics on pneu- monia transmission modeling. European Journal of Pure and Applied Mathematics, 18(2):5901–5901, 2025. [29] Sayed Saber and Abdullah Alahmari. Mathematical insights into zoonotic disease spread: Application of the milstein method. European Journal of Pure and Applied Mathematics, 18(2):5881–5881, 2025. [30] Mohammed Althubyani, Haroon DS Adam, Ahmad Alalyani, Nidal E Taha, Khdija O Taha, Rasmiyah A Alharbi, and Sayed Saber. Understanding zoonotic disease spread with a fractional order epidemic model. Scientific Reports, 15(1):13921, 2025. [31] Sana Abdulkream Alharbi, Mohamed A. Abdoon, Rania Saadeh, Reima Daher Alsemiry, Reem Allogmany, Mohammed Berir, and Fathelrhman EL Guma. Modeling and analysis of visceral leishmaniasis dynamics using fractional-order operators: A comparative study. Mathematical Methods in the Applied Sciences, 47(12):9918–9937, 2024. [32] Fathelrhman EL Gumaa, Mohamed A Abdoon, Ahmad Qazza, Rania Saadeh, Mo- hammed Ali Arishi, and Abdoelnaser M Degoot. Analyzing the impact of control strategies on visceralleishmaniasis: a mathematical modeling perspective. European Journal of Pure and Applied Mathematics, 17(2):1213–1227, 2024. [33] Rania Saadeh, Mohamed A Abdoon, Ahmad Qazza, Mohammed Berir, Fathel- rhman EL Guma, Naseam Al-Kuleab, and Abdoelnaser M Degoot. Mathematical modeling and stability analysis of the novel fractional model in the caputo derivative operator: A case study. Heliyon, 10(5), 2024. [34] FE Guma, Ossama M Badawy, AG Musa, Badawi Osman Mohammed, Mohamed A Abdoon, Mohammed Berir, and Salih Yousuf Mohamed Salih. Risk factors for death among covid-19 patients admitted to isolation units in gedaref state, eastern sudan: a retrospective cohort study. Journal of Survey in Fisheries Sciences, 10(3s):712–722, 2023. [35] Salem Mubarak Alzahrani and FE Guma. Improving seasonal influenza forecast- ing using time series machine learning techniques. Journal of Information Systems Engineering and Management, 9(4):30195, 2024. [36] Mawada Ali, Fathelrhman EL Guma, Ahmad Qazza, Rania Saadeh, Nahaa E Al- subaie, Mohammed Althubyani, and Mohamed A Abdoon. Stochastic modeling of influenza transmission: Insights into disease dynamics and epidemic management. Partial Differential Equations in Applied Mathematics, 11:100886, 2024. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 33 of 34 [37] Sana Abdulkream Alharbi, Mohamed A Abdoon, Abdoelnaser M Degoot, Reima Da- her Alsemiry, Reem Allogmany, Fathelrhman EL Guma, and Mohammed Berir. Mathematical modeling of influenza dynamics: A novel approach with sveihr and fractional calculus. International Journal of Biomathematics, page 2450147, 2025. [38] Yung-Gyung Kang and Jeong-Man Park. A stochastic susceptible-infected-susceptible epidemic model with stratonovich processes. Journal of the Korean Physical Society, 84(2):158–163, 2024. [39] Xuerong Mao, Glenn Marion, and Eric Renshaw. Environmental brownian noise suppresses explosions in population dynamics. Stochastic Processes and their Appli- cations, 97(1):95–110, 2002. [40] Michael J Panik. Stochastic Differential Equations: An Introduction with Applications in Population Dynamics Modeling. John Wiley & Sons, 2017. [41] Aadil Lahrouz and Lahcen Omari. Extinction and stationary distribution of a stochas- tic sirs epidemic model with non-linear incidence. Statistics & Probability Letters, 83(4):960–968, 2013. [42] programme/surveillance-and monitoring/fluid. https://www.who.int/teams/global influenza-560, volume Available. cited 28 Jul 2022, 2022. [43] Gongzheng Yao, Di Zhang, and Yingbo Liu. A bayesian non-parametric approach for estimating covid-19’s vaccine effectiveness in a stochastic epidemic model. Informatics in Medicine Unlocked, 42:101329, 2023. [44] Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin. Bayesian Data Analysis. CRC press, 2013. [45] Christian P Robert and George Casella. Monte Carlo Statistical Methods. Springer Science & Business Media, 2004. [46] Mawada Ali, Salem Mubarak Alzahrani, Rania Saadeh, Mohamed A Abdoon, Ahmad Qazza, Naseam Al-kuleab, and Fathelrhman EL Guma. Modeling covid-19 spread and non-pharmaceutical interventions in south africa: A stochastic approach. Scientific African, 24:e02155, 2024. [47] Nahaa E. Alsubaie, Fathelrhman EL Guma, Kaouther Boulehmi, Naseam Al-kuleab, and Mohamed A. Abdoon. Improving influenza epidemiological models under caputo fractional-order calculus. Symmetry, 16(7):929, 2024. [48] Fathelrhman EL Guma, Abdelaziz GM Musa, Fatimah Dhafer Alkhathami, Rania Saadehm, and Ahmad Qazza. Prediction of visceral leishmaniasis incidences utilizing machine learning techniques. In 2023 2nd International Engineering Conference on Electrical, Energy, and Artificial Intelligence (EICEEAI), pages 1–6. IEEE, 2023. [49] Mohamed A Abdoon. Fractional derivative approach for modeling chaotic dynamics: Applications in communication and engineering systems. In International Confer- ence on Mathematical Modelling, Applied Analysis and Computation, pages 82–95. Springer, 2025. [50] Faeza Hasan, Mohamed A Abdoon, Rania Saadeh, Mohammed Berir, and Ahmad Qazza. A new perspective on the stochastic fractional order materialized by the exact solutions of allen-cahn equation. International Journal of Mathematical, Engineering and Management Sciences, 8(5):912, 2023. R. Saadeh et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6379 34 of 34 [51] Reem Allogmany, Nada A Almuallem, Reima Daher Alsemiry, and Mohamed A Ab- doon. Exploring chaos in fractional order systems: A study of constant and variable- order dynamics. Symmetry, 17(4):605, 2025. [52] DK Almutairi, Dalal M AlMutairi, Nidal E Taha, Mohammed E Dafaalla, and Mo- hamed A Abdoon. Variable-fractional-order nosé–hoover system: Chaotic dynamics and numerical simulations. Fractal and Fractional, 9(5):277, 2025. [53] Diaa Diaa Eldin Elgezouli, Mohamed Abdoon, Samir Brahim Belhaouari, and Dalal Khalid Almutairi. A novel fractional edge detector based on generalized frac- tional operator. European Journal of Pure and Applied Mathematics, 17(2):1009– 1028, 2024.