# Analysis of a Tuberculosis Infection Model considering the Influence of Saturated Recovery (Treatment)

Authors: Fatima Sulayman; Farah Aini Abdullah
Journal: Complexity
 (2021)

---

## Abstract

Tuberculosis (TB) is a serious global health threat that is caused by the bacteriumMycobacterium tuberculosis, is extremely infectious, and has a high mortality rate. In this paper, we developed a model of TB infection to predict the impact of saturated recovery. The existence of equilibrium and its stability has been investigated based on the effective reproduction number RC. The model displays interesting dynamics, including backward bifurcation and Hopf bifurcation, which further results in the stable disease-free and stable endemic equilibria to be coexisting. Bifurcation analysis demonstrates that the saturation parameter is accountable for the phenomenon of backward bifurcation. We derive a condition that guarantees that the model is globally asymptotically stable using the Lyapunov function approach to global stability. The numerical simulation also reveals that the extent of saturation of TB infection is the mechanism that is fuelling TB disease in the population.


---

## Body

## 1. Introduction

Modeling and simulation are significant decision tools for controlling human diseases [1, 2]. Although every disease demonstrates its own series of biological characteristics, the models should be adjusted to every case in order to tackle real-world scenarios [3, 4].Tuberculosis (TB) kills more people every year than any other infectious disease, including HIV and malaria, making it one of the most serious global health challenges [5]. Despite the fact that TB is a curable and preventable disease, it is still one of the top ten causes of mortality worldwide [5, 6]. TB is usually caused by the bacterium Mycobacterium tuberculosis. In 2018, there were 10 million new cases of TB reported, with 8.6% of those living with HIV [5]. To be infected with TB, a person just has to inhale a few TB bacteria [5, 6]. The infection may also affect other areas of the body, like the kidneys, brain, skin, and spine, in some cases. The TB infection is not new to the community, yet it has a long history of presence since antiquity in China, Egypt, and India [7]. Tuberculosis is transferred from one infected individual to another by germ-infested air droplets. When a person with active tuberculosis germs spits, talks, coughs, or sneezes, the TB germs are propelled into the air [5]. TB can be contracted effectively from relatives and companions. An individual becomes contaminated effectively when they inhale a couple of TB germs. Approximately ten percent of TB-infected people develop active TB infection, while the other ninety percent remain latent [5]. TB infection is not transmitted by people who are latently sick. Individuals with TB may be identified by skin or blood testing. People with immunocompromised diseases are more susceptible to TB (HIV and diabetic patients). A cough that lasts longer than three weeks, fever, coughing up blood, chest pain, fatigue, weight loss, exhaustion, and night sweats are all symptoms of active tuberculosis [5, 8]. TB is a transferable disease, which can be treated with medication therapy [5, 8].Treatment is a significant and viable technique that can be used within a population to control the spread of infectious diseases. In modeling techniques, the role of treatment function addresses the likelihood of treatment against the disease at a given time for each infected person. In standard epidemic models, the treatment rate is proportional to the infective’s population size; however, in general, the rate of recovery is determined by medical resources such as medications, vaccinations, hospital beds, isolation areas, and the effectiveness of the treatments. It is understood in practice that societies have limited medical services for the treatment of infections, so offering these services places considerable pressure on facilities for public health. It is therefore more realistic to consider treatment functions of the saturated form, which appear to be large as the number of infectives increases [9–11]. Wang and Ruan [12] proposed an SIR model with a constant treatment as follows:(1)hIτ,I>0,0,I=0,that modeled a treatment capacity limit. The following piecewise linear treatment function was considered by Wang and Ruan [12]:(2)hIτI,0≤I≤I0,τI0,I>I0,where I0 is the infectious threshold at which the healthcare system gains capacity; that is, treatment rises linearly with I until capacity is achieved, and then it takes its maximal value, τI. This appears to be more realistic than the standard linear function. Eckalbar and Eckalbar [13] built the following SIR model using a quadratic treatment function:(3)TI=maxτI−gI2,0,r,g>0.Apart from that, we are aware that delaying treatment for infective persons has a negative impact on treatment efficiency. Zhang and Liu [9] employed a saturated treatment function that was both continuous and differentiable, hI=τI/1+kI, where τ>0,τ corresponds to the cure rate, and the quantity k represents the delay in treatment.It has been noted that saturation in recovery or treatment function allows multiple endemic equilibria to occur whenR0 changes [9, 10, 14–17]. For R0<1, the presence of backward bifurcation (the existence of endemic equilibrium along with disease-free equilibrium) implies that only R0<1 is insufficient to eradicate the disease from the community [9, 14, 16, 17].In this situation, an endemic equilibrium is established along with the stable disease-free equilibrium (bistability of equilibria when R0<1). Moreover, some researchers have observed the presence of limit cycles alongside the stable disease-free equilibrium when R0<1 [16, 17]. Subsequently, if R0<1, the occurrence of a stable endemic equilibrium and oscillations in the infective population indicate that maximum efforts are needed to eliminate or monitor the disease [14, 16, 17].Moreover, whenR0>1, some different scenarios, for example, Hopf bifurcation and the occurrence of several endemic equilibria, also have been observed as a result of saturated recovery or treatment [12, 14, 15, 17]. For example, Cui et al. [14] proposed and investigated an SIS model with saturated recovery rate and discovered the occurrence of oscillatory behavior in the population via Hopf bifurcation because of the treatment capacity limitations. Wang and Ruan [12] presented and analyzed an SIR model with a constant treatment rate of infective. They noticed that when R0>1, the system experiences different bifurcations including saddle-node bifurcation and Hopf bifurcation. Global stability is observed without removal rate, while oscillatory behavior for constant removal rate was discovered. The behaviour of SEIR epidemic model with saturated treatment function is investigated by [18], a backward bifurcation that results in bistability emergence, and numerical result work suggested that they should improve the medical conditions to control the epidemic [19]. To fully understand the influence of delayed treatment on transmission of the dynamics of disease, Zhang et al. [20] considered an SEIR model approach with saturated incidence and saturated treatment method. Their finding recommends that to eliminate the disease, they should raise the effectiveness and extend the capacity of treatment. In other words, our medical technology and contribution of more drugs and beds provide timely treatment to patients.In this paper, we incorporated saturated recovery which was excluded in previous four-dimensionalSEIR-like models with exogenous infection used to study the phenomenon of backward bifurcation (see [21, 22]).


## 2. Construction of the Model

We examined the transmission dynamics of TB by employing a nonlinear ODE system ofSιEιIιRι type system, which was developed in [21, 22]. Our aim is to better understand the interplay between exogenous reinfection and saturated treatment (recovery) rate in determining the outbreaks of TB. The whole population Nιt is classified into four (9) compartments, Sι, Eι, Iι, and Rι, which, respectively, denote susceptible, exposed, infected, and treated individuals. A normal mass action known as density-dependent incidence form interaction for the TB epidemic is examined and the population is considered to also be homogeneously mixed. The Sι,Eι,Iι,Rι process of TB spreading is shown in Figure 1.Figure 1
Schematic diagram showing the dynamics of TB infection with saturated recovery.


### 2.1. The Model with Saturated Recovery Rate

The nonlinear model that is considered in this study consists of a system of ordinary differential equationsODEs:(4)dSιdt=Λ−μSι−βSιIι,dEιdt=βSιIι−pβEιIι−κ+μEι,dIιdt=pβEιIι+κEι−μ+τ+δIι−cI1+ϖIι,dRιdt=cIι1+ϖIι+τI−μRι,with the initial conditions given by the following equation:(5)Sι0=Sι0>0,Eι0=Eι0>0,Iι0=Iι>0,Rι0=Rι0>0.In this epidemiological system, the susceptible compartment is increased by recruiting individuals, either by immigration or birth, into the population at the constant rateΛ. The term μ is the natural death rate. The exposed compartment becomes infectious at the rate of κE and progresses to actively infected state. Exogenous reinfection can enable individuals that are already infected to develop active TB. This happens when they contract new infection from infectious individuals at a constant rate pβEI. Infected individuals are decreased in number because of treatment, which is given at the rate τI, cI/1+ϖI. δ respectively. The term cI/1+ϖI is the saturation form of treatment (recovery) to the infected populations with the recovery rate c, and the parameter ϖ determines the magnitude of the impact of delaying the treatment of infected individuals.


## 2.1. The Model with Saturated Recovery Rate

The nonlinear model that is considered in this study consists of a system of ordinary differential equationsODEs:(4)dSιdt=Λ−μSι−βSιIι,dEιdt=βSιIι−pβEιIι−κ+μEι,dIιdt=pβEιIι+κEι−μ+τ+δIι−cI1+ϖIι,dRιdt=cIι1+ϖIι+τI−μRι,with the initial conditions given by the following equation:(5)Sι0=Sι0>0,Eι0=Eι0>0,Iι0=Iι>0,Rι0=Rι0>0.In this epidemiological system, the susceptible compartment is increased by recruiting individuals, either by immigration or birth, into the population at the constant rateΛ. The term μ is the natural death rate. The exposed compartment becomes infectious at the rate of κE and progresses to actively infected state. Exogenous reinfection can enable individuals that are already infected to develop active TB. This happens when they contract new infection from infectious individuals at a constant rate pβEI. Infected individuals are decreased in number because of treatment, which is given at the rate τI, cI/1+ϖI. δ respectively. The term cI/1+ϖI is the saturation form of treatment (recovery) to the infected populations with the recovery rate c, and the parameter ϖ determines the magnitude of the impact of delaying the treatment of infected individuals.


## 3. Basic Properties of the TB Model with Saturated Recovery

Following the technique in [23, 24] (and furthermore the proof in [24]), it is not hard to demonstrate that when all model parameters remain nonnegative, the state variables Sι,Eι,Iι,andRι are all positive for all time t given they are initially positive.Theorem 1.
The region(6)Φ=Sι,Eι,Iι,Rι∈R+4,Nι≤Λμ,is positively invariant and attracts all nonnegative solutions of model system (4).Proof.
As mentioned in the study of [23], considering the nonlinear system of model (4), we take the first equation and let λ=βIι.(7)dSιdt=Λ−λ+μ,that is,(8)dSιdt≥−λ+μS,(9)⟹dSιSι≥λ+μdt.
Integrating (9) by separation of variables gives(10)∫dSιSι≥∫λ+μdt,and therefore(11)Sιt≥S0e−λ+μt≥0.
Similarly, it can also been shown thatEιt>0, Iιt>0, Rιt>0, for all t>0.


### 3.1. Invariant Region

Theorem 2.
The solution of TB model system (4) is enclosed in the region Φ subset ∈R+6, given by(12)Φ=Sι,Eι,Iι,Rι∈R+4,Nι≤Λμ,for the initial conditions (5) in Φ.


## 3.1. Invariant Region

Theorem 2.
The solution of TB model system (4) is enclosed in the region Φ subset ∈R+6, given by(12)Φ=Sι,Eι,Iι,Rι∈R+4,Nι≤Λμ,for the initial conditions (5) in Φ.


## 4. Analysis of Disease-Free EquilibriumDFE, P0, and Effective Basic Reproduction Number RC

The disease-free equilibriumDFE state, P0, of system model (4) is obtained by setting the right-hand side of (4) to zero which is given by(13)P0=Sι,Eι,Iι,RιΛμ,0,0,0.The basic reproduction number of system model (4) was calculated using the next generation operator approach (see [25, 26]), where it is obtained to be(14)RC=βκΛμμ+κμ+δ+τ+c,where RC is the effective reproduction number of system model (4). The following outcomes would guaranty the stability of P0, when the threshold quantity RC<1 ascertains that the system will be free from TB infection under a threshold condition. A significant result is the following theorem.Theorem 3.
The disease-free equilibrium state,P0, of model system (4) is locally asymptotically stable (LAS) when RC<1 and unstable if RC>1.Proof.
Consider model system (4), and the Jacobian of system (4) at P0 is given by(15)ΣP0=−μ0−βΛμ00−κ+μβΛμ00κ−μ+δ+τ+cRC−1000τ+c−μ.
Then, the root characteristic equation ofΣP0 is −μ, −κ+μ, −μ, and −μ+δ+τ+cRC−1. If RC<1, all the characteristic roots are negative and thus P0 is locally asymptotically stable and unstable when RC>1.


### 4.1. Global Stability ofP0 State

In order to ensure that TB is independent of the initial size of the subpopulations of model system (4), it is necessary to show that P0 is globally asymptotically stable (GAS) [27–30]. Several approaches, such as the Lyapunov theorem, [27, 29, 31, 32] and global stability theorem have been used to prove the global stability of the disease free equilibrium, P0. Here, we shall apply the method that was introduced in [32].Theorem 4.
The disease-free equilibriumP0 of model system (4) is global asymptotically stable in Ω if RC<1 and conditions H1 and H2 are satisfied as given in [32].Proof.
The two conditionsH1 and H2 in Castillo-Chavez global stability theorem must be satisfied for RC<1 to ascertain the global stability of the P0. Write system model (4) in the following form:(16)dXdt=FX,Z,dZdt=GX,Z;GX,Z=0.
Here, the componentsX=Sι∗,Rι∗ and Z=Eι∗,Iι∗ where X∈R2 denotes uninfected population and Z∈R2 denotes the infected population. The disease-free equilibrium is defined by(17)P0=X∗,0,where(18)X∗=Λμ,0.
For conditionH1 (the global asymptotical stability of X∗) to be fulfilled, we have(19)dXdt=FX,0Λ−μSι∗μRι∗,written in the form of linear differential equations as follows:(20)dSι∗dt=Λ−μSι∗,dRι∗dt−μRι.
Solving linear differential equation (20) yields(21)Sι∗t=Λμ+Sι∗0e−μt,Rι∗t=Rι∗0e−μt.
Now, clearlyRι∗t⟶0 and Sι∗t⟶Λ/μ as t⟶∞, regardless of the value of Sι∗0 and Rι∗0.
Hence,X∗=Λ/μ,0 is globally asymptotically stable.
Next, applying second condition of the theorem proposed in [32], that is, G⌢X,Z=AZ−GX,Z, G⌢X,Z≥0, we have(22)GX,Z=G1X,ZG2X,Z=βSιIι−μ+κ+pβIιEικEι−μ+δ+τ+cIι1+ϖIι,A=∂G∂ZX∗,0−μ+κ0κ−μ+δ+τ+cϖ.
This is certainly anM- matrix (the off-diagonal elements of A are nonnegative) [32].
Since(23)GX,Z=−μ+κEι∗0κE∗−μ+δ+τ+cϖIι∗,then,(24)G^X,Z=AZ−GX,Z=−μ+κ0κ−μ+δ+τ+cϖEι∗Iι∗−−μ+κEι∗0κEι∗−μ+δ+τ+cϖIι∗=−μ+κEι∗0κEι∗−μ+δ+τ+cϖIι∗−−μ+κEι∗0κEι∗−μ+δ+τ+cϖIι∗=00,so that(25)G^X,Z=AZ−GX,Z=00.
SinceG^X1,X2≥0, we have showed that condition H2 of the theorem proposed in [32, 33] is satisfied by our model. This shows that regardless of the initial population of infected, TB still can be controlled.


### 4.2. Endemic Equilibria and Bifurcation Analysis

To find endemic equilibriumSι∗,Eι∗,Iι∗,Rι∗ as before, we set the right-hand side of system model (4) to zero. Solving for Sι∗, Eι∗, and Rι∗ from first, third, and last equation of (4), respectively, we have(26)Sι∗=Λμ+βIι∗,(27)Eι∗=M+c/1+ϖIι∗Iι∗p∗βIι∗,where M=μ+τ+δ and(28)Rι∗=τ+c/1+ϖIι∗Iι∗μ.By substituting and simplifying equations (26)–(28), we get the following equation in Iι∗:(29)FIι∗=A˜Iι3+B˜Iι2+C˜Iι+D˜=0,where(30)A˜=ϖpβ2M,B˜=p−ϖΛ+M+cβ+ϖMp+1μ+κβ,C˜=−pϖΛβ2+p+1M+cμ+κ−ϖΛ+M+cβ+μϖMκ+μ,D˜=κ+μμ+τ+δ+c1−RC.Clearly,A˜>0. When RC>1, then D˜<0 and thus F0<0. Furthermore, as as Iι⟶∞, F∞>0. Accordingly, from the continuity of F0, at least one positive Iι exists, such that FIι=0, and thus there will be at least one endemic equilibrium of the model system (4). Notice that if RC>1 and B˜,C˜ are positive, then model system (4) has a unique endemic equilibrium defined as P∗=Sι∗,Eι∗,Iι∗,Rι∗ from Descartes’ rule of sign. The number of possible positive real roots of the cubic polynomial (29) depends on the signs of B˜,C˜, and D˜. We, therefore, propose the following results.Theorem 5.
The TB model system (4) has(1)
A unique endemic equilibrium whenRC>1.(2)
One or more than one endemic equilibrium whenB˜<0,C˜>0,andRC<1.(3)
No endemic equilibrium whenRC<1.Now we will look at the polynomial in (29) to see how it influences the existence of endemic equilibrium.Ifϖ=0, the polynomial (29) reduces to a linear equation in Iι∗, with a unique solution(31)Iι∗=−D˜C˜,which is positive if and only if RC>1. As a result, we can now say that if ϖ=0, then a unique endemic equilibrium exist if RC>1 and there cannot be an endemic equilibrium when RC<1. To put it in another way, if the treatment of certain infected outpatients is not delayed, when the effective reproduction number is smaller than unity, an endemic steady state solution is not possible, ignoring the possibility of a backward bifurcation in this scenario.


### 4.3. Local Stability of Endemic EquilibriumP∗

The Jacobian matrix of (4) at P∗ is given as(32)ΣP∗=Υ0χ0βIι−κ+μ+pβIι−pβEι+βSι00pβIι+κΨ000Ξ−μ.Here(33)Υ=−μ+βIι,χ=−βSι,Ψ=pβEι−M−c1+ϖIι+cϖIι1+ϖIι2,Ξ=τ+c1+ϖIι−cϖIι1+ϖIι2.We obtain the characteristic equation after some row and column operations ofΣP∗ given by(34)λ4+A1λ3+A2λ2+A3λ+A4=0,where(35)A1=pβIι−Ψ+κ+2μ−Υ,A2=ppEι−SιIιβ2+κEι+μ−Ψ−ΥpIι−Sικβ+μ−Ψ−Υκ+μ2+−2Ψ−2Yμ+ΨΥ,A3=ppEι−SιμIι+pEι−SιΥ+χIιβ2+pIιΥ+κEι−ΨIιp−κSι+−κEι+Ψιp+κSιΥ−κχIιβ+−Υ+Ψμ2+2Ψ−κΥ−κΨμ+κΨΥ,A4=μpIιpEι+Sιβ2+PΨIι−κpEι−Sιβ+Ψκ+μΥ−βχIιpβIι+κ.Thus, all the roots of the characteristic equationΣP∗ are with negative real parts if A1>0, A3>0, A4>0, and A1A2A3>A32+A12A4. Therefore, from Routh–Hurwitz criterion [34], the endemic equilibrium state P∗ is locally asymptotically stable if these conditions are true. The following result can be obtained.Theorem 6.
IfRC>1 and B˜ and D˜ are positive, then the unique endemic equilibrium P∗ of system model (4) is locally asymptotically stable provided A3>0,A4>0,andA1A2A3>A32+A12A4.


### 4.4. Analysis of Backward Bifurcation

The existence and stability of endemic equilibrium are determined by investigating the possibility of backward or forward bifurcation due to the existence of endemic equilibrium. To investigate the possibility of backward or forward bifurcation of model system (4), we use the center manifold theory [35, 36]. Let the bifurcation parameter be β=β∗.Firstly, we obtained the bifurcation parameter atR0=1, and thus(36)βκΛμκ+μμ+δ+τ+c=1,and therefore(37)β=β∗=μκ+μμ+δ+γ+cΛκ.To investigate the use of center manifold theory in [35], it is convenient to make simplification and transform the variables on model system (4). This is done by rewriting our system model (4).Let(38)x1=Sι,x2=Eι,x3=Iι,x4=Rι.Moreover, by using the vector notationV=x1,x2,x3,x4T, model system (4) can be restated in the form of dV/dt=f1,f2,f3,x4T as follows:(39)dx1dt=Λ−βx1x3−μx1=f1,dx2dt=βx1x3−pβx2x3−κ+μx2=f2,dx3dt=κx2+pβx2x3−μ+δ+τx3−cx31+ϖx3=f3,dx4dt=τx3+cx31+ϖx3−μx4=f4.We next show that the Jacobian matrix of (4) at the point P0,β∗ has a simple zero eigenvalue, that is,(40)ΣP0,β∗=−μ0−β∗Λμ00−κ+μβ∗Λμ00κ−μ+δ+τ+c000τ+c−μ.Withβ=β∗, system model (39) has a simple zero real part eigenvalue, and all other eigenvalues are negative (i.e., has a hyperbolic equilibrium point). The center manifold theory [35] can therefore be used to examine the dynamics of the system (39) close to β=β∗. It is practicable to derive the right eigenvectors of ΣP0,β∗ represented by ω=ω1,ω2,ω3,ω4T, where(41)ω1=βΛμμ+κ,ω2=μ+τ+δ+cκ,ω3=1,ω4=τ+cμ.Conversely, the left eigenvectors represented asυ=υ1,υ2,υ3,υ4T are easy to obtain, where(42)υ1=0,υ2=κ2μκ+τδ+c,κ+μ2μκ+τδ+c,υ4=0.


### 4.5. Computation ofa and b

Hence, associated bifurcation coefficients are denoted bya and b, respectively; as explained in [35], when bifurcation coefficients a and b are both negative, then the system undergoes a backward bifurcation; otherwise, forward bifurcation will occur.It is convenient to find botha and b defined by [35](43)a=∑k,i,j=14vkwiwj∂2fk∂xi∂xjP0,β∗,(44)b=∑k,i=14vkwi∂2fk∂xi∂βP0,β∗.Taking into account model system (39) and examining a and b only, which are not equal to zero derivatives, for the terms(45)∂2fk∂xi∂xjP0,β∗,∂2fk∂xi∂β∗P0,β∗,we get(46)a=2υ1ω1wω∂2f1∂x1∂x4P0,β∗+2υ2ω2ω4∂2f2∂x2∂x4P0,β∗+2υ3ω1ω4∂2f3∂x1∂x4P0,β∗+2υ3ω2ω4∂2f3∂x3∂x4P0,β∗+2υ3ω3ω4∂2f3∂x3∂x4P0,β∗+2υ4ω3ω4∂2f4∂x3∂x4P0,β∗,b=2υ1υ4∂2f1∂x4∂βP0,β∗+υ3ω2ω4∂2f3∂x4∂βP0,β∗.Now considering (43) and (44), we obtain(47)a=22μ+κ2μ+τ+δ+c−μ+κμ+τ+δ+cκΛμ−μ+τ+δ+c+cϖ,b=κ2μ+κ+τ+δ+cΛμ.As indicated by the result shown in [35], our system experiences backward bifurcation at β=β∗, only when both a and b are positive at P0,β∗. Obviously, b is consistently positive. Hence, the positivity of a offers the threshold circumstance for the phenomenon of backward bifurcation.


### 4.6. Analysis of Hopf Bifurcation

Hopf bifurcation of system (4) around endemic equilibrium point will be considered in this section. We employ the criterion that was derived in [37] to show that Hopf bifurcation exists. Let P∗∗=Sι∗,Eι∗,Iι∗,Rι∗ be an endemic equilibrium of system (4) when RC>1 and β be a bifurcation parameter for which P∗∗ is locally stable when β<βo and unstable when β>βo. At β=βo, P∗∗ loses its stability and periodic solutions emerge. The characteristic equation at endemic equilibrium equivalent to model system (4) at P∗∗ is given by(48)X4+A1X3+A2X2+A3X+A4=0.It was also observed that the component of the endemic equilibriumP∗∗ smoothly relies on β, and the corresponding characteristic equation is given as(49)X4+A1βX3+A2βX2+A3βX+A4β=0.To this place, the coefficient is a smooth function ofβ. Following the approach in [38], a Hopf bifurcation occurs when the following two conditions are satisfied:The Jacobian matrixJP∗∗,βo has a pair of purely imaginary eigenvalues and the remaining eigenvalues have negative real parts.dReA/dβ|β=βo≠0.According to [37], the coefficients of the characteristic equation have to meet the prerequisite for a pair of purely imaginary eigenvalues, which areA4>0, A4>0 and A3A2−A4A1>0.Δ2=A3A2A1−A32−A4A12=0.To complete the discussion, we need to drive the transversality condition (ii). For this reason, we let±iα be a pair of purely imaginary eigenvalues corresponding to βo. Here, differentiating characteristic equation (48) with respect to β, we get(50)4X3+3A1X2+2A2X+A3dXdβ+X3dA1dβ+X2dA2dβ+XdA3dβ+dA4dβ=0.Further, we get(51)dXdβ−1=−4X3+3A1X2+2A2X+A3X3dA1/dβ+X2dA2/dβ+XdA3/dβ+dA4/dβ.The following relation will be used.(52)signdReXdβX=iα,Δ2=0=signRedXdβ−1X=iα,Δ2=0=signRe−A3−3A1α2+2A2α−4α3idA4/dβ−dA2/dβα2+dA3/dβα−dA1/dβα3i=signΥ.Here(53)Υ=−A3−3A1α2dA4/dβ−dA2/dβα2+2A2α−4α3dA3/dβα−dA1/dβα3dA4/dβ−dA2/dβα22+dA3/dβα−dA1/dβα32.IfA3−3A1α2dA4/dβ−dA2/dβα2+2A2α−4α3dA3/dβα−dA1/dβα3<0, then signdReX/dββo,α>0. We then summarize this result in the following equation:(54)A3−3A1α2dA4dβ−dA2dβα2+2A2α−4α3dA3dβα−dA1dβα3<0.Therefore, the transversality condition holds, and the system undergoes the phenomena of Hopf bifurcation atβo=α. Thus, transmission rate βo crosses its critical value, βo=α, and individuals start oscillating around endemic equilibrium point P∗∗. This ends the proof of the theorem.


### 4.7. Global Stability of Endemic Equilibrium WhenRC>1

To prove the global stability of endemic equilibrium stateP∗∗ of system model (4), we consider the case when exogenous reinfection parameter p=0. We employed the Lyapunov function method and LaSalle’s invariance principles [39, 40]. System model (4) with p=0 therefore reduces to(55)dSιdt=Λ−μSι−βSιIι,dEιdt=βSιIι−l1Eι,dIιdt=κEι−l2Iι−fIι,dRιdt=fIι+τI−μRι,with(56)l1=μ+κ,l2=μ+τ+δ,fIι=cIι1+ϖIι.Furthermore, define(57)Γo=Sι,Eι,Iι,Rι≥Γ:Eι=Iι=Rι=0.The following result is obtained.Theorem 7.
IfRC>1, then the endemic equilibrium state P∗∗ given by (55) is GAS in Γ\Γo.Proof.
Consider the model without exogenous reinfection given by (55), and RC>1, so that the endemic equilibrium P∗∗ of the model (55) exists. In addition, consider the following Lyapunov function method of this type which has been utilized in the study of epidemiology and ecology, such as in [29, 41–43], for model system (55) comprising the first three equations of (55), defined by the following. Define the Lyapunov function(58)USι,Eι,Iι=Sι−Sι∗∗−Sι∗∗lnSSι∗∗+Eι−Eι∗∗−Eι∗∗lnEEι∗∗+l1κIι−Iι∗∗−Iι∗∗lnIIι∗∗,and the derivatives of U are given by(59)U´=S´ι−Sι∗∗SιS´+E´−Eι∗∗EιE´ι+l1κI´ι−Iι∗∗IιI´ι.
Filling in the derivatives and using (55) in U′, we have(60)=Λ−βSιIι−μSι−Sι∗∗SιΛ−βSιIι−μSι+βSιIι−l1Eι−Eι∗∗EιβSιIι−l1Eι+l1κκEι−l2Iι−Iι∗∗IικEι−l2Iι,(61)U´=Λ1−Sι∗∗Sι−μSι1−Sι∗∗Sι+βSι∗∗−l1l2κ−Eι∗∗βSιIιEι+l1Eι∗∗−l1Iι∗∗EιIι+l1l2κ.
At the endemic steady state, it can be observed that(62)βSι∗∗=μ+κμ+τ+δ+fκ=l1l2κ.
Using expression (62) in equation (61) yields(63)U´=Λ1−Sι∗∗Sι−μSι1−Sι∗∗Sι−Eι∗∗βSιIιEι+Eι∗∗l1−l1Iι∗∗EιIι+βSι∗∗Iι∗∗,U´=βSι∗∗Iι∗∗+μSι∗∗1−Sι∗∗Sι−μSι1−Sι∗∗Sι−Eι∗∗βSιIιEι+Eι∗∗l1−l1Iι∗∗EιIι+βSι∗∗Iι∗∗,U´=μSι∗∗2−Sι∗∗Sι−SιSι∗∗+βSι∗∗Iι∗∗−βSι∗∗2Sι−βSιIιEι∗∗Eι+l1Eι∗∗−l1Iι∗∗EιIι+βSι∗∗Iι∗∗.
We also observed from (55) at steady state that(64)l1Eι∗=βSι∗∗Iι∗∗.
Using expression (64) in equation (63) gives(65)U´=μSι∗∗2−Sι∗∗Sι−SιSι∗∗+βSι∗∗Iι∗∗3−Sι∗∗Sι−SιIιEι∗∗Sι∗∗Iι∗∗Eι−EIι∗∗IιEι∗∗.
We can simplify the first term in the last equation of (65) as follows.(66)2−Sι∗∗Sι−SιSι∗∗=2SιSι∗∗−Sι∗∗2−Sι2SιSι∗∗=−Sι−Sι∗∗2Sι−Sι∗∗≤0.
Finally, it follows that since the arithmetic mean exceeds the geometric mean,(67)3−Sι∗∗Sι−SιIιEι∗∗Sι∗∗Iι∗∗Eι−EIι∗∗IιEι∗∗≤0.
Accordingly, it follows from this thatU´≤0 when RC> with U´=0 if and only if Sι=Sι∗∗, Eι=Eι∗∗, and Iι=Iι∗∗. Thus, U is a Lyapunov function for system model (55) in Γ∖Γ0. Thus, the largest compact invariant set where U´=0 is a singleton Sι,Eι,Iι=Sι∗∗,Eι∗∗,Iι∗∗, and therefore it follows that by [39], Sιt⟶Sι∗∗, Eιt⟶Eι∗∗, and Iιt⟶Iι∗∗ as t⟶∞.


## 4.1. Global Stability ofP0 State

In order to ensure that TB is independent of the initial size of the subpopulations of model system (4), it is necessary to show that P0 is globally asymptotically stable (GAS) [27–30]. Several approaches, such as the Lyapunov theorem, [27, 29, 31, 32] and global stability theorem have been used to prove the global stability of the disease free equilibrium, P0. Here, we shall apply the method that was introduced in [32].Theorem 4.
The disease-free equilibriumP0 of model system (4) is global asymptotically stable in Ω if RC<1 and conditions H1 and H2 are satisfied as given in [32].Proof.
The two conditionsH1 and H2 in Castillo-Chavez global stability theorem must be satisfied for RC<1 to ascertain the global stability of the P0. Write system model (4) in the following form:(16)dXdt=FX,Z,dZdt=GX,Z;GX,Z=0.
Here, the componentsX=Sι∗,Rι∗ and Z=Eι∗,Iι∗ where X∈R2 denotes uninfected population and Z∈R2 denotes the infected population. The disease-free equilibrium is defined by(17)P0=X∗,0,where(18)X∗=Λμ,0.
For conditionH1 (the global asymptotical stability of X∗) to be fulfilled, we have(19)dXdt=FX,0Λ−μSι∗μRι∗,written in the form of linear differential equations as follows:(20)dSι∗dt=Λ−μSι∗,dRι∗dt−μRι.
Solving linear differential equation (20) yields(21)Sι∗t=Λμ+Sι∗0e−μt,Rι∗t=Rι∗0e−μt.
Now, clearlyRι∗t⟶0 and Sι∗t⟶Λ/μ as t⟶∞, regardless of the value of Sι∗0 and Rι∗0.
Hence,X∗=Λ/μ,0 is globally asymptotically stable.
Next, applying second condition of the theorem proposed in [32], that is, G⌢X,Z=AZ−GX,Z, G⌢X,Z≥0, we have(22)GX,Z=G1X,ZG2X,Z=βSιIι−μ+κ+pβIιEικEι−μ+δ+τ+cIι1+ϖIι,A=∂G∂ZX∗,0−μ+κ0κ−μ+δ+τ+cϖ.
This is certainly anM- matrix (the off-diagonal elements of A are nonnegative) [32].
Since(23)GX,Z=−μ+κEι∗0κE∗−μ+δ+τ+cϖIι∗,then,(24)G^X,Z=AZ−GX,Z=−μ+κ0κ−μ+δ+τ+cϖEι∗Iι∗−−μ+κEι∗0κEι∗−μ+δ+τ+cϖIι∗=−μ+κEι∗0κEι∗−μ+δ+τ+cϖIι∗−−μ+κEι∗0κEι∗−μ+δ+τ+cϖIι∗=00,so that(25)G^X,Z=AZ−GX,Z=00.
SinceG^X1,X2≥0, we have showed that condition H2 of the theorem proposed in [32, 33] is satisfied by our model. This shows that regardless of the initial population of infected, TB still can be controlled.


## 4.2. Endemic Equilibria and Bifurcation Analysis

To find endemic equilibriumSι∗,Eι∗,Iι∗,Rι∗ as before, we set the right-hand side of system model (4) to zero. Solving for Sι∗, Eι∗, and Rι∗ from first, third, and last equation of (4), respectively, we have(26)Sι∗=Λμ+βIι∗,(27)Eι∗=M+c/1+ϖIι∗Iι∗p∗βIι∗,where M=μ+τ+δ and(28)Rι∗=τ+c/1+ϖIι∗Iι∗μ.By substituting and simplifying equations (26)–(28), we get the following equation in Iι∗:(29)FIι∗=A˜Iι3+B˜Iι2+C˜Iι+D˜=0,where(30)A˜=ϖpβ2M,B˜=p−ϖΛ+M+cβ+ϖMp+1μ+κβ,C˜=−pϖΛβ2+p+1M+cμ+κ−ϖΛ+M+cβ+μϖMκ+μ,D˜=κ+μμ+τ+δ+c1−RC.Clearly,A˜>0. When RC>1, then D˜<0 and thus F0<0. Furthermore, as as Iι⟶∞, F∞>0. Accordingly, from the continuity of F0, at least one positive Iι exists, such that FIι=0, and thus there will be at least one endemic equilibrium of the model system (4). Notice that if RC>1 and B˜,C˜ are positive, then model system (4) has a unique endemic equilibrium defined as P∗=Sι∗,Eι∗,Iι∗,Rι∗ from Descartes’ rule of sign. The number of possible positive real roots of the cubic polynomial (29) depends on the signs of B˜,C˜, and D˜. We, therefore, propose the following results.Theorem 5.
The TB model system (4) has(1)
A unique endemic equilibrium whenRC>1.(2)
One or more than one endemic equilibrium whenB˜<0,C˜>0,andRC<1.(3)
No endemic equilibrium whenRC<1.Now we will look at the polynomial in (29) to see how it influences the existence of endemic equilibrium.Ifϖ=0, the polynomial (29) reduces to a linear equation in Iι∗, with a unique solution(31)Iι∗=−D˜C˜,which is positive if and only if RC>1. As a result, we can now say that if ϖ=0, then a unique endemic equilibrium exist if RC>1 and there cannot be an endemic equilibrium when RC<1. To put it in another way, if the treatment of certain infected outpatients is not delayed, when the effective reproduction number is smaller than unity, an endemic steady state solution is not possible, ignoring the possibility of a backward bifurcation in this scenario.


## 4.3. Local Stability of Endemic EquilibriumP∗

The Jacobian matrix of (4) at P∗ is given as(32)ΣP∗=Υ0χ0βIι−κ+μ+pβIι−pβEι+βSι00pβIι+κΨ000Ξ−μ.Here(33)Υ=−μ+βIι,χ=−βSι,Ψ=pβEι−M−c1+ϖIι+cϖIι1+ϖIι2,Ξ=τ+c1+ϖIι−cϖIι1+ϖIι2.We obtain the characteristic equation after some row and column operations ofΣP∗ given by(34)λ4+A1λ3+A2λ2+A3λ+A4=0,where(35)A1=pβIι−Ψ+κ+2μ−Υ,A2=ppEι−SιIιβ2+κEι+μ−Ψ−ΥpIι−Sικβ+μ−Ψ−Υκ+μ2+−2Ψ−2Yμ+ΨΥ,A3=ppEι−SιμIι+pEι−SιΥ+χIιβ2+pIιΥ+κEι−ΨIιp−κSι+−κEι+Ψιp+κSιΥ−κχIιβ+−Υ+Ψμ2+2Ψ−κΥ−κΨμ+κΨΥ,A4=μpIιpEι+Sιβ2+PΨIι−κpEι−Sιβ+Ψκ+μΥ−βχIιpβIι+κ.Thus, all the roots of the characteristic equationΣP∗ are with negative real parts if A1>0, A3>0, A4>0, and A1A2A3>A32+A12A4. Therefore, from Routh–Hurwitz criterion [34], the endemic equilibrium state P∗ is locally asymptotically stable if these conditions are true. The following result can be obtained.Theorem 6.
IfRC>1 and B˜ and D˜ are positive, then the unique endemic equilibrium P∗ of system model (4) is locally asymptotically stable provided A3>0,A4>0,andA1A2A3>A32+A12A4.


## 4.4. Analysis of Backward Bifurcation

The existence and stability of endemic equilibrium are determined by investigating the possibility of backward or forward bifurcation due to the existence of endemic equilibrium. To investigate the possibility of backward or forward bifurcation of model system (4), we use the center manifold theory [35, 36]. Let the bifurcation parameter be β=β∗.Firstly, we obtained the bifurcation parameter atR0=1, and thus(36)βκΛμκ+μμ+δ+τ+c=1,and therefore(37)β=β∗=μκ+μμ+δ+γ+cΛκ.To investigate the use of center manifold theory in [35], it is convenient to make simplification and transform the variables on model system (4). This is done by rewriting our system model (4).Let(38)x1=Sι,x2=Eι,x3=Iι,x4=Rι.Moreover, by using the vector notationV=x1,x2,x3,x4T, model system (4) can be restated in the form of dV/dt=f1,f2,f3,x4T as follows:(39)dx1dt=Λ−βx1x3−μx1=f1,dx2dt=βx1x3−pβx2x3−κ+μx2=f2,dx3dt=κx2+pβx2x3−μ+δ+τx3−cx31+ϖx3=f3,dx4dt=τx3+cx31+ϖx3−μx4=f4.We next show that the Jacobian matrix of (4) at the point P0,β∗ has a simple zero eigenvalue, that is,(40)ΣP0,β∗=−μ0−β∗Λμ00−κ+μβ∗Λμ00κ−μ+δ+τ+c000τ+c−μ.Withβ=β∗, system model (39) has a simple zero real part eigenvalue, and all other eigenvalues are negative (i.e., has a hyperbolic equilibrium point). The center manifold theory [35] can therefore be used to examine the dynamics of the system (39) close to β=β∗. It is practicable to derive the right eigenvectors of ΣP0,β∗ represented by ω=ω1,ω2,ω3,ω4T, where(41)ω1=βΛμμ+κ,ω2=μ+τ+δ+cκ,ω3=1,ω4=τ+cμ.Conversely, the left eigenvectors represented asυ=υ1,υ2,υ3,υ4T are easy to obtain, where(42)υ1=0,υ2=κ2μκ+τδ+c,κ+μ2μκ+τδ+c,υ4=0.


## 4.5. Computation ofa and b

Hence, associated bifurcation coefficients are denoted bya and b, respectively; as explained in [35], when bifurcation coefficients a and b are both negative, then the system undergoes a backward bifurcation; otherwise, forward bifurcation will occur.It is convenient to find botha and b defined by [35](43)a=∑k,i,j=14vkwiwj∂2fk∂xi∂xjP0,β∗,(44)b=∑k,i=14vkwi∂2fk∂xi∂βP0,β∗.Taking into account model system (39) and examining a and b only, which are not equal to zero derivatives, for the terms(45)∂2fk∂xi∂xjP0,β∗,∂2fk∂xi∂β∗P0,β∗,we get(46)a=2υ1ω1wω∂2f1∂x1∂x4P0,β∗+2υ2ω2ω4∂2f2∂x2∂x4P0,β∗+2υ3ω1ω4∂2f3∂x1∂x4P0,β∗+2υ3ω2ω4∂2f3∂x3∂x4P0,β∗+2υ3ω3ω4∂2f3∂x3∂x4P0,β∗+2υ4ω3ω4∂2f4∂x3∂x4P0,β∗,b=2υ1υ4∂2f1∂x4∂βP0,β∗+υ3ω2ω4∂2f3∂x4∂βP0,β∗.Now considering (43) and (44), we obtain(47)a=22μ+κ2μ+τ+δ+c−μ+κμ+τ+δ+cκΛμ−μ+τ+δ+c+cϖ,b=κ2μ+κ+τ+δ+cΛμ.As indicated by the result shown in [35], our system experiences backward bifurcation at β=β∗, only when both a and b are positive at P0,β∗. Obviously, b is consistently positive. Hence, the positivity of a offers the threshold circumstance for the phenomenon of backward bifurcation.


## 4.6. Analysis of Hopf Bifurcation

Hopf bifurcation of system (4) around endemic equilibrium point will be considered in this section. We employ the criterion that was derived in [37] to show that Hopf bifurcation exists. Let P∗∗=Sι∗,Eι∗,Iι∗,Rι∗ be an endemic equilibrium of system (4) when RC>1 and β be a bifurcation parameter for which P∗∗ is locally stable when β<βo and unstable when β>βo. At β=βo, P∗∗ loses its stability and periodic solutions emerge. The characteristic equation at endemic equilibrium equivalent to model system (4) at P∗∗ is given by(48)X4+A1X3+A2X2+A3X+A4=0.It was also observed that the component of the endemic equilibriumP∗∗ smoothly relies on β, and the corresponding characteristic equation is given as(49)X4+A1βX3+A2βX2+A3βX+A4β=0.To this place, the coefficient is a smooth function ofβ. Following the approach in [38], a Hopf bifurcation occurs when the following two conditions are satisfied:The Jacobian matrixJP∗∗,βo has a pair of purely imaginary eigenvalues and the remaining eigenvalues have negative real parts.dReA/dβ|β=βo≠0.According to [37], the coefficients of the characteristic equation have to meet the prerequisite for a pair of purely imaginary eigenvalues, which areA4>0, A4>0 and A3A2−A4A1>0.Δ2=A3A2A1−A32−A4A12=0.To complete the discussion, we need to drive the transversality condition (ii). For this reason, we let±iα be a pair of purely imaginary eigenvalues corresponding to βo. Here, differentiating characteristic equation (48) with respect to β, we get(50)4X3+3A1X2+2A2X+A3dXdβ+X3dA1dβ+X2dA2dβ+XdA3dβ+dA4dβ=0.Further, we get(51)dXdβ−1=−4X3+3A1X2+2A2X+A3X3dA1/dβ+X2dA2/dβ+XdA3/dβ+dA4/dβ.The following relation will be used.(52)signdReXdβX=iα,Δ2=0=signRedXdβ−1X=iα,Δ2=0=signRe−A3−3A1α2+2A2α−4α3idA4/dβ−dA2/dβα2+dA3/dβα−dA1/dβα3i=signΥ.Here(53)Υ=−A3−3A1α2dA4/dβ−dA2/dβα2+2A2α−4α3dA3/dβα−dA1/dβα3dA4/dβ−dA2/dβα22+dA3/dβα−dA1/dβα32.IfA3−3A1α2dA4/dβ−dA2/dβα2+2A2α−4α3dA3/dβα−dA1/dβα3<0, then signdReX/dββo,α>0. We then summarize this result in the following equation:(54)A3−3A1α2dA4dβ−dA2dβα2+2A2α−4α3dA3dβα−dA1dβα3<0.Therefore, the transversality condition holds, and the system undergoes the phenomena of Hopf bifurcation atβo=α. Thus, transmission rate βo crosses its critical value, βo=α, and individuals start oscillating around endemic equilibrium point P∗∗. This ends the proof of the theorem.


## 4.7. Global Stability of Endemic Equilibrium WhenRC>1

To prove the global stability of endemic equilibrium stateP∗∗ of system model (4), we consider the case when exogenous reinfection parameter p=0. We employed the Lyapunov function method and LaSalle’s invariance principles [39, 40]. System model (4) with p=0 therefore reduces to(55)dSιdt=Λ−μSι−βSιIι,dEιdt=βSιIι−l1Eι,dIιdt=κEι−l2Iι−fIι,dRιdt=fIι+τI−μRι,with(56)l1=μ+κ,l2=μ+τ+δ,fIι=cIι1+ϖIι.Furthermore, define(57)Γo=Sι,Eι,Iι,Rι≥Γ:Eι=Iι=Rι=0.The following result is obtained.Theorem 7.
IfRC>1, then the endemic equilibrium state P∗∗ given by (55) is GAS in Γ\Γo.Proof.
Consider the model without exogenous reinfection given by (55), and RC>1, so that the endemic equilibrium P∗∗ of the model (55) exists. In addition, consider the following Lyapunov function method of this type which has been utilized in the study of epidemiology and ecology, such as in [29, 41–43], for model system (55) comprising the first three equations of (55), defined by the following. Define the Lyapunov function(58)USι,Eι,Iι=Sι−Sι∗∗−Sι∗∗lnSSι∗∗+Eι−Eι∗∗−Eι∗∗lnEEι∗∗+l1κIι−Iι∗∗−Iι∗∗lnIIι∗∗,and the derivatives of U are given by(59)U´=S´ι−Sι∗∗SιS´+E´−Eι∗∗EιE´ι+l1κI´ι−Iι∗∗IιI´ι.
Filling in the derivatives and using (55) in U′, we have(60)=Λ−βSιIι−μSι−Sι∗∗SιΛ−βSιIι−μSι+βSιIι−l1Eι−Eι∗∗EιβSιIι−l1Eι+l1κκEι−l2Iι−Iι∗∗IικEι−l2Iι,(61)U´=Λ1−Sι∗∗Sι−μSι1−Sι∗∗Sι+βSι∗∗−l1l2κ−Eι∗∗βSιIιEι+l1Eι∗∗−l1Iι∗∗EιIι+l1l2κ.
At the endemic steady state, it can be observed that(62)βSι∗∗=μ+κμ+τ+δ+fκ=l1l2κ.
Using expression (62) in equation (61) yields(63)U´=Λ1−Sι∗∗Sι−μSι1−Sι∗∗Sι−Eι∗∗βSιIιEι+Eι∗∗l1−l1Iι∗∗EιIι+βSι∗∗Iι∗∗,U´=βSι∗∗Iι∗∗+μSι∗∗1−Sι∗∗Sι−μSι1−Sι∗∗Sι−Eι∗∗βSιIιEι+Eι∗∗l1−l1Iι∗∗EιIι+βSι∗∗Iι∗∗,U´=μSι∗∗2−Sι∗∗Sι−SιSι∗∗+βSι∗∗Iι∗∗−βSι∗∗2Sι−βSιIιEι∗∗Eι+l1Eι∗∗−l1Iι∗∗EιIι+βSι∗∗Iι∗∗.
We also observed from (55) at steady state that(64)l1Eι∗=βSι∗∗Iι∗∗.
Using expression (64) in equation (63) gives(65)U´=μSι∗∗2−Sι∗∗Sι−SιSι∗∗+βSι∗∗Iι∗∗3−Sι∗∗Sι−SιIιEι∗∗Sι∗∗Iι∗∗Eι−EIι∗∗IιEι∗∗.
We can simplify the first term in the last equation of (65) as follows.(66)2−Sι∗∗Sι−SιSι∗∗=2SιSι∗∗−Sι∗∗2−Sι2SιSι∗∗=−Sι−Sι∗∗2Sι−Sι∗∗≤0.
Finally, it follows that since the arithmetic mean exceeds the geometric mean,(67)3−Sι∗∗Sι−SιIιEι∗∗Sι∗∗Iι∗∗Eι−EIι∗∗IιEι∗∗≤0.
Accordingly, it follows from this thatU´≤0 when RC> with U´=0 if and only if Sι=Sι∗∗, Eι=Eι∗∗, and Iι=Iι∗∗. Thus, U is a Lyapunov function for system model (55) in Γ∖Γ0. Thus, the largest compact invariant set where U´=0 is a singleton Sι,Eι,Iι=Sι∗∗,Eι∗∗,Iι∗∗, and therefore it follows that by [39], Sιt⟶Sι∗∗, Eιt⟶Eι∗∗, and Iιt⟶Iι∗∗ as t⟶∞.


## 5. Numerical Simulation and Discussion

To confirm the feasibility of theoretical findings regarding stability and bifurcation of model system (4), we employed numerical continuation package XPPAUT [44] to perform one parameter bifurcation analysis and Matlab to perform numerical experiment. Motivated by some epidemiological literature of TB model [21, 45], model system (4) was validated by choosing a set of parameter values as illustrated in TB relevant literature, and the estimated parameters used are captured in each respective figure.Figure2 demonstrates the population of exposed individuals as the transmission rate β is varied from 0 to 2. The y-axis defines the exposed individuals and x-axis defines the transmission rate beta β. The red and black lines represent stable and unstable steady states, respectively. The green (blue, respectively) cycles illustrated the occurrence of stable and unstable limit cycles in this epidemiological system. In general, this epidemiological system demonstrates rich bifurcation dynamics which are epidemiologically important. Several dynamical behaviors of the system are observed as β changes. As demonstrated in Figure 2, the system is identified by Sub⋅H,BP,LP, and LPC, which corresponded to Hopf bifurcation, transcritical bifurcation, saddle-node bifurcation, and saddle-node bifurcation of cycles. van Voorn [46] stated that “when an unstable limit cycle is born, while the equilibrium becomes stable, we have a subcritical Hopf bifurcation.” Hence, the bifurcation is subcritical. Indeed, Figure 2 reveals that for the bifurcation point, the stable equilibrium (red line) becomes unstable (black line). In addition, stable and unstable limit cycles also emerge. As the bifurcation parameter β increases, the unstable limit cycles collide with stable limit cycles via a saddle-node bifurcation of cycles LPC.Figure 2
One-parameter bifurcation diagram of exposed individualsE with transmission rate β as a bifurcation parameter. The fixed parameters are μ=0.18, Λ=4, p=0.5, κ=0.03, c = 0.02, b = 0.2, δ=0.08, and τ=2. The Hopf bifurcation points occurred at β=1.035. Sub⋅H corresponds to subcritical Hopf bifurcation, LPC corresponds to saddle-node bifurcation of cycles, LP corresponds to saddle-node bifurcation, BP corresponds to transcritical bifurcation, SLC corresponds to stable limit cycle, and ULC corresponds to unstable limit cycle.Figure3 depicts a codimension-two bifurcation graph as β and τ are changed. To examine the interaction of subcritical Hopf, saddle-node, and transcritical bifurcations, we conducted a codimension-two bifurcation analysis of TB model (4) as we varied the transmission rate β and recovery parameter τ, and the result can be seen in Figure 3. There is indeed a codimension-2 point, black curve with red point, referring to cycle saddle-node bifurcation. This kind of bifurcation frequently occurs as a result of a collision between stable limit cycles (green curve) and unstable limit cycles (blue curve). In general, the codimension-2 bifurcation point serves as an organizing center, separating the parameter space into various curves with distinct qualitative findings: red curve corresponds to saddle-node bifurcation, blue curve corresponds to Hopf bifurcation, and cyan curve corresponds to transcritical bifurcation as demonstrated in previous result (that is, Figure 2, which illustrates the same qualitative mechanisms as Figure 3).Figure 3
Two-parameter bifurcation diagram of model system (9) in the τ,β-plane. The fixed parameters are μ=0.18, Λ=4, p=0.5, κ=0.03, c = 0.02, b = 0.2, δ=0.08, and τ=2. Red curve denotes saddle-node bifurcation, blue curve denotes Hopf bifurcation, cyan curve denotes transcritical bifurcation, and black curve at (red point) denotes saddle-node bifurcation of cycles.In Figures4 and 5, we set the parameter values as Λ=4, μ=0.18, δ=0.08, κ=0.03, p=0.5, c=0.02, and ϖ=0.1, and β is varied around the critical value βo=1.7. The threshold quantity RC=2.3670>1, and the eigenvalues are (−0.2163698605 + 1.349717265I, −3.652541848, −0.2163698605–1.349717265I) and the model system endemic equilibrium is given as (2.53865, 3.9028, 0.820965). Here, P∗ is found to be unstable. The endemic equilibrium P∗ is stable for β<βo and unstable for β>βo. We discover that the condition for Hopf bifurcation is fulfilled resulting in the presence of periodic oscillations around P∗. The graphs for the results obtained are given in Figures 4 and 5.Figure 4
The oscillatory simulation of susceptible individuals (blue colour), exposed individuals (green colour), and infected individuals (red colour). Parameters areβ=1.7, Λ=4, μ=0.18, δ=0.08, κ=0.03, p=0.5, c=0.02, and ϖ=0.1 which give RC=2.3670>1. The oscillatory solution is unstable for some initial conditions. Those used to produce the figure are S(0) = 2, E(0) = 3, and I(0) = 1.5.Figure 5
Phase trajectory for model system (4). The solution approaches periodic orbit around endemic equilibrium when RC>1. (a) Phase diagram of S, I- plane. (b) Phase diagram of E, I-plane.
(a)(b)Furthermore, we consider another set of representative parameters asΛ=4, μ=0.18, δ=0.04, κ=0.03, p=2, c=0.02, ϖ=0.1, and β=1.7. The threshold quantity RC=2.4093>1, and the equilibrium is (2.11674, 0.997134, 1.0057) while the eigenvalues are given as (−0.231843, −0.231843± 5.08514I). We can conclude that P∗ is stable since all of the eigenvalues are negative or negative real parts (see Figure 6).Figure 6
The trajectories and phase diagram ofS, I-plane of model system (4). (a) Profile of S, I with respect to time when p=0. (b) Phase trajectory for S, I-plane when p=2.
(a)(b)Figure7 depicts the phase diagram of Sι and Iι planes with different values of saturated recovery c and with μ=0.18, Λ=4, p=0.5, κ=0.03, b = 0.1, and τ=2. When c=0.04, the endemic equilibrium is unstable, and as we increase the value of c to 0.1, the endemic equilibrium becomes stable; furthermore, we observed the appearance of two (5) stable limit cycles. We observed that the stability endemic equilibrium state changes from stable to unstable, and this kind of bifurcation is called subcritical Hopf bifurcation (see [46, 47]).Figure 7
The phase diagram ofS and I planes with different values of saturated recovery c. Here, there are two limit cycles. The inner one is unstable and the outer one is stable.Figure8 illustrates a typical bifurcation graphs of system model (4). In this diagram, the recruitment rate is changed, whereas other parameters are kept constant. The set of representative parameter values which result in Figure 8 are μ=0.18, δ=0.12, κ=0.03, β=1.6, p=0.12, c=0.9, ϖ=0.11, and 3.5≤Λ≤5. Figure 8(a) demonstrates the scenario of forward bifurcation that when RC<1, the TB-free equilibrium is globally asymptotically stable, while if RC>1, the TB infection can persist. Conversely, as we observed from Figure 8(b), increasing the value of the parameter ϖ from 0.1 to 0.25, the TB infection can persist once defined for the range of RC values below unity which suggests the occurrence of the phenomenon of backward bifurcation. This suggests that diminishing RC below one will not be fundamentally enough for the eradication of TB epidemic from the population. When RC is adequately reduced such that RC<RC∗, the positive equilibrium does not exist anymore and TB epidemic will cease to develop and will gradually fall from its generally high endemic to the TB-free equilibrium. From Figure 8(b), we observed that when RC∗<RC<1, there are stable endemic equilibrium, an unstable endemic equilibrium, and a stable TB-free equilibrium, whereas when RC<1, there is just one steady endemic equilibrium.Figure 8
(a) Forward bifurcation and (b) backward bifurcation of endemic equilibria where the infected individuals are plotted as a function ofRC. The blue dash line represents the stable equilibrium while the red dash line represents the unstable equilibrium. (a) shows forward bifurcation at the critical value RC=1, and the parameter values used include μ=0.18, δ=0.12, κ=0.03, β=1.6, p=0.12, c=0.9, ϖ=0.1, and 3.5≤Λ≤5, while (b) represents backward bifurcation and parameters used are the same as in (a) except ϖ=0.25.
(a)(b)Figure9 illustrates the backward bifurcation diagram of system model (4) as the effective basic reproduction number RC varies against the infected population I. Figures 9(a)–9(d) demonstrate that increasing the value of ϖ contributes to the extension of the region of bistability, while reducing the value of ϖ brings about contraction of the bistability region. The TB eradication quantity also known as critical reproduction number shifts from left to the right when ϖ decreases and vice versa. Here, high value of ϖ means insufficient treatment for a huge population of TB infections; therefore, we prefer a scenario where there will be consistent TB epidemics within the population despite the fact that RC<1.Figure 9
Demonstrations of the impact of varying the parameter that accounts for delaying the treatment of infected individuals. Here the parameters are the same as in Figure8 except c = 1.2 while ϖ is shown on top of the figure. (a) Forward bifurcation. (b) Backward bifurcation. (c) Forward bifurcation. (d) Backward bifurcation.
(a)(b)(c)(d)


## 6. Conclusions

Dynamical models of TB infection have been analyzed by numerous authors [21, 22, 48–50]. However, they all proposed linear recovery function to investigate transmission dynamic behavior of their model. However, we discovered that recovery from TB infection depends upon numerous factors such as antibiotic treatment and the number of hospital beds. This leads to nonlinearity in the number of recoveries. The nonlinear recovery function rate principle was then introduced to reveal some insight into the eradication of TB disease by examining the effect of antibiotic treatment and the number of hospital beds.The main goal of this study is to investigate the qualitative dynamics of the TB infection incorporating saturated recovery (treatment) of the formcIι/1+ϖIι into the model proposed in [21, 22] which gives a more realistic model. On the other hand, complicated dynamic behaviors, including oscillations, can be caused by saturated recovery [9, 10, 14–17].Typically in epidemiology, the criterion for the persistence and extinction of an infection is vital. The effective reproduction number is computed and the conditions for the local stability of equilibria and the existence of backward bifurcation and Hopf bifurcation are determined. Under the condition ofRC<1, the infection disappears completely from Theorem 3. But if RC>1, then according to Theorem 3, there will be endemic disease. In reality, the epidemic cannot be eliminated except that the basic reproduction number is reduced under a lower level to such an extent that RC∗<RC<1 (see Figure 8(b)). The mathematical analysis of the TB model system (1) has shown two important mathematical phenomena: bistability and periodicity (oscillatory), which have been observed in some other infectious diseases. In the bistability phenomenon, which usually comes with a conversation of backward bifurcation (see [21] and the references therein), the system exhibits multiple endemic equilibria despite RC<1. In such situation, a stable endemic equilibrium has been shown to compete with a stable disease-free equilibrium. From these results, it has been confirmed that decreasing RC<1 cannot control the spread of TB infection in the population. The second phenomenon undergoes oscillatory behavior under some certain conditions. The occurrence of limit cycles supports the behavior stated in several studies on the dynamics of some contagious diseases such as whooping cough, measles, rubella, and so on [51–54]. Numerical findings indicate that ϖ which is the saturation parameter is accountable for backward bifurcation. Inability to intervene before TB infection has accumulated in the population will prompt a circumstance where a TB epidemic exists despite the fact that the RC<1. Improving existing medical technologies and channeling adequate resources in medicines can essentially encourage early intervention by ensuring that individuals infected with TB get treatment immediately.



---

*Source: 1805651-2021-10-23.xml*