Fractal–fractional model for analysis of tuberculosis infection using Ethiopia incidence data

Abstract

In this study, a fractional‐order SE{I}_{u}ITR model was proposed to examine the transmission of tuberculosis (TB) in the community. The proposed model was rigorously examined for well‐posedness. The basic reproduction number was also derived, and the disease‐free equilibrium was obtained. The local and global stabilities of the proposed model were also examined. The proposed model was fitted with the actual TB incidence in Ethiopia, while numerical simulations were carried out in MATLAB to validate the results obtained in this study. Local and global sensitivity analyses were carried out to determine the key parameters that influence the transmission of TB in the community. Phase‐plane portraits and phase‐plane diagrams demonstrate the dynamic process of the state variables. The results indicate that as the fractional‐order increases, the number of infections decreases due to the memory effect of early detection and treatment, while as the order decreases, the number of infections increases. Early detection and treatment have a profound effect on recovery and prevalence.

Graphical abstract

Keywords

convolution theorem / Euler fractional method / fixed point theory / fractional‐order / Laplace transform / numerical simulation / tuberculosis

1 INTRODUCTION

Tuberculosis (TB) is one of the most dangerous infectious diseases caused by the bacterium Mycobacterium tuberculosis [1–5]. One of the most important features of TB is its latent form, where people who have been in contact with the bacteria or are carriers may not even know their status.

The primary infection stage of TB is when the immune cells in the body recognize and fight off the invading bacteria. In some instances, the immune system is able to destroy the bacteria entirely [6,7]. However, some of the bacteria may survive and live inside the body. One of the biggest problems in the fight against TB is that many people infected with the disease do not show symptoms during the early stages of infection. When symptoms do appear, they are usually mild and may include influenza‐like symptoms such as low‐grade fever, tiredness, and coughing. After the primary infection, a person may develop latent TB infection, where the immune system is able to hold off the bacteria by creating a barrier in the lungs. During this period, the bacteria remain alive but inactive, and no clinical symptoms are observed. Active TB develops when the immune system fails to suppress the infection, allowing the bacteria to proliferate and cause disease in the lungs or other parts of the body.

Active TB may occur shortly after initial infection or emerge months or even years later due to reactivation of latent infection. The symptoms of active TB typically worsen gradually and may include persistent cough, hemoptysis, chest pain, fever, chills, night sweats, loss of weight, loss of appetite, fatigue, and general malaise [8,9].

The symptoms of this type of TB depend on the organs that are affected and may include the kidneys, liver, central nervous system, heart, lymph nodes, bones and joints, skin, blood vessels, and larynx. The symptoms include fever, night sweats, loss of weight, fatigue, pain, and loss of organ function. In children, especially infants, the symptoms of TB can be very different from those of adults. In adolescents, the symptoms may resemble those of adults, while in younger children, the symptoms may include fever, failure to gain weight, developmental delays, and neurological problems due to TB meningitis [7,8].

Mathematical modeling is a key component in the study of complex phenomena in the real world by formulating them in a mathematical framework [10]. In the field of epidemiology, mathematical models have been extensively employed to describe the dynamics of disease transmission and assess the efficacy of control measures [11]. Traditional integer‐order dynamical systems have been employed to model TB transmission, including the role of exogenous reinfection and endogenous reactivation of recovered or pseudo‐recovered individuals [3,12]. Nevertheless, these models have limitations in describing the memory and hereditary properties of TB, which is associated with long latent periods and delayed treatment responses.

Fractional‐order models have been recognized as useful tools for implementing memory effects and have been extensively employed to improve the accuracy and realism of infectious disease models [13–15]. Specifically, the non‐singular memory kernels of the Caputo–Fabrizio (CF) and Atangana–Baleanu–Caputo (ABC) operators make it possible to implement memory effects without encountering singularities, which is advantageous for modeling long‐term dynamics of infectious diseases.

The fractional‐order \alpha measures the intensity of memory activation by incorporating the entire history of disease evolution into the current state of the system. Lower values of α indicate higher memory effects, which represent important epidemiological processes such as disease latency, waning immunity, and pathogen reactivation [16,17]. These properties are particularly important in the case of TB, in which treatment compliance and late diagnosis are key factors in the persistence of the disease.

Recent research has emphasized the importance of awareness, reinfection, and effective intervention strategies in controlling TB transmission [9,12,13,16–22]. However, TB still constitutes a major health problem worldwide, especially in developing nations. It is in light of these challenges that this paper proposes a CF fractional TB model that includes both informed and uninformed compartments to explore the effect of memory activation on TB transmission based on incidence data from Ethiopia.

2 RESULTS

2.1 Non‐negativity of solutions

In this section, we prove that all state variables of the system are non‐negative under the fractional operator.

Theorem 1.(Non‐negativity of solutions). The solutions of the TB system of Equation (118) with anticipated initial conditions behave non‐negative for all t\ge 0.

Proof The proof is based on the generalized mean value theorem for the CF fractional derivative, which states that if a function reaches a minimum at a point, then the fractional derivative of the function at that point is non‐positive.

First to show the non‐negativity of S\left(t\right), suppose that there is a first time t_1>0 such that

For system of Equation (118),

\begin{array}{}\text{(2)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }S\left(t_1\right)=\mathrm{\Lambda }>0.\end{array}

However, since S\left(t\right) attains its minimum at t_1, the generalized fractional mean value theorem implies

which is a contradiction. Thus,

\begin{array}{}\text{(4)}& S\left(t\right)\ge 0,\mathit{\forall }t\ge 0.\end{array}

Again to show the non‐negativity of E\left(t\right) assume there exists t_2 such that

\begin{array}{}\text{(5)}& E\left(t_2\right)=0.\end{array}

Then,

\begin{array}{}\text{(6)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }E\left(t_2\right)=\beta S\left(t_2\right)\left(I_{\mathrm{u}}\left(t_2\right)+\theta I\left(t_2\right)\right)\ge 0.\end{array}

Thus, E\left(t\right) cannot become negative.

To show the non‐negativity of I_{\mathrm{u}}\left(t_3\right), if I_{\mathrm{u}}\left(t_3\right)=0, then

\begin{array}{}\text{(7)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}\left(t_3\right)=\kappa E\left(t_3\right)\ge 0.\end{array}

Thus,

\begin{array}{}\text{(8)}& I_{\mathrm{u}}\left(t\right)\ge 0.\end{array}

Again to show the non‐negativity of I\left(t\right), if I\left(t_4\right)=0, then

\begin{array}{}\text{(9)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }I\left(t_4\right)=\varphi I_{\mathrm{u}}\left(t_4\right)+\omega E\left(t_4\right)\ge 0.\end{array}

Thus,

\begin{array}{}\text{(10)}& I\left(t\right)\ge 0.\end{array}

Also, to show the non‐negativity of T\left(t\right), if T\left(t_5\right)=0, then

\begin{array}{}\text{(11)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }T\left(t_5\right)=\rho I\left(t_5\right)\ge 0.\end{array}

Thus,

\begin{array}{}\text{(12)}& T\left(t\right)\ge 0.\end{array}

To show the non‐negativity of R\left(t\right), if R\left(t_6\right)=0, then

\begin{array}{}\text{(13)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }R\left(t_6\right)=\eta T\left(t_6\right)\ge 0.\end{array}

Gives R(t) ≥ 0 and satisfies the non‐negativity condition.

Therefore, the solutions remain biologically feasible and this completes the proof.

2.2 Invariant region using Mittag‐Leffler function

Theorem 2.(Invariant region). Consider system of Equation (118) with non‐negative initial conditions. Then the region

\begin{array}{}\text{(14)}& \mathrm{\Omega }=\left\{\left(S,E,I_{\mathrm{u}},I,T,R\right)\mathit{\in }{\mathbb{R}}_{+}^6:0\mathit{\le }N\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu }\right\},\end{array}

is positively invariant.

Proof Summing left and right of system of Equation (118), we obtain

\begin{array}{}\text{(15)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }N\left(t\right)=\mathrm{\Lambda }-\mu N\left(t\right).\end{array}

Since I_{\mathrm{u}},I,T\ge 0 by non‐negativity. Consider the auxiliary scalar equation

\begin{array}{}\text{(16)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }X\left(t\right)=\mathrm{\Lambda }-\mu X\left(t\right),X\left(t\right),X\left(0\right)=N\left(0\right)\ge 0.\end{array}

Applying the Laplace transform to Equation (16):

\begin{array}{}\text{(17)}& \mathcal{L}\left\{{}_0{}^{\text{CF}}{D}_t^{\alpha }X\left(t\right)\right\}=\frac{s\bar{X}\left(s\right)-N\left(0\right)}{\left(1-\alpha \right)s+\alpha }=\frac{\mathrm{\Lambda }}{s}-\mu \bar{X}\left(s\right).\end{array}

Solving for \bar{X}\left(s\right):

\begin{array}{}\text{(18)}& \bar{X}\left(s\right)=\frac{N\left(0\right)+\mathrm{\Lambda }\left(1-\alpha \right)/\mu }{s\left(1+\mu \left(1-\alpha \right)\right)+\mu \alpha }+\frac{\alpha \mathrm{\Lambda }/s}{s\left(1+\mu \left(1-\alpha \right)\right)+\mu \alpha }.\end{array}

Using partial fraction decomposition and the Laplace transform of the Mittag‐Leffler function, the inverse transform gives

\begin{array}{}\text{(19)}& X\left(t\right)=\frac{\mathrm{\Lambda }}{\mu }+\left(N\left(0\right)-\frac{\mathrm{\Lambda }}{\mu }\right)E_1\left(-\frac{\mu \alpha }{1+\mu \left(1-\alpha \right)}t\right),\end{array}

where E_1\left(-x\right) is the Mittag‐Leffler function of order one. Since for x\ge 0, we have

\begin{array}{}\text{(20)}& 0\le X\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu },\mathit{\forall }t\ge 0.\end{array}

By the fractional comparison principle, since

\begin{array}{}\text{(21)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }N\left(t\right)\le {}_0{}^{\text{CF}}{D}_t^{\alpha }X\left(t\right),N\left(0\right)=\mathrm{X}\left(0\right).\end{array}

It follows that

\begin{array}{}\text{(22)}& N\left(t\right)\le X\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu },\mathit{\forall }t\ge 0.\end{array}

Hence, the region

\begin{array}{}\text{(23)}& \mathrm{\Omega }=\left\{\left(S,E,I_{\mathrm{u}},I,T,R\right)\mathit{\in }{R}_{+}^6:0\le N\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu }\right\},\end{array}

is positively invariant.

2.3 Existence and uniqueness of solutions

One of the known conditions for the well‐posedness of the model is the existence and uniqueness of the solution. Consider the CF fractional epidemiological system of Equation (118) with non‐negative initial conditions. The existence and uniqueness of the solutions can be guaranteed by the Banach fixed‐point theorem.

Theorem 3.(Existence and uniqueness). Let

\begin{array}{}\text{(24)}& S\left(0\right),E\left(0\right),I_{\mathrm{u}}\left(0\right),I\left(0\right),T\left(0\right),R\left(0\right)\ge 0,\end{array}

and assume all model parameters are non‐negative constants. Then there exists a unique solution

\begin{array}{}\text{(25)}& \left(S\left(t\right),E\left(t\right),I_{\mathrm{u}}\left(t\right),T\left(t\right),R\left(t\right)\right)\mathit{\in }C\left[0,\tau \right],{\mathbb{R}}_{+}^6,\end{array}

for some \tau >0.

Proof We write system of Equation (118) in the abstract form

\begin{array}{}\text{(26)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }\mathit{X}\left(t\right)=\mathit{F}\left(t,\mathit{X}\left(t\right)\right),\mathit{X}\left(0\right)={\mathit{X}}_0,\end{array}

where

\begin{array}{}\text{(27)}& \mathit{X}\left(t\right)=\left(\begin{array}{c}S\left(t\right)\\ E\left(t\right)\\ I_{\mathrm{u}}\left(t\right)\\ I\left(t\right)\\ T\left(t\right)\\ R\left(t\right)\end{array}\right),{\mathit{X}}_0=\left(\begin{array}{c}S\left(0\right)\\ E\left(0\right)\\ I_{\mathrm{u}}\left(0\right)\\ I\left(0\right)\\ T\left(0\right)\\ R\left(0\right)\end{array}\right),\end{array}

and F is the vector of right‐hand side functions of the system of Equation (118). By definition of the CF integral, this system of equations can be re‐written in the equivalent integral equation

\begin{array}{}\text{(28)}& \mathit{z}\left(t\right)={\mathit{z}}_0+\frac{1-\alpha }{\alpha }{\int }_0^t\mathit{F}\left(\tau ,\mathit{z}\left(\tau \right)\right)d\tau .\end{array}

Define an operator or function \varkappa on the Banach space

by

\begin{array}{}\text{(30)}& \left(\varkappa \mathit{Y}\right)\left(t\right)={\mathit{Y}}_0+\frac{1-\alpha }{\alpha }{\int }_0^t\mathit{F}\left(\tau ,\mathit{Y}\left(\tau \right)\right)d\tau .\end{array}

Since \mathit{F} is Lipschitz continuous with respect to \mathit{z} in any bounded subset of {\mathbb{R}}_{+}^6, there exists a constant L>0 such that

\begin{array}{}\text{(31)}& \Vert \mathit{F}\left(t,{\mathit{Y}}_1\left(t\right)\right)-\left(\mathit{F}\left(t,{\mathit{Y}}_2\left(t\right)\right)\right)\Vert \le L\Vert {\mathit{Y}}_1\left(t\right)-{\mathit{Y}}_2\left(t\right)\Vert ,\forall t\in \left[0,\tau \right].\end{array}

Consequently, for all {\mathit{Y}}_{\mathbf{1}},{\mathit{Y}}_{\mathbf{2}}\mathit{\in }\mathcal{B}, we have

\begin{array}{}\text{(31)}& \Vert \left(\varkappa {\mathit{Y}}_1\right)\left(t\right)-\left(\varkappa {\mathit{Y}}_2\right)\left(t\right)\Vert \le \frac{1-\alpha }{\alpha }{\int }_0^t\Vert \mathit{F}\left(\tau ,{\mathit{Y}}_1\left(\tau \right)\right)-\left(\mathit{F}\left(\tau ,{\mathit{Y}}_2\left(\tau \right)\right)\right)\Vert d\tau \end{array}

\begin{array}{}\text{(32)}& \le \frac{1-\alpha }{\alpha }L{\int }_0^t\Vert {\mathit{Y}}_1\left(\tau \right)-{\mathit{Y}}_2\left(\tau \right)\Vert d\tau ,\end{array}

\begin{array}{}\text{(33)}& \le \frac{1-\alpha }{\alpha }L\tau {\Vert {\mathit{Y}}_1-{\mathit{Y}}_2\Vert }_{\infty },\end{array}

by choosing appropriately, that \tau >0

the operator \varkappa is a contraction operator on \mathcal{B}. According to the Banach fixed point theorem, a contraction operator on a complete metric space has a unique fixed point. Therefore, the contraction operator \varkappa on \mathcal{B} possesses a distinct fixed point {\mathit{z}}^{\ast }\mathit{\in }\mathcal{B} satisfying

\begin{array}{}\text{(35)}& {\mathit{z}}^{\ast }=\mathit{\varkappa }{\mathit{z}}^{\ast }.\end{array}

This fixed point corresponds to the unique solution to the system of Equation (118) on the interval \left[0,t\right].

2.4 Disease‐free equilibrium

The disease‐free equilibrium (DFE) (X_0) represents a state where there is no disease in the population [14,15,18,23–27]. This means that the populations in the exposed and infectious classes are zero, that is, E=I_{\mathrm{u}}=I=T=0. At equilibrium, the CF fractional derivatives vanish. Setting

\begin{array}{}\text{(36)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }S={}_0{}^{\text{CF}}{D}_t^{\alpha }E={}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}={}_0{}^{\text{CF}}{D}_t^{\alpha }I={}_0{}^{\text{CF}}{D}_t^{\alpha }T={}_0{}^{\text{CF}}{D}_t^{\alpha }R=0.\end{array}

We obtain the algebraic system of equation

\begin{array}{}\text{(37)}& 0=\mathrm{\Lambda }-\mu S,0=-\mu R.\end{array}

Therefore, the solution of the preceding system of equation yields

\begin{array}{}\text{(38)}& X_0=\left(S,E,I_{\mathrm{u}},I,T,R\right)=\left(\frac{\mathrm{\Lambda }}{\mu },0,0,0,0,0\right).\end{array}

2.5 Basic reproduction number

In a completely susceptible population, the expected number of secondary infections caused by one infected person is defined as the basic reproduction number R_0, which is calculated as the spectral radius of the next generation matrix [28–34]. For the model, let’s consider \mathit{F} the matrix of the new infection terms and \mathit{V} the transition matrix:

\begin{array}{}\text{(39)}& \mathit{F}=\left(\begin{array}{c}\beta S\left(I_{\mathrm{u}}+\theta I\right)\\ 0\\ 0\\ 0\end{array}\right),\mathit{V}=\left(\begin{array}{c}\left(\omega +\kappa +\mu \right)E\\ -\kappa E+\left(\varphi +\delta +\mu \right)I_{\mathrm{u}}\\ -\varphi I_{\mathrm{u}}-\omega E+\left(\rho +\delta +\mu \right)I\\ -\rho I+\left(\eta +\delta +\mu \right)T\end{array}\right).\end{array}

The Jacobian matrices of vectors f and v are the matrices \mathit{F} and \mathit{V} respectively and given by

\begin{array}{}\text{(40)}& \mathit{F}=\left(\begin{array}{cccc}0& \beta S& \beta \theta S& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),\mathit{V}=\left(\begin{array}{cccc}\omega +\kappa +\mu & 0& 0& 0\\ -\kappa & \varphi +\delta +\mu & 0& 0\\ -\omega & -\varphi & \rho +\delta +\mu & 0\\ 0& 0& -\rho & \eta +\delta +\mu \end{array}\right).\end{array}

Computation of R_0=\mathit{F}{\mathit{V}}^{-1} gives

\begin{array}{}\text{(41)}& R_0=\frac{\beta S_0}{\omega +\kappa +\mu }\left(\frac{\kappa }{\varphi +\delta +\mu }+\frac{\varphi }{\rho +\delta +\mu }+\frac{\omega }{\rho +\delta +\mu }\right).\end{array}

In the R_0, substituting S_0 with \frac{\mathrm{\Lambda }}{\mu } gives,

\begin{array}{}\text{(42)}& R_0=\frac{\beta \mathrm{\Lambda }}{\mu \left(\omega +\kappa +\mu \right)}\left(\frac{\kappa }{\varphi +\delta +\mu }+\frac{\varphi }{\rho +\delta +\mu }+\frac{\omega }{\rho +\delta +\mu }\right).\end{array}

This formula expresses the expected number of secondary infections in a completely susceptible population. If , the infection is not sustained and will die out. If R_0>1, then the disease can invade the population.

2.6 Endemic equilibrium in terms of R_0

In the state of equilibrium, the CF fractional derivatives are zero:

\begin{array}{}\text{(43)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }S={}_0{}^{\text{CF}}{D}_t^{\alpha }E={}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}={}_0{}^{\text{CF}}{D}_t^{\alpha }I={}_0{}^{\text{CF}}{D}_t^{\alpha }T={}_0{}^{\text{CF}}{D}_t^{\alpha }R=0.\end{array}

Setting the right‐hand sides of the system of Equation (118) to zero, we have the following algebraic system of equations:

\begin{array}{}\text{(44)}& \begin{array}{ll}0=& \mathrm{\Lambda }-\beta S\left(I_{\mathrm{u}}+\theta I\right)-\mu S,\\ 0=& \beta S\left(I_{\mathrm{u}}+\theta I\right)-\left(\omega +\kappa +\mu \right)E,\\ 0=& \kappa E-\left(\varphi +\delta +\mu \right)I_{\mathrm{u}},\\ 0=& \varphi I_{\mathrm{u}}+\omega E-\left(\rho +\delta +\mu \right)I,\\ 0=& \rho I-\left(\eta +\delta +\mu \right)T,\\ 0=& \eta T-\mu R.\end{array}\end{array}

The last three equations of Equation (118) can be expressed in terms of the variable E as

\begin{array}{}\text{(45)}& I_{\mathrm{u}}=\frac{\kappa }{\varphi +\delta +\mu }E,I=\frac{\varphi I_{\mathrm{u}}+\omega E}{\rho +\delta +\mu }=\frac{\varphi \kappa +\omega \left(\varphi +\delta +\mu \right)}{\left(\varphi +\delta +\mu \right)\left(\rho +\delta +\mu \right)}E,T=\frac{\rho }{\eta +\delta +\mu }I,R=\frac{\eta }{\mu }T.\end{array}

Substituting these expressions for variables I_{\mathrm{u}} and I into the first two equations of Equation (118), we get

\begin{array}{}\text{(46)}& \beta S\left(I_{\mathrm{u}}+\theta I\right)=\left(\omega +\kappa +\mu \right)E.\end{array}

Using the expression for R_0:

\begin{array}{}\text{(47)}& R_0=\frac{\beta S_0}{\omega +\kappa +\mu }\left(\frac{\kappa }{\varphi +\delta +\mu }+\frac{\varphi }{\rho +\delta +\mu }+\frac{\omega }{\rho +\delta +\mu }\right),S_0=\frac{\mathrm{\Lambda }}{\mu },\end{array}

the susceptible population at endemic equilibrium, it can be written as

\begin{array}{}\text{(48)}& S=\frac{S_0}{R_0}.\end{array}

Thus, the endemic equilibrium \left({\mathit{X}}_{\mathit{e}}\right) is given by

\begin{array}{}\text{(49)}& {\mathit{X}}_{\mathit{e}}=\left(S,E,I_{\mathrm{u}},I,T,R\right)=\left(\frac{S_0}{R_0},E,\frac{\kappa }{\varphi +\delta +\mu }E,\frac{\varphi \kappa +\omega \left(\varphi +\delta +\mu \right)}{\left(\varphi +\delta +\mu \right)\left(\rho +\delta +\mu \right)}E,\frac{\rho }{\eta +\delta +\mu }I,\frac{\eta }{\mu }T\right),\end{array}

where E can be computed from

\begin{array}{}\text{(50)}& E=\frac{\mu \left(R_0-1\right)S_0}{\omega +\kappa +\mu +\beta \left(\frac{\kappa }{\varphi +\delta +\mu }+\theta \frac{\varphi \kappa +\omega \left(\varphi +\delta +\mu \right)}{\left(\varphi +\delta +\mu \right)\left(\rho +\delta +\mu \right)}\right)S_0}.\end{array}

2.7 Local stability of equilibria

Theorem 4.(Local stability of TB‐free steady state). Let {\mathit{X}}_{\mathbf{0}}=\left(S,E,I_{\mathrm{u}},I,T,R\right)=\left(S_0,0,0,0,0,0\right) the DFE with S_0=\frac{\mathrm{\Lambda }}{\mu }. Then,

•If , then the DFE is locally asymptotically stable.

•If R_0>1, then the DFE is unstable.

Proof Let {\mathit{X}}_{\mathit{I}}=\left(\begin{array}{c}E\\ I_{\mathrm{u}}\\ I\\ T\end{array}\right) the TB‐infected compartments and we construct the Lyapunov function,

\begin{array}{}\text{(51)}& V=E+I_{\mathrm{u}}+I+T,\end{array}

which is non‐negative and has an equilibrium value of 0 at the DFE. We calculate the CF fractional derivative of V with respect to time t

\begin{array}{}\text{(52)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }V={}_0{}^{\text{CF}}{D}_t^{\alpha }E+{}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}+{}_0{}^{\text{CF}}{D}_t^{\alpha }I+{}_0{}^{\text{CF}}{D}_t^{\alpha }T.\end{array}

By substituting the system equations and evaluating at S\mathit{\le }{S}_0, we can write that

By applying the standard Lyapunov theory for fractional‐order systems, we can conclude that, since {}_0{}^{\text{CF}}{D}_t^{\alpha }V\mathit{\le }0 and V is positive definite, the DFE is locally asymptotically stable for , and unstable if R_0>1.

Theorem 5.(Local stability of the endemic equilibrium). Let {\mathit{X}}_{\mathit{e}} be the endemic equilibrium expressed in terms of R_0>1. Then the endemic equilibrium is locally asymptotically stable.

Proof Let the deviation from the endemic equilibrium be given by

\begin{array}{}\text{(54)}& \check{\mathit{X}}=\mathit{X}-{\mathit{X}}_{\mathit{e}}.\end{array}

Define Lyapunov function to show the local stability state

\begin{array}{}\text{(55)}& V=\sum _i{\check{X}}_i-{\check{X}}_i\ln \frac{{\check{X}}_i+X_{e,i}}{X_{e,i}},\end{array}

where the sum is over all compartments i=S,E,I_{\mathrm{u}},I,T,R. This function is positive definite and vanishes at \mathit{X}={\mathit{X}}_{\mathit{e}}. Differentiating along the trajectories of system Equation (118) and using the relation between compartments at the endemic equilibrium expressed in terms of R_0, we obtain

\begin{array}{}\text{(56)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }V\mathit{\le }0,\end{array}

equality if and only if \check{\mathit{X}}=\mathbf{0}. For fractional‐order systems, the endemic equilibrium {\mathit{X}}_{\mathit{e}} is thus locally asymptotically stable whenever R_0>1, according to Lyapunov’s direct technique.

2.8 Global stability of the DFE

Theorem 6.(Global stability of DFE). If {\mathit{X}}_{\mathit{e}}=\left(S,E,I_{\mathrm{u}},I,T,R\right)=\left(S_0,0,0,0,0,0\right) is the DFE with S_0=\frac{\mathrm{\Lambda }}{\mu }, then the TB‐free steady state is globally asymptotically stable in the feasible region of

\begin{array}{}\text{(57)}& \mathrm{\Omega }=\left\{\left(S,E,I_{\mathrm{u}},I,T,R\right)\mathit{\in }{\mathbb{R}}_{+}^6:S+E+I_{\mathrm{u}}+I+T+R\mathit{\le }{S}_0\right\},\end{array}

the basic reproduction number R_0\mathit{\le }1.

Proof To study the global behavior of the infection, we focus on the infected population given by the vector

\begin{array}{}\text{(58)}& {\mathit{X}}_{\mathit{I}}=\left(\begin{array}{c}E\\ I_{\mathrm{u}}\\ I\\ T\end{array}\right).\end{array}

To prove the stability, we define the Lyapunov function

\begin{array}{}\text{(59)}& V=E+I_{\mathrm{u}}+I+T.\end{array}

This function is positive definite in the feasible region \Omega and only attains zero at the DFE when E=I_{\mathrm{u}}=I=T=0.

Taking the CF fractional derivative with respect to time t along the solutions of the system of equations, we obtain

\begin{array}{}\text{(60)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }V={}_0{}^{\text{CF}}{D}_t^{\alpha }E+{}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}+{}_0{}^{\text{CF}}{D}_t^{\alpha }I+{}_0{}^{\text{CF}}{D}_t^{\alpha }T.\end{array}

Replacing the system equations for the infected population, we have

\begin{array}{}\text{(61)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }V=\beta S\left(I_{\mathrm{u}}+\theta I\right)-\left(\omega +\kappa +\mu \right)E-\left(\phi +\delta +\mu \right)I_{\mathrm{u}}-\left(\rho +\delta +\mu \right)I-\left(\eta +\delta +\mu \right)T\le \beta S_0\left(I_{\mathrm{u}}+\theta I\right)-\left(\omega +\kappa +\mu \right)E-\left(\phi +\delta +\mu \right)I_{\mathrm{u}}-\left(\rho +\delta +\mu \right)I-\left(\eta +\delta +\mu \right)T,\end{array}

where we have used the fact that S\mathit{\le }{S}_0 in the feasible region. Using the definition of the basic reproduction number,

\begin{array}{}\text{(62)}& R_0=\frac{\beta S_0}{\omega +\kappa +\mu }\left(\frac{\kappa }{\varphi +\delta +\mu }+\frac{\varphi }{\rho +\delta +\mu }+\frac{\omega }{\rho +\delta +\mu }\right),\end{array}

we obtain

\begin{array}{}\text{(63)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }V\mathit{\le }\left(\omega +\kappa +\mu \right)\left(R_0-1\right)E.\end{array}

Thus, when R_0\mathit{\le }1, we have {}_0{}^{\text{CF}}{D}_t^{\alpha }V\mathit{\le }0 and equality occurs only when E=I_{\mathrm{u}}=I=T=0, that is, at the DFE.

Using LaSalle’s invariance principle for fractional systems, all trajectories in \Omega will converge to the largest invariant set where {}_0{}^{\text{CF}}{D}_t^{\alpha }V=0, which is the DFE.

The global asymptotic stability of the DFE {\mathit{X}}_{\mathbf{0}} is guaranteed if R_0\mathit{\le }1.

2.9 Global stability of TB model endemic equilibrium

Theorem 7.If the assumption R_0\mathit{\ge }1 holds, then the TB model endemic equilibrium

\begin{array}{}\text{(64)}& {\mathit{X}}_{\mathit{e}}=\left(S_e,E_e,I_{\mathrm{u},e},I_e,T_e,R_e\right),\end{array}

of the TB dynamic transmission fractional‐order system with CF fractional derivative guarantees a globally asymptotically stable behavior in the feasible region

\begin{array}{}\text{(65)}& \Omega =\left\{\left(S,E,I_{\mathrm{u}},I,T,R\right)\mathit{\in }{\mathbb{R}}_{+}^6\right\}.\end{array}

Proof In order to show the guarantee of the global stability of the system, we consider the following Lyapunov function

\begin{array}{}\text{(66)}& V=\left(S-S_e-S_e\ln \frac{S}{S_e}\right)+\left(E-E_e-E_e\ln \frac{E}{E_e}\right)+\left(I_{\mathrm{u}}-I_{\mathrm{u},e}-I_{\mathrm{u},e}\ln \frac{I_{\mathrm{u}}}{I_{\mathrm{u},e}}\right)+\left(I-I_e-I_e\ln \frac{I}{I_e}\right)+\left(T-T_e-T_e\ln \frac{T}{T_e}\right)+\left(R-R_e-R_e\ln \frac{R}{R_e}\right).\end{array}

Accordingly, the function satisfies,

\begin{array}{}\text{(67)}& V\mathit{\ge }0,\end{array}

and

\begin{array}{}\text{(68)}& V=0,\text{if}\text{and}\text{only}\text{if}\mathit{X}={\mathit{X}}_{\mathit{e}}.\end{array}

Since the function

\begin{array}{}\text{(69)}& f\left(x\right)=x-x_e-x_e\ln \left(\frac{x}{x_e}\right),\end{array}

is non‐negative for all x > 0.

Applying the CF fractional derivative along system trajectories,

\begin{array}{}\text{(70)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }V=\left(1-\frac{S_e}{S}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }S+\left(1-\frac{E_e}{E}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }E+\left(1-\frac{I_{\mathrm{u},e}}{I_{\mathrm{u}}}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}+\left(1-\frac{I_e}{I}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }I+\left(1-\frac{T_e}{T}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }T+\left(1-\frac{R_e}{R}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }R.\end{array}

On the other hand, from the equation,

\begin{array}{}\text{(71)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }E=\beta S\left(I_{\mathrm{u}}+\theta I\right)-\left(\omega +\kappa +\mu \right)E.\end{array}

At equilibrium,

\begin{array}{}\text{(72)}& \beta S_e\left(I_{\mathrm{u},e}+\theta I_e\right)=\left(\omega +\kappa +\mu \right)E_e.\end{array}

Thus,

\begin{array}{}\text{(73)}& \left(1-\frac{E_e}{E}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }E\mathit{\le }0,\end{array}

with equality holds if and only if

\begin{array}{}\text{(74)}& E=E_e.\end{array}

Again, for I_{\mathrm{u}},

\begin{array}{}\text{(75)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}=\kappa E-\left(\varphi +\delta +\mu \right)I_{\mathrm{u}},\end{array}

and using

\begin{array}{}\text{(76)}& \kappa E_e=\left(\varphi +\delta +\mu \right)I_{\mathrm{u},e}.\end{array}

Hence,

\begin{array}{}\text{(77)}& \left(1-\frac{I_{\mathrm{u},e}}{I_{\mathrm{u}}}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}\mathit{\le }0.\end{array}

Also for I,

\begin{array}{}\text{(78)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }I=\varphi I_{\mathrm{u}}+\omega E-\left(\rho +\delta +\mu \right)I,\end{array}

and

\begin{array}{}\text{(79)}& \varphi I_{\mathrm{u},e}+\omega E_e=\left(\rho +\delta +\mu \right)I_e.\end{array}

Thus,

\begin{array}{}\text{(80)}& \left(1-\frac{I_e}{I}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }I\mathit{\le }0.\end{array}

Also for T,

\begin{array}{}\text{(81)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }T=\rho I-\left(\eta +\delta +\mu \right)T,\end{array}

and

\begin{array}{}\text{(82)}& \rho I_e=\left(\eta +\delta +\mu \right)T_e,\end{array}

thus,

\begin{array}{}\text{(83)}& \left(1-\frac{T_e}{T}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }T\mathit{\le }0.\end{array}

But,

\begin{array}{}\text{(84)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }S=\mathrm{\Lambda }-\beta S\left(I_{\mathrm{u}}+\theta I\right)-\mu S.\end{array}

Considering the equilibrium condition,

\begin{array}{}\text{(85)}& \mathrm{\Lambda }=\beta S_e\left(I_{\mathrm{u},e}+\theta I_e\right)+\mu S_e,\end{array}

yields,

\begin{array}{}\text{(86)}& \left(1-\frac{S_e}{S}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }S\mathit{\le }0.\end{array}

Again,

\begin{array}{}\text{(87)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }R=\eta T-\mu R,\end{array}

and

\begin{array}{}\text{(88)}& \eta T_e=\mu R_e.\end{array}

Thus,

\begin{array}{}\text{(89)}& \left(1-\frac{R_e}{R}\right){}_0{}^{\text{CF}}{D}_t^{\alpha }R\mathit{\le }0.\end{array}

Summing all components yields,

\begin{array}{}\text{(90)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }V\mathit{\le }0,\end{array}

with equality holds if and only if

\begin{array}{}\text{(91)}& \left(S,E,I_{\mathrm{u}},I,T,R\right)=\left(S_e,E_e,I_{\mathrm{u},e},I_e,T_e,R_e\right).\end{array}

Therefore, by LaSalle’s invariance principle for fractional‐order systems, the endemic equilibrium is globally asymptotically stable whenever

\begin{array}{}\text{(92)}& R_0>1.\end{array}

2.10 Numerical scheme of model

The fractional TB model was solved numerically using the fractional Euler method adapted for the CF fractional derivative. This approach is well‐suited for epidemiological models incorporating memory effects, as it accounts for past states of the system, providing a more realistic representation of disease progression [35–37].

Let the fractional‐order be \alpha \mathit{\in }(0,1]. Consider the following general fractional differential equation:

\begin{array}{}\text{(93)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }X\left(t\right)=F\left(t,X\left(t\right)\right),X\left(0\right)=X_0.\end{array}

Using a fixed step size h, the CF fractional derivative can be approximated by the numerical fractional Euler scheme that is simple and suitable for the analysis of various epidemiological models [35,36]:

\begin{array}{}\text{(94)}& X_{n+1}=X_n+\frac{\left(1-\alpha \right)}{M\left(\alpha \right)}F\left(t_n,X_n\right)+\frac{\alpha h}{M\left(\alpha \right)}F\left(t_n,X_n\right),\end{array}

where M\left(\alpha \right) is the normalization function of the derivative. This method was iteratively applied to the entire compartments of the TB model, including S,E,I_{\mathrm{u}},I,T, and R. This method also guarantees the positivity and feasibility of the solution variables. The method was implemented in MATLAB to efficiently simulate the fractional system. The fractional weight \alpha offers flexibility to interpolate between memoryless behavior \left(\alpha =1\right) and memory‐dependent behavior , in order to match long‐term epidemiological data.

2.11 Consistency, stability, and convergence of the fractional Euler scheme

In this subsection, we examine the numerical properties of the CF fractional Euler scheme for the TB model. Let the fractional TB model be written in compact form as

\begin{array}{}\text{(95)}& {}_0{}^{\text{CF}}{D}_t^{\alpha }X\left(t\right)=F\left(t,X\left(t\right)\right),X\left(0\right)=X_0,\end{array}

where X\left(t\right)=\left(S,E,I_{\mathrm{u}},I,T,R\right) and F represents the nonlinear transmission dynamics. The CF fractional Euler scheme is given by

\begin{array}{}\text{(96)}& X_{n+1}=X_n+\frac{1-\alpha }{M\left(\alpha \right)}F\left(t_n,X_n\right)+\frac{\alpha h}{M\left(\alpha \right)}\sum _{k=0}^{n}F\left(t_k,X_k\right),\end{array}

where h is the time step and M\left(\alpha \right) is the normalization function. Let’s say the function F meets the Lipschitz condition:

\begin{array}{}\text{(97)}& \Vert F\left(t,X_1\right)-F\left(t,X_2\right)\Vert \mathit{\le }L\Vert X_1-X_2\Vert .\end{array}

2.11.1 Consistency

Theorem 8.The Caputo–Fabrizio fractional Euler method is consistent with the Caputo–Fabrizio TB model, and the local truncation error \left({\tau }_n\right) meets the requirements of

\begin{array}{}\text{(98)}& {\tau }_n=O\left(h\right).\end{array}

Proof Consistency is established by comparing the exact solution with the numerical scheme. Using Taylor expansion of the exact solution,

\begin{array}{}\text{(99)}& X\left(t_{n+1}\right)=X\left(t_n\right)+h{X}^{\mathit{\prime }}\left(t_n\right)+O\left(h^2\right).\end{array}

From the definition of the CF fractional derivative,

\begin{array}{}\text{(100)}& X^{\mathit{\prime }}\left(t_n\right)=F\left(t_n,X_n\right).\end{array}

Substituting into the Taylor expansion gives

\begin{array}{}\text{(101)}& X\left(t_{n+1}\right)=X\left(t_n\right)+hF\left(t_n,X_n\right)+O\left(h^2\right).\end{array}

The numerical scheme gives

\begin{array}{}\text{(102)}& X_{n+1}=X_n+hF\left(t_n,X_n\right)+O\left(h\right).\end{array}

Thus, the local truncation error is

\begin{array}{}\text{(103)}& {\tau }_n=X\left(t_{n+1}\right)-X_{n+1}=O\left(h\right).\end{array}

Thus, the method is consistent of order one.

2.11.2 Stability

Theorem 9.If the function F satisfies the Lipschitz condition

\begin{array}{}\text{(104)}& \Vert F\left(t,X_1\right)-F\left(t,X_2\right)\Vert \mathit{\le }L\Vert X_1-X_2\Vert ,\end{array}

then the CF fractional Euler scheme is stable for sufficiently small step size h.

Proof Assume that two numerical solutions with different initial conditions are represented by X_n and {\tilde{X}}_n. Following that, the description of an error is

\begin{array}{}\text{(105)}& e_n=X_n-{\tilde{X}}_n.\end{array}

From the numerical scheme,

\begin{array}{}\text{(106)}& e_{n+1}=e_n+\frac{1-\alpha }{M\left(\alpha \right)}\left[F\left(t_n,X_n\right)-F\left(t_n,{\tilde{X}}_n\right)\right]+\frac{\alpha h}{M\left(\alpha \right)}\sum _{k=0}^n\left[F\left(t_k,X_k\right)-F\left(t_k,{\tilde{X}}_k\right)\right].\end{array}

Taking norms and using Lipschitz condition,

\begin{array}{}\text{(107)}& \Vert e_{n+1}\Vert \mathit{\le }\left(1+Ch\right)\Vert e_n\Vert .\end{array}

Applying this inequality recursively,

\begin{array}{}\text{(108)}& \Vert e_n\Vert \mathit{\le }\Vert e_0\Vert e^{C{t}_n}.\end{array}

Thus, the numerical solution remains bounded and the scheme is stable.

2.11.3 Convergence

Theorem 10.The Caputo–Fabrizio fractional Euler scheme converges to the exact solution of the TB model, and the global error satisfies

\begin{array}{}\text{(109)}& \underset{0\mathit{\le }n\mathit{\le }N}{\max }\Vert X\left(t_n\right)-X_n\Vert =O\left(h\right).\end{array}

Proof Define the global error

\begin{array}{}\text{(110)}& E_n=X\left(t_n\right)-X_n.\end{array}

Using consistency,

\begin{array}{}\text{(111)}& E_{n+1}=E_n+h\left[F\left(t_n,X\left(t_n\right)\right)-F\left(t_n,X_n\right)\right]+O\left(h^2\right).\end{array}

Using Lipschitz condition,

\begin{array}{}\text{(112)}& \Vert E_{n+1}\Vert =\left(1+Lh\right)\Vert E_n\Vert +C{h}^2.\end{array}

Applying discrete Gronwall inequality,

\begin{array}{}\text{(113)}& \Vert E_n\Vert \mathit{\le }Ch.\end{array}

Thus,

\begin{array}{}\text{(114)}& \underset{h\to 0}{\lim }\Vert X\left(t_n\right)-X_n\Vert =0.\end{array}

Therefore, the fractional Euler method converges to the exact solution.

2.12 Curve fitting and estimation of the parameters

To verify the effectiveness of the presented SE{I}_{\mathrm{u}}ITR fractional model, the TB incidence data for Ethiopia from 2000 to 2023 were utilized [38]. The parameters of the model were determined by fitting the model’s incidence curve to the real data using the MATLAB programming language. The fitting was conducted by minimizing the difference between the TB incidence data and the model’s output using a least squares method. The goodness of fit of the model was determined by computing the root mean square error, which is given by

\begin{array}{}\text{(115)}& \text{RMSE}=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(I_{\text{model}}\left(t_i\right)-I_{\text{data}}\left(t_i\right)\right)}^2,}\end{array}

where I_{\text{model}}\left(t_i\right) is the TB incidence predicted by the model at time t_i, I_{\text{data}}\left(t_i\right) is the real TB incidence data at time t_i, and n is the total number of data points. Figure 1 illustrates the comparison of the observed TB incidence and the result obtained from the model simulation. The results obtained demonstrate a good comparison and thus verify the model’s ability to simulate the TB transmission in the Ethiopian population. The estimation of the parameters ensures reliable quantitative estimates of the model’s key parameters associated with TB transmission and intervention strategies. The estimated model is then employed in the sensitivity analysis, phase‐plane analysis, and numerical simulations. The estimated model is then employed in projecting the TB incidence in the next 10 years under existing intervention strategies and thus emphasizes the need for enhanced intervention strategies to reduce TB prevalence in the region.

2.13 Sensitivity analysis on a local and global level

2.13.1 Local sensitivity analysis

Local sensitivity analysis is a measure of how a small change in each parameter affects the model [39]. The normalized forward sensitivity index for a parameter p with respect to R₀ is given by

\begin{array}{}\text{(116)}& {\mathrm{{\rm Y}}}_p^{R_0}=\frac{{\mathit{\partial }R}_0}{\mathit{\partial }p}\times \frac{p}{R_0}.\end{array}

Table 1 describes model parameters, values, and sources used in the TB transmission model. Estimated parameters were obtained by fitting the model to TB incidence data, while fixed parameters were adopted from standard epidemiological literature.

Local sensitivity indices of the basic reproduction number of R_0 are described in Table 2.

Global sensitivity analysis partial rank correlation coefficients (PRCC) of the basic reproduction number R_0 with respect to model parameters is described in Table 3. Baseline values are taken from MATLAB estimation, and parameter ranges are used for Latin hypercube sampling (LHS).

Figure 2 shows a plot of the local sensitivity indices for all parameters in the model. From the results, it is clear that the transmission rate parameter \beta and the modification parameter \theta are directly proportional to R_0, implying that a small change in these parameters results in a significant increase in disease transmission. Conversely, the parameter \omega associated with interventions is highly negative, which shows that it is a key parameter in controlling R_0.

2.13.2 Global sensitivity analysis

On the other hand, global sensitivity analysis is used to investigate the overall influence of each parameter across its entire range, considering possible interactions between the parameters. To obtain the global sensitivity indices, the LHS method with PRCC is used. Figure 3 shows the results for the global sensitivity analysis. From the results, it is evident that β and θ are the most important parameters in increasing R₀, while ω is the most sensitive parameter in reducing R_0. This shows that control strategies focusing on reducing transmission and changing behavior are important in controlling the spread of TB. From the above discussions, it is evident that the results from the local and global sensitivity analyses are important in controlling the spread of TB in a public health intervention context.

2.14 Numerical simulations of the population dynamics

The numerical simulations of the population dynamics are described in the subsequent figures (Figures 4–15). Moreover, the initial conditions and time of simulations are visible on the figures.

3 DISCUSSION

This section provides an interpretation of the results of the suggested mathematical model for fractional‐order TB transmission. The discussion of results entails the model fitting, sensitivity analysis, reproduction number determination, fractional order dynamics, and phase plane dynamics, where every figure is analyzed in light of the epidemiological model. Figure 1 depicts the data on the TB incidence in relation to the model predictions of incidence from 2000 to 2023. From the results, it can be noted that there is good agreement between the model predictions and the observed data, hence suggesting that the estimated parameters effectively model the transmission dynamics of the disease. The figure below provides an illustration of the local sensitivity indices for the basic reproduction number R_0 with respect to the model parameters. As can be seen from the results, \beta and \mathrm{\Lambda } have the highest positive indices which imply that they contribute significantly in increasing R_0 and promoting transmission of the infection.

On the other hand, parameters like \mu ,\delta , and \omega have negative sensitivity indices which imply that they contribute to the decrease in transmission and control of the disease. The figure below also represents the PRCC global sensitivity indices analysis for the basic reproduction number R_0. As expected from the local analysis, \beta and \mathrm{\Lambda } have very high positive correlations with R_0 while \mu ,\delta , and \omega exhibit strong negative correlations.

In Figure 4, we present the behavior of the susceptible individuals with respect to time. It can be seen that the number of susceptible is decreasing in time due to the propagation of disease. Fractional order leads to slower decrease as memory effect reduces the rate of infection. In Figure 5, the behavior of the exposed group is shown. Initially, the exposed group rises with an increase in infection rates. Then it starts declining with an advancement towards either infectious or another compartment. Smaller values for α result in more persistence of the exposed group, showing delayed progression of TB infection. Figure 6 demonstrates the dynamics of the unreported infectious I_{\mathrm{u}}. The obtained findings clearly demonstrate the effect of the fractional order on the magnitude of the unreported infections. It can be stated that a higher level of fractional order significantly decreases their occurrence and makes them appear later. Unreported infection plays a crucial role in TB transmission within the population. In Figure 7 the dynamics of active infection is presented.

The changes in the treatment population can be viewed in Figure 8. In view of the outcomes presented above, the number of treatments grows with increasing diagnosis and inclusion of patients into the process. However, a smaller fraction order leads to lower growth rates due to delays associated with the detection and treatment process. The change in the recovered population is depicted in Figure 9. As is seen from the outcomes provided above, recovery steadily increases because of treatment procedures performed. Recovery is characterized by delayed growth rates due to memory effects impacting disease development and treatment.

Figure 10 depicts the phase plane for susceptible S and infectious I populations. As can be seen from the outcome obtained, an increase in infection rate leads to a decrease in the susceptible population. This clearly indicates that the transmission of TB infection takes place as was expected from the outset. The relationship between the exposed E and infectious I populations is illustrated in Figure 11. As is evident, there is a smooth passage from latent infection to active disease state indicating stable and biological consistent dynamics. Phase plane between I and I_{\mathrm{u}} populations is shown in Figure 12. The simulation shows both continue increasing for a while and both decrease after a long period. Hence, an unaware and aware population should be considered to control the dynamics of TB infection in Ethiopia. Phase plane between T and I populations is shown in Figure 13. It is evident from the outcomes that an increase in T results in a decrease in I, thereby highlighting the role of the effectiveness of treatment intervention in disease control. Phase plane between R and I is provided in Figure 14. It is noted that an increase in recovery occurs when infection is minimized over time. Phase plane between R_0 and the parameters \beta and \theta is shown in Figure 15. It is observed that an increase in either of the two parameters causes an increase in R_0; hence the significance of these transmission parameters in disease persistence.

Finally, it is evident that the fractional‐order TB model proposed gives an effective way of modeling and studying TB dynamics. In addition, it has been observed that the use of memory effects in the model significantly enhances disease modeling, and sensitivity and phase plane analysis reveal the significant disease control parameters.

4 CONCLUSION

In this manuscript, a fractional‐order SE{I}_{\mathrm{u}}ITR mathematical model was proposed and used to explore the transmission of TB. The basic properties of the model, including the existence, positivity, and boundedness of solutions, were proved to guarantee the feasibility of the model. The basic reproduction number R_0 was calculated and used to establish the threshold condition for disease persistence and elimination. The model was successfully used to fit the TB incidence data from Ethiopia, and it was found that there was a good agreement between the model and real data, which confirmed the accuracy and reliability of the model. The sensitivity analysis was carried out and showed that the transmission rate \beta and the modification parameter \theta are directly proportional to R_0 and play a significant role in disease transmission, and the parameter \omega associated with intervention is the most sensitive parameter in reducing R_0 and thus is significant in disease control through strategies such as early diagnosis, treatment, and awareness. Additionally, the fractional‐order analysis showed that increasing the fractional‐order decreases the level of infection and improves stability, as a result of the memory effect included in the fractional derivative. Stability of the DFE was proved for , which confirms that TB can be controlled if appropriate intervention strategies are put in place. Future work can be carried out by adding other epidemiological parameters, which include vaccination, reinfection, and drug resistance, and developing optimal control strategies to determine how to control TB. Stochastic, spatial, and co‐infection factors, which include TB and lung cancer progression, can be included to make the model more realistic. This can be done in future work to make the model more applicable in supporting public health decisions and improving TB control strategies.

5 METHODS

The current fractional model TB transmission consists of

1.Susceptible (S): This compartment encloses TB free individuals with a chance to be infected in the future.

2.Exposed (E): This class includes persons who have been in contact with the TB bacterium but are not yet infectious and thus cannot spread the disease.

3.Undetected (unaware) infected (I_{\mathrm{u}}): This compartment represents persons who are infected with TB but are not yet aware of their status.

4.Detected (aware) infected (I): This class comprises persons who are TB positive but not on treatment.

5.On treatment (T): This compartment includes persons who are receiving treatment for TB infection.

6.Recovered (R): This class comprises persons recovered.

The total population size at time t, represented by N\left(t\right), is expressed as

\begin{array}{}\text{(117)}& N\left(t\right)=S\left(t\right)+E\left(t\right)+I_{\mathrm{u}}\left(t\right)+I\left(t\right)+T\left(t\right)+R\left(t\right).\end{array}

The explanation of the state variables and parameters used in the model development is presented in Tables 4 and 5, respectively.

In the construction of the schematic diagram of the TB dynamic transmission system, it is assumed that the disease induced death rate {\delta }_1={\delta }_2={\delta }_3=\delta. The fractional dynamics of each variable on Figure 16 TB model is given by

\begin{array}{}\text{(118)}& \begin{array}{l}{}_0{}^{\text{CF}}{D}_t^{\alpha }S\left(t\right)=\mathrm{\Lambda }-\beta S\left(t\right)\left(I_{\mathrm{u}}\left(t\right)+\theta I\left(t\right)\right)-\mu S\left(t\right),\\ {}_0{}^{\text{CF}}{D}_t^{\alpha }E\left(t\right)=\beta S\left(t\right)\left(I_{\mathrm{u}}\left(t\right)+\theta I\left(t\right)\right)-\left(\omega +\kappa +\mu \right)E\left(t\right),\\ {}_0{}^{\text{CF}}{D}_t^{\alpha }{I}_{\mathrm{u}}\left(t\right)=\kappa E\left(t\right)-\left(\varphi +\delta +\mu \right)I_{\mathrm{u}}\left(t\right),\\ {}_0{}^{\text{CF}}{D}_t^{\alpha }I\left(t\right)=\varphi I_{\mathrm{u}}\left(t\right)+\omega E\left(t\right)-\left(\rho +\delta +\mu \right)I\left(t\right),\\ {}_0{}^{\text{CF}}{D}_t^{\alpha }T\left(t\right)=\rho I\left(t\right)-\left(\eta +\delta +\mu \right)T\left(t\right),\\ {}_0{}^{\text{CF}}{D}_t^{\alpha }R\left(t\right)=\eta T\left(t\right)-\mu R\left(t\right),\end{array}\end{array}

with initial conditions

\begin{array}{}\text{(119)}& S\left(0\right)\ge 0,E\left(0\right)\mathit{\ge }0,I_{\mathrm{u}}\left(0\right)\mathit{\ge }0,I\left(0\right)\mathit{\ge }0,T\left(0\right)\mathit{\ge }0,R\left(0\right)\mathit{\ge }0.\end{array}

References

[1] Das DK, Khajanchi S, Kar TK. Influence of multiple re‐infections in tuberculosis transmission dynamics: a mathematical approach. In: 2019 8th international conference on modeling simulation and applied optimization (ICMSAO). IEEE; 2019. p. 1–5.

[2] Das DK, Khajanchi S, Kar TK. Transmission dynamics of tuberculosis with multiple re‐infections. Chaos Solitons Fractals. 2020;130:109450.

[3] Khajanchi S, Das DK, Kar TK. Dynamics of tuberculosis transmission with exogenous reinfections and endogenous reactivation. Phys Stat Mech Appl. 2018;497:52–71.

[4] Jabbari A, Lotfi M, Kheiri H, Khajanchi S. Mathematical analysis of the dynamics of a fractional‐order tuberculosis epidemic in a patchy environment under the influence of re‐infection. Math Methods Appl Sci. 2023;46(17):17798–17817.

[5] Ojo MM, Peter OJ, Goufo EFD, Panigoro HS, Oguntolu FA. Mathematical model for control of tuberculosis epidemiology. J Appl Math Comput. 2023;69(1):69–87.

[6] Ayinla AY, Othman WAM, Rabiu M. A mathematical model of the tuberculosis epidemic. Acta Biotheor. 2021;69(3):225–255.

[7] Kudryavtsev I, Starshinova A, Rubinstein A, Kulpina A, Ling H, Zhuang M, et al. The role of the immune response in developing tuberculosis infection: from latent infection to active tuberculosis. Front Tuberc. 2024;2:1438406.

[8] Omeje D, Acheneje G, Odiba P, Bolaji B. Modeling the transmission dynamics of the co‐infection of malaria and tuberculosis with optimal control strategies and cost‐benefit analysis. Research Square [Preprint]. 2024.

[9] Peter OJ, Abidemi A, Fatmawati F, Ojo MM, Oguntolu FA. Optimizing tuberculosis control: a comprehensive simulation of integrated interventions using a mathematical model. Math Model Numer Simul Appl. 2024;4(3):238–255.

[10] Ledder G. Mathematical modeling for epidemiology and ecology, Vol. 50. Berlin/Heidelberg: Springer; 2023. p. 223–276.

[11] Arenas A, Cota W, Gómez‐Gardeñes J, Gómez S, Granell C, Matamalas JT, et al. Modeling the spatiotemporal epidemic spreading of COVID‐19 and the impact of mobility and social distancing interventions. Phys Rev X. 2020;10(4):041055.

[12] Das DK, Khajanchi S, Kar TK. The impact of the media awareness and optimal strategy on the prevalence of tuberculosis. Appl Math Comput. 2020;366:124732.

[13] Jajarmi A, Baleanu D. A new fractional analysis on the interaction of HIV with CD4⁺ T‐cells. Chaos Solitons Fractals. 2018;113:221–229.

[14] Panwar VS, Uduman PS, Gómez‐Aguilar JF. Mathematical modeling of coronavirus disease COVID‐19 dynamics using CF and ABC non‐singular fractional derivatives. Chaos Solitons Fractals. 2021;145:110757.

[15] Peter OJ, Yusuf A, Oshinubi K, Oguntolu FA, Lawal JO, Abioye AI, et al. Fractional order of pneumococcal pneumonia infection model with Caputo Fabrizio operator. Results Phys. 2021;29:104581.

[16] Moore EJ, Sirisubtawee S, Koonprasert S. A Caputo–Fabrizio fractional differential equation model for HIV/AIDS with treatment compartment. Adv Differ Equ. 2019;2019(1):200.

[17] Unaegbu EN, Onah IS, Oyesanya MO. A fractional order HIV/AIDS model using Caputo‐Fabrizio operator. Afr J Infect Dis. 2021;15(2 Suppl):1–18.

[18] Ghosh S, Mondal J, Khajanchi S. Effect of media awareness in the spread of infectious diseases: S. Appl Math Comput. 2025;71(3):3833–3860.

[19] Sulayman F, Abdullah FA, Mohd MH. An SVEIRE model of tuberculosis to assess the effect of an imperfect vaccine and other exogenous factors. Mathematics. 2021;9(4):327.

[20] Malek A, Hoque A. Mathematical model of tuberculosis with seasonality, detection, and treatment. Inform Med Unlocked. 2024;49:101536.

[21] Farman M, Amilo D, Ghannam M, Nisar KS, Hafez M. Stability and optimizing the treatment control of tuberculosis model via numerical approach. Results Control Optim. 2026;22:100650.

[22] Hamam H, Awan AU, Medani M, Alroobaea R, Bukhari SAHS, Fathima D. Theoretical analysis of colonic crypt and colorectal cancer model through Caputo–Fabrizio fractional derivative. Mod Phys Lett B. 2024;38(36):2450365.

[23] Farman M, Akgül A, Ahmad A, Imtiaz S. Analysis and dynamical behavior of fractional‐order cancer model with vaccine strategy. Math Methods Appl Sci. 2020;43(7):4871–4882.

[24] Owolabi KM, Pindza E. A nonlinear epidemic model for tuberculosis with Caputo operator and fixed point theory. Healthc Anal. 2022;2:100111.

[25] Mollah S, Biswas S, Khajanchi S. Impact of awareness program on diabetes mellitus described by fractional‐order model solving by homotopy analysis method. Ricerche Matemat. 2024;73(5):2359–2384.

[26] Khajanchi S, Bera S, Kar TK. An optimal control problem for HTLV‐I infection model. Optim Control Appl Methods. 2025;46(2):798–810.

[27] Sardar M, Khajanchi S, Biswas S. Stochastic dynamics of a nonlinear tumor‐immune competitive system. Nonlinear Dyn. 2025;113(5):4395–4423.

[28] Malek A, Hoque A. Model development and prediction of COVID‐19 pandemic in Bangladesh with nonlinear incident. Iran J Sci. 2023;47(1):249–258.

[29] Peter OJ, Oguntolu FA, Ojo MM, Olayinka Oyeniyi A, Jan R, Khan I. Fractional order mathematical model of monkeypox transmission dynamics. Phys Scr. 2022;97(8):084005.

[30] Ullah S, Khan MA, Farooq M, Gul T. Modeling and analysis of tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan. Math Comput Simulat. 2019;165:181–199.

[31] Wojtak W, Silva CJ, Torres DF. Uniform asymptotic stability of a fractional tuberculosis model. Math Model Nat Phenom. 2018;13(1):9.

[32] van den Driessche P, Zou X. Modeling relapse in infectious diseases. Math Biosci. 2007;207(1):89–103.

[33] Van den Driessche P, Watmough J. Reproduction numbers and sub‐threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. 2002;180(1–2):29–48.

[34] Diekmann O, Heesterbeek H, Britton T. Mathematical tools for understanding infectious disease dynamics, Vol. 7. Princeton, NJ: Princeton University Press; 2013.

[35] Caputo M, Fabrizio M. A new definition of fractional derivative without singular kernel. Prog Fract Differ Appl. 2015;1(2):73–85.

[36] Diethelm K, Ford NJ. The analysis of fractional differential equations. J Math Anal Appl. 2002;265(2):229–248.

[37] Podlubny I. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Vol. 198. San Diego, California: Elsevier; 1998.

[38] World Bank. Tuberculosis incidence (per 100,000 people)–Ethiopia. Washington, D.C.: World Bank Group; 2024.

[39] Chitnis N, Hyman JM, Cushing JM. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull Math Biol. 2008;70(5):1272–1296.

[40] World Health Organization Staff. Global tuberculosis report 2013. World Health Organization; 2013.

Rights & permissions

2026 The Author(s). Quantitative Biology published by John Wiley & Sons Australia, Ltd on behalf of Higher Education Press.