Next Article in Journal
Data Augmentation for Deep Learning-Based Speech Reconstruction Using FOC-Based Methods
Next Article in Special Issue
An Efficient Petrov–Galerkin Scheme for the Euler–Bernoulli Beam Equation via Second-Kind Chebyshev Polynomials
Previous Article in Journal
Topology of Locally and Non-Locally Generalized Derivatives
Previous Article in Special Issue
Investigating a Nonlinear Fractional Evolution Control Model Using W-Piecewise Hybrid Derivatives: An Application of a Breast Cancer Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Discrete Fractional-Order Modeling of Recurrent Childhood Diseases Using the Caputo Difference Operator

1
Department of Mathematics, College of Science, University of Ha’il, Ha’il 55473, Saudi Arabia
2
School of Science, Monash University, Selangor 47500, Malaysia
3
Department of Mathematics, College of Science, Qassim University, Buraydah 51452, Saudi Arabia
4
Department of Mathematics, Turabah University College, Taif University, Taif 21944, Saudi Arabia
5
Department of Mathematics, Faculty of Science, University of Tabuk, Tabuk 71491, Saudi Arabia
6
Department of Mathematics, Faculty of Science, Islamic University of Madinah, Medina 42351, Saudi Arabia
7
Department of Mathematics, College of Science and Humanities, Prince Sattam Bin Abdulaziz University, Al-Kharj 11942, Saudi Arabia
*
Authors to whom correspondence should be addressed.
Fractal Fract. 2025, 9(1), 55; https://doi.org/10.3390/fractalfract9010055
Submission received: 20 December 2024 / Revised: 11 January 2025 / Accepted: 16 January 2025 / Published: 20 January 2025

Abstract

:
This paper presents a new SIRS model for recurrent childhood diseases under the Caputo fractional difference operator. The existence theory is established using Brouwer’s fixed-point theorem and the Banach contraction principle, providing a comprehensive mathematical foundation for the model. Ulam stability is demonstrated using nonlinear functional analysis. Sensitivity analysis is conducted based on the variation of each parameter, and the basic reproduction number ( R 0 ) is introduced to assess local stability at two equilibrium points. The stability analysis indicates that the disease-free equilibrium point is stable when R 0 < 1 , while the endemic equilibrium point is stable when R 0 > 1 and otherwise unstable. Numerical simulations demonstrate the model’s effectiveness in capturing realistic scenarios, particularly the recurrent patterns observed in some childhood diseases.

1. Introduction

Epidemiology is a prominent research field within the biological sciences that examines health, illness, and related factors at the population level. The term “epidemiology” is derived from three Greek words, epi (upon), demos (public), and logos (study), indicating that this field specifically focuses on human populations [1]. Currently, the study of infectious diseases is a rich area for applying mathematical research in biological sciences. Infectious diseases contribute significantly to morbidity and mortality rates in both human and non-human populations. Consequently, the ability to utilize formal models to reduce these burdens has piqued the interest of both mathematicians and biologists [2]. Many researchers have dedicated their efforts to various established mathematical models to explore the spread of infectious diseases within the realm of mathematical biology [3].
Childhood diseases such as measles, rotavirus, varicella, mumps, respiratory syncytial virus (RSV), norovirus, enterovirus, influenza, smallpox, chickenpox, rubella, and polio are prevalent and pose significant risks as epidemic infectious diseases affecting young individuals. These diseases can spread rapidly among children aged five and below due to frequent contact with their peers in school or other social settings. Vaccination remains an effective control strategy against these childhood diseases [4,5]. Although numerous infectious and non-infectious conditions affect children, it is impractical to address all of them in this study. Instead, we highlight some common viral and bacterial infections, as well as immunological disorders such as ear infections, glue ear, croup, and hand, foot, and mouth disease, which tend to affect children more severely than adults. The distinctive characteristics of childhood diseases make them particularly well suited for mathematical modeling, prompting significant research attention in this area [3]. For instance, a SIR model was developed in [6] to track the time-dependent dynamics of a childhood disease in the presence of a preventive vaccine.
Incorporating reinfection into the SIR model for childhood diseases is crucial for accurately capturing the real-world dynamics of these diseases. Many childhood illnesses do not confer lifelong immunity, leading to a potential for reinfection after a certain period. This phenomenon is observed in diseases such as pertussis (whooping cough) [7], rotavirus [8], influenza [9], respiratory syncytial virus (RSV) [10], and norovirus [11], where immunity gradually wanes over time. Additionally, vaccination schedules that include multiple and booster doses highlight the necessity for renewed immunity. Incorporating reinfection into the model enhances our understanding of outbreaks, assists in planning vaccination strategies, and provides insight into the long-term behavior of diseases in populations with high birth rates.
This paper introduces an enhanced SIR model for childhood diseases by incorporating a reinfection rate, building on the foundational model presented in [12]. The assumptions of the model are as follows:
  • Newborns enter the susceptible class at a constant birth rate Λ , with a fraction λ v vaccinated at birth.
  • Individuals in the population have a constant natural death rate μ , which affects all compartments equally.
  • The disease spreads via contact between susceptible ( S ) and infected ( I ) individuals, with a transmission rate β .
  • Infected individuals recover at rate σ and move to the recovered class ( R ) but lose immunity over time at rate γ , returning to the susceptible class.
  • Individuals die at a natural rate μ , and infected individuals have an additional disease-related mortality rate δ .
With the above assumptions, we consider the following SIR model of ordinary differential equations:
S ( t ) = ( 1 λ v ) Λ β S I Λ S + γ R ,   S ( 0 ) 0 , I ( t ) = β S I ( σ + μ + δ ) I , I ( 0 ) 0 , R ( t ) = λ v Λ + σ I ( μ + γ ) R , R ( 0 ) 0 ,
where Λ is the constant birth rate (new susceptible individuals), μ is the natural death rate, β is the average contact rate leading to infection, σ is the recovery rate from infection, λ v is the proportion of the population vaccinated at birth, γ is the rate at which recovered individuals lose immunity and return to the susceptible class, and δ is the death rate due to disease. The schematic diagram of the proposed model (1) is given in Figure 1.
Fractional calculus (FC) tools are highly emerging in today’s science, as they have been proven to be the best tools for modeling real-world problems due to their memory properties, which are not held by classical calculus tools. The history of FC is as old as that of classical calculus [13]. These tools find widespread applications across diverse scientific disciplines, including, but not limited to, physics and polymer technology [14], robotics [15], electrical circuits [16], bioengineering [17], electrochemistry [18], fluid mechanics [19], electrodynamics of complex media [20], optical fibers [21], control theory [22], thermodynamics [23], viscoelasticity [24], aerodynamics [25], capacitor theory [26], disease modeling [27], blood flow [28], and the fitting of experimental data [29], as well as many other areas [30,31].
Several definitions of fractional-order operators have been proposed in the literature, each with its own benefits and drawbacks [32]. One such operator is known as a fractional-order discrete difference operator, which is a natural extension of the integer-order difference equation. Typically, real-world data exist in discrete form; therefore, these operators perform very well when studying phenomena in the real world under these operators, as they provide output in the discrete form. They have a wide range of applications, including image processing [33], game theory [34], ecology [35], neural networks [36], epidemiology [37], economics [38], and electrical engineering [39], among others. These operators are defined as follows:
Definition 1 
([37]). For function f ( t ) , the ϕ-Caputo fractional difference operator is represented by
C Δ a ϕ f ( t ) = Δ a ( m ϕ ) Δ m f ( t ) = 1 Γ ( m ϕ ) s = a t ( m ϕ ) ( t s 1 ) m ϕ 1 Δ m f ( s ) ,
where t N a + m ϕ , m = ϕ + 1 , and ϕ > 0 with ϕ N . The falling factorial function, denoted by ( t 1 s ) ( m ϕ 1 ) , and the m -th integer difference operator, represented as Δ m f ( s ) , can be described as follows:
( t s 1 ) ( m ϕ 1 ) = Γ ( t s ) Γ ( t s m + ϕ + 1 ) ,
and
Δ m f ( t ) = Δ ( Δ m 1 f ( t ) ) = k = 0 m m k ( 1 ) m k f ( t + k ) , t N a .
Remark 1. 
For m = 1 , the ϕ-Caputo fractional difference operator is defined as
C Δ a ϕ f ( t ) = Δ a ( 1 ϕ ) Δ f ( t ) = 1 Γ ( 1 ϕ ) s = a t ( 1 ϕ ) ( t s 1 ) ϕ Δ f ( s ) ,
where t N a + 1 ϕ .
Definition 2 
([37]). Let f ( t ) be a function defined on N a with the mapping f : N a R . In this context, the ϕ-th fractional sum of f is defined as follows:
Δ a ϕ f ( t ) = 1 Γ ( ϕ ) s = a t ϕ ( t s 1 ) ϕ 1 f ( s ) ,
where ϕ > 0 , t N a + ϕ = { a + ϕ , a + ϕ + 1 , } is a time scale and a R .
Lemma 1 
([40]). Let ϕ > 0 and f be defined on N a , then
Δ ϕ C Δ ϕ f ( t ) = f ( t ) k = 0 m 1 ( t a ) ( k ) k ! Δ k f ( a ) = f ( t ) + c 0 + c 1 t + + c m t m 1 ,
where m denotes the smallest integer satisfying m ϕ , and c i R , i = 1 , 2 , , m 1 .
Lemma 2 
([40]). The following equality holds:
s = 0 t ϕ ( t s 1 ) ( ϕ 1 ) = Γ ( t + 1 ) ϕ Γ ( t ϕ + 1 ) .
This identity simplifies the summation in terms of the Gamma function.
Given the significance of the discrete fractional-order derivative, which also encompasses the ordinary order case, the discrete fractional-order model from Equation (1) can be expressed in the following form:
C Δ ϕ S ( t ) = ( 1 λ v ) Λ β S ( t + ϕ 1 ) I ( t + ϕ 1 ) Λ S ( t + ϕ 1 ) + γ R ( t + ϕ 1 ) , C Δ ϕ I ( t ) = β S ( t + ϕ 1 ) I ( t + ϕ 1 ) ( σ + μ + δ ) I ( t + ϕ 1 ) , C Δ ϕ R ( t ) = λ v Λ + σ I ( t + ϕ 1 ) ( μ + γ ) R ( t + ϕ 1 ) ,
where t [ 0 , b ] N 0 . So the proposed model (9) can be expressed in the following four equivalent equations:
S ( t ) = S ( 0 ) + 1 Γ ( ϕ ) s = 1 ϕ t ϕ ( t s 1 ) ϕ 1 ( ( 1 λ v ) Λ β S ( t + ϕ 1 ) I ( t + ϕ 1 ) Λ S ( t + ϕ 1 ) + γ R ( t + ϕ 1 ) ) , I ( t ) = I ( 0 ) + 1 Γ ( ϕ ) s = 1 ϕ t ϕ ( t s 1 ) ϕ 1 ( β S ( t + ϕ 1 ) I ( t + ϕ 1 ) ( σ + μ + δ ) × I ( t + ϕ 1 ) ) , R ( t ) = R ( 0 ) + 1 Γ ( ϕ ) s = 1 ϕ t ϕ ( t s 1 ) ϕ 1 ( λ v Λ + σ I ( t + ϕ 1 ) μ R ( t + ϕ 1 ) γ R ( t + ϕ 1 ) ) .
On the other hand, existence theory is one of the fundamental components of mathematical analysis, developed with the help of fixed-point theorems. The concept of fixed-point theory becomes valuable when solving problems in science and engineering involving nonlinear functional equations. These equations can be effectively addressed by transforming them into fixed-point problems. The fixed-point theory has proven to be a valuable tool in a wide range of mathematical analysis problems, including split feasibility problems, variational inequality problems, nonlinear optimization problems, the theory of fractals, equilibrium problems, complementarity problems, and selection and matching problems, as well as problems related to establishing the existence of solutions for integral and differential equations [27].
The study of stability in functional and differential equations, including exponential, Lyapunov, and Ulam–Hyers stability, is a significant focus in mathematical analysis. Ulam–Hyers stability, which connects exact and approximate solutions, was introduced by Ulam in 1940 and further developed by Hyers in 1941 for linear functional equations. Rassias extended this concept in 1978 to linear mappings, inspiring further exploration of its applications to ordinary differential equations, fractional differential equations, functional equations in multiple variables, first-order linear differential equations, advection–reaction diffusion systems, and fields like biology and economics [30].
Therefore, this paper will explore the existence theory and Ulam stability in detail, along with a comprehensive analysis that includes the non-negativity of the solutions, equilibrium points, basic reproduction number, stability of equilibrium points, sensitivity analysis, scenario analysis, numerical solutions, and the dynamics of the proposed model (9) using Caputo fractional-order difference operators.
The remainder of the paper is organized as follows: In Section 2, the existence theory for the solution of model (9) is established. Section 3 provides the Ulam stability of the proposed model (9). Section 4 addresses parameter estimation. Section 5 then provides a comprehensive basic analysis of the proposed model (9). In Section 6, the numerical solution is presented along with a discussion. Finally, Section 7 offers the concluding remarks of the study.

2. Existence and Uniqueness

Let us consider the functions of model (10) as follows:
f 1 ( t + ϕ 1 , S , I , R ) = ( 1 λ v ) Λ β S ( t + ϕ 1 ) I ( t + ϕ 1 ) Λ S ( t + ϕ 1 ) + γ R ( t + ϕ 1 ) , f 2 ( t + ϕ 1 , S , I , R ) = β S ( t + ϕ 1 ) I ( t + ϕ 1 ) ( σ + μ + δ ) I ( t + ϕ 1 ) , f 3 ( t + ϕ 1 , S , I , R ) = λ v Λ + σ I ( t + ϕ 1 ) ( μ + γ ) R ( t + ϕ 1 ) ,
then Equation (10) can be written in the following form:
f ( t ) = f ( 0 ) + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 G ( s + ϕ 1 , f ( s + ϕ 1 ) ) ,
where
f ( t ) = S ( t ) I ( t ) R ( t ) , f ( 0 ) = S 0 I 0 R 0 , G ( t + ϕ 1 , f ( t + ϕ 1 ) ) = f 1 ( t + ϕ 1 , S , I , R ) f 2 ( t + ϕ 1 , S , I , R ) f 3 ( t + ϕ 1 , S , I , R ) .
To establish the existence and uniqueness of solutions for the system (10), it is necessary to utilize the subsequent fixed-point theorems.
Theorem 1 
(see [40] (Banach contraction mapping principle)). The unique fixed-point property holds for a contraction mapping defined on a complete metric space, meaning that such a mapping has precisely one fixed point.
Theorem 2 
(see [40] (Brouwer fixed-point theorem)). If Φ : C R n C R n is a continuous mapping defined on a nonempty, bounded, closed, and convex set C, then Φ has at least one fixed point in C.
For the existence theory, we define Banach space as J = B 1 × B 2 × B 3 , where B i , ( i = 1 , 2 , 3 ) is the set of all real sequences ( S , I , R ) = { S ( t ) , I ( t ) , R ( t ) } t = ϕ 1 b + ϕ with the norm f = ( S , I , R = sup t [ ϕ 1 , b + ϕ ] N ϕ 1 [ | S ( t ) | + | I ( t ) | + | R ( t ) | ] . Then, J is a Banach space.
Let Φ : J J be an operator such that
( Φ f ) ( t ) = f ( 0 ) + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 G ( s + ϕ 1 , f ( s + ϕ 1 ) ) ,
If f ( t ) is a fixed point of the operator Φ , then it is evident that f ( t ) is a solution of (9).
To simplify computations, define C ϕ = max Γ ( t + 1 ) Γ ( ϕ + 1 ) Γ ( t ϕ + 1 ) = Γ ( b + ϕ + 1 ) Γ ( ϕ + 1 ) Γ ( b + 1 ) .
Theorem 3. 
Assume that
( A 1 ) :
There is a bounded function K : [ ϕ 1 , b + ϕ ] N ϕ 1 R such that | G ( t , f ) | K ( t ) | f | for all f J .
Then, the system (9) admits at least one solution on J as long as K < C ϕ , where K = max { K ( t ) : t [ ϕ 1 , b + ϕ ] N ϕ 1 }
Proof. 
Let M > 0 , and define the set W = f | [ ϕ 1 , b + ϕ ] N ϕ 1 R , f M . To establish this theorem, it suffices to verify that Φ maps W into W. For any f ( t ) W , we have
| ( Φ f ) ( t ) | K ( t ) Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 | f ( s + ϕ 1 ) | , K ( t ) f Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 , C ϕ K M .
Therefore, if | ( Φ f ) ( t ) | < M , then ( Φ f ) M , indicating that Φ maps W to itself. By the Brouwer fixed-point Theorem 2, Φ must have at least one fixed point that is a solution to (9). □
Theorem 4. 
Assume the following:
( A 2 ) :
There is a constant L > 0 such that | G ( t , f ) G ( t , g ) | L | f g | for every t [ ϕ 1 , b + ϕ ] N ϕ 1 for all f , g J .
Then, the system (9) has a unique solution on J as long as L < C ϕ .
Proof. 
Let f , g J , then for each t [ ϕ 1 , b + ϕ ] N ϕ 1 , we have
| ( Φ f ) ( t ) ( Φ g ) ( t ) | 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 | G ( s + ϕ 1 , f ( s + ϕ 1 ) ) G ( s + ϕ 1 , g ( s + ϕ 1 ) ) | , L Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 | f ( s + ϕ 1 ) ) g ( s + ϕ 1 ) ) | , L f g Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 , C ϕ L f g .
Hence, ( Φ f ) ( Φ g )     f g . This means that Φ is a contraction, and hence, by the Banach fixed-point Theorem 1, there exists a unique fixed point of Φ which is also a unique solution of the system (9). Therefore, the proof is complete. □

3. Ulam Stability

This section addresses the Ulam stability of the proposed model (9). Consider Equation (9) along with Equations (11) and (13), subject to the inequality
| C Δ ϕ f ( t ) G ( s + ϕ 1 , f ( s + ϕ 1 ) ) | ϵ , t [ 0 , b ] N 0 .
Definition 3 
([40]). Equation (9) is considered Ulam–Hyers stable if, for every ϵ > 0 , there exists a constant C f > 0 such that for any function f J satisfying inequality (15), a solution f ˜ J to (9) can be found, satisfying the following condition:
| f ( t ) f ˜ ( t ) | C f ϵ , t [ ϕ 1 , b + ϕ ] N ϕ 1 .
Equation (9) is generalized Ulam–Hyers stable if the function φ ( ϵ ) is substituted for the constant C f ϵ in inequality (15), where φ ( ϵ ) C ( R + , R + ) and φ ( 0 ) = 0 .
Remark 2. 
A function f J is a solution to inequality (15) if and only if there exists a function Υ : [ ϕ 1 , b + ϕ ] N ϕ 1 R , such that the following hold:
(i): 
| Υ ( s + ϕ 1 ) | ϵ , t [ 0 , b ] N 0 ,
(ii): 
C Δ ϕ f ( t ) = G ( s + ϕ 1 , f ( s + ϕ 1 ) ) + Υ ( s + ϕ 1 ) , t [ 0 , b ] N 0 .
Theorem 5. 
Assume that the hypothesis ( A 2 ) is satisfied. Let g J be a solution of inequality (15), and f J be a solution of the model (9). The proposed model (9) is Ulam–Hyers stable if L < C ϕ 1 .
Proof. 
Let f J be the solution of Equation (9), and let g J be the solution of inequality (15). Based on Remark 2(ii), it follows that
f ( t ) = f ( 0 ) + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 G ( s + ϕ 1 , f ( s + ϕ 1 ) ) + Υ ( s + ϕ 1 ) .
Using Remark 2(i), for t [ ϕ 1 , b + ϕ ] N ϕ 1 } , one can easily obtain the following:
f ( t ) f ( 0 ) + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 G ( s + ϕ 1 , f ( s + ϕ 1 ) )     C ϕ ϵ .
Next, for t [ ϕ 1 , b + ϕ ] N ϕ 1 , it follows that
| f ( t ) g ( t ) | = | f ( t ) g ( 0 ) + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 G ( s + ϕ 1 , g ( s + ϕ 1 ) ) | , C ϕ ϵ + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 × | G ( s + ϕ 1 , f ( s + ϕ 1 ) ) G ( s + ϕ 1 , g ( s + ϕ 1 ) ) | , C ϕ ϵ + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 L f g .
This leads to
f g C ϕ 1 L C ϕ ϵ = C f ϵ .
Thus, if C f = C ϕ 1 L C ϕ > 0 , it follows from Definition 3 that the solution to model (9) is Ulam–Hyers stable. □
Remark 3. 
By selecting φ ( ϵ ) = C f ϵ , where φ ( 0 ) = 0 , it can be concluded that the solution of model (9) satisfies the criteria for generalized Ulam–Hyers stability.

4. Parameters Estimation

The initial conditions for the proposed model are set with the susceptible population at S ( 0 ) = 0.92 , the infected population at I ( 0 ) = 0.06 , and the recovered population at R ( 0 ) = 0.02 , following the values provided in [12]. The birth rate, Λ = 0.01 [12], is kept constant across all cases to stabilize population influx, isolating the effects of disease parameters for a clearer analysis of disease dynamics without birth rate fluctuations.
  • Case 1. This scenario features a moderate infection rate ( β = 0.6 ) and recovery rate ( σ = 0.3 ), representing a disease that can spread within the population but lacks an exceedingly high transmission rate or fatality impact ( μ = 0.015 ). With a disease-induced mortality rate of δ = 0.01 , the risk of death is low but present. The moderate immunity loss rate ( γ = 0.02 ) suggests that individuals may become susceptible again, though not immediately, pointing to diseases that have limited reinfection rates and maintain a relatively stable population effect. This case could represent diseases where reinfection is uncommon, and overall population health remains steady.
  • Case 2. This case models a highly infectious disease scenario characterized by a higher infection rate ( β = 0.8 ) alongside a substantial recovery rate ( σ = 0.4 ), reflecting a rapidly spreading disease with a reasonable chance of recovery. However, the mortality rate ( δ = 0.03 ) is elevated, underscoring the disease’s severity. The increased vaccination rate ( λ v = 0.3 ) suggests a need to control transmission more aggressively. The parameter μ = 0.02 represents the natural death rate, while γ = 0.01 . This setup mirrors diseases such as measles in low-vaccination communities or pertussis, where high transmission and mortality risks necessitate vigilant infection control measures.
  • Case 3. This case represents a recurrent disease, characterized by a heightened immunity loss rate ( γ = 0.04 ), leading to individuals quickly losing immunity and re-entering the susceptible population. With an infection rate of β = 0.7 , a low mortality rate ( δ = 0.005 ), and a moderate recovery rate ( σ = 0.35 ), the disease poses minimal risk of death but persists due to frequent reinfections. The low natural death rate ( μ = 0.01 ) and vaccination rate ( λ v = 0.1 ) reflect minimal demographic impact and vaccination efforts. Diseases like RSV or rotavirus, where reinfections are common but rarely fatal, exemplify this case, highlighting the ongoing healthcare burden posed by recurrent infections.

5. Basic Analysis of the Model

This section discusses the non-negativity of solutions, equilibrium points, basic reproduction number, stability, sensitivity analysis, and scenario analysis.

5.1. Non-Negativity of the Solutions

In this population dynamics model, all solutions of system (9) are required to stay non-negative. The theorem below demonstrates the non-negativity of the solutions for system (9).
Theorem 6. 
All solutions of system (9) with non-negative initial values remain non-negative.
Proof. 
Given system (9) and the non-negativity of the initial values S ( 0 ) 0 , I ( 0 ) 0 , and R ( 0 ) 0 , it follows that
S ( t ) = S ( 0 ) + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 f 1 ( t + ϕ 1 , S , I , R ) 0 , I ( t ) = I ( 0 ) + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 f 2 ( t + ϕ 1 , S , I , R ) 0 , R ( t ) = R ( 0 ) + 1 Γ ( ϕ ) s = 0 t ϕ ( t s 1 ) ϕ 1 f 3 ( t + ϕ 1 , S , I , R ) 0 ,
where f 1 , f 2 and f 3 are defined as in (11). Thus, all solutions of system (9) remain non-negative. □

5.2. Disease-Free Equilibrium

A disease-free equilibrium describes a condition within a population or community where there are no cases of a specific disease. In this state, any introduction of the infection will not lead to an outbreak because the factors required for the disease’s propagation are not present. Reaching a disease-free equilibrium is often a primary objective in public health initiatives, which may include strategies like vaccination, effective treatments, and other preventive actions aimed at eliminating the disease from the community. To determine this equilibrium, the fixed points are identified by equating the right-hand side of the system in Equation (9) to zero as follows:
( 1 λ v ) Λ β S ( t + ϕ 1 ) I ( t + ϕ 1 ) Λ S ( t + ϕ 1 ) + γ R ( t + ϕ 1 ) = 0 , β S ( t + ϕ 1 ) I ( t + ϕ 1 ) ( σ + μ + δ ) I ( t + ϕ 1 ) = 0 , λ v Λ + σ I ( t + ϕ 1 ) ( μ + γ ) R ( t + ϕ 1 ) = 0 .
In the absence of disease, it is assumed that I = R = 0 . Under these circumstances, the disease-free equilibrium point, denoted as E 0 is E 0 = ( S 0 , I 0 , R 0 ) = ( γ + μ μ λ v γ + μ , 0 , Λ λ v γ + μ ) .

5.3. Endemic Equilibrium

An endemic equilibrium represents a condition in which a disease persists within a population, maintaining a stable prevalence over time. In this scenario, the infection continues to transmit among individuals, leading to a steady number of cases. For the proposed model (9), the endemic equilibrium point is defined as E = ( S , I , R ) , where
S = δ + μ + σ β , I = Λ β γ + β μ γ δ γ μ γ σ δ μ μ 2 μ σ β μ λ v β γ δ + γ μ + δ μ + μ 2 + μ σ , R = β Λ σ δ Λ σ Λ μ σ Λ σ 2 + β δ Λ λ v + β Λ μ λ v β γ δ + γ μ + δ μ + μ 2 + μ σ .
The next subsection will present the derivation of the basic reproduction number ( R 0 ) .

5.4. Basic Reproduction Number ( R 0 )

The basic reproduction number ( R 0 ) is a crucial metric for assessing the potential spread of an epidemic. Mathematically, R 0 is defined as the spectral radius, which is the largest absolute eigenvalue of the next-generation matrix [41]. Using the next-generation matrix method, the following formula for R 0 is obtained for the proposed model:
R 0 = β S 0 ( μ + σ + δ ) = β ( γ + μ μ λ v ) ( γ + μ ) ( μ + σ + δ ) .
Equation (22) defines the basic reproduction number for the proposed model (9). The following subsection will discuss the stability analysis of the equilibrium points.

5.5. Stability Analysis of Fixed Points

Theorem 7. 
The disease-free equilibrium state E 0 is locally asymptotically stable when R 0 < 1 , and unstable when R 0 1 .
Proof. 
The Jacobian matrix J at a point ( S , I , R ) of the proposed model (9) is defined as follows:
J = β I Λ β S γ β I δ μ σ + β S 0 0 σ γ μ .
The characteristic equation of the Jacobian matrix (23) at the disease-free equilibrium point is det ( J ( E 0 ) λ I ) = 0 , where I is the identity matrix and λ is the eigenvalue. Then, one can obtain the following eigenvalues:
λ 1 = Λ , λ 2 = ( γ + μ ) , λ 3 = δ + μ + σ β ( γ + μ μ λ v ) γ + μ .
The term δ + μ + σ β ( γ + μ μ λ v ) γ + μ must be positive for stability, which corresponds to β ( γ + μ μ λ v ) ( γ + μ ) ( μ + σ + δ ) < 1 or equivalently, R 0 < 1 , which guarantees that all eigenvalues have negative real parts, ensuring stability at the disease-free equilibrium. □
Theorem 8. 
The disease-endemic equilibrium state E is locally asymptotically stable if R 0 1 .
Proof. 
The characteristic equation for the Jacobian matrix (23) at the disease-endemic equilibrium point E = ( 0.5416 ,   0.0256 ,   0.0428 ) for Case 1 is det ( J ( E ) λ I ) = 0 . Solving this equation yields the following eigenvalues:
λ 1 = 0.0222148 0.0684798 i , λ 2 = 0.0222148 + 0.0684798 i , λ 3 = 0.0159657 .
Similarly, the equilibrium points for Cases 2 and 3 are E = ( 0.5625 ,   0.0075 ,   0.1 ) and E = ( 0.5214 ,   0.0539 ,   0.02 ) , respectively. The corresponding eigenvalues for E = ( 0.5625 ,   0.0075 ,   0.1 ) are
λ 1 = 0.0214827 , λ 2 = 0.0122587 0.0500303 i , λ 3 = 0.0122587 + 0.0500303 i
and for E = ( 0.5214 ,   0.0539 ,   0.02 ) are
λ 1 = 0.0436216 0.115549 i , λ 2 = 0.0436216 + 0.115549 i , λ 3 = 0.0105215 .
This analysis confirms that the disease-endemic equilibrium states for all cases are stable, as the real parts of all eigenvalues are negative. □

5.6. R 0 Sensitivity Analysis

Analyzing the sensitivity of the basic reproduction number ( R 0 ) to its defining parameters is essential for understanding each parameter’s influence on the spread or control of infection. The sensitivity index of R 0 with respect to any given parameter ( ) can be expressed as
S ( ) = ( ) R 0 × R 0 ( ) .
The sensitivity values provided in Table 1 highlight the magnitude and direction of each parameter’s impact on R 0 , shedding light on which factors contribute most significantly to the infection dynamics.
These sensitivity values are depicted in Figure 2. The parameter most sensitive to increasing R 0 across all three cases is the infection coefficient ( β ) . A sensitivity value of 1.0 for β suggests that a small change in this parameter will cause a significant increase in R 0 , underscoring the strong influence of β on the basic reproduction number.
The second most sensitive parameter with a positive impact on R 0 is the rate at which recovered individuals lose immunity ( γ ) , across all three cases. However, its effect varies slightly in each case, with sensitivity values of 0.039 , 0.083 , and 0.016 for cases 1, 2, and 3, respectively. The sensitivity values of 0.085 , 0.127 , and 0.043 for the natural death rate ( μ ) indicate that small changes in μ produce only minor changes in R 0 , with the negative sign signifying an inverse relationship. This means that an increase in μ will slightly decrease R 0 .
For the recovery rate ( σ ) , the sensitivity values are 0.923 , 0.888 , and 0.958 for cases 1, 2, and 3, respectively, indicating that small changes in σ will lead to significant decreases in R 0 . The negative sign again signifies an inverse relationship, where an increase in σ strongly reduces R 0 . Similarly, the sensitivity values of 0.030 , 0.066 , and 0.013 for the disease-induced death rate ( δ ) suggest that a small change in δ will result in a minimal decrease in R 0 . Finally, the sensitivity values for the proportion of the population vaccinated ( λ v ) are 0.068 , 0.250 , and 0.020 for cases 1, 2, and 3, respectively, indicating that changes in λ v will lead to minimal decreases in R 0 . The negative sign reflects an inverse relationship, where increasing λ v slightly reduces R 0 .

5.7. Scenario Analysis

In addition to analyzing the sensitivity of R 0 to individual parameters, examining the effect of the parameter pairs offers deeper understanding into the complex interactions within the model. By assessing these pairwise parameter influences, we can better understand how multiple factors jointly affect the disease dynamics, particularly in scenarios where single-parameter variations do not fully capture the model’s sensitivity. Figure 3, Figure 4 and Figure 5 depict the contour plots of R 0 as influenced by combinations of parameter pairs, illustrating the basic reproduction number’s response in Case 1, Case 2, and Case 3, respectively.

6. Numerical Solution and Discussion

The model in Equation (9) is used to present the numerical results. The numerical approach, shown in Equation (25) and derived by substituting s + ϕ = j into system (10), is applied to compute the numerical solutions. The values for λ v , Λ , β , σ , μ , σ , and γ are taken from Table 2, and the initial values for the susceptible ( S ), infected ( I ), and recovered ( R ) populations are provided in Section 4:
S ( n ) = S ( 0 ) + 1 Γ ( ϕ ) j = 1 n Γ ( n j + ϕ ) Γ ( n j + 1 ) ( ( 1 λ v ) Λ β S ( j 1 ) I ( j 1 ) Λ S ( j 1 ) + γ R ( j 1 ) ) , I ( n ) = I ( 0 ) + 1 Γ ( ϕ ) j = 1 n Γ ( n j + ϕ ) Γ ( n j + 1 ) ( β S ( j 1 ) I ( j 1 ) ( σ + μ + δ ) × I ( j 1 ) ) , R ( n ) = R ( 0 ) + 1 Γ ( ϕ ) j = 1 n Γ ( n j + ϕ ) Γ ( n j + 1 ) ( λ v Λ + σ I ( j 1 ) μ R ( j 1 ) γ R ( j 1 ) ) .
  • Case 1. Figure 6 illustrates the model’s dynamics for Case 1 parameter values. In Figure 6a with ϕ = 0.7 , the susceptible population ( S ) declines gradually as children become infected, and then rises slowly as immunity fades, bringing individuals back to susceptibility. This slower change is due to a moderate infection rate ( β = 0.6 ) and a moderate immunity loss rate ( γ = 0.02 ). The infected population ( I ) shows moderate peaks, reflecting a stable pattern influenced by the balanced shifts in the susceptible population. Meanwhile, the recovered population ( R ) increases steadily after outbreaks and gradually declines as immunity wanes, demonstrating stable recovery and immunity loss dynamics.
In Figure 6b with ϕ = 0.8 , shifts in the susceptible population ( S ) become quicker. The decrease during outbreaks happens faster, and the increase due to waning immunity accelerates, indicating a more responsive dynamic than at ϕ = 0.7 . The infected population ( I ) exhibits a sharper peak, responding more quickly to changes in the susceptible population. Likewise, the recovered population ( R ) experiences more noticeable rises after outbreaks, followed by faster declines as immunity is lost, indicating heightened fluctuation and a more dynamic recovery process.
In Figure 6c with ϕ = 0.9 , changes in the susceptible population ( S ) are even more rapid. The drop during outbreaks is significant, and the increase as immunity fades is swift, reflecting a high level of responsiveness. The infected population ( I ) displays a very sharp peak, showing rapid response to the infectious spread. The recovered population ( R ) experiences substantial increases post-outbreaks and quick decreases as immunity is lost, revealing marked fluctuation and a highly dynamic recovery pattern.
Finally, in Figure 6d with ϕ = 1.0 , the changes in the susceptible population ( S ) occur rapidly, with steep declines during outbreaks followed by sharp increases as immunity wanes, reflecting a highly responsive scenario. The infected population ( I ) displays significant fluctuations, with pronounced peaks and valleys that correspond to changes in the susceptible population. These variations highlight the dynamic interplay between infection and reinfection. The recovered population ( R ) increases quickly after outbreaks but decreases sharply as immunity fades, showing rapid shifts in recovery and reinfection dynamics, with observable oscillations that gradually stabilize over time.
From the behaviors depicted in Figure 6, as ϕ progresses from 0.7 to 1.0, the speed and intensity of changes in the susceptible, infected, and recovered populations increase. Lower ϕ values result in smoother shifts and quicker stabilization, while higher values lead to faster and more pronounced fluctuations.
Figure 6. Case 1: The relative population dynamics of childhood disease for various values of ϕ , where λ v = 0.15 , Λ = 0.01 , β = 0.6 , σ = 0.3 ; μ = 0.015 , δ = 0.01 and γ = 0.02 .
Figure 6. Case 1: The relative population dynamics of childhood disease for various values of ϕ , where λ v = 0.15 , Λ = 0.01 , β = 0.6 , σ = 0.3 ; μ = 0.015 , δ = 0.01 and γ = 0.02 .
Fractalfract 09 00055 g006
Figure 7 presents phase-plane portraits showing the interactions among susceptible ( S ), infected ( I ), and recovered ( R ) populations in the proposed disease model (9) of Case 1 for various values of the fractional-order ϕ (0.7, 0.8, 0.9, and 1.0).
For ϕ = 0.7 , as shown in Figure 7a–c, the dynamics exhibit closed trajectories, indicating stable oscillatory behavior around an equilibrium point. This implies that the system stabilizes with sustained but limited oscillations. When ϕ = 0.8 , as seen in Figure 7d–f, the trajectories indicate mild, stable oscillations around the equilibrium, suggesting that the populations quickly stabilize, leading to a steady state with minimal long-term fluctuations.
At ϕ = 0.9 , displayed in Figure 7g–i, the phase trajectories develop inward spirals, indicating dampened oscillations that take longer to settle into equilibrium, suggesting prolonged cycles of infection and recovery. Finally, when ϕ = 1.0 , as shown in Figure 7j–l, the trajectories exhibit more pronounced inward spirals, indicating sustained oscillatory dynamics with extended fluctuations before reaching equilibrium. As ϕ approaches 1, the system exhibits more persistent oscillations, reflecting prolonged dynamics before stabilization.
  • Case 2. Figure 8 illustrates the dynamics of a highly infectious childhood disease with a high infection rate ( β = 0.8 ) and substantial recovery rate ( σ = 0.4 ). The figure includes four subplots, each showing a different fractional-order value for ϕ : 0.7, 0.8, 0.9, and 1.0. These subplots demonstrate how varying ϕ influences the susceptible ( S ), infected ( I ), and recovered ( R ) populations over time, while considering the role of immunity loss ( γ = 0.01 ).
In Figure 8a with ϕ = 0.7 , the susceptible population ( S ) declines gradually during outbreaks, driven by the high infection rate ( β = 0.8 ), while immunity loss ( γ = 0.01 ) replenishes the susceptible class at a moderate rate. The infected population ( I ) displays moderate peaks, consistent with a relatively slower disease spread under fractional-order dynamics. The recovered population ( R ) grows steadily due to the high recovery rate ( σ = 0.4 ) but decreases slowly as individuals lose immunity.
As ϕ increases to 0.8 (Figure 8b), the dynamics of ( S ) become sharper, with more pronounced declines during outbreaks and faster rebounds due to the combined effects of immunity loss ( γ = 0.01 ) and the fractional memory effect. The infected population ( I ) peaks more sharply, reflecting faster transmission dynamics, while the recovered population ( R ) begins to exhibit visible oscillations, driven by individuals transitioning back to the susceptible class.
For ϕ = 0.9 (Figure 8c), the susceptible population ( S ) undergoes even more rapid shifts, with steep declines during outbreaks and quicker rebounds. The infected population ( I ) reaches higher and more pronounced peaks due to accelerated disease spread. The recovered population ( R ) experiences significant oscillations, highlighting the interplay between recovery, immunity loss ( γ ), and heightened disease dynamics as the system approaches classical-order behavior.
At ϕ = 1.0 (Figure 8d), the dynamics are most intense, reflecting the absence of fractional-order memory effects. The susceptible population ( S ) drops steeply during outbreaks, with immunity loss ( γ = 0.01 ) rapidly replenishing it. The infected population ( I ) displays very sharp peaks, indicating immediate and pronounced responses, while the recovered population ( R ) rises and falls abruptly, showcasing substantial fluctuations driven by the faster interplay between recovery and immunity loss.
Figure 8. Case 2: The relative population dynamics of childhood disease for various values of ϕ , where λ v = 0.3 , Λ = 0.01 , β = 0.8 , σ = 0.4 ; μ = 0.02 , δ = 0.03 and γ = 0.01 .
Figure 8. Case 2: The relative population dynamics of childhood disease for various values of ϕ , where λ v = 0.3 , Λ = 0.01 , β = 0.8 , σ = 0.4 ; μ = 0.02 , δ = 0.03 and γ = 0.01 .
Fractalfract 09 00055 g008
As ϕ increases from 0.7 to 1.0, as shown in Figure 8, the system transitions from smoother fractional-order dynamics to more volatile and intense classical-order behavior. The role of γ becomes increasingly significant, as immunity loss drives faster replenishment of the susceptible population and amplifies oscillations in the recovered population. This model reflects a highly infectious childhood disease with high infection and recovery rates, moderate mortality ( δ = 0.03 ), and robust vaccination efforts ( λ v = 0.3 ).
Figure 9 displays phase-plane plots of the susceptible ( S ), infected ( I ), and recovered ( R ) populations for Case 2 across different fractional-order values of ϕ (0.7, 0.8, 0.9, and 1.0). At ϕ = 0.7 , the trajectories show rapid convergence with stable, closed loops, indicating the quick stabilization of infection and recovery cycles. As ϕ increases to 0.8, the system exhibits more prolonged oscillations, reflecting slower convergence and longer stabilization periods.
With ϕ = 0.9 , oscillations become more intricate, suggesting even slower convergence, with recurring cycles that delay stabilization. At ϕ = 1.0 , the trajectories display persistent inward spirals and loops, indicating extended oscillations and a marked delay in reaching equilibrium. This progression demonstrates that higher ϕ values result in more pronounced oscillatory behavior and slower stabilization in the population dynamics.
  • Case 3. The dynamics of a recurrent disease with high immunity loss ( γ = 0.04 ), a moderate infection rate ( β = 0.7 ), a low mortality rate ( δ = 0.005 ), and a moderate recovery rate ( σ = 0.35 ) are illustrated in Figure 10. This figure includes four subplots, each corresponding to a different value of the fractional order ϕ : 0.7, 0.8, 0.9, and 1.0. The following parameters remain fixed: vaccination rate ( λ v = 0.1 ), natural death rate ( μ = 0.01 ), and recruitment rate ( Λ = 0.01 ).
For ϕ = 0.7 , as shown in Figure 10a, the susceptible population ( S ) gradually decreases as individuals become infected. However, due to the high immunity loss rate ( γ = 0.04 ), the susceptible class is replenished over time, resulting in a relatively smooth transition. The infected population ( I ) stabilizes at a low level, reflecting the moderate infection rate, while the recovered population ( R ) increases slowly but does not accumulate significantly due to frequent immunity loss, leading to stable but low levels of recovery.
As ϕ increases to 0.8, the dynamics of the system become more pronounced as shown in Figure 10b. The susceptible population ( S ) declines more rapidly during outbreaks and stabilizes at a lower equilibrium level compared to ϕ = 0.7 . The infected population ( I ) reaches higher peak levels before stabilizing, and the recovered population ( R ) grows more quickly, though still constrained by the immunity loss rate, resulting in a higher equilibrium level than at ϕ = 0.7 .
Figure 10. Case 3: The relative population dynamics of childhood disease for various values of ϕ , where λ v = 0.1 , Λ = 0.01 , β = 0.7 , σ = 0.35 ; μ = 0.01 , δ = 0.005 and γ = 0.04 .
Figure 10. Case 3: The relative population dynamics of childhood disease for various values of ϕ , where λ v = 0.1 , Λ = 0.01 , β = 0.7 , σ = 0.35 ; μ = 0.01 , δ = 0.005 and γ = 0.04 .
Fractalfract 09 00055 g010
At ϕ = 0.9 , the transitions between population classes become sharper as depicted in Figure 10c. The susceptible population ( S ) declines significantly during outbreaks before stabilizing around 0.5. The infected population ( I ) exhibits more pronounced peaks, reflecting an increased intensity of infection. Meanwhile, the recovered population ( R ) oscillates more dramatically and stabilizes at a higher level compared to ϕ = 0.8 , although still impacted by the ongoing immunity loss.
With ϕ = 1.0 , the dynamics of the disease reach their fastest pace as shown in Figure 10d. The susceptible population ( S ) drops sharply during outbreaks and then stabilizes around 0.52. The infected population ( I ) shows a brief, sharp peak before rapidly stabilizing at a low level. The recovered population ( R ) increases initially, then stabilizes close to the level of the susceptible population. The higher fractional order leads to faster transitions between population classes and quicker convergence to equilibrium.
As ϕ increases from 0.7 to 1.0, the dynamics of the system become faster and more pronounced. Higher fractional orders result in quicker equilibrium and more abrupt transitions between susceptible, infected, and recovered populations. In contrast, lower values of ϕ correspond to smoother and slower transitions, reflecting more gradual changes in population levels. The model captures the interplay between moderate infection, high immunity loss, and recovery, simulating recurrent diseases like RSV or rotavirus, which exhibit low mortality but persist due to frequent reinfections.
The phase-plane plots in Figure 11 depict the interactions among susceptible ( S ), infected ( I ), and recovered ( R ) populations of Case 3 for varying fractional-order values of ϕ (0.7, 0.8, 0.9, and 1.0) in the disease model (9). At ϕ = 0.7 , the trajectories converge directly to equilibrium, indicating quick stabilization. With ϕ = 0.8 , closed loops appear, reflecting mild oscillations that stabilize quickly.
As ϕ increases to 0.9, the trajectories form inward spirals, showing prolonged dampened oscillations, suggesting a slower path to equilibrium with potential recurring outbreaks. At ϕ = 1.0 , sustained inward spirals indicate extended oscillations, modeling dynamics with delayed stabilization and higher susceptibility to reinfection.
Figure 12 compares the MATLAB (R2023b) ode45 solver with the discrete fractional-order numerical method (25) for ϕ = 1.0 . Since the exact solutions are not available, the ode45 solutions are used as the reference. Figure 12a presents the comparison for Case 1, Figure 12b illustrates the comparison for Case 2, and Figure 12c depicts the comparison for Case 3. In all cases, the numerical solutions from the discrete fractional-order method align well with those from the ode45 solver, demonstrating a good level of agreement.
This indicates that the discrete fractional-order method (25) provides accurate approximations in the absence of exact solutions. These results collectively highlight the effectiveness and reliability of the method for solving fractional-order differential equations and suggest that it could serve as a viable alternative to traditional solvers like ode45 in similar applications.

7. Conclusions

This study presents a novel SIRS model for recurrent childhood diseases based on a discrete fractional-order operator. The discrete fractional approach is crucial, as it enables capturing the memory effects and non-locality often present in biological processes and most of the disease data are in discrete form, thereby providing a more accurate reflection of disease dynamics over time. The model incorporates a reinfection rate to account for the dynamics of immunity loss in recovered individuals. This comprehensive framework enhances understanding of the transmission and recovery processes in childhood diseases.
The existence theory is established using Brower’s fixed-point theorem and the Banach contraction principle, providing a robust mathematical foundation. The Ulam stability of the model’s solution is also demonstrated, ensuring its reliability under small perturbations. Sensitivity analysis further highlights the impact of parameter variations, and the basic reproduction number R 0 is derived to assess local stability at the disease-free and endemic equilibria. The analysis shows that the disease-free equilibrium remains stable for R 0 < 1 , while the endemic equilibrium becomes stable for R 0 > 1 and unstable otherwise.
Numerical simulations demonstrate the model’s effectiveness in capturing realistic disease patterns, underscoring the value of discrete fractional-order models for predicting long-term trends and outbreak potential. Additionally, fractional orders show greater stability compared to integer-order models. This work contributes meaningfully to disease management policies for childhood illnesses, offering a refined mathematical framework to support public health decision-making. Future directions for this research include investigating the proposed model under alternative differential operators to assess its behavior and applicability within varied mathematical frameworks. Additionally, the model can be extended to integrate critical parameters of childhood diseases, enhancing its effectiveness in analyzing disease dynamics and supporting intervention strategies.

Author Contributions

Conceptualization, Y.A.M.; Software, B.M.; Formal analysis, M.R., A.A. and N.H.E.E.; Investigation, A.A. and N.H.E.E.; Writing—original draft, Z.A.; Writing—review and editing, Y.A.M. and K.A.; Project administration, K.A.; Funding acquisition, M.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data are contained within the article.

Acknowledgments

This research was funded by the Deanship of Graduate Studies and Scientific Research at Qassim University, grant number (QU-APC-2025).

Conflicts of Interest

Authors declare that there are no conflicts of interest with this paper.

References

  1. Martcheva, M. An Introduction to Mathematical Epidemiology; Springer: New York, NY, USA, 2015; Volume 61. [Google Scholar]
  2. Roberts, M.; Heesterbeek, J. Mathematical Models in Epidemiology; Springer: New York, NY, USA, 2019; Volume 69. [Google Scholar]
  3. Momani, S.; Kumar, R.; Srivastava, H.M.; Kumar, S.; Hadid, S. A chaos study of fractional SIR epidemic model of childhood diseases. Results Phys. 2021, 27, 104422. [Google Scholar] [CrossRef]
  4. Darti, I.; Suryanto, A. Dynamics of a SIR Epidemic Model of Childhood Diseases with a Saturated Incidence Rate: Continuous Model and Its Nonstandard Finite Difference Discretization. Mathematics 2020, 8, 1459. [Google Scholar] [CrossRef]
  5. Alter, S.J.; Bennett, J.S.; Koranyi, K.; Kreppel, A.; Simon, R. Common Childhood Viral Infections. Curr. Probl. Pediatr. Adolesc. Health Care 2015, 45, 21–53. [Google Scholar] [CrossRef] [PubMed]
  6. Makinde, O.D. Adomian decomposition approach to a SIR epidemic model with constant vaccination strategy. Appl. Math. Comput. 2007, 184, 842–848. [Google Scholar] [CrossRef]
  7. Platt, L.; Thun, M.; Harriman, K.; Winter, K. A Population-Based Study of Recurrent Symptomatic Bordetella pertussis Infections in Children in California, 2010–2015. Clin. Infect. Dis. 2017, 65, 2099–2104. [Google Scholar] [CrossRef] [PubMed]
  8. Crawford, S.; Ramani, S.; Tate, J.; Parashar, U.D.; Svensson, L.; Hagbom, M.; Franco, M.A.; Greenberg, H.B.; O’Ryan, M.; Kang, G.; et al. Rotavirus infection. Nat. Rev. Dis. Primers 2017, 3, 17083. [Google Scholar] [CrossRef]
  9. Streng, A.; Prifert, C.; Weissbrich, B.; Sauerbrei, A.; Krumbholz, A.; Schmidt-Ott, R.; Liese, J.G. Similar severity of influenza primary and re-infections in pre-school children requiring outpatient treatment due to febrile acute respiratory illness: Prospective, multicentre surveillance study (2013–2015). BMC Infect. Dis. 2022, 22, 12. [Google Scholar] [CrossRef] [PubMed]
  10. Kutsaya, A.; Teros-Jaakkola, T.; Kakkola, L.; Toivonen, L.; Peltola, V.; Waris, M.; Julkunen, I. Prospective clinical and serological follow-up in early childhood reveals a high rate of subclinical RSV infection and a relatively high reinfection rate within the first 3 years of life. Epidemiol. Infect. 2016, 144, 1622–1633. [Google Scholar] [CrossRef] [PubMed]
  11. Shah, M.P.; Hall, A.J. Norovirus Illnesses in Children and Adolescents. Infect. Dis. Clin. N. Am. 2018, 32, 103–118. [Google Scholar] [CrossRef]
  12. Selvam, A.G.M.; Vianny, D.A.; Jacintha, M. Stability in a fractional order SIR epidemic model of childhood diseases with discretization. J. Phys. Conf. Ser. 2018, 1139, 012009. [Google Scholar] [CrossRef]
  13. de Oliveira, E.C.; Machado, J.A.T. A Review of Definitions for Fractional Derivatives and Integral. Math. Probl. Eng. 2014, 2014, 238459. [Google Scholar] [CrossRef]
  14. Hilfer, R. Applications of Fractional Calculus in Physics; World Scientific: Singapore, 2000. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  16. Das, S. Application of Generalized Fractional Calculus in Electrical Circuit Analysis and Electromagnetics. In Functional Fractional Calculus; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar] [CrossRef]
  17. Magin, R.L. Fractional calculus in bioengineering. Crit. Rev. Biomed. Eng. 2004, 32, 1–104. [Google Scholar] [CrossRef] [PubMed]
  18. Oldham, K. Fractional differential equations in electrochemistry. Adv. Eng. Softw. 2010, 41, 9–12. [Google Scholar] [CrossRef]
  19. Ali, Z.; Nia, S.N.; Rabiei, F.; Shah, K.; Tan, M.K. A semi-analytical approach for the solution of time-fractional Navier-Stokes equation. Adv. Math. Phys. 2021, 2021, 5547804. [Google Scholar] [CrossRef]
  20. Tarasov, V.E. Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
  21. Alqahtani, A.M.; Akram, S.; Ahmad, J.; Aldwoah, K.A.; ur Rahman, M. Stochastic wave solutions of fractional Radhakrishnan-Kundu-Lakshmanan equation arising in optical fibers with their sensitivity analysis. J. Opt. 2024. [Google Scholar] [CrossRef]
  22. Matušu, R. Application of fractional order calculus to control theory. Int. J. Math. Model. Methods Appl. Sci. 2011, 5, 1162–1169. [Google Scholar]
  23. Meilanov, R.P.; Magomedov, R.A. Thermodynamics in Fractional Calculus. J. Eng. Phys. Thermophy 2014, 87, 1521–1531. [Google Scholar] [CrossRef]
  24. Meral, F.; Royston, T.; Magin, R. Fractional calculus in viscoelasticity: An experimental study. Commun. Nonlinear Sci. Numer. Simul. 2010, 15, 939–945. [Google Scholar] [CrossRef]
  25. Abdulwahhab, O.W.; Abbas, N.H. A new method to tune a fractional-order PID controller for a twin rotor aerodynamic system. Arab. J. Sci. Eng. 2017, 42, 5179–5189. [Google Scholar] [CrossRef]
  26. Jiang, Y.; Zhang, B.; Shu, X.; Wei, Z. Fractional-order autonomous circuits with order larger than one. J. Adv. Res. 2020, 25, 217–225. [Google Scholar] [CrossRef] [PubMed]
  27. Ali, Z. Theoretical and Computational Study of Fractional-Order Mathematical Models for Infectious Diseases. Ph.D. Thesis, Monash University, Clayton, VIC, Australia, 2023. [Google Scholar] [CrossRef]
  28. Shah, N.A.; Vieru, D.; Fetecau, C. Effects of the fractional order and magnetic field on the blood flow in cylindrical domains. J. Magn. Magn. Mater. 2016, 409, 10–19. [Google Scholar] [CrossRef]
  29. Can, N.H.; Jafari, H.; Ncube, M.N. Fractional calculus in data fitting. Alex. Eng. J. 2020, 59, 3269–3274. [Google Scholar] [CrossRef]
  30. Madani, Y.A.; Rabih, M.N.A.; Alqarni, F.A.; Ali, Z.; Aldwoah, K.A.; Hleili, M. Existence, Uniqueness, and Stability of a Nonlinear Tripled Fractional Order Differential System. Fractal Fract. 2024, 8, 416. [Google Scholar] [CrossRef]
  31. Adel, M.; Khader, M.M.; Algelany, S.; Aldwoah, K. An Accurate Approach to Simulate the Fractional Delay Differential Equations. Fractal Fract. 2023, 7, 671. [Google Scholar] [CrossRef]
  32. Ali, Z.; Rabiei, F.; Shah, K.; Majid, Z.A. Dynamics of SIR mathematical model for COVID-19 outbreak in Pakistan under fractal-fractional derivative. Fractals 2021, 29, 2150120. [Google Scholar] [CrossRef]
  33. Ostalczyk, P. Discrete Fractional Calculus: Applications in Control and Image Processing; World Scientific: Singapore, 2015; Volume 4. [Google Scholar]
  34. Al-Khedhairi, A.; Elsadany, A.A.; Elsonbaty, A. On the dynamics of a discrete fractional-order cournot-bertrand competition duopoly game. Math. Probl. Eng. 2022, 2022, 8249215. [Google Scholar] [CrossRef]
  35. Elsonbaty, A.; Elsadany, A.A. On discrete fractional-order Lotka-Volterra model based on the Caputo difference discrete operator. Math. Sci. 2023, 17, 67–79. [Google Scholar] [CrossRef]
  36. Huang, L.-L.; Park, J.H.; Wu, G.-C.; Mo, Z.-W. Variable-order fractional discrete-time recurrent neural networks. J. Comput. Appl. Math. 2020, 370, 112633. [Google Scholar] [CrossRef]
  37. He, Z.-Y.; Abbes, A.; Jahanshahi, H.; Alotaibi, N.D.; Wang, Y. Fractional-Order Discrete-Time SIR Epidemic Model with Vaccination: Chaos and Complexity. Mathematics 2022, 10, 165. [Google Scholar] [CrossRef]
  38. Tarasov, V.E.; Tarasova, V.V. Long and short memory in economics: Fractional-order difference and differentiation. arXiv 2016, arXiv:1612.07903. [Google Scholar] [CrossRef]
  39. Peng, Y.; He, S.; Sun, K. Chaos in the discrete memristor-based system with fractional-order difference. Results Phys. 2021, 24, 104106. [Google Scholar] [CrossRef]
  40. Chen, F.; Zhou, Y. Existence and Ulam Stability of Solutions for Discrete Fractional Boundary Value Problem. Discret. Dyn. Nat. Soc. 2013, 2013, 459161. [Google Scholar] [CrossRef]
  41. van den Driessche, P. Reproduction Numbers of Infectious Disease Models. Infect. Dis. Model. 2017, 2, 288–303. [Google Scholar] [CrossRef]
Figure 1. Flowchart illustrating the dynamics of model (1).
Figure 1. Flowchart illustrating the dynamics of model (1).
Fractalfract 09 00055 g001
Figure 2. The sensitivity of R 0 to each parameter in the model.
Figure 2. The sensitivity of R 0 to each parameter in the model.
Fractalfract 09 00055 g002
Figure 3. Case 1: Contour plots of the contributions of pairs of parameters to R 0 .
Figure 3. Case 1: Contour plots of the contributions of pairs of parameters to R 0 .
Fractalfract 09 00055 g003
Figure 4. Case 2: Contour plots of the contributions of pairs of parameters to R 0 .
Figure 4. Case 2: Contour plots of the contributions of pairs of parameters to R 0 .
Fractalfract 09 00055 g004
Figure 5. Case 3: Contour plots of the contributions of pairs of parameters to R 0 .
Figure 5. Case 3: Contour plots of the contributions of pairs of parameters to R 0 .
Fractalfract 09 00055 g005
Figure 7. Case 1: Phase-plane behavior for various values of ϕ , where λ v = 0.15 , Λ = 0.01 , β = 0.6 , σ = 0.3 ; μ = 0.015 , δ = 0.01 and γ = 0.02 .
Figure 7. Case 1: Phase-plane behavior for various values of ϕ , where λ v = 0.15 , Λ = 0.01 , β = 0.6 , σ = 0.3 ; μ = 0.015 , δ = 0.01 and γ = 0.02 .
Fractalfract 09 00055 g007
Figure 9. Case 2: Phase-plane behavior for various values of ϕ , where λ v = 0.3 , Λ = 0.01 , β = 0.8 , σ = 0.4 ; μ = 0.02 , δ = 0.03 and γ = 0.01 .
Figure 9. Case 2: Phase-plane behavior for various values of ϕ , where λ v = 0.3 , Λ = 0.01 , β = 0.8 , σ = 0.4 ; μ = 0.02 , δ = 0.03 and γ = 0.01 .
Fractalfract 09 00055 g009
Figure 11. Case 3: Phase-plane behavior for various values of ϕ , where λ v = 0.1 , Λ = 0.01 , β = 0.7 , σ = 0.35 ; μ = 0.01 , δ = 0.005 and γ = 0.04 .
Figure 11. Case 3: Phase-plane behavior for various values of ϕ , where λ v = 0.1 , Λ = 0.01 , β = 0.7 , σ = 0.35 ; μ = 0.01 , δ = 0.005 and γ = 0.04 .
Fractalfract 09 00055 g011
Figure 12. Comparison of the MATLAB ode45 solver with the discrete fractional-order method (25) (for ϕ = 1.0 ).
Figure 12. Comparison of the MATLAB ode45 solver with the discrete fractional-order method (25) (for ϕ = 1.0 ).
Fractalfract 09 00055 g012
Table 1. Sensitivity of R 0 to each parameter.
Table 1. Sensitivity of R 0 to each parameter.
Parameters S λ v S β S σ S μ S δ S γ
Case 1 0.068 1.0 0.923 0.085 0.030 0.039
Case 2 0.250 1.0 0.888 0.127 0.066 0.083
Case 3 0.020 1.0 0.958 0.043 0.013 0.016
Table 2. Parameters for each case study with subcases.
Table 2. Parameters for each case study with subcases.
Parameters λ v Λ β σ μ δ γ
Case 10.150.010.60.30.0150.010.02
Case 20.30.010.80.40.020.030.01
Case 30.10.010.70.350.010.0050.04
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.

Share and Cite

MDPI and ACS Style

Madani, Y.A.; Ali, Z.; Rabih, M.; Alsulami, A.; Eljaneid, N.H.E.; Aldwoah, K.; Muflh, B. Discrete Fractional-Order Modeling of Recurrent Childhood Diseases Using the Caputo Difference Operator. Fractal Fract. 2025, 9, 55. https://doi.org/10.3390/fractalfract9010055

AMA Style

Madani YA, Ali Z, Rabih M, Alsulami A, Eljaneid NHE, Aldwoah K, Muflh B. Discrete Fractional-Order Modeling of Recurrent Childhood Diseases Using the Caputo Difference Operator. Fractal and Fractional. 2025; 9(1):55. https://doi.org/10.3390/fractalfract9010055

Chicago/Turabian Style

Madani, Yasir A., Zeeshan Ali, Mohammed Rabih, Amer Alsulami, Nidal H. E. Eljaneid, Khaled Aldwoah, and Blgys Muflh. 2025. "Discrete Fractional-Order Modeling of Recurrent Childhood Diseases Using the Caputo Difference Operator" Fractal and Fractional 9, no. 1: 55. https://doi.org/10.3390/fractalfract9010055

APA Style

Madani, Y. A., Ali, Z., Rabih, M., Alsulami, A., Eljaneid, N. H. E., Aldwoah, K., & Muflh, B. (2025). Discrete Fractional-Order Modeling of Recurrent Childhood Diseases Using the Caputo Difference Operator. Fractal and Fractional, 9(1), 55. https://doi.org/10.3390/fractalfract9010055

Article Metrics

Back to TopTop