EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6082 ISSN 1307-5543 – ejpam.com Published by New York Business Global Advanced Fractional Reaction-Diffusion Modeling for Spatio-Temporal Dynamics of Poliovirus Transmission with Disability Outcomes and Vaccination Impacts Kamel Guedri1, Rahat Zarin2,∗, Basim M. Makhdoum1,3, Hatoon A. Niyazi4 1 Mechanical Engineering Department, College of Engineering and Architecture, Umm Al-Qura University, P.O. Box 5555, Makkah 21955, Saudi Arabia 2 Department of Mathematics, Faculty of Science, King Mongkut’s University of Technology, Thonburi (KMUTT), Bangkok 10140, Thailand 3 King Salman Center for Disability Research, Riyadh 11614, Saudi Arabia 4 Department of Clinical Microbiology and Immunology, Faculty of Medicine, King Abdulaziz University, Jeddah 21589, Saudi Arabia Abstract. This work presents a novel fractional reaction-diffusion model to analyze the spatio- temporal dynamics of poliovirus transmission. Polio, a highly contagious viral infection that pri- marily affects children under five and can lead to permanent disability, spreads through fecal-oral and airborne transmission, often exacerbated by poor sanitation and environmental conditions. Traditional polio models, predominantly based on ordinary differential equations (ODEs), have assumed spatial uniformity an oversimplification of real-world scenarios influenced by popula- tion density, environmental heterogeneity, and mobility. Our study extends classical models by incorporating spatial heterogeneity and memory effects using Caputo fractional derivatives and reaction-diffusion dynamics. The model divides the population into seven epidemiological com- partments: Susceptible (S), Vaccinated (V), Exposed (E), Non-paralytic Infected (Np), Paralytic Infected (P), Recovered (R), and Post-paralytic (A), with corresponding diffusion terms to capture spatial mobility. The inclusion of fractional derivatives accounts for the memory-dependent nature of disease progression, offering a more realistic depiction of poliovirus dynamics. Key contributions include proving the existence and uniqueness of positively bounded solutions, identifying equilib- rium points, and assessing their local and global stability using the basic reproduction number (R0) and Lyapunov functions under fractional dynamics. Sensitivity analysis is conducted to find the most sensitive parameters using the direct differentiation method. Sensitivity analysis highlights critical parameters influencing disease propagation, while numerical simulations validate the theo- retical findings. Graphical results demonstrate the impact of fractional order (α) and vaccination on disease spread, illustrating how memory effects influence convergence to steady states. This fractional reaction-diffusion framework provides valuable insights into poliovirus transmission, of- fering a robust tool for predicting outbreaks and guiding effective intervention strategies in diverse populations. 2020 Mathematics Subject Classifications: 35R11, 92D30, 35K57, 65M70 ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6082 Email address: rahat.zari@mail.kmutt.ac.th (R. Zarin) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 2 of 39 Key Words and Phrases: Fractional reaction-diffusion, polio virus, sensitivity analysis, disabil- ity outcomes, global stability, numerical simulation 1. Introduction Poliomyelitis, commonly referred to as polio, is a viral infection that primarily affects the nervous system. Children under the age of five are considered particularly vulner- able, although individuals of all age groups can contract the disease. Polio is typically transmitted through the ingestion of material contaminated with feces from an infected person. It may also spread through airborne droplets released when an infected individual coughs or sneezes. The virus primarily resides in the throat and intestinal tract of the host. Communities where open defecation is prevalent are at higher risk of polio outbreaks. Cul- tural traditions, a lack of access to sanitation facilities, or both, often contribute to open defecation practices. Additionally, improper waste disposal and unsanitary environmental conditions can facilitate the persistence and spread of the virus within a community. Once the virus enters the body, it targets the nervous system, including the brain. The incuba- tion period for polio ranges from as few as four days to up to 35 days [1]. Currently, there is no known cure for this disease. While many individuals infected with the poliovirus do not exhibit symptoms, they remain capable of transmitting the virus to others. In cases where symptoms do appear, they can include fever, sore throat, headache, vomiting, and fatigue. Severe cases may progress to paralysis, resulting in reduced reflexes, weakened or deformed limbs, and other motor impairments. Recovery from polio is possible, but some individuals may experience a recurrence of symptoms years later, a condition known as post-polio syndrome (PPS). PPS typically affects polio survivors 15 to 40 years after their initial recovery. The exact prevalence and mechanisms of PPS remain uncertain, with proposed explanations including overuse of surviving nerve cells and potential brain damage. The first documented polio outbreak in the United States occurred in Vermont in 1894, resulting in 132 cases. Since that time, significant efforts have been made by medi- cal researchers to address the devastating effects of the virus. In 1988, the World Health Assembly adopted a resolution aimed at eradicating poliomyelitis globally, leading to the establishment of the Global Polio Eradication Initiative (GPEI). This initiative has sig- nificantly reduced the prevalence and incidence of polio in many parts of the world [2]. Despite these efforts, countries such as Afghanistan, Nigeria, and Pakistan continue to report cases of the disease. As long as the poliovirus persists in even one individual, there is a risk of further transmission. Polio remains a global health concern due to its ability to spread silently over several weeks and its potential to travel long distances via land, air, or sea. Vaccination plays a crucial role in helping individuals develop immunity against the poliovirus. However, administering a vaccine during the virus’s incubation period can have harmful effects, as it may increase the viral load in the individual. Agarwal and Bhadauria [3] observed that receiving an inactivated polio vaccine (IPV) during the incubation phase K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 3 of 39 could lead to a progression from the incubation stage to paralysis. For this reason, it is recommended to conduct appropriate screening before administering the vaccine. Compartmental models address refined numerical structures intended to give an item- ized comprehension of the transmission elements of infectious diseases inside human pop- ulation [4–9] . The starting points of these models can be followed back to 1760, when Daniel Bernoulli fostered the principal numerical compartmental model [10]. From that point forward, various compartmental models have been planned to portray the spread of different scourges and irresistible specialists, frequently using old style whole number request subsidiaries [11, 12]. A few compartmental models explicitly address the transmis- sion elements of poliomyelitis. For example, Duque-Marin et al. [13] fostered a model that catches the elements of polio contaminations inside a populace, consolidating two kinds of immunizations, delineation by age gatherings, and transient impacts. Browne et al. [14] proposed a model that records for sickness elements across interconnected locales, in- corporating variables like occasional varieties, ecological repositories, and district explicit intermittent heartbeat immunization techniques. Dénes et al. [15] presented a compart- mental model inspecting the expected effect of unvaccinated people relocating into regions with low immunization inclusion. Various extra investigations have investigated different elements of polio transmission elements, adding to a more profound comprehension of the infection’s spread and control systems [16–18]. In recent years, reaction-diffusion models have been widely utilized to depict the spatial and temporal dynamics of infection transmission, both in vitro and in vivo [19–23]. These models normally expect that all elements engaged with the disease interaction like target cells, contaminated cells, and viral particles follow old style Fickian dissemination with steady dispersion rates. Notwithstanding, trial discoveries, for example, those detailed in [24], propose that the dissemination properties of cells and infections can fluctuate alto- gether contingent upon the encompassing tissues or natural circumstances. For instance, with regards to poliovirus, its portability and determination are intensely affected by the digestive system, brain tissues, or ecological repositories, for example, water and soil, which are normal transmission pathways. To investigate the effect of spatial heterogene- ity on viral elements, scientists have researched response dispersion models custom fitted to explicit infections. For example, Pankavich and Parkinson [25] broke down a response dissemination model for HIV contamination, representing different dispersion rates among target and tainted cells, and concentrated on the worldwide security of the sickness free balance. Likewise, Cai et al. [26] inspected the spatial elements of the Zika infection, consolidating particular dissemination coefficients for different specialists. Applying this way to deal with poliovirus, figuring out the heterogeneous dissemination of viral particles through various mediums like feces, water, or air is critical for precisely displaying its spatial spread. Fractional derivatives (FDs), representing derivatives of non-integer order, have emerged as a vital component in mathematical models (MMs) due to their ability to incorporate K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 4 of 39 memory and hereditary effects. These properties are particularly advantageous for mod- eling processes where the current state depends on historical influences, such as in in- fectious disease dynamics and ecological systems. Unlike classical derivatives, FDs offer greater flexibility and a higher degree of freedom, enabling a more precise representation of complex system dynamics. Consequently, fractional-order models (FOMs) have been extensively applied in various domains, including the study of infectious diseases [27–33], physical processes [34], control theory [35, 36], biology [37], viscoelastic materials [38], and engineering [39]. FOMs possess unique features such as the ability to model anomalous diffusion, ac- count for intrinsic memory effects, capture fractal structures, fit data more effectively, and handle multiscale phenomena. These characteristics make FOMs a powerful framework for incorporating memory and hereditary properties into infectious disease models, surpass- ing the capabilities of integer-order models, which lack non-local interactions and memory effects. The growing interest in fractional calculus (FC) has led to the development of multiple definitions of FDs, each with distinct advantages. Among these, the Caputo and Riemann–Liouville derivatives remain the most commonly used [40, 41]. Fractional-order models based on these derivatives have demonstrated superior accuracy in describing dis- ease progression compared to classical models [42]. Traditional ordinary differential equa- tions (ODEs) often fail to adequately capture the discontinuous nature of disease spread, whereas fractional-order models effectively address these limitations [43]. The use of frac- tional derivatives, particularly the Caputo FD, in mathematical epidemiology has been shown to improve the realism and accuracy of infectious disease dynamics [40, 44]. The novelty of this research lies in extending the classical ordinary differential equa- tion (ODE)-based polio transmission models to a more comprehensive fractional reaction- diffusion framework. While traditional ODE models effectively describe temporal dy- namics, they assume spatial uniformity, an oversimplification in the context of real-world disease transmission influenced by environmental heterogeneity and human mobility. By incorporating diffusion terms into the model, this study accounts for spatial heterogeneity, reflecting the localized impact of population density, sanitation conditions, and vaccina- tion efforts on poliovirus dynamics. Furthermore, replacing classical time derivatives with Caputo fractional derivatives introduces memory effects and non-local temporal interac- tions into the model, which are crucial for accurately capturing the history-dependent nature of disease progression. This combined fractional reaction-diffusion framework of- fers significant advantages over fractional ODEs by integrating both spatial and temporal complexity, providing a more realistic and precise depiction of poliovirus spread. The proposed model not only enhances our understanding of spatially heterogeneous disease dynamics but also offers a valuable tool for informing effective intervention strategies and predicting outbreak scenarios in diverse populations. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 5 of 39 2. Model Formulation The population under study is categorized into seven distinct compartments based on the epidemiological state of individuals. These compartments are: susceptible (S(t)), vac- cinated (V (t)), exposed (E(t)), infected non-paralytic (Np(t)), infected paralytic (P (t)), recovered (R(t)), and post-paralytic (A(t)). The total population size is expressed as: N = S + V + E +Np + P +R+A. Susceptible individuals increase in number due to recruitment through births and immi- gration at a constant rate Λ. They are exposed to the polio virus via contact with envi- ronmental contamination originating from sources such as infected individuals (E,Np, P ). The rate of exposure to the virus is represented by: λ = βcτ (E +Np + P ) N , where β denotes the likelihood of transmission from an infectious source, c represents the average contact rate with contaminated fecal waste per unit time, and τ > 1 is a factor accounting for enhanced environmental contamination due to poor sanitation practices, such as open defecation. Individuals in the susceptible compartment (S) transition to the exposed compartment upon infection, driven by the force of infection λ. Additionally, it is assumed that a proportion of the susceptible population is vaccinated at a rate δ. Once vaccinated, individuals gain temporary immunity and eventually move to the recovered compartment (R) at a recovery rate γ. However, vaccine failure may occur. This can happen if the body either: • Fails to fight the virus introduced by vaccination and becomes infected. • Successfully fights the introduced virus but does not gain immunity, becoming sus- ceptible to infection upon contact with the virus source. In such cases, individuals with vaccination failure move from the vaccinated class into the exposed class with a probability of (1− δ). The natural and polio-induced death rates are µ and σ, respectively. The rates of progression of exposed individuals to the non-paralytic and paralytic stages are α1 and α2, respectively. Recovery rates for individuals in the exposed, non-paralytic, and paralytic classes are ω, ω1, and ω2, respectively. Individuals in the recovered class may progress to the post-paralytic class at rate ρ, while a fraction of them may lose immunity and re-enter the susceptible class at a rate ρ1. Under these assumptions, and following the flow diagram in Figure 1, the dynamics of polio infection are governed by the following system of differential equations, as proposed K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 6 of 39 in [1]:  dS dt = Λ+ ρ1R− (λ+ δ + µ)S, dV dt = δS − (1− δ)λV − (γ + µ)V, dE dt = λS + (1− δ)λV − (α1 + α2 + ω + µ)E, dNp dt = α1E − (σ + µ+ ω1)Np, dP dt = α2E − (σ + µ+ ω2)P, dR dt = γV + ωE + ω1Np + ω2P − (ρ+ ρ1 + µ)R, dA dt = ρR− (σ + µ)A. (1) 2.1. Fractional Reaction-Diffusion Polio Model The classical models for the dynamics of polio virus infection predominantly focus on temporal changes in the population compartments (e.g., [1–3]). These models are often based on lumped parameters, assuming a homogeneously mixed population. However, this assumption oversimplifies the real-world scenario where spatial factors, such as population mobility, environmental heterogeneity, and local climatic conditions, significantly influence disease transmission dynamics. Incorporating Spatial Heterogeneity The baseline model (1) assumes a homogeneous population where individuals are uniformly distributed across the domain. This implies that the infection spreads uniformly without spatial variation. However, in practice, disease spread can vary significantly due to: • Environmental factors (e.g., climate and sanitation conditions), • Population density, • Localized interventions (e.g., vaccination campaigns). To address these complexities, we introduce spatial dependence into the model by adding a diffusion term to represent the mobility of individuals within the domain [45–48]. The K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 7 of 39 resulting reaction-diffusion system is as follows: ∂S(t, x) ∂t − d1∆S(t, x) = Λ + ρ1R(t, x)− (λ+ δ + µ)S(t, x), ∂V (t, x) ∂t − d2∆V (t, x) = δS(t, x)− (1− δ)λV (t, x)− (γ + µ)V (t, x), ∂E(t, x) ∂t − d3∆E(t, x) = λS(t, x) + (1− δ)λV (t, x)− (α1 + α2 + ω + µ)E(t, x), ∂Np(t, x) ∂t − d4∆Np(t, x) = α1E(t, x)− (σ + µ+ ω1)Np(t, x), ∂P (t, x) ∂t − d5∆P (t, x) = α2E(t, x)− (σ + µ+ ω2)P (t, x), ∂R(t, x) ∂t − d6∆R(t, x) = γV (t, x) + ωE + ω1Np(t, x) + ω2P (t, x)− (ρ+ ρ1 + µ)R(t, x), ∂A(t, x) ∂t − d7∆A(t, x) = ρR(t, x)− (σ + µ)A(t, x), where di (i = 1, 2, . . . , 7) are the diffusion coefficients corresponding to each compartment, and ∆ is the Laplacian operator capturing spatial mobility. Incorporating Fractional Dynamics Recent advancements in mathematical modeling emphasize the utility of fractional calcu- lus for capturing memory effects and non-local dynamics in complex systems. Fractional derivatives allow for more accurate modeling of biological processes, as they incorporate the history-dependent nature of the disease spread and population interactions. To further generalize the model, we replace the classical time derivative with the Caputo fractional derivative of order α (0 < α ≤ 1) and retain the diffusion term. The resulting fractional reaction-diffusion system is given by: C 0 D α t S(t, x)− d1∆S(t, x) = Λ + ρ1R(t, x)− (λ+ δ + µ)S(t, x), C 0 D α t V (t, x)− d2∆V (t, x) = δS(t, x)− (1− δ)λV (t, x)− (γ + µ)V (t, x), C 0 D α t E(t, x)− d3∆E(t, x) = λS(t, x) + (1− δ)λV (t, x)− (α1 + α2 + ω + µ)E(t, x), C 0 D α t Np(t, x)− d4∆Np(t, x) = α1E(t, x)− (σ + µ+ ω1)Np(t, x), C 0 D α t P (t, x)− d5∆P (t, x) = α2E(t, x)− (σ + µ+ ω2)P (t, x), C 0 D α t R(t, x)− d6∆R(t, x) = γV (t, x) + ωE(t, x) + ω1Np(t, x) + ω2P (t, x)− (ρ+ ρ1 + µ)R(t, x), C 0 D α t A(t, x)− d7∆A(t, x) = ρR(t, x)− (σ + µ)A(t, x), (2) Here, C 0 D α t denotes the Caputo fractional derivative of order α. The combination of fractional derivatives and diffusion terms enhances the model’s capability to capture both temporal and spatial heterogeneity in disease dynamics. The initial conditions (ICs) for the above model are as follows: S(x, 0) = ψ1(x), V (x, 0) = ψ2(x), E(x, 0) = ψ3(x), Np(x, 0) = ψ4(x), P(x, 0) = ψ5(x), R(x, 0) = ψ6(x), A(x, 0) = ψ7(x), x ∈ Ω. (3) K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 8 of 39 No-flux boundary conditions (BCs) for the model (2) ∂ ∂n S(x, t) = ∂ ∂n V (x, t) = ∂ ∂n E(x, t) = ∂ ∂n Np(x, t) = 0, ∂ ∂n P (x, t) = ∂ ∂n R(x, t) = ∂ ∂n A(x, t) = 0, t ≥ 0, x ∈ ∂Ω (4) where ∂Ω represent the smooth boundary with a bounded domain Ω ⊂ R, and the homo- geneous Neumann boundary conditions mean that no individual crosses the boundary ∂Ω. The coefficients d1,..., d7 are positive constants which control the movement or migration of classes of population, respectively. S E P Np V R A Λ λ δ γ α2 (1− δ)λ ω ω2 ω1 ρ µ µ µ+ σ µ+ σ µ+ σ µ µ ρ1 α1 Figure 1: Flow Diagram of Polio Virus Dynamics 3. Preliminaries We recall the definitions of the Caputo Fractional Partial Derivative (CFPD), Laplace Transform (LT), and Mittag-Leffler function (see, e.g., [49–51]). Definition 1 (Riemann-Liouville Fractional Integral). Let g(t) ∈ L1 (R+). The Riemann- Liouville fractional integral is defined as: ℑα t g(t) = 1 Γ(α) ∫ t 0 (t− ζ)α−1g(ζ) dζ, t > 0, α > 0, (5) K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 9 of 39 with the special case: ℑ0 t g(t) = g(t). (6) Definition 2 (Caputo Fractional Partial Derivative). Let g(x, t) ∈ ACn ([0,+∞),R+). The Caputo fractional derivative of order α is defined as: cDα t g(x, t) = { 1 Γ(n−α) ∫ t 0 (t− ζ)n−α−1 ∂ ng(x,ζ) ∂ζn dζ, n− 1 < α < n, ∂ng(x,t) ∂tn , α = n ∈ N. (7) Definition 3 (Laplace Transform of CFPD). Let G(s) be the Laplace Transform of the function g(t). Then, the LT of CFPD is given by: L{Dα t g(x, t); s} = sαG(x, s)− n−1∑ i=0 sα−i−1g(i)(x, 0), (8) where α ∈ (n− 1, n], n ∈ N. Definition 4 (Mittag-Leffler Function). The two-parameter Mittag-Leffler function Ma,b(x) is defined as: Ma,b(x) = ∞∑ m=0 xm Γ(am+ b) , x ∈ R, a > 0, b > 0. (9) The following properties of Ma,b(x) hold [49]: Ma,b(x) = xMa,a+b(x) + 1 Γ(b) , (10) L [ tb−1Ma,b (±κta) ] = sa−b sa ∓ κ . (11) Proposition 1 (Green’s Formula [50]). Let D be a domain of Rn and v(x) its exterior normal. Then, for two regular functions v and w, Green’s formula is expressed as:∫ D (∆v)w dx = − ∫ D ∇v · ∇w dx+ ∫ ∂D ∂v ∂n w dσ, (12) where ∂D is the boundary of region D. We now introduce the following definition of the Lyapunov function and some lemmas [49–51] to prove the stability of equilibrium points. Definition 5 (Lyapunov Function). Let Θ(Ξ) be a neighborhood of Ξ. A real-valued differentiable function V defined on Θ(Ξ) is said to be a Lyapunov function for system ( 2) if: 1. V(Ξ) = 0, 2. V (Ξ1) > 0 in Θ(Ξ), ∀Ξ ̸= Ξ1. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 10 of 39 Lemma 1. Let Ω be a positively invariant subset of Π, where Π ⊆ Rn. If y : Π → R is continuously differentiable and y(x) > 0, and CDα t y(x(t)) ≤ 0 in Ω, then: • The set F contains all points in Ω where CDα t y(x(t)) = 0, • Σ is the largest invariant set in F , • Every bounded solution starting in Ω → Σ as t→ ∞. Lemma 2. Let y(t) ∈ R+ be a continuously differentiable function. For t ≥ 0 and α ∈ (0, 1): cDα t [ y∗Φ ( y(t) y∗ )] ≤ ( 1− y∗ y(t) ) cDα t y(t), y∗ ∈ R+, (13) where Φ(y) = (y − 1− ln(y)) ≥ 0 for any y > 0. 4. Qualitative Analysis of Model Dynamics In this section, we analyze the behavior of the proposed fractional-order model (2) to understand its general characteristics. These analyses allow us to establish the positivity and boundedness of the solutions, describe solution trends by evaluating equilibrium points (EPs), investigate their stability, and perform sensitivity analysis. 4.1. Existence, Positivity, and Boundedness of Solutions To ensure that the proposed fractional-order reaction-diffusion model (2) reflects bio- logical reality, we analyze the existence, positivity, and boundedness of its solutions. Let X = C(0̄,R) be a Banach space equipped with the usual norms. The system (2) can be rewritten in a compact form as:{ C 0 D α t λ(t, x)−Aλ(t, x) = F (t, x), λ(0, x) = λ0, (14) where λ = (S, V,E,Np, P,R,A) T , λ0 = (S0, V0, E0, Np,0, P0, R0, A0) T , and Aλ(t, x) = (d1∆S, d2∆V, d3∆E, d4∆Np, d5∆P, d6∆R, d7∆A) T . Here, A : D(A) ⊂ X7 → X7 is a linear diffusion operator with domain: D(A) = {λ ∈ X7 : ∆λ ∈ X7, ∂λ ∂n = 0R7 for x ∈ ∂0}. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 11 of 39 The nonlinear function F : [0, T ]×X7 → X7 is defined as: F (λ(t, x)) =  Λ + ρ1R− (λ+ δ + µ)S δS − (1− δ)λV − (γ + µ)V λS + (1− δ)λV − (α1 + α2 + ω + µ)E α1E − (σ + µ+ ω1)Np α2E − (σ + µ+ ω2)P γV + ωE + ω1Np + ω2P − (ρ+ ρ1 + µ)R ρR− (σ + µ)A  . Theorem 1. The problem (14) has a unique positive solution for all α ∈ (0, 1]. Proof. The operator A is linear and densely defined on D(A), and the function F is Lipschitz continuous in X. By the fixed-point theorem for fractional differential equations (see [52]), there exists a unique solution λ(t, x) for the problem (14). Furthermore, the positivity of the initial conditions λ0 ≥ 0 ensures that all components of the solution remain non-negative due to the structure of F , which reflects biological constraints. Hence, the solution is positive for all t > 0 and x ∈ 0. Theorem 2. The solution of the fractional-order system (2) is bounded on 0 × [0,+∞) for all t ≥ 0. Proof. The total population at time t is given by: N(t) = ∫ 0 [S(t, x) + V (t, x) + E(t, x) +Np(t, x) + P (t, x) +R(t, x) +A(t, x)] dx. (15) Adding the equations in (2) and integrating over 0, we have:∫ 0 7∑ i=1 C 0 D α t λi dx = ∫ 0 7∑ i=1 di∆λi dx+ ∫ 0 7∑ i=1 Fi dx. Using Green’s formula and the boundary conditions, we get:∫ 0 7∑ i=1 di∆λi dx = 0. Hence: C 0 D α t N(t) + ηN(t) ≤ Λ∥0∥. Taking the Laplace Transform of the above inequality, we obtain: sαL{N} − sα−1N(0) + ηL{N} ≤ Λ s . Solving for L{N}, we get: L{N} ≤ Λ s(sα + η) + sα−1N(0) sα + η . K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 12 of 39 Using the inverse Laplace Transform and the properties of the Mittag-Leffler function, we have: N(t) ≤ max { Λ η ,N(0) } [ηtαMα,α+1(−ηtα) +Mα,1(−ηtα)] . Since Mα,α+1 and Mα,1 are bounded for all t, it follows that: N(t) ≤ Λ η , ensuring boundedness of the solution. This completes the proof. 5. Equilibria of the Model Model (2) exhibits two types of equilibrium states. The first is the polio-free equilib- rium, expressed as: ε0 = ( Λk5k6 ξ , δΛk6 k1k5k7 − δγρ1 , 0, 0, 0, σγΛ k1k5k7 − δγρ1 , ρσγΛ k6(k1k5k7 − δγρ1) ) , along with a second equilibrium, referred to as the polio-endemic equilibrium Ξ∗. Utilizing the next-generation matrix approach as described by Shuai and van den Driessche [53], the transmission and transition matrices for model (2) at the equilibrium point E0 are determined as follows: F = βcr(S0+(1−δ)V0) N0 βcr(S0+(1−δ)V0) N0 βcr(S0+(1−δ)V0) N0 0 0 0 0 0 0  , where S0, V0, and N0 are the values of S, V , and N at E0, and V = −k2 0 0 α1 −k3 0 α2 0 −k4  . By applying the approach outlined in [53], the basic reproduction number is derived as the spectral radius of the next-generation matrix FV−1. Consequently, the basic reproduction number is expressed as: R0 = βcrk5 (k1 + δ(1− δ)) (α1k4 + α2k3 + k3k4) ηk2k3k4 , where η = k1k5 + δ ( k5 + γ + ργ k6 ) . k1 = γ + µ, k2 = α1 + α2 + µ+ ω, k3 = σ + µ+ ω1, k5 = ρ+ ρ1 + µ, k6 = σ + µ, k7 = δ + µ. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 13 of 39 The polio-endemic equilibrium is represented as Ξ∗ = ( S∗, V ∗, E∗, N∗ p , P ∗, R∗, A∗), with the components defined as follows: S∗ = Λk5 (k5 − Γ1ρ1) (λ∗ + k7) , V ∗ = k5δΛ (k5 − Γ1ρ1) (λ∗ + k7) ((1− δ)λ∗ + k1) , E∗ = λ∗k5Λ (k5 − Γ1ρ1) (λ∗ + k7) k2 ( 1 + δ(1− δ) (1− δ)λ∗ + k1 ) , N∗ p = λ∗k5Λα1 (k5 − Γ1ρ1) (λ∗ + k7) k2k3 ( 1 + δ(1− δ) (1− δ)λ∗ + k1 ) , P ∗ = λ∗k5Λα2 (k5 − Γ1ρ1) (λ∗ + k7) k2k4 ( 1 + δ(1− δ) (1− δ)λ∗ + k1 ) , A∗ = ρΓ1Λ k5 (k6 (λ∗ + k7)− ρ1Γ1) , R∗ = Γ1Λ k6 (λ∗ + k7)− ρ1Γ1 . Here, Γ1 is given by Γ1 = 1 (x∗ + k1) [ x∗ (ωk6k4 + α1k6 + α2k5) + δ (1− δ) (1− δ)x∗ + k1 + γ ((1− δ)x∗ + k1) ] . Substituting the appropriate variables into the expression for λ and simplifying shows that at equilibrium, x∗ satisfies the following equation: x∗ ( Φ2(x ∗)2 +Φ1x ∗ +Φ0 ) = 0, (16) where Φ2 = (1− δ) [ k5 − k6 k4 ( γ + δ k1 ) − ωk6 k4k3 (α1k4 + α2k5) ] , Φ1 = k4k5 (1− δ) [1− k1(k6(1− δ))]− ωk6k4k5 k3k4k5 (α1k4 + α2k5) , Φ0 = k4k5 [γ + ρR∗ − ρ1R ∗] . Equation (16) describes the equilibrium states of model (2), where the condition x∗ = 0 represents the polio-free equilibrium. The non-trivial solution x∗ = √ − Φ1 2Φ2 + √ Φ2 1 4Φ2 2 − Φ0 Φ2 represents the polio-endemic equilibrium. The number of endemic equilibria in model (1) is determined by the positive roots of Eq. (16). These positive roots depend on the signs of the coefficients Φ0, Φ1, and Φ2, and each positive root corresponds to an endemic K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 14 of 39 equilibrium point. It can be shown that Φ2 > 0, meaning the number of positive roots (and therefore the number of endemic equilibria) is influenced only by the signs of Φ0 and Φ1. Using Descartes’ rule of signs, the following conclusions can be drawn: Theorem 3. For model ( 2), the following conditions apply: 1. No endemic equilibria exist (i.e., only x∗ = 0 is valid) if Φ0 > 0. 2. A unique endemic equilibrium exists when Φ0 < 0. 3. Two endemic equilibria exist if Φ1 < 0 and Φ0 > 0. The next section focuses on analyzing the stability of the model’s equilibrium states. 6. Global Stability of the Endemic Equilibrium Theorem 4. The Endemic Equilibrium Ξ∗ is globally asymptotically stable (GAS) if R0 > 1. Proof. Assuming the Lyapunov function as follows: V2(t, x) = ∫ Ω [ S∗Φ ( S S∗ ) + E∗Φ ( E E∗ ) + η + δ δ N∗ pΦ ( Np N∗ p ) + . . . ] dx, where Φ(z) = z − 1− ln(z) satisfies Φ(z) ≥ 0, with equality only if z = 1. Using Lemma 2 and the system of equations, we calculate: C 0 D α t V2(t, x) = ∫ Ω [ C 0 D α t ( S∗Φ ( S S∗ )) + C 0 D α t ( E∗Φ ( E E∗ )) + η + δ δ C 0 D α t ( N∗ pΦ ( Np N∗ p )) + . . . ] dx. Using the equation for S: C 0 D α t S = Λ+ ρ1R− (λ+ δ + µ)S − d1∆S, we compute: C 0 D α t ( S∗Φ ( S S∗ )) = S∗ ( 1− S∗ S ) [Λ + ρ1R− (λ+ δ + µ)S − d1∆S] . Expanding and simplifying, we get: C 0 D α t ( S∗Φ ( S S∗ )) ≤ ∫ Ω [ d1∆S − d1S ∗∆S S + source terms ] dx. Using the equation for E: C 0 D α t E = λS + (1− δ)λV − (α1 + α2 + ω + µ)E − d3∆E, K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 15 of 39 we compute: C 0 D α t ( E∗Φ ( E E∗ )) = E∗ ( 1− E∗ E ) [λS + (1− δ)λV − (α1 + α2 + ω + µ)E − d3∆E] . This leads to: C 0 D α t ( E∗Φ ( E E∗ )) ≤ ∫ Ω [ d3∆E − d3E ∗∆E E + source terms ] dx. Using the equation for Np: C 0 D α t Np = α1E − (σ + µ+ ω1)Np − d4∆Np, we compute: C 0 D α t ( N∗ pΦ ( Np N∗ p )) = N∗ p ( 1− N∗ p Np ) [α1E − (σ + µ+ ω1)Np − d4∆Np] . Simplifying: C 0 D α t ( N∗ pΦ ( Np N∗ p )) ≤ ∫ Ω [ d4∆Np − d4N ∗ p ∆Np Np + source terms ] dx. Adding the contributions of all terms: C 0 D α t V2(t, x) ≤ − ∫ Ω [ d1S ∗ |∇S|2 S2 + d3E ∗ |∇E|2 E2 + d4N ∗ p |∇Np|2 N2 p + . . . ] dx. The terms involving (1−X∗/X) simplify using equilibrium relations:∫ Ω [ΛS∗ − ρ1R ∗ + . . . ] dx, which are all non-positive due to the equilibrium relations and the properties of Φ(z). Apply Fractional LaSalle Invariance Principle Since C 0 D α t V2(t, x) ≤ 0 and equality holds only when S = S∗, E = E∗, . . ., the largest invariant set is {Ξ∗}. By the fractional LaSalle Invariance Principle, Ξ∗ is GAS when R0 > 1. 7. Global Stability of the Disease-Free Equilibrium To analyze the global stability of the Disease-Free Equilibrium (DFE) of the model at ε0, we proceed with the following steps. The Disease-Free Equilibrium of the model is given by: ε0 = (S∗, V ∗, 0, 0, 0, R∗, A∗) , where: S∗ = Λk5k6 ξ , V ∗ = δΛk6 ξ , R∗ = σγΛ ξ , A∗ = ρσγΛ k6ξ , and ξ = k1k5k7 − δγρ1. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 16 of 39 Theorem 5. The Disease-Free Equilibrium ε0 of the system is globally asymptotically stable if R0 < 1. Proof. We define the Lyapunov function: V(t, x) = ∫ Ω [ S∗Φ ( S S∗ ) + E∗Φ ( E E∗ ) +N∗ pΦ ( Np N∗ p ) + P ∗Φ ( P P ∗ )] dx, where Φ(z) = z − 1− ln(z), which satisfies Φ(z) ≥ 0 and is equal to 0 only if z = 1. Here, E∗ = N∗ p = P ∗ = 0 at the DFE (ε0). Time Derivative of V, using the system of equations, compute the fractional derivative of V: C 0 D α t V(t, x) = ∫ Ω [ C 0 D α t ( S∗Φ ( S S∗ )) +C 0 D α t ( E∗Φ ( E E∗ )) +C 0 D α t ( N∗ pΦ ( Np N∗ p )) +C 0 D α t ( P ∗Φ ( P P ∗ ))] dx. Evaluate Each Compartment. The equation for S is: C 0 D α t S = Λ+ ρ1R− (λ+ δ + µ)S − d1∆S. Substitute into the Lyapunov derivative: C 0 D α t ( S∗Φ ( S S∗ )) = S∗ ( 1− S∗ S )[ Λ + ρ1R− (λ+ δ + µ)S − d1∆S ] . At the DFE (S = S∗): C 0 D α t ( S∗Φ ( S S∗ )) ≤ −d1S∗ |∇S|2 S2 . The equation for E is: C 0 D α t E = λS + (1− δ)λV − (α1 + α2 + ω + µ)E − d3∆E. Substitute into the Lyapunov derivative: C 0 D α t ( E∗Φ ( E E∗ )) = E∗ ( 1− E∗ E )[ λS + (1− δ)λV − (α1 + α2 + ω + µ)E − d3∆E ] . At the DFE (E = 0): C 0 D α t ( E∗Φ ( E E∗ )) ≤ −d3E∗ |∇E|2 E2 . The equation for Np is: C 0 D α t Np = α1E − (σ + µ+ ω1)Np − d4∆Np. Substitute into the Lyapunov derivative: C 0 D α t ( N∗ pΦ ( Np N∗ p )) = N∗ p ( 1− N∗ p Np )[ α1E − (σ + µ+ ω1)Np − d4∆Np ] . K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 17 of 39 At the DFE (Np = 0): C 0 D α t ( N∗ pΦ ( Np N∗ p )) ≤ −d4N∗ p |∇Np|2 N2 p . The equation for P is: C 0 D α t P = α2E − (σ + µ+ ω2)P − d5∆P. Substitute into the Lyapunov derivative: C 0 D α t ( P ∗Φ ( P P ∗ )) = P ∗ ( 1− P ∗ P )[ α2E − (σ + µ+ ω2)P − d5∆P ] . At the DFE (P = 0): C 0 D α t ( P ∗Φ ( P P ∗ )) ≤ −d5P ∗ |∇P |2 P 2 . Combine the contributions of all compartments: C 0 D α t V(t, x) ≤ − ∫ Ω [ d1S ∗ |∇S|2 S2 + d3E ∗ |∇E|2 E2 + d4N ∗ p |∇Np|2 N2 p + d5P ∗ |∇P |2 P 2 ] dx. All terms are non-positive, with equality only when S = S∗, E = 0, Np = 0, and P = 0. The basic reproduction number R0 is derived using the next-generation matrix method: R0 = βcrk5 (k1 + δ(1− δ)) (α1k4 + α2k3 + k3k4) ηk2k3k4 , where: η = k1k5 + δ ( k5 + γ + ργ k6 ) . If R0 < 1, the disease dies out, and the Lyapunov derivative is strictly negative. By the fractional LaSalle Invariance Principle, the DFE is globally asymptotically stable when R0 < 1. 8. Sensitivity Analysis Determining the parameters that help reduce the spread of infectious diseases is car- ried out through sensitivity analysis. Forward sensitivity analysis is considered a vital component of disease modeling, though its computation becomes tedious for complex bio- logical models. Sensitivity analysis of R0 has received much attention from ecologists and epidemiologists. Definition 6. The normalized forward sensitivity index of R0 that depends differentiably on a parameter Ω is defined as K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 18 of 39 SΩ = Ω R0 ∂R0 ∂Ω . Three methods are commonly used to calculate sensitivity indices: (i) by direct differenti- ation, (ii) by a Latin hypercube sampling method, and (iii) by linearizing the system (2) and then solving the obtained set of linear algebraic equations. We will apply the direct differentiation method as it provides analytical expressions for the indices. The indices not only show us the influence of various aspects associated with the spread of infectious diseases but also provide important information regarding the comparative change be- tween R0 and different parameters. Consequently, it helps in developing effective control strategies. Table 1 shows that the parameters β, c, r, ρ1, γ, and γ1 positively influence the reproduction number R0. This implies that an increase or decrease in these parameters by 10% will proportionally increase or decrease R0 by 10%, 10%, 10%, 2.2384%, 7.283%, and 9.347%, respectively. On the other hand, the indices for parameters ρ, µ, δ, σ, ω, ω1, and ω2 show that increasing their values by 10% will decrease R0 by 0.25475%, 0.46834%, 5.4797%, 7.725%, 0.10323%, 0.034112%, and 0.0031115%, respectively. Parameters α1 and α2 also have a negative impact, with R0 decreasing by 0.56475% and 0.48046%, respectively, for a 10% increase in these parameters. Figs. 3(a, b, c, d, e, f) and Figs. 4(a, b, c, d, e, f) depict the sensitivity of various parameters. These results underline the importance of prioritizing interventions targeting the most sensitive parameters to effectively control the disease. Parameter SIndex Value Parameter SIndex Value β Sβ 1.00 c Sc 1.00 r Sr 1.00 ρ Sρ -0.025475 ρ1 Sρ1 0.22384 µ Sµ -0.046834 γ Sγ -0.084718 δ Sδ -0.54797 σ Sσ -0.7725 ω Sω -0.010323 ω1 Sω1 -0.0034112 ω2 Sω2 -0.00031115 α1 Sα1 -0.056475 α2 Sα2 -0.048046 Table 1: Sensitivity indices of the reproduction number R0 against mentioned parameters. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 19 of 39 Sensitivity Analysis of R 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 c r 1 1 2 1 2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Figure 2: Sensitivity bar chart related to model parameters. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 20 of 39 (a) R0 vs σ and µ Contour Plot of R0 for vs 0.0026639 0.0028231 0.0028231 0.0029823 0.0029823 0.0031415 0.0031415 0.0033007 0.0033007 0.0034599 0.0034599 0.00361910.0037783 0.1 0.12 0.14 0.16 0.18 0.2 0.05 0.06 0.07 0.08 0.09 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (b) R0 vs σ and µ (c) R0 vs σ and β Contour Plot of R0 for vs 1. 64 38 1. 64 38 2. 18 26 2. 18 26 2. 72 14 2. 72 14 3. 26 02 3. 26 02 3. 79 9 3. 79 9 4. 33 78 4. 33 78 4. 87 66 4. 87 66 5. 41 53 5. 95 41 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.002 0.004 0.006 0.008 0.01 1 2 3 4 5 6 7 8 (d) R0 vs σ and β (e) R0 vs µ and β Contour Plot of R0 for vs 0. 08 59 15 0. 08 59 15 0. 17 02 7 0. 17 02 7 0. 25 46 2 0. 25 46 2 0. 33 89 8 0. 33 89 8 0. 42 33 3 0. 42 33 3 0. 50 76 9 0. 50 76 9 0. 59 20 4 0. 59 20 4 0. 67 63 9 0. 76 07 5 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.02 0.04 0.06 0.08 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (f) R0 vs µ and β Figure 3: 3D sensitivity analysis profiles showing the impact of various parameter pairs (σ, β, µ) on the basic reproduction number R0. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 21 of 39 (a) R0 vs ω1 and α2 Contour Plot of R0 for 1 vs 2 0.002902 0. 00 31 04 9 0.0031049 0. 00 33 07 8 0.0033078 0. 00 35 10 7 0. 00 35 10 7 0. 00 37 13 5 0. 00 37 13 5 0. 00 39 16 4 0. 00 39 16 4 0.0041193 0.0041193 0.0043222 0.0043222 0.0045251 0.1 0.2 0.3 0.4 0.5 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 2 0 0.2 0.4 0.6 0.8 1 (b) R0 vs ω1 and α2 (c) R0 vs µ and α1 Contour Plot of R0 for 1 vs 0.00407750.004185 0.0042925 0.0042925 0.0044 0.0044 0.0045075 0.0045075 0.004615 0.004615 0.00472260.0048301 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1 0.005 0.01 0.015 0.02 0.025 0.03 0 0.2 0.4 0.6 0.8 1 (d) R0 vs µ and α1 (e) R0 vs α2 and α1 Contour Plot of R0 for 2 vs 1 0.004333 0.004364 0.004364 0.004395 0.004395 0.0044261 0.0044261 0.0044571 0.0044571 0.00448810.0045192 0.4 0.5 0.6 0.7 0.8 0.9 1 2 0.4 0.5 0.6 0.7 0.8 0.9 1 1 0 0.2 0.4 0.6 0.8 1 (f) R0 vs α2 and α1 Figure 4: 3D sensitivity analysis profiles showing the impact of various parameter pairs (ω1, α1, α2) on the basic reproduction number R0. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 22 of 39 9. Numerical Scheme In this section, we present the numerical simulations of the fractional-order reaction diffusion model (2). The fractional-order derivatives are approximated using a forward finite difference scheme, while the diffusion operator is discretized using a centered finite difference approach. The domain is given by U = [0, L], with temporal and spatial step sizes ht = T N and hx = L n , respectively. Here, ti = iht for i = 0, 1, . . . , N and xj = jhx for j = 0, 1, . . . , n. For simplicity, the approximate solutions for S(ti, xj), V (ti, xj), E(ti, xj), Np(ti, xj), P (ti, xj), R(ti, xj), and A(ti, xj) are denoted as Si j , V i j , E i j , N i p,j , P i j , R i j , and A i j , respectively. We can write, for example, the approximation of C 0 D α t S(ti, xj) and ∆S(ti, xj), respectively, as follows: C 0 D α t S(ti, xj) ≈ 1 Γ(2− α) i∑ l=0 (l + 1)1−α − l1−α hαt (Si+1−l j − Si−l j ), and ∆S(ti, xj) ≈ Si j+1 − 2Si j + Si j−1 h2x . Using the discretization, we derive the following numerical schemes for each compartment of the model: Si+1 j = Si j − i∑ l=1 (l + 1)1−α − l1−α Γ(2− α)hαt (Si+1−l j − Si−l j ) + dS Γ(2− α)hαt h2x (Si j+1 − 2Si j + Si j−1) +Γ(2− α)hαt [ Λ + ρ1R i j − (λ+ δ + µ)Si j ] , V i+1 j = V i j − i∑ l=1 (l + 1)1−α − l1−α Γ(2− α)hαt (V i+1−l j − V i−l j ) + dV Γ(2− α)hαt h2x (V i j+1 − 2V i j + V i j−1) +Γ(2− α)hαt [ δSi j − (1− δ)λV i j − (γ + µ)V i j ] , Ei+1 j = Ei j − i∑ l=1 (l + 1)1−α − l1−α Γ(2− α)hαt (Ei+1−l j − Ei−l j ) + dE Γ(2− α)hαt h2x (Ei j+1 − 2Ei j + Ei j−1) +Γ(2− α)hαt [ λSi j + (1− δ)λV i j − (α1 + α2 + ω + µ)Ei j ] , N i+1 p,j = N i p,j− i∑ l=1 (l + 1)1−α − l1−α Γ(2− α)hαt (N i+1−l p,j −N i−l p,j )+dNp Γ(2− α)hαt h2x (N i p,j+1−2N i p,j+N i p,j−1) +Γ(2− α)hαt [ α1E i j − (σ + µ+ ω1)N i p,j ] , P i+1 j = P i j − i∑ l=1 (l + 1)1−α − l1−α Γ(2− α)hαt (P i+1−l j − P i−l j ) + dP Γ(2− α)hαt h2x (P i j+1 − 2P i j + P i j−1) K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 23 of 39 +Γ(2− α)hαt [ α2E i j − (σ + µ+ ω2)P i j ] , Ri+1 j = Ri j − i∑ l=1 (l + 1)1−α − l1−α Γ(2− α)hαt (Ri+1−l j −Ri−l j ) + dR Γ(2− α)hαt h2x (Ri j+1 − 2Ri j +Ri j−1) +Γ(2− α)hαt [ γV i j + ωEi j + ω1N i p,j + ω2P i j − (ρ+ ρ1 + µ)Ri j ] , Ai+1 j = Ai j − i∑ l=1 (l + 1)1−α − l1−α Γ(2− α)hαt (Ai+1−l j −Ai−l j ) + dA Γ(2− α)hαt h2x (Ai j+1 − 2Ai j +Ai j−1) +Γ(2− α)hαt [ ρRi j − (σ + µ)Ai j ] . This numerical scheme allows for a detailed investigation of the spatiotemporal behavior of the model under various parameter configurations and initial conditions. 10. Numerical Results and Discussion In this section, we analyze the numerical results obtained from the fractional-order reaction-diffusion model. The simulations are conducted using the parameter values listed in Table ??. For the spatial domain, we consider the one-dimensional interval 0 ≤ x ≤ L, while for the temporal domain, we have 0 ≤ t ≤ T . The diffusion coefficients used in the model are specified as follows: DS = 0.01, DV = 0.2, DE = 0.2, DNp = 0.2, DP = 0.01, DR = 0.2, and DA = 0.1 (all in km per day). The numerical solutions for the com- partments S, V , E, Np, P , R, and A demonstrate complex spatiotemporal behaviors under varying parameter configurations. These solutions reveal the interactions among susceptible, exposed, vaccinated, and infected individuals over time and space. The choice of diffusion coefficients significantly influences the spatial spread of the disease. Higher diffusion coefficients, such as DE , DNp , and DR, result in faster spatial propagation of the disease, leading to more homogeneous distributions of infected individuals across the domain. Conversely, smaller diffusion coefficients, such as DS and DP , create localized variations, maintaining higher concentrations of susceptible and infected individuals in specific regions. The recovery rates (ω, ω1, ω2) and the of progression rates (α1, α2) decide the elements of changes between various compartments. Higher recovery rates (ω1, ω2) decline the span people spend in the infected class, in this way decreasing the spread of the illness. Interest- ingly, higher movement rates (α1, α2) speed up the development of uncovered people into tainted compartments, prompting an expanded infection trouble. These elements highlight the significance of recuperation and movement boundaries in forming the overall disease trajectory. To approve the model’s precision, we simulated different situations utilizing different initial conditions and parameter values. The outcomes, illustrated in Figures, portray the temporal evolution of the compartments over the time interval 0 ≤ t ≤ T . The simulations result show that the model effectively catches fundamental elements of illness elements, including beginning dramatic development, top contamination levels, and K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 24 of 39 inevitable decay because of recuperation and inoculation endeavors. The transaction of boundary values, for example, immunization rate (δ) and disappearing resistance rate (ρ1), was basic in deciding the elements of the powerless and recuperated populaces over time. The results of our simulations propose a few significant experiences. Expanding inocu- lation rates (δ) significantly reduces the quantity of susceptible people, which by implica- tion brings down the disease rate. Moreover, the waning immunity rate (ρ1) assumes a vital part in deciding the drawn out dynamics of the recovered compartment, as higher upsides of ρ1 lead to quicker changes of recuperated people back to the powerless class. Moreover, the boundary c, addressing the typical pace of contact with waste, firmly impacts the underlying spread of the illness. Reducing c through improved sanitation measures can ef- fectively control the outbreak, highlighting the importance of environmental interventions in disease management. Table 2: Descriptions of model parameters and their values used for simulation. Parameter Description Value Ref. µ Natural death rate 8.75× 10−3 [1] σ Polio-induced death rate 0.125 [1] β Transmission probability 0.002 [1] ρ Rate of development of PPS 0.030 [1] Λ Recruitment rate 1000 [1] ω Recovery rate of exposed individuals 0.010 [1] ρ1 Waning rate of post-recovery immunity 0.330 [1] γ Recovery rate of vaccinated individuals 0.140 [1] δ Vaccination rate of susceptibles 0.378 [1] ω2 Recovery rate in paralytic class 10−4 [1] ω1 Recovery rate in non-paralytic class 0.001 [1] α2 Progression rate to paralytic class 0.450 [1] α1 Progression rate to non-paralytic class 0.500 [1] c Contact rate with faecal waste 0.500 [1] K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 25 of 39 (a) 0 50 100 150 200 250 300 350 400 Time in days 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 S (t ) - S us ce pt ib le P op ul at io n 109 Susceptible Population over Time for Different = 0.7 = 0.75 = 0.8 = 0.85 = 0.9 (b) Figure 5: Simulation results of the susceptible population for different fractional orders α. (a) 0 100 200 300 400 500 600 700 800 Time in days 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 S (t ) - S us ce pt ib le P op ul at io n 109 Susceptible Population over Time for Different = 0.8 = 0.85 = 0.9 = 0.95 = 1 (b) Figure 6: Long term behavior of susceptible population for different fractional orders α. The simulation results in Figure 5 illustrate the impact of fractional-order derivatives and spatial heterogeneity on the dynamics of the susceptible population S(t, x) in the fractional-order reaction-diffusion polio model. The fractional order α introduces memory effects, influencing the depletion rate of susceptibles. Lower values of α (e.g., α = 0.7) correspond to stronger memory effects, accelerating the depletion due to the cumulative influence of historical interactions. Higher α values (e.g., α = 0.9) result in a slower de- pletion, reflecting behavior closer to classical models with limited memory effects. The three-dimensional plot shows significant spatial heterogeneity, where regions with higher K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 26 of 39 initial densities of susceptibles experience faster depletion due to localized transmission dy- namics. These findings underscore the importance of spatially targeted interventions, such as vaccination campaigns and improved sanitation, particularly in high-density regions. The interplay of diffusion and memory effects highlights the need to account for both mo- bility and historical interactions in control strategies, emphasizing the utility of fractional models in capturing the complex temporal and spatial dynamics of polio outbreaks. (a) 0 50 100 150 200 250 300 350 400 Time in days 0 0.5 1 1.5 2 2.5 3 3.5 V (t ) - V a c c in a te d P o p u la ti o n 108 Vaccinated Population over Time for Different = 0.7 = 0.75 = 0.8 = 0.85 = 0.9 (b) Figure 7: Simulation results of the vaccinated population for different fractional orders α. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 27 of 39 (a) 0 100 200 300 400 500 600 700 800 Time in days 0 0.5 1 1.5 2 2.5 3 3.5 4 V (t ) - V a c c in a te d P o p u la ti o n 108 Vaccinated Population over Time for Different = 0.8 = 0.85 = 0.9 = 0.95 = 1 (b) Figure 8: Long term behavior of vaccinated population for different fractional orders α. The dynamics of the vaccinated population V (t, x), depicted in Figure 7, are signifi- cantly influenced by fractional-order memory effects and spatial heterogeneity. Lower α values (e.g., α = 0.7) lead to a prior and more keen top in the vaccinated individuals, reflecting quick take-up driven by more grounded memory impacts. Conversely, higher α values (e.g., α = 0.9) bring about deferred and continuous immunization elements, similar to old style models. The three-layered plot uncovers spatial changeability in im- munization inclusion, with longer times expected to accomplish immersion in certain areas. This stresses the requirement for geologically designated endeavors to accomplish uniform immunization inclusion. The cooperation among dissemination and memory impacts pro- poses that areas with high populace thickness or more prominent openness to the infection accomplish quicker immunization tops however may likewise encounter more fast down- falls because of fading resistance or antibody disappointment. These outcomes highlight the worth of fragmentary request models in improving vaccination methodologies for polio control in heterogeneous populations. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 28 of 39 (a) 0 50 100 150 200 250 300 350 400 Time in days 0 2 4 6 8 10 12 14 16 18 E (t ) - E x p o s e d P o p u la ti o n 107 Exposed Population over Time for Different = 0.7 = 0.75 = 0.8 = 0.85 = 0.9 (b) Figure 9: Simulation results of the exposed population for different fractional orders α. (a) 0 100 200 300 400 500 600 700 800 Time in days 0 2 4 6 8 10 12 14 16 18 E (t ) - E x p o s e d P o p u la ti o n 107 Exposed Population over Time for Different = 0.8 = 0.85 = 0.9 = 0.95 = 1 (b) Figure 10: Long term behavior of exposed population for different fractional orders α. Figure 10 illustrates the dynamics of the exposed population E(t, x), reflecting the influence of memory effects and spatial heterogeneity on disease progression. Smaller K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 29 of 39 α values (e.g., α = 0.7) result in earlier and sharper peaks, driven by enhanced mem- ory effects, which amplify historical exposure and infection dynamics. Larger α values (e.g., α = 0.9) lead to slower and later peaks, reflecting weaker memory effects. Spatial heterogeneity, shown in the three-dimensional plot, highlights areas with concentrated transmission potential. Diffusion redistributes exposed individuals, amplifying localized outbreaks in regions with high susceptibility or insufficient vaccination. These results em- phasize the importance of region-specific interventions to mitigate exposure and control disease progression, showcasing the utility of fractional-order models in designing effective strategies. (a) 0 50 100 150 200 250 300 350 400 Time in days 0 1 2 3 4 5 6 7 8 9 N p (t ) - In fe c te d n o n -p a ra ly ti c P o p u la ti o n 107Infected non-paralytic Population over Time for Different = 0.7 = 0.75 = 0.8 = 0.85 = 0.9 (b) Figure 11: Simulation results of the infected non-paralytic population for different frac- tional orders α. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 30 of 39 (a) 0 100 200 300 400 500 600 700 800 Time in days -1 0 1 2 3 4 5 6 7 8 9 N p (t ) - In fe c te d n o n -p a ra ly ti c P o p u la ti o n 107Infected non-paralytic Population over Time for Different = 0.8 = 0.85 = 0.9 = 0.95 = 1 (b) Figure 12: Long term behavior of infected non-paralytic population for different fractional orders α. The infected non-paralytic population Np(t, x), shown in Figure 12, highlights the effects of fractional orders and spatial factors on disease dynamics. Smaller α values (e.g., α = 0.7) produce earlier and sharper peaks, reflecting accelerated disease progression due to memory effects. Conversely, larger α values (e.g., α = 0.9) lead to delayed peaks and gradual transitions. The spatial plot underscores localized transmission influenced by high population densities and environmental factors, necessitating targeted public health measures. Fractional-order models provide critical insights into managing non-paralytic infections and mitigating transmission risks. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 31 of 39 (a) 0 50 100 150 200 250 300 350 400 Time in days 0 0.5 1 1.5 2 2.5 3 3.5 4 P (t ) - In fe c te d P a ra ly ti c P o p u la ti o n 107 Infected Paralytic Population over Time for Different = 0.7 = 0.75 = 0.8 = 0.85 = 0.9 (b) Figure 13: Simulation results of the infected paralytic population for different fractional orders α. (a) 0 100 200 300 400 500 600 700 800 Time in days -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 P (t ) - In fe c te d P a ra ly ti c P o p u la ti o n 107 Infected Paralytic Population over Time for Different = 0.8 = 0.85 = 0.9 = 0.95 = 1 (b) Figure 14: Long term behavior of infected paralytic population for different fractional orders α. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 32 of 39 The dynamics of the infected paralytic population P (t, x), shown in Figure 14, reveal critical insights into severe polio manifestations. Smaller α values (e.g., α = 0.7) yield faster progression to paralysis, while larger α values (e.g., α = 0.9) result in delayed but more prolonged peaks. Spatial heterogeneity highlights regions with concentrated severe cases, driven by diffusion and clustering effects. These results emphasize the need for targeted medical interventions and enhanced sanitation in high-incidence regions, demon- strating the model’s value in understanding paralytic polio dynamics. (a) 0 50 100 150 200 250 300 350 400 Time in days 0 5 10 15 R (t ) - R e c o v e re d P o p u la ti o n 108 Recovered Population over Time for Different = 0.7 = 0.75 = 0.8 = 0.85 = 0.9 (b) Figure 15: Simulation results of the recovered population for different fractional orders α. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 33 of 39 (a) 0 100 200 300 400 500 600 700 800 Time in days 0 2 4 6 8 10 12 14 16 18 R (t ) - R e c o v e re d P o p u la ti o n 108 Recovered Population over Time for Different = 0.8 = 0.85 = 0.9 = 0.95 = 1 (b) Figure 16: Long term behavior of recovered population for different fractional orders α. The recovered population R(t, x), depicted in Figure 16, increases consistently over time across all α values. Smaller α values (e.g., α = 0.7) result in slower recovery rates, while larger α values (e.g., α = 0.9) show rapid recovery and earlier saturation. Spatial analysis highlights regional disparities influenced by vaccination and intervention efforts. These results guide the allocation of resources to optimize recovery outcomes and reduce disease burdens. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 34 of 39 (a) 0 50 100 150 200 250 300 350 400 Time in days 0 0.5 1 1.5 2 2.5 A (t ) - P o s t- p a ra ly ti c o r D is a b le d P o p u la ti o n 107 Post-paralytic Population over Time for Different = 0.7 = 0.75 = 0.8 = 0.85 = 0.9 (b) Figure 17: Simulation results of the post-paralytic population for different fractional orders α. (a) 0 100 200 300 400 500 600 700 800 Time in days 0 2 4 6 8 10 12 A (t ) - P o s t- p a ra ly ti c o r D is a b le d P o p u la ti o n 107 Post-paralytic Population over Time for Different = 0.8 = 0.85 = 0.9 = 0.95 = 1 (b) Figure 18: Long term behavior of post-paralytic population for different fractional orders α. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 35 of 39 Figure 18 highlights the dynamics of the post-paralytic population A(t, x), representing individuals with long-term effects of polio. Smaller α values (e.g., α = 0.7) lead to slower growth, while larger values (e.g., α = 0.9) show rapid accumulation. Spatial heterogeneity emphasizes regions with higher densities of disabled individuals, necessitating targeted healthcare and rehabilitation efforts. These results underscore the importance of fractional models in addressing chronic disease burdens and guiding long-term interventions. 11. Conclusion This study presents a novel fractional reaction-diffusion model tailored to the spatio- temporal dynamics of poliovirus transmission, integrating the effects of vaccination and post-paralytic outcomes. Unlike traditional models, which often assume spatial homo- geneity and classical derivatives, this research incorporates Caputo fractional derivatives to model memory-dependent disease progression and heterogeneous diffusion terms to cap- ture population mobility and environmental influences. The model divides the population into seven distinct compartments, addressing the intricate transitions from susceptibility to post-paralytic disability. Rigorous mathematical analysis demonstrates the existence, positivity, and boundedness of solutions, ensuring the biological realism of the model. Stability assessments of equilibrium points, using tools such as Lyapunov functions and reproduction number (R0), reveal conditions for both the eradication and persistence of poliovirus. Sensitivity analysis identifies critical factors, such as the fractional order (α) and vaccination rates, that strongly influence disease dynamics. Numerical simulations provide practical insights into how fractional memory effects alter the pace of disease spread and recovery. Lower values of α correlate with accelerated depletion of susceptibles and earlier peaks in infected populations, underscoring the importance of incorporating memory effects in outbreak predictions. Furthermore, the spatial distribution of disease burden, as revealed by the reaction-diffusion framework, highlights the necessity of region- specific public health interventions, particularly in high-density or sanitation-compromised areas. This study’s fractional reaction-diffusion approach significantly enhances the pre- dictive power of mathematical models in epidemiology, offering actionable insights into poliovirus control. By addressing both spatial and temporal complexities, it provides a robust tool for optimizing vaccination campaigns, improving sanitation strategies, and mitigating long-term disability outcomes, paving the way for more effective public health responses to poliovirus outbreaks. Acknowledgements The authors extend their appreciation to the King Salman center For Disability Re- search for funding this work through Research Group no KSRG-2024-200. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 36 of 39 Funding The authors extend their appreciation to the King Salman center For Disability Re- search for funding this work through Research Group no KSRG-2024-200. Conflict of interest The authors declare that they have no conflict of interest. Data Availability Statement “All data generated or analyzed during this study are included in this article”. References [1] C. S. Bornaa, B. Seidu, and O. D. Makinde. Mathematical analysis of the impact of vaccination and poor sanitation on the dynamics of poliomyelitis. International Journal of Nonlinear Sciences and Numerical Simulation, 24(1):161–169, 2023. [2] J. Radboud D. Tebbens, M. A. Pallansch, et al. A dynamic model of poliomyelitis outbreaks: learning from the past to help inform the future. American Journal of Epidemiology, 162(4):223–241, 2005. [3] M. Agarwal and A. S. Bhadauria. Modeling spread of polio with the role of vaccina- tion. Applications and Applied Mathematics: An International Journal, 6(2):552–571, 2011. [4] Y. Sabbar and A. A. Raezah. The influence of independent jumps on the dynamics of a perturbed SIRS epidemic model with altered behavior. International Journal of Dynamics and Control, 13(1):1–16, 2025. [5] Y. Sabbar. Exploring threshold dynamics of a behavioral epidemic model featuring two susceptible classes and second-order jump-diffusion. Chaos, Solitons & Fractals, 186:115216, 2024. [6] A. Din, Y. Sabbar, and P. Wu. A novel stochastic Hepatitis B virus epidemic model with second-order multiplicative α-stable noise and real data. Acta Mathematica Scientia, 44(2):752–788, 2024. [7] Y. Sabbar and A. A. Raezah. Modeling mosquito-borne disease dynamics via stochas- tic differential equations and generalized tempered stable distribution. AIMS Math- ematics, 9(8):22454–22485, 2024. [8] K. S. Nisar and Y. Sabbar. Long-run analysis of a perturbed HIV/AIDS model with antiretroviral therapy and heavy-tailed increments performed by tempered stable Lévy jumps. Alexandria Engineering Journal, 78:498–516, 2023. [9] A. Khan, Y. Sabbar, and A. Din. Stochastic modeling of the Monkeypox 2022 epi- demic with cross-infection hypothesis in a highly disturbed environment. Mathemat- ical Biosciences and Engineering, 19(12):13560–13581, 2022. K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 37 of 39 [10] K. Dietz and J. A. P. Heesterbeek. Daniel Bernoulli’s epidemiological model revisited. Mathematical Biosciences, 180(1-2):1–21, 2002. [11] Fred Brauer. Compartmental models in epidemiology. Mathematical Epidemiology, 1945:19–79, 2008. [12] I. Nali and A. Dénes. Global dynamics of a within-host model for Usutu virus. Computation, 11(11):226, 2023. [13] E. Duque-Maŕın, J. G. Vergaño-Salazar, I. Duarte-Gandica, and K. Vilches. Mathe- matical modelling of some poliomyelitis vaccination and migration scenarios in Colom- bia. Journal of Physics: Conference Series, 1160:012021, 2019. [14] C. J. Browne, R. J. Smith, and L. Bourouiba. From regional pulse vaccination to global disease eradication: insights from a mathematical model of poliomyelitis. Jour- nal of Mathematical Biology, 71(1):215–253, 2015. [15] A. Dénes and L. Székely. Global dynamics of a mathematical model for the possible re-emergence of polio. Mathematical Biosciences, 293:64–74, 2017. [16] F. A. Alrawajeh, F. M. Allehiany, A. Raza, S. A. M. Abdelmohsen, T. N. Cheema, M. Rafiq, and M. Mohsin. Bio-inspired computational methods for the polio virus epidemic model. Computers, Materials & Continua, 72(2):2351–2374, 2022. [17] A. Raza, D. Baleanu, Z. U. Khan, M. Mohsin, N. Ahmad, M. Rafiq, and P. Anwar. Stochastic analysis for the dynamics of a poliovirus epidemic model, 2023. Preprint, no journal specified. [18] M. Agarwal and A. S. Bhadauria. Modeling spread of polio with the role of vaccina- tion. Applications and Applied Mathematics: An International Journal, 6(2):552–571, 2011. [19] H. Miao, Z. Teng, X. Abdurahman, and Z. Li. Global stability of a diffusive and delayed virus infection model with general incidence function and adaptive immune response. Journal of Computational and Applied Mathematics, 37(12):3780–3805, 2018. [20] N. Ahmed, J. E. Macias-Diaz, N. Shahid, A. Raza, and M. Rafiq. On the computa- tional simulation of a temporally non-local and nonlinear diffusive epidemic model of disease transmission. International Journal of Modern Physics C, 2024. [21] E. Fadhal, M. M. Al-Shamiri, M. W. Yasin, S. M. H. Ashfaq, N. Ahmed, A. Raza, and M. Rafiq. Spatio-temporal dynamics of the vector-borne plant disease model. Modeling Earth Systems and Environment, 10(6):6417–6430, 2024. [22] M. W. Yasin, N. Ahmed, M. S. Iqbal, A. Raza, M. Rafiq, E. M. T. Eldin, and I. Khan. Spatio-temporal numerical modeling of stochastic predator-prey model. Scientific Reports, 13(1):1990, 2023. [23] R. Zarin. A robust study of dual variants of SARS-CoV-2 using a reaction-diffusion mathematical model with real data from the USA. The European Physical Journal Plus, 138(11):1011, 2023. [24] X. Ren, Y. Tian, L. Liu, and X. Liu. A reaction–diffusion within-host HIV model with cell-to-cell transmission. Journal of Mathematical Biology, 76(7):1831–1872, 2018. [25] S. Pankavich and C. Parkinson. Mathematical analysis of an in-host model of viral dynamics with spatial heterogeneity. Discrete and Continuous Dynamical Systems - K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 38 of 39 Series B, 21(4):1237–1257, 2016. [26] Y. Cai, Z. Ding, and W. Wang. Transmission dynamics of Zika virus with spatial structure—a case study in Rio de Janeiro, Brazil. Physica A: Statistical Mechanics and its Applications, 514:729–740, 2019. [27] Rashid Jan and Asif Jan. MSGDTM for solution of fractional order dengue disease model. International Journal of Science and Research, 6(3):1140–1144, 2017. [28] I. Ahmad, I. Ali, Rashid Jan, S. A. Idris, and M. Mousa. Solutions of a three- dimensional multi-term fractional anomalous solute transport model for contamina- tion in groundwater. PLoS ONE, 18(12):e0294348, 2023. [29] R. Zarin, A. Khan, and P. Kumar. Fractional-order dynamics of Chagas-HIV epidemic model with different fractional operators. AIMS Mathematics, 7(10):18897–18924, 2022. [30] M. C. Bahi, S. Bahramand, Rashid Jan, S. Boulaaras, H. Ahmad, and R. Guefaifia. Fractional view analysis of sexual transmitted human papilloma virus infection for public health. Scientific Reports, 14(1):3048, 2024. [31] Rashid Jan, Mona Alsulami, and N. N. A. Razak. Modeling the response of the im- mune system to HIV-tumor interaction via a fractional framework. European Journal of Pure and Applied Mathematics, 18(1):5670, 2025. [32] Rashid Jan, S. Boulaaras, A. Alharbi, and N. N. A. Razak. Nonlinear dynamics of a zoonotic disease with control interventions through fractional derivative. European Journal of Pure and Applied Mathematics, 17(4):3781–3800, 2024. [33] I. Ahmad, Rashid Jan, N. N. A. Razak, A. Khan, and T. Abdeljawad. Numerical investigation of the dynamical behavior of hepatitis B virus via Caputo-Fabrizio frac- tional derivative. European Journal of Pure and Applied Mathematics, 18(1):5509, 2025. [34] A. E. Matouk and I. Khan. Complex dynamics and control of a novel physical model using nonlocal fractional differential operator with singular kernel. Journal of Ad- vanced Research, 24:463–474, 2020. [35] A. Al-khedhairi, A. E. Matouk, and I. Khan. Chaotic dynamics and chaos control for the fractional-order geomagnetic field model. Chaos, Solitons & Fractals, 128:390– 401, 2019. [36] R. P. Agarwal, D. Baleanu, J. J. Nieto, D. F. M. Torres, and Y. Zhou. A survey on fuzzy fractional differential, and optimal control nonlocal evolution equations. Journal of Computational and Applied Mathematics, 339:3–29, 2018. [37] A. E. Matouk and B. Lahcene. Chaotic dynamics in some fractional predator-prey models via a new Caputo operator based on the generalised Gamma function. Chaos, Solitons & Fractals, 166:112946, 2023. [38] R. Almeida, D. Tavares, and D. F. M. Torres. The variable-order fractional calculus of variations. Springer, Cham, 2019. [39] L. Debnath. Recent applications of fractional calculus to science and engineering. In- ternational Journal of Mathematics and Mathematical Sciences, 2003(54):3413–3442, 2003. [40] M. Awadalla, J. Alahmadi, K. R. Cheneke, and S. Qureshi. Fractional optimal control K. Guedri et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6082 39 of 39 model and bifurcation analysis of human syncytial respiratory virus transmission dynamics. Fractal and Fractional, 8(1):44, 2024. [41] K. Diethelm. A fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlinear Dynamics, 71(4):613–619, 2013. [42] S. Rosa and D. F. M. Torres. Optimal control of a fractional order epidemic model with application to human respiratory syncytial virus infection. Chaos, Solitons & Fractals, 117:142–149, 2018. [43] A. B. Salati and M. Shamsi. Direct transcription methods based on fractional integral approximation formulas for solving nonlinear fractional optimal control problems. Communications in Nonlinear Science and Numerical Simulation, 67:334–350, 2019. [44] S. Qureshi and A. Yusuf. Fractional derivatives applied to MSEIR problems: com- parative study with real world data. The European Physical Journal Plus, 134(4):171, 2019. [45] R. Zarin. Artificial neural network-based approach for simulating influenza dynamics: a nonlinear SVEIR model with spatial diffusion. Engineering Analysis with Boundary Elements, 176:106230, 2025. [46] R. Zarin and U. W. Humphries. Analyzing spatial diffusion and vaccination strategies in malaria epidemics: a numerical approach. Modeling Earth Systems and Environ- ment, 11(3):1587–1605, 2025. [47] K. Guedri, R. Zarin, A. Khan, A. Khan, B. M. Makhdoum, and H. A. Niyazi. Model- ing hepatitis B transmission dynamics with spatial diffusion and disability potential in the chronic stage. AIMS Mathematics, 10(1):1322–1349, 2025. [48] K. Guedri, R. Zarin, A. Zeb, B. M. Makhdoum, H. A. Niyazi, and A. Khan. A numerical study of HIV/AIDS transmission dynamics and the onset of long-term disability in chronic infection. The European Physical Journal Plus, 139(12):1089, 2024. [49] Kai Diethelm. The analysis of fractional differential equations: an application- oriented exposition using operators of Caputo type. 2010. [50] S. A. Khuri. A Laplace decomposition algorithm applied to a class of nonlinear differential equations. Journal of Applied Mathematics, 1(4):141–155, 2001. [51] Roland Duduchava. The Green formula and layer potentials. Integral Equations and Operator Theory, 41(2):127–178, 2001. [52] M. M. El-Borai. Some probability density and fundamental solutions of fractional evolution equations. Chaos, Solitons & Fractals, 14(3):433–440, 2002. [53] Zhisheng Shuai and P. van den Driessche. Global stability of infectious disease models using Lyapunov functions. SIAM Journal on Applied Mathematics, 73(4):1513–1532, 2013.