Epidemiological Modeling of Pneumococcal Pneumonia: Insights from ABC Fractal-Fractional Derivatives

1  Introduction

Pneumococcal pneumonia, caused by Streptococcus pneumoniae, remains a significant global health issue, particularly affecting young children and the elderly. Despite medical advancements, it continues to contribute notably to morbidity and mortality worldwide. Traditional mathematical models, including classical and integer-order differential equations, have been employed to study disease dynamics, but these models often fail to account for nonlinear, long-memory effects and hereditary factors that influence disease transmission, immunity, and treatment responses. This lung disease, which can affect individuals of all ages, can manifest with varying severity, with some forms being non-airborne due to inhalation of harmful substances. While most pneumococcal toxins are harmless, certain strains can be extremely dangerous, leading to brain damage and hearing loss. Pneumococcal meningitis is the most severe form of illness, primarily affecting children under five, though it can also impact adults. The bacteria involved in pneumococcal infections can spread through the bloodstream, and while the mortality rate in children under five is around 1%, the elderly have a significantly higher mortality risk, with death rates around 5%. Pneumococcal pneumonia remains one of the leading causes of death in the elderly population.

Pneumococcal pneumonia, a major public health concern, has been extensively modeled using classical and fractional-order differential equations. However, existing models primarily rely on integer-order or Caputo fractional derivatives, which have limitations in capturing complex nonlocal memory effects, hereditary properties, and fractal-like transmission dynamics observed in real-world epidemiological data. These traditional models assume that disease progression depends only on current states, failing to incorporate long-term dependencies and memory effects, which are crucial for understanding infectious disease behavior. To address these limitations, we introduce a novel fractal-fractional Susceptible-Carrier-Infected-Recovered model formulated with the Atangana-Baleanu in Caputo (ABC) sense. Unlike classical or Caputo-based fractional models, the ABC derivative incorporates a non-singular Mittag-Leffler kernel, allowing for a more realistic representation of disease transmission.

A mathematical model of bacteremic Pneumococcal Pneumonia in young children (5-year-old) was developed by Ong’ala et al. [1] to describe the occurrence of this disease. They analyzed the transmission rates and paths between the carriers and the infected class, using the bifurcation theory and the stability of equilibrium points as a means of analyzing the transmission rates. Mochan et al. [2] provided a dynamic ordinary differential equation model of the host immunological response to bacterial Pneumococcal Pneumonia infection in Murine strains. Drusano et al. [3] investigated the efficacy of granulates in preventing bacterial growth and concluded that antibiotics play no part in this process. Ndelwa et al. [4] mathematically represented the dynamic features for the transmission of Pneumococcal Pneumonia, including screening and medicine, and analyzed the transmission and consequences. Kosasih et al. [5] analyzed a mathematical model of cough sounds using wavelet-based crackling detection to rapidly identify bacterial Pneumococcal Pneumonia in young patients. In 2016, César et al. [6] used a mathematical model for pediatric asthma and Pneumococcal Pneumonia over a large population. Based on the work of Marchello et al. [7], it was shown in 2016 that atypical bacterial infections were mostly responsible for the spread of respiratory illnesses such as coughing, bronchitis, and chronic obstructive pulmonary disease (COPD). In 2017, Cheng et al. [8] presented a dynamical mathematical model of the influenza A virus and Streptococcus pneumoniae. Kosasih and Abeyratne [9] provided a simple mathematical model illuminating the clinical analysis of measures for the diagnosis of childhood Pneumococcal Pneumonia and discussed the most common reasons why young infants in low-income areas of the world get Pneumococcal Pneumonia.

In 2018, a co-infection model for Pneumococcal Pneumonia and typhoid was presented by Tilahun et al. [10,11], and their defining connection in the face of a cure and therapeutic approaches was mathematically examined. With the use of mathematical properties of cough sounds, Raj et al. [12] examined the categorization of asthma and Pneumococcal Pneumonia in low-resource communities. Kizito and Tumwiine [13] produced a mathematical model that shows how microorganisms limit the spread of Pneumococcal Pneumonia. Vaccine formulation and treatment dynamics were also examined. Mbabazi et al. [14] looked into a nonlinear mathematical model that describes influenza A virus and Pneumococcal Pneumonia inside the host. A model of Pneumococcal Pneumonia-meningitis coinfection was developed by Tilahun [15], making use of some theorems and ordinary differential equations. Various strategies for eradicating diseases were outlined in detail. Dynamic mathematical models of Pneumococcal Pneumonia were evaluated by Diah and Aziz [16]. Pneumococcal Pneumonia risk was calculated using a different method. In 2019, Tilahun [17] developed a mathematical model of Pneumococcal Pneumonia and bacterial meningitis. A mathematical model of pneumococcal Pneumococcal Pneumonia with temporal delays was developed by Mbabazi et al. [18]. The effectiveness of the model was examined by Otoo et al. [19]. Models of respiratory diseases were presented in graphic form by Zephaniah et al. [20]. A coronavirus Pneumococcal Pneumonia outbreak was reported in Wuhan, China, by Ming et al. [21]. Based on clinical testing, Jung et al. [22] confirmed the findings and identified the pathogen responsible. A mathematical model of Pneumococcal Pneumonia-HIV co-infection was developed by Wafula et al. [23] based on anti-Pneumococcal Pneumonia treatments. Moreover, Oluwatobi and Erinle-Ibrahim, Ref. [24] examined the effects of Pneumococcal Pneumonia while explaining the presence of fundamental and effective reproduction numbers. In Sayed Saber et al. [25], the transmission dynamics of Pneumococcal Pneumonia were investigated using a fractional mathematical susceptible-vaccinated-carrier-infected-recovered model. Similar result on the Pneumonia diseases has been studied in [26].

Fractional-order models, particularly those using Caputo fractional derivatives, have been introduced to incorporate memory-dependent interactions into epidemiological dynamics. Fractal fractional dynamics are observed in various natural systems, including excitation-relaxation systems and natural oscillatory systems. Atangana and Qureshi [27] explore the modeling of chaotic attractors using fractal-fractional operators, focusing on their applicability to dynamical systems in the context of fractional calculus, see also [28–32]. Khalid et al. [33] propose the Modified Minimal Model for diabetes treatment strategies [34] and further investigate the use of fractal-fractional derivatives in diabetes modeling. Almutairi et al. [35] focus on the use of the fractal-fractional Atangana-Baleanu operator in modeling pneumonia, with detailed stability, statistical, and numerical analyses. Abodayeh et al. [36] highlight the use of efficient computational techniques in pneumonia modeling. Boukhouima et al. [37] study a fractional-order HIV infection model with functional response and cure rate, see also [38–40].

The application of fractional calculus in modeling dynamical systems has gained significant attention in recent years due to its ability to capture memory effects and complex dynamics more accurately than classical integer-order models. Several studies have demonstrated the effectiveness of fractional-order systems in various fields, including bioengineering, finance, and control theory [41–45]. For instance, Ref. [41] provides a comprehensive overview of fractional-order modeling and control applications, highlighting their advantages in representing real-world phenomena with hereditary properties. In the context of biological systems, fractional-order models have been particularly useful in describing glucose-insulin interactions. For example, Refs. [42,44] developed fractional-order models to analyze the dynamics of glucose-insulin regulation, demonstrating improved agreement with experimental data compared to traditional models. Similarly, epidemiological models have benefited from fractional calculus, as seen in [43,45], where fractional-order differential equations were used to model the spread of COVID-19, providing insights into transmission dynamics and stability.

Saber and Alahmari [46] discussed the modeling and mathematical analysis of zoonotic diseases, see also [47–49] presented additional advancements in zoonotic disease modeling. Althubyani and Saber [50] also contributed significantly to the understanding of complex dynamical systems. Several studies have demonstrated the effectiveness of fractional-order systems in various fields, including bioengineering, finance, and control theory [51].

The stability analysis of fractional-order systems has also been a key area of research. Theoretical frameworks such as those presented in [52,53] have been extended to fractional-order systems, enabling rigorous stability assessments. Additionally, chaos control in fractional-order chaotic systems has been explored in works such as [54], where the Burke-Shaw system was analyzed using fractal-fractional operators. The study of fractional-order systems is deeply rooted in both classical and modern mathematical frameworks.

Building on these theoretical underpinnings, Petras [55] systematized the modeling, analysis, and simulation of fractional-order nonlinear systems, bridging abstract mathematical theory with practical applications. Fractional calculus has further been applied to chaotic systems, including the Chua system [56], the Lorenz system [57], the Lü system [58], Burke-Shaw [59] and the Rössler system [60]. Recent studies, such as [61], have investigated chaos control in the fractional Newton-Leipnik system using different fractional derivatives. Further applications include disease modeling, where fractional derivatives have been used to study smoking epidemic model [62] and Chikungunya transmission [63] and typhoid fever dynamics [64]. The theoretical foundations of fractal-fractional differential equations have also been explored in [65], providing a mathematical basis for their application in complex systems. Additionally, fractal-fractional operators have been employed to analyze chaotic dynamics, as demonstrated in [66–68].

In this study, we present a fractal-fractional Atangana-Baleanu Caputo (ABC) derivative-order model that describes pneumococcal pneumonia transmission dynamics. The objective is to establish the existence and uniqueness of solutions and analyze the qualitative characteristics of the model. To capture actual disease behavior, we employ the fractal-fractional derivative of the Atangana-Baleanu equation and conduct a local and global stability analysis of steady states. Specifically, the basic reproductive number ℛ0 is derived to demonstrate the stability of the disease-free steady state when ℛ0<1. The existence of a positive (endemic) steady state that is both locally and globally asymptotically stable is ensured if ℛ0>1. To solve the fractional-order system numerically, we utilize the fractal-fractional Atangana-Baleanu numerical method. The proposed system offers a foundational framework for studying, analyzing, and computing solutions to various epidemiological illness models. The numerical simulations validate the theoretical findings, reinforcing the model’s biological relevance.

Addressing the limitations of existing models, this study introduces a fractal-fractional Susceptible-Carrier-Infected-Recovered model using the ABC derivative. The ABC operator, with its non-singular Mittag-Leffler kernel, enhances epidemiological forecast accuracy by ensuring smooth transitions between disease states. Unlike previous models, this approach accounts for nonlocal effects and memory-dependent disease transmission, making it a more effective tool for studying pneumococcal pneumonia outbreaks. We establish the theoretical validity of the model by proving the existence and uniqueness of solutions. We conduct local and global stability assessments, and demonstrate the impact of memory and fractal structures through numerical simulations. Furthermore, we conduct a comprehensive sensitivity analysis using 3D visualizations to identify the most influential parameters affecting disease spread. By integrating a fractal-fractional calculus into epidemiological modeling, this study significantly enhances the predictive capability of disease transmission models, offering a more effective framework for intervention strategies.

To overcome the shortcomings of existing models, this study introduces a fractal-fractional Susceptible-Carrier-Infected-Recovered model using the Atangana-Baleanu in Caputo (ABC) derivative. The ABC operator, which is based on a non-singular Mittag-Leffler kernel, provides a more realistic representation of disease dynamics than previous fractional models. This approach ensures smoother transitions between disease states and enhances epidemiological forecast accuracy, making it a valuable tool for analyzing the long-term behavior of pneumococcal pneumonia outbreaks. We extend traditional pneumococcal pneumonia models by incorporating fractal-fractional operators, which account for nonlocal effects and memory-dependent disease transmission. We establish the existence and uniqueness of solutions, ensuring theoretical validity. We conduct local and global stability assessments to determine the conditions under which the disease persists or is eradicated. Unlike prior models, we employ numerical simulations based on the ABC operator to demonstrate the effects of memory and fractal structure on disease transmission. We conducted a comprehensive sensitivity analysis using 3D visualizations to identify the most influential parameters affecting disease spread. By integrating fractal-fractional calculus into epidemiological modeling, this study enhances the accuracy and predictive capability of disease transmission models, paving the way for more effective intervention strategies.

The novelty of our study lies in the integration of fractal and fractional dynamics. This better reflects real-world epidemic patterns, particularly in cases where transmission pathways are irregular and heterogeneous. This research applies fractal-fractional calculus to generalize traditional epidemiological models. Unlike classical differential equation models, which assume instantaneous effects and local interactions, the proposed approach accounts for long-memory behavior and nonlocal influences. This allows a more accurate representation of pneumococcal pneumonia dynamics, where past infections and immunity play a significant role in shaping future disease spread. Another crucial innovation is the incorporation of nonlinear, memory-dependent interactions into the modeling process. Classical epidemiological models often rely on simplified assumptions about disease transmission, ignoring the cumulative effects of immunity, delayed treatment response, and persistent infections. By integrating memory-dependent dynamics, this study provides a more comprehensive framework for understanding Pneumococcal Pneumonia progression and evaluating potential intervention strategies.

Mathematically, the study establishes the existence and uniqueness of solutions to the proposed model, ensuring its validity and applicability. Additionally, stability properties–including both local and global stability–are rigorously examined to determine the conditions under which disease persistence or eradication occurs. These theoretical foundations are essential for understanding the model’s behavior and implications for real-world disease control. To further illustrate the effectiveness of the proposed model, numerical simulations were conducted to analyze the impact of long-memory effects on disease spread. These simulations provide valuable insights into the transmission dynamics of pneumococcal pneumonia and demonstrate the advantages of fractal-fractional modeling in capturing complex epidemiological patterns. By comparing the results with classical models, the study highlights the significant improvements in accuracy and predictive capability offered by the fractal-fractional approach.

In previous studies on pneumococcal pneumonia dynamics, the limitations have primarily involved the use of classical integer-order models, which fail to capture the memory effects and fractal characteristics inherent in real-world biological systems. These models often assume homogeneous populations and instantaneous changes in disease progression, thereby oversimplifying the complex dynamics of reinfection, immunity decay, and bacterial carriage. This study aims to address these gaps by employing the Atangana-Baleanu Caputo (ABC) fractional operator, which incorporates memory effects and fractional-order parameters to provide a more comprehensive and accurate representation of pneumococcal pneumonia transmission dynamics.

The remainder of the paper is structured as follows: Section 2 provides preliminary definitions and a mathematical background. Section 3 introduces the pneumococcal pneumonia Susceptible-Carrier-Infected-Recovered model formulation with fractal-fractional derivatives. Sections 4 and 5 discuss the qualitative properties of the model, including stability analyses. Section 6 presents a sensitivity analysis of the basic reproduction number ℛ0. Sections 7 and 8 present a numerical scheme and simulations to validate the theoretical findings. Finally, Sections 9 and 10 conclude with a discussion and conclusion on the implications of the results and potential future research directions.

2  Preliminary Definitions

Definition 1([51]). The Gamma function Γ(x) is a continuous extension of the factorial function, defined for real and complex numbers except for non-positive integers. It is given by the following integral representation:

Γ(x)=∫0∞tx−1e−tdt,for Re(x)>0.

Definition 2([51]). The Mittag-Leffler function is a generalization of the exponential function, while the gamma function is a generalization of the factorial function. The Mittag-Leffler function is defined as a power series expansion

Eσ1(z)=∑s=0∞zsΓ(σ1s+1),σ1>0,σ1∈R,z∈C,

where Eσ1 is the single-parameter Mittag-Leffler function.

The two-parameter generalization of the Mittag-Leffler function, which plays an important role in the fractional calculus, is defined by the series expansion

Eσ1,β(z)=∑s=0∞zsΓ(σ1s+β),σ1>0,β>0,σ1,β∈R,z∈C.

When β=1, we denote Eσ1(z)=Eσ1,1(z) to be a one-parameter MittagLeffler function. When both σ1,β are real and positive, the above series converges for values of z, and the Mittag-Leffler function is an entire function.

Definition 3([69,70]). Consider the fractal differentiable on (a,b) of order 0<σ2≤1 for ϕ∈C((a,b),R). The following is a fractal-fractional derivative operator for t in the Atangana-Baleanu setting:

𝒟0,tσ1,σ2FFPϕ(t)=ℏ(σ1)1−σ1ddtσ2∫0tϕ(s)Eσ1[−σ11−σ1(t−s)σ1]ds,

where Eσ1,1(t)=∑k=0∞tkΓ(kq+1)>0 is the Mittag Leffler function with the normalization function ℏ(σ1) is given by: ℏ(σ1)=1−σ1+σ1Γ(σ1), and dh(s)dsσ2=limt→st(t)−t(s)tσ2−ςσ2.

Definition 4([69,70]). The fractal-fractional integration operator are provided by

I0,tσ1,σ2FFPϕ(t)=σ1σ2ℏ(σ1)Γ(σ1)∫0tsσ2−1ϕ(s)(t−s)σ1−1ds+σ2(1−σ1)ϑσ2−1ℏ(σ1)ϕ(t).

The parameter σ1 in the Atangana-Baleanu fractal-fractional (ABC-FF) derivative represents the fractional order of differentiation, while σ2 is associated with the fractal dimension. The distinction is crucial, as σ1 controls the memory effect and smoothness of the function, whereas σ2 determines the complexity of the fractal structure embedded within the model.

When σ1=1, the AB-FF operator reduces to a classical derivative. Specifically:

•   The Mittag-Leffler kernel simplifies to an exponential function, removing nonlocal memory effects.

•   The fractional component vanishes, resulting in an integer-order derivative.

•   The fractal properties disappear, causing the model to behave as a standard differential equation.

However, if σ1=1, the entire formulation must be carefully examined to ensure that singularities do not arise. The authors should explicitly state under what conditions the equation remains well-defined at this limit.

By integrating memory effects and fractal structures, our ABC-based model provides:

•   More accurate predictions of outbreak dynamics, helping policymakers design better intervention strategies.

•   Improved understanding of reinfection risks and long-term immunity, guiding vaccination campaigns.

•   A more realistic assessment of treatment impact over time, aiding in antibiotic stewardship programs.

To help health practitioners grasp the importance of ABC operators, a simple analogy can be introduced:

Traditional models treat disease progression like flipping a switch–either you are infected or you are not. However, pneumococcal pneumonia leaves a lasting imprint on the body, much like a fading memory that still influences decisions. The ABC operator helps us mathematically capture this fading memory, allowing for a more realistic and accurate disease model.

The conditions σ1>0 and σ2≤1 here are imposed to ensure the mathematical validity and the physical meaningfulness of the fractal-fractional derivative operator in the Atangana-Baleanu setting. The reasoning behind these constraints is as follows:

•   σ1>0: The parameter σ1 controls the fractional order of the derivative in the Mittag-Leffler kernel and must be positive to ensure that the fractional derivative is well-defined. A negative or zero value for σ1 would lead to an undefined or singular behavior, particularly in terms of the Mittag-Leffler function, which could break the mathematical formulation and cause instability.

•   σ2≤1: The parameter σ2 defines the order of the fractional derivative with respect to the temporal variable t. The condition σ2≤1 ensures that the fractional order does not exceed unity, which keeps the physical model within realistic bounds. If σ2>1, the operator would involve derivatives of higher order, potentially resulting in a non-physical behavior that doesn’t align with the intended modeling of complex systems (such as diseases or biological processes).

When σ1=1 in the Atangana-Baleanu fractal-fractional (ABC-FF) derivative, the operator reduces to a classical derivative. If σ1=1:

•   The Mittag-Leffler kernel Eσ1 simplifies to an exponential function, removing the nonlocal memory effects of the operator.

•   The AB-FF derivative reduces to the classical derivative, losing its fractional and fractal properties.

•   The model no longer retains memory effects and behaves like a standard first-order differential equation.

•   The fractional-order nature of the derivative disappears, meaning the system no longer retains past information.

•   The fractal nature is lost, causing the system to behave as a regular differential equation without the complexity introduced by fractal geometry.

•   The model transitions back to a standard integer-order system, eliminating the advantages of fractional calculus in capturing long-memory effects and nonlocal interactions.

Setting σ1=1 removes both the fractional and fractal effects, essentially converting the AB-FF derivative into a classical derivative. The model would no longer benefit from nonlocality or memory effects, making it behave like a standard epidemiological model.

3  Formulation of the Fractal-Fractional Model

The Susceptible-Carrier-Infected-Recovered model proposed by Abodayeh et al. in [36] describes the dynamics of a population divided into four compartments: The susceptible, who is at risk of acquiring Pneumococcal Pneumonia infection, is characterized by S(t), C(t): represents individuals who carry Pneumococcal Pneumonia bacteria and may transmit the disease, I(t): indicates who is infected and at risk, R(t): indicates the number of individuals who have recovered from Pneumococcal Pneumonia. The Susceptible-Carrier-Infected-Recovered model assumes a well-mixed population where individuals transition between compartments based on rates of transmission, recovery, and mortality. The model captures the dynamics of both symptomatic and asymptomatic individuals, reflecting the complexity of disease spread in the population. It is important for the authors to explicitly state the biological motivations behind these transitions and to ensure that the parameters reflect real-world processes, especially in the context of disease dynamics and public health interventions.

For any given time t, the parameters and variables describing pneumonia are as follows:

•   S(t): The number of individuals susceptible to pneumonia, at risk of infection.

•   C(t): The number of carriers who harbor the pneumonia bacteria and can spread the infection.

•   I(t): The number of infected individuals capable of transmitting pneumonia.

•   R(t): The number of individuals who have recovered after receiving pneumonia treatment.

•   μ: The natural mortality rate per capita.

•   Π: The recruitment rate of new individuals into the susceptible population per capita.

•   θ: The proportion of susceptible individuals who become carriers.

•   σ: The disease-induced mortality rate per capita.

•   β: The per capita recovery rate of carriers.

•   σ1: The transmission rate of pneumonia to susceptible individuals.

•   τ: The per capita recovery rate of infected individuals.

•   π: The rate at which carriers develop symptoms.

•   η: The rate at which treated individuals return to being susceptible.

•   γ: The vaccination rate of susceptible individuals.

•   ω: The rate of vaccinated treated individuals.

•   δ: The transmission coefficient within the carrier subgroup.

•   p: The probability that a contact leads to infection.

•   k: The contact rate.

•   ℓ: Transmission rate. This is the rate at which individuals become infected (either as carriers or symptomatic) from contact with others in the population.

The governing equations of the model are as follows:

dS(t)dt=Π−ℓS(t)−μS(t)+ηR(t),dC(t)dt=ℓθS(t)−(π+β+μ)C(t),dI(t)dt=ℓ(1−θ)S(t)+πC(t)−(μ+τ+σ)I(t),dR(t)dt=βC(t)+τI(t)−(μ+η)R(t).(1)

Here, ℓ=δ(I(t) + ϖC(t))N∗ is the transmission rate, and the initial conditions are:

S(0)=S0≥0,     C(0)=C0≥0,     I(0)=I0≥0,     R(0)=R0≥0.

Assumptions and Motivations:

•   Susceptible Population (S(t)):

Assumption: Individuals who are susceptible to the disease are denoted by S(t). This group is replenished by a recruitment rate Π, which represents new individuals entering the population (e.g., through birth or immigration).

–   Motivation: In many epidemiological models, the susceptible population increases due to natural birth or migration. The loss of susceptible individuals occurs due to transmission, natural mortality, and movement into the recovered compartment.

Carriers (C(t)):

–   Assumption: Carriers are individuals who are infected but do not exhibit symptoms. These individuals may transmit the disease to others and are important in the spread of the infection.

–   Motivation: The rate at which susceptible individuals become carriers is given by ℓθS(t), where ℓ is the transmission rate and θ is the proportion of susceptible individuals who transition to the carrier state rather than becoming symptomatic. This represents the asymptomatic carriers in the population.

–   Loss of carriers: Carriers can either develop symptoms and transition to the infected group or be removed from the population due to natural mortality, recovery, or vaccination. The parameter π+β+μ represents the combined rate of loss from the carrier compartment.

Infected Population (I(t)):

–   Assumption: Infected individuals are those who exhibit symptoms and contribute to the disease spread. They can transition to recovery or death.

–   Motivation: The term ℓ(1−θ)S(t) represents the rate at which susceptible individuals become symptomatic. The term ℓC(t) accounts for the infection spread from carriers who become symptomatic.

–   Loss of infected individuals: Infected individuals can recover with rate τ, die from the disease with rate σ, or experience natural mortality μ. The total loss rate of infected individuals is μ+τ+σ, which combines these factors.

Recovered Population (R(t)):

–   Assumption: Recovered individuals are those who have overcome the infection, either by recovering naturally or through vaccination. They can no longer transmit the disease.

–   Motivation: The term βC(t) represents the recovery of carriers, while τI(t) represents the recovery of infected individuals. The loss of recovered individuals occurs due to natural mortality with rate μ and the rate at which carriers lose their symptoms (η).

To obtain the model with the fractal-fractional derivative from the classical model of differential equations, you need to replace the standard first-order derivatives with the generalized fractal-fractional derivative 𝒟tσ1,σ2FFP, which generalizes the differentiation process. Here’s how you can proceed: To transform the system into a fractal-fractional system, replace each ordinary derivative ddt with the corresponding fractal-fractional derivative 𝒟tσ1,σ2FFP. Apply the fractal-fractional derivative operator to each equation in the classical model.

𝒟tσ1,σ2FFPS(t)=Π−ℓS(t)−μS(t)+ηR(t),𝒟tσ1,σ2FFPC(t)=ℓθS(t)−(π+β+μ)C(t),𝒟tσ1,σ2FFPI(t)=ℓ(1−θ)S(t)+πC(t)−(μ+τ+σ)I(t),𝒟tσ1,σ2FFPR(t)=βC(t)+τI(t)−(μ+η)R(t).(2)

This is the model with the fractal-fractional derivatives, denoted as model (2). The transition from model (1) to model (2) is made by replacing the usual derivative ddt with the fractal-fractional derivative 𝒟tσ1,σ2FFP in each equation. Finally, the main transformation is the replacement of classical derivatives with the fractal-fractional derivatives in each equation, which introduces memory effects and nonlocality into the system dynamics. The system is now described by equations that account for complex time dependencies rather than relying on instantaneous rates of change, as in the classical case.

The ABC operator has become a crucial tool in fractal-fractional calculus, offering significant advantages for modeling complex systems with memory and hereditary properties, particularly in biological and epidemiological contexts. Its key benefits include incorporating memory effects through a non-singular Mittag-Leffler kernel, which realistically represents historical dependencies in disease dynamics, such as incubation periods and immune responses. The operator’s nonlocality ensures a comprehensive account of past states, surpassing classical derivatives in capturing infection progression and intervention effectiveness over time. Additionally, its flexibility across different orders allows for adaptable modeling of varying complexity levels, while its non-singular kernel enhances numerical stability, preventing singularities in simulations. The Atangana-Baleanu in Caputo (ABC) operator generalizes classical and fractional derivatives, making it a versatile framework for analyzing dynamic systems, with demonstrated success in modeling epidemic dynamics such as Pneumococcal Pneumonia. Traditional integer-order Susceptible-Carrier-Infected-Recovered models assume Markovian processes with instantaneous transitions, which fail to capture pneumococcal pneumonia’s memory effects, heterogeneous transmission, and long-tail recovery distributions. In contrast, the fractional-order Susceptible-Carrier-Infected-Recovered model using the AB operator in the Caputo sense addresses these limitations by incorporating memory effects into disease progression, allowing for non-exponential transmission and recovery rates, and improving data fitting to clinical observations. The Mittag-Leffler kernel further refines this approach by providing a smooth decay of memory effects, avoiding singularities, and ensuring a more accurate representation of pneumococcal pneumonia transmission and control dynamics.

This is the first study that applies the Atangana-Baleanu in Caputo (ABC) operators to model the dynamics of pneumococcal pneumonia. The novelty of this approach appears by incorporating fractal-fractional derivatives with a non-singular Mittag-Leffler kernel, which provides a more realistic representation of disease transmission compared to previous models using integer-order or Caputo derivatives.

4  Positivity, Existence and Uniqueness of the Model Solutions

By R+, we refer to the set of all positive real numbers, Ω={(S,C,I,R)∈R+4:S≥0,C≥0,I≥0,R≥0,  max (|S|,|C|,|I|,|R|)≤N}.

Theorem 1. Model (3.2) has non-negative solutions if and only if θ<1.

Proof. One has

𝒟tσ1,σ2FFPS(t)|S=0=Π+δR(t)>0,𝒟tσ1,σ2FFPC(t)|C=0=δθI(t)S(t)N+δϖθC(t)S(t)N>0,𝒟tσ1,σ2FFPI(t)|I=0=δϖC(t)N(1−θ)S(t)+πC(t)>0,𝒟tσ1,σ2FFPR(t)|R=0=βC(t)+τI(t)>0.

To understand why Π+δR(t)>0 here, we need to consider the nature of the terms involved: If Π represents a positive rate, such as a birth or recruitment rate, it is positive by definition. Since R(t)>0 by assumption, and δ>0, the term δR(t) is also positive. Thus, the sum Π+δR(t) is positive because it is the sum of two positive quantities: Π and δR(t). As a result, we conclude that:

Π+δR(t)>0.

Based on Lemmas 5 and 6 in [37], all of the solutions to the model (2) have a semi-positive outcome. ▀

Model (2) can be stated as because the integration is differentiable.

𝒟0,tσ1,σ2ABRS(t)=σ2tσ2−1F1(t,S,C,I,R),𝒟0,tσ1,σ2ABRC(t)=σ2tσ2−1F2(t,S,C,I,R),𝒟0,tσ1,σ2ABRI(t)=σ2tσ2−1F3(t,S,C,I,R),𝒟0,tσ1,σ2ABRR(t)=σ2tσ2−1F4(t,S,C,I,R).

with F1,F2,F3,F4 defined as:

F1=Π−ℓS(t)−μS(t)+ηR(t),F2=ℓθS(t)−(π+β+μ)C(t),F3=ℓ(1−θ)S(t)+πC(t)−(μ+τ+σ)I(t),F4=βC(t)+τI(t)−(μ+η)R(t).(3)

Consider

𝒟0,tσ1,σ2ABRΘ(t)=σ2tσ2−1⨿(t,Θ(t)),

with Θ(0)=Θ0. With Θ(t)=(S(t),C(t),I(t)),Θ(0)=(S(0),C(0),I(0)), and using fractional integral and 𝒟0,tσ1,σ2ABR in place of 𝒟0,tσ1,σ2ABR, we have

Θ(t)=Θ(0)+σ2tσ2−1(1−σ1)k(σ1)⨿(t,Θ(t))+σ1σ2Γ(σ1)k(σ1)×∫0tεσ2−1(t−ε)σ2−1⨿(ε,Θ(ε))dε.

Banach space J=C([0, K]×R3,R) is used to find existence and uniqueness, where J=C[0,K] having

‖Θ‖=maxt∈[0,K][|S(t)|+|C(t)|+|I(t|)].

Establishing an operator Ψ:J→J, as

Ψ(Θ)(t)=Θ(0)+σ2tσ2−1(1−σ1)k(σ1)⨿(t,Θ(t))+σ1σ2Γ(σ1)k(σ1)∫0tεσ2−1(t−ε)σ2−1×⨿(ε,Θ(ε))dε.

Assuming a nonlinear function ⨿(t,Θ(t)) meets the Lipchitz function and growth, then

(H1) For all Θ∈J, there are positive constants 𝒞φ>0 and 𝒮φ satisfies

|⨿(t,Θ(t))|≤𝒞φ|Θ(t)|+𝒮φ.

(H2) For all Θ,Θ¯∈J, there is a positive constant Tφ>0 satisfies

|⨿(t,Θ(t))−⨿(t,Θ¯(t))|≤Tφ|Θ(t)−Θ¯(t)|.

Theorem 2. Assume (H1) and (H2) are true. Let ⨿:[0,K]×J→J be a continuous function. As a result, there is only one solution to the model.293

Proof. There is no difference between ⨿ and Ψ since ⨿ is permanently fixed. Assume the following:

Φ={Θ∈J:‖Θ‖≤L,L>0}.

For any Θ∈J, then we get

|Ψ(Θ)|=maxt∈[0,K]|Θ(0)+vt(1−σ1)k(σ1)×⨿(t,Θ(t))σ1σ2Γ(σ1)k(σ1)×∫0tεσ2−1(t−ε)σ2−1⨿(ε,Θ(ε))dε|≤Θ(0)+vKσ2−1(1−σ1)k(σ1)(𝒞φ‖Θ‖+𝒮φ)+σ1σ2Γ(σ1)k(σ1)(𝒞φ‖Θ‖+𝒮φ)×Kσ1+v−1Φ(σ1,v)≤L.

To ensure that the homogeneity of the function is guaranteed, the operator Ψ must have an input function Φ(σ1,v). Assume t1,t2≤K for Ψ equicontinuity. Take into consideration the following points:

|Ψ(Θ)(t2)−Ψ(Θ)(t1)|=|σ2t2σ2−1(1−σ1)k(σ1)⨿(t2,Θ(t2))+σ1σ2Γ(σ1)k(σ1)∫0t2εσ2−1(t2−ε)σ2−1⨿(ε,Θ(ε))dε−σ2t1σ2−1(1−σ1)k(σ1)⨿(t1,Θ(t1))σ1σ2Γ(σ1)k(σ1)∫0t1εσ2−1(t1−ε)σ2−1⨿(ε,Θ(ε))dε|−σ2t2σ2−1(1−σ1)k(σ1)(𝒞φ‖Θ‖+𝒮φ)+σ1σ2Γ(σ1)k(σ1)(𝒞φ‖Θ‖+𝒮φ)t2σ1+σ2−1Φ(σ1,σ2)−σ2t1σ2−1(1−σ1)k(σ1)(𝒞φ‖Θ‖+𝒮φ)+σ1σ2Γ(σ1)k(σ1)(𝒞φ‖Θ‖+𝒮φ)t1σ1+σ2−1Φ(σ1,σ2).

Consequently, as t1→t2, ‖Ψ(Θ)(t2)−Ψ(Θ)(t1)‖→0. Since Ψ is continuous, the theory of Arzela-Ascoli is fully continuous. ▀

Theorem 3. Suppose that (H1) and (H2) is true. If Ξ<1 where

Ξ=(vkσ2−1(1−σ1)k(σ1)+σ1σ2Γ(σ1)k(σ1)Kσ1+v−1Φ(σ1,v))Tφ,

Then the solution is unique, if it exists.

Proof. For Θ,Θ¯∈J, we have

|Ψ(Θ)−Ψ(Θ¯)|=maxt∈[0,K]|σ2tσ2−1(1−σ1)k(σ1)[⨿(t,Θ(t))−⨿(t,Θ(t))]+σ1σ2Γ(σ1)k(σ1)∫0tεσ2−1(t−ε)σ2−1dε×(⨿(ε,Θ(ε))−⨿(ε,Θ¯(ε)))|≤[vKσ2−1(1−σ1)k(σ1)+σ1σ2Γ(σ1)k(σ1)Kσ1+v−1Φ(σ1,v)]‖Θ−Θ¯‖≤Ξ‖Θ−Θ¯‖.

5  Stability Analysis

The total population at any time t is defined as:

N(t)=S(t)+C(t)+I(t)+R(t).

To derive the equation for N(t), we sum the given differential equations:

dN(t)dt=dS(t)dt+dC(t)dt+dI(t)dt+dR(t)dt.

Using the given system of equations, we have:

dN(t)dt=(Π−ℓS(t)−μS(t)+ηR(t))+(ℓθS(t)−(π+β+μ)C(t))+(ℓ(1−θ)S(t)+ℓC(t)−(μ+τ+σ)I(t))+(βC(t)+τI(t)−(μ+η)R(t)).

The resulting equation is:

dN(t)dt=Π−μN(t)−σI(t).

At equilibrium (dN(t)dt=0), we have:

0=Π−μN∗−σI∗.

Rearranging the equation gives:

N∗=Π−σI∗μ.

The total population at equilibrium, N∗, depends on the recruitment rate Π, the natural death rate μ, and the disease-induced death rate σI∗. The term σI∗ represents the loss of population due to infection-induced deaths.

When I=C=0, one obtains the following equation:

0=Π−μS(t)+ηR(t),0=−(μ+η)R(t).

Implies E1=(Πμ,0,0,0). At E1, as in by [71], the next-generation matrix theory allows us to deduce the reprduction number ℛ0 as:

ℛ0=(ϖθ(μ+τ+σ)+(1−θ)((πϖ+(π+β+μ)))δS0N∗(π+β+μ)(μ+τ+σ).

Lemma 1 ([36]). The Pneumococcal Pneumonia infection-free equilibrium E1 is locally asymptotically stable in Ω if ℛ0<1, but it is unstable if ℛ0>1.

Theorem 4. The disease-free equilibrium E1 of the fractal-fractional system is globally asymptotically stable in Ω.

Proof. Consider the Lyapunov function:

L1(S,C,I,R)=(S−S0−S0ln⁡SS0).

From [53], applying the fractal-fractional derivative 𝒟tσ1,σ2FFP, we have:

𝒟tσ1,σ2FFPL1(S,C,I,R)≤(S−S0S)𝒟tσ1,σ2FFPS.

Substituting the fractal-fractional system equations, we get:

𝒟tσ1,σ2FFPL1(S,C,I,R)=(S−S0S)(Π−ℓS−μS+ηR).

At E1, we have:

𝒟tσ1,σ2FFPL1(S,C,I,R)≤(S−S0)(ΠS−ℓ−μ+ηRS)=(S−S0)(ΠS+ηRS−ΠS0−ηRS0)=(S−S0)(−ΠSS0(S−S0)−ηRSS0(S−S0))=−ΠSS0(S−S0)2−ηRSS0(S−S0)2.

Since this expression is negative for all (S,C,I,R)∈Ω, it follows from LaSalle’s invariance principle [52] that E1 is globally asymptotically stable in Ω.

For I∗>0, and taking

{𝒟tσ1,σ2FFPS(t)=0,     𝒟tσ1,σ2FFPC(t)=0,     𝒟tσ1,σ2FFPI(t)=0,     𝒟tσ1,σ2FFPR(t)=0.

With α1=θπ+(1−θ)(π+β+μ)(μ+τ+σ), one obtains the endemic equilibrium E2=(S∗,C∗,I∗,R∗), with

S∗=N∗(π+β+μ)α1+ϖ,C∗=(μ+η)(Π(α1+ϖ)−μN∗(π+β+μ))(α1+ϖ)(δ(π+β+μ)(μ+η)−η(σ2+τα1)),I∗=α1(μ+η)(Π(α1+ϖ)−μN∗(π+β+μ))(α1+ϖ)(δ(π+β+μ)(μ+η)−η(σ2+τα1)),R∗=(σ2+τα1)(Π(α1+ϖ)−μN∗(π+β+μ))(α1+ϖ)(δ(π+β+μ)(μ+η)−η(σ2+τα1)).

Lemma 2 ([36]). For ℛ0>1, a unique endemic equilibrium point E2 exists; otherwise, no endemic equilibrium exists.

Lemma 3 ([36]). If R0>1, E2 is stable locally asymptotically in Ω.

Theorem 5. The endemic equilibrium E2 of the fractal-fractional system is globally asymptotically stable in Ω.

Proof. Consider the Lyapunov function:

L2(S,C,I,R)=(S−S∗−S∗ln⁡SS∗)+(C−C∗−C∗ln⁡CC∗)+(I−I∗−I∗ln⁡II∗)+(R−R∗−R∗ln⁡RR∗).

From [53], applying the fractal-fractional derivative 𝒟tσ1,σ2FFP, we have:

𝒟tσ1,σ2FFPL2(S,C,I,R)≤(S−S∗S)𝒟tσ1,σ2FFPS+(C−C∗C)𝒟tσ1,σ2FFPC+(I−I∗I)𝒟tσ1,σ2FFPI+(R−R∗R)𝒟tσ1,σ2FFPR.

Substituting the fractal-fractional system equations, we get:

𝒟tσ1,σ2FFPL2(S,C,I,R)=−ΠN∗SS∗(S−S∗)2−ηRN∗SS∗(S−S∗)2−ℓθSN∗CC∗(C−C∗)2−ℓ(1−θ)SCN∗II∗(I−I∗)2−πCII∗(I−I∗)2−βCRR∗(R−R∗)2−τIRR∗(R−R∗)2.

Since this expression is negative for all (S,C,I,R)∈Ω, it follows from LaSalle’s invariance principle [52] that E2 is globally asymptotically stable in Ω.

6  Sensitivity Analysis of the Basic Reproduction Number ℛ0

The sensitivity index of ℛ0 with respect to a given parameter σ is defined as:

Γσℛ0=∂ℛ0∂σ×σℛ0.

The partial derivatives of ℛ0 with respect to the key parameters are given below.

•   For ϖ:

∂ℛ0∂ϖ=θ(μ+τ+σ)+(1−θ)πδS0N(μ+τ+σ)(μ+β+π).

Γϖℛ0=∂ℛ0∂ϖ×ϖℛ0.

•   For θ:

∂ℛ0∂θ=ϖ(μ+τ+σ)−πϖδS0N(μ+τ+σ)(μ+β+π).

Γθℛ0=∂ℛ0∂θ×θℛ0.

•   For τ:

∂ℛ0∂τ=ϖθN(μ+τ+σ)2(μ+β+π).

Γτℛ0=∂ℛ0∂τ×τℛ0.

•   For μ:

∂ℛ0∂μ=−ϖθ(μ+τ+σ)+(1−θ)(πϖ+δS0)N(μ+τ+σ)2(μ+β+π).

Γμℛ0=∂ℛ0∂μ×μℛ0.

•   For σ:

∂ℛ0∂σ=−ϖθN(μ+τ+σ)2(μ+β+π).

Γσℛ0=∂ℛ0∂σ×σℛ0.

•   For π:

∂ℛ0∂π=−(1−θ)ϖδS0N(μ+τ+σ)(μ+β+π)2.

Γπℛ0=∂ℛ0∂π×πℛ0.

•   For δ:

∂ℛ0∂δ=(1−θ)(πϖ+δS0)N(μ+τ+σ)(μ+β+π).

Γδℛ0=∂ℛ0∂δ×δℛ0.

The Basic Reproduction Number ℛ0 for the given model is computed for the Disease-Free Equilibrium (DFE) and Endemic Equilibrium (EE) as follows:

ℛ0=0.0547     (DFE),     ℛ0=0.0684     (EE).

The sensitivity indices for the parameters are computed as follows:

Γϖℛ0=0.00018,     Γθℛ0=−1.2882,     Γμℛ0=−0.0829,     Γτℛ0=1.539×10−5

Γσℛ0=−1.273×10−5,     Γπℛ0=−0.0128,     Γδℛ0=0.5542.

These sensitivity indices indicate the effect of parameter changes on the basic reproduction number ℛ0. For example, the parameter θ has a strong negative impact on ℛ0, while δ has a positive impact, indicating that increasing δ increases the disease transmission potential.

The sensitivity index quantifies how changes in each parameter of the model affect the basic reproduction number ℛ0, which represents the average number of secondary infections produced by a single infected individual in a completely susceptible population. A positive sensitivity index indicates that an increase in the corresponding parameter leads to an increase in ℛ0, while a negative value means that an increase in the parameter leads to a decrease in ℛ0.

Fig. 1 presents the sensitivity index analysis. The objective is to help readers quickly ascertain the parameters that significantly influence the spread of the disease in the population. By analyzing the sensitivity of model parameters, we can determine which factors have the most substantial impact on disease transmission dynamics. Parameters such as the transmission rate, recovery rate, and vaccination rate exhibit the highest sensitivity, indicating their crucial role in controlling disease spread. These sensitivity analyses offer valuable insights for public health interventions by identifying which parameters should be targeted for effective disease mitigation strategies.

Figure 1: Sensitivity index analysis using 3D visualizations

7  Numerical Scheme of the Fractal-Fractional Pneumonia Model

Atangana-Baleanu fractal-fractional operators are implemented via Lagrangian piecewise interpolation for the proposed model. As in [72], consider system (3.2) in this case

𝒟0,tσ1,σ2FF-ABΨ(t)=Π(t,Ψ(t)),

The Antangana-Baleanu integral gives us

ϑ(t)=Ψ(0)+1−σ1ℏ(σ1)Π(t,Ψ(t))+σ1ℏ(σ1)Γ(σ1)∫0t(t−ξ)σ1−1ξσ2−1Π(ξ,Ψ(ξ))dξ.

Replacing t by tn+1 we have

Ψn+1=Ψ(0)+1−σ1ℏ(σ1)Π(tn+1,Ψ(t))+σ1ℏ(σ1)Γ(σ1)∫0tn+1(tn+1−ξ)σ1−1ξσ2−1Π(ξ,Ψ(ξ))dξ.

Application of the two-step Lagrange polynomial yields

Π(t,(y,Ψ(t))=(y−tξ−1)Π(t,(tξ,Ψ(tξ))tξ−tξ−1−(y−tξ)Π(tξ−1,Ψ(tξ−1)tξ−tξ−1=Π(t,(tξ,Ψ(tξ))(x−tξ−1)tξ−tξ−1−Π(tξ−1,Ψ(tξ−1)(y−tξ)tξ−tξ−1=Π(t,(tξ,Ψ(tξ))(y−tξ−1)h−Π(tξ−1,Ψ(tξ−1)(y−tξ)h.