The fractal-fractional Atangana-Baleanu operator for pneumonia disease: stability, statistical and numerical analyses

1.   Introduction

Infections of the lungs, such as pneumonia, have a variety of causes. The prevalence of this disease is increasing in all age groups and is a major medical concern. Several researchers are working on mathematical models that describe disease spread and optimal control problems in epidemics because they are highly interesting. As a result, these models play a critical role in predicting the effects of epidemics and diseases on areas and populations, as well as the environment. Researchers have presented models for modeling pneumonia dynamics based on a review of the literature; see, e.g., ^{[1,2,3,4,5,6,7,8]}. Based on a mathematical analysis of pneumonia and typhoid characteristics, Tilahun et al. ^[6] proposed a coinfection model. Tilahun et al. ^[7] used ordinary differential equations and a few theorems to model pneumonia and meningitis coinfections in 2018.

Since 1970, infectious disease dynamics has emerged as an interdisciplinary field. Epidemiology studies disease spread. Modeling diseases and their effects on humans is described in ^[5]. Fractional and fractal calculus are combined here. In engineering, physics, biology, and biomedicine, fractal-fractional operators are widely used to model real-world processes. Comparable to classical models, fractional order integrals and fractional derivatives are more precise than classical models. In fractional derivatives, there are three different types of operators: Riemann-Liouville and Caputo, Caputo-Fabrizio and Atangana-Baleanu, which are connected to power laws, exponential decay laws, and extended Mittag-Leffler functions ^{[9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52]}.

The concept of fractal-fractional order integration and differentiation was developed by Atangana with two orders, i.e., one fractional and the other fractal ^[46,54]. Besides, fractal differentiation is equivalent to classical differentiation if the fractal order tends to 1. Fractal behaviors are investigated through the use of these combined operators. Several researchers have shown that fractal-fractional operators better capture real-world mathematics ^{[55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79]}. These include but are not limited to, for instance, the HIV/AIDS model ^[57], Leishmania model ^[58], tuberculosis model ^[59], Q fever model ^[60], hepatitis C virus model ^[61], AH1N1/09 virus model ^[62] and tobacco smoking model ^[63].

This study investigates the formulation of the Tilahun et al. mathematical model of pneumonia transmission dynamics by using fractal-fractional derivatives in the Atangana-Baleanu sense. The model is improved by assuming five unknowns and using this Caputo, and Atangana-Baleanu-type fractional derivatives. The study aims to investigate and compare the solutions to this system, which unique as compared to other studies. The authors construct schemes for this system, using the fractal-fractional Atangana-Baleanu operator to prove the existence, uniqueness, non-negativity and boundedness of solutions. Three levels of stability are established: local, global and Hyers-Ulam. Sensitivity analysis is conducted to assess the impact of parameters on initial disease transmission. The study finds that $\mathscr{R}_{0}$ is most sensitive to pneumonia disease density. The results show that the scheme's method is effective and suitable for the system defined by Caputo and Atangana-Baleanu fractional derivatives. Simulations using Matlab show that the schemes method is suitable for both types of problems and has approximate solutions that are close to the exact solution. The study also discusses other fractional operators and numerically verifies their mathematical findings for the proposed model's dynamical behavior.

2.   Preliminaries

2.1.   Fractal-fractional pneumonia model

Definition 1 (^[76,77]). Consider the fractal to be differentiable on $(a, b)$ of order $0 < \tau_{2} \leq 1$ for $\phi \in C((a, b), \mathbb{R})$. The following is a fractal-fractional derivative operator for $\mathtt{t}$ in the Atangana-Baleanu setting:

where, $\hbar(\tau_{1}) = 1-\tau_{1}+\frac{\tau_{1}}{\Gamma(\tau_{1})}$, and $\frac{dh(s)}{ds^{\tau_{2}}} = \lim _{t \rightarrow s} \frac{t(\mathtt{t})-t(s)}{\mathtt{t}^{\tau_{2}}-\varsigma^{\tau_{2}}}$.

Definition 2 (^[76,77]). The fractal-fractional integration operator is given by

There are five populations in the pneumonia model: susceptible $(\operatorname{x})$, vaccinated $(\operatorname{y})$, carrier $(\operatorname{z})$, infected $(\operatorname{u})$ and recovered $(\operatorname{v})$. The total human population, denoted by \(N\), can be expressed as \(N = \operatorname{x} + \operatorname{y} + \operatorname{z} + \operatorname{u} + \operatorname{v}\). Therefore, our suggested fractal-fractional pneumonia model in the sense of the Atangana-Baleanu derivative looks like this:

subject to $\operatorname{x}(0)\geq 0, \; \operatorname{y}(0)\geq 0, \; \operatorname{z}(0)\geq 0, \; \operatorname{u}(0)\geq 0, \; \operatorname{v}(0)\geq 0.$

All of the positive parameters are listed in Table 1. When people get infected, they either join the carrier class $\operatorname{z}$ or the infectious class $\operatorname{u}$ based on a probability of $1-\varrho$. Let $\varUpsilon$ be the transmission coefficient for the carrier. Infection force is defined as $\lambda = \operatorname{a} \operatorname{z}+ \operatorname{b} \operatorname{u}$, where $a = \frac{k\tau \varUpsilon}{N}$ represents the carrier compartment transmission and $b = \frac{k\tau}{N}$ represents the infective compartment transmission. The population is $\operatorname{N} = \operatorname{x}(\mathtt{t})+\operatorname{y}(\mathtt{t})+\operatorname{z}(\mathtt{t})+\operatorname{u}(\mathtt{t})+\operatorname{v}(\mathtt{t})$. If $\varUpsilon > 1$, carriers are more likely to infect susceptibles than infectious individuals. Both carriers and infectives have the same chance of spreading when $\varUpsilon = 1$. Nevertheless, if $\varUpsilon < 1$, the infective has a higher chance of contacting the susceptible.

2.2.   Existence and uniqueness

The matrix form of (1.2) is given by:

Define the Banach space $\mathbb{U} = X^5$, where $X = C(I, \mathbb{R})$ is subject to the norm

Let $\|\operatorname{x} \| \leq \eta_{1}$, $\|\operatorname{y} \| \leq \eta_{2}$, $\|\operatorname{z} \| \leq \eta_{3}$, $\|\operatorname{u} \| \leq \eta_{4}$ and $\|\operatorname{v} \| \leq \eta_{5}$ for some constants $\eta_{1}, \eta_{2}, \eta_{3}, \eta_{4}, \eta_{5} > 0$.

Rewrite (2.1) as follows

where

Applying fractional integrals, we get

$\varLambda(\mathtt{t}, \Psi(\mathtt{t}))$ must satisfy these Lipschitz and growth conditions.

Theorem 1. For each $\Psi_{1}, \Psi_{2} \in \mathscr{B}$, $\exists$ a constant $\operatorname{A} > 0$ that satisfies

where $\operatorname{A} = \text{max}\left\{\omega_1, \omega_2, \omega_3, \omega_4, \omega_5\right\},$ with $\omega_1 = \mu+2\lambda +2\lambda \varrho+2\vartheta$, $\omega_2 = \mu+2\phi+2\varepsilon\lambda +2\lambda \varepsilon\varrho$, $\omega_3 = 2\chi+\mu+2\beta, \omega_4 = 2\mathtt{q}\eta+\mu+\alpha+2\eta, \omega_5 = 2\delta+\mu.$

Proof. For each $\Psi_{1}, \Psi_{2} \in \mathscr{B}$, one obtains

Therefore, $\varLambda(\mathtt{t}, \Psi(\mathtt{t})$ satisfies the Lipschitz condition. □

Theorem 2. There are constants $\boldsymbol{\operatorname{z}}_{\Psi} > 0$ and $\mathbb{M}_{\Psi}$ satisfies the following, for each $\Psi$ in $\mathscr{B}$,

So, there is at least one solution to the suggested model.

Proof. To begin with, we demonstrate that the operator $\varLambda$ stated in (2.2) is totally continuous. Due to the continuous nature of $\Psi$, $N$ is also continuous. □

Theorem 3. Assume that (2.4) is true; then,

So, it has a unique solution.

Proof. For $\Psi_{1}, \Psi_{2}$ in $\mathscr{B}$, we obtain

Due to this, $\mathscr{H}$ is a contraction and there is only one solution to the model. □

The non-negativity and boundedness of the solutions of the system (2.1) in the fractional case have been proved in ^[8].

3.   Model analysis

An infection-free equilibrium is $E^0 = (\operatorname{x}_0, \operatorname{y}_0, 0, 0, 0) = \Big(\frac{\pi\aleph_{1}}{\mu}, \frac{\pi\aleph_{2}}{\mu}, 0, 0, 0\Big)$, with $\aleph_{1} = \frac{\mu+\phi-p\mu}{\mu+\phi+\vartheta}$ and $\aleph_{2} = \frac{\vartheta+p\mu}{\mu+\phi+\vartheta}$ for $\operatorname{z} = \operatorname{u} = \operatorname{v} = 0$. The endemic equilibrium, given as $E^\ast = (\operatorname{x}^\ast, \operatorname{y}^\ast, \operatorname{z}^\ast, \operatorname{u}^\ast, \operatorname{v}^\ast)$ for $\operatorname{u} > 0, \operatorname{z} > 0$ was acquired by applying the following:

Thus, $E^\ast = (\operatorname{x}^\ast, \operatorname{y}^\ast, \operatorname{z}^\ast, \operatorname{u}^\ast, \operatorname{v}^\ast) = (\operatorname{E}_5+\operatorname{E}_6\operatorname{u}^\ast, \operatorname{E}_3+\operatorname{E}_4\operatorname{E}_5+\operatorname{E}_4\operatorname{E}_6\operatorname{u}^\ast, \operatorname{E}_1\operatorname{u}^\ast, \operatorname{u}^\ast, \operatorname{E}_2\operatorname{u}^\ast),$ with

The basic reproduction number is given by ^[7,8],

3.1.   Stability of an infection-free equilibrium $E_0$

Lemma 1. If $\mathscr{R}_{0} = 1$, $E_0$ is locally asymptotically stable for model (2.1), and unstable if $\mathscr{R}_{0} > 1$. Moreover, model (2.1) has a globally asymptotically stable $E_0$.

Proof. The first part follows as in ^[8]. We present a positive definite Lyapunov function:

One obtains

At $E_0$, one obtains

Thus, ${ }^\mathtt{FF-AB} \mathscr{D}_{0, \mathtt{t}}^{\tau_{1}, \tau_{2}} L_1 < 0$ for all $(\operatorname{x}, \operatorname{y}, \operatorname{z}, \operatorname{u}, \operatorname{v})\in \Lambda$. Moreover, ${ }^\mathtt{FF-AB} \mathscr{D}_{0, \mathtt{t}}^{\tau_{1}, \tau_{2}} L_1 = 0$ implies that $\operatorname{x} = \operatorname{x}_0$, $\operatorname{y} = \operatorname{y}_0$, $\operatorname{z} = \operatorname{z}_0$, $\operatorname{u} = \operatorname{u}_0$ and $\operatorname{v} = \operatorname{v}_0$. So, $\{E_0\}$ is the only set satisfying that ${ }^\mathtt{FF-AB} \mathscr{D}_{0, \mathtt{t}}^{\tau_{1}, \tau_{2}} L_1 = 0$. □

3.2.   Stability of an endemic equilibrium $E^\ast$

Lemma 2. $E^\ast$ exists when $\mathscr{R}_{0} > 1$; otherwise, there is no endemic equilibrium.

Proof. The following characteristics are required for a disease to be endemic: ${ }^\mathtt{FF-AB} \mathscr{D}_{0, \mathtt{t}}^{\tau_{1}, \tau_{2}} \operatorname{z}(\mathtt{t}) > 0$ and ${ }^\mathtt{FF-AB} \mathscr{D}_{0, \mathtt{t}}^{\tau_{1}, \tau_{2}} \operatorname{u}(\mathtt{t}) > 0$, that is,

Given (3.1), based on the first inequality,

Then,

Because $\frac{(\operatorname{x}(\mathtt{t})+\varepsilon \operatorname{y}(\mathtt{t}))}{N} \leqslant1$, one obtains

As a result of the second inequality of (3.1),

Then,

Using the fact that $\frac{(\operatorname{x}(\mathtt{t})+\varepsilon \operatorname{y}(\mathtt{t}))}{N} \leqslant1$, one obtains

Substituting (3.2) into (3.3), one obtains

After rearranging and canceling $\operatorname{u}(\mathtt{t})$, one gets

Thus, $\mathscr{R}_{0} > 1$ creates a unique endemic equilibrium. □

Lemma 3 (^[8]). $E^\ast$ is locally asymptotically stable for $\mathscr{R}_{0} > 1$. Moreover, $E^\ast$ is globally asymptotically stable.

Proof. Define

Thus, one obtains

Thus, ${ }^\mathtt{FF-AB} \mathscr{D}_{0, \mathtt{t}}^{\tau_{1}, \tau_{2}} L_2 < 0$ for all $(\operatorname{x}, \operatorname{y}, \operatorname{z}, \operatorname{u}, \operatorname{v})\in \Lambda$. Furthermore, ${ }^\mathtt{FF-AB} \mathscr{D}_{0, \mathtt{t}}^{\tau_{1}, \tau_{2}} L_2 = 0$ implies that $\operatorname{x} = \operatorname{x}^*$, $\operatorname{y} = \operatorname{y}^*$, $\operatorname{z} = \operatorname{z}^*$, $\operatorname{u} = \operatorname{u}^*$, and $\operatorname{v} = \operatorname{v}^*$. Therefore, according to Theorem 5, $E^*$ is globally asymptotically stable. □

3.3.   Hyers-Ulam stability

The Hyers-Ulam stability has been motivated by the work done in ^[80,81].

Definition 3. The constants $\zeta_i > 0$, for $i \in \mathbb{N}_{1}^{5}$ must meet the following conditions for every $\zeta_i > 0$, $i \in \mathbb{N}_{1}^{5}$, for model (2.1) to have Hyers-Ulam stability:

In the model (2.1), an approximation is $\Big(\operatorname{x}_1(\mathtt{t}), \operatorname{y}_1(\mathtt{t}), \operatorname{z}_1(\mathtt{t}), \operatorname{u}_1(\mathtt{t}), \operatorname{v}_1(\mathtt{t})\Big)$, which satisfies the following:

so that

Theorem 4. If (3.1) is true, then model (2.1) has Hyers-Ulam stability.

Proof.

Then,

Similarly, one obtains

Hence, the results follows. □

3.4.   Sensitivity analysis

According to the parameters of our model, the following equation yields the sensitivity index of $\mathscr{R}_0$:

where $\omega$ is a value from Table 1. Table 1 lists the sensitivity indices of $\mathscr{R}_0$. It is easy to verify that

The sensitivity index of each parameter in the model obtained as in Table 1 by applying (3.5). The sensitivity indexes of Table 1 indicate that $\mathscr{R}_0$ increases as the parameters $\chi$, $\mu$, $a$, $b$ and $\epsilon$ are increased. In contrast, the values of other parameters are fixed. Based on these indices, it appears that disease endemicity has increased. In contrast, when the parameters $\beta$, $\eta$, $q$, $\alpha$ and $\varrho$ are decreased while the rest of the parameters are maintained, $\mathscr{R}_0$ decreases.

4.   Numerical procedures

Atangana-Baleanu fractal-fractional operators are implemented via Lagrangian piecewise interpolation for the proposed model.

4.1.   Adams-Bashforth-Moulton method application

As in ^[82,83], consider system (2.2) in the following case:

subject to ceil function $n = [\tau_{1}]$ and for $t\in [0, T]$, $0 < \tau_{1}\leq1$ with ${ }^{\operatorname{z}}\mathscr{D}_{0, \mathtt{t}}^{\tau_{1}} \Psi(0) = \Psi^{(\kappa)}_0$, $\kappa = 0, 1, 2, ..., n-1$. Volterra's integral equation of system (2.3) is given by

It is easy to reconstruct Eq (4.1) by using the product rule for rectangles,

where $\mathbb{A}_{\kappa, n+1}$ is given by

Let $\{t_n = nh: n = -k, -k +1, ..., -1, 0, 1, ..., N\}$, with $h = T/N$. Then, (4.1) can be discretized as follows:

The predicted value $\Psi^p_h(t_{n+1})$ is determined as follows:

where

According to (4.2), (2.1) is as follows:

where

4.2.   Application of the fractal-fractional Atangana-Baleanu scheme

The Antangana-Baleanu integral gives us

Replacing $\mathtt{t}$ with $\mathtt{t}_{n+1}$ we have

Application of the two-step Lagrange polynomial yields

By using the Lagrange polynomial to solve the given problem, we obtain

Now, solving the integral we get

Now, replacing the value of $\varLambda(y, \Psi(\mathtt{t}))$, we get

As a result, the numerical scheme above rewritten as follows:

4.3.   Numerical simulations and discussion

The above analyses are displayed in Figures 2–16, which display the time series of the model (2.1) under the following initial conditions: $\operatorname{x}(0) = 8200, \operatorname{y}(0) = 2800, \operatorname{z}(0) = 200, \operatorname{u}(0) = 210, \operatorname{v}(0) = 200$. Figures 2–6 show that the time series of the model (2.1) with different trajectories of infected states tends to zero whenever $\operatorname{R}_0 = 0.0083 < 1$. The proposed model was simulated for approximately 100 days for different fractal fractional-order values $\tau_1$ and $\tau_2$. According to these parameters, $E_0 = (0.1192, 0.2961, 0, 0, 0)$ is asymptomatically stable. As predicted, the solutions of (2.1) converge to the unique disease-free equilibrium $E_0$. The biological implication is that we need to bring $\operatorname{R}_0$ to below 1 to ensure a reduction of the disease in the country.

Figures 2–6 show the results of the fractal-fractional Atangana-Baleanu and the Adams-Bashforth-Moulton methods for pneumonia transmission, with (a) $\tau_1 = 0.7$, $\tau_2 = 0.95$, (b) $\tau_1 = 0.75$, $\tau_2 = 0.95$, (c) $\tau_1 = 0.85$, $\tau_2 = 0.95$ (d) $\tau_1 = 0.95$, $\tau_2 = 0.95$ and (e) $\tau_1 = 1$, $\tau_2 = 0.95$. Comparisons between the ordinary differential system, ABC fractal, and fractional derivative, can also be seen in Figures 2–6.

Figures 2–6 show the influence of varying $\tau_{1}$ between 0.7 and 1 with a fixed $\tau_{2} = 0.95$ on model dynamics. The black curve in each of these figures represents the numerical results of model (2.1) when the fractional order is equal to 1. From the results of Figures 2–6, it follows that the variation of the fractional parameter has a great impact on the quantitative dynamics of the model. Indeed, in Figure 5, the classes of infected humans peak after 10 years and decrease according to the decrease of the fractional parameter $\tau_{1}$.

Lemmas 2 and 3 are validated numerically in Figures 2–6. It is clear that varying the fractional order parameter $\tau_{1}$ does not influence the model dynamics whenever $\mathscr{R}_{0} = 0.0083 < 1$. Indeed, whatever the value of $\tau_{1}$, the infected compartments tend to zero asymptotically whenever $\mathscr{R}_{0} = 0.0083 < 1$. This validates the fact that the pneumonia-free equilibrium of the fractional model is globally asymptotically stable whenever $\mathscr{R}_{0} = 0.0083 < 1$.

Figures 7–11 show the phase plots ($\operatorname{x}-\operatorname{y}-\operatorname{z}-\operatorname{u}-\operatorname{v}$) for different values, i.e., (a) $\tau_1 = 0.7$, $\tau_2 = 0.95$, (b) $\tau_1 = 0.75$, $\tau_2 = 0.95$, (c) $\tau_1 = 0.85$, $\tau_2 = 0.95$ (d) $\tau_1 = 0.95$, $\tau_2 = 0.95$ and (e) $\tau_1 = 1$, $\tau_2 = 0.95$.

Comparison between ordinary differential system, ABC, ABC fractal fractional derivative schemes can be seen in Figures 12–16.

5.   Conclusions

A fractal fractional-order mathematical model based on the Atangana-Baleanu operator was constructed to describe pneumonia transmission in a population. The Caputo operator was used to analyze the dynamics of the virus, and a fractal fractional derivative was used to maximize the number of recovered populations. For models of infectious diseases, the vaccination rate coefficient is considered as a control to reduce the disease burden. It is important to prove the existence of optimal control, characterize the optimal control, prove the uniqueness of optimal control and compute the optimal control numerically. The model was subjected to dynamic analysis, and the results show that the rate at which susceptible individuals contract an infectious disease is the most significant parameter. In this simulation, the value of $\mathscr{R}_{0} = 0.0083 < 1$, which is smaller than 1. As you can see, disease spread is controlled, and the number of infected people is reduced to zero. We also see that each function tends to its equilibrium point, and that the equilibrium point becomes stable as the system approaches its equilibrium point. We also note that the total number of susceptible humans decreases rapidly according to the increase of the fractional parameter (Figure 2).

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research, for the financial support for this research under grant number (37622-BSRC-FFT-2023) during the academic year 1445 AH/2023 AD.

Conflict of interest

The authors declare that they have no conflict of interest.