[picture] | Computer Modeling in Engineering & Sciences | [picture]

DOI: 10.32604/cmes.2024.053681

ARTICLE

Modeling the Dynamics of Tuberculosis with Vaccination, Treatment, and Environmental Impact: Fractional Order Modeling

Muhammad Altaf Khan¹,*, Mahmoud H. DarAssi², Irfan Ahmad³, Noha Mohammad Seyam⁴ and Ebraheem Alzahrani⁵

1 Faculty of Natural and Agricultural Sciences, University of the Free State, Bleomfontein, 9300, South Africa
2 Department of Basic Sciences, Princess Sumaya University for Technology, Amman, 11941, Jordan
3 Department of Clinical Laboratory Sciences, College of Applied Medical Science, King Khalid University, Abha, 62529, Saudi Arabia
4 Mathematical Sciences Department, College of Applied Sciences, Umm Al-Qura University, Makkah, 24381, Saudi Arabia
5 Department of Mathematics, Faculty of Science, King Abdulaziz University, P.O. Box 80203, Jeddah, 21589, Saudi Arabia
Corresponding Author: Muhammad Altaf Khan. Email: altafdir@gmail.com
Received: 07 May 2024; Accepted: 18 July 2024; Published: 27 September 2024

Abstract: A mathematical model is designed to investigate Tuberculosis (TB) disease under the vaccination, treatment, and environmental impact with real cases. First, we introduce the model formulation in non-integer order derivative and then, extend the model into fractional order derivative. The fractional system’s existence, uniqueness, and other relevant properties are shown. Then, we study the stability analysis of the equilibrium points. The disease-free equilibrium (DFE) 𝒟0 is locally asymptotically stable (LAS) when ℛv<1. Further, we show the global asymptotical stability (GAS) of the endemic equilibrium (EE) 𝒟∗ for ℛv>1 and 𝒟0 for ℛv≤1. The existence of bifurcation analysis in the model is investigated, and it is shown the system possesses the forward bifurcation phenomenon. Sensitivity analysis has been performed to determine the sensitive parameters that impact ℛv. We consider the real TB statistics from Khyber Pakhtunkhwa in Pakistan and parameterized the model. The computed basic reproduction number obtained using the real cases is ℛ0≈3.6615. Various numerical results regarding disease elimination of the sensitive parameters are shown graphically.

Keywords: Tuberculosis; real data; stability analysis; parameter estimations; discussion

1  Introduction

A Mycobacterium tuberculosis (MTB) infection is the reason for tuberculosis (TB), a persistent contagious disease that frequently spreads through the air in the form of droplets [1]. Pulmonary tuberculosis is the name given to TB that mostly affects the lungs, however, it can affect any organ in the body [2]. Before the COVID-19 pandemic in 2019, TB ranked as the world’s 13th most common cause of mortality and the primary cause of death from single-source infections. Even though TB prevention and mitigation have come a long way over the last 20 years, 10.6 million new cases of the disease were reported globally in 2021. Based on estimates, Pakistan ranked fifth among B high-burden nations globally, accounting for 61% of the TB burden in the World Health Organization (WHO) Eastern Mediterranean Region, with 510,000 new cases of TB arising yearly and roughly 15,000 drug-resistant cases developing annually [3]. According to the WHO in 2002, approximately 1.3 million death cases have been reported, and 10.6 million infected cases are reported worldwide [4]. There is an urgent need to make more control and prevention measures to conduct research on TB and highlight this issue to make more efforts for its resolution.

Several recent mathematical studies have been conducted on TB dynamics in the past. For example, the authors in [5] introduced a mathematical formulation for TB dynamics and discussed the disease trends in India and other nations. The work in [6], provided the development of the age structure modeling of TB controlling under dots. Further, different case studies have investigated related to TB cures in different areas of the world. In [7], the authors conducted a study, to examine the TB dynamics incorporated by the features of slow and fast progression within the SEIR (Susceptible, Exposed, Infected, Recovered) model. The recurrence features in the TB model formulation have been done by [8], where the authors considered an XLTiTn R model and divided the treatment into infectious and non-infectious classes. They studied the time-dependent uncertainty and sensitivity analysis quantitatively, to understand the TB when the treatment effect is absent. The Bacillus Calmette-Guerin (BCG) for vaccinating susceptible people has been incorporated in the study [9], where the susceptible people are divided into unvaccinated and vaccinated individuals. Similarly, the exposed people with a history of no pre-exposure vaccination, and those with a history of pre-exposure vaccination are divided. They showed that treatment combined with vaccination is more effective than vaccination alone. The treatment and vaccination impact on the TB disease dynamics in terms of mathematical formulation has been completed in [10], where the global stability and optimal control results are obtained regarding the disease curtail in the population. A mathematical model for TB infection, considering the potential public health impacts of new vaccines and their predictions for countries with high incidence rates, is discussed in [11]. The model is subdivided into six groups: vaccinated uninfected, unvaccinated uninfected, unvaccinated latent, vaccinated latent, infected, and recovered classes. The features of complete and incomplete treatment in the transmission dynamics of TB have been formulated in terms of a mathematical model [12], where the SEIR model in terms of age structure has been studied. A mathematical model for TB with drug resistance by studying the high endemic countries of the Asia-Pacific area [13]. The authors considered various results regarding vaccine efficacy in partial and temporary, and reinfection during the latency period. A mathematical model in non-integer derivative describing the TB disease is discussed in [14], where the real data is considered and various results regarding the disease curtail are presented. TB and HIV coinfection models have been studied in [15] where the details analysis of each model and their co-dynamics are studied. The TB model is considered in the form of an SEIR type. The dynamics of TB with relapse and treatment under the age-structured model is analyzed in [16]. They considered the infected data of TB in China from 2017 to 2018 to obtain the optimal values of the parameters using the Grey Wolf Optimizer algorithm. A non-autonomous mathematical model for TB with seasonality and age structure has been studied in [17]. The findings show that immunizing vulnerable people over the age of sixty-five as well as those between the ages of 20 and 24 are significantly more efficient at lowering the overall incidence of TB and that each enhanced vaccination approach, screening approach, and therapeutic the approach results in significant decreases in the prevalence of TB per 100,000 people in comparison with present methods, and that combining all three methods is even more successful. A tuberculosis mathematical study with seasonality and without seasonality is discussed in [18]. They studied the model without seasonality and showed that it undergoes backward bifurcation while the model with seasonality provided the data fitting to the model. A mathematical model to study TB using the cases in China is given in [19] where the authors used the Chinese data on TB from 2005 to 2021. Further, they designed an optimal control problem using four controls to obtain the optimal control results. In [20], the authors used the control theory to study the middle eastern respiratory syndrome (MERS) coronavirus disease model.

Vaccination is considered a powerful tool for the prevention and control of diseases. The authors in [21] studied the COVID-19 spread in China with vaccination in some percentage. The authors in [22] studied the model with the Omicron strain and the Delta variant strain. To prevent TB disease, administer the BCG vaccine, which is made from a suspension of attenuated Bacillus bovis. It boosts macrophage action, enhances tumor cell killing, activates T lymphocytes, and improves cellular immunity. The two French scientists, Calmette and Guerin, are the first to develop it [23,24]. For a long time, it was considered that the vaccine BCG does not protect against the mycobacterium tuberculosis infection, however, it protects against progression of the active disease [25]. A vaccination model in [26] is formulated to study TB infection with relapse incomplete treatment, and slow and past progression.

Scientists and researchers are always exploring ways to find the best possible treatment against the disease for it curtail. Presently, the treatment available for TB is drug treatment. The TB treatment follows five principles and with many constraints, treatment for an individual frequently falls short of the intended outcomes. The incomplete treatment causes some of the infected people to have drug-resistant infections [19]. The latent individuals, as a control strategy to design an optimal control problem for TB, have been studied in [27]. Various studies indicate the survival of the TB bacteria outside the host in favorable conditions for a long time and even alcohol with 75% cannot be able to vanish it. Healthy individuals can easily be infected with the pathogens in the environment that is available to live freely. Thus, the concept of environmental study in TB modeling has been incorporated by many researchers [28]. The contaminated environment effect includes handkerchiefs, doorknobs, towels, and toys used by infected individuals. The authors conducted research on TB with environmental class in their mathematical modeling in China [29]. Some more related work in applied mathematics area are provided here in the work [30,31].

An expanded form of integer-order calculus is called fractional calculus. In the last few decades, it has been widely applied in various scientific and technological fields, such as electronic circuit analysis, control theory, heat transfer, and fluid dynamics [32]. The fractional derivative is primarily used because it incorporates memory, a feature shared by the majority of biological systems. Moreover, the fractional order derivative has a nonlocal characteristic. Additionally, the dynamical systems’ stability zone is expanded by fractional-order derivatives. It has been heavily utilized in the past several years to address numerous biological issues. The difficulties are modeled using differential equations of fractional order. In this connection, we may cite various research in which non-integer derivatives are considered to simulate and examine the transmission of different types of illnesses other than TB (see [33,34]). For example, a fractional study to examine TB disease has been considered in [35]. A mathematical model of TB in Caputo-Fabrizio derivative is explored in [36]. The COVID-19 mathematical study in different fractional kernels is given in [37]. A Zika virus model under fractional differential equation modeling is studied in [38]. There is more related research on fractional derivatives, see [39–42].

We design a mathematical model in the present analysis for tuberculosis infection with vaccination, environmental contamination, and therapies that significantly influence the TB disease’s dynamics. In the previous works mentioned above none of us used the concept of fractional derivative with real data of KP, in Pakistan to study TB, with vaccination, environmental impacts, and treatment. We considered a perfect vaccination and provided the results when there was an increase in vaccine and less waning rate of vaccine. We use real data and provide a more detailed analysis of the model with qualitative and numerical results. The rest of the paper has been organized in detail section-wise as follows: The basic definitions and the construction of the mathematical model are given in Section 2. The fundamental analysis of the model and its existence and uniqueness is shown in Section 3 while bifurcation and the local stability analysis are discussed in Section 4. Section 5studies the global stability analysis of the proposed model. Parameters estimation is shown in Section 6 while the detailed numerical results are given in Section 7. Results are summarized finally in Section 8.

2  Model Construction and Basic of Fractional Formulae

We shall provide some important formulae that will be used in the modeling of the problem. We follow the definitions from [32,43].

Definition 1. A Caputo derivative of order q of a function g∈Ck([0,∞],R), where q∈(k−1,k] is described by

Dqg(t)=1Γ(k−q)∫0t(t−ϕ)k−q−1g(k)(ϕ)dϕ,t>0.(1)

Definition 2. We can define the fractional Riemann–Liouville of the function g:R+→R shown by

Iqg(t)=1Γ(q)∫0t(t−ϕ)q−1g(ϕ)dϕ,(2)

where Γ(q) represents Euler Gamma function.

2.1 Model Construction

The present section formulates a mathematical model for tuberculosis infection, with treatment, vaccination, and environmental contamination. The model consists of seven classes where the population of humans is distributed among six different categories, such as susceptible people, S(t) (people who are vulnerable to TB infection but are not vaccinated), vaccinated people, V(t) (people who are vaccinated against the TB infection), exposed people, E(t) (people who are exposed to the disease but not yet infected), infected people, I(t) (people infected with TB disease), treated people, T(t) (people who are treated after being identified as TB infected), and those recovered people by R(t) (those who have recovered from the sickness). The compartment W(t) demonstrates the pathogen abundance in the polluted environment. The contaminated environment may include the bedding, toys, handkerchiefs, towels, etc. The total population of the human compartments can be represented by N(t)=S(t)+V(t)+E(t)+I(t)+T(t)+R(t). Currently, the BCG vaccine available in the market does not provide full prevention against mycobacterium tuberculosis, but it can give some protection to people who are vaccinated. As a result, in this research, we assume that inoculated individuals move to the exposed condition following infection. The TB model with vaccination and the combination of the environment contamination compartment takes the following shape in terms of the evolutionary differential equation:

{dSdt=Π−κ1SIN−κ2SWN−(μ+α)S+θV,dVdt=αS−μV−θV,dEdt=κ1SIN+κ2SWN−(μ+β+r1)E,dIdt=βE+ϕR−(μ+d+τ+r2)I,dTdt=τI−(r3+d+μ)T,dRdt=r1E+r2I+r3T−(ϕ+μ)R,dWdt=ψI−bW,(3)

where the initial conditions subject to the system (3) are given by

S(0)≥S0, V(0)≥V0, E(0)≥E0, I(0)≥I0, T(0)≥T0, R(0)≥R0, W(0)≥W0.(4)

In the TB system (3), the parameters Π and μ denote respectively the birth rate and natural mortality rate in each human class. The rate at which healthy humans are infected after close contact with the infected person by the rate κ1. The rate κ2 measures the contact of healthy people becoming infected through indirect exposure to the contaminated environment. Healthy people are being vaccinated through the rate α while the vaccine waning rate is given by θ. The exposed individuals become infected at a rate β, and the infected people are treated at a rate τ. The disease-related death of the TB infected and treated individuals is given by d. The recovery rate of the exposed, infected, and treated individuals is given by r1,r2, and r3, respectively. At a rate ϕ, the recovered people are infected again and join infected compartment I(t). The parameter ψ is the shedding rate of the virus from infected humans, and b is the clearance rate of pathogens in the environment.

{DtqS=Π−κ1SIN−κ2SWN−(μ+α)S+θV,DtqV=αS−μV−θV,DtqE=κ1SIN+κ2SWN−(μ+β+r1)E,DtqI=βE+ϕR−(μ+d+τ+r2)I,DtqT=τI−(r3+d+μ)T,DtqR=r1E+r2I+r3T−(ϕ+μ)R,DtqW=ψI−bW.(5)

3  Model Analysis

We shall carry out some important results regarding the fractional system (5) in the present portion.

Theorem 1. All the associated solutions to the system (5) are uniformly bounded and nonnegative.

Proof. It follows from the model (5) we get

DtqS|S=0 =Π+θV≥0,DtqV|V=0 =αS≥0,DtqE|E=0 =κ1SIN−E+κ2SWN−E≥0,DtqI|I=0 =βE+ϕR≥0,DtqT|T=0 =τI≥0,DtqR|R=0 =r1E+r2I+r3T≥0,DtqW|W=0 =ψI≥0.

So, it follows from the work given in [44] that all the solution will remain always in R+7. Assume that N(t)=S(t)+V(t)+E(t)+I(t)+T(t)+R(t), and adding the equations of the model (5) except the W equation, we have

DtqN(t) =Π−μN−d(I+T), ≤Π−μN.

Now consider the result from [45], we obtain

N(t)≤Πμ+(−Πμ+N(0))Eq(−μtq),

where Eq is the Mittag-Leffler function. So N(t)→Π/μ for t→∞, and so 0<N(t)≤Π/μ. We also use the W equation of the model (5) and obtain

DtqW=ψI−bW≤ψΠμ−bW,

DtqW+bW≤ψΠμ.

Again using the result from [45], we have

W(t)≤(W(0)−ψΠbμ)Eq(−btq)+ψΠbμ.

Hence, W(t)→ψΠ/bμ when t→∞, and thus 0<W(t)≤ψΠ/bμ. Therefore, all the solutions starting in R+7 are restricted in the region Θ1×Θ2, where

Θ1={(S,V,E,I,T,R)|0≤N≤Πμ},andΘ2={W|0≤W≤ψΠbμ}.

■

The existence of a unique solution is given in the theorem below.

Theorem 2. There exists a unique solution associated to the system (5).

Proof. Denote J(t)=(S,V,E,I,T,R,W)=(j1,j2,j3,j4,j5,j6,j7)T. Then, system (5) follows the form given by

DtqJ(t)=𝒬1J(t)+j1𝒬2J(t)+𝒬3,(6)

where

𝒬1=(−z1θ00000α−z20000000−z3000000β−z40ϕ0000τ−z50000r1r2r3−z60000ψ00−b),𝒬2=(000−κ1N00−κ2N0000000000κ1N00κ2N0000000000000000000000000000)

𝒬3=(Π000000).

Denoting Θ(t,J(t))=𝒬1J(t)+j1𝒬2J(t)+𝒬3. Now,

‖Θ(t,J(t))−Θ(t,P(t))‖ =‖(𝒬1J(t)+j1𝒬2J(t)+𝒬3)−(𝒬1P(t)+j1𝒬2P(t)+𝒬3)‖, =‖𝒬1(J(t)−P(t))+j1𝒬2(J(t)−P(t))‖, =‖(𝒬1+j1𝒬2)(J(t)−P(t))‖, ≤L‖J(t)−P(t)‖,

where L=max(𝒬1+j1𝒬2), and ∥.∥ represents the usual Euclidean norm. So, Θ(t,J(t)) holds Lipschitz condition thus, follows from [46], the fractional system (5) exists and has a unique solution. ■

4  Analysis of the Equilibrium Points

First, we obtain the disease-free equilibrium point of the model (5) in the following by solving the equations of the model (5) at a steady state by setting:

DtqS=0,DtqV=0,DtqE=0,DtqI=0,DtqT=0,DtqR=0,DtqW=0.

We get the following equilibrium point called the disease-free equilibrium is given by

𝒟0=(S0,V0,0,0,0,0,0),

where

S0=Π(θ+μ)μ(α+θ+μ), V0=αΠμ(α+θ+μ).

To study the analysis of the equilibrium points, first, we need to obtain the expression for the basic reproduction number ℛ0 using the concept given in [47]. According to [47], we have the following expressions:

ℱ=(κ1SIN+κ2SWN0000),𝒱=((μ+β+r1)E−βE−ϕR+(μ+d+τ+r2)I−τI+(r3+d+μ)T−r1E−r2I−r3T+z6R−ψI+bW).

Further, we have the following matrices:

F=(0κ1S0S0+V000κ2S0S0+V000000000000000000000),V=(z30000−βz40−ϕ00−τz500−r1−r2−r3z600−ψ00b).

Using the concept of spectral radius of FV−1, we obtain

ℛv =κ2r1ψϕS0bz3z4z6(S0+V0)+βκ2ψS0bz3z4(S0+V0)+κ1r1ϕS0z3z4z6(S0+V0)+βκ1S0z3z4(S0+V0)+r3τϕz4z5z6+r2ϕz4z6, =ℛ1+ℛ2+ℛ3+ℛ4+ℛ5+ℛ6,

where z1=(α+μ),z2=(θ+μ),z3=(β+μ+r1), z4=(d+μ+r2+τ),z5=(d+μ+r3), z6=(μ+ϕ). Now, we obtain the expression of the threshold quantity known as the basic reproduction number ℛ0, which can be derived from the vaccine-induced reproduction number ℛv, given by

ℛ0=(bκ1+κ2ψ)(r1ϕ+βz6)bz3z4z6+ϕ(r3τ+r2z5)z4z5z6.

We give the comparison of the basic reproduction number ℛ0 and ℛv in Fig. 1. We can see that vaccine reproduction number ℛv decreases the basic reproduction number ℛ0 effectively.

Figure 1: Comparison of ℛ0 vs. ℛv, values have been taken from Table 1

Endemic Equilibria

The finding of the expression of the endemic equilibria is important in disease models as it determines the number of possible equilibrium points in the presence of infection. Further, the existence of unique or multiple endemic equilibria determines the possibility of the bifurcation which may be of different types, such as forward or backward. We use the notation 𝒟∗ to denote the endemic equilibrium of the system (5) and obtain in the following:

𝒟∗=(S∗,V∗,E∗,I∗,T∗,R∗,W∗)

where

S∗ =Πz2−αθ+λ∗z2+z1z2,V∗ =αΠ−αθ+λ∗z2+z1z2,E∗ =λ∗S∗z3,I∗ =z3ϕR∗+βλS∗z3z4,T∗ =τI∗z5,R∗ =r1E∗+r2I∗+r3T∗z6,W∗ =ψI∗b.

Now, putting the above equations into

λ∗=κ1I∗+κ2W∗S∗+V∗+E∗+I∗+T∗+R∗

and then simplifying, we obtain the following:

AoI∗+A1=0,

where

A0 =bβz2(r3τ+r2z5)+br1z2(τϕ+z5(z4+ϕ))+bβz2z6(τ+z5)+bz2z4z5z6(1−ℛ5−ℛ6),A1 =bz3z4z5z6(α+z2)(1−ℛv).

Here, A0>0 obviously based on the fact ℛ5<ℛv, and ℛ6<ℛv. Also, ℛv<1 ensures that the coefficient A1>0. So, I∗=−A1/A0 ensures the existence of endemic equilibria depends on ℛv>1. Further, with the linear expression of the endemic equilibria, there is no possibility of backward bifurcation in the fractional system (5).

Theorem 3. The fractional system (5) at 𝒟0 is LAS if ℛ0<1.

Proof. At 𝒟0, we obtain the Jacobian matrix

J=(−z1θ0−κ1S0S0+V000−κ2S0S0+V0α−z20000000−z3κ1S0S0+V000κ2S0S0+V000β−z40ϕ0000τ−z50000r1r2r3−z60000ψ00−b).

The characteristics equation related to J is

λ7+f1λ6+f2λ5+f3λ4+f4λ3+f5λ2+f6λ+f7=0,(7)

where

f1 =b+z1+z2+z3+z4+z5+z6,f2 =z6(d+μ+τ)+μr2+μ(α+z2)+z3z4(1−ℛ4),f3 =b(z3+z4)z5+z6(bz5+z3(b+z5))+(z1+z2)z7+z4z6(b+z1+z2)(1−ℛ6) +μ(b+z3+z4+z5+z6)(α+z2)+bz3z4(1−ℛ2−ℛ4)+z3z4z6(1−ℛ3−ℛ4−ℛ6) +z3z4(z1+z2+z5)(1−ℛ4)+z4z5z6(1−ℛ5−ℛ6),f4 =z3z4z6{b(1−ℛ1−ℛ2−ℛ3−ℛ4−ℛ6)+z5(1−ℛ3−ℛ4−ℛ5−ℛ6)} +(−αθ+z2z5+z1(z2+z5))z3z4(1−ℛ4)+z4z5z6(b+z1+z2)(1−ℛ5−ℛ6) +(bz5+z4(b+z5)+z6(b+z5)+z3(b+z5+z6))(z1z2−αθ) +(z4z6(−αθ+z2(b+z3)+z1(b+z2+z3)))(1−ℛ6) +bz3z4((z2+z5)(1−ℛ2−ℛ4)+(1−ℛ2)z1) +b(z1+z2)(z3+z4)z5+z6(b(z1+z2)z3+z5(b(z1+z2)+z3(b+z1+z2))),f5 =bz3z4(z1+z2+z5)z6(1−ℛ1−ℛ2−ℛ3−ℛ4−ℛ6)+b(z1+z2)z3z4z5(1−ℛ2−ℛ4) +(1−ℛ4)z3z4z5(z1z2−αθ)+bz3z4(1−ℛ2−ℛ4)(z1z2−αθ)+(z1+z2)z3z5z6z4 ×(1−ℛ3−ℛ4−ℛ5−ℛ6)+z3z6z4(1−ℛ3−ℛ4−ℛ6)((z1z2−αθ)) +(bz5(z4+z6)+z3(bz6+z5(b+z6)))(z1z2−αθ)+z5z4z6(1−ℛ5−ℛ6) ×(b(z1+z2)+(z1z2−αθ))−br3τz3ϕ+b(1−ℛ6)z4z6(z1z2−αθ)+b(z1+z2)z5z6z3,f6 =−Abr2ϕ+Abz3z4z6(1−ℛ1−ℛ2−ℛ3−ℛ4)+Abz3z4z5(1−ℛ2−ℛ4) +Abz4z5z6(1−ℛ5−ℛ6)+Abz3z5z6+Az3z4z5z6(1−ℛ3−ℛ4−ℛ5−ℛ6) +b(z1+z2)z3z4z5z6(1−ℛ1−ℛ2−ℛ3−ℛ4−ℛ6),f7 =bz3z4z5z6(μ(α+z2))(1−ℛv),

where

z7=bz5+z4(b+z5)+z6(b+z5)+z3(b+z5+z6)

and z1z3−αθ>0. So all the coefficients f1,…,f7 are positive whenever ℛv<1. Further, the Routh-Hurwitz conditions can be implemented easily to ensure that the Jacobian matrix exhibits eigenvalues with negative real parts. Therefore, the fractional system of TB with vaccination is locally asymptotically stable at 𝒟0 whenever ℛv<1. ■

4.1 Backward Bifurcation Analysis

We present the backward bifurcation analysis of the model (5) when q=1 (as the bifurcation analysis is only studied at the the steady state of disease-free equilibrium (DFE)). It is important in epidemic models where the infection elimination occurs for ℛv<1. However, in this case, there exists yet another stable endemic equilibrium point for ℛv<1, indicating that disease elimination is dependent not solely on ℛv but also on the initial infection level. A Jacobian matrix related to the model (5) has an eigenvalue, zero at DFE 𝒟0 when ℛv=1. Here, we established the results for the existence of the backward bifurcation in the model (5) using the method presented in [48]. To do this, we assume that y1=S, y2=V, y3=E, y4=I, y5=T, y6=R, y7=W, then, the system (5) takes the form below:

{dy1dt=Π−κ1y1y4N−κ2y1y7N−z1y1+θy2=h1,dy2dt=αy1−z2y2=h2,dEdt=κ1y1y4N+κ2y1y7N−z3y3=h3,dy4dt=βy3+ϕy6−z4y4=h4,dy5dt=τy4−z5y5=h5,dy6dt=r1y3+r2y4+r3y5−z6y6=h6,dy7dt=ψy3−by7=h7,(8)

where N=y1+y2+y3+y4+y5+y6. Taking ℛv=1, then the solving for κ1 as a bifurcation parameter, we get κ1=κ1∗=−κ2ψb−z3(α+z2)(r3τϕ+z5(r2ϕ−z4z6))z2z5(r1ϕ+βz6). The Jacobian matrix of the system (8) at 𝒟0 when κ1=κ1∗ is given by

J∗=(−z1θ0−κ1∗S0S0+V000−κ2S0S0+V0α−z20000000−z3κ1∗S0S0+V000κ2S0S0+V000β−z40ϕ0000τ−z50000r1r2r3−z60000ψ00−b).

J∗ has an eigenvalue ‘zero’ while the other eigenvalues contain negative real parts. So, we can apply the method in [48] to our model (5). Next, we find that the matrix J∗ has right eigenvectors G=[g1,g2,…,g7]T and the left eigenvectors V=[v1,v2,…,v7], given by

g1 =g4(z2z3(r3τϕ+z5(r2ϕ−z4z6)))z5(z1z2−αθ)(r1ϕ+βz6), g2=−αg4(z3(r3τϕ+z5(r2ϕ−z4z6)))z5(αθ−z1z2)(r1ϕ+βz6),g3 =−g4(r3τϕ+z5(r2ϕ−z4z6))z5(r1ϕ+βz6), g4=g4>0,g5=g4τz5,g6 =g4(βr3τ+z5(βr2+r1z4))z5(r1ϕ+βz6), g7=g4ψb,

and

v1=v2=0,v3=v4(r1ϕ+βz6)z3z6, v4=v4>0, v5=r3v4ϕz5z6, v6=v4ϕz6, v7=κ2v4z2(r1ϕ+βz6)bz3z6(α+z2).

According to [48], the coefficients are

a¯=∑k,i,j=17vkgigj∂2hk∂xi∂xj(0,0)

b¯=∑i,k=17vkgi∂2hk∂xi∂κ1(0,0)

need to be computed using the second-order partial derivatives obtained at 𝒟0. Obtaining the partial derivatives and then inserting into the coefficients a¯, and b¯, we have

a¯=−2g42μv4z2(θ+z1)[l1(r3τ+r2z5)+l3z5]bΛz3z52z6(α+z2)2,

and

b¯=g4v3z2α+z2>0,

where l1=bκ1(z5(β−ϕ)−τϕ)+κ2ψz5(β−ϕ), l2=bκ1(β+z4)(τ+z5)+κ2ψ(βτ+z5(β+z4)), l3=r1(bκ1+κ2ψ)(τϕ+z5(z4+ϕ))+l2z6. We observe that b¯>0 and the coefficient a¯ is also negative that a¯<0, and this can be verified through the values of the parameters in Table 1, that is a¯=−3.006667∗10−8<0, which is a clear indication of the non-existence of the backward bifurcation in the system (5). According to the result in [48], there exists a forward bifurcation in the given system as shown in Fig. 2.

Figure 2: The plot describes the existence of the forward bifurcation in the model (3)

5  Global Stability

We use the Lemma in the proof of global stability.

Lemma 1. [49]. Consider that a continuous and derivable function f(t)∈R+. Then, for any time t, we have t≥t0,

Dtq(f(t)−f∗(t)−ln⁡f(t)f∗(t))≤1−f∗(t)f(t)Dtqf(t), f∗∈R+, q∈(0,1).(9)

Theorem 4. The fractional system (5) at 𝒟0 is globally asymptotical stable whenever ℛv≤1.

Proof. We consider the following Lyapunov function:

U(t)=g1E+g2I+g3T+g4R+g5W,(10)

where the positive constants g1,…,g5 might have values at a later time. Further, the time differentiation of the system (10) and then utilizing the equations of the system (5), we arrive at the following result:

DtqU(t)=g1DtqE+g2DtqI+g3DtqT+g4DtqR+g5DtqW.

DtqU(t) =g1[κ1SIN+κ2SWN−z3E]+g2[βE+ϕR−z4I]+g3[τI−z5T] +g4[r1E+r2I+r3T−z6R]+g5[ψI−bW], =[g2β−z3g1+g4r1]E+[g1κ1S0S0+V0−g2z4+g3τ+g4r2+g5ψ]I +[g4r3−g3z5]T+[g2ϕ−g4z6]R +[g1κ2S0S0+V0−g5b]W, ≤[g2β−z3g1+g4r1]E+[g1κ1S0S0+V0−g2z4+g3τ+g4r2+g5ψ]I +[g4r3−g3z5]T+[g2ϕ−g4z6]R+[g1κ2S0S0+V0−g5b]W.(11)

Now, consider the values assigned to the constants g1,…,g5, which are given by, g1=r1ϕ+βz6, g2=z3z6, g3=ϕr3z3z5, g4=z3ϕ and g5=κ2S0(r1ϕ+βz6)b(S0+V0). Now, using these values in the last equation of (11) and after simplifications, we obtain

DtqU≤z3z4z6 (ℛv−1)I.

It follows that DtqU=0 whenever ℛv=1, and DtqU if ℛv<1. The result follows that there exists singleton set 𝒟0 which is the maximal compact invariant set in {(S,V,E,I,T,R,W)∈Θ1×Θ2:DtqU=0}. Thus, the LaSalle Invariance Principle implies that every solution of the system (5) with (4) in Θ1×Θ2 tends to 𝒟0 whenever t→∞ for ℛ0≤1. ■

5.1 Global Stability Endemic Equilibrium (EE)

The following assumptions are made based on steady-state 𝒟∗ of the model (5): We assume S/N≤1, so S≤N.

{Π=κ1S∗I∗+κ2S∗W∗+z1S∗−θV∗,αS∗=z2V∗,κ1S∗I∗+κ2S∗W∗=z3E∗,βE∗+ϕR∗=z4I∗,τI∗=z5T∗,r1E∗+r2I∗+r3T∗=z6R∗,ψI∗=bW∗.

Theorem 5. The infected equilibrium 𝒟∗ is GAS if ℛv>1.

Proof. Let us consider the Lyapunov function in the structure given below:

K(t) =K1(S(t))+θz2K2V(t)+K3E(t)+κ1I∗S∗βE∗K4I(t)+κ1S∗I∗ϕR∗z5βE∗T∗K5(T(t)) +κ2S∗W∗r2I∗K6R(t)+κ2S∗bK7W(t),(12)

where

Kj=m(t)−m∗−m∗ln⁡m(t)m∗,

where j=1,…,7 and m=(S,V,E,I,T,R,W). Now taking the fractional derivative of K(t) and applying Lemma 1, we have

DtqK(t) ≤(1−S∗S)DtqS(t)+θz2(1−V∗V)DtqV(t)+(1−E∗E)DtqE(t) +κ1I∗S∗βE∗(1−I∗I)DtqI(t)+κ1S∗I∗ϕR∗z5βE∗T∗(1−T∗T)DtqT(t) +κ2S∗W∗r2I∗(1−R∗R)DtqR(t)+κ2S∗b(1−W∗W)DtqW(t).(13)

Computing the terms on the right side of Eq. (13) one by one below:

(1−S∗S)DtqS =(1−S∗S)[Π−κ1SI−κ2SW−z1S+θV], =(1−S∗S)[κ1S∗I∗+κ2S∗W∗+z1S∗−θV∗−κ1SI−κ2SW−z1S+θV], =κ1S∗I∗(1−S∗S−SIS∗I∗+II∗)+κ2S∗W∗(1−S∗S−SWS∗W∗+WW∗) +z1S∗(2−S∗S−SS∗)+θV∗(−1+S∗S+VV∗−VS∗SV∗), ≤κ1S∗I∗(1−S∗S−SIS∗I∗+II∗)+κ2S∗W∗(1−S∗S−SWS∗W∗+WW∗) +θV∗(−1+S∗S+VV∗−VS∗SV∗),

θz2(1−V∗V)DtqV =θz2(1−V∗V)[αS−z2V], =θ(1−V∗V)[V∗S∗S−V], =θV∗(1−VV∗+SS∗−SV∗VS∗),

(1−E∗E)DtqE =(1−E∗E)[κ1SI+κ2SW−z3E], =(1−E∗E)[κ1SI+κ2SW−(κ1S∗I∗+κ2S∗W∗)EE∗], =k1S∗I∗(1−EE∗+SIS∗I∗−SIE∗S∗I∗E)+κ2S∗W∗(1−EE∗+SWS∗W∗−SWE∗ES∗W∗),

κ1S∗I∗βE∗(1−I∗I)DtqI =κ1S∗I∗βE∗(1−I∗I)[βE+ϕR−z4I], =κ1S∗I∗βE∗(1−I∗I)[βE+ϕR−(βE∗+ϕR∗)I∗I], =κ1S∗I∗(1−II∗+EE∗−EI∗IE∗)+κ1S∗I∗βE∗ϕR∗(1−II∗+RR∗−RI∗IR∗),

κ1S∗I∗ϕR∗z5βE∗T∗(1−T∗T)DtqT =κ1S∗I∗ϕR∗z5βE∗T∗(1−T∗T)[τI−z5T], =κ1S∗I∗ϕR∗βE∗T∗(1−T∗T)[T∗TI−T], =κ1S∗I∗ϕR∗βE∗(1−TT∗+II∗−IT∗TI∗),

κ2S∗W∗r2I∗(1−R∗R)DtqR =κ2S∗W∗r2I∗(1−R∗R)[r1E+r2I+r3T−z6R], =κ2S∗W∗r2I∗(1−R∗R)[r1E+r2I+r3T−r1E∗R∗R−r2I∗R∗R−r3T∗R∗R], =κ2S∗W∗r2I∗r1E∗(1−RR∗+EE∗−ER∗RE∗)+κ2S∗W∗(1−RR∗+II∗−IR∗RI∗) +κ2S∗W∗r2I∗r3T∗(1−RR∗+TT∗−TR∗RT∗),

κ2S∗b(1−W∗W)DtqW =κ2S∗b(1−W∗W)[ψI−bW], =κ2S∗(1−W∗W)[II∗W∗−W], =κ2S∗W∗(1−WW∗+II∗−IW∗WI∗).

We submit the above into Eq. (12) and simplifying, we obtain,

DtqK(t) =κ1S∗I∗(3−S∗S−EI∗IE∗−SIE∗S∗I∗E)+κ2S∗W∗(4−S∗S−RR∗−EE∗−IW∗WI∗−IR∗RI∗−SWE∗ES∗W∗) +κ1S∗I∗ϕR∗βE∗(2−TT∗−IT∗TI∗+RR∗−RI∗IR∗)+θV∗(S∗S−VS∗SV∗+SS∗−SV∗VS∗) +κ2S∗W∗r2I∗r1E∗(1−RR∗+EE∗(1−R∗R))+κ2S∗W∗r2I∗r3T∗(1−RR∗+TT∗(1−R∗R)).

Here, the arithmetic mean is greater or equal to the geometric mean, so we get,

(3−S∗S−EI∗IE∗−SIE∗S∗I∗E)≤0,(4−S∗S−RR∗−EE∗−IW∗WI∗−IR∗RI∗−SWE∗ES∗W∗)≤0.

If (2−TT∗−IT∗TI∗+RR∗−RI∗IR∗)≤0, (S∗S−VS∗SV∗+SS∗−SV∗VS∗)≤0, (1−RR∗+EE∗(1−R∗R))≤0, (1−RR∗+TT∗(1−R∗R))≤0, then the DtqK(t)≤0. Thus, the system (5) is globally asymptotically stable if ℛv>1.

■

6  Estimation of Parameters

The estimation of the parameters of the model (5) will be investigated for q=1 in the absence of vaccination using the cases reported TB in Khyber Pakhtunkhwa in Pakistan for the period 2002–2017 [50]. We utilize the nonlinear least square fitting method to estimate the model parameters where some of the parameters as the birth and recruitment rates are considered from literature or estimated from existing equations of the considered model while the parameters other than birth and death are fitted to the model. The life expectancy of an individual in Pakistan was 67.7 years in 2017. So, we write the natural death rate can be obtained as the inverse of life expectancy given by 1/67.7 years [51]. In 2002, the total population of Khyber Pakhtunkhwa was 30,523,371, so, we assume the initial size of the population is N(0)=30,523,371 [52]. We assume that the total population of Pakistan from 2002–2017 will not change significantly. Therefore, in our numerical simulation, the total population N(t) is assumed to be constant, i.e., Π=μ∗N(0) which describes the birth of the susceptible people. According to the reported data, the initial cases in 2002 were reported as 8010, so we consider I(0)=8010 as the initial TB-infected cases. The remaining initial values of the variables involved in the model are computed as S(0)=N(0)−E(0)−8010−T(0)−R(0)=30,075,361, where we assume that there is no treatment and recovery from infection on the disease starting period, that is T(0)=R(0)=0, while E(0)=440,000 has been adjusted to the data fitting. The population of the contaminated environment is subjected to the data fitting is W(0)=100. The method is explained in the following follows [53] to get the results:

i)   Create the objective function that measures how much the real and projected values vary from one another. The stated designed function is the sum of squared residuals between the model’s forecast and the corresponding real data cases.

ii)   First, provide the initial values of the parameters that need estimations. It might be an initial guess or a prior knowledge base.

iii)   Using the initial values of the parameters shown in step (ii), the model is simulated to obtain the predictions of the model.

iv)   Applying the “lsqcurvefit” optimization technique to see whether the objective function can be minimized.

v)   Setting convergence criteria to end the process of iterative estimating. It is based on reaching a predetermined threshold for the relevant objective function, an appropriate number of iterations, or the smallest modification in parameter estimations.

vi)   Set the original parameter estimations back to repeat steps (i–v) until an acceptable agreement across the model simulation and the actual data is attained if the estimated parameter values do not fulfill the convergence condition or do not match the real data curve sufficiently.

The model in the absence of vaccination provides a reasonable fitting which is shown in Fig. 3 while the obtained fitted/estimated parameters are shown briefly in Table 1. For the parameter values listed in Table 1, we computed ℛ0, which is ℛ0≈3.6615 while vaccine basic reproduction number is ℛv≈2.83209.

Figure 3: Model fitting with real data. The circle shows real TB data while the bold line is the model solution. (a) model vs. data fitting (b) data vs. model long-time predictions (c) data vs. fractional model q=0.99

6.1 Sensitivity Analysis

6.1.1 Local Sensitivity Analysis

Here, we study the sensitivity analysis of the parameters involved in ℛ0. Sensitivity analysis is important in epidemic modeling because it provides information about the parameters that how much they are sensitive to ℛ0. Such sensitive parameters are then considered in focus for the eradication of disease as a control. We first study the local sensitivity analysis of the parameters in ℛ0. The formula in [54] that used to get the sensibility of the parameters involved in ℛ0 is given by

∏cℛ0=∂ℛ0∂c×cℛ0,(14)

where c is a general parameter of ℛ0. Using the the formula above, we obtained the sensitivity indices, which are shown in Table 2. We see from Table 2 that the sensitive parameters are κ1, μ, r2, etc., that enhance the basic reproduction number. We give some parameters in terms of ℛ0 and show it graphically in Figs. 4 and 5.

Figure 4: The plot represents the parameters κ1 and τ in terms of ℛ0 (a) 3D contour (b) 2D contour

Figure 5: The plot represents the parameters κ1 and r2 in terms of ℛ0 (a) 3D contour (b) 2D contour

6.1.2 Global Sensitivity Analysis

This section studies the sensitivity and uncertainty of the basic reproduction ℛ0 in relation to each parameter, which is given in Fig. 6. We do this by using the Partial rank correlation coefficient (PRCC) technique and the Latin hypercube sampling (LHS) technique [55]. It follows from the Fig. 6 that the positive correlated parameters with ℛ0 are respectively, κ1, κ2, β, ϕ, r1, r2, r3 and ψ while the parameters that are negative correlated are μ, τ, b and d. The absolute value of the PRCC greater or equal to 0.5 represents a high correlation between the output and input parameters, we can see a moderate correlation when the value lies in the range of 0.2 and 0.4 while a weak correlation exists for the range 0 to 0.2. Here, the parameters μ, b, d, κ2, β, ϕ, ψ are correlated highly with the basic reproduction number ℛ0.

Figure 6: The plot indicates the sensitivity of ℛ0 for there parameter

7  Numerical Results

We perform the numerical simulation in the present section by considering the parameter values shown in Table 1. The model (5) is solved numerically by the method explained in [46,56]. According to the nature of the data, we consider the units in simulation per year. We present graphical results for the fractional model (5) under different fractional orders. Figs. 7 and 8 show the dynamics of the model using different fractional order q. We note, that the solutions converge to the DFE equilibrium when the order q decreases.

Figure 7: Numerical simulation of the healthy, exposed, and infected compartments (a–c) respectively for different fractional order q

Figure 8: Numerical simulation of the treated, recovered and contaminated environmental compartments (a–c) respectively for different fractional order q

We see from Fig. 9 the impact of the κ1 on the dynamics of exposed, infected and environmental populations. By decreasing the contact between healthy and infected populations, the number of future cases in the community decreases. Adequate ventilation and natural light are essential safety measures to reduce the risk of TB infection. Without proper ventilation, TB can remain in the air for hours. Ultraviolet (UV) light kills bacteria, so exposure to natural light helps eliminate the bacteria. It is also important to maintain good hygiene by covering the mouth and nose when sneezing or coughing to prevent the transmission of TB bacteria.

Figure 9: The plot shows the dynamics of exposed, infected and environmental population for various values of κ1 when q=0.97. (a) exposed population, (b) infected population, and (c) environmentally contaminated population

Fig. 10 represents the comparison of vaccination vs. no vaccination on the dynamics of exposed, infected, and environmental populations. It can be seen that the vaccine impacts greatly disease dynamics and reduces future cases in society. When increasing the vaccine rate α=0.02 and fixing the waning rate of vaccine θ=0.01, there is a clear difference between the two models with and without vaccination, see Fig. 11. The vaccination impact on the exposed and infected population by varying α, see the result in Fig. 12. Increasing the vaccination to individuals will decrease the number of cases in the future in the country. Infants in Pakistan are immunized with the BCG vaccination, which can prevent serious TB in children. The BCG vaccination, which prevents tuberculosis has been used extensively since 1921. In Pakistan, newborns are usually given a single dose of the BCG vaccination.

Figure 10: The plot shows the dynamics of exposed, infected and environmental population with and without vaccination. (a) to (c) represent the exposed, infected, and environmental population

Figure 11: The plot shows the dynamics of exposed and infected population for α=0.01 and θ=0.01 Subfigures (a) to (b) represent the exposed and the infected population

Figure 12: The plot shows the dynamics of exposed and infected populations with the variation in α. (a) to (b) represent the exposed and the infected population

We show the results for the treatment parameter τ on the exposed and infected population, see Figs. 13 and 14. We can see that by increasing the treatment for the infected and exposed populations, there is a decrease in the population of exposed and infected individuals. The TB program in Pakistan provides 1700 TB care institutions in both commercial and public sectors of the nation. In 2002, more than 424,000 infected persons with TB were treated successfully. If this process continues, a significant decrease in TB cases will be observed in the coming years.

Figure 13: The treatment impact on the dynamics of exposed and infected population when varying τ, q=1. (a) exposed population (b) infected population

Figure 14: The treatment impact on the dynamics of exposed and infected population when varying τ, q=0.93. (a) exposed population (b) infected population

8  Conclusion

In this study, we investigated tuberculosis infection with treatment, vaccination, and environmental impact using non-integer order derivatives. The basic modeling is done initially in integer case and later the model is extended to non-integer order based on the definition of Caputo. First, the formulation of the model has been obtained and discussed in detail. Then, we extended the model into a fractional order differential equation using the concept of Caputo derivative. The existence and uniqueness of the fractional model have been obtained. The stability of the equilibrium points is obtained and discussed based on ℛ0. We found that the disease-free equilibrium 𝒟0 is locally asymptotically stable for ℛv<1. The equilibrium point 𝒟0 is stable globally asymptotically when ℛv≤1. We proved that there exists a unique endemic equilibrium when ℛv>1. The existence of the unique endemic equilibrium ensures the impossibility of the backward bifurcation existence in the model where the DFE coexists with EE. We also proved the impossibility of the backward bifurcation in the given model, hence, there exists only the forward bifurcation analysis and the only focus to reduce the cases is ℛ0. The non-existence of the backward bifurcation suggests that the infection cases can be minimized by focusing on the sensitive parameters involved in ℛ0. The model is stable globally asymptotically when ℛv>1. The sensitivity analysis is performed for the basic reproduction number and for obtaining the results graphically. We considered the local and global sensitivity analysis and determined the most sensitive parameters that impact the basic reproduction number.

We considered the real data of TB in Khyber Pakhtunkhwa in Pakistan and experimented using the nonlinear least square method. The estimated parameters obtained from the experiment that provided reasonable fitting to the data have been used in numerical simulation and obtained graphical results regarding disease elimination. We observed that vaccinating more and more individuals and reducing the contact between healthy and infected people can minimize well the future cases of TB in the country. The treatment provides important results regarding the disease control of TB in the country. Increasing the treatment for the long-term program, the number of TB cases has decreased significantly in the country. Vaccination impact on the disease has a significant impact on the TB cases in Pakistan. The results obtained in the numerical section show that vaccination has a great impact on the reduction of TB cases. So, vaccination, treatment, and controlling environmental contamination will significantly decrease TB cases in the country shortly.

The results obtained in the work are based on the formulation of the model with their parameters and the available data. The time domain considered for the data is per year which may have different impacts with other time units. The model results rely on the real data used and the accuracy is limited to the quality and quantitative data which is different from the other data used in the model. The demographic parameters may vary from region to region and thus will have different results. The environmental factors are also different from region to region and hence will impact the result.

Acknowledgement: The authors are thankful to the Deanship of Research and Graduate Studies, King Khalid University, Abha, Saudi Arabia, for financially supporting this work through the Large Research Group Project under Grant No. R.G.P.2/507/45.

Funding Statement: This work is supported by King Khalid University.

Author Contributions: Muhammad Altaf Khan, Mahmoud H. DarAssi, Irfan Ahmad; data collection: Noha Mohammad Seyam; Ebraheem Alzahrani; Muhammad Altaf Khan, analysis and interpretation of results: Muhammad Altaf Khan, Mahmoud H. DarAssi and Ebraheem Alzahrani; draft manuscript preparation: Muhammad Altaf Khan, Mahmoud H. DarAssi. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: The data used in the manuscript is available from the corresponding author on a reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Tuberculosis; 2023. Available from: https://www.who.int/news-room/fact-sheets/detail/tuberculosis. [Accessed 2023].

2. George M. The challenge of culturally competent health care: applications for asthma. Heart & Lung. 2001;30(5):392–400. doi:10.1067/mhl.2001.118364.

3. Pakistan. Tuberculosis. Available from: https://www.emro.who.int/pak/programmes/stop-tuberculosis.html. [Accessed 2023].

4. World Health Organization. Available from: https://www.who.int/en/news-room/factsheets/detail/tuberculosis. [Accessed 2023].

5. Waaler H, Geser A, Andersen S. The use of mathematical models in the study of the epidemiology of tuberculosis. Am J Public Health and the Nations Health. 1962;52(6):1002–13. doi:10.2105/AJPH.52.6.1002.

6. Dye C, Garnett GP, Sleeman K, Williams BG. Prospects for worldwide tuberculosis control under the WHO DOTS strategy. Lancet. 1998;352(9144):1886–91. doi:10.1016/S0140-6736(98)03199-7.

7. Song B, Castillo-Chavez C, Aparicio JP. Tuberculosis models with fast and slow dynamics: the role of close and casual contacts. Math Bio. 2002;180(1–2):187–205. doi:10.1016/S0025-5564(02)00112-8.

8. Porco TC, Blower SM. Quantifying the intrinsic transmission dynamics of tuberculosis. Theor Popul Biol. 1998;54(2):117–32. doi:10.1006/tpbi.1998.1366.

9. Bhunu C, Garira W, Mukandavire Z, Magombedze G. Modelling the effects of pre-exposure and post-exposure vaccines in tuberculosis control. J Theor Biol. 2008;254(3):633–49. doi:10.1016/j.jtbi.2008.06.023.

10. Yang Y, Tang S, Xlaohong R, Zhao H, Guo C. Global stability and optimal control for a tuberculosis model with vaccination and treatment. Discrete Conti Dyn Syst Ser B. 2016;21(3):1–26. doi:10.3934/dcdsb.

11. Ziv E, Daley CL, Blower S. Potential public health impact of new tuberculosis vaccines. Emerg Infect Dis. 2004;10(9):1529–35. doi:10.3201/eid1009.030921.

12. Ren S. Global stability in a tuberculosis model of imperfect treatment with age-dependent latency and relapse. Math Biosci Eng. 2017;14(5&6):1337–60. doi:10.3934/mbe.2017069.

13. Trauer JM, Denholm JT, McBryde ES. Construction of a mathematical model for tuberculosis transmission in highly endemic regions of the Asia-Pacific. J Theor Biol. 2014;358(33):74–84. doi:10.1016/j.jtbi.2014.05.023.

14. Khan MA, Ullah S, Farooq M. A new fractional model for tuberculosis with relapse via Atangana-Baleanu derivative. Chaos Soliton Fract. 2018;116(2):227–38. doi:10.1016/j.chaos.2018.09.039.

15. Bhunu C, Garira W, Mukandavire Z. Modeling HIV/AIDS and tuberculosis coinfection. Bull Math Biol. 2009;71(7):1745–80. doi:10.1007/s11538-009-9423-9.

16. Guo ZK, Xiang H, Huo HF. Analysis of an age-structured tuberculosis model with treatment and relapse. J Math Biol. 2021;82(5):1–37. doi:10.1007/s00285-021-01595-1.

17. Xue L, Jing S, Wang H. Evaluating strategies for tuberculosis to achieve the goals of WHO in China: a seasonal age-structured model study. Bull Math Biol. 2022;84(6):61. doi:10.1007/s11538-022-01019-1.

18. Bowong S, Kurths J. Modeling and analysis of the transmission dynamics of tuberculosis without and with seasonality. Nonlin Dyn. 2012;67(3):2027–51. doi:10.1007/s11071-011-0127-y.

19. Kang TL, Huo HF, Xiang H. Dynamics and optimal control of tuberculosis model with the combined effects of vaccination, treatment and contaminated environments. Math Biosci Eng. 2024;21(4):5308–34. doi:10.3934/mbe.2024234.

20. Fatima B, Yavuz M, ur Rahman M, Al-Duais FS. Modeling the epidemic trend of middle eastern respiratory syndrome coronavirus with optimal control. Math Biosci Eng. 2023;20(7):11847–74. doi:10.3934/mbe.2023527.

21. Li T, Guo Y. Modeling and optimal control of mutated COVID-19 (Delta strain) with imperfect vaccination. Chaos Soliton Fract. 2022;156(4):111825. doi:10.1016/j.chaos.2022.111825.

22. Guo Y, Li T. Modeling the competitive transmission of the Omicron strain and Delta strain of COVID-19. J Math Anal Appl. 2023;526(2):127283. doi:10.1016/j.jmaa.2023.127283.

23. Hawn TR, Day TA, Scriba TJ, Hatherill M, Hanekom WA, Evans TG, et al. Tuberculosis vaccines and prevention of infection. Microb Mol Biol Rev. 2014;78(4):650–71. doi:10.1128/MMBR.00021-14.

24. Bin Sayeed MS, James SL, Abate D, Abate KH, Abay SM, Abbafati C, et al. Global, regional, and national incidence, prevalence, andyears lived with disability for 354 diseases and injuries for 195 countries and territories, 1990–2017: a systematicanalysis for the Global Burden of Disease Study 2017. Global Health Metrics. 2020;396(10258):1204–22. doi:10.1016/S0140-6736(18)32279-7.

25. White PJ, Garnett GP. Mathematical modelling of the epidemiology of tuberculosis. In: Modelling parasite transmission and control. 2010. vol. 673, p. 127–40. doi:10.1007/978-1-4419-6064-1_9.

26. Li Y, Liu X, Yuan Y, Li J, Wang L. Global analysis of tuberculosis dynamical model and optimal control strategies based on case data in the United States. Appl Math Comput. 2022;422:126983.

27. Jiang Q, Liu Z, Wang L, Tan R. A tuberculosis model with early and late latency, imperfect vaccination, and relapse: an application to China. Math Methods Appl Sci. 2023;46(9):10929–46.

28. Xu A, Wen ZX, Wang Y, Wang WB. Prediction of different interventions on the burden of drug-resistant tuberculosis in China: a dynamic modelling study. J Glob Antimicrob Resist. 2022;29:323–30.

29. Cai Y, Zhao S, Niu Y, Peng Z, Wang K, He D, et al. Modelling the effects of the contaminated environments on tuberculosis in Jiangsu. China J Theor Biol. 2021;508:110453.

30. Sk N, Mondal B, Thirthar AA, Alqudah MA, Abdeljawad T. Bistability and tristability in a deterministic prey-predator model: transitions and emergent patterns in its stochastic counterpart. Chaos Soliton Fract. 2023;176:114073.

31. Nyaberi H, Mutuku W, Malonza D, Gachigua G, Alworah G. An optimal control model for Coffee Berry Disease and Coffee Leaf Rust co-infection. J Math Anal Model. 2024;5(1):1–25.

32. Podlubny I. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. 1st ed. Elsevier; 1998 Oct 21. vol. 198.

33. Majee S, Jana S, Das DK, Kar T. Global dynamics of a fractional-order HFMD model incorporating optimal treatment and stochastic stability. Chaos Soliton Fract. 2022;161:112291.

34. Naik PA, Zu J, Owolabi KM. Global dynamics of a fractional order model for the transmission of HIV epidemic with optimal control. Chaos Soliton Fract. 2020;138:109826.

35. Ullah S, Khan MA, Farooq M. A fractional model for the dynamics of TB virus. Chaos Soliton Fract. 2018;116:63–71.

36. Ullah S, Khan MA, Farooq M, Hammouch Z, Baleanu D. A fractional model for the dynamics of tuberculosis infection using Caputo-Fabrizio derivative, Discrete Cont Dyn-S, 2020;13(3):975–93. doi:10.3934/dcdss.2020057.

37. Kumar S, Chauhan R, Momani S, Hadid S. Numerical investigations on COVID-19 model through singular and non-singular fractional operators. Numer Methods Partial Diff Equ. 2024;40(1):e22707. doi:10.1002/num.22707.

38. Zafar ZUA, Khan MA, Inc M, Akgül A, Asiri M, Riaz MB, et al. The analysis of a new fractional model to the Zika virus infection with mutant. Heliyon. 2024;10(1):e23390. doi:10.1016/j.heliyon.2023.e23390.

39. Manikandan S, Gunasekar T, Kouidere A, Venkatesan K, Shah K, Abdeljawad T. Mathematical modelling of HIV/AIDS treatment using caputo-fabrizio fractional differential systems. Qual Theory Dyn Syst. 2024;23(4):1–29. doi:10.1007/s12346-024-01005-z.

40. Kalra P, Malhotra N. Modeling and analysis of fractional order logistic equation incorporating additive allee effect. Contemp Math. 2024;380–401. doi:10.37256/cm.5120243183.

41. Alla Hamou A, Azroul E, Bouda S, Guedda M. Mathematical modeling of HIV transmission in a heterosexual population: incorporating memory conservation. Model Earth Syst Environ. 2024;10(1):393–416. doi:10.1007/s40808-023-01791-6.

42. Paul S, Mahata A, Karak M, Mukherjee S, Biswas S, Roy B. A fractal-fractional order Susceptible-Exposed-Infected-Recovered (SEIR) model with Caputo sense. Healthc Anal. 2024;5(14):100317. doi:10.1016/j.health.2024.100317.

43. Almeida R, Brito da Cruz AM, Martins N, Monteiro MTT. An epidemiological MSEIR model described by the Caputo fractional derivative. Int J Dyn Control. 2019;7(2):776–84. doi:10.1007/s40435-018-0492-1.

44. Odibat ZM, Shawagfeh NT. Generalized Taylors formula. Appl Math Comput. 2007;186(1):286–93.

45. Li HL, Zhang L, Hu C, Jiang YL, Teng Z. Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge. J Appl Math Comput. 2017;54(1–2):435–49. doi:10.1007/s12190-016-1017-8.

46. Diethelm K, Ford NJ. Analysis of fractional differential equations. J Math Anal Appl. 2002;265(2):229–48. doi:10.1006/jmaa.2000.7194.

47. 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. doi:10.1016/S0025-5564(02)00108-6.

48. Castillo-Chavez C, Song B. Dynamical models of tuberculosis and their applications. Math Biosci Eng. 2004;1(2):361–404. doi:10.3934/mbe.2004.1.361.

49. Vargas-De-León C. Volterra-type Lyapunov functions for fractional-order epidemic systems. Comm Nonlin Sci Numer Simul. 2015;24(1–3):75–85. doi:10.1016/j.cnsns.2014.12.013.

50. National TB Control Program Pakistan (NTP). Available from: http://www.ntp.gov.pk/national_data.php. [Accessed 2023].

51. Pakistan Bureau of Statistics. Pakistans 6th census: population of major cities census. 2017. Available from: http://www.pbscensus.gov.pk/sites/default/files/population_of_major_cities_census_2017.pdf. [Accessed 2023].

52. World Health Organization. Who country cooperation strategic. 2016. Available from: http://apps.who.int/iris/bitstream/10665/136607/1/ccsbrief_pak_en.pdf. [Accessed 2023].

53. Pal KK, Sk N, Rai RK, Tiwari PK. Examining the impact of incentives and vaccination on COVID-19 control in India: addressing environmental contamination and seasonal dynamics. Eur Phy J Plus. 2024;139(3):225. doi:10.1140/epjp/s13360-024-04997-4.

54. 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–96. doi:10.1007/s11538-008-9299-0.

55. Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theor Biol. 2008;254(1):178–96. doi:10.1016/j.jtbi.2008.04.011.

56. Gu Y, Khan MA, Hamed Y, Felemban BF. A comprehensive mathematical model for SARS-CoV-2 in Caputo derivative. Fractal Fract. 2021;5(4):271. doi:10.3390/fractalfract5040271.

[picture] | Copyright © 2024 The Authors. Published by Tech Science Press. This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.