Analysing of Tuberculosis in Turkey through SIR, SEIR and BSEIR Mathematical Models

ABSTRACT Since mathematical models play a key role in investigating the dynamics of infectious diseases, many mathematical models for these diseases are developed. In this paper, it is aimed to determine the dynamics of Tuberculosis (TB) in Turkey, how much it will affect the future and the impact of vaccine therapy on the disease. For this purpose, three mathematical models (SIR, SEIR and BSEIR) in the literature are considered for the case of Turkey. The model parameters are obtained with TB reported data from 2005 to 2015 by using the least square method. The obtained results revealed that the basic reproduction ratio for all three models is less than 1. Moreover, the stability analysis of the models and sensitivity analysis of the model parameters are presented and discussed. Finally, the accuracy of results for all three models is compared and the effect of the vaccination rate is discussed.


Introduction
During human history, many epidemics such as tuberculosis, influenza, SARS, MERS and Ebola have affected the population in many respect such as health, politics and economy. Therefore, scientists and governments have tried to keep epidemics under control. Recently, the current outbreak of the coronavirus disease (COVID-19) re-exposes the importance of epidemic researches and development of the mathematical models to describe the behaviour of epidemics. Among epidemic diseases, Tuberculosis (TB) which is a chronic infectious disease caused by Mycobacterium (MTB) usually affects the lungs but can also other organs like the brain, kidneys, gastrointestinal tract, bone, lymph nodes, etc. The bacteria that cause TB are spread through contaminated air released during the coughing of TB patients. Compared with other diseases caused by a single infectious agent, TB is the second leading lethal disease all over the world, especially in Asia and Africa. According to the World Health Organization (WHO) data, 1.8 million people died of the disease and 10.4 million fell ill in 2015. This shows that TB is a danger to human health and affects economic and social life negatively. Therefore, the current state of the disease should be understood and control programmes should be planned in order to prevent the spread of the disease.
Mathematical models have an important role in the planning of TB control programs. Modelling is a useful tool to understand the dynamics of an epidemic that would help to prevent spreading. Besides, models contribute to predict the future of epidemic and control of the disease. The first study on the mathematical modelling of the spread of disease was proposed by Bernoulli in 1766 [1]. Since the 20th century, interest in studying mathematical modelling of infectious diseases has increased. The earlier related studies can be found in [2][3][4][5][6][7][8][9].
A significant improvement to the SIR model is the addition of the exposed group which is infected but not infectious, called the Susceptible -Exposed -Infected -Recovered (SEIR) model [25][26][27][28][29]. The role of the season on the transmission of an epidemic was first investigated by Aron and Schwartz [30] using the SEIR model. Li et al. [31] studied the global dynamics of the SEIR model in the case of variable total population size. Newton and Reiter [32] developed an SEIR model for observing the behaviour of dengue fever. After that, the studies on TB began to model by the SEIR models. Chavez and Feng [33] focused on four models to understand the disease transmission dynamics of TB. Röst and Wu [34] proposed a new SEIR model in which the infectivity depends on age. Dontwi et al. [35] described the spreading of TB in Amansie West District Ghana by using the standard SEIR model. Yali Yang et al. [36] evaluated the cost of control strategies by using an SEIR model. Side et al. [37] proposed a SIR and an SEIR models for TB and analysed these models. Zhang et al. [38] set up a new mathematical model for TB in China using the data from January 2005 to December 2012. Xu et al. [39] proposed a mathematical model to investigate the control and precautions in Guangdong of China. Besides, many authors interested in global stability of TB models [40][41][42][43].
One of the most important factors in preventing and controlling the spread of tuberculosis is vaccination. Since 1921, the Bacillus Calmette-Guérin (BCG) vaccine remains the most widely used vaccine for the prevention of TB. Evidence has shown that BCG has a protective efficacy of about 75% in preventing some serious types of TB (e.g. meningitis) in children [44]. BCG is currently administered to newborns of high-risk populations as part of the World Health Organization (WHO) Expanded Programme on Immunization (EPI) [45]. Therefore, mathematical modelling of TB including vaccination has gained more importance to make more accurate predictions. In literature, there are many studies on this topic for various epidemics [46][47][48][49][50]. Specifically for TB, Liu et al. [51] proposed a mixed vaccination strategy that is the combination of constant vaccination and pulse vaccination. Egbetade and Ibrahim [52] set up a new mathematical model incorporating treatment, migration and vaccination. Rangkuti et al. [53] explained the spread of TB in North Sumatera Indonesia using VSEIR, which was created by adding the vaccination parameter to the SEIR model. Egonmwan et al. [54] formulated a mathematical model that incorporates vaccination of newborn children and older susceptible individuals into the transmission dynamics of TB in a population.
Most of the epidemic diseases are modelled by a system of nonlinear ordinary differential equations with respect to model parameters and state variables. The main problem in these models is to determine the model parameters that describe the behaviour of the model in a realistic way. Generally, these parameters are adjusted by using nonlinear optimization techniques. Obtaining the parameters is important to calculate the basic reproduction number, R 0 , which represents the expected number of new cases of infectious generated by an infected individual. R 0 is a key to understand how fast the disease will spread and the impact of control strategies. If R 0 > 1, disease breaks out into epidemics, but R 0 < 1, disease dies out.
TB is a lethal disease and is still struggled hard to control for some countries. Even though TB is controlled in Turkey, the geopolitical position of the country and letting in too many immigrants always has put it at risk. Therefore, the development of researches and prevention strategies should be continued. Although the studies [55,56] have been successfully used to describe the spreading of TB in Turkey, in these models, the assumption of absence of some factors such as the birth and the death rate, exposed individuals and prevention, lead to unrealistic estimations. Therefore, in this paper, TB disease in Turkey has been analysed by three epidemiological models (a modified SIR, an SEIR, and a BSEIR) to obtain more realistic predictions. Our aim is to investigate and discuss the characteristics of all three models regarding TB in Turkey, the information they provide and the situations in which the models are used. Besides, one of our main interests is remarking the effect of vaccination on the propagation of virus. In the modified SIR model, the population has been characterized by a death rate and a birth rate equal to the death rate. Through using the SEIR model, it has been included a latent or incubation time which means a certain time for an infected individual to become infectious. The control of the spread of TB through BCG vaccination therapy has been investigated with the BSEIR model. Besides, the model parameters have been determined by fitting the real data reported by the WHO [57][58][59][60][61][62][63][64][65][66][67]. The obtained results concluded that the basic reproduction number for all three models is R 0 < 1. This means that the disease in Turkey is not alarming, but since the number of patients does not decrease dramatically, the control strategies should be continued. The stability analysis of the disease-free equilibrium and endemic equilibrium points is investigated for all three models. Moreover, the sensitivity analysis of the model parameters has been performed and their simulation results have been plotted. In addition, the %95 confidence interval estimate graphs have been presented. The computations show that the predictions produced by all three models approximate the real data very well. Besides, it is emphasized that the model including the vaccination factor (the BSEIR) presents more realistic approaches. Furthermore, based on our literature review, such a comprehensive study is not worked through TB for Turkey so far.
The paper is organized as follows: In Section 2, the models, the modified SIR, the SEIR, and the BSEIR, for describing the dynamics of Tuberculosis are introduced and the equilibria and stability analysis of the models are studied. In addition, the sensitivity analysis of the parameters of the models is investigated. Then, the numerical simulations of all three models are presented and discussed how this is fitted against the data available in Turkey in Section 3. Besides, the predictions of the three models are compared and the effect of vaccination is remarked. Finally, it is concluded the paper in Section 4.

The modified SIR model
Since the basic SIR model proposed by Kermack and McKendrick neglects any natural deaths and births, some authors modified the model by considering a population by a death and a birth rates. Thus, the nonlinear system of differential equations including the death and the birth rates for TB is given by where β, γ, μ, and b are the transmission rate, the recovery rate, the birth rate and the death rate, respectively. In this model, it is assumed that the death and the birth rates are equal to each other. This model comprises three subgroups: SðtÞ represents the number of noninfected but liable to infection in the population at time t, IðtÞ is the number of the infected individuals who can transmit the disease to the susceptible individuals in the population and RðtÞ is the recovered individuals who have been immunized and they are not able to infect again. The total population which is homogeneous and isolated is denoted by N ¼ SðtÞ þ IðtÞ þ RðtÞ for all t. The transmission coefficient β estimates the probability of getting the disease in contact between a susceptible and an infectious individuals. This model does not take into account age, vaccination or waning immunity. Since β and γ are the transition rate, are stated as probabilities, their range is 0 < β; γ < 1.
Equilibrium points and stability. The dimensionless form of the system (1) is presented by where RðtÞ ¼ 1 À SðtÞ À IðtÞ. The third equation in system (1) can be omitted because the first two equations do not depend on it. To find the equilibrium points, the system (2) should be equated to zero: Then the equilibrium points are obtained, namely disease-free equilibrium (DFE) and endemic equilibrium (EE), respectively: The Jacobian matrices for the system (2) are evaluated at the DFE and the EE as follows: In order to the equilibrium points of the system (2) are stable, the eigenvalues must be negative or have a negative real part [68]. The eigenvalues can be obtained by solving the characteristic equation λ 2 À TrðJÞλ þ detðJÞ ¼ 0. Since the eigenvalues must satisfy the rule stated above, detðJÞ > 0 and TrðJÞ < 0 must be proven. and This includes: The DFE D 1 is locally asymptotically stable if and only if β γþμ < 1, otherwise unstable.
The EE D 2 is locally asymptotically stable if and only if β γþμ > 1, otherwise unstable. Here, β γþμ is a threshold value, called the basic reproduction number R 0 , which has an important role in determining dynamics of the disease and is used for prevention strategies. It helps to decide whether the disease spreads or eradicates. If the basic reproduction number is less than 1, then the disease will not spread in the population and dies out. Otherwise, the disease breaks out into an epidemic. Therefore, R 0 is a crucial significance in mathematical and biological aspects.

The SEIR model
The SEIR epidemic model is a generalization of the SIR model. This model on the spread of TB consists of four compartments, namely susceptible ðSðtÞÞ, Exposed ðEðtÞÞ, Infected ðIðtÞÞ and Recovered ðRðtÞÞ, can be interpreted as follows: with subject to Sð0Þ � 0; Eð0Þ � 0; Ið0Þ � 0; and Rð0Þ � 0; where a new compartment EðtÞ which denotes the individuals who are infected but the symptoms of the disease are not yet visible and the parameter P is the rate at which the exposed individuals become infective so that 1 P is the mean latent period. The transmission rate, the recovery rate, the birth and death rates are represented by β, γ, μ and b, respectively. This model assumes that the death and the birth rates are equal to each other. The total population which is homogeneous and isolated is denoted by N ¼ SðtÞ þ EðtÞ þ IðtÞ þ RðtÞ for all t. In this model, the other factors such as age, sex, vaccination, etc. are neglected.
Equilibrium points and stability. The dimensionless form of the system (4) is presented by where RðtÞ ¼ 1 À SðtÞ À IðtÞ À EðtÞ. The fourth equation in system (4) can be omitted because the first three equations do not depend on it. To find the equilibrium points, the system (5) should be equated to zero: Then the equilibrium points are obtained, namely DFE and EE, respectively: Here, Þ is the basic reproduction number [69]. The threshold quantity ðR 0 Þ can be interpreted as the product of the contact rate ðβÞ and average fraction 2 2þμ � � surviving the incubation period and average infectious period 1 γþμ � � [35]. The Jacobian matrices for the system (5) are evaluated at the DFE and the EE as follows: In order to the equilibrium points of the system (5) are stable, the eigenvalues must be negative or have a negative real part [68]. The eigenvalues can be obtained by solving the characteristic equation corresponding to JðD 1 Þ and JðD 2 Þ.
i. Stability analysis of DFE Characteristic equation for D 1 is Then the eigenvalues are ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi λ 1 and λ 2 are negative, clearly. λ 3 < 0 only for R 0 less than 1. Hence, D 1 is locally asymptotically stable if and only if R 0 < 1, otherwise unstable.

ii. Stability analysis of EE
Characteristic equation for D 2 is Methods that allow determining whether all roots have negative real parts and ensuring the stability of the system without solving the characteristic equation are of great importance. The Routh-Hurwitz criterion, which contains the necessary and sufficient conditions for the stability of the system, is one of these methods. Theorem 2.1. (Routh-Hurwitz Stability Criterion) [70] Given the polynomial, where the coefficients a i are the real constants, i ¼ 1; . . . ; n, define the n Hurwitz matrices using the coefficient a i of the characteristics polynomial: The Routh-Hurwitz criteria for which n ¼ 3 taken can be simplified as follows From characteristic equation (9), Hence, a 3 > 0 if and only if R 0 > 1: Hence, it can be concluded that: The EE D 2 is locally asymptotically stable if and only if R 0 > 1, otherwise unstable.

The BSEIR model
In this subsection, to investigate the effect of vaccination on the behaviour of TB, BSEIR model proposed by Liu et al. [51] can be described by the dynamic system of differential equations as follows: where the parameters p ð0 < p < 1Þ and Λ denote the fraction of the newborns vaccinated successfully and the recruitment rate, respectively. Since the positive effect of BCG vaccination is limited, the vaccinated successfully individuals become susceptible again by the rate k. The transmission rate, the recovery rate, the natural death rate and the disease-induced death rate are represented by β, γ, μ and d, respectively. The model comprises five subgroups: BCG vaccinated BðtÞ, Susceptible SðtÞ, Exposed EðtÞ, Infected IðtÞ and Recovered RðtÞ. In addition to the SEIR model, the BSEIR model includes the BCG vaccinated subgroup which denotes the vaccinated newborns successfully. In the BCG protection period, they will not get infected even if they contact infected individuals when the vaccination provides immunity to all of them. The total population, N ¼ BðtÞ þ SðtÞ þ EðtÞ þ IðtÞ þ RðtÞ, changes with respect to time since the natural death rate and the birth rate do not have taken equally. Besides, since the vaccination is thought to prevent for only 10 to 15 years and the natural death of children is about %1, the natural death rate in subgroup BðtÞ is neglected. Equilibrium points and stability. To find the equilibrium points, the system (11) should be equated to zero: Then the equilibrium points are obtained, namely DFE and EE, respectively: where i � 2 ¼ βPΛÀ NμðPþμÞðμþdþγÞ βðPþμÞðμþdþγÞ . We can obtain the basic reproduction number R 0 by using the next-generation matrix [71]: where ρ denotes the spectral radius and the matrices F and V are given by The Jacobian matrix for the system (11) is evaluated at the EE as follows: Solving detðJ À λIÞ ¼ 0, the characteristic equation is Roots of the characteristic equations are λ 1 ¼ À k, λ 2 ¼ À μ and λ 3 ¼ À μ are negative and other two roots satisfies the following quadratic equation Clearly, if a 1 > 0 and a 2 > 0, two roots of (14) will be negative. Hence, D 1 is locally asymptotically stable if and only if R 0 < 1, otherwise unstable.

Numerical solution and estimation of model parameters
The models in the present work have been numerically solved using software MATLAB R2019b by the command 'ode15s'. The reason for choosing this command is the stiffness of differential equation systems that generate the SIR, the SEIR and the BSEIR models. Since stiff problems are characterized by significantly different magnitudes of variation rates, the command 'ode45ʹ may struggle to solve the system. The model parameters have been adjusted based on TB incidence data taken from the WHO Global Tuberculosis Report [57][58][59][60][61][62][63][64][65][66][67] from 2005 to 2015. Some of the parameters have been obtained from the literature [72,73], as seen in Table 1 (15) it has been used the command 'nlinfit' which solves nonlinear regression problems through the Levenberg-Marquardt algorithm in MATLAB R2019b.

Sensitivity analysis
Sensitivity analysis provides to analyse the effect of each parameter on disease transmission and prevalence. It is commonly used to determine the robustness of model predictions to parameter values because of errors in data collection and assumed parameters. It is important to determine parameters that have a high impact on R 0 and should be targeted by intervention strategies. Therefore, here, the partial derivatives of the basic reproduction number R 0 with respect to model parameters p, k, β and γ have been calculated. Since @R 0 @k > 0, increasing k means increasing R 0 . It shows that the number of infected will increase faster. Parameter p is the fraction of vaccinated successfully and @R 0 @p < 0, and parameter γ is the recovery rate and @R 0 @γ < 0, also. Hence, it can be said that the total number of infected populations can be reduced by increasing the parameters p and γ. Besides, being @R 0 @β > 0, it is meant that the infection can be reduced by decreasing the transmission rate β.
To estimate the relative change in a variable when parameters change, sensitivity indices should be calculated. The calculation of these indicates has been executed by means of the following definition.
Definition. The normalized forward sensitivity index of R 0 , which is differentiable with respect to a given parameter σ, is defined by The obtained sensitivity indices of the basic reproduction number R 0 for the baseline model parameters calculated by the formula (16) are presented in Table 3.

Results and discussion
In this section, to illustrate the numerical results for all three models, initial data have been taken as the following values. Since the epidemic data from 2005 to 2015 have been taken into account, the total initial population has been accepted as Nð0Þ ¼ 68;860;540 as the same as the reported population of Turkey in 2005 [72]. The initial infected population was given in the report [57] as Ið0Þ ¼ 20;535. Eð0Þ has been assumed by comparing the rate of the infected to exposed individuals in literature and adapting to data in Turkey. Sð0Þ ¼ N À Eð0Þ À Ið0Þ À Rð0Þ; Sð0Þ ¼ N À Bð0Þ À Eð0Þ À Ið0Þ À Rð0Þ: As indicated in Section 2.2, the parameters β; γ; P; p; k are fitted from real data by using the minimization function (15). The fitted parameter values and initial data have been listed in Table 1.
The epidemic models generated by the systems (1), (4) and (11) have been solved by MATLAB R2019b using command 'ode15s'. This command uses numerical backward differentiation formulas with a maximum order of k ¼ 5 using a quasi constant step method. Since the timelines of real incidence data are according to years, the timelines of obtained numerical solutions have been taken by specifying the interval of integration as a vector of years in the ode15s.
( Figure 1) presents the bar diagram of the per year total infection of real cases and model predictions for all three models. The obtained results revealed that the model predictions approximate the real data very well.
In (Figures (2-4)), it is plotted graphs to show the behaviours of all groups of individuals for all three models. The variation of SðtÞ, IðtÞ and RðtÞ for the modified SIR model is displayed in (Figure 2). As seen in the figure, the population of infected with TB has decreased and hence, the populations of recovered have increased. With this increase, a decrease in the number of susceptible populations is observed. Therefore, it is concluded that the TB disease is under control. Also, the behaviour of the compartments of the SEIR model given by (Figure 3) is similar to the SIR model. Additionally, the exposed population has been decreasing since the infected individuals have been decreasing. However, if the absence of control strategies such as treatment, vaccination, etc., the infected individuals in the population will increase, and therefore, TB becomes an epidemic. The behaviour of the compartments of the BSEIR model has been simulated in ( Figure  4). Unlike the SIR and the SEIR models, the BSEIR model has included the vaccination rate, and the birth and the death rates are not equal to each other in this model. For the prevention of TB, the BCG vaccine is widely used, which is 80% effective in preventing TB and the duration is about 10 years [75]. Therefore, vaccination is one of the treatments for TB patients and hence, an assumption for the transfer of a proportion of the susceptible population to the vaccination class is considered. As can be seen in (Figure 4), the BðtÞ class is increasing. However, since the population is not constant, SðtÞ class shows an increase. The fact that population growth does not have a negative effect on the number of patients due to the vaccination factor. Since the high recovery rate and the effect of the vaccine, the number of infected decreases.
The behaviour of infected individuals for all three models is plotted in ( Figure 5). As displayed in the figure, the reported real incidence data and the obtained predictions for all three models are in good agreement. From the figure, it can be seen the superiority of the BSEIR model when compared to the SIR and the SEIR models. The factors that make the BSEIR model more realistic are the active TB vaccine is included in the model and the population is not constant. Since BCG vaccination is used in Turkey for a long time, when taking into account the vaccination rate has been revealed more realistic  hðtÞ ¼ jI c ðtÞ À I d ðtÞj I d ðtÞ ; (17) where I c ðtÞ and the I d ðtÞ are the model prediction and the corresponding data at time t, respectively. The reported infected data, model predictions and relative errors have been listed in (Table 3). In the calculation, it is used the reported infected data from 2005 to 2015 to compare the model predictions. According to values in (Table 3)    for all three models approach the reported infected data 2016-2018 very well. Hence, it can be said that the model solution is reliable to predict the future behaviour of TB. The error percentages between 2005 and 2018 have been calculated as 0.83%, 0.83%, and 0.45% for the SIR, the SEIR, and the BSEIR models, respectively, with the help of the L 2 norm. As indicated above, the obtained results show that the BSEIR model is more accurate than the SIR and the SEIR model. On the other hand, the epidemic models which are discussed and our methodology have some limitations. The main limitations are to not assuming age-structure, emigration, seasonality, all individuals in the population have the same immune response and the same reaction to the disease for all three models. In addition to these drawbacks, the SIR model does not take into account the special features introduced by the presence of a large set of exposed individuals and, therefore, in the case of too much-exposed individuals, this makes poor predictions. Looking at the limitations of our methodology, the main limitation for all three models is to have limited data for predicting. We used the reported infected data from 2005 to 2015 since data for 2004 and before were not properly listed according to years. Besides, some initial conditions were assumed according to the examples in the literature. This uncertainty has affected the prediction results directly. Therefore, to observe the errors of the predictions caused by the assumed parameters, the sensitivity analysis has been performed. In addition to these disadvantages, fitting the parameters to the reported data has required a bit of computational work.
The dynamical behaviour of the model is defined by the basic reproduction number, R 0 , which is the average number of new infectious caused by a single infective. Since R 0 depends on the model parameters, to investigate the effect of the parameters on disease transmission, the parameters p, k, β and γ were subjected to sensitivity analysis. (Figures (6-8)) show the effect of parameter variation on the disease for the SIR, the SEIR and the BSEIR models, respectively. Each figure displays the number of infected individuals with respect to parameter values in (Table 1) and the corresponding curves with a specific parameter increase of 2%. Clearly from (Figures (6-8) the decrease in the total number of patients is possible by increasing the recovery rate γ or decreasing the transmission parameter β. Considering the BSEIR model in Fig. 8, a certain increase in the number of patients is observed when the rate of successful vaccination p decreases. At the same time, the number of patients decreases when the parameter k, which refers to the rate of loss of effect of the vaccine, decreases. Besides, sensitivity indicates have been given in (Table 2). According to values in (Table 2), the parameter has a positive sign which means that R 0 increases with the parameter. While the parameter has a negative sign means that R 0 decreases for higher values of the parameter. Further, this implies that decreasing (or increasing) of the parameter β by 10% will decrease (or increase) the basic reproductive number by 10%. Also, increasing parameter γ by 10% will decrease the basic reproduction number by 9.946%. Finally, a statistical evaluation of the confidence of the predictions of the SIR, the SEIR, and the BSEIR models given by the %95 confidence interval estimate has been plotted in (Figure 9).     Figure 9. Confidence intervals for the SIR, the SEIR and the BSEIR models.

Conclusion
In this research, the transmission dynamics of TB in Turkey have been analysed and discussed. For this purpose, three epidemic models (the SIR, the SEIR, and the BSEIR) have investigated by describing the spread of the epidemic in a certain population. By presenting the modified SIR model in which the population is fixed and the death and birth rates are equal has been started. After that, the SEIR model which includes the exposed group has been discussed. Finally, a more comprehensive model included the vaccination rate, the BSEIR, has been investigated and compared with the other models.
The results showed that the obtained predictions by the BSEIR model more accurate than the other models. Since BCG vaccination is used in Turkey for a long time, the BSEIR model reveals more realistic predictions for the dynamics of the disease. The model parameters have been determined through the least square method by fitting the reported infected data. Besides, the mathematical analysis for all three models in terms of stability analysis and the importance of the reproduction number R 0 have been discussed. As indicated before, the basic reproduction number has a key role in classifying the dynamical behaviour of the models. According to calculations, R 0 < 1 has been calculated for all three models. It means that the status of TB in Turkey does not lead to an epidemic. However, if necessary precautions do not take, the infected people may increase, because the disease has not eradicated yet. Besides, the sensitivity analysis of the model parameters is performed and discussed in detail. Finally, a statistical interpretation has been added by the confidence interval estimate. It is believed that this study helps to describe the course of TB disease and determine the precautions strategy.