A Comparative Study of a Fractional‐Order Rabies Model With MCMC Estimation Using Singular and Nonsingular Operators

Abstract

Rabies remains a significant public health concern, particularly in regions with high dog‐mediated transmission, and understanding its dynamics is crucial for effective control strategies. This study investigates the transmission dynamics of rabies by developing a deterministic human‐dog model extended to fractional‐order derivatives, incorporating three operators: Caputo, Caputo–Fabrizio (CF), and Atangana–Baleanu–Caputo (ABC), to capture memory and hereditary effects. Model parameters were estimated from field data using the Markov chain Monte Carlo (MCMC) method, and the effective reproduction number, ${\mathcal{R}}_e$, was derived via a graph‐theory approach. Mathematical analysis establishes the positivity, boundedness, and stability of solutions. Comparative simulations indicate that fractional‐order models capture slower disease progression compared to classical integer‐order systems, with the ABC operator producing the most conservative epidemic projections, reflecting realistic epidemic inertia. The study highlights the critical impact of vaccination, culling, and postexposure prophylaxis (PEP) in controlling rabies. The novelty of this work lies in the comprehensive comparison of different fractional‐order operators within the same modeling framework, providing new insights into the role of memory effects in rabies transmission and guiding more effective intervention strategies.

1. Introduction

Rabies is a highly fatal viral disease caused by the rabies virus, primarily affecting mammals, with significant public health implications, particularly in Africa and Asia, where it remains endemic [1]. Annually, rabies results in approximately 59,000 human deaths globally, with domestic dogs identified as the primary vectors for transmission to humans [2–4]. The zoonotic nature of rabies underscores its threat, as rabid dogs frequently bite humans, leading to severe neurological symptoms and almost certain mortality once clinical signs appear [5, 6]. This situation is exacerbated by the close interactions between humans and dogs, especially in regions with high pet ownership, which increases the risk of zoonotic disease transmission [2, 4].

The urgency for effective control strategies against rabies is required, given its high impact on public health and the economy. Vaccination emerges as the most reliable method for preventing rabies, not only in dogs but also in humans, particularly in rabies-endemic areas [3, 4]. Vaccination programs targeting domestic dogs are crucial, as they significantly reduce the incidence of rabies in human populations [5, 6]. Furthermore, sustained international commitment is essential for the global elimination of rabies, focusing on comprehensive vaccination strategies and public awareness campaigns to mitigate the spread of the disease [2, 5].

Mathematical models play a crucial role in understanding the transmission dynamics of infectious diseases by providing frameworks to analyze and predict disease spread, evaluate control strategies, and inform public health interventions. These models simplify complex real-world phenomena, allowing researchers and planners to estimate key epidemiological tools like reproduction numbers, which are essential for forecasting disease trajectories and assessing the impact of interventions [7]. They also facilitate the exploration of transmission uncertainties and the effectiveness of various control measures, as demonstrated in studies of metapopulation dynamics where individual movement between patches influences infection spread [8]. Despite their limitations, mathematical models remain invaluable tools for guiding healthcare resources and shaping policy decisions in the face of infectious disease outbreaks [9].

Classical models in biological systems, particularly those utilizing ordinary differential equations (ODEs), often struggle to encapsulate the complexities inherent in real-world phenomena such as disease progression and immune responses. While ODEs can effectively describe certain dynamics, they typically fail to account for critical factors like memory and hereditary effects, which are essential in understanding the nuanced behavior of biological systems [10]. Moreover, the limitations of ODEs have prompted a shift towards more sophisticated modeling techniques, including agent-based models and microscopic simulations, which better capture the emergent behaviors and inhomogeneities present in biological systems [11]. This evolution reflects a growing recognition of the need for models that can incorporate the complexities of biological realities and, thereby enhancing predictive accuracy.

Fractional calculus extends classical calculus by allowing derivatives and integrals to be of noninteger order, providing a more generalized framework for mathematical analysis and simulation. This concept, which dates back to the correspondence between Leibniz and L’Hopital in the late 17th century, has evolved significantly over the centuries with contributions from numerous mathematicians such as Laplace, Fourier, and Riemann, among others [12]. The essence of fractional calculus lies in its ability to model systems with memory, nonlocality, and hereditary properties, which are not adequately captured by traditional integer-order calculus [12]. This makes it particularly suitable for modeling complex systems where past states influence current dynamics, such as in the case of infectious disease modeling. For instance, fractional-order derivatives have been applied to extend classical epidemiological models to incorporate memory effects, thereby providing a more realistic representation of disease transmission dynamics [13]. The flexibility and depth offered by fractional calculus make it a powerful tool in various fields, particularly in scenarios where traditional calculus falls short in capturing the intricacies of real-world phenomena [14].

The study of fractional operators reveals two major types: singular kernel operators and nonsingular kernel operators. Singular kernel operators, such as the Caputo fractional derivative, exhibit power–law memory and are characterized by singularities in their kernels, which can complicate numerical stability and convergence in certain applications [15, 16]. In contrast, nonsingular kernel operators, including the Caputo–Fabrizio (CF) and Atangana–Baleanu derivatives, are designed to provide finite memory effects, enhancing numerical stability and realism in modeling biological systems [17, 18]. These nonsingular operators utilize functions like the Mittag–Leffler function, which facilitate easier computation and better performance in solving fractional differential equations, making them preferable for practical applications [19]. Overall, the choice between these operators significantly impacts the modeling of systems with memory effects.

The modeling of rabies using fractional order derivatives, specifically the Caputo, CF, and Atangana–Baleanu derivatives, has gained traction in recent research due to their ability to capture the complexities of disease dynamics more effectively than traditional integer-order models. Aydogan et al. [20] utilized the CF derivative to establish a unique solution for rabies dynamics, demonstrating the impact of fractional derivatives on the stability and behavior of infected populations. Zainab et al. [21] developed a model using the Atangana–Baleanu derivative, emphasizing its predictive capabilities for livestock health management. Furthermore, Zarin et al. [22] explored optimal control strategies within a fractional framework, highlighting the effectiveness of these models in minimizing rabies spread through various control measures. Garg and Chauhan [23] focused on the stability of solutions using a fixed-point approach, reinforcing the importance of fractional modeling in understanding rabies transmission. Collectively, these studies illustrate the advantages of fractional calculus in epidemiological modeling, offering deeper insights into disease control and prevention strategies.

Recent studies have highlighted the significance of fractional-order derivatives in modeling complex dynamical systems, particularly in epidemiology and ecological interactions. The CF derivative and its generalizations have been applied to electrical circuits [24] and epidemiological models [25, 26], demonstrating their ability to capture memory and hereditary effects more effectively than classical integer-order derivatives. Fractional models with singular and nonsingular kernels have also been employed in host–parasitoid [27], predator–prey–pathogen [28], and COVID-19 transmission studies [29, 30], showing improved fidelity in simulating real-world dynamics. Recent works further emphasize the importance of nonsingular kernels for incorporating vaccination efficacy and intervention strategies [26]. In parallel, fractional calculus has been leveraged to explore soliton wave propagation, Lie symmetry analysis, and invariant structures in physical systems [31–38], providing robust mathematical frameworks for analyzing nonlinear dynamical behavior. Collectively, these studies motivate the current work in employing both singular and nonsingular fractional operators to analyze rabies transmission, offering a comparative perspective that captures memory effects and improves predictive modeling.

Markov chain Monte Carlo (MCMC) estimation plays a crucial role in modeling infectious diseases by enabling the analysis of complex and stochastic epidemic models that account for uncertainties and missing data inherent in real-world scenarios. MCMC methods facilitate the estimation of model parameters from partially observed data, which is essential given that complete observation of epidemic events is often impractical [39, 40]. These techniques allow for more realistic assumptions in data analysis, improving the accuracy of parameter estimates and forecasts, particularly when incorporating population-level variations in transmission dynamics [41]. Furthermore, MCMC approaches, such as the Metropolis–Hastings algorithm and Gibbs sampler, have demonstrated effectiveness in inferring parameters from diverse datasets, including those from household outbreaks of diseases like measles and influenza [42]. Overall, MCMC estimation enhances the robustness of epidemic modeling, providing valuable insights for public health interventions and outbreak management.

While numerous mathematical models have examined rabies transmission and control through vaccination, most have relied on classical integer-order differential equations, which often fail to capture the memory and hereditary properties of disease progression and host immune response. Fractional-order models, incorporating singular and nonsingular operators such as the Caputo, CF, and Atangana–Baleanu derivatives, provide a more flexible and accurate framework for modeling such dynamics. Despite their potential, a systematic comparative study evaluating the effects of these different fractional operators on rabies transmission and vaccination outcomes remains lacking. To address this gap, the present study aims to assess the role of memory effects in shaping rabies dynamics and optimizing control strategies. In order to calibrate the model accurately, we employ the MCMC method, a Bayesian approach that generates samples from the posterior distribution of parameters, allowing incorporation of prior knowledge and quantification of uncertainty. By using MCMC, robust parameter estimates with credible intervals are obtained, ensuring that the stochastic variability inherent in rabies transmission is adequately captured and that model projections are closely aligned with observed data.

2. Preliminaries for Fractional Operators

This study considers three types of fractional operators: Caputo (C), Caputo–Fabrizio (CF), and Atangana–Baleanu–Caputo (ABC). The fractional derivative of order θ is denoted by ${}^i{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)$, and the corresponding fractional integral by ${}^i{\mathcal{I}}_t^{\theta }\mathcal{X}\left(t\right)$, where i = C, CF, or ABC refers to the respective operator. These operators are designed to capture memory and hereditary properties in dynamical systems, with varying types of kernels (singular or nonsingular).

$$

Statement

Definition 1 Gamma function [43].

The Gamma function Γ(⋅) is defined for x > 0 by: 1 $\begin{array}{}\Gamma \left(x\right)={\int }_0^{\mathrm{\infty }}{t}^{x-1}{e}^{-t}dt,\end{array}$ and generalizes the factorial function via Γ(n) = (n − 1)! for natural numbers n.

$$

Statement

Definition 2 Mittag–Leffler function [43].

The one-parameter Mittag–Leffler function is defined for $z\in \mathbb{C}$, θ > 0, as: 2 $\begin{array}{}{E}_{\theta }\left(z\right)=\sum _{k=0}^{\mathrm{\infty }}\frac{z^k}{\Gamma \left(\theta k+1\right)}.\end{array}$

$$

Statement

Definition 3 Caputo fractional derivative [44].

The Caputo fractional derivative of order θ > 0 for a function $\mathcal{X}\in C^n\left[0,t\right]$ is defined as: 3 $\begin{array}{}{}^C{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)=\frac{1}{\Gamma \left(n-\theta \right)}{\int }_0^t\frac{{\mathcal{X}}^{\left(n\right)}\left(\varphi \right)}{{\left(t-\varphi \right)}^{\theta -n+1}}d\varphi ,\end{array}$ where n = ⌈θ⌉. This operator includes a singular kernel and is widely used due to its compatibility with classical initial conditions.

$$

Statement

Definition 4 Caputo fractional integral [44].

Let $\mathcal{X}:{\mathbb{R}}^{+}\to \mathbb{R}$ be a locally integrable function. The Caputo fractional integral of order θ > 0 is given by: 4 $\begin{array}{}{}^C{\mathcal{I}}_t^{\theta }\mathcal{X}\left(t\right)=\frac{1}{\Gamma \left(\theta \right)}{\int }_0^t{\left(t-\varphi \right)}^{\theta -1}\mathcal{X}\left(\varphi \right)d\varphi .\end{array}$

The Caputo derivative uses a singular power–law kernel (t − ϕ)^−β, which captures long-memory effects, but may lead to numerical instability. To overcome this, regular kernels have been proposed, such as in the CF and ABC models.

$$

Statement

Definition 5 CF fractional derivative [45].

The CF derivative of order 0 < θ < 1 for a function $\mathcal{X}\in AC\left(\left[0,T\right]\right)$ is defined as: 5 $\begin{array}{}{}^{CF}{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)=\frac{M\left(\theta \right)}{1-\theta }{\int }_0^t{\mathcal{X}}^{\prime }\left(\varphi \right)e^{-\frac{\theta }{1-\theta }\left(t-\varphi \right)}d\varphi ,\end{array}$ where M(θ) is a normalization constant satisfying M(0) = M(1) = 1.

$$

Statement

Definition 6 CF fractional integral [45].

The CF fractional integral of order 0 < θ < 1 for a function $\mathcal{X}\in L^1\left(\left[0,T\right]\right)$ is given by: 6 $\begin{array}{}{}^{CF}{\mathcal{I}}_t^{\theta }\mathcal{X}\left(t\right)=\frac{1-\theta }{M\left(\theta \right)}\mathcal{X}\left(t\right)+\frac{\theta }{M\left(\theta \right)}{\int }_0^t\mathcal{X}\left(\varphi \right)d\varphi .\end{array}$

$$

Statement

Definition 7 Atangana–Baleanu–Caputo (ABC) fractional derivative [46].

The ABC fractional derivative of order 0 < θ < 1 for a function $\mathcal{X}\in H^1\left(\left[0,T\right]\right)$ is defined as: 7 $\begin{array}{}{}^{\text{ABC}}{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)=\frac{B\left(\theta \right)}{1-\theta }{\int }_0^t{\mathcal{X}}^{\prime }\left(\varphi \right)E_{\theta }\left(-\frac{\theta }{1-\theta }{\left(t-\varphi \right)}^{\theta }\right)d\varphi ,\end{array}$ where B(θ) is a normalization constant with B(0) = B(1) = 1, and E_θ(⋅) is the one-parameter Mittag–Leffler function as defined in Definition 2.

$$

Statement

Definition 8 Atangana–Baleanu–Caputo (ABC) fractional integral [46].

The ABC fractional integral of order 0 < θ < 1 for a function $\mathcal{X}\in L^1\left(\left[0,T\right]\right)$ is given by: 8 $\begin{array}{}{}^{\text{ABC}}{\mathcal{I}}_t^{\theta }\mathcal{X}\left(t\right)=\frac{1-\theta }{B\left(\theta \right)}\mathcal{X}\left(t\right)+\frac{\theta }{B\left(\theta \right)\Gamma \left(\theta \right)}{\int }_0^t{\left(t-\varphi \right)}^{\theta -1}\mathcal{X}\left(\varphi \right)d\varphi .\end{array}$

$$

Statement

Definition 9 Fractional initial value problem (FIVP).

A general FIVP is expressed as: 9 $\begin{array}{}{}^i{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)=f\left(t,\mathcal{X}\left(t\right)\right),\mathcal{X}\left(0\right)={\mathcal{X}}_0,\end{array}$ where ${}^i{\mathcal{D}}_t^{\theta }$ denotes the fractional derivative of Caputo, CF, or ABC type and f is a continuous function describing the system dynamics.

3. Model Formulation

In Section 3, we present an extension of the deterministic mathematical model for the transmission dynamics of rabies, developed by Ndendya et al. [47], by incorporating two interacting populations: humans and domestic dogs.

3.1. Model Assumptions

The model is built on the following biological and epidemiological assumptions:

• i.

  Rabies transmission occurs only from infectious dogs to humans and dogs, meaning that there is no human-to-human or human-to-dog transmission, which reflects the known zoonotic nature of rabies in endemic areas where dogs are the primary reservoir.

• ii.

  Only susceptible dogs are vaccinated at a constant rate and move to a protected class, from which it provide permanent immunity.

• ii.

  No recovery occurs in either humans or dogs; infection with rabies is invariably fatal. Once an individual becomes infectious, they do not recover, but ultimately die due to rabies-induced mortality.

3.2. Model Description

Based on these assumptions, the model describes the dynamics of rabies in two interconnected populations: humans and domestic dogs. The human population is subdivided into three compartments: susceptible S_h, exposed E_h, and infectious I_h. Individuals enter the human population at rate π_h, and susceptible humans may become exposed through contact with infectious dogs. The force of infection for humans is given by ${\lambda }_{\text{h}}=\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}}{N_{\text{d}}}$, where σ ∈ [0, 1] represents the protective effect of public health education, and N_d is the total dog population. The parameter ω represent the exposed humans who do not develop rabies due to effectiveness of postexposure prophylaxis (PEP) and return to the susceptible class. Also exposed humans progress to the infectious stage at rate τ_h or die naturally at rate ɛ_h, while infectious humans exit the system due to natural or rabies-induced death.

The dog population is structured into susceptible S_d, exposed E_d, infectious I_d, and vaccinated V_d compartments. Susceptible dogs are recruited at rate π_d, and may become exposed upon contact with infectious dogs, with infection occurring at rate ${\lambda }_{\text{d}}={\beta }_{\text{d}}\frac{I_{\text{d}}}{N_{\text{d}}}$. Susceptible dogs may also be vaccinated at rate γ and move to the vaccinated class, or die naturally at rate ɛ_d. Exposed dogs become infectious at rate τ_d, and infectious dogs are removed through natural death, rabies-induced mortality, or culling at rate μ. Because the vaccine provide full immunity that means vaccinated dogs exit the system only through natural death. Based on the assumptions and model description, the conceptual framework is illustrated in Figure 1, the corresponding system of differential equations is presented in Equation (10) and parameter description in Table 1.

[IMAGE OMITTED. SEE PDF]

Table 1 Model parameter values, descriptions, and sources.

Parameter	Description	Value	Source
πh	Human recruitment (birth/immigration) rate	50	[47]
πd	Dog recruitment rate	40	[47]
βh	Transmission rate from infectious dogs to humans	0.05372	Estimated
βd	Dog-to-dog transmission rate	0.12250	Estimated
τh	Progression rate from exposed to infectious in humans	0.05430	Estimated
τd	Progression rate from exposed to infectious in dogs	0.12000	Estimated
γ	Vaccination rate of susceptible dogs	0.01984	Estimated
αh	Rabies-induced death rate in humans	0.01144	Estimated
αd	Rabies-induced death rate in dogs	0.05230	Estimated
μ	Culling/removal rate of infectious dogs	0.05153	Estimated
ɛh	Natural death rate of humans	0.0001	[47]
ɛd	Natural death rate of dogs	0.8	[47]
ω	Rate of loss of latency in humans (progress to susceptible)	0.21334	[47]
σ	Public health education efficacy	0.8826	[47]

10 $\begin{array}{}\begin{array}{c}\left\{\begin{array}{cc}\frac{d{S}_{\text{h}}}{dt}\hfill & ={\pi }_{\text{h}}+\omega E_{\text{h}}-\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}}{N_{\text{d}}}{S}_{\text{h}}-{\varepsilon }_{\text{h}}{S}_{\text{h}}\\ \frac{d{E}_{\text{h}}}{dt}\hfill & =\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}}{N_{\text{d}}}{S}_h-\omega E_{\text{h}}-{\tau }_{\text{h}}{E}_{\text{h}}-{\varepsilon }_{\text{h}}{E}_{\text{h}}\\ \frac{d{I}_{\text{h}}}{dt}\hfill & ={\tau }_{\text{h}}{E}_{\text{h}}-\left({\varepsilon }_{\text{h}}+{\alpha }_{\text{h}}\right)I_{\text{h}}\\ \frac{d{S}_{\text{d}}}{dt}\hfill & ={\pi }_{\text{d}}-{\beta }_{\text{d}}\frac{I_{\text{d}}}{N_{\text{d}}}{S}_{\text{d}}-\gamma S_{\text{d}}-{\varepsilon }_{\text{d}}{S}_{\text{d}}\\ \frac{d{E}_{\text{d}}}{dt}\hfill & ={\beta }_{\text{d}}\frac{I_{\text{d}}}{N_{\text{d}}}{S}_{\text{d}}-{\tau }_{\text{d}}{E}_{\text{d}}-{\varepsilon }_{\text{d}}{E}_{\text{d}}\\ \frac{d{I}_{\text{d}}}{dt}\hfill & ={\tau }_{\text{d}}{E}_{\text{d}}-\left({\varepsilon }_{\text{d}}+{\alpha }_{\text{d}}+\mu \right)I_{\text{d}}\\ \frac{d{V}_{\text{d}}}{dt}\hfill & =\gamma S_{\text{d}}-{\varepsilon }_{\text{d}}{V}_{\text{d}}\end{array}\right.\end{array},\end{array}$ with initial conditions satisfy: 11 $\begin{array}{}\begin{array}{c}{S}_{\text{h}}\left(0\right)>0000000.,E_{\text{h}}\left(0\right)\ge ,I_{\text{h}}\left(0\right)>,S_{\text{d}}\left(0\right)>,E_{\text{d}}\left(0\right)\ge ,I_{\text{d}}\left(0\right)>,\text{and\hspace{0.17em}}{V}_{\text{d}}\left(0\right)\ge \end{array}\end{array}$

3.3. Rabies Fractional-Order Model

Here, we developed a fractional-order rabies transmission model based on the deterministic system (Equation (10)), employing the Caputo, CF, and ABC derivative to capture memory effects inherent in disease dynamics. The model accounts for interactions among humans and dogs. 12 $\begin{array}{}\begin{array}{c}\left\{\begin{array}{cc}{}_0^i{\mathcal{D}}_t^{\theta }{S}_{\text{h}}\hfill & ={\pi }_{\text{h}}^{\theta }+{\omega }^{\theta }{E}_{\text{h}}-\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }\frac{I_{\text{d}}}{N_{\text{d}}}{S}_{\text{h}}-{\varepsilon }_{\text{h}}^{\theta }{S}_{\text{h}},\\ {}_0^i{\mathcal{D}}_t^{\theta }{E}_{\text{h}}\hfill & =\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }\frac{I_{\text{d}}}{N_{\text{d}}}{S}_{\text{h}}-{\omega }^{\theta }{E}_{\text{h}}-{\tau }_{\text{h}}^{\theta }{E}_{\text{h}}-{\varepsilon }_{\text{h}}^{\theta }{E}_{\text{h}},\\ {}_0^i{\mathcal{D}}_t^{\theta }{I}_{\text{h}}\hfill & ={\tau }_{\text{h}}^{\theta }{E}_{\text{h}}-\left({\varepsilon }_{\text{h}}^{\theta }+{\alpha }_{\text{h}}^{\theta }\right)I_h,\\ {}_0^i{\mathcal{D}}_t^{\theta }{S}_{\text{d}}\hfill & ={\pi }_{\text{d}}^{\theta }-{\beta }_{\text{d}}^{\theta }\frac{I_{\text{d}}}{N_{\text{d}}}{S}_{\text{d}}-{\gamma }^{\theta }{S}_{\text{d}}-{\varepsilon }_{\text{d}}^{\theta }{S}_{\text{d}},\\ {}_0^i{\mathcal{D}}_t^{\theta }{E}_{\text{d}}\hfill & ={\beta }_{\text{d}}^{\theta }\frac{I_{\text{d}}}{N_{\text{d}}}{S}_{\text{d}}-{\tau }_{\text{d}}^{\theta }{E}_{\text{d}}-{\varepsilon }_{\text{d}}^{\theta }{E}_{\text{d}},\\ {}_0^i{\mathcal{D}}_t^{\theta }{I}_{\text{d}}\hfill & ={\tau }_{\text{d}}^{\theta }{E}_{\text{d}}-\left({\varepsilon }_{\text{d}}^{\theta }+{\alpha }_{\text{d}}^{\theta }+{\mu }^{\theta }\right)I_{\text{d}},\\ {}_0^i{\mathcal{D}}_t^{\theta }{V}_{\text{d}}\hfill & ={\gamma }^{\theta }{S}_{\text{d}}-{\varepsilon }_{\text{d}}^{\theta }{V}_{\text{d}},\end{array}\right.\end{array},\end{array}$ with initial conditions satisfying: 13 $\begin{array}{}\begin{array}{c}{S}_{\text{h}}\left(0\right)>0000000,E_{\text{h}}\left(0\right)\ge ,I_{\text{h}}\left(0\right)>,S_{\text{d}}\left(0\right)>,E_{\text{d}}\left(0\right)\ge ,I_{\text{d}}\left(0\right)>,V_{\text{d}}\left(0\right)\ge ,\text{and\hspace{0.17em}}\theta \in \left(01,\right],\end{array}\end{array}$ where i denotes C, CF, and ABC, which stand for the Caputo, CF, and ABC derivatives, respectively.

4. Mathematical Analysis of the Fractional-Order Model

Section 4 presents a rigorous mathematical analysis of the fractional-order rabies transmission model. We establish fundamental properties to ensure the model’s epidemiological validity and mathematical well-posedness.

4.1. Positivity and Boundedness of Solutions

We employ the methodology developed by Agarwal et al. [48] to demonstrate the positivity and boundedness of solutions for the fractional rabies model. Consider the feasible region: 14 $\begin{array}{}\Lambda ={\Lambda }_{\text{h}}\times {\Lambda }_{\text{d}}\subset {\mathbb{R}}_{+}^7,\end{array}$ where the human and dog population subspaces are, respectively, defined as: 15 $\begin{array}{}\begin{array}{c}{\Lambda }_{\text{h}}=\left\{\left(S_{\text{h}}\left(t\right),E_{\text{h}}\left(t\right),I_{\text{h}}\left(t\right)\right)\in {\mathbb{R}}_{+}^3:N_{\text{h}}\left(t\right)\le \frac{{\pi }_{\text{h}}^{\theta }}{{\varepsilon }_{\text{h}}^{\theta }}\right\},\end{array}\end{array}$ 16 $\begin{array}{}\begin{array}{cc}{\Lambda }_{\text{d}}& =\left\{\left(S_{\text{d}}\left(t\right),V_{\text{d}}\left(t\right),E_{\text{d}}\left(t\right),I_{\text{d}}\left(t\right)\right)\in {\mathbb{R}}_{+}^4:N_{\text{d}}\left(t\right)\le \frac{{\pi }_{\text{d}}^{\theta }}{{\varepsilon }_{\text{d}}^{\theta }}\right\}.\end{array}\end{array}$

$$

Statement

Lemma 1 Positive invariance.

The closed set Λ is positively invariant for the fractional system (Equation (12)).

$$

Statement

Proof.

The proof proceeds in three main steps:

Step 1: Human population dynamics: The total human population satisfies: 17 $\begin{array}{}{}_0^i{D}_t^{\theta }{N}_{\text{h}}\left(t\right)={\pi }_{\text{h}}^{\theta }-{\varepsilon }_{\text{h}}^{\theta }{N}_{\text{h}}\left(t\right).\end{array}$

Applying the Laplace transform to Equation (17) yields: 18 $\begin{array}{}s\mathcal{L}\left\{N_{\text{h}}\left(t\right)\right\}-s^{\theta -1}{N}_{\text{h}}\left(0\right)=\frac{{\pi }_{\text{h}}^{\theta }}{s}-{\varepsilon }_{\text{h}}^{\theta }\mathcal{L}\left\{N_{\text{h}}\left(t\right)\right\}.\end{array}$

Solving for $\mathcal{L}\left\{N_{\text{h}}\left(t\right)\right\}$ gives: 19 $\begin{array}{}\mathcal{L}\left\{N_{\text{h}}\left(t\right)\right\}=\frac{{\pi }_{\text{h}}^{\theta }}{s\left(s+{\varepsilon }_{\text{h}}^{\theta }\right)}+\frac{s^{\theta -1}{N}_{\text{h}}\left(0\right)}{s+{\varepsilon }_{\text{h}}^{\theta }}.\end{array}$

The inverse Laplace transform yields the solution: 20 $\begin{array}{}{N}_{\text{h}}\left(t\right)=\frac{{\pi }_{\text{h}}^{\theta }}{{\varepsilon }_{\text{h}}^{\theta }}\left[1-E_{\theta }\left(-{\varepsilon }_{\text{h}}^{\theta }t\right)\right]+N_{\text{h}}\left(0\right)E_{\theta }\left(-{\varepsilon }_{\text{h}}^{\theta }t\right),\end{array}$ where E_θ(z) is the Mittag–Leffler function.

Step 2: Dog population dynamics

Similarly, the dog population follows: 21 $\begin{array}{}{N}_{\text{d}}\left(t\right)=\frac{{\pi }_{\text{d}}^{\theta }}{{\varepsilon }_{\text{d}}^{\theta }}\left[1-E_{\theta }\left(-{\varepsilon }_{\text{d}}^{\theta }t\right)\right]+N_{\text{d}}\left(0\right)E_{\theta }\left(-{\varepsilon }_{\text{d}}^{\theta }t\right).\end{array}$

Step 3: Invariance of Λ

From Equations (20) and (21), we observe that:

• •

  If $N_{\text{h}}\left(0\right)\le \frac{{\pi }_{\text{h}}^{\theta }}{{\varepsilon }_{\text{h}}^{\theta }}$, then $N_{\text{h}}\left(t\right)\le \frac{{\pi }_{\text{h}}^{\theta }}{{\varepsilon }_{\text{h}}^{\theta }}$ for all t > 0.

• •

  If $N_{\text{d}}\left(0\right)\le \frac{{\pi }_{\text{d}}^{\theta }}{{\varepsilon }_{\text{d}}^{\theta }}$, then $N_{\text{d}}\left(t\right)\le \frac{{\pi }_{\text{d}}^{\theta }}{{\varepsilon }_{\text{d}}^{\theta }}$ for all t > 0.

Thus, the set Λ is positively invariant under the system dynamics.

4.2. Existence and Uniqueness of Solutions

Section 4.2 establishes the existence and uniqueness of solutions for the fractional-order rabies model. Following Edward [44], we employ the Banach contraction principle to analyze the system’s well-posedness. The model in Caputo form can be expressed as: 22 $\begin{array}{}{}_0^i{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)=\Psi \left(t,\mathcal{X}\left(t\right)\right),t\ge 0,\mathcal{X}\left(0\right)=x_0\in {\mathbb{R}}_{+}^7,\end{array}$ where $\Psi :{\mathbb{R}}_{+}\times {\mathbb{R}}_{+}^7\to {\mathbb{R}}^7$ represents the nonlinear vector field.

$$

Statement

Lemma 2 Existence and uniqueness conditions [44].

The system (Equation (22)) admits a unique solution on [0, ∞) if Ψ satisfies:

• i.

  Continuity in t for all $\mathcal{X}\in {\mathbb{R}}_{+}^7$.

• ii.

  Lipschitz continuity in $\mathcal{X}$ with continuous partial derivatives $\frac{\partial \Psi }{\partial \mathcal{X}}$.

• iii.

  Linear growth condition: $‖\Psi \left(t,\mathcal{X}\right)‖\le p_1+p_2‖\mathcal{X}‖$ for constants p₁, p₂ > 0.

$$

Statement

Proof.

We reformulate the Model (12) in vector form: 23 $\begin{array}{}{}_0^i{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)={\mathcal{P}}_2+{\mathcal{P}}_1\mathcal{X}\left(t\right)+\mathcal{N}\left(\mathcal{X}\left(t\right)\right),\end{array}$ where:

• •

  State vector: $\mathcal{X}\left(t\right)={\left(S_{\text{h}},E_{\text{h}},I_{\text{h}},S_{\text{d}},V_{\text{d}},E_{\text{d}},I_{\text{d}}\right)}^T.$

• •

  Initial condition: ${\mathcal{X}}_0={\left({\pi }_{\text{h}}^{\theta },00000,,,{\pi }_{\text{d}}^{\theta },,\right)}^T.$

• •

  Nonlinear terms: $\mathcal{N}\left(\mathcal{X}\left(t\right)\right)={\left(\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }{S}_{\text{h}}\frac{I_{\text{d}}}{N_{\text{h}}},00000,,{\beta }_{\text{d}}^{\theta }{S}_{\text{d}}\frac{I_{\text{d}}}{N_{\text{d}}},,,\right)}^T.$

The system matrices are: $\begin{array}{}{\mathcal{P}}_1=\left(\begin{array}{ccccccc}-{\varepsilon }_{\text{h}}^{\theta }& {\omega }^{\theta }& 0& 0& 0& 0& 0\\ 0& -{\varepsilon }_{\text{h}}^{\theta }-{\tau }_{\text{h}}^{\theta }-{\omega }^{\theta }& 0& 0& 0& 0& 0\\ 0& {\tau }_{\text{h}}^{\theta }& -{\alpha }_{\text{h}}^{\theta }-{\varepsilon }_{\text{h}}^{\theta }& 0& 0& 0& 0\\ 0& 0& 0& -{\gamma }^{\theta }-{\varepsilon }_{\text{d}}^{\theta }& 0& 0& 0\\ 0& 0& 0& {\gamma }^{\theta }& -{\varepsilon }_{\text{d}}^{\theta }& 0& 0\\ 0& 0& 0& 0& 0& -{\varepsilon }_{\text{d}}^{\theta }-{\tau }_{\text{d}}^{\theta }& 0\\ 0& 0& 0& 0& 0& {\tau }_{\text{d}}^{\theta }& -{\alpha }_{\text{d}}^{\theta }-{\varepsilon }_{\text{d}}^{\theta }-{\mu }^{\theta }\end{array}\right),\end{array}$ $\begin{array}{}{\mathcal{P}}_2=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right),\mathcal{N}\left(\mathcal{X}\right)=\left(\begin{array}{c}-\frac{\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }{S}_{\text{h}}}{N_{\text{d}}}\\ \frac{\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }{S}_{\text{h}}}{N_{\text{d}}}\\ 0\\ -\frac{{\beta }_{\text{d}}^{\theta }{S}_{\text{d}}}{N_{\text{d}}}\\ 0\\ \frac{{\beta }_{\text{d}}^{\theta }{S}_{\text{d}}}{N_{\text{d}}}\\ 0\end{array}\right).\end{array}$

To verify Lemma 2’s conditions:

• i.

  The continuity of Ψ follows immediately from the smoothness of polynomial functions.

• ii.

  The Jacobian $\frac{\partial \Psi }{\partial \mathcal{X}},$ exists and is continuous for all $\mathcal{X}\in {\mathbb{R}}_{+}^7.$

• iii.

  The growth condition holds with: $\begin{array}{}‖\Psi \left(t,\mathcal{X}\right)‖\le ‖{\mathcal{P}}_2‖+‖{\mathcal{P}}_1‖‖\mathcal{X}‖+{\beta }_h^{\theta }‖S_h‖‖I_d‖+{\beta }_d^{\theta }‖S_d‖‖I_d‖.\end{array}$

Thus, by the Banach fixed-point theorem, there exists a unique solution $\mathcal{X}\left(t\right)\in C\left([0,\mathrm{\infty }\right);{\mathbb{R}}_{+}^7)$ to the fractional-order system.

4.3. Effective Reproduction Number ${\mathcal{R}}_e$

The rabies-free equilibrium (RFE) is characterized by the absence of infection in both human and dog populations. This state is given by: $\begin{array}{}\xi =\left(S_{\text{h}}^0,0000,,S_{\text{d}}^0,V_{\text{d}}^0,,\right)=\left(\frac{{\pi }_{\text{h}}^{\theta }}{{\varepsilon }_{\text{h}}^{\theta }},0000,,\frac{{\pi }_{\text{d}}}{{\gamma }^{\theta }+{\varepsilon }_{\text{d}}^{\theta }},\frac{{\gamma }^{\theta }{\pi }_{\text{d}}^{\theta }}{{\varepsilon }_{\text{d}}^{\theta }\left({\gamma }^{\theta }+{\varepsilon }_{\text{d}}^{\theta }\right)},,\right).\end{array}$

The effective reproduction number ${\mathcal{R}}_e$ quantifies disease transmissibility under control measures, representing the expected number of secondary cases produced by a typical infected individual during its infectious period [1]. We compute ${\mathcal{R}}_e$ using the Graph-Theoretic (GT) method [49], which involves:

• i.

  Decomposing the system into new infections ($\mathcal{F}$) and transfer terms ($\mathcal{V}$).

• ii.

  Constructing the next generation matrix $F{\mathcal{R}}_e^{-1}-V.$

• iii.

  Applying GT–based Gaussian elimination to determine the spectral radius.

For our model, the infection and transition terms are: $\begin{array}{}\mathcal{F}=\left(\begin{array}{c}\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }\frac{I_{\text{d}}{S}_{\text{h}}}{N_{\text{d}}}\\ 0\\ {\beta }_{\text{d}}^{\theta }\frac{I_{\text{d}}{S}_{\text{d}}}{N_{\text{d}}}\\ 0\end{array}\right),\mathcal{V}=\left(\begin{array}{c}{\omega }^{\theta }{E}_{\text{h}}+{\tau }_{\text{h}}^{\theta }{E}_{\text{h}}+{\varepsilon }_{\text{h}}^{\theta }{E}_{\text{h}}\\ -{\tau }_{\text{h}}^{\theta }{E}_{\text{h}}+\left({\varepsilon }_{\text{h}}^{\theta }+{\alpha }_{\text{h}}^{\theta }\right)I_{\text{h}}\\ {\tau }_{\text{d}}^{\theta }{E}_{\text{d}}+{\varepsilon }_{\text{d}}^{\theta }{E}_{\text{d}}\\ -{\tau }_{\text{d}}^{\theta }{E}_{\text{d}}+\left({\varepsilon }_{\text{d}}^{\theta }+{\alpha }_{\text{d}}^{\theta }+{\mu }^{\theta }\right)I_{\text{d}}\end{array}\right).\end{array}$ Evaluated at DFE, the Jacobian matrices become: $\begin{array}{}F=\left(\begin{array}{cccc}0& 0& 0& \frac{\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }{S}_{\text{h}}^0}{N_{\text{d}}}\\ 0& 0& 0& 0\\ 0& 0& 0& \frac{{\beta }_{\text{d}}^{\theta }{S}_{\text{d}}^0}{N_{\text{d}}}\\ 0& 0& 0& 0\end{array}\right),V=\left(\begin{array}{cccc}{\varepsilon }_{\text{h}}^{\theta }+{\tau }_{\text{h}}^{\theta }+{\omega }^{\theta }& 0& 0& 0\\ -{\tau }_{\text{h}}^{\theta }& {\alpha }_{\text{h}}^{\theta }+{\varepsilon }_{\text{h}}^{\theta }& 0& 0\\ 0& 0& {\varepsilon }_{\text{d}}^{\theta }+{\tau }_{\text{d}}^{\theta }& 0\\ 0& 0& -{\tau }_{\text{d}}^{\theta }& {\alpha }_{\text{d}}^{\theta }+{\varepsilon }_{\text{d}}^{\theta }+{\mu }^{\theta }\end{array}\right).\end{array}$ Then, the matrix expression $F{\mathcal{R}}_e^{-1}-V$ can be written as follows: 24 $\begin{array}{}F{\mathcal{R}}_e^{-1}-V=\left(\begin{array}{cccc}-\left({\varepsilon }_{\text{h}}^{\theta }+{\tau }_{\text{h}}^{\theta }+{\omega }^{\theta }\right)& 0& 0& \frac{\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }{S}_{\text{h}}^0}{N_{\text{d}}}{\mathcal{R}}_e^{-1}\\ {\tau }_{\text{h}}^{\theta }& -\left({\alpha }_{\text{h}}^{\theta }+{\varepsilon }_{\text{h}}^{\theta }\right)& 0& 0\\ 0& 0& -\left({\varepsilon }_{\text{d}}^{\theta }+{\tau }_{\text{d}}^{\theta }\right)& \frac{{\beta }_{\text{d}}^{\theta }{S}_{\text{d}}^0}{N_{\text{d}}}{\mathcal{R}}_e^{-1}\\ 0& 0& {\tau }_{\text{d}}^{\theta }& -\left({\alpha }_{\text{d}}^{\theta }+{\varepsilon }_{\text{d}}^{\theta }+{\mu }^{\theta }\right)\end{array}\right).\end{array}$

4.3.1. GT Derivation

The matrix $F{\mathcal{R}}_e^{-1}-V$ admits the digraph representation shown in Figure 2, where nodes correspond to infected compartments and edges represent transition pathways. Through sequential node elimination [49, 50], we reduce the system into single node as shown in Figure 3.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

The effective reproduction number emerges as the solution to: $\begin{array}{}-\left({\varepsilon }_{\text{d}}^{\theta }+{\tau }_{\text{d}}^{\theta }\right)+\frac{{\pi }_{\text{d}}^{\theta }{\beta }_{\text{d}}^{\theta }{\tau }_{\text{d}}^{\theta }{\mathcal{R}}_e^{-1}}{N_{\text{d}}\left({\gamma }^{\theta }+{\varepsilon }_{\text{d}}^{\theta }\right)\left({\alpha }_{\text{d}}^{\theta }+{\varepsilon }_{\text{d}}^{\theta }+{\mu }^{\theta }\right)}=0,⟹{\mathcal{R}}_e=\frac{{\pi }_{\text{d}}^{\theta }{\beta }_{\text{d}}^{\theta }{\tau }_{\text{d}}^{\theta }}{N_{\text{d}}\left({\gamma }^{\theta }+{\varepsilon }_{\text{d}}^{\theta }\right)\left({\varepsilon }_{\text{d}}^{\theta }+{\tau }_{\text{d}}^{\theta }\right)\left({\alpha }_{\text{d}}^{\theta }+{\varepsilon }_{\text{d}}^{\theta }+{\mu }^{\theta }\right)}.\end{array}$

4.4. Local Stability Analysis of the RFE

To examine the local stability properties of the RFE, we employ eigenvalue analysis following the approach of Edward [44]. The stability criterion is established through the following theorem.

$$

Statement

Theorem 1 Local stability condition for fractional-order systems.

Consider the initial value problem (IVP) for a fractional-order dynamical system: $\begin{array}{}{}_0^C{\mathcal{D}}_t\mathbf{x}\left(t\right)=\mathbf{F}\left(t,\mathbf{x}\left(t\right)\right),t\ge 0,\mathbf{x}\left(0\right)={\mathbf{x}}_0\in {\mathbb{R}}_{+}^7,\end{array}$ where α^θ ∈ (0, 1] denotes the fractional order, ${}_0^C{\mathcal{D}}_t$ is the Caputo fractional derivative operator, and $\mathbf{F}:{\mathbb{R}}^{+}\times {\mathbb{R}}_{+}^7\to {\mathbb{R}}^7$ defines the system’s vector field. The RFE is locally asymptotically stable if all eigenvalues χ_i (i = 1, …, 7) of the Jacobian matrix evaluated at the RFE satisfy: $\begin{array}{}\left|\text{arg}\left({\chi }_i\right)\right|>\frac{\pi \theta }{2}.\end{array}$

$$

Statement

Proof.

The Jacobian matrix ${\mathcal{J}}_{\xi }$ is derived by computing the partial derivatives of the model system (Equation (10)) with respect to each state variable (S_h, E_h, I_h, S_d, V_d, E_d, I_d) at the RFE ξ: $\begin{array}{}{\mathcal{J}}_{\xi }=\left(\begin{array}{ccccccc}-{\varepsilon }_{\text{h}}^{\theta }& {\omega }^{\theta }& 0& 0& 0& 0& -\frac{\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }{S}_{\text{h}}^0}{N_{\text{d}}}\\ 0& -{\varepsilon }_{\text{h}}^{\theta }-{\tau }_{\text{h}}^{\theta }-{\omega }^{\theta }& 0& 0& 0& 0& \frac{\left(1-\sigma \right){\beta }_{\text{h}}^{\theta }{S}_{\text{h}}^0}{N_{\text{d}}}\\ 0& {\tau }_{\text{h}}^{\theta }& -{\alpha }_{\text{h}}^{\theta }-{\varepsilon }_{\text{h}}^{\theta }& 0& 0& 0& 0\\ 0& 0& 0& -{\gamma }^{\theta }-{\varepsilon }_{\text{d}}^{\theta }& 0& 0& -\frac{{\beta }_{\text{d}}^{\theta }{S}_{\text{d}}^0}{N_{\text{d}}}\\ 0& 0& 0& {\gamma }^{\theta }& -{\varepsilon }_{\text{d}}^{\theta }& 0& 0\\ 0& 0& 0& 0& 0& -{\varepsilon }_{\text{d}}^{\theta }-{\tau }_{\text{d}}^{\theta }& \frac{{\beta }_{\text{d}}^{\theta }{S}_{\text{d}}^0}{N_{\text{d}}}\\ 0& 0& 0& 0& 0& {\tau }_{\text{d}}^{\theta }& -{\alpha }_{\text{d}}^{\theta }-{\varepsilon }_{\text{d}}^{\theta }-{\mu }^{\theta }\end{array}\right).\end{array}$ The characteristic equation $\text{det}\left({\mathcal{J}}_{\xi }-\lambda I\right)=0$ yields the following explicit eigenvalues: $\begin{array}{}{\lambda }_1=-{\varepsilon }_{\text{d}}^{\theta },{\lambda }_2=-{\gamma }^{\theta }-{\varepsilon }_{\text{d}}^{\theta },{\lambda }_3=-{\varepsilon }_{\text{h}}^{\theta },{\lambda }_4=-{\alpha }_{\text{h}}^{\theta }-{\varepsilon }_{\text{h}}^{\theta },\text{\hspace{0.17em}and\hspace{0.17em}}{\lambda }_5=-{\varepsilon }_{\text{h}}^{\theta }-{\tau }_{\text{h}}^{\theta }-{\omega }^{\theta }.\end{array}$

The remaining eigenvalues originate from the reduced matrix P governing the dog population dynamics: $\begin{array}{}\mathcal{P}=\left(\begin{array}{cc}-{\varepsilon }_{\text{d}}^{\theta }-{\tau }_{\text{d}}^{\theta }& \frac{{\beta }_{\text{d}}^{\theta }{S}_{\text{d}}^0}{N_{\text{d}}}\\ {\tau }_{\text{d}}^{\theta }& -{\alpha }_{\text{d}}^{\theta }-{\varepsilon }_{\text{d}}^{\theta }-{\mu }^{\theta }\end{array}\right).\end{array}$

Following Ndendya et al. [1], the stability of $\mathcal{P}$ requires:

• (i)

  Negative trace: $\text{Tr}\left(\mathcal{P}\right)=-\left({\alpha }_{\text{d}}^{\theta }+2{\varepsilon }_{\text{d}}^{\theta }+{\tau }_{\text{d}}^{\theta }+{\mu }^{\theta }\right)<0$.

• (ii)

  Positive determinant: $\text{det}\left(\mathcal{P}\right)=N_{\text{d}}\left({\gamma }^{\theta }+{\varepsilon }_{\text{d}}^{\theta }\right)\left({\varepsilon }_{\text{d}}^{\theta }+{\tau }_{\text{d}}^{\theta }\right)\left({\alpha }_{\text{d}}^{\theta }+{\varepsilon }_{\text{d}}^{\theta }+{\mu }^{\theta }\right)\left(1-{\mathcal{R}}_e\right)>0$.

Condition (i) is always satisfied, while Condition (ii) holds precisely when ${\mathcal{R}}_e<1$. Consequently:

• •

  All eigenvalues satisfy $\left|\text{arg}\left({\chi }_i\right)\right|>\frac{\pi \theta }{2}$ when ${\mathcal{R}}_e<1$.

• •

  The RFE ξ is locally asymptotically stable if ${\mathcal{R}}_e<1$.

• •

  The equilibrium becomes unstable when ${\mathcal{R}}_e\ge 1$.

4.5. Hyers–Ulam Stability of the Generalized Fractional Rabies Model

This section establishes the Hyers–Ulam stability of the generalized fractional-order rabies Model (12), where the operator ${}_0^i{\mathcal{D}}_t^{\theta }$ denotes either the Caputo, CF, or ABC fractional derivative depending on the kernel type. The compact form of the system can be written as: 25 $\begin{array}{}{}_0^i{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)=H\left(t,\mathcal{X}\left(t\right)\right),t\ge 0,\mathcal{X}\left(0\right)={\mathcal{X}}_0\in {\mathbb{R}}_{+},\end{array}$ where $H:{\mathbb{R}}_{+}\times {\mathbb{R}}^n\to {\mathbb{R}}^n$ is a continuous vector-valued function defined by the right-hand side of system (Equation (12)).

Consider the perturbation inequality: 26 $\begin{array}{}{‖{}_0^i{\mathcal{D}}_t^{\theta }\mathcal{X}\left(t\right)-H\left(t,\mathcal{X}\left(t\right)\right)‖}_1<\varepsilon ,\forall t\in \left[0,T\right].\end{array}$

$$

Statement

Definition 10.

A function $z\in \mathbb{E}$ is said to be an approximate solution of Equation (25), if there exists a perturbation $\zeta \left(t\right)\in \mathbb{E}$ such that: 27 $\begin{array}{}\left|\zeta \left(t\right)\right|\le \varepsilon ,\forall t\in \left[0,T\right],\end{array}$ and 28 $\begin{array}{}{}_0^i{\mathcal{D}}_t^{\theta }z\left(t\right)=H\left(t,z\left(t\right)\right)+\zeta \left(t\right),\forall t\in \left[0,T\right].\end{array}$ Applying the appropriate fractional integral corresponding to ${}_0^i{\mathcal{D}}_t^{\theta }$ to both sides of Equation (28), we obtain: 29 $\begin{array}{}z\left(t\right)=z_0+{\mathcal{I}}_t^{\theta }H\left(t,z\left(t\right)\right)+{\mathcal{I}}_t^{\theta }\zeta \left(t\right),\end{array}$ where ${\mathcal{I}}_t^{\theta }$ denotes the corresponding fractional integral operator (e.g., Riemann–Liouville, CF, or Atangana–Baleanu integral) depending on the choice of kernel.

$$

Statement

Definition 11.

Consider the deviation from the exact mild solution: 30 $\begin{array}{}{‖z\left(t\right)-z_0-{\mathcal{I}}_t^{\theta }H\left(t,z\left(t\right)\right)‖}_1.\end{array}$ Substituting Equation (29) into Equation (30) and applying Equation (27) yields: 31 $\begin{array}{}\begin{array}{cc}{‖z\left(t\right)-z_0-{\mathcal{I}}_t^{\theta }H\left(t,z\left(t\right)\right)‖}_1& ={‖{\mathcal{I}}_t^{\theta }\zeta \left(t\right)‖}_1\le \varepsilon \cdot {\mathcal{K}}_{\theta }\left(t\right),\end{array}\end{array}$ where ${\mathcal{K}}_{\theta }\left(t\right)$ is a constant depending on the kernel and order θ. For instance, $\begin{array}{}{\mathcal{K}}_{\theta }\left(t\right)=\left\{\begin{array}{cc}\frac{t^{\theta }}{\Gamma \left(\theta +1\right)}& \text{(Caputo)}\\ \frac{1-\exp \left(-\frac{\theta t}{1-\theta }\right)}{\theta }& \text{(Caputo–Fabrizio)}\\ \frac{1-E_{\theta }\left(-\theta t^{\theta }\right)}{\theta }& \text{(ABC)}\end{array}\right.,\end{array}$ and E_θ(⋅) is the Mittag–Leffler function.

$$

Statement

Definition 12.

Let z(t) and $\mathcal{X}\left(t\right)$ be the perturbed and exact solutions, respectively. If there exists a constant M > 0 such that for every ɛ > 0: 32 $\begin{array}{}{‖z\left(t\right)-\mathcal{X}\left(t\right)‖}_1\le M\varepsilon ,\forall t\in \left[0,T\right],\end{array}$ then the Model (25) is said to be Hyers––Ulam stable.

$$

Statement

Theorem 2.

Assume $H\left(t,\mathcal{X}\right)$ satisfies a Lipschitz condition with constant κ > 0, that is, $\begin{array}{}{‖H\left(t,{\mathcal{X}}_1\right)-H\left(t,{\mathcal{X}}_2\right)‖}_1\le \kappa {‖{\mathcal{X}}_1-{\mathcal{X}}_2‖}_1.\end{array}$ Then, the generalized fractional-order Model (25) is Hyers–Ulam stable on [0, T], if the condition: 33 $\begin{array}{}{\mathcal{K}}_{\theta }\left(T\right)\kappa <1,\end{array}$ is satisfied.

$$

Statement

Proof.

Let $z_0={\mathcal{X}}_0$. Then, from the perturbed solution (Equation (29)) and the mild form of the exact solution: $\begin{array}{}\mathcal{X}\left(t\right)={\mathcal{X}}_0+{\mathcal{I}}_t^{\theta }H\left(t,\mathcal{X}\left(t\right)\right),\end{array}$ the difference satisfies: $\begin{array}{}{‖z\left(t\right)-\mathcal{X}\left(t\right)‖}_1\le {‖z\left(t\right)-z_0-{\mathcal{I}}_t^{\theta }H\left(t,z\left(t\right)\right)‖}_1+{‖{\mathcal{I}}_t^{\theta }\left[H\left(t,z\left(t\right)\right)-H\left(t,\mathcal{X}\left(t\right)\right)\right]‖}_1\le \varepsilon {\mathcal{K}}_{\theta }\left(T\right)+\kappa {\mathcal{K}}_{\theta }\left(T\right){‖z\left(t\right)-\mathcal{X}\left(t\right)‖}_1.\end{array}$ Rearranging terms, we get: $\begin{array}{}{‖z\left(t\right)-\mathcal{X}\left(t\right)‖}_1\le \frac{\varepsilon {\mathcal{K}}_{\theta }\left(T\right)}{1-\kappa {\mathcal{K}}_{\theta }\left(T\right)}=M\varepsilon ,\end{array}$ provided that $\kappa {\mathcal{K}}_{\theta }\left(T\right)<1$. Hence, the model is Hyers–Ulam stable.

5. Parameter Estimation and Numerical Simulations

In this study, the MCMC method was employed to estimate key epidemiological parameters of the rabies model using available field data. MCMC is a powerful Bayesian statistical approach that generates samples from the posterior distribution of parameters by constructing a Markov chain whose equilibrium distribution is the desired posterior. This approach is particularly advantageous in epidemiological modeling where parameters are often uncertain or difficult to measure directly. By using MCMC, we are able to incorporate prior knowledge, account for uncertainty, and obtain robust parameter estimates with credible intervals. This allows for a more reliable calibration of the model compared to classical point-estimation techniques. The adoption of MCMC in this work ensures that the stochastic variability and uncertainty inherent in rabies transmission dynamics are adequately captured, thereby improving the model’s predictive power and aligning projections more closely with observed data.

5.1. Parameter Estimation

Parameter estimation was carried out using rabies incidence data reported in Ulanga District, Morogoro Region, Tanzania, from 2006 to 2009 (Figure 4).

[IMAGE OMITTED. SEE PDF]

To capture long-term transmission patterns and improve model realism, synthetic data were added to extend the dataset over a 20-year period. We employed the MCMC method, a Bayesian framework that enables sampling from the posterior distributions of the parameters conditional on observed data. The estimated parameters were then integrated into the fractional-order rabies transmission model to perform numerical simulations using an appropriate scheme for fractional differential equations. Table 1 summarizes the key model parameters, including their biological interpretations, mean values estimated via MCMC, and relevant literature sources used to inform the model structure.

Figures 5 and 6 present the model fitting results for the exposed and infected human and dog populations over a 20-year period, along with associated uncertainty quantification. Each subplot shows the model-predicted trajectory (solid blue line) overlaid with simulated observations (red circles) and shaded confidence intervals (C.I.) at 50%, 90%, 95%, and 99% levels, respectively. The exposed human population (Figure 5a demonstrates a gradual exponential decline, indicating effective reduction in early-stage infections, with the majority of the data points falling within the 90%–99% confidence bands, suggesting good model reliability. In contrast, the infected human population (Figure 5b exhibits a transient rise followed by a decline, capturing the delayed progression from exposure to active infection. Similarly, the exposed and infected dog populations (Figure 6a,b follow downward trends, reflecting potential impacts of control strategies such as vaccination or culling. Overall, the close alignment between the model fit and observed data across all compartments supports the robustness of the model in capturing rabies transmission dynamics and host responses overtime.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Figure 7 displays trace plots of MCMC samples for eight estimated parameters: β_h, β_d, τ_h, τ_d, γ, α_h, α_d, and μ. Each plot illustrates the sampled values of a single parameter against the MCMC iteration number, ranging up to 2 × 10⁴. These plots serve as crucial diagnostics for assessing MCMC chain convergence and mixing, indicating that the samples are fluctuating randomly around a stable mean without any discernible trends or long-range dependence, suggesting thorough exploration of the posterior distribution. Visually, all trace plots in this figure exhibit these desirable characteristics, implying that the MCMC simulation successfully converged and generated reliable samples from the posterior distributions of these epidemiological transmission and progression parameters. Figure 8 presents a pairwise scatter plot showing the posterior correlations among various estimated transmission and progression parameters from a MCMC analysis. Each blue dot in the individual scatter plots represents a sample from the joint posterior distribution of the two parameters being compared on the respective axes. The diagonal elements display the names of the parameters (β_h, β_d, τ_h, τ_d, γ, α_h, α_d, and μ), while the off-diagonal plots in the lower triangular part of the matrix show the relationships between different parameter pairs. The spread and orientation of the point clouds within each scatter plot reveal the nature and strength of the correlation: for instance, an upward-sloping cloud indicates a positive correlation (e.g., α_h vs. γ), a downward-sloping cloud suggests a negative correlation and a widely dispersed and nontrending cloud indicates a weak or no linear correlation. This visualization is crucial for understanding dependencies and identifiability issues among the estimated parameters.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

5.2. Numerical Simulation

To investigate the dynamic behavior of the proposed human–dog rabies model under fractional-order formulations, we conducted a series of numerical simulations using the Caputo, CF, and ABC operators. These simulations aim to illustrate how memory effects, introduced via fractional derivatives, influence the trajectory of rabies transmission across human and dog populations. The numerical schemes were implemented for a simulation horizon of 20 years to observe long-term epidemic outcomes.

5.3. Numerical Schemes for Rabies Fractional-Order Model

Section 5.3 presents the numerical solution of the rabies model under different derivative operators: the classical (integer-order), Caputo, CF, and Atangana–Baleanu in Caputo sense. Each scheme reflects different memory properties of the system and is suited for various biological interpretations. Let Δt denote the time step, α ∈ (0, 1) the fractional order, and t_a = aΔt the current time level.

5.3.1. Caputo Fractional-Derivative Scheme

The Caputo fractional-derivative captures long-term memory effects via a power–law kernel. It is defined as [51]: $\begin{array}{}{}^C{D}_t^{\alpha }y\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\int }_0^t\frac{dy\left(\tau \right)}{d\tau }\frac{d\tau }{{\left(t-\tau \right)}^{\alpha }},\end{array}$ with 0 < α < 1. Discretizing using the L1 method yields: 34 $\begin{array}{}\begin{array}{cc}{S}_{\text{h}}^{a+1}& =S_{\text{h}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{\mu =1}^a{w}_{a-\mu +1}\cdot \left[{\pi }_{\text{h}}+\omega E_{\text{h}}^{\mu }-\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}^{\mu }}{N_{\text{d}}^{\mu }}{S}_{\text{h}}^{\mu }-{\varepsilon }_{\text{h}}{S}_{\text{h}}^{\mu }\right],\\ E_{\text{h}}^{a+1}& =E_{\text{h}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{\mu =1}^a{w}_{a-\mu +1}\cdot \left[\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}^{\mu }}{N_{\text{d}}^{\mu }}{S}_{\text{h}}^{\mu }-\omega E_{\text{h}}^{\mu }-{\tau }_{\text{h}}{E}_{\text{h}}^{\mu }-{\varepsilon }_{\text{h}}{E}_{\text{h}}^{\mu }\right],\\ I_{\text{h}}^{a+1}& =I_{\text{h}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{\mu =1}^a{w}_{a-\mu +1}\cdot \left[{\tau }_{\text{h}}{E}_{\text{h}}^{\mu }-\left({\varepsilon }_{\text{h}}+{\alpha }_{\text{h}}\right)I_{\text{h}}^{\mu }\right],\\ S_{\text{d}}^{a+1}& =S_{\text{d}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{\mu =1}^a{w}_{a-\mu +1}\cdot \left[{\pi }_{\text{d}}-{\beta }_{\text{d}}\frac{I_{\text{d}}^{\mu }}{N_{\text{d}}^{\mu }}{S}_{\text{d}}^{\mu }-\gamma S_{\text{d}}^{\mu }-{\varepsilon }_{\text{d}}{S}_{\text{d}}^{\mu }\right],\\ E_{\text{d}}^{a+1}& =E_{\text{d}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{\mu =1}^a{w}_{a-\mu +1}\cdot \left[{\beta }_{\text{d}}\frac{I_{\text{d}}^{\mu }}{N_{\text{d}}^{\mu }}{S}_{\text{d}}^{\mu }-{\tau }_{\text{d}}{E}_{\text{d}}^{\mu }-{\varepsilon }_{\text{d}}{E}_{\text{d}}^{\mu }\right],\\ I_{\text{d}}^{a+1}& =I_{\text{d}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{\mu =1}^a{w}_{a-\mu +1}\cdot \left[{\tau }_{\text{d}}{E}_{\text{d}}^{\mu }-\left({\varepsilon }_{\text{d}}+{\alpha }_{\text{d}}+\mu \right)I_{\text{d}}^{\mu }\right],\\ V_{\text{d}}^{a+1}& =V_{\text{d}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{\mu =1}^a{w}_{a-\mu +1}\cdot \left[\gamma S_{\text{d}}^{\mu }-{\varepsilon }_{\text{d}}{V}_{\text{d}}^{\mu }\right],\end{array}\end{array}$ where the weights are: $\begin{array}{}{w}_k=k^{1-\alpha }-{\left(k-1\right)}^{1-\alpha }.\end{array}$

5.3.2. CF Fractional Derivative Scheme

The CF fractional derivative avoids singular kernels using exponential decay memory [52]: $\begin{array}{}{}^{CF}{D}_t^{\alpha }y\left(t\right)=\frac{M\left(\alpha \right)}{1-\alpha }{\int }_0^t{e}^{-\frac{\alpha }{1-\alpha }\left(t-\tau \right)}\frac{dy\left(\tau \right)}{d\tau }d\tau ,\end{array}$ with M(α) such that M(0) = M(1) = 1. The time-discrete form (forward Euler-like) becomes: 35 $\begin{array}{}\begin{array}{cc}{S}_{\text{h}}^{a+1}& =S_{\text{h}}^a+\frac{\Delta t}{1-\frac{\alpha }{2}}\left[{\pi }_{\text{h}}+\omega E_{\text{h}}^a-\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}^a}{N_{\text{d}}^a}{S}_{\text{h}}^a-{\varepsilon }_{\text{h}}{S}_{\text{h}}^a\right],\\ E_{\text{h}}^{a+1}& =E_{\text{h}}^a+\frac{\Delta t}{1-\frac{\alpha }{2}}\left[\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}^a}{N_{\text{d}}^a}{S}_{\text{h}}^a-\omega E_{\text{h}}^a-{\tau }_{\text{h}}{E}_{\text{h}}^a-{\varepsilon }_{\text{h}}{E}_{\text{h}}^a\right],\\ I_{\text{h}}^{a+1}& =I_{\text{h}}^a+\frac{\Delta t}{1-\frac{\alpha }{2}}\left[{\tau }_{\text{h}}{E}_{\text{h}}^a-\left({\varepsilon }_{\text{h}}+{\alpha }_{\text{h}}\right)I_{\text{h}}^a\right],\\ S_{\text{d}}^{a+1}& =S_{\text{d}}^a+\frac{\Delta t}{1-\frac{\alpha }{2}}\left[{\pi }_{\text{d}}-{\beta }_{\text{d}}\frac{I_{\text{d}}^a}{N_{\text{d}}^a}{S}_{\text{d}}^a-\gamma S_{\text{d}}^a-{\varepsilon }_{\text{d}}{S}_{\text{d}}^a\right],\\ E_{\text{d}}^{a+1}& =E_{\text{d}}^a+\frac{\Delta t}{1-\frac{\alpha }{2}}\left[{\beta }_{\text{d}}\frac{I_{\text{d}}^a}{N_{\text{d}}^a}{S}_{\text{d}}^a-{\tau }_{\text{d}}{E}_{\text{d}}^a-{\varepsilon }_{\text{d}}{E}_{\text{d}}^a\right],\\ I_{\text{d}}^{a+1}& =I_{\text{d}}^a+\frac{\Delta t}{1-\frac{\alpha }{2}}\left[{\tau }_{\text{d}}{E}_{\text{d}}^a-\left({\varepsilon }_{\text{d}}+{\alpha }_{\text{d}}+\mu \right)I_{\text{d}}^a\right],\\ V_{\text{d}}^{a+1}& =V_{\text{d}}^a+\frac{\Delta t}{1-\frac{\alpha }{2}}\left[\gamma S_{\text{d}}^a-{\varepsilon }_{\text{d}}{V}_{\text{d}}^a\right].\end{array}\end{array}$

5.3.3. Atangana–Baleanu in Caputo Sense Scheme

The ABC derivative uses the nonlocal and nonsingular Mittag–Leffler kernel. It is defined as [46]: $\begin{array}{}{}^{ABC}{D}_t^{\alpha }y\left(t\right)=\frac{B\left(\alpha \right)}{1-\alpha }{\int }_0^t{E}_{\alpha }\left(-\frac{\alpha }{1-\alpha }{\left(t-\tau \right)}^{\alpha }\right)\frac{dy\left(\tau \right)}{d\tau }d\tau ,\end{array}$ where $\begin{array}{}B\left(\alpha \right)=1-\alpha +\frac{\alpha }{\Gamma \left(\alpha \right)}.\end{array}$

The discrete scheme becomes: 36 $\begin{array}{}\begin{array}{cc}{S}_{\text{h}}^{a+1}& =S_{\text{h}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{B\left(\alpha \right)}\sum _{\mu =1}^a\left[{\pi }_{\text{h}}+\omega E_{\text{h}}^{\mu }-\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}^{\mu }}{N_{\text{d}}^{\mu }}{S}_{\text{h}}^{\mu }-{\varepsilon }_{\text{h}}{S}_{\text{h}}^{\mu }\right],\\ E_{\text{h}}^{a+1}& =E_{\text{h}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{B\left(\alpha \right)}\sum _{\mu =1}^a\left[\left(1-\sigma \right){\beta }_{\text{h}}\frac{I_{\text{d}}^{\mu }}{N_{\text{d}}^{\mu }}{S}_{\text{h}}^{\mu }-\omega E_{\text{h}}^{\mu }-{\tau }_{\text{h}}{E}_{\text{h}}^{\mu }-{\varepsilon }_{\text{h}}{E}_{\text{h}}^{\mu }\right],\\ I_{\text{h}}^{a+1}& =I_{\text{h}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{B\left(\alpha \right)}\sum _{\mu =1}^a\left[{\tau }_{\text{h}}{E}_{\text{h}}^{\mu }-\left({\varepsilon }_{\text{h}}+{\alpha }_{\text{h}}\right)I_{\text{h}}^{\mu }\right],\\ S_{\text{d}}^{a+1}& =S_{\text{d}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{B\left(\alpha \right)}\sum _{\mu =1}^a\left[{\pi }_{\text{d}}-{\beta }_{\text{d}}\frac{I_{\text{d}}^{\mu }}{N_{\text{d}}^{\mu }}{S}_{\text{d}}^{\mu }-\gamma S_{\text{d}}^{\mu }-{\varepsilon }_{\text{d}}{S}_{\text{d}}^{\mu }\right],\\ E_{\text{d}}^{a+1}& =E_{\text{d}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{B\left(\alpha \right)}\sum _{\mu =1}^a\left[{\beta }_{\text{d}}\frac{I_{\text{d}}^{\mu }}{N_{\text{d}}^{\mu }}{S}_{\text{d}}^{\mu }-{\tau }_{\text{d}}{E}_{\text{d}}^{\mu }-{\varepsilon }_{\text{d}}{E}_{\text{d}}^{\mu }\right],\\ I_{\text{d}}^{a+1}& =I_{\text{d}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{B\left(\alpha \right)}\sum _{\mu =1}^a\left[{\tau }_{\text{d}}{E}_{\text{d}}^{\mu }-\left({\varepsilon }_{\text{d}}+{\alpha }_{\text{d}}+\mu \right)I_{\text{d}}^{\mu }\right],\\ V_{\text{d}}^{a+1}& =V_{\text{d}}^0+\frac{{\left(\Delta t\right)}^{\alpha }}{B\left(\alpha \right)}\sum _{\mu =1}^a\left[\gamma S_{\text{d}}^{\mu }-{\varepsilon }_{\text{d}}{V}_{\text{d}}^{\mu }\right].\end{array}\end{array}$

5.4. Comparative Simulation Results

The comparative simulation for fractional and classical models was conducted for 20 years. The analysis of human and dog population dynamics under different fractional operators (Caputo, CF, and ABC) with θ = 0.75 reveals significant differences in disease progression patterns compared to the classical integer-order model (θ = 1). Figure 9 illustrates the human population dynamics (susceptible S_h, exposed E_h, and infected I_h). The classical model (integer-order, θ = 1) exhibits the fastest growth and decline in infectious populations. Among the fractional models, ABC (θ = 0.75) exhibits the slowest disease evolution, indicating a stronger memory effect and delayed transmission peak.

[IMAGE OMITTED. SEE PDF]

Figure 10 presents the dog population dynamics, where the fractional models show delayed peaks in exposed and infectious compartments compared to the integer-order case. The ABC derivative exhibits smoother transitions between compartments, reflecting its Mittag–Leffler kernel’s ability to capture both short- and long-term memory effects. These differences highlight how fractional operators can provide more nuanced representations of disease transmission, particularly in accounting for latency periods and intervention delays.

[IMAGE OMITTED. SEE PDF]

5.5. Impact of Varying Intervention Parameter Rates on the Rabies Transmission Dynamics Using the ABC Fractional-Derivative With Order θ = 0.75

The analysis of intervention strategies using the ABC fractional derivative (θ = 0.75) provides critical insights for rabies control. Figure 11 demonstrates the effect of PEP rate ω on human populations. Higher PEP rates (ω = 0.5) significantly reduce both exposed and infectious human cases compared to lower rates (ω = 0.1), emphasizing the importance of rapid treatment after exposure. Figure 12 examines vaccination rates (γ) in dog populations. Increasing vaccination from γ = 0.1 to γ = 0.9 substantially decreases exposed and infectious dog cases, with the higher rate achieving near-elimination of dog rabies within 10 years. This nonlinear response suggests threshold effects in vaccination coverage. Figure 13 shows as the culling rate μ increases, both E_d and I_d decrease sharply. At μ = 0.9, the infected dog population nearly vanishes within the first 5 years.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

6. Discussion

This study presents a comprehensive exploration of rabies transmission dynamics using fractional-order models involving Caputo, CF, and ABC derivatives. The comparative simulations under a fixed fractional order θ = 0.75 revealed that the model embedded with the ABC derivative produced the most restrained epidemic progression across all compartments, particularly in the exposed and infectious classes of both humans and dogs.

The observed suppression of peak infection levels in the ABC model is attributable to its incorporation of a nonsingular and nonlocal Mittag–Leffler kernel, which imparts a smoother memory effect compared to the power–law memory in the Caputo model. This property enables the ABC model to better replicate the delayed response and long-term interactions in disease dynamics, as often seen in real-world zoonotic transmission systems such as rabies [46]. The CF derivative, although also nonsingular, offered intermediate dynamics, while the classical (integer-order) model overestimated disease intensity due to its lack of memory representation.

In agreement with our findings, Aydogan et al. [20] demonstrated that applying the CF operator in a rabies framework led to more stable trajectories and reduced peaks, highlighting the benefits of memory-aware operators. Similarly, Zainab et al. [21] used the ABC operator to model rabies in animal and reported enhanced predictive capabilities, particularly in forecasting long-term endemic levels. Garg and Chauhan [23] supported the theoretical robustness of fractional-order frameworks through fixed-point stability analyses, reinforcing their practical and analytical advantages over classical models.

What distinguishes the present study is the integration of real-world parameter calibration using MCMC techniques with a full comparative analysis of three fractional operators within the same epidemiological framework. This unified approach permits a consistent comparison of model performance under identical biological and computational settings, thereby improving the validity of the conclusions. The incorporation of parameter uncertainty further enhances the realism and applicability of the predictions.

Moreover, the intervention-focused simulations using the ABC model provide actionable insights. Varying the PEP rate ω clearly showed that higher values significantly reduce the burden of human exposure and infection. Likewise, increasing the dog vaccination rate γ and the culling rate μ markedly decreased the prevalence of infection among dogs. These findings affirm that combined intervention strategies are far more effective than any single measure, echoing previous studies emphasizing multi-pronged approaches to rabies control [53].

Overall, the results suggest that fractional-order modeling, especially with the ABC operator, not only enhances the theoretical fidelity of the model but also improves its practical utility in designing and evaluating rabies control programs. The memory effect embedded within such operators aligns well with the prolonged incubation, delayed response to vaccination, and lingering effects of transmission observed in rabies dynamics.

7. Conclusion

This study has demonstrated the potential of fractional-order models in accurately capturing the complex and delayed dynamics of rabies transmission between humans and dogs. By integrating Caputo, CF, and ABC derivatives into a unified framework, we observed that fractional derivatives-particularly the ABC operator-can more realistically reflect the memory and hereditary effects of rabies transmission and control interventions. Numerical simulations using MCMC-estimated parameters revealed that the ABC model provides the most conservative estimate of infection spread, resulting in reduced peaks and slower progression in both human and dog populations. These results underscore the importance of including memory effects in epidemiological modeling, especially for diseases with noninstantaneous incubation and prolonged infectious periods. The analysis of intervention strategies further highlighted that effective rabies control requires a multifaceted approach, including high coverage of dog vaccination, rapid PEP for humans, and consistent culling of infected dogs. When implemented in combination, these measures significantly reduced the disease burden in both host populations under the ABC fractional framework.

Future studies should extend this work by considering spatial heterogeneity and stochastic effects, particularly in rural or resource-limited settings where rabies persists. Incorporating time-varying control measures, cost-effectiveness analysis, and real-time data assimilation through data-driven fractional models could further enhance the utility of this modeling framework. Moreover, the exploration of variable-order fractional derivatives and hybrid models may offer even greater flexibility and accuracy in understanding rabies dynamics and informing public health interventions.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest.

Author Contributions

All authors contributed equally to the preparation and publication of this work.

Funding

The authors did not receive any funding for this research.

Acknowledgments

We confirm that no AI software or tools were used at any stage of the preparation of this manuscript. All writing, editing, data analysis, and figure preparation were performed solely by the authors.

References

1 Ndendya J. Z., Leandry L., and Kipingu A. M., A Next-Generation Matrix Approach Using Routh-Hurwitz Criterion and Quadratic Lyapunov Function for Modeling Animal Rabies With Infective Immigrants, Healthcare Analytics. (2023) 4, https://doi.org/10.1016/j.health.2023.100260, 100260.

2 Burgos-Cáceres S., Canine Rabies: A Looming Threat to Public Health, Animals. (2011) 1, no. 4, 326–342, https://doi.org/10.3390/ani1040326, 2-s2.0-84879916425.

3 Kumar A., Bhatt S., Kumar A., and Rana T., Canine Rabies: An Epidemiological Significance, Pathogenesis, Diagnosis, Prevention, and Public Health Issues, Comparative Immunology, Microbiology and Infectious Diseases. (2023) 97, https://doi.org/10.1016/j.cimid.2023.101992, 101992.

4 Abraham S., Ravindran J., Abishaw N., Sandam N. P., Thimmareddy P., and Govindaraju G., Review on Rabies and Vaccines, International Journal of Current Microbiology and Applied Sciences. (2017) 6, no. 12, 2064–2085, https://doi.org/10.20546/ijcmas.2017.612.237.

5 Gu Y., Rabies: The Scientific Basis and Its Public Threat, Highlights in Science, Engineering and Technology. (2023) 45, 91–96, https://doi.org/10.54097/hset.v45i.7316.

6 Ibrahim S., Audu S., Usman A., and Kaltugo B., Rabies in a 6-Week Old Bunaji-Bull Calf in Zaria: A Case Report, Journal of Microbes and Microbiology Techniques. (2017) 1, no. 1, 101.

7 Sun L., Leung G. M., and Cook A. R., FisherD., Modelling to Inform Strategy, Infectious Disease Emergencies: Preparedness and Response, 2025, 1st edition, NUS Press, https://emergencies.pubpub.org/pub/modelling.

8 Mutuguta J. W., Rotich T., and Chepkwony I., Mathematical Modeling of Infectious Disease Transmission Dynamics in a Metapopulation, 2019, Technical report http://repository.mut.ac.ke:8080/xmlui/handle/123456789/4380.

9 Wei Y., Guan J., Zhao Y., and Chen F., YeD. Q., Progress and Perspective of Transmission Dynamics Models in Prevention and Control of Infectious Diseases, Progress in China Epidemiology, 2022, Springer, https://doi.org/10.1007/978-981-19-2199-5_22.

10 Louzoun Y., Solomon S., Atlan H., and Cohen I. R., Modeling Complexity in Biology, Physica A: Statistical Mechanics and Its Applications. (2001) 297, no. 1-2, 242–252, https://doi.org/10.1016/S0378-4371(01)00201-1, 2-s2.0-0035426487.

11 Louzoun Y., The Evolution of Mathematical Immunology, Immunological Reviews. (2007) 216, no. 1, 9–20, https://doi.org/10.1111/j.1600-065X.2006.00495.x, 2-s2.0-33947283687.

12 Abdelsheed I. G., Fractional Calculus: Numerical Methods and SIR Models, 2016, PhD thesis.

13 Skwara U., Mozyrska D., Aguiar M., and Stollenwerk N., Applications of Fractional Calculus to Mathematical Modelling of Vector-Borne Diseases, 2021.

14 Lopes A. M. and Chen L., Fractional Order Systems and Their Applications, Fractal and Fractional. (2022) 6, no. 7, https://doi.org/10.3390/fractalfract6070389, 389.

15 Bas E., Acay B., and Ozarslan R., Fractional Models With Singular and Non-Singular Kernels for Energy Efficient Buildings, Chaos: An Interdisciplinary Journal of Nonlinear Science. (2019) 29, no. 2, https://doi.org/10.1063/1.5082390, 2-s2.0-85061361037.

16 Atangana A., AtanganaA., Fractional Operators and Their Applications, Fractional Operators With Constant and Variable Order With Application to Geo-Hydrology, 2018, Academic Press, 79–112, https://doi.org/10.1016/B978-0-12-809670-3.00005-9.

17 Hattaf K., On Some Properties of the New Generalized Fractional Derivative With Non-Singular Kernel, Mathematical Problems in Engineering. (2021) 2021, no. 1, https://doi.org/10.1155/2021/1580396, 1580396.

18 Ercan A., Comparative Analysis for Fractional Nonlinear Sturm-Liouville Equations With Singular and Non-Singular Kernels, AIMS Mathematics. (2022) 7, no. 7, 13325–13343, https://doi.org/10.3934/math.2022736.

19 Yavuz M., Characterizations of Two Different Fractional Operators Without Singular Kernel, Mathematical Modelling of Natural Phenomena. (2019) 14, no. 3, https://doi.org/10.1051/mmnp/2018070, 2-s2.0-85062022411, 302.

20 Aydogan S. M., Baleanu D., Mohammadi H., and Rezapour S., On the Mathematical Model of Rabies by Using the Fractional Caputo-Fabrizio Derivative, Advances in Difference Equations. (2020) 2020, no. 1, https://doi.org/10.1186/s13662-020-02798-4, 382.

21 Zainab M., Boulaaras S., Aslam A., Shafiq S., Hussain T., and Ozair M., Study of Fractional Order Rabies Transmission Model via Atangana-Baleanu Derivative, Scientific Reports. (2024) 14, no. 1, https://doi.org/10.1038/s41598-024-77282-0, 25875.

22 Zarin R., Ahmed I., Kumam P., Zeb A., and Din A., Fractional Modeling and Optimal Control Analysis of Rabies Virus Under the Convex Incidence Rate, Results in Physics. (2021) 28, https://doi.org/10.1016/j.rinp.2021.104665, 104665.

23 Garg P. and Chauhan S. S., Stability Analysis of a Solution for the Fractional-Order Model on Rabies Transmission Dynamics Using a Fixed-Point Approach, Mathematical Methods in the Applied Sciences. (2024) https://doi.org/10.1002/mma.9903.

24 Alshabanat A., Jleli M., Kumar S., and Samet B., Generalization of Caputo-Fabrizio Fractional Derivative and Applications to Electrical Circuits, Frontiers in Physics. (2020) 8, https://doi.org/10.3389/fphy.2020.00064, 64.

25 Khan M. A., Ullah S., and Kumar S., A Robust Study on 2019-Ncov Outbreaks Through Non-Singular Derivative, The European Physical Journal Plus. (2021) 136, no. 2, https://doi.org/10.1140/epjp/s13360-021-01159-8, 168.

26 Dhar B., Gupta P. K., and Sajid M., Solution of a Dynamical Memory Effect Covid-19 Infection System With Leaky Vaccination Efficacy by Non-Singular Kernel Fractional Derivatives, Mathematical Biosciences and Engineering. (2022) 19, no. 5, 4341–4367, https://doi.org/10.3934/mbe.2022201.

27 Kumar S., Kumar A., Samet B., and Dutta H., A Study on Fractional Host-Parasitoid Population Dynamical Model to Describe Insect Species, Numerical Methods for Partial Differential Equations. (2021) 37, no. 2, 1673–1692, https://doi.org/10.1002/num.22603.

28 Ghanbari B. and Kumar S., A Study on Fractional Predator-Prey-Pathogen Model With Mittag-Leffler Kernel-Based Operators, Numerical Methods for Partial Differential Equations. (2024) 40, no. 1, https://doi.org/10.1002/num.22689, e22689.

29 Kumar S., Chauhan R. P., Momani S., and Hadid S., Numerical Investigations on Covid-19 Model Through Singular and Non-Singular Fractional Operators, Numerical Methods for Partial Differential Equations. (2024) 40, no. 1, https://doi.org/10.1002/num.22707, e22707.

30 Kumar S., Kumar R., Momani S., and Hadid S., A Study on Fractional Covid-19 Disease Model by Using Hermite Wavelets, Mathematical Methods in the Applied Sciences. (2023) 46, no. 7, 7671–7687, https://doi.org/10.1002/mma.7065.

31 Faridi W. A., Wazwaz A.-M., Mostafa A. M., Myrzakulov R., and Umurzakhova Z., The Lie Point Symmetry Criteria and Formation of Exact Analytical Solutions for Kairat-II Equation: Paul-Painlevé Approach, Chaos, Solitons & Fractals. (2024) 182, https://doi.org/10.1016/j.chaos.2024.114745, 114745.

32 Asghar U., Asjad M. I., Faridi W. A., and Akgül A., Novel Solitonic Structure, Hamiltonian Dynamics and Lie Symmetry Algebra of Biofilm, Partial Differential Equations in Applied Mathematics. (2024) 9, https://doi.org/10.1016/j.padiff.2024.100653, 100653.

33 Tipu G. H., Faridi W. A., Alshehri M., and Yao F., The Propagation of Hyper-Geometric Optical Soliton Waves in Plasma Dynamics for the Integrable Zhanbota-Iia Equation, Modern Physics Letters B. (2025) 39, no. 22, https://doi.org/10.1142/S0217984925500769, 2550076.

34 Faridi W. A., Jhangeer A., Riaz M. B., Asjad M. I., and Muhammad T., The Fractional Soliton Solutions of the Dynamical System of Equations for Ion Sound and Langmuir Waves: A Comparative Analysis, Scientific Reports. (2024) 14, no. 1, https://doi.org/10.1038/s41598-024-73983-8, 30473.

35 Faridi W. A., Bakar M. Abu, Jhangeer A., Tawfiq F., Myrzakulov R., and Naizagarayeva A., The Generalized Soliton Wave Structures and Propagation Visualization for Akbota Equation, Zeitschrift für Naturforschung A. (2024) 79, no. 12, 1075–1091, https://doi.org/10.1515/zna-2024-0120.

36 Jhangeer A., Faridi W. A., and Alshehri M., Soliton Wave Profiles and Dynamical Analysis of Fractional Ivancevic Option Pricing Model, Scientific Reports. (2024) 14, no. 1, https://doi.org/10.1038/s41598-024-74770-1, 23804.

37 Faridi W. A., Iqbal M., and Ramzan B., et al.The Formation of Invariant Optical Soliton Structures to Electric-Signal in the Telegraph Lines on Basis of the Tunnel Diode and Chaos Visualization, Conserved Quantities: Lie Point Symmetry Approach, Optik. (2024) 305, https://doi.org/10.1016/j.ijleo.2024.171785, 171785.

38 Faridi W. A. and AlQahtani S. A., The Formation of Invariant Exact Optical Soliton Solutions of Landau-Ginzburg-Higgs Equation via Khater Analytical Approach, International Journal of Theoretical Physics. (2024) 63, no. 2, https://doi.org/10.1007/s10773-024-05559-1, 31.

39 Höhle M. and Jørgensen E., Estimating Parameters for Stochastic Epidemics, 2002, The Royal Veterinary and Agricultural University.

40 Boni M. F., Mølbak K., and Krogfelt K. A., Inferring the Time of Infection from Serological Data, Handbook of Infectious Disease Data Analysis, 2019, Chapman and Hall/CRC, 287–303.

41 Li M., Dushoff J., and Bolker B. M., Fitting Mechanistic Epidemic Models to Data: A Comparison of Simple Markov Chain Monte Carlo Approaches, Statistical Methods in Medical Research. (2018) 27, no. 7, 1956–1967, https://doi.org/10.1177/0962280217747054, 2-s2.0-85047953488.

42 O’Neill P. D., Balding D. J., Becker N. G., Eerola M., and Mollison D., Analyses of Infectious Disease Data From Household Outbreaks by Markov Chain Monte Carlo Methods, Journal of the Royal Statistical Society Series C: Applied Statistics. (2000) 49, no. 4, 517–542, https://doi.org/10.1111/1467-9876.00210, 2-s2.0-0034354606.

43 Podlubny I., Fractional Differential Equations: an Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications, 1999, Elsevier.

44 Edward S., A Fractional Order Model for the Transmission Dynamics of Shigellosis, Heliyon. (2024) 10, no. 10, https://doi.org/10.1016/j.heliyon.2024.e31242, e31242.

45 Caputo M. and Fabrizio M., A New Definition of Fractional Derivative Without Singular Kernel, Progress in Fractional Differentiation & Applications. (2015) 1, no. 2, 73–85.

46 Atangana A. and Baleanu D., New Fractional Derivatives With Nonlocal and Non-Singular Kernel: Theory and Application to Heat Transfer Model, Thermal Science. (2016) 20, no. 2, 763–769, https://doi.org/10.2298/TSCI160111018A, 2-s2.0-84971631213.

47 Ndendya J. Z., Mwasunda J. A., Edward S., and Shaban N., A Fractional-Order Model for Rabies Transmission Dynamics Using the Atangana-Baleanu-Caputo Derivative and Mcmc Methods, Scientific African. (2025) 29, https://doi.org/10.1016/j.sciaf.2025.e02800, e02800.

48 Agarwal R., Airan P., and Agarwal R. P., Exploring the Landscape of Fractional-Order Models in Epidemiology: A Comparative Simulation Study, Axioms. (2024) 13, no. 8, https://doi.org/10.3390/axioms13080545, 545.

49 de Camino-Beck T., Lewis M. A., and van den Driessche P., A Graph-Theoretic Method for the Basic Reproduction Number in Continuous Time Epidemiological Models, Journal of Mathematical Biology. (2009) 59, no. 4, 503–516, https://doi.org/10.1007/s00285-008-0240-9, 2-s2.0-70349941248.

50 Oswald S., Mureithi E., Tsanou B., Chapwanya M., Mashoto K., and Kahesa C., Mcmc-Driven Mathematical Modeling of the Impact of Hpv Vaccine Uptake in Reducing Cervical Cancer, Scientific African. (2025) 28, https://doi.org/10.1016/j.sciaf.2025.e02633, e02633.

51 Diethelm K., Ford N. J., and Freed A. D., A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations, Nonlinear Dynamics. (2002) 29, no. 1–4, 3–22, https://doi.org/10.1023/A:1016592219341, 2-s2.0-0036650479.

52 Caputo M. and Fabrizio M., A New Definition of Fractional Derivative Without Singular Kernel, Progress in Fractional Differentiation and Applications. (2015) 1, no. 2, 1–13.

53 World Health Organization, Zero by 30: The Global Strategic Plan to End Human Deaths From Dog-Mediated Rabies by 2030, 2018, https://www.who.int/publications/i/item/9789241513838.