Mathematical Analysis of TB Model with Vaccination and Saturated Incidence Rate

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 RouthHurwitz criterion is used to prove the local stability of the endemic equilibrium when R0 > 1. The model is illustrated using parameters applicable to Ethiopia. A variety of numerical simulations are carried out to illustrate our main results.


Introduction
Tuberculosis (TB) is a contagious disease caused by Mycobacterium 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][7][8][9][10]). Another kind of incidence rate proposed by Capasso and Serio [11] is the saturation incidence rate βSI/ð1 + bIÞ. Here, βI 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.

Construction of the Model
Based on the disease status of individuals, we divide the total population NðtÞ into four subclasses, namely, susceptible SðtÞ, high-risk latently infected EðtÞ, infectious (or active TB) IðtÞ, and low-risk latent LðtÞ. 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][18][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 βSð1 − εÞ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 the E -class and I-class is denoted by α and r, respectively (vi) If treatment is administered for the I-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 − pÞrI, 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 to I-class, and we denote this mortality rate as δ (xi) We further assume that all parameters to be used in this model are nonnegative Based on the above assumptions, we describe the dynamics of TB by the following system of ordinary differential equations (ODEs) with four compartments.
with the initial conditions S 0 , E 0 , I 0 , L 0 ≥ 0. The complete transfer flow of the model parameters is shown in Figure 1. The time unit "t" has been considered in "years." 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 data S 0 , E 0 , I 0 , and L 0 be nonnegative. Then, the solution set Ω = fSðtÞ, EðtÞ, IðtÞ, LðtÞg is nonnegative for all t ≥ 0: Proof. If we let ðð1 − ε + θεÞβIðtÞÞ/ð1 + bIðtÞÞ − μ = HðtÞ, then it follows from the first equation of the model (1) that Multiplying both sides of (2) by exp f Ð t 0 HðτÞdτg gives  Abstract and Applied Analysis By the product rule of the derivative, we have Hence, Integrating both sides of (5) gives Then, Similarly, it can be shown that EðtÞ, IðtÞ, and LðtÞ are nonnegative for all time t ≥ 0.
This completes the proof.

Invariant
Region. The TB model (1) will be studied in a biologically feasible region as given below.
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 Then, Substituting the values from the model (1), we get Thus, integrating both sides and taking as t ⟶ ∞, we get 0 ≤ NðtÞ ≤ Λ/μ. Therefore, the feasible solution set for the model is given by

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 a disease-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.

Disease-Free Equilibrium and the Basic Reproduction
Number. It is easy to check that model (1) always has a disease-free equilibrium (DFE), P * 0 = ðS * 0 , 0, 0, 0Þ, where S * 0 = Λ/μ. Now, we introduce the basic reproduction number R 0 , 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 R 0 by using the next-generation matrix method given in [21]. Let X = ðE, I, LÞ T , then it follows from the system (1) that where

Abstract and Applied Analysis
Evaluating the derivatives of F and M at DFE leads to the following matrices F and Then, R 0 is the dominant eigenvalue of the matrix FM −1 . Thus, we have 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.
Proof. To establish the global stability of the disease-free equilibrium, we construct the following Lyapunov function: Note that T is nonnegative and TðP * 0 Þ = 0. Differentiation with respect to t and using (1) leads to Since βSI/ð1 + bIÞ ≤ βSI, we have Using SðtÞ ≤ Λ/μ and after simplification, we have Finally, with some rearrangement, we obtain Since all parameters of the model are nonnegative, it follows that _ T ≤ 0 for R 0 ≤ 1 and _ T = 0 only if I = 0. Hence, T is a Lyapunov function on Ω and the largest compact invariant subset of fðSðtÞ, EðtÞ, IðtÞ, LðtÞÞ ∈ Ω : T = 0g is the singleton P * 0 . La Salle's invariant principle [22] implies that P * 0 is globally asymptotically stable in Ω. Proof. The EE of the model (1) at P * = ðS * , E, I * , L * Þ is obtained by setting the right-hand side of the system (1) equal to zero, we get

The Endemic Equilibrium Point (EE)
Clearly, from the above, we can conclude that a unique positive EE exists whenever R 0 > 1.

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). If R 0 > 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: 4 Abstract and Applied Analysis 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 * : The associated characteristic equation of J * is given by Clearly, a 1 , a 2 , a 3 , and a 1 a 2 − a 3 are positive whenever R 0 > 1. Thus, the Routh-Hurwitz criterion is satisfied. Therefore, the endemic equilibrium point P * of the system (23) is locally asymptotically stable for R 0 > 1. This completes the proof.

Sensitivity Analysis of R 0
In this subsection, we investigate the effect of the various model parameters on the basic reproduction number, R 0 . For this purpose, we used the normalized forward sensitivity index also known as elasticity [24]. If we let R 0 to be a differentiable function with respect to π, then the normalized for-ward sensitivity index of R 0 with respect to the parameter π is given by This normalized sensitivity index measures the relative change of R 0 with respect to π.
To study the sensitivity of R 0 , we choose the parameters β, α, r, ε, and p. 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.
The positive (negative) values indicate a positive (negative) correlation with R 0 , 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 5 Abstract and Applied Analysis total infected population. The effect of the transmission rate ðβÞ has the largest influence on R 0 while the effect of successful treatment rate for the I class ðpÞ has the smallest impact to R 0 . 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.

Numerical Simulation and Discussion
To illustrate some of the analytic results obtained above, we give numerical simulations using the parameter values in Table 1. 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. Figure 2 presents the dynamics of the model for different initial conditions. It shows that only the susceptible population (S * = 8:7 × 10 7 ) 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 R 0 ≤ 1, which supports our analytical results stated in Theorem 3. Similarly, in Figure 3, for R 0 > 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 4 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, R 0 , is varied. In Figure 4(a), the susceptible individuals are being depleted rapidly as R 0 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, R 0 , and all become large as it becomes large.
Generally, Figures 2-4 shows the role of reproduction number to determine the dynamics of TB. It is shown that if R 0 < 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. Figure 5 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.

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, highrisk 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

Abstract and Applied Analysis
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 as R 0 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

Abstract and Applied Analysis
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 p=0.2 p=0.5 p=0.9 Figure 9: Simulation of high-risk exposed and active TB-infected population with different values of p. 9 Abstract and Applied Analysis 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.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that they have no conflicts of interest.