# Mathematical Analysis of TB Model with Vaccination and Saturated Incidence Rate

Authors: Ashenafi Kelemu Mengistu; Peter J. Witbooi
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

The model system of ordinary differential equations considers two classes of latently infected individuals, with different risk of becoming infectious. The system has positive solutions. By constructing a Lyapunov function, it is proved that if the basic reproduction number is less than unity, then the disease-free equilibrium point is globally asymptotically stable. The Routh-Hurwitz criterion is used to prove the local stability of the endemic equilibrium whenR0>1. The model is illustrated using parameters applicable to Ethiopia. A variety of numerical simulations are carried out to illustrate our main results.


---

## Body

## 1. Introduction

Tuberculosis (TB) is a contagious disease caused byMycobacterium tuberculosis and is transmitted from person to person through the air. Even though TB disease can be controlled and cured, it continues to be a major health problem in the world and is now one of the top ten causes of death and the leading cause of a single infectious agent (above HIV/AIDS). In 2019, approximately 10 million people developed TB and 1.4 million died [1]. The prevention of TB involves vaccination with Bacillus Calmette-Guérin (BCG), screening of those at high risk, early detection, and treatment of cases [2, 3].BCG vaccination is one of the strategies to control the transmission of TB. In most TB endemic countries BCG vaccination is recommended for tuberculosis prevention and is usually given shortly after birth to prevent TB in infants [3, 4]. The incidence rate plays a very important role in the research of epidemiological models. In many epidemic models, bilinear incidence rate βSI as in the classical Kermack and McKendrick model [5] is frequently used. Many models of TB also use this type of incidence rate (see [6–10]). Another kind of incidence rate proposed by Capasso and Serio [11] is the saturation incidence rate βSI/1+bI. Here, β𝐼 measures the infection force when the disease is entering a fully susceptible population, and 1/1+bI measures the inhibition effect from the behavioral change of susceptible individuals when their number increases or from the crowding effect of the infective individuals [12, 13]. The saturation incidence rate is more reasonable than the bilinear saturating incidence rate when we need to include the behavioral change and crowding effect of the infective individuals, in order to curb the contact rate [13, 14].Mathematical modeling fulfills a significant role to examine, explain, and predict the dynamics of infectious disease transmission, including tuberculosis [8, 15, 16]. Ongoing research is aimed at developing more realistic mathematical models for investigating the transmission dynamics of infectious diseases. One of the main issues in mathematical epidemiology is the study of the asymptotic behavior of epidemic models, and for this purpose, we need to analyze steady states and their stability [17]. In this paper, we present and analyze a basic tuberculosis mathematical model with a saturated incidence rate and TB vaccination.


## 2. Construction of the Model

Based on the disease status of individuals, we divide the total populationNt into four subclasses, namely, susceptible St, high-risk latently infected Et, infectious (or active TB) It, and low-risk latent Lt. For the model framework, we consider the following assumptions:
(i)
The rate at which new individuals enter into the susceptible class due to birth is denoted byΛ(ii)
Following [17–19], the disease transmission rate is considered as βI/1+bI, where b is the saturation rate and β is the maximum contact rate between susceptible and infected individuals(iii)
The BCG vaccine will be given to the susceptible population of a rateεS, where 0≤ε≤1. The BCG vaccine efficacy in preventing adults from TB is incomplete, with an average efficacy of 50% [20]. This shows that vaccinated individuals may still be vulnerable to bacteria. Hence, it is reasonable to classify the vaccinated and nonvaccinated individuals into a single class with a different chance of infection. It is assumed that the chance of being infected by the bacteria for the vaccinated and nonvaccinated population is θβεSI/1+bI and βS1−εI/1+bI, respectively, where θ, with 0≤θ≤1, is the loss of immunity for a vaccinated person(iv)
k is the rate at which individuals in the high-risk latent class E becomes infective(v)
It is assumed that treatment will be administered for both high-risk latent and active TB-infected classes. The treatment rate for theE-class and I-class is denoted by α and r, respectively(vi)
If treatment is administered for theI-class with a rate r, then some of them will complete their treatment and recover completely at a rate 0≤p≤1 and move to the L-class. Others, 1−prI, will not be cured and will remain vulnerable to the bacteria and move to the E-class(vii)
Patients who have completed anti-TB treatment will recover, but these individuals may remain latent because the TB bacteria stay dormant in the host body. Accordingly, we classify these individuals as low-risk latent(viii)
After being cured, it is assumed that some of the recovered (low-risk latent) individuals can be reinfected with the bacteria and become high-risk latently infected with the rateσ(ix)
We assume that the natural death rate (any death which is not due to TB) is the same for all classes and is denoted byμ(x)
Death due to the TB disease will happen only toI-class, and we denote this mortality rate as δ(xi)
We further assume that all parameters to be used in this model are nonnegativeBased on the above assumptions, we describe the dynamics of TB by the following system of ordinary differential equations (ODEs) with four compartments.(1)dSdt=Λ−1−ε+θεβSI1+bI−μS,dEdt=1−ε+θεβSI1+bI+1−prΙ+σL−k+α+μE,dIdt=kE−μ+r+δΙ,dLdt=prΙ+αE−μ+σL,Nt=St+Et+It+Lt,with the initial conditions S0,E0,I0,L0≥0.The complete transfer flow of the model parameters is shown in Figure1. The time unit “t” has been considered in “years.”Figure 1
Flow diagram of the TB transmission model.


### 2.1. Basic Properties of the Model

#### 2.1.1. Positivity of Solutions

Since the model system (1) involves the human population, it is necessary to prove that all its associated variables and parameters are nonnegative, so we have the following theorem.Theorem 1.
Let the initial dataS0, E0, I0, and L0 be nonnegative. Then, the solution set Ω=St,Et,It,Lt is nonnegative for all t≥0.Proof.
If we let1−ε+θεβIt/1+bIt−μ=Ht, then it follows from the first equation of the model (1) that
(2)dStdt+HtSt=Λ.Multiplying both sides of (2) by exp∫0tHτdτ gives
(3)dStdtexp∫0tHτdτ+HtStexp∫0tHτdτ=Λexp∫0tHτdτ.By the product rule of the derivative, we have(4)dStdtexp∫0tHτdτ+HtStexp∫0tHτdτ=ddtStexp∫0tHτdτ.Hence,(5)ddtStexp∫0tHτdτ=Λexp∫0tHτdτ.Integrating both sides of (5) gives
(6)Stexp∫0tHτdτ−S0=Λ∫0texp∫0τHudu.Then,(7)St=S0exp−∫0tHτdτ+Λ∫0texp∫0τHuduexp−∫0tHτdτ≥0.Similarly, it can be shown thatEt, It, and Lt are nonnegative for all time t≥0.This completes the proof.


#### 2.1.2. Invariant Region

The TB model (1) will be studied in a biologically feasible region as given below.
(8)Ω=St,Et,It,Lt∈ℝ+4∣0≤St+Et+It+Lt≤Λμ.Lemma 2.
The regionΩ⊂ℝ+4 is positively invariant for the model (1) with nonnegative initial conditions in ℝ+4.
Proof. For this model, the total population is given by(9)Nt=St+Et+It+Lt.Then,(10)dNtdt=dStdt+dEtdt+dItdt+dLtdt.Substituting the values from the model (1), we get
(11)dNtdt=Λ−μNt−δΙt≤Λ−μNt.Thus, integrating both sides and taking ast⟶∞, we get 0≤Nt≤Λ/μ.Therefore, the feasible solution set for the model is given by(12)Ω=St,Et,It,Lt∈ℝ+4∣0≤St+Et+It+Lt≤Λμ.


## 2.1. Basic Properties of the Model

### 2.1.1. Positivity of Solutions

Since the model system (1) involves the human population, it is necessary to prove that all its associated variables and parameters are nonnegative, so we have the following theorem.Theorem 1.
Let the initial dataS0, E0, I0, and L0 be nonnegative. Then, the solution set Ω=St,Et,It,Lt is nonnegative for all t≥0.Proof.
If we let1−ε+θεβIt/1+bIt−μ=Ht, then it follows from the first equation of the model (1) that
(2)dStdt+HtSt=Λ.Multiplying both sides of (2) by exp∫0tHτdτ gives
(3)dStdtexp∫0tHτdτ+HtStexp∫0tHτdτ=Λexp∫0tHτdτ.By the product rule of the derivative, we have(4)dStdtexp∫0tHτdτ+HtStexp∫0tHτdτ=ddtStexp∫0tHτdτ.Hence,(5)ddtStexp∫0tHτdτ=Λexp∫0tHτdτ.Integrating both sides of (5) gives
(6)Stexp∫0tHτdτ−S0=Λ∫0texp∫0τHudu.Then,(7)St=S0exp−∫0tHτdτ+Λ∫0texp∫0τHuduexp−∫0tHτdτ≥0.Similarly, it can be shown thatEt, It, and Lt are nonnegative for all time t≥0.This completes the proof.


### 2.1.2. Invariant Region

The TB model (1) will be studied in a biologically feasible region as given below.
(8)Ω=St,Et,It,Lt∈ℝ+4∣0≤St+Et+It+Lt≤Λμ.Lemma 2.
The regionΩ⊂ℝ+4 is positively invariant for the model (1) with nonnegative initial conditions in ℝ+4.
Proof. For this model, the total population is given by(9)Nt=St+Et+It+Lt.Then,(10)dNtdt=dStdt+dEtdt+dItdt+dLtdt.Substituting the values from the model (1), we get
(11)dNtdt=Λ−μNt−δΙt≤Λ−μNt.Thus, integrating both sides and taking ast⟶∞, we get 0≤Nt≤Λ/μ.Therefore, the feasible solution set for the model is given by(12)Ω=St,Et,It,Lt∈ℝ+4∣0≤St+Et+It+Lt≤Λμ.


## 2.1.1. Positivity of Solutions

Since the model system (1) involves the human population, it is necessary to prove that all its associated variables and parameters are nonnegative, so we have the following theorem.Theorem 1.
Let the initial dataS0, E0, I0, and L0 be nonnegative. Then, the solution set Ω=St,Et,It,Lt is nonnegative for all t≥0.Proof.
If we let1−ε+θεβIt/1+bIt−μ=Ht, then it follows from the first equation of the model (1) that
(2)dStdt+HtSt=Λ.Multiplying both sides of (2) by exp∫0tHτdτ gives
(3)dStdtexp∫0tHτdτ+HtStexp∫0tHτdτ=Λexp∫0tHτdτ.By the product rule of the derivative, we have(4)dStdtexp∫0tHτdτ+HtStexp∫0tHτdτ=ddtStexp∫0tHτdτ.Hence,(5)ddtStexp∫0tHτdτ=Λexp∫0tHτdτ.Integrating both sides of (5) gives
(6)Stexp∫0tHτdτ−S0=Λ∫0texp∫0τHudu.Then,(7)St=S0exp−∫0tHτdτ+Λ∫0texp∫0τHuduexp−∫0tHτdτ≥0.Similarly, it can be shown thatEt, It, and Lt are nonnegative for all time t≥0.This completes the proof.


## 2.1.2. Invariant Region

The TB model (1) will be studied in a biologically feasible region as given below.
(8)Ω=St,Et,It,Lt∈ℝ+4∣0≤St+Et+It+Lt≤Λμ.Lemma 2.
The regionΩ⊂ℝ+4 is positively invariant for the model (1) with nonnegative initial conditions in ℝ+4.
Proof. For this model, the total population is given by(9)Nt=St+Et+It+Lt.Then,(10)dNtdt=dStdt+dEtdt+dItdt+dLtdt.Substituting the values from the model (1), we get
(11)dNtdt=Λ−μNt−δΙt≤Λ−μNt.Thus, integrating both sides and taking ast⟶∞, we get 0≤Nt≤Λ/μ.Therefore, the feasible solution set for the model is given by(12)Ω=St,Et,It,Lt∈ℝ+4∣0≤St+Et+It+Lt≤Λμ.


## 3. Equilibria

The long-term behavior of a model is quite important. When left for a long time with no external interference, we want to know whether there will eventually be a stable state to which the system converges. If the disease consistently vanishes, we say that the system converges to adisease-free equilibrium. It is also possible that the disease persists, and the class sizes tend to converge to a single stable point, an endemic equilibrium point.


### 3.1. Disease-Free Equilibrium and the Basic Reproduction Number

It is easy to check that model (1) always has a disease-free equilibrium (DFE), P0∗=S0∗,0,0,0, where S0∗=Λ/μ. Now, we introduce the basic reproduction number R0, which is defined as the expected average number of new TB infections caused by a single infected individual during his/her infected time when in contact with a completely susceptible population. The basic reproduction number plays a key role in determining the global dynamics of the disease in an epidemiological study. It helps us to know whether the disease will spread and persist or whether it will eventually vanish from the population. If the value of the reproduction number is less than one, it indicates that the disease is not spreading in the community, whereas when its value is greater than one, then the disease can invade the community, i.e., the disease becomes endemic and will produce deaths and perhaps an outbreak. We obtained R0 by using the next-generation matrix method given in [21].LetX=E,I,LT, then it follows from the system (1) that
(13)dXdt=F−M,where
(14)F=1−ε+θεβSI1+bI00,M=−1−prΙ−σL+k+α+μE−kE+μ+r+δΙ−prΙ−αE+μ+σL.Evaluating the derivatives ofF and M at DFE leads to the following matrices F and
(15)F=0βΛ1−ε+θεμ0000000,M=k+α+μ−1+pr−σ−kr+δ+μ0−α−prμ+σ.Then,R0 is the dominant eigenvalue of the matrix FM−1.Thus, we have(16)R0=kβΛ1−ε+θεμ+σμμr+δ+μα+μ+σ+kprμ+δ+μμ+σ.


### 3.2. Global Stability of the DFE

Global asymptotic stability of an equilibrium point means that starting from any initial state, the system will eventually converge to that particular equilibrium point.Theorem 3 [22].
For model (1), the disease-free equilibrium point P0∗ is globally asymptotically stable if R0≤1.Proof.
To establish the global stability of the disease-free equilibrium, we construct the following Lyapunov function:(17)T=μ+σE+kμ+σ+μα+μ+σkI+σL.Note thatT is nonnegative and TP0∗=0. Differentiation with respect to t and using (1) leads to
(18)Ṫ=μ+σĖ+kμ+σ+μα+μ+σkİ+σL̇=μ+σ1−ε+θεβSI1+bI+1−prΙ+σL−k+α+μΕ+kμ+σ+μα+μ+σkkE−μ+r+δΙ+σprΙ+αΕ−μ+σL.SinceβSI/1+bI≤βSI, we have
(19)Ṫ≤μ+σ1−ε+θεβSI+1−prΙ+σL−k+α+μΕ+kμ+σ+μα+μ+σkkE−μ+r+δΙ+σprΙ+αΕ−μ+σL.UsingSt≤Λ/μ and after simplification, we have
(20)Ṫ≤βΛμ+σ1−ε+θεμ+μ+σ1−pr−kμ+σ+μα+μ+σkμ+r+δ+σprI.Finally, with some rearrangement, we obtain(21)Ṫ≤μr+δ+μα+μ+σ+kprμ+δ+μμ+σkR0−1I.Since all parameters of the model are nonnegative, it follows thatṪ≤0 for R0≤1and Ṫ=0 only if I=0. Hence,T is a Lyapunov function on Ω and the largest compact invariant subset of St,Et,It,Lt∈Ω:T=0 is the singleton P0∗. La Salle’s invariant principle [22] implies that P0∗ is globally asymptotically stable in Ω.


### 3.3. The Endemic Equilibrium Point (EE)

Lemma 4. (existence of EE).
IfR0>1, then the model (1) has a unique positive endemic equilibrium point.Proof.
The EE of the model (1) atP∗=S∗,E,I∗,L∗ is obtained by setting the right-hand side of the system (1) equal to zero, we get
(22)S∗=μr+δ+μα+μ+σ+kprμ+kδ+bΛ+μμ+σkβ1−ε+θεμ+σ+kbμμ+σ,Ε∗=μr+δ+μkβ1−ε+θε+bμR0−1,I∗=μβ1−ε+θε+bμR0−1,L∗=μkrp+αr+δ+μkβ1−ε+θε+bμμ+σR0−1.Clearly, from the above, we can conclude that a unique positive EE exists wheneverR0>1.


### 3.4. Local Stability of the EE

Local asymptotic of an equilibrium point stability means that if the system is very near to that equilibrium point, then it will necessarily converge to the equilibrium state. It also means that minor perturbations will not disrupt this convergence.Theorem 5 (local stability at EE).
IfR0>1 and σ=0, then the system (1) is locally asymptotically stable about the endemic equilibrium point P∗.Proof.
Forσ=0, the variable L will only appear in the fourth equation of the model (1), hence the system can be reduced to the three-dimensional system:
(23)dSdt=Λ−1−ε+θεβSI1+bI−μS,dEdt=1−ε+θεβSI1+bI+1−prΙ−k+α+μΕ,dIdt=kE−μ+r+δΙ.By applying the Routh-Hurwitz criterion of stability [23] and by denoting ψ=1−ε+θεβ for the system (23), we have the following Jacobian matrix at P∗:
(24)J∗=−μ−βψI∗1+bI∗0−βψS∗1+bI∗2βψI∗1+bI∗−k−α−μr−pr+βψS∗1+bI∗20k−r−δ−μ.The associated characteristic equation ofJ∗ is given by
(25)λ3+a1λ2+a2λ+a3=0.where
(26)a1=k+r+α+δ+3μ+βψI∗1+bI∗,a2=kβ2Λk+r+α+δ+2μψ2+bμα+μr+δ+μ+kpr+δ+μ2R0−1+kβΛψk+r+α+δμ+2μ2ψβα+μr+δ+μ+kpr+δ+bΛ+μψ,a3=μbμ+βψα+μr+δ+μ+kpr+δ+μ2R0−1βα+μr+δ+μ+kpr+δ+bΛ+μψ,a1a2−a3=R0−1μα+μr+δ+μ+kpr+δ+μ2α+μr+δ+μ+kpr+δ+bΛ+μbμ+βψβα+μr+δ+μ+kpr+δ+bΛ+μ2ψ+Akβ2Λk+r+α+δ+2μψ2+bAR0−1μα+μr+δ+μ+kpr+δ+μ2βα+μr+δ+μ+kpr+δ+bΛ+μ2ψ+bAkβk+r+α+δΛψ+2μ2ψβα+μr+δ+μ+kpr+δ+bΛ+μ2ψ,where
(27)A=k2pr+δ+bΛ+μ+α+μr+δ+μr+α+δ+2μ+krα+δ+bΛ+2μ+βΛψ+k2αδ+δ2+bαΛ+bδΛ+2αμ+4δμ+3bΛμ+3μ2+prr+α+δ+2μ.Clearly,a1, a2,a3, and a1a2−a3 are positive whenever R0>1. Thus, the Routh-Hurwitz criterion is satisfied. Therefore, the endemic equilibrium point P∗ of the system (23) is locally asymptotically stable for R0>1. This completes the proof.


## 3.1. Disease-Free Equilibrium and the Basic Reproduction Number

It is easy to check that model (1) always has a disease-free equilibrium (DFE), P0∗=S0∗,0,0,0, where S0∗=Λ/μ. Now, we introduce the basic reproduction number R0, which is defined as the expected average number of new TB infections caused by a single infected individual during his/her infected time when in contact with a completely susceptible population. The basic reproduction number plays a key role in determining the global dynamics of the disease in an epidemiological study. It helps us to know whether the disease will spread and persist or whether it will eventually vanish from the population. If the value of the reproduction number is less than one, it indicates that the disease is not spreading in the community, whereas when its value is greater than one, then the disease can invade the community, i.e., the disease becomes endemic and will produce deaths and perhaps an outbreak. We obtained R0 by using the next-generation matrix method given in [21].LetX=E,I,LT, then it follows from the system (1) that
(13)dXdt=F−M,where
(14)F=1−ε+θεβSI1+bI00,M=−1−prΙ−σL+k+α+μE−kE+μ+r+δΙ−prΙ−αE+μ+σL.Evaluating the derivatives ofF and M at DFE leads to the following matrices F and
(15)F=0βΛ1−ε+θεμ0000000,M=k+α+μ−1+pr−σ−kr+δ+μ0−α−prμ+σ.Then,R0 is the dominant eigenvalue of the matrix FM−1.Thus, we have(16)R0=kβΛ1−ε+θεμ+σμμr+δ+μα+μ+σ+kprμ+δ+μμ+σ.


## 3.2. Global Stability of the DFE

Global asymptotic stability of an equilibrium point means that starting from any initial state, the system will eventually converge to that particular equilibrium point.Theorem 3 [22].
For model (1), the disease-free equilibrium point P0∗ is globally asymptotically stable if R0≤1.Proof.
To establish the global stability of the disease-free equilibrium, we construct the following Lyapunov function:(17)T=μ+σE+kμ+σ+μα+μ+σkI+σL.Note thatT is nonnegative and TP0∗=0. Differentiation with respect to t and using (1) leads to
(18)Ṫ=μ+σĖ+kμ+σ+μα+μ+σkİ+σL̇=μ+σ1−ε+θεβSI1+bI+1−prΙ+σL−k+α+μΕ+kμ+σ+μα+μ+σkkE−μ+r+δΙ+σprΙ+αΕ−μ+σL.SinceβSI/1+bI≤βSI, we have
(19)Ṫ≤μ+σ1−ε+θεβSI+1−prΙ+σL−k+α+μΕ+kμ+σ+μα+μ+σkkE−μ+r+δΙ+σprΙ+αΕ−μ+σL.UsingSt≤Λ/μ and after simplification, we have
(20)Ṫ≤βΛμ+σ1−ε+θεμ+μ+σ1−pr−kμ+σ+μα+μ+σkμ+r+δ+σprI.Finally, with some rearrangement, we obtain(21)Ṫ≤μr+δ+μα+μ+σ+kprμ+δ+μμ+σkR0−1I.Since all parameters of the model are nonnegative, it follows thatṪ≤0 for R0≤1and Ṫ=0 only if I=0. Hence,T is a Lyapunov function on Ω and the largest compact invariant subset of St,Et,It,Lt∈Ω:T=0 is the singleton P0∗. La Salle’s invariant principle [22] implies that P0∗ is globally asymptotically stable in Ω.


## 3.3. The Endemic Equilibrium Point (EE)

Lemma 4. (existence of EE).
IfR0>1, then the model (1) has a unique positive endemic equilibrium point.Proof.
The EE of the model (1) atP∗=S∗,E,I∗,L∗ is obtained by setting the right-hand side of the system (1) equal to zero, we get
(22)S∗=μr+δ+μα+μ+σ+kprμ+kδ+bΛ+μμ+σkβ1−ε+θεμ+σ+kbμμ+σ,Ε∗=μr+δ+μkβ1−ε+θε+bμR0−1,I∗=μβ1−ε+θε+bμR0−1,L∗=μkrp+αr+δ+μkβ1−ε+θε+bμμ+σR0−1.Clearly, from the above, we can conclude that a unique positive EE exists wheneverR0>1.


## 3.4. Local Stability of the EE

Local asymptotic of an equilibrium point stability means that if the system is very near to that equilibrium point, then it will necessarily converge to the equilibrium state. It also means that minor perturbations will not disrupt this convergence.Theorem 5 (local stability at EE).
IfR0>1 and σ=0, then the system (1) is locally asymptotically stable about the endemic equilibrium point P∗.Proof.
Forσ=0, the variable L will only appear in the fourth equation of the model (1), hence the system can be reduced to the three-dimensional system:
(23)dSdt=Λ−1−ε+θεβSI1+bI−μS,dEdt=1−ε+θεβSI1+bI+1−prΙ−k+α+μΕ,dIdt=kE−μ+r+δΙ.By applying the Routh-Hurwitz criterion of stability [23] and by denoting ψ=1−ε+θεβ for the system (23), we have the following Jacobian matrix at P∗:
(24)J∗=−μ−βψI∗1+bI∗0−βψS∗1+bI∗2βψI∗1+bI∗−k−α−μr−pr+βψS∗1+bI∗20k−r−δ−μ.The associated characteristic equation ofJ∗ is given by
(25)λ3+a1λ2+a2λ+a3=0.where
(26)a1=k+r+α+δ+3μ+βψI∗1+bI∗,a2=kβ2Λk+r+α+δ+2μψ2+bμα+μr+δ+μ+kpr+δ+μ2R0−1+kβΛψk+r+α+δμ+2μ2ψβα+μr+δ+μ+kpr+δ+bΛ+μψ,a3=μbμ+βψα+μr+δ+μ+kpr+δ+μ2R0−1βα+μr+δ+μ+kpr+δ+bΛ+μψ,a1a2−a3=R0−1μα+μr+δ+μ+kpr+δ+μ2α+μr+δ+μ+kpr+δ+bΛ+μbμ+βψβα+μr+δ+μ+kpr+δ+bΛ+μ2ψ+Akβ2Λk+r+α+δ+2μψ2+bAR0−1μα+μr+δ+μ+kpr+δ+μ2βα+μr+δ+μ+kpr+δ+bΛ+μ2ψ+bAkβk+r+α+δΛψ+2μ2ψβα+μr+δ+μ+kpr+δ+bΛ+μ2ψ,where
(27)A=k2pr+δ+bΛ+μ+α+μr+δ+μr+α+δ+2μ+krα+δ+bΛ+2μ+βΛψ+k2αδ+δ2+bαΛ+bδΛ+2αμ+4δμ+3bΛμ+3μ2+prr+α+δ+2μ.Clearly,a1, a2,a3, and a1a2−a3 are positive whenever R0>1. Thus, the Routh-Hurwitz criterion is satisfied. Therefore, the endemic equilibrium point P∗ of the system (23) is locally asymptotically stable for R0>1. This completes the proof.


## 4. Sensitivity Analysis ofR0

In this subsection, we investigate the effect of the various model parameters on the basic reproduction number,R0. For this purpose, we used the normalized forward sensitivity index also known as elasticity [24]. If we let R0 to be a differentiable function with respect to π, then the normalized forward sensitivity index of R0 with respect to the parameter π is given by
(28)ΓπR0=∂R0∂π×πR0.This normalized sensitivity index measures the relative change ofR0 with respect to π.To study the sensitivity ofR0, we choose the parameters β, α, r, ε, andp. Using the above formula (28) and taking the parameters’ value in Table 1, we calculate the sensitivity indices of the basic reproduction number with respect to these parameters.Table 1
List of parameter values used for the model simulations.
ParametersDescriptionValueN0Initial total population7.25×107S0The initial number of susceptible individuals3.85×107E0The initial number of high-risk latent individuals1.19×107I0The initial number of infectious individuals3.73×105L0The initial number of low-risk latent individuals2.18×107ΛRecruitment rate1.4×106βThe transmission coefficient1.646×10−7εVaccination coverage rate0.715θLose of protection for vaccination0.5μThe natural death rate0.016kProgression rate fromE to I0.023rThe treatment rate ofI0.546pSuccessful treatment rate ofI0.832αTreatment rate ofE0.153δTB induced death rate0.17σThe relapse rate0.0013bSaturation rate0.0004 (assuming)The positive (negative) values indicate a positive (negative) correlation withR0, whereas the magnitude determines the importance of parameter [25]. From Table 2, we can see that β has a positive sensitivity index value, implying that a positive change on this parameter will increase the number of total infected population (E+I). In contrast α, r, ε, and p have negative sensitivity index values, and thus, raising these parameters will consequently decrease the number of the total infected population. The effect of the transmission rate β has the largest influence on R0 while the effect of successful treatment rate for the I class p has the smallest impact to R0. Further, this implies that increasing the value of β by 10% will increase the basic reproduction number by 10%. On the other hand, the parameters α, r, ε, and p have negative influence; hence, increasing these parameters by 10% will decrease the basic reproduction number by 8.015%, 7.402%, 5.564%, and 0.748%, respectively.Table 2
Sensitivity indices ofR0 with respect to model parameters.
ParametersDescriptionSensitivity index ofR0βThe transmission coefficient+1αTreatment rate ofE−0.801497rThe treatment rate ofI−0.74021εVaccination coverage rate−0.55642pSuccessful treatment rate ofI−0.0747728


## 5. Numerical Simulation and Discussion

To illustrate some of the analytic results obtained above, we give numerical simulations using the parameter values in Table1. In our previous study [7], we estimated the value of parameters of the TB model based on Ethiopian data, and in the current study, we used those parameters’ values. In particular, this means that the simulations are applicable to Ethiopia.Figure2 presents the dynamics of the model for different initial conditions. It shows that only the susceptible population (S∗=8.7×107) persists but the high-risk latent E, infective I, and the low-risk latent L decline to zero. It shows that system (1) is globally asymptotically stable at the DFE whenever R0≤1, which supports our analytical results stated in Theorem 3. Similarly, in Figure 3, for R0>1, the solution curves of the model are plotted by varying the initial values of the compartments, and it tends to the endemic equilibrium point. This confirms that P∗ is locally asymptotically stable, supporting the conclusion of Theorem 5.Figure 2
The plot shows the convergence of the solution of the model to the disease-free equilibrium pointP0∗=8.7×107,0,0,0 for different initial values of S,E,I,L when β=1.65×10−8 and R0=0.164<1.Figure 3
The plot shows for different initial values ofS,E,I,L, the solution of the model converges to the endemic equilibrium point P∗=8.5×107,0.15×106,0.048×105,1.5×106 forβ=3.02×10−7 and R0=3.0014>1.Figure4 is a graphical representation of the components of the endemic equilibrium point, P∗, and shows the changes in the susceptible individuals, high-risk latent, infectious, and low-risk latent classes as the reproduction number, R0, is varied. In Figure 4(a), the susceptible individuals are being depleted rapidly as R0 becomes large. Figures 4(b)–4(d) show that the number of high-risk latent, infectious, and low-risk latent individuals has linear relationships with the reproduction number, R0, and all become large as it becomes large.Figure 4
Graphs showing the behavior of the susceptibleS, high-risk latent E, infectious I, and low-risk latent L classes as the reproduction number gets larger.
(a)(b)(c)(d)Generally, Figures2–4 shows the role of reproduction number to determine the dynamics of TB. It is shown that if R0<1, the system appears in a disease-free state, that is, TB will be eradicated ultimately from the population. On the other hand, when the reproduction number is higher than unity, the system will persist in the endemic state, that is, TB will spread in the population.Figure5 explains the effect of the transmission coefficient β on the number of high-risk exposed and active TB population. It shows that an increase in the TB transmission rate increases the number of both high-risk exposed and active TB population. Figure 6 shows the simulation of active TB and high-risk exposed populations for different values of α. Thus, it can be seen that increasing the treatment rate for high-risk latent TB patients helps to reduce the total number of both high-risk exposed and active TB individuals. From Figure 7, by increasing the treatment rate for active TB patients, the total number of high-risk exposed and active TB individuals decreases. Similarly, Figure 8 shows the simulation of active TB and high-risk exposed populations for different values of ε. As we can see from the figure, increasing the coverage of the BCG vaccine for newborns can reduce the total number of the infected population. Finally, in Figure 9, the effect of successful treatment rate p is shown. The total infected individuals decrease when p increases.Figure 5
Simulation of high-risk exposed and active TB-infected population with different values ofβ.Figure 6
Simulation of high-risk exposed and active TB-infected population with different values ofα.Figure 7
Simulation of high-risk exposed and active TB-infected population with different values ofr.Figure 8
Simulation of high-risk exposed and active TB-infected population with different values ofε.Figure 9
Simulation of high-risk exposed and active TB-infected population with different values ofp.


## 6. Conclusion and Future Directions

In this paper, we investigate the dynamical behavior of a tuberculosis disease model with vaccination and a saturated incidence rate. In our model, we have divided the total population into four compartments, namely, susceptible, high-risk latent, infective, and low-risk latent population. For the proposed model, two equilibria, the disease-free equilibrium and endemic equilibrium, are derived. We have found the basic reproduction number, which helps us to determine the dynamical behavior of the model. We obtain that the disease-free equilibrium point of the model (1) is globally asymptotically stable when R0≤1. On the other hand, if R0>1, then the endemic equilibrium point is locally asymptotically stable.From the numerical simulations, the number of high-risk exposed class and active TB-infected population has a linear relationship with the reproduction number and all become large asR0 becomes large. On the other hand, the basic reproduction number depends on the transmission rate β, treatment rate of high-risk latent α, treatment rate of active TB r, vaccination coverage rate ε, and successful treatment rate of active TB p. Therefore, it is important to identify the best strategies to control and eradicate the disease. Using sensitivity analysis in this study, we found that the first effective way to control the spread of tuberculosis in Ethiopia is to minimize contact between TB-infected and susceptible individuals. The second important strategy is to increase access to treatment for latently infected individuals. Therefore, we recommend that the following strategies be implemented first and foremost to control tuberculosis and eradicate it from Ethiopia. These strategies include early detection and isolation of infectious people, conducting health campaigns and educating the community, screening of high-risk exposed individuals, and the treatment of latent TB. At the same time, the BCG vaccine plays an important role in the prevention of the disease, and immunization of the children should be continued.Many TB model studies use the bilinear incidence rate (see [7, 26, 27]), while other studies have used the saturated incidence rate but have not explicitly considered the BCG vaccine in their model (e.g., [18]). Some models combine saturated incidence rate and vaccination [28], but their numerical analysis does not apply to real data. In our study, we developed a mathematical model that included vaccination and saturated incidence rates, and in the numerical analysis of the study, we used the parameters that we fitted based on the data from Ethiopia. This allows us to make the model unique and more realistic. However, due to the complex nature of tuberculosis, we can make the model more realistic if we include other parameters such as social distance and wearing facemasks. In our next work, we will improve this model by incorporating these TB prevention methods.



---

*Source: 6669997-2020-12-04.xml*