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

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∗.