Theoretical and Numerical Analysis of the SIR Model and Its Symmetric Cases with Power Caputo Fractional Derivative

Abstract

This paper introduces a novel fractional Susceptible-Infected-Recovered ($\mathbb{SIR}$) model that incorporates a power Caputo fractional derivative (PCFD) and a density-dependent recovery rate. This enhances the model’s ability to capture memory effects and represent realistic healthcare system dynamics in epidemic modeling. The model’s utility and flexibility are demonstrated through an application using parameters representative of the COVID-19 pandemic. Unlike existing fractional $\mathbb{SIR}$ models often limited in representing diverse memory effects adequately, the proposed PCFD framework encompasses and extends well-known cases, such as those using Caputo–Fabrizio and Atangana–Baleanu derivatives. We prove that our model yields bounded and positive solutions, ensuring biological plausibility. A rigorous analysis is conducted to determine the model’s local stability, including the derivation of the basic reproduction number (${\mathbf{R}}_0$) and sensitivity analysis quantifying the impact of parameters on ${\mathbf{R}}_0$. The uniqueness and existence of solutions are guaranteed via a recursive sequence approach and the Banach fixed-point theorem. Numerical simulations, facilitated by a novel numerical scheme and applied to the COVID-19 parameter set, demonstrate that varying the fractional order significantly alters predicted epidemic peak timing and severity. Comparisons across different fractional approaches highlight the crucial role of memory effects and healthcare capacity in shaping epidemic trajectories. These findings underscore the potential of the generalized PCFD approach to provide more nuanced and potentially accurate predictions for disease outbreaks like COVID-19, thereby informing more effective public health interventions.

1. Introduction

The study of infectious diseases is a critical area in epidemiology, where mathematical models serve as essential tools for understanding disease spread dynamics, predicting outbreaks, and evaluating the effectiveness of interventions. The Susceptible-Infected-Recovered ($\mathbb{SIR}$) model, a cornerstone of epidemiological modeling, classifies populations into susceptible ($\mathbb{S}$), infected ($\mathbb{I}$), and recovered ($\mathbb{R}$) compartments. Its simplicity and versatility have allowed for extensive adaptations to incorporate real-world factors. For example, Marinov [1] employed an adaptive $\mathbb{SIRV}$ model with time-dependent rates to analyze the dynamics of COVID-19, integrating data from various national contexts to capture transmission variations and to evaluate vaccination strategies. Similarly, Balderrama et al. [2] explored optimal control strategies for a $\mathbb{SIR}$ epidemic model under quarantine limitations, highlighting the trade-offs between quarantine stringency and economic or social disruption. Further illustrating the adaptability of the $\mathbb{SIR}$ model, El Khalifi [3] investigated an extended $\mathbb{SIR}$ model with gradually waning immunity, acknowledging individual heterogeneity in immune system responses and demonstrating the impact of duration of immunity on long-term disease prevalence. These studies highlight the ability of the $\mathbb{SIR}$ model to address real-world complexities, while also suggesting limitations in capturing the influence of past events on current disease dynamics.

However, traditional $\mathbb{SIR}$ models often simplify disease dynamics by assuming instantaneous interactions and neglecting historical factors, such as the influence of past infection rates on current immunity levels. To address these limitations, fractional-order derivatives, which extend the concept of differentiation to non-integer orders [4,5,6,7], have emerged as a powerful tool. In contrast to integer-order derivatives, fractional derivatives inherently incorporate memory effects and long-range interactions, leading to more accurate representations of biological processes like disease transmission and recovery. This characteristic has led to several applications in the context of $\mathbb{SIR}$ models, which offer potentially more accurate predictions of epidemic spread. For example, Alqahtani [8] analyzed a fractional-order $\mathbb{SIR}$ model that incorporates the capacity of the healthcare system, showing an improved fit of the model to the observed infection data compared to its integer-order counterpart. Kim [9] introduced a normalized time-fractional $\mathbb{SIR}$ model using a novel fractional derivative designed to improve understanding of the influence of fractional-order on epidemiological dynamics and disease prediction accuracy, specifically demonstrating a reduction in prediction error when forecasting peak infection rates. Riabi et al. [10] investigated a fractional $\mathbb{SIR}$ epidemic model with the Atangana–Baleanu–Caputo operator, utilizing the homotopy perturbation method to obtain a series solution and demonstrating an increased number of immunized individuals compared to other methods. Alazman et al. [11] introduced a diffusion component into a fractional $\mathbb{SIR}$ model and analyzed its impact using a general fractional derivative, illustrating the effects of diffusion on the model’s dynamics, revealing how diffusion can alter the spatial distribution of the infected population, a feature absent in traditional models. Beyond epidemiology, fractional calculus shows promise in areas like modeling drug delivery in pharmacokinetics, where non-local tissue interactions affect drug distribution, and in neuroscience for capturing memory effects in neuronal signaling. However, challenges remain in the widespread adoption of fractional-order models, including the computational cost of solving fractional differential equations and the difficulty in directly interpreting fractional-order parameters in terms of underlying biological mechanisms. Nevertheless, the ability of fractional calculus to capture memory and non-local effects makes it a valuable tool for enhancing the realism and predictive power of models in diverse biological systems exhibiting delayed responses or cumulative effects.

The utility of fractional calculus extends far beyond epidemiological modeling, with applications across various scientific fields. These applications leverage fractional calculus’s ability to capture complex dynamics and memory effects often observed in real-world phenomena. Examples include physics and polymer technology, where fractional calculus aids in modeling complex material behaviors [12], and electrical circuits, enabling the incorporation of fractional-order elements for enhanced circuit representation [13]. In bioengineering, fractional calculus is utilized to model biological processes [14], while in robotics, it facilitates the design of fractional-order controllers for improved performance [15]. Its utility extends to fluid mechanics, where it helps model non-Newtonian fluid behavior [16], and electrodynamics of complex media, aiding in describing materials with memory effects [17]. Control theory benefits from fractional-order controllers, offering advantages over traditional methods [18], and, as discussed previously, disease modeling leverages fractional calculus to capture memory effects and non-local interactions in epidemiological models [19]. This widespread applicability underscores the power of fractional calculus in capturing complex behaviors not adequately represented by traditional integer-order calculus.

Within epidemiological modeling, the application of fractional calculus to the $\mathbb{SIR}$ framework has been explored using various fractional derivatives, including Caputo, Caputo–Fabrizio, and Atangana–Baleanu. These derivatives incorporate memory effects and non-local interactions, leading to more realistic representations of the spread of infectious diseases [20,21,22]. The Caputo fractional derivative is well suited for systems with well-defined initial conditions, while the Caputo–Fabrizio derivative is useful for systems with less-defined initial states. The Atangana–Baleanu derivative, with its non-singular kernel, provides advantages in modeling complex dynamics with crossovers. However, many studies focus on specific fractional derivatives, potentially limiting the exploration of generalized operators that can encompass a wider range of behaviors and simulate diseases with diverse memory characteristics. The selection of an appropriate fractional derivative is a crucial consideration, as it can significantly influence the model’s properties, such as stability and the existence of solutions. While fractional $\mathbb{SIR}$ models have been applied to various diseases, a more thorough investigation of their qualitative and quantitative properties is still needed. Many previous studies have concentrated on numerical simulations, often lacking in-depth exploration of the theoretical foundations, such as the boundedness, positivity, and stability of solutions. Furthermore, rigorous comparisons between different fractional derivatives are often absent.

To address these gaps, we extend the classical $\mathbb{SIR}$ framework [23] by considering the incidence rate as ${\displaystyle \frac{2\beta \mathbb{SI}}{N}}$, which suggests a closed population with density-dependent interactions influenced by the total population size N ($N=\mathbb{S}+\mathbb{I}+\mathbb{R}$). Furthermore, we incorporate a $\delta \mathbb{R}$ term, representing the rate at which recovered individuals lose immunity and return to the susceptible compartment. Crucially, we employ a Power Caputo fractional derivative (PCFD) [24], which generalizes well-known fractional derivatives like Caputo–Fabrizio [25], Atangana–Baleanu [26], weighted Atangana–Baleanu [27], and weighted Hattaf fractional derivatives [28]. The PCFD provides a flexible and adaptable modeling framework capable of capturing diverse memory and non-local effects within disease dynamics. This work primarily focuses on the theoretical development, mathematical analysis, and numerical simulation of the PCFD SIR model to demonstrate its properties and potential. We use parameters representative of the COVID-19 pandemic to demonstrate the behavior of our novel PCFD SIR model. The construction of this paper is as follows: In Section 2, we present the construction of the proposed fractional $\mathbb{SIR}$ model. In Section 3, we recall the necessary mathematical foundations, detailing the power Caputo fractional derivative. Section 4 presents a rigorous qualitative analysis that demonstrates the boundedness and positivity of the solutions, which is critical to ensuring biological feasibility and applicability in the real world. In this section, we also investigate the stability of the disease-free equilibrium (DFE), deriving the basic reproduction number ${\mathbf{R}}_0$, a key epidemiological parameter for assessing the potential for disease spread. In addition, sensitivity analysis identifies influential parameters, revealing factors affecting disease transmission and recovery. In Section 5, we introduce a two-step Lagrange interpolation polynomial-based numerical method for approximating solutions to the fractional $\mathbb{SIR}$ model. Then, in Section 4.2, we explore symmetric model cases, including Caputo–Fabrizio, Atangana–Baleanu, and weighted Hattaf. Finally, Section 7 provides the biological interpretation of our results and conclusions.

2. Mathematical $\mathbb{SIR}$ Model

In this section, we extend the classical $\mathbb{SIR}$ model [23], given by

$\left\{\begin{array}{c}{\displaystyle \frac{d}{dz}}\mathbb{S}\left(z\right)=\Lambda -{\displaystyle \frac{\beta \mathbb{S}\mathbb{I}}{\mathbb{S}+\mathbb{I}}}-\mu \mathbb{S},\hfill \\ {\displaystyle \frac{d}{dz}}\mathbb{I}\left(z\right)={\displaystyle \frac{\beta \mathbb{S}\mathbb{I}}{\mathbb{S}+\mathbb{I}}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\gamma \mathbb{I}-\mu \mathbb{I},\hfill \\ {\displaystyle \frac{d}{dz}}\mathbb{R}\left(z\right)=\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\mu \mathbb{R},\hfill \end{array}\right.$

by considering the incidence rate ${\displaystyle \frac{2\beta \mathbb{S}\mathbb{I}}{N}}$. This suggests a closed population where interactions are density-dependent and influenced by the total population size N ($N=\mathbb{S}+\mathbb{I}+\mathbb{R}$). We assume that the total population size N is constant throughout the duration of the epidemic, which is justified by the relatively short timescale of the epidemic compared to demographic processes such as birth and death. Also, we add $\delta \mathbb{R}$ term, which represents the rate at which recovered individuals lose immunity and return to the susceptible compartment. We incorporate a power Caputo fractional-order dynamics and various important parameters related to disease transmission, recovery, and mortality. The fractional derivatives introduce memory effects and non-local interactions into the dynamics. The following model strength lies in its adaptability according to power parameters p in ${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}$, allowing for the exploration of diverse scenarios related to intervention and disease:

(1)$\left\{\begin{array}{c}{}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{S}\left(z\right)=\Lambda -{\displaystyle \frac{2\beta \mathbb{S}\mathbb{I}}{N}}-\mu \mathbb{S}+\delta \mathbb{R},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{I}\left(z\right)={\displaystyle \frac{2\beta \mathbb{S}\mathbb{I}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{R}\left(z\right)=\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R},\hfill \end{array}\right.$

with initial conditions $\mathbb{S}\left(0\right)\ge 0,\mathbb{I}\left(0\right)\ge 0$ and $\mathbb{R}\left(0\right)\ge 0.$ The definitions of parameters are presented in Table 1. The $\mathbb{SIR}$ model (1) contains three equations as follows:

Susceptible Population Dynamics ($\mathbb{S}$)

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{S}\left(z\right)=\Lambda -{\displaystyle \frac{2\beta \mathbb{S}\mathbb{I}}{N}}-\mu \mathbb{S}+\delta \mathbb{R},$

where $\Lambda$ represents the birth rate or influx of new susceptible individuals into the population. The term ${\displaystyle \frac{2\beta \mathbb{S}\mathbb{I}}{N}}$ represents the rate at which susceptible individuals become infected, $\beta$ is the transmission rate, $\mathbb{S}$ is the number of susceptible individuals, $\mathbb{I}$ is the number of infected individuals, N is the total population size ($N=\mathbb{S}+\mathbb{I}+\mathbb{R}$), assumed constant (as shown in Figure 1).

This form modifies the incidence rate ${\displaystyle \frac{\beta \mathbb{S}\mathbb{I}}{N}}$ to account for saturation effects, emphasizing that the infection rate depends on the number of individuals in $\mathbb{S}$ and $\mathbb{I}$ without relying on the recovered population $\mathbb{R}$. The factor of ‘2’ accounts for the increased contact rates within households of two individuals, where the probability of transmission is higher [29]. The term $\mu \mathbb{S}$ represents the number of susceptible individuals dying per unit time, where $\mu$ is the natural death rate. The term $\delta \mathbb{R}$ represents the rate at which recovered individuals lose immunity and return to the susceptible compartment.

Infected Population Dynamics ($\mathbb{I}$)

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{I}\left(z\right)={\displaystyle \frac{2\beta \mathbb{S}\mathbb{I}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I},$

The term ${\displaystyle \frac{2\beta \mathbb{S}\mathbb{I}}{N}}$ represents the same infection process as in the susceptible equation. The term $\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}$ represents the rate at which infected individuals recover through healthcare intervention and leave the infected compartment. Here, ${\alpha }_0$ is the baseline recovery rate attributable to healthcare intervention, ${\alpha }_1$ is the maximum recovery rate when healthcare resources are sufficient, b is a constant modulating the recovery rate based on the infected population [30]. The fraction ${\displaystyle \frac{b}{b+\mathbb{I}}}$ models the effect of healthcare resource constraints: When $\mathbb{I}$ is small, the recovery rate approaches ${\alpha }_1$, when $\mathbb{I}$ is large, the recovery rate asymptotically approaches ${\alpha }_0$. This reflects the real-world scenario in which a surge in infections can overwhelm healthcare systems, reducing the quality and availability of care for each infected individual. The term $\gamma \mathbb{I}$ represents the rate at which infected individuals die from infection. The term $\mu \mathbb{I}$ represents the number of infected individuals who die per unit time from natural causes.

Recovered Population Dynamics ($\mathbb{R}$)

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{R}\left(z\right)=\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}.$

The term $\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}$ represents the rate at which infected individuals recover through healthcare intervention and enter the recovered compartment, as described in the infected population equation. The term $\left(\mu +\delta \right)\mathbb{R}$ accounts for the removal of recovered individuals due to: Natural death ($\mu$), loss of immunity ($\delta$), causing individuals to re-enter the susceptible compartment.

3. Basic Concepts

Definition 1

([24]). For $\zeta \in \left[0,1\right),$ with $\psi ,p>0,$ and $\mathbb{X}\in H^1\left(a,b\right),$ where $H^1\left(a,b\right)$ is Sobolev space. The PCFD of order ζ, of a function $\mathbb{X}$ w.r.t the weight function $w,$ $0<w\in C^1\left(\left[a,b\right]\right)$, is defined by

(2)${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{X}\left(z\right)={\displaystyle \frac{\mathbb{PC}\left(\zeta \right)}{1-\zeta }}{\displaystyle \frac{1}{w\left(z\right)}}{\int }_a^z{}^p\mathbb{E}_{\psi ,1}\left(-{\displaystyle \frac{\zeta }{1-\zeta }}{\left(z-s\right)}^{\psi }\right){\left(w\mathbb{X}\right)}^{\prime }\left(s\right)ds,$

where ${}^p\mathbb{E}_{\psi ,1}$ is the Power Mittag-Leffler function given by

${}^p\mathbb{E}_{\psi ,1}\left(s\right)=\sum _{n=0}^{+\infty }{\displaystyle \frac{{\left(slnp\right)}^n}{\Gamma \left(kn+l\right)}},s\in \mathbb{C},and k,l,p>0,$

and $\mathbb{PC}\left(\zeta \right)$ is the normalization positive function satisfying $\mathbb{PC}\left(0\right)=\mathbb{PC}\left(1\right)=1.$

Definition 2

([24]). The Power fractional integral with order ζ, of a function $\mathbb{X},$ w.r.t the non-decreasing weight function $w,$$0<w\in C^1\left(\left[a,b\right]\right)$, is defined by

${}_{a }{}^{\mathbb{PC}}\mathbf{I}_{z,w}^{\zeta ,\psi ,p}\mathbb{X}\left(z\right)={\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}\mathbb{X}\left(z\right)+lnp{{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}}^{\mathbb{RL}}{\mathbf{I}}_{a,w}^{\psi }\mathbb{X}\left(z\right),$

where ${}^{\mathbb{RL}}\mathbf{I}_{a,w}^{\psi }\mathbb{X}\left(z\right)$ is the standard weighted R–L fractional integral of order ψ given by

${}^{\mathbb{RL}}\mathbf{I}_{a,w}^{\psi }\mathbb{X}\left(z\right)={\displaystyle \frac{1}{\Gamma \left(\psi \right)}}{\displaystyle \frac{1}{w\left(z\right)}}{\int }_a^z{\left(z-s\right)}^{\psi -1}\left(w\mathbb{X}\right)\left(s\right)ds.$

Remark 1.

The power Caputo fractional derivative, as given by Definition 1, generalizes some fractional derivatives as follows:

(1) If $p=e,w\left(z\right)=1$ and $\zeta =\psi .$ Then, the Definition 1 reduced to the following definition of ABC fractional derivative [26]

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,1}^{\zeta ,e}\mathbb{X}\left(z\right)={\displaystyle \frac{\mathbb{PC}\left(\zeta \right)}{1-\zeta }}{\int }_a^z{}^e\mathbb{E}_{\zeta ,1}\left(-{\displaystyle \frac{\zeta }{1-\zeta }}{\left(z-s\right)}^{\zeta }\right){\left(\mathbb{X}\right)}^{\prime }\left(s\right)ds.$

(2) If $p=e,w\left(z\right)=1$ and $\psi =1.$ Then, the Definition 1 reduced to the following definition of CF fractional derivative [25]

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,1}^{\zeta ,1,e}\mathbb{X}\left(z\right)={\displaystyle \frac{\mathbb{PC}\left(\zeta \right)}{1-\zeta }}{\int }_a^z{}^e\mathbb{E}_{1,1}\left(-{\displaystyle \frac{\zeta }{1-\zeta }}\left(z-s\right)\right){\left(\mathbb{X}\right)}^{\prime }\left(s\right)ds.$

(3) If $p=e$ and $\zeta =\psi .$ Then, the Definition 1 reduced to the following definition of weighted ABC fractional derivative [27]

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,e}\mathbb{X}\left(z\right)={\displaystyle \frac{\mathbb{PC}\left(\zeta \right)}{1-\zeta }}{\displaystyle \frac{1}{w\left(z\right)}}{\int }_a^z{}^p\mathbb{E}_{\zeta ,1}\left(-{\displaystyle \frac{\zeta }{1-\zeta }}{\left(z-s\right)}^{\zeta }\right){\left(w\mathbb{X}\right)}^{\prime }\left(s\right)ds.$

(4) If $p=e.$ Then, the Definition 1 reduced to the following definition of weighted generalized Hattaf fractional derivative [28]

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,e}\mathbb{X}\left(z\right)={\displaystyle \frac{\mathbb{PC}\left(\zeta \right)}{1-\zeta }}{\displaystyle \frac{1}{w\left(z\right)}}{\int }_a^z{}^p\mathbb{E}_{\psi ,1}\left(-{\displaystyle \frac{\zeta }{1-\zeta }}{\left(z-s\right)}^{\psi }\right){\left(w\mathbb{X}\right)}^{\prime }\left(s\right)ds.$

4. Behavioral Characteristics of the $\mathbb{SIR}$ Model (1)

This section undertakes a rigorous analysis of the $\mathbb{SIR}$ model defined in model (1), focusing on its fundamental properties. We will establish the boundedness and positivity of the model solutions, determine the stability of the DFE, derive the basic reproduction number ${\mathit{R}}_0$, and quantify the sensitivity of ${\mathit{R}}_0$ to variations in model parameters, thereby revealing the key drivers of disease transmission.

4.1. Analysis of Solution Boundedness

Theorem 1.

The $\mathbb{SIR}$ PCFD Model (1) yields solutions $\left(\mathbb{S},\mathbb{I},\mathbb{R}\right)$ that are guaranteed to be physically and mathematically feasible, lying within the region $\Omega ,$ where

$\Omega =\left\{\left(\mathbb{S},\mathbb{I},\mathbb{R}\right);\mathbb{S}+\mathbb{I}+\mathbb{R}\le {\displaystyle \frac{\Lambda }{\mu }}\right\}.$

Proof. 

The critical condition for the model to be biologically and mathematically feasible is that the total population size must be bounded. Thus, we have

$\begin{array}{ccc}\hfill {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}N(z& =& {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{S}\left(z\right){+}_a^{\mathbb{PC}}{\mathbf{D}}_{z,w}^{\zeta ,\psi ,p}\mathbb{I}\left(z\right){+}_a^{\mathbb{PC}}{\mathbf{D}}_{z,w}^{\zeta ,\psi ,p}\mathbb{R}\left(z\right)\hfill \\ & =& \Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}+{\displaystyle \frac{2\beta \mathbb{SI}}{N}}\hfill \\ & & -\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\hfill \\ & & +\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\hfill \\ & =& \Lambda -\mu N\left(z\right)-\gamma \mathbb{I},\hfill \end{array}$

where

$N\left(z\right)=\mathbb{S}\left(z\right)+\mathbb{I}\left(z\right)+\mathbb{R}\left(z\right).$

Clearly

$\Lambda -\mu N\left(z\right)-\gamma \mathbb{I}\le \Lambda -\mu N\left(z\right).$

This implies that

(3)${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}N\left(z\right)\le \Lambda -\mu N\left(z\right).$

Applying the Laplace transform of PCFD [24] on both sides of (3), we obtain

$\mathcal{L}\left(w{\left(z\right)}_a^{\mathbb{PC}}{\mathbf{D}}_{z,w}^{\zeta ,\psi ,p}N\left(z\right)\right)\le \mathcal{L}\left[w\left(z\right)\left(\Lambda -\mu N\left(z\right)\right)\right].$

This implies that,

$\begin{array}{ccc}\hfill N\left(z\right)& \le & {\displaystyle \frac{\Lambda }{\mu }}+{\displaystyle \frac{\mathbb{PC}\left(\zeta \right)w\left(a\right)}{\left[\mathbb{PC}\left(\zeta \right)-\left(1-\zeta \right)\mu \right]w\left(z\right)}}\left(N\left(0\right)-{\displaystyle \frac{\Lambda }{\mu }}\right){}^p\mathbb{E}_{\psi ,1}\left({\displaystyle \frac{\zeta \mu }{\left[\mathbb{PC}\left(\zeta \right)-\left(1-\zeta \right)\mu \right]}}{z}^{\zeta }\right)\hfill \\ & & -{\displaystyle \frac{\mathbb{PC}\left(\zeta \right)\mu }{\left[\mathbb{PC}\left(\zeta \right)-\left(1-\zeta \right)\mu \right]w\left(z\right)}}{}^p\mathbb{E}_{\psi ,1}\left({\displaystyle \frac{\zeta \mu }{\left[\mathbb{PC}\left(\zeta \right)-\left(1-\zeta \right)\mu \right]}}{z}^{\zeta }\right)\ast w^{\prime }\left(z\right).\hfill \end{array}$

Consequently, $N\left(z\right)$ bounded by ${\displaystyle \frac{\Lambda }{\mu }}.$ According to the fact $N\left(z\right)=\mathbb{S}\left(z\right)+\mathbb{I}\left(z\right)+\mathbb{R}\left(z\right)$, we deduce that $(\mathbb{S},\mathbb{I},\mathbb{R})$ are bounded in $\Omega$. This means ensuring the biological feasibility of the model. This boundedness result ensures that the model predicts realistic population sizes and prevents unbounded growth. □

4.2. Nonnegativity of Solutions

Theorem 2.

The $\mathbb{SIR}$ PCFD Model (1), with initial conditions ($\mathbb{S}\left(0\right)>0,\mathbb{I}\left(0\right)>0$, and $\mathbb{R}\left(0\right)>0$), guarantees non-negative solutions for all time.

Proof. 

To biological relevance of our model, all state variables must remain non-negative. We will mathematically justify this by showing that each state variable is strictly positive for all $z\in \left[a,T\right]$. This justification will be based on the analysis of the third equation in model (1), presented below:

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{R}\left(z\right)=\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}.$

Then, we have

(4)${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{R}\left(z\right)\ge -\left(\mu +\delta \right)\mathbb{R}.$

Applying Laplace transform of PCFD [24] on both sides of (4), we have

$\mathcal{L}\left[w{\left(z\right)}_a^{\mathbb{PC}}{\mathbf{D}}_{z,w}^{\zeta ,\psi ,p}\mathbb{R}\left(z\right)\right]\ge -\left(\mu +\delta \right)\mathcal{L}\left[w\left(z\right)\mathbb{R}\left(z\right)\right]\left(s\right).$

Thus, we obtain

$\mathbb{R}\left(z\right)\ge {\displaystyle \frac{\mathbb{PC}\left(\zeta \right)w\left(a\right)\mathbb{R}\left(0\right)}{\left[\mathbb{PC}\left(\zeta \right)-\left(1-\zeta \right)\left(\mu +\delta \right)\right]w\left(z\right)}}{}^p\mathbb{E}_{\psi ,1}\left({\displaystyle \frac{\zeta \left(\mu +\delta \right)}{\left[\mathbb{PC}\left(\zeta \right)-\left(1-\zeta \right)\left(\mu +\delta \right)\right]}}{z}^{\zeta }\right).$

Since $\mathbb{R}\left(0\right)>0$ and $0{\le }^p{E}_{\psi ,1}\le 1$, we determine that $\mathbb{R}\left(z\right)$ is positive for all $z\in \left[a,T\right]$. Using an analogous argument, we can show that $\mathbb{S}$ and $\mathbb{I}$ are also positive, thereby establishing the model’s biological feasibility. Consequently, the population compartments ($\mathbb{S}$, $\mathbb{I}$, $\mathbb{R}$) are guaranteed to remain non-negative, reflecting the biological reality that populations cannot have negative sizes. □

4.3. Disease-Free Equilibrium Point (DFE)

In the context of a $\mathbb{SIR}$ model, an equilibrium point is a state where the system is not changing. Mathematically, this means that the time derivatives of all the state variables are equal to zero. In this case, we are looking for values of $\mathbb{S}$, $\mathbb{I}$, and $\mathbb{R}$ where:

•. ${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{S}\left(z\right)=0$ (The rate of change of susceptible individuals is zero).

•. ${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{I}\left(z\right)=0$ (The rate of change of infected individuals is zero).

•. ${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{R}\left(z\right)=0$ (The rate of change of recovered individuals is zero).

Thus, we obtain

$\left\{\begin{array}{c}\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}=0,\hfill \\ {\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}=0,\hfill \\ \left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}=0.\hfill \end{array}\right.$

By solving the above equilibrium equations, one can easily obtain the DFE point, ${\ell }_0$, for model (1) as follows:

${\ell }_0=\left(\mathbb{S}\left(0\right),\mathbb{I}\left(0\right),\mathbb{R}\left(0\right)\right)=\left({\displaystyle \frac{\Lambda }{\mu }},0,0\right).$

4.4. Basic Reproduction Number

To derive the basic reproduction number, ${\mathit{R}}_0$, we first consider the equations governing the dynamics of the susceptible ($\mathbb{S}$) and infected ($\mathbb{I}$) compartments, which are given by:

$\begin{array}{c}{}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{S}\left(z\right)=\Lambda -{\displaystyle \frac{2\beta \mathbb{S}\mathbb{I}}{N}}-\mu \mathbb{S}+\delta \mathbb{R},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\mathbb{I}\left(z\right)={\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}.\hfill \end{array}$

The Disease-Free Equilibrium (DFE) is ${\ell }_0=(\mathbb{S}=\Lambda /\mu ,\mathbb{I}=0,\mathbb{R}=0)$. We assume a constant population size, $N=\mathbb{S}+\mathbb{I}+\mathbb{R}$, and that at equilibrium, the birth rate balances the death rate, i.e., $\Lambda \approx \mu N$. This system of equations can be expressed in a compact form as:

${}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,p}\left[\begin{array}{c}\mathbb{S}\left(z\right)\\ \mathbb{I}\left(z\right)\end{array}\right]=F\left(z\right)-V\left(z\right),$

where $F\left(z\right)$ represents the rate of new infections and $V\left(z\right)$ represents the rate of transfer out of the infected population. These are defined as:

$F\left(z\right)=\left[\begin{array}{c}-{\displaystyle \frac{2\beta \mathbb{SI}}{N}}\\ {\displaystyle \frac{2\beta \mathbb{SI}}{N}}\end{array}\right],$

and

$V\left(z\right)=\left[\begin{array}{c}\mu \mathbb{S}-\Lambda +\delta \mathbb{R}\\ \left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}+\left(\gamma +\mu \right)\mathbb{I}\end{array}\right],$

The Jacobian matrices of $F\left(z\right)$ and $V\left(z\right)$, evaluated at the DFE (${\ell }_0=(\mathbb{S}=\Lambda /\mu$, $\mathbb{I}=0$, $\mathbb{R}=0)$), denoted by $\mathcal{F}$ and $\mathcal{V}$, respectively, are:

$\mathcal{F}=\left[\begin{array}{cc}0& -{\displaystyle \frac{2\beta \Lambda }{N\mu }}\\ 0& {\displaystyle \frac{2\beta \Lambda }{N\mu }}\end{array}\right],\mathcal{V}=\left[\begin{array}{cc}\mu & 0\\ 0& {\alpha }_1+\gamma +\mu \end{array}\right].$

Using the fact that ${\mathit{R}}_0$ is the spectral radius of the next generation matrix ${\mathcal{FV}}^{-1}$, and substituting $\Lambda =\mu N$, the basic reproduction number ${\mathit{R}}_0$ for model (1) is given by

${\mathit{R}}_0={\displaystyle \frac{2\beta }{{\alpha }_1+\gamma +\mu }}.$

4.5. Stability Analysis

The stability of the Disease-Free Equilibrium (DFE) is critical in epidemiology. A locally stable DFE prevents sustained epidemics, as pathogens diminish due to low reproductive capacity. An unstable DFE risks outbreaks, even with minimal pathogen introduction. Understanding DFE stability guides public health interventions, enabling targeted control strategies to prevent disease spread.

Theorem 3.

The DFE of the $\mathbb{SIR}$ model (1) exhibits local asymptotic stability for ${\mathit{R}}_0<1$, whereas values of ${\mathit{R}}_0>1$ leads to instability in this equilibrium.

Proof. 

The model (1) is linearized at the no-disease equilibrium to examine its local stability. This procedure leads to the Jacobian matrix, $J\left({\ell }_0\right)$, which governs the dynamics of the linearized model.

$J^{\left[0\right]}=\left[\begin{array}{ccc}-\mu & -2\beta & 0\\ 0& 2\beta -\left({\alpha }_1+\gamma +\mu \right)& 0\\ 0& {\alpha }_1& -\mu \end{array}\right].$

The eigenvalues of the above matrix are ${\lambda }_1=-\mu$ (this eigenvalue has multiplicity 2) and ${\lambda }_2=2\beta -\left({\alpha }_1+\gamma +\mu \right)$. For the DFE to be locally asymptotically stable, all eigenvalues of the Jacobian matrix evaluated at the equilibrium must have strictly negative real parts. Since ${\lambda }_1=-\mu$, where $\mu$ is a positive parameter (death rate), it is always negative. Therefore, the stability is determined by the sign of ${\lambda }_2$. The DFE is stable if ${\lambda }_2<0$, that is: $2\beta -\left({\alpha }_1+\gamma +\mu \right)<0$. This can be rearranged as $2\beta <\left({\alpha }_1+\gamma +\mu \right).$ Dividing both sides by $\left({\alpha }_1+\gamma +\mu \right)$, we obtain

${\displaystyle \frac{2\beta }{{\alpha }_1+\gamma +\mu }}<1.$

This means that the DFE of the model (1) exhibits local asymptotic stability for ${\mathit{R}}_0<1$.

If ${\lambda }_2>0$, then the DFE is unstable. This corresponds to $2\beta -\left({\alpha }_1+\gamma +\mu \right)>0.$ Which, following the same steps as above leads to ${\displaystyle \frac{2\beta }{{\alpha }_1+\gamma +\mu }}>1.$ Thus, the DFE is locally asymptotically stable when ${\mathit{R}}_0<1$, and unstable when ${\mathit{R}}_0>1$. □

4.6. Sensitivity Analysis

This section is devoted to the application of sensitivity analysis of the basic reproduction number, ${\mathit{R}}_0$ with the model parameters. The derived indices elucidate the significance of individual parameters in the context of disease emergence and transmission processes. Furthermore, this sensitivity analysis serves to gauge the model’s resilience to alterations in parameter values. The following formula is used to ascertain the sensitivity indices:

${\mathbf{SEN}}_{\ell }^{{\mathit{R}}_0}={\displaystyle \frac{\ell }{{\mathit{R}}_0}}\left[{\displaystyle \frac{\partial {\mathit{R}}_0}{\partial \ell }}\right].$

Applying the above formula, we obtain the sensitivity indices of the parameters as follows:

•. ${\mathbf{SEN}}_{\beta }^{{\mathit{R}}_0}=1,$

•. ${\mathbf{SEN}}_{{\alpha }_1}^{{\mathit{R}}_0}=-0.3556,$

•. ${\mathbf{SEN}}_{\gamma }^{{\mathit{R}}_0}=-0.1531,$

•. ${\mathbf{SEN}}_{\mu }^{{\mathit{R}}_0}=-0.6436.$

The sensitivity analysis indicates that controlling the transmission rate $\beta$ is essential to mitigate disease spread, as the basic reproduction number (${\mathit{R}}_0$) exhibits the highest sensitivity to this parameter. While increasing the infected rate $(\gamma ,{\alpha }_0,{\alpha }_1,{\alpha }_2)$ also helps lower ${\mathit{R}}_0$ by increasing the recovery rate, its impact is comparatively less pronounced. The natural death rate ($\mu$) also influence ${\mathit{R}}_0$, though indirectly. Consequently, interventions that directly target transmission remain the most effective for controlling the disease within the model’s framework, followed by strategies that enhance recovery. Figure 2 presented the sensitivity of ${\mathit{R}}_0$ to each parameter in the model.

4.7. Scenario Analysis

Beyond the isolated impact of individual parameters on ${\mathit{R}}_0$, investigating the interplay between parameter pairs unlocks a more nuanced comprehension of the model’s complex dynamics. This pairwise analysis reveals how synergistic effects and countervailing influences between parameters jointly shape disease transmission, a phenomenon often obscured when considering single-parameter variations alone. The 3D contour plots shown in Figure 3 effectively visualize the responses of ${\mathit{R}}_0$ across these multifaceted parameter landscapes, thus illuminating the sensitivity of the basic reproduction number to changes in joint parameters. In addition, Figure 4 shows the contour plots of ${\mathit{R}}_0$ as combinations of pairs of two parameters, illustrating the basic reproduction number.

4.8. Lipschitz Property

By Lemma 4 in [31], we can convert the PCFD model (1) as the following equivalent integral equations:

(5)$\left[\begin{array}{c}\mathbb{S}\left(z\right)\\ \mathbb{I}\left(z\right)\\ \mathbb{R}\left(z\right)\end{array}\right]={\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}\left[\begin{array}{c}{\mathbb{S}}_0\\ {\mathbb{I}}_0\\ {\mathbb{R}}_0\end{array}\right]{+}_0^{\mathbb{PC}}{\mathbf{I}}_{z,w}^{\zeta ,\psi ,p}\left[\begin{array}{c}{\mathbb{K}}_1\left(z,\mathbb{S}\right)\\ {\mathbb{K}}_2\left(z,\mathbb{I}\right)\\ {\mathbb{K}}_3\left(z,\mathbb{R}\right)\end{array}\right],$

where

$\begin{array}{ccc}\hfill {\mathbb{K}}_1\left(z,\mathbb{S}\right)& =& \Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R},\hfill \\ \hfill {\mathbb{K}}_2\left(z,\mathbb{I}\right)& =& {\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I},\hfill \\ \hfill {\mathbb{K}}_3\left(z,\mathbb{R}\right)& =& \left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}.\hfill \end{array}$

Theorem 4.

Let $\mathbb{S},\mathbb{I},\mathbb{R},\widehat{\mathbb{S}},\widehat{\mathbb{I}},\widehat{\mathbb{R}}$ be continuous functions in $L^1\left[0,1\right]$. Define positive constants $x_1,x_2$ and $x_3$ such that

$∥\mathbb{S}∥=\underset{𝚤\in \mathcal{J}}{max}\left|\mathbb{S}(𝚤)\right|<x_1,∥\mathbb{I}∥=\underset{𝚤\in \mathcal{J}}{max}\left|\mathbb{I}(𝚤)\right|<x_2,∥\mathbb{R}∥=\underset{𝚤\in \mathcal{J}}{max}\left|\mathbb{R}(𝚤)\right|<x_3.$

Then, the following kernels

$\begin{array}{ccc}\hfill {\mathbb{K}}_1\left(z,\mathbb{S}\right)& =& \Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R},\hfill \\ \hfill {\mathbb{K}}_2\left(z,\mathbb{I}\right)& =& {\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I},\hfill \\ \hfill {\mathbb{K}}_3\left(z,\mathbb{R}\right)& =& \left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R},\hfill \end{array}$

satisfy Lipschitz conditions with Lipschitz constant $k={max}_{i=1}^3\left\{{\mathcal{L}}_{{\mathbb{K}}_i}\right\}>0$, such that

$\begin{array}{ccc}\hfill {\mathcal{L}}_{{\mathbb{K}}_1}& =& \left({\displaystyle \frac{2\beta x_2}{N}}+\mu \right),\hfill \\ \hfill {\mathcal{L}}_{{\mathbb{K}}_2}& =& \left({\displaystyle \frac{2\beta x_1}{N}}+\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]+\left(\gamma +\mu \right)\right),\hfill \\ \hfill {\mathcal{L}}_{{\mathbb{K}}_3}& =& \left(\mu +\delta \right).\hfill \end{array}$

Proof. 

For ${\mathbb{K}}_1\left(z,\mathbb{S}\right)=\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R},$ let $\mathbb{S},\widehat{\mathbb{S}}\in L^1\left[0,1\right].$ Thus,

$\begin{array}{ccc}\hfill ∥{\mathbb{K}}_1\left(z,\mathbb{S}\right)-{\mathbb{K}}_1\left(z,\widehat{\mathbb{S}}\right)∥& =& ∥\left(\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}\right)-\left(\Lambda -{\displaystyle \frac{2\beta \widehat{\mathbb{S}}\mathbb{I}}{N}}-\mu \widehat{\mathbb{S}}+\delta \mathbb{R}\right)∥\hfill \\ & \le & {\displaystyle \frac{2\beta ∥\mathbb{I}∥}{N}}∥\left(\mathbb{S}-\widehat{\mathbb{S}}\right)∥+\mu ∥\left(\mathbb{S}-\widehat{\mathbb{S}}\right)∥\hfill \\ & \le & \left({\displaystyle \frac{2\beta x_2}{N}}+\mu \right)∥{\mathbb{S}}_1-{\widehat{\mathbb{S}}}_1∥.\hfill \end{array}$

Put ${\mathcal{L}}_{{\mathbb{K}}_1}=\left({\displaystyle \frac{2\beta x_2}{N}}+\mu \right)>0.$ Thus, we get

$∥{\mathbb{K}}_1\left(z,{\mathbb{S}}_1\right)-{\mathbb{K}}_1\left(z,{\widehat{\mathbb{S}}}_1\right)∥\le {\mathcal{L}}_{{\mathbb{K}}_1}∥{\mathbb{S}}_1-{\widehat{\mathbb{S}}}_1∥.$

To further demonstrate the concept, we can obtain the following:

$∥{\mathbb{K}}_2\left(z,{\mathbb{S}}_2\right)-{\mathbb{K}}_2\left(z,{\widehat{\mathbb{S}}}_2\right)∥\le {\mathcal{L}}_{{\mathbb{K}}_2}∥{\mathbb{S}}_2-{\widehat{\mathbb{S}}}_2∥,$

and

$∥{\mathbb{K}}_3\left(z,\mathbb{R}\right)-{\mathbb{K}}_3\left(z,\widehat{\mathbb{R}}\right)∥\le {\mathcal{L}}_{{\mathbb{K}}_5}∥\mathbb{R}-\widehat{\mathbb{R}}∥.$

Let

$k={max}_{i=1}^3\left\{{\mathcal{L}}_{{\mathbb{K}}_i}\right\}>0.$

Thus, the kernels ${\mathbb{K}}_i,i=1,2,3$ are Lipschitz continuous with a Lipschitz constant $k>0.$ □

4.9. Existence of Solution via Recursive Sequences

In this section, we aim to prove the existence of a solution to the following model using a recursive sequence approach. We will use the contraction mapping theorem to show that the sequence converges to a unique solution. By (5) the solution of the model (1) is given by

$\left[\begin{array}{c}\mathbb{S}\left(z\right)\\ \mathbb{I}\left(z\right)\\ \mathbb{R}\left(z\right)\end{array}\right]={\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}\left[\begin{array}{c}{\mathbb{S}}_0\\ {\mathbb{I}}_0\\ {\mathbb{R}}_0\end{array}\right]+{\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}\left[\begin{array}{c}{\mathbb{K}}_1\left(z,\mathbb{S}\right)\\ {\mathbb{K}}_2\left(z,\mathbb{I}\right)\\ {\mathbb{K}}_3\left(z,\mathbb{R}\right)\end{array}\right]+lnp{{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}}^{\mathbb{RL}}{\mathbf{I}}_{a,w}^{\psi }\left[\begin{array}{c}{\mathbb{K}}_1\left(z,\mathbb{S}\right)\\ {\mathbb{K}}_2\left(z,\mathbb{I}\right)\\ {\mathbb{K}}_3\left(z,\mathbb{R}\right)\end{array}\right].$

Let’s represent the given system in a compact operator form. Define:

$\mathbf{X}\left(z\right)=\left[\begin{array}{c}\mathbb{S}\left(z\right)\\ \mathbb{I}\left(z\right)\\ \mathbb{R}\left(z\right)\end{array}\right],{\mathbf{X}}_0=\left[\begin{array}{c}{\mathbb{S}}_0\\ {\mathbb{I}}_0\\ {\mathbb{R}}_0\end{array}\right],$

and

$\mathbf{F}(z,\mathbf{X}\left(z\right))=\left[\begin{array}{c}{\mathbb{K}}_1\left(z,\mathbb{S}\right)\\ {\mathbb{K}}_2\left(z,\mathbb{I}\right)\\ {\mathbb{K}}_3\left(z,\mathbb{R}\right)\end{array}\right].$

Thus, the original system can be written as:

$\mathbf{X}\left(z\right)={\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}{\mathbf{X}}_0+{\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}\mathbf{F}(z,\mathbf{X}\left(z\right))+lnp{{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}}^{\mathbb{RL}}{\mathbf{I}}_{a,w}^{\psi }\mathbf{F}(z,\mathbf{X}\left(z\right)).$

Define the operator

$\mathbf{H}\left(\mathbf{X}\left(z\right)\right)={\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}{\mathbf{X}}_0+{\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}\mathbf{F}(z,\mathbf{X}\left(z\right))+lnp{{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}}^{\mathbb{RL}}{\mathbf{I}}_{a,w}^{\psi }\mathbf{F}(z,\mathbf{X}\left(z\right)).$

We define a recursive sequence of vector functions $\left\{{\mathbf{X}}_n\left(z\right)\right\},n=0,1,2,\cdots$ as follows:

${\mathbf{X}}_{n+1}\left(z\right)={\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}{\mathbf{X}}_0+{\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}{\mathbf{X}}_0+{\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}\mathbf{F}(z,{\mathbf{X}}_n\left(z\right))+lnp{{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}}^{\mathbb{RL}}{\mathbf{I}}_{a,w}^{\psi }\mathbf{F}(z,{\mathbf{X}}_n\left(z\right)).$

Theorem 5.

Assume that $w\left(z\right)\ne 0$ for all z in the considered interval $[0,T]$ and the components of $\mathbf{F}(z,\mathbf{X}(z\left)\right)$ are continuous and bounded for all $\mathbf{X}$ and z in the interval $[0,T]$. Then, the model (1) possesses a solution provided that:

$k\left[{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}+{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}{\displaystyle \frac{T^{\psi }}{\Gamma (\psi +1)}}\right]<1,$

where k is the Lipschitz constant defined in Theorem 4.

Proof. 

Let us define the operator $\mathbf{H}:\mathbf{C}\left(\left[0,T\right]\right)\to \mathbf{C}\left(\left[0,T\right]\right)$ as follows

$\mathbf{H}\left(\mathbf{X}\left(z\right)\right)={\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}{\mathbf{X}}_0+{\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}\mathbf{F}(z,\mathbf{X}\left(z\right))+lnp{{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}}^{\mathbb{RL}}{\mathbf{I}}_{a,w}^{\psi }\mathbf{F}(z,\mathbf{X}\left(z\right)).$

For all $\mathbf{X},\mathbf{Y}\in \mathbf{C}\left(\left[0,T\right]\right)$ and $z\in \left[0,T\right],$ we have

$\begin{array}{cc}\hfill ∥\mathbf{H}\left(\mathbf{X}\right)-\mathbf{H}\left(\mathbf{Y}\right)∥& \le {\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}∥\mathbf{F}(z,\mathbf{X}(z\left)\right)-\mathbf{F}(z,\mathbf{Y}(z\left)\right)∥\hfill \\ \hfill & +lnp{{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}}^{\mathbb{RL}}{\mathbf{I}}_{a,w}^{\psi }∥\mathbf{F}(z,\mathbf{X}(z\left)\right)-\mathbf{F}(z,\mathbf{Y}(z\left)\right)∥,\hfill \end{array}$

By Theorem 4, $\mathbf{F}(z,\mathbf{X}(z\left)\right)$ satisfies the Lipschitz condition $k={max}_{i=1}^3\left\{{\mathcal{L}}_{{\mathbb{K}}_i}\right\}>0$. Thus, we obtain that

$\begin{array}{cc}\hfill ∥\mathbf{H}\left(\mathbf{X}\right)-\mathbf{H}\left(\mathbf{Y}\right)∥& \le {\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}k∥\mathbf{X}-\mathbf{Y}∥\hfill \\ \hfill & +lnp{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}k{∥\mathbf{X}-\mathbf{Y}∥}^{\mathbb{RL}}{\mathbf{I}}_{0,w}^{\psi }\left(1\right)z\hfill \\ \hfill & \le k\left[{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}+{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}{\displaystyle \frac{T^{\psi }}{\Gamma (\psi +1)}}\right]∥\mathbf{X}-\mathbf{Y}∥.\hfill \end{array}$

Thus, by the Banach fixed point theorem, we have that $\mathbf{H}$ is a contraction operator. Since $\mathbf{H}$ is a contraction mapping, the sequence $\left\{{\mathbf{X}}_n\left(z\right)\right\}$ converges to a limit, which we denote by $\mathbf{X}\left(z\right)$. This means that:

$\underset{n\to \infty }{lim}{\mathbf{X}}_n\left(z\right)=\mathbf{X}\left(z\right).$

And since ${\mathbf{X}}_{n+1}\left(z\right)=\mathbf{H}\left({\mathbf{X}}_n\left(z\right)\right)$, taking limits as $n\to \infty$, and since $\mathbf{H}$ is a continuous operator. We have that:

$\underset{n\to \infty }{lim}{\mathbf{X}}_{n+1}\left(z\right)=\underset{n\to \infty }{lim}\mathbf{H}\left({\mathbf{X}}_n\left(z\right)\right)=\mathbf{H}\left(\underset{n\to \infty }{lim}{\mathbf{X}}_n\left(z\right)\right).$

and

$\mathbf{X}\left(z\right)=\mathbf{H}\left(\mathbf{X}\right(z\left)\right).$

That means the limit $\mathbf{X}\left(z\right)$ is a fixed point of the operator $\mathbf{H}$. Therefore, $\mathbf{X}\left(z\right)$ satisfies:

$\mathbf{X}\left(z\right)={\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}{\mathbf{X}}_0+{\displaystyle \frac{1-\zeta }{\mathbb{PC}\left(\zeta \right)}}\mathbf{F}(z,\mathbf{X}\left(z\right))+lnp{{\displaystyle \frac{\zeta }{\mathbb{PC}\left(\zeta \right)}}}^{\mathbb{RL}}{\mathbf{I}}_{a,w}^{\psi }\mathbf{F}(z,\mathbf{X}\left(z\right)).$

Thus, a recursive sequence of functions ${\mathbf{X}}_n\left(z\right)$ approaches the solution. This sequence converges to a unique function $\mathbf{X}\left(z\right)$, which represents the solution to the given system, according to the contraction mapping theorem. Therefore, a solution exists for the given system. □

5. Numerical Scheme with Power Caputo Fractional Derivative

We will now introduce a numerical method, based on the two-step Lagrange interpolation polynomial [32], to approximate the solution of model (1). This approach is chosen for its ability to achieve a balance between computational efficiency and accuracy in approximating solutions to systems of ordinary differential equations. The two-step nature of the method allows for the inclusion of previous solution values, improving the approximation in each iteration, while the use of a Lagrange interpolation polynomial ensures that the approximation fits the known solution points well. From (5), the solution of (1) is given by

$\begin{array}{ccc}\hfill \mathbb{S}\left(z\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z\right)}}{\mathbb{S}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left(\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|\Gamma \left(\psi \right)w\left(z\right)}}{\int }_a^z{(z-s)}^{\psi -1}w\left(s\right)\left(\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}\right)ds,\end{array}\right.\hfill \\ & & \\ \hfill \mathbb{I}\left(z\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left({\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|\Gamma \left(\psi \right)w\left(z\right)}}{\int }_a^z{(z-s)}^{\psi -1}w\left(s\right)\times \\ \left({\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)ds,\end{array}\right.\hfill \\ & & \\ \hfill \mathbb{R}\left(z\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(0\right)}{w\left(z\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}+\gamma \mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|\Gamma \left(\psi \right)w\left(z\right)}}{\int }_a^z{(z-s)}^{\psi -1}w\left(s\right)\left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right).\end{array}\right.\hfill \end{array}$

Let $z_m=a+mh$ with $m\in \mathbb{N}$ and h are the discretization step. One has

$\begin{array}{ccc}\hfill \mathbb{S}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z_m\right)}}{\mathbb{S}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left(\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\int }_a^{z_{m+1}}{(z_{m+1}-s)}^{\psi -1}w\left(s\right)\left(\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}\right)ds,\end{array}\right.\hfill \\ & & \\ \hfill \mathbb{I}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z_m\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left({\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\int }_a^{z_{m+1}}{(z_{m+1}-s)}^{\psi -1}w\left(s\right)\times \\ \left({\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)ds,\end{array}\right.\hfill \\ & & \\ \hfill \mathbb{R}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(0\right)}{w\left(z_m\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\int }_a^{z_{m+1}}{(z_{m+1}-s)}^{\psi -1}w\left(s\right)\left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right)ds,\end{array}\right.\hfill \end{array}$

which yields

$\begin{array}{ccc}\hfill \mathbb{S}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z_m\right)}}{\mathbb{S}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left(\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\sum }_{l=0}^m{\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}w\left(s\right)\times \\ \left(\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}\right)ds,\end{array}\right.\hfill \\ & & \\ \hfill \mathbb{I}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z_m\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left({\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\sum }_{l=0}^m{\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}w\left(s\right)\times \\ \left({\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)ds,\end{array}\right.\hfill \\ & & \\ \hfill \mathbb{R}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(0\right)}{w\left(z_m\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}+\gamma \mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\sum }_{l=0}^m{\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}w\left(s\right)\times \\ \left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}+\gamma \mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right)ds.\end{array}\right.\hfill \end{array}$

By Lagrange interpolation polynomial through the points $\left(z_{l-1},\mathbb{S}\left(z_{l-1}\right),\mathbb{I}\left(z_{l-1}\right),\mathbb{R}\left(z_{l-1}\right)\right)$ and $\left(z_l,\mathbb{S}\left(z_l\right),\mathbb{I}\left(z_l\right),\mathbb{R}\left(z_l\right)\right),l=1,2,3,\cdots ,m$ and $h=z_{l-1}-z_l,$ we obtain

(6)$\begin{array}{ccc}\hfill \mathbb{S}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z_m\right)}}{\mathbb{S}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left(\Lambda -{\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\mu \mathbb{S}+\delta \mathbb{R}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\sum }_{l=1}^m\left[{\displaystyle \frac{w\left(l-1\right)\left(\Lambda -\frac{2\beta \mathbb{SI}}{N}-\mu \mathbb{S}+\delta \mathbb{R}\right)\left(l-1\right)}{h}}\times \right.\\ {\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}(z_l-s)ds+{\displaystyle \frac{w\left(l\right)\left(\Lambda -\frac{2\beta \mathbb{SI}}{N}-\mu \mathbb{S}+\delta \mathbb{R}\right)\left(l\right)}{h}}\times \\ \left.{\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}(s-z_{l-1})ds\right],\end{array}\right.\hfill \end{array}$

(7)$\begin{array}{ccc}\hfill \mathbb{I}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z_m\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left({\displaystyle \frac{2\beta \mathbb{SI}}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\sum }_{l=1}^m\left[{\displaystyle \frac{w\left(l-1\right)\left(\frac{2\beta \mathbb{SI}}{N}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)\left(l-1\right)}{h}}\times \right.\\ {\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}(z_l-s)ds+{\displaystyle \frac{w\left(l\right)\left(\frac{2\beta \mathbb{SI}}{N}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I}\right)\left(l\right)}{h}}\times \\ \left.{\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}(s-z_{l-1})ds\right]\end{array}\right.\hfill \end{array}$

(8)$\begin{array}{ccc}\hfill \mathbb{R}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(0\right)}{w\left(z_m\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}\left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right)\left(m\right)\\ +{\displaystyle \frac{|ln(p\left)\right|\zeta }{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma \left(\psi \right)w\left(z_m\right)}}{\sum }_{l=1}^m\left[{\displaystyle \frac{w\left(l-1\right)\left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right)\left(l-1\right)}{h}}\times \right.\\ {\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}(z_l-s)ds+{\displaystyle \frac{w\left(l\right)\left(\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}\right)\left(l\right)}{h}}\times \\ \left.{\int }_{z_l}^{z_{l+1}}{(z_{l+1}-s)}^{\psi -1}(s-z_{l-1})ds\right]\end{array}\right.\hfill \end{array}$

Furthermore, we have

(9)${\int }_{z_l}^{z_{l+1}}{\left(z_{m+1}-s\right)}^{\psi -1}\left(z_l-s\right)ds={\displaystyle \frac{h^{\psi +1}}{\psi \left(\psi +1\right)}}\left[{\left(m-l\right)}^{\psi }\left(m-l+1+\psi \right)-{\left(m-l+1\right)}^{\psi +1}\right],$

and

(10)${\int }_{z_l}^{z_{l+1}}{\left(z_{m+1}-s\right)}^{\psi -1}\left(s-z_{l-1}\right)ds={\displaystyle \frac{h^{\psi +1}}{\psi \left(\psi +1\right)}}\left[\begin{array}{c}{\left(m-l+1\right)}^{\psi }\left(m-l+2+\psi \right)\\ -{\left(m-l\right)}^{\psi }\left(m-l+2+2\psi \right)\end{array}\right].$

Thus, by (9) and (10), the Equations (6)–(8) becomes as follows

(11)$\begin{array}{ccc}\hfill \mathbb{S}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z_m\right)}}{\mathbb{S}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}{\wp }_1\left(z_m,\mathbb{S}\left(z_m\right)\right)\left(m\right)\\ +{\displaystyle \frac{|ln\left(p\right)|h^{\psi }}{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma (\psi +2)w\left(z_m\right)}}{\sum }_{l=1}^m\left[w\left(l-1\right){\wp }_1\left(z_{l-1},\mathbb{S}\left(z_{l-1}\right)\right){\mathcal{A}}_{m,l}^{\psi }\right.\\ \left.+w\left(l\right){\wp }_1\left(z_l,\mathbb{S}\left(z_l\right)\right)\left(l\right){\mathcal{B}}_{m,l}^{\psi }\right],\end{array}\right.\hfill \end{array}$

(12)$\begin{array}{ccc}\hfill \mathbb{I}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(a\right)}{w\left(z_m\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}{\wp }_2\left(z_m,\mathbb{I}\left(z_m\right)\right)\\ +{\displaystyle \frac{|ln\left(p\right)|h^{\psi }}{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma (\psi +2)w\left(z_m\right)}}{\sum }_{l=1}^m\left[w\left(l-1\right){\wp }_2\left(z_{l-1},\mathbb{I}\left(z_{l-1}\right)\right){\mathcal{A}}_{m,l}^{\psi }\right.\\ \left.+w\left(l\right){\wp }_2\left(z_l,\mathbb{I}\left(z_l\right)\right)\left(l\right){\mathcal{B}}_{m,l}^{\psi }\right]\end{array}\right.\hfill \end{array}$

(13)$\begin{array}{ccc}\hfill \mathbb{R}\left(z_{m+1}\right)& =& \left\{\begin{array}{c}{\displaystyle \frac{w\left(0\right)}{w\left(z_m\right)}}{\mathbb{I}}_0+{\displaystyle \frac{1-\zeta }{\left|\mathbb{PC}\right(\zeta \left)\right|}}{\wp }_3\left(z_m,\mathbb{R}\left(z_m\right)\right)\\ +{\displaystyle \frac{|ln\left(p\right)|h^{\psi }}{\left|\mathbb{PC}\left(\zeta \right)\right|\Gamma (\psi +2)w\left(z_m\right)}}{\sum }_{l=1}^m\left[w\left(l-1\right){\wp }_3\left(z_{l-1},\mathbb{R}\left(z_{l-1}\right)\right){\mathcal{A}}_{m,l}^{\psi }\right.\\ \left.+w\left(l\right){\wp }_3\left(z_l,\mathbb{R}\left(z_l\right)\right)\left(l\right){\mathcal{B}}_{m,l}^{\psi }\right]\end{array}\right.\hfill \end{array}$

where

$\begin{array}{ccc}\hfill {\wp }_1\left(z,\mathbb{S}\left(z\right)\right)& =& \Lambda -{\displaystyle \frac{2\beta \mathbb{S}\left(z\right)\mathbb{I}\left(z\right)}{N}}-\mu \mathbb{S}\left(z\right)+\delta \mathbb{R}\left(z\right),\hfill \\ \hfill {\wp }_2\left(z,\mathbb{I}\left(z\right)\right)& =& {\displaystyle \frac{2\beta \mathbb{S}\left(z\right)\mathbb{I}\left(z\right)}{N}}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}\left(z\right)-\left(\gamma +\mu \right)\mathbb{I}\left(z\right),\hfill \\ \hfill {\wp }_3\left(z,\mathbb{R}\left(z\right)\right)& =& \left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]\mathbb{I}\left(z\right)-\left(\mu +\delta \right)\mathbb{R}\left(z\right),\hfill \\ \hfill {\mathcal{A}}_{m,l}^{\psi }& =& {\left(m-l\right)}^{\psi }\left(m-l+1+\psi \right)-{\left(m-l+1\right)}^{\psi +1},\hfill \\ \hfill {\mathcal{B}}_{m,l}^{\psi }& =& {\left(m-l+1\right)}^{\psi }\left(m-l+2+\psi \right)-{\left(m-l\right)}^{\psi }\left(m-l+2+2\psi \right).\hfill \end{array}$

6. SIR Model on COVID-19

A key strength of this model lies in its enhanced capabilities to simulate a range of real-world infectious disease scenarios. Specifically, the model incorporates the $\delta \mathbb{R}$ term (where $\delta$ represents the rate of immunity loss), thereby enabling the capture of diseases where protection following infection is not lifelong, such as influenza. In such cases, the model can effectively investigate the initial propagation of novel strains within a susceptible population, and providing insights into the effectiveness of early intervention strategies. Furthermore, the PCFD employs a flexible framework for capturing diverse memory and non-local effects within disease dynamics according to its power parameter p, and generalizes well-known fractional derivatives. Moreover, a density-dependent recovery rate, represented mathematically by $\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right){\displaystyle \frac{b}{b+\mathbb{I}}}\right]$, accounts for the impact of healthcare resource limitations, a feature particularly relevant for simulating outbreaks where access to medical care significantly influences outcomes. In this section, we illustrate the application and behavior of the SIR model using parameters representative of the COVID-19 pandemic (Table 2) to provide a concrete real-world example and motivate the use of this advanced fractional framework.

Here, we consider $z\in \left[0,1000\right]$, and the values of parameters as in Table 2 with initial conditions $\left({\mathbb{S}}_0,{\mathbb{I}}_0,{\mathbb{R}}_0\right)=\left(90,40,30\right).$ The complete code of simulations is provided in a GitHub repository via the link: https://github.com/Almalahi/COMPLETE-CODE-SIR-MODEL (accessed on 29 March 2025).

By these values, with PCFD model (1), we have Figure 5, Figure 6 and Figure 7 present a graphical depiction of the $\mathbb{S}$, $\mathbb{I}$ and $\mathbb{R}$ with different fractional order of the PCFD model (1).

These visualizations offer a direct view into the dynamic interplay of the three epidemiological classes in different cases, illustrating their temporal evolution.

7. Symmetric Cases of Model (1)

The fractional derivative employed within model (1) offers a high degree of generalization, encompassing several symmetric cases contingent upon the specific choices of its parameters $\zeta ,\psi$, the fractional derivative’s power, p, and the weighting function, $w\left(z\right)$. In the ensuing subsections, we will explore and analyze simulations of these distinct symmetric scenarios, for comparison, to highlight the flexibility of the PCFD and illustrate the versatility and richness of the fractional model using the COVID-19 representative parameter set.

$\mathbb{SIR}$ COVID-19 Model with Caputo–Fabrizio Fractional Approach

If $w\left(z\right)=1,p=e,\psi =1.$ Then, the model (1) reduce to the Caputo–Fabrizio fractional COVID-19 model given by

(14)$\left\{\begin{array}{c}{}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,1}^{\zeta ,1,e}\mathbb{S}\left(z\right)=\Lambda -\frac{2\beta \mathbb{SI}}{N}-\mu \mathbb{S}+\delta \mathbb{R},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,1}^{\zeta ,1,e}\mathbb{I}\left(z\right)=\frac{2\beta \mathbb{SI}}{N}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,1}^{\zeta ,1,e}\mathbb{R}\left(z\right)=\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}.\hfill \end{array}\right.$

With the same parameter values in Table 2, the graphs of approximate solutions in case of Caputo–Fabrizio model (14) are given as follows:

Figure 8 and Figure 9 provide a detailed graphical representation of the classes $\mathbb{S}$, $\mathbb{I}$, and $\mathbb{R}$, with $w\left(z\right)=1,p=e,\psi =1.$ in different fractional order of the Caputo–Fabrizio fractional model (14).

$\mathbb{SIR}$ COVID-19 Model with Atangana–Baleanu Fractional Approach

If $w\left(z\right)=1,p=e,\zeta =\psi .$ Then, the model (1) reduce to the Atangana–Baleanu fractional COVID-19 model given by

(15)$\left\{\begin{array}{c}{}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,1}^{\zeta ,\zeta ,e}\mathbb{S}\left(z\right)=\Lambda -\frac{2\beta \mathbb{SI}}{N}-\mu \mathbb{S}+\delta \mathbb{R},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,1}^{\zeta ,\zeta ,e}\mathbb{I}\left(z\right)=\frac{2\beta \mathbb{SI}}{N}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,1}^{\zeta ,\zeta ,e}\mathbb{R}\left(z\right)=\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}.\hfill \end{array}\right.$

The graphs of approximate solutions of Atangana–Baleanu fractional model (15) are given as follows:

•. Figure 10 provide a detailed graphical depiction of $\mathbb{S}$, $\mathbb{I}$, and $\mathbb{R}$ populations as simulated by the Atangana–Baleanu fractional model (15) with $w\left(z\right)=1,p=e,\zeta =\psi$ and $\zeta =0.3,0.35,0.4,0.45$.

$\mathbb{SIR}$ COVID-19 Model with Weighted Atangana–Baleanu Fractional Approach

If $p=e,\zeta =\psi .$ Then, the model (1) reduce to the weighted Atangana–Baleanu fractional model given by

(16)$\left\{\begin{array}{c}{}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\zeta ,e}\mathbb{S}\left(z\right)=\Lambda -\frac{2\beta \mathbb{SI}}{N}-\mu \mathbb{S}+\delta \mathbb{R},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\zeta ,e}\mathbb{I}\left(z\right)=\frac{2\beta \mathbb{SI}}{N}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\zeta ,e}\mathbb{R}\left(z\right)=\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}.\hfill \end{array}\right.$

The graphs of approximate solutions of weighted Atangana–Baleanu fractional model (16) are given as follows:

•. Figure 11 and Figure 12 provide a detailed graphical depiction of $\mathbb{S}$, $\mathbb{I}$, and $\mathbb{R}$ populations as simulated by the weighted Atangana–Baleanu fractional model (16) with $w\left(z\right)=z+1$, $p=e,\zeta =\psi$ with different fractional order.

$\mathbb{SIR}$ COVID-19 Model with Weighted Generalized Hattaf Fractional Approach

If $p=e.$ Then, the model (1) reduce to the weighted generalized Hattaf fractional model given by

(17)$\left\{\begin{array}{c}{}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,e}\mathbb{S}\left(z\right)=\Lambda -\frac{2\beta \mathbb{SI}}{N}-\mu \mathbb{S}+\delta \mathbb{R},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,e}\mathbb{I}\left(z\right)=\frac{2\beta \mathbb{SI}}{N}-\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\gamma +\mu \right)\mathbb{I},\hfill \\ {}_{a }{}^{\mathbb{PC}}\mathbf{D}_{z,w}^{\zeta ,\psi ,e}\mathbb{R}\left(z\right)=\left[{\alpha }_0+\left({\alpha }_1-{\alpha }_0\right)\frac{b}{b+\mathbb{I}}\right]\mathbb{I}-\left(\mu +\delta \right)\mathbb{R}.\hfill \end{array}\right.$

The graphs of approximate solutions of weighted generalized Hattaf fractional model (17) are given as follows:

Figure 13 and Figure 14 provide a detailed graphical depiction of $\mathbb{S}$, $\mathbb{I}$, and $\mathbb{R}$ classes as simulated by the weighted generalized Hattaf fractional model (17) with $w\left(z\right)=z+1,p=e$ with different fractional order.

Table 3, Table 4, Table 5, Table 6 and Table 7 address the comparison of fractional models and the standard integer-order SIR model against the “Actual COVID-19 Trend” characteristics.

8. Discussion and Biological Interpretation

This study introduced a novel SIR model incorporating a generalized PCFD and applied it using parameters representative of the COVID-19 pandemic (Table 2) to provide a concrete real-world example and illustrate the potential of this advanced fractional framework. The simulation results, presented in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 and summarized comparatively in Table 3, Table 4, Table 5, Table 6 and Table 7, offer significant insights into how fractional calculus, particularly the flexible PCFD approach, can capture diverse epidemic dynamics. This section discusses the biological interpretation of these findings, evaluating the performance of different fractional derivatives against the standard integer-order model and benchmark characteristics derived from actual COVID-19 trends, thereby addressing the need to demonstrate the model’s relevance and potential advantages through illustrative simulations.

Our analysis reveals that the choice of fractional derivative and its associated parameters ($\zeta ,p,\psi ,w\left(z\right)$) profoundly influences the predicted epidemic trajectory, even when using the same underlying parameter set (Table 2). This is clearly demonstrated in Table 4, Table 5, Table 6 and Table 7, where different symmetric cases of the PCFD yield markedly different alignments with the benchmark COVID-19 trend. This highlights the importance of selecting an appropriate modeling framework and tuning its parameters carefully for specific applications. The fractional order, $\zeta$, is particularly influential, primarily modulating the “memory” embedded in the system—how strongly past events influence present dynamics. The ability to adjust this memory effect via $\zeta$ is key to potentially achieving improved alignment with real-world data. For instance, compared to the baseline integer-order model which showed only moderate alignment (Table 4, Table 5, Table 6 and Table 7), specific fractional models like Caputo–Fabrizio ($\zeta =0.85$, Table 4) and Weighted Hattaf ($\zeta =0.95$, Table 7) demonstrated high alignment, successfully capturing the peak timing, magnitude, and duration characteristics of the benchmark trend much more closely. This improved fit, as summarized in the overall comparison (Table 3), suggests that the memory effects implicit in these specific fractional orders better represent the underlying dynamics of the illustrative COVID-19 scenario than the standard derivative. Conversely, other fractional derivatives like Atangana–Baleanu (Table 5) and Weighted Atangana–Baleanu (Table 6) showed low alignment for the tested parameters, emphasizing that simply using any fractional derivative does not guarantee superiority.

The PCFD model, by its generalized nature encompassing these various forms, allows for tuning these elements, offering enhanced flexibility to potentially match specific disease characteristics more accurately than restrictive models. We now examine the behavior of each population compartment, interpreting the simulation results (Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14) in light of the comparative evaluation Tables:

•. Susceptible Population ($\mathbb{S}$): As expected, $\mathbb{S}$ initially declines in all simulations. However, the rate of decline and subsequent recovery or stabilization varies significantly, impacting the overall epidemic duration and alignment score. Models achieving high alignment (Table 4 and Table 7) exhibit $\mathbb{S}$ dynamics consistent with the benchmark epidemic duration (130–140 days), showing significant depletion by the peak infection time (e.g., S ≈ 220 for CF, S ≈ 200 for WGH at peak) and partial recovery towards the end (S ≈ 200 for CF, S ≈ 180 for WGH). In contrast, models with lower alignment, such as Atangana–Baleanu (Table 5), show dynamics (e.g., S ≈ 150 at end) reflecting the much longer predicted epidemic duration (∼250 days). The diversity in S-dynamics across Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 illustrates the PCFD framework’s capacity to represent varied scenarios, including those with potentially faster (e.g., Figure 11, WAB) or slower (e.g., Figure 10, AB) susceptibility changes compared to the benchmark.

•. Infected Population ($\mathbb{I}$): The dynamics of the $\mathbb{I}$ compartment are central to the comparative evaluation. The benchmark trend showed a peak around 50–55 days with a relative magnitude of ∼39–42 individuals. The integer-order model predicted an earlier (35 days) and higher (48 individuals) peak (Table 4, Table 5, Table 6 and Table 7). Significantly, the Caputo–Fabrizio ($\zeta =0.85$, Table 4) and Weighted Hattaf ($\zeta =0.95$, Table 7) models closely matched the benchmark peak time (∼40/∼50 days) and magnitude (∼39/∼40 individuals). This successful replication highlights the potential of these fractional approaches (summarized in Table 3). In contrast, the Atangana–Baleanu model (Table 5) predicted an early peak (∼35 days), and the Weighted Atangana–Baleanu model (Table 6) predicted either a very early peak (∼22 days for $\zeta =0.80$) or a lower peak magnitude (∼30 for $\zeta =0.45$), both failing to align well with the benchmark $\mathbb{I}$ curve characteristics. The PCFD’s ability to generalize allows it, in principle, to capture dynamics ranging from the well-aligned cases (like WGH) to the less aligned ones (like AB), depending on the chosen parameters ($p,\psi ,w\left(z\right),\zeta$). The modulation of peak characteristics via the fractional definition, combined with the density-dependent recovery term [30], is crucial for realistic simulation.

•. Recovered Population ($\mathbb{R}$): The accumulation of the $\mathbb{R}$ population reflects the progression towards the end of the epidemic wave. In models with high alignment (Table 4 and Table 7), the $\mathbb{R}$ curve rises steadily and approaches its plateau around the benchmark duration of 130–140 days (reaching R ≈ 160 for CF, R ≈ 180 for WGH). This contrasts sharply with models showing poor duration alignment, like Atangana–Baleanu (Table 5), where the $\mathbb{R}$ population continues to rise significantly beyond 140 days, reaching R ≈ 210 only around 250 days. The diverse shapes of the R curves in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 again showcase the flexibility conferred by the fractional derivative choice, influencing factors like apparent recovery speed and the final proportion recovered within a given timeframe, relevant to understanding immunity accumulation [33].

In conclusion, this section explicitly addressed the need for demonstrating the real-world relevance and motivation of the proposed PCFD SIR model through comparative evaluation (Table 3). By applying the model using COVID-19 representative parameters (Table 2) and comparing the outcomes against a benchmark trend (Table 4, Table 5, Table 6 and Table 7), we have shown that specific fractional derivatives generalized by the PCFD (namely Caputo–Fabrizio and Weighted Hattaf under the tested conditions) can offer superior alignment compared to the standard integer-order model. The primary motivation for using the generalized PCFD framework lies in its inherent flexibility to capture a wider spectrum of dynamics—particularly varying memory effects influencing transmission, peak characteristics, and recovery patterns—than is possible with standard integer-order models or single fixed fractional derivatives. The results clearly indicate that the choice of fractional derivative significantly impacts predicted epidemic dynamics, and careful selection or fitting is crucial. The comparative tables strongly suggest that the PCFD approach offers a valuable and adaptable tool for exploring.

9. Conclusions

This study introduced and analyzed a novel fractional Susceptible-Infected-Recovered (SIR) model incorporating PCFD and a density-dependent recovery rate reflecting healthcare capacity constraints. We proved solutions’ boundedness and positivity, analysed the stability of the disease-free equilibrium, derived an explicit formula for the basic reproduction number ($R_0$), and conducted a sensitivity analysis. The analysis confirms the biological plausibility of the model and reveals the dominant influence of the transmission rate ($\beta$) on $R_0$. Numerical simulations vividly demonstrate the significant impact of the fractional order ($\zeta$) on crucial epidemic characteristics, such as peak timing and severity. This highlights a core strength and challenge of fractional modeling: the choice of derivative and its associated parameters—the fractional order $\zeta$, the power parameters p and $\psi$, and the weighting function $w\left(z\right)$—collectively determine the type and strength of memory embedded within the model. These choices profoundly influence predicted epidemic dynamics in simulation, modulating how past events shape current infection rates, recovery processes, and mortality, thereby substantially altering projections of epidemic spread, peak characteristics, and overall duration. For instance, lower fractional orders generally emphasize longer-term historical dynamics, while higher orders prioritize more recent events. The PCFD framework’s generality, encompassing specific derivatives like Caputo–Fabrizio, Atangana–Baleanu, and generalized Hattaf (including weighted variants) as special cases, offers significant flexibility. However, this underscores the critical importance of selecting or fitting these fractional parameters appropriately for specific disease contexts, as different choices lead to distinct predictions in model outputs. Furthermore, the model’s inclusion of detailed recovery pathways (both dependent on and independent of healthcare intervention) and an infection-induced death rate enhances its realism in representing diverse disease outcomes and the impact of healthcare systems. The comparison Table 3, Table 4, Table 5, Table 6 and Table 7 demonstrates that models incorporating fractional derivatives—particularly the Caputo–Fabrizio (CF) and weighted generalized Hattaf (WGH) cases—yield predictions that are more consistent with observed data in terms of peak infection timing and total case count. The advantage of fractional-order derivatives over classical models lies in their inherent ability to capture memory and hereditary properties of the infection dynamics. This allows the model to account for the influence of historical infection rates on current behavior—something integer-order models fundamentally lack. As seen in our simulations and comparative analysis, fractional models adjust more effectively to real-world outbreak patterns, thereby offering superior descriptive and predictive power. One of the unique strengths of our approach is the use of the Power Caputo Fractional Derivative (PCFD), which serves as a unifying operator encompassing various well-known fractional derivatives as special cases. This flexibility not only provides a broader mathematical foundation but also allows the model to be calibrated based on specific memory kernels suited to different types of epidemics. Such generality enhances the model’s adaptability across a spectrum of diseases with varying temporal characteristics.. Future work will focus on extending this model to incorporate spatial dynamics and age-structured populations, as well as calibrating and validating the model against specific real-world epidemiological time-series data, aiming to further enhance its utility for detailed epidemic forecasting and control.

Author Contributions

Conceptualization, M.S.A. and M.A.; methodology, M.A.; software, M.A.; formal analysis, M.S.; investigation, A.S.A. and F.A.A.; writing—original draft preparation, M.A.; writing—review and editing, M.S., K.A. and B.Y.; visualization, F.A.A.; supervision, M.A. and K.A.; project administration, K.A.; funding acquisition, M.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

All data used in this work are contained within the article.

Acknowledgments

The Researchers would like to thank the Deanship of Graduate Studies and Scientific Research at Qassim University for financial support (QU-APC-2025).

Conflicts of Interest

There are no conflicts of interest declared by the authors.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Schematic diagram of the modified $\mathbb{SIR}$ model.

Figure 2 Sensitivity indices of $\mathbb{SIR}$ model parameters.

Figure 3 Contour 3D plots for illuminating the basic reproduction number ${\mathit{R}}_0$ sensitivity to joint parameter changes.

Figure 4 Contour 2D plots for illuminating the basic reproduction number ${\mathit{R}}_0$ sensitivity to joint parameter changes.

Figure 5 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes with $p=10$, $\psi =2$ and $\zeta =0.3,0.35,0.4,0.45$ of the power-law Caputo fractional model.

Figure 6 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes with $p=100$, $\psi =2.5$ and $\zeta =0.8,0.85,0.9,0.95$ of the power-law Caputo fractional model.

Figure 7 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes with $p=10$, $\psi =2.5$ and $\zeta =0.8,0.85,0.9,0.95$ of the power-law Caputo fractional model.

Figure 8 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes as simulated by the Caputo–Fabrizio fractional model (14) with $w\left(z\right)=1,p=e,\psi =1.$ and $\zeta =0.8,0.85,0.9,0.95$.

Figure 9 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes as simulated by the Caputo–Fabrizio fractional model (14) with $w\left(z\right)=1,p=e,\psi =1.$ and $\zeta =0.3,0.35,0.4,0.45$.

Figure 10 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes as simulated by the Atangana–Baleanu fractional model (15) with $w\left(z\right)=1,p=e,\zeta =\psi$ and $\zeta =0.3,0.35,0.4,0.45$.

Figure 11 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes as simulated by the weighted Atangana–Baleanu fractional model (16) with $w\left(z\right)=z+1,p=e,\zeta =\psi$ and $\zeta =0.8,0.85,0.9,0.95$.

Figure 12 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$, and Recovered $\mathbb{R}$ classes as simulated by the weighted Atangana–Baleanu fractional model (16) with $p=e,\zeta =\psi$.

Figure 13 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes as simulated by the weighted generalized Hattaf fractional model (17) with $w\left(z\right)=z+1,p=e$ and and $\zeta =0.8,0.85,0.9,0.95$.

Figure 14 Graphical depiction of the Susceptible $\mathbb{S}$, Infected $\mathbb{I}$ and Recovered $\mathbb{R}$ classes as simulated by the weighted generalized Hattaf fractional model (17) with $w\left(z\right)=z+1,p=e$ and and $\zeta =0.3,0.35,0.4,0.45$.

Table 1

Description of Model Parameters.

Parameter	Description	Units
$\Lambda$	Recruitment rate	Individual/Time
$\beta$	Transmission rate	(Individual · Time${)}^{-1}$
$\mu$	Natural death rate	${\mathrm{Time}}^{-1}$
$\delta$	Rate of loss of immunity	${\mathrm{Time}}^{-1}$
${\alpha }_0$	Baseline recovery rate attributable to healthcare intervention	${\mathrm{Time}}^{-1}$
${\alpha }_1$	Maximum recovery rate when healthcare resources are sufficient	${\mathrm{Time}}^{-1}$
b	Density dependence influence parameter	Individual
$\gamma$	Infection-induced death rate	${\mathrm{Time}}^{-1}$
N	Total population size	Individual

Table 2

Values of Model Parameters.

Parameter	Value	Units	Ref.
$\Lambda$	1.75	Individual/Time	[29]
$\beta$	0.01	(Individual · Time${)}^{-1}$	[33]
$\mu$	0.005	${\mathrm{Time}}^{-1}$	[29]
$\delta$	0.04	${\mathrm{Time}}^{-1}$	[34]
${\alpha }_0$	0.2	${\mathrm{Time}}^{-1}$	[29]
${\alpha }_1$	0.21	${\mathrm{Time}}^{-1}$	[23]
b	0.3	Individual	[29]
$\gamma$	0.2	${\mathrm{Time}}^{-1}$	[35]

Table 3

Comparative Evaluation Between Classical and power Fractional SIR Models vs. Actual COVID-19 Data (Duration for $\zeta =0.85$ adjusted based on visual inspection of Figure 6 and Figure 7).

Model Type	Peak Inf. Time (days)	Max Infected Individuals	Epidemic Duration (days)	Alignment with Real Data
Integer-order SIR [23,29]	35	48	∼90	Moderate
Fractional SIR ($\zeta =0.85$)	50	40	∼150	High
Actual COVID-19 Trend [36]	50–55	∼39–42	∼130–140	—

Table 4

Caputo–Fabrizio SIR Model ($w\left(z\right)=1,p=e,\psi =1.$) vs. Actual COVID-19 Data.

Model Type	Peak Inf. Time (days, approx.)	Max Infected (Imax, approx.)	Epidemic Duration (days, approx.)	Alignment with Real Data
Integer-order SIR [23,29]	35	48	∼90	Moderate
Fractional SIR (Caputo–Fabrizio, $\zeta =0.85$) Figure 8 and Figure 9	∼40	∼39	∼130	High
Actual COVID-19 Trend [36]	50–55	∼39–42	∼130–140	—

Table 5

Atangana–Baleanu SIR Model ($w\left(z\right)=1,p=e,\zeta =\psi .$) vs. Actual COVID-19 Data.

Model Type	Peak Inf. Time (days, approx.)	Max Infected (Imax, approx.)	Epidemic Duration (days, approx.)	Alignment with Real Data
Integer-order SIR [23,29]	35	48	∼90	Moderate
Fractional SIR (Atangana–Baleanu, $\zeta =0.45$) Figure 10	∼35	∼38	∼250	Low to Moderate
Actual COVID-19 Trend [36]	50–55	∼39–42	∼130–140	—

Table 6

Weighted Atangana–Baleanu SIR Model ($w\left(z\right)=z+1,p=e,\zeta =\psi$) vs. Actual COVID-19 Data.

Model Type	Peak Inf. Time (days, approx.)	Max Infected (Imax, approx.)	Epidemic Duration (days, approx.)	Alignment with Real Data
Integer-order SIR [23,29]	35	48	∼90	Moderate
Fractional SIR (Weighted AB, $\zeta =0.80$) Figure 11 and Figure 12	∼22	∼40	∼110	Low
Fractional SIR (Weighted AB, $\zeta =0.45$) Figure 11 and Figure 12	∼45	∼30	∼250	Low
Actual COVID-19 Trend [36]	50–55	∼39–42	∼130–140	—

Table 7

Weighted Generalized Hattaf SIR Model ($w\left(z\right)=z+1,p=e$) vs. Actual COVID-19 Data.

Model Type	Peak Inf. Time (days, approx.)	Max Infected (Imax, approx.)	Epidemic Duration (days, approx.)	Alignment with Real Data
Integer-order SIR [23,29]	35	48	∼90	Moderate
Fractional SIR (Weighted Hattaf, $\zeta =0.95$) Figure 13 and Figure 14	∼50	∼40	∼130	High
Actual COVID-19 Trend [36]	50–55	∼39–42	∼130–140	—

References

1. Marinov, T.T.; Marinova, R.S. Adaptive SIR model with vaccination: Simultaneous identification of rates and functions illustrated with COVID-19. Sci. Rep.; 2022; 12, 15688. [DOI: https://dx.doi.org/10.1038/s41598-022-20276-7]

2. Balderrama, R.; Peressutti, J.; Pinasco, J.P.; Vazquez, F.; Vega, C.S.D.L. Optimal control for a SIR epidemic model with limited quarantine. Sci. Rep.; 2022; 12, 12583. [DOI: https://dx.doi.org/10.1038/s41598-022-16619-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35869150]

3. El Khalifi, M.; Britton, T. SIRS epidemics with individual heterogeneity of immunity waning. J. Theor. Biol.; 2022; 587, 111815. [DOI: https://dx.doi.org/10.1016/j.jtbi.2024.111815]

4. Podlubny, I. Fractional Differential Equations; Academic Press: San Diego, CA, USA, 1999.

5. Samko, S.G.; Kilbas, A.A.; Marichev, O.I. Fractional Integrals and Derivatives; Gordon & Breach: Yverdon, Switzerland, 1993.

6. Magin, R.L. Fractional Calculus in Bioengineering; Begell House: Redding, CA, USA, 2006.

7. Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations. North-Holland Mathematics Studies; Elsevier: Amsterdam, The Netherlands, 2006.

8. Alqahtani, R.T. Mathematical model of SIR epidemic system (COVID-19) with fractional derivative: Stability and numerical analysis. Adv. Differ. Equ.; 2021; 2021, 2. [DOI: https://dx.doi.org/10.1186/s13662-020-03192-w] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33424955]

9. Kim, J. Influence of fractional order on the behavior of a normalized time-fractional SIR model. Mathematics; 2024; 12, 3081. [DOI: https://dx.doi.org/10.3390/math12193081]

10. Riabi, L.; Hamdi Cherif, M.; Cattani, C. An Efficient Approach to Solving the Fractional SIR Epidemic Model with the Atangana–Baleanu-Caputo Fractional Operator. Fractal Fract.; 2023; 7, 618. [DOI: https://dx.doi.org/10.3390/fractalfract7080618]

11. Alazman, I.; Mishra, M.N.; Alkahtani, B.S.; Dubey, R.S. Analysis of infection and diffusion coefficient in an sir model by using generalized fractional derivative. Fractal Fract.; 2024; 8, 537. [DOI: https://dx.doi.org/10.3390/fractalfract8090537]

12. Hilfer, R. Applications of Fractional Calculus in Physics; World Scientific: Singapore, 2000.

13. Das, S. Application of generalized fractional calculus in electrical circuit analysis and electromagnetics. Functional Fractional Calculus; Springer: Berlin/Heidelberg, Germany, 2011; pp. 387-436.

14. Oldham, K.B. Fractional differential equations in electrochemistry. Adv. Eng. Softw.; 2010; 41, pp. 9-12. [DOI: https://dx.doi.org/10.1016/j.advengsoft.2008.12.012]

15. Chávez-Vázquez, S.; Gómez-Aguilar, J.F.; Lavín-Delgado, J.E.; Escobar-Jiménez, R.F.; Olivares-Peregrino, V.H. Applications of fractional operators in robotics: A review. J. Intell. Robot. Syst.; 2022; 104, 63. [DOI: https://dx.doi.org/10.1007/s10846-022-01597-1]

16. Ali, Z.; Nia, S.N.; Rabiei, F.; Shah, K.; Tan, M.K. A Semianalytical Approach to the Solution of Time-Fractional Navier-Stokes Equation. Adv. Math. Phys.; 2021; 2021, 5547804. [DOI: https://dx.doi.org/10.1155/2021/5547804]

17. Tarasov, V.E. Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media; Springer Science and Business Media: Berlin, Germany, 2011.

18. Matušů, R. Application of fractional-order calculus to control theory. Int. J. Math. Model. Methods Appl. Sci.; 2011; 5, pp. 1162-1169.

19. Ali, Z. Theoretical and Computational Study of Fractional-Order Mathematical Models for Infectious Diseases. Doctoral Dissertation; Monash University: Melbourne, Australia, 2023.

20. Aldwoah, K.A.; Almalahi, M.A.; Shah, K. Theoretical and numerical simulations on the hepatitis B virus model through a piecewise fractional order. Fractal Fract.; 2023; 7, 844. [DOI: https://dx.doi.org/10.3390/fractalfract7120844]

21. Alraqad, T.; Almalahi, M.A.; Mohammed, N.; Alahmade, A.; Aldwoah, K.A.; Saber, H. Modeling Ebola Dynamics with a Φ-Piecewise Hybrid Fractional Derivative Approach. Fractal Fract.; 2024; 8, 596. [DOI: https://dx.doi.org/10.3390/fractalfract8100596]

22. Saber, H.; Almalahi, M.A.; Albala, H.; Aldwoah, K.; Alsulami, A.; Shah, K.; Moumen, A. Investigating a Nonlinear Fractional Evolution Control Model Using W-Piecewise Hybrid Derivatives: An Application of a Breast Cancer Model. Fractal Fract.; 2024; 8, 735. [DOI: https://dx.doi.org/10.3390/fractalfract8120735]

23. Alshammari, F.S.; Khan, M.A. Dynamic behaviors of a modified SIR model with nonlinear incidence and recovery rates. Alex. Eng. J.; 2021; 60, pp. 2997-3005. [DOI: https://dx.doi.org/10.1016/j.aej.2021.01.023]

24. Lotfi, E.M.; Zine, H.; Torres, D.F.; Yousfi, N. The power fractional calculus: First definitions and properties with applications to power fractional differential equations. Mathematics; 2022; 10, 3594. [DOI: https://dx.doi.org/10.3390/math10193594]

25. Caputo, M.; Fabrizio, M. A new definition of fractional derivative without singular kernel. Prog. Fract. Differ. Appl.; 2015; 1, pp. 73-85.

26. Atangana, A.; Baleanu, D. New fractional derivatives with non-local and non-singular kernel: Theory and application to heat transfer model. Therm. Sci.; 2016; 20, pp. 763-769. [DOI: https://dx.doi.org/10.2298/TSCI160111018A]

27. Al-Refai, M. On weighted Atangana–Baleanu fractional operators. Adv. Differ. Equ.; 2020; 2020, 3. [DOI: https://dx.doi.org/10.1186/s13662-019-2471-z]

28. Hattaf, K. A new generalized definition of fractional derivative with non-singular kernel. Computation; 2020; 8, 49. [DOI: https://dx.doi.org/10.3390/computation8020049]

29. Chawla, S.R.; Ahmad, S.; Khan, A.; Albalawi, W.; Nisar, K.S.; Ali, H.M. Stability analysis and optimal control of a generalized SIR epidemic model with harmonic mean type of incidence and nonlinear recovery rates. Alex. Eng. J.; 2024; 97, pp. 44-60. [DOI: https://dx.doi.org/10.1016/j.aej.2024.04.017]

30. Shan, C.; Zhu, H. Bifurcations and complex dynamics of an SIR model with the impact of the number of hospital beds. J. Differ. Equ.; 2014; 257, pp. 1662-1688. [DOI: https://dx.doi.org/10.1016/j.jde.2014.05.030]

31. Zitane, H.; Torres, D.F. A class of fractional differential equations via power non-local and non-singular kernels: Existence, uniqueness and numerical approximations. Phys. D Nonlinear Phenom.; 2024; 457, 133951. [DOI: https://dx.doi.org/10.1016/j.physd.2023.133951]

32. Gassem, F.; Almalahi, M.; Osman, O.; Muflh, B.; Aldwoah, K.; Kamel, A.; Eljaneid, N. Nonlinear Fractional Evolution Control Modeling via Power Non-Local Kernels: A Generalization of Caputo–Fabrizio, Atangana–Baleanu, and Hattaf Derivatives. Fractal. Fract.; 2025; 9, 104.

33. Lipsitch, M.; Krammer, F.; Regev-Yochay, G.; Lustig, Y.; Balicer, R.D. SARS-CoV-2 breakthrough infections in vaccinated individuals: Measurement, causes and impact. Nat. Rev. Immunol.; 2022; 22, pp. 57-65. [DOI: https://dx.doi.org/10.1038/s41577-021-00662-4]

34. Cromer, D.; Steain, M.; Reynaldi, A.; Schlub, T.E.; Wheatley, A.K.; Juno, J.A.; Kent, S.J.; Triccas, J.A.; Khoury, D.S.; Davenport, M.P. Neutralizing antibody titres as predictors of protection against SARS-CoV-2 variants and the impact of boosting: A meta-analysis. Lancet Microbe; 2022; 3, pp. e52-e61. [DOI: https://dx.doi.org/10.1016/S2666-5247(21)00267-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34806056]

35. Verity, R.; Okell, L.C.; Dorigatti, I.; Winskill, P.; Whittaker, C.; Imai, N.; Ferguson, N.M. Estimates of the severity of COVID-19 disease. Lancet Infect. Dis.; 2020; 20, pp. 669-677. [DOI: https://dx.doi.org/10.1016/S1473-3099(20)30243-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32240634]

36. World Health Organization. WHO COVID-19 Dashboard-WHO Data. Available online: https://data.who.int/dashboards/covid19/cases (accessed on 23 March 2025).