Hasnain Ahmed¹, Md. Farooq Hasan¹, Mohammad Monzu Uddin², Mohammad Saifuddin^1
1Department of Mathematics and Statistics, Bangladesh University of Business and Technology, Dhaka, Bangladesh.
²Department of Public Health, Daffodil International University, Dhaka, Bangladesh.
DOI: 10.4236/ajcm.2025.154020   PDF    HTML   XML   105 Downloads   661 Views   Citations

Abstract

Tuberculosis (TB) is one of the top ten most common infectious diseases in the world. It is also a serious public health problem in Bangladesh, with the number of patients increasing at an alarming rate. In this paper, we are interested in learning more about Bangladesh’s dynamics of tuberculosis, and for this purpose, we proposed and analyzed a four-compartmental model namely the SEIR model. According to the model, the disease-free equilibrium point is asymptotically stable (unstable), and the endemic equilibrium point is unstable (asymptotically stable) if the basic reproduction number is less than (greater than) unity. We analytically calculate the basic reproduction number, and based on the number, we perform the stability analysis of the model at the equilibrium points to understand the epidemic and endemic cases. In our model, the parameters play a significant role in controlling the spread of TB which are generated from TB-reported data from 2016 to 2021 in Bangladesh using the least squares method. The obtained results indicate that the basic reproduction number is greater than unity and this means that TB disease is increasing at an alarming rate every year, and the control strategies should be changed through an alternative system. Moreover, the stability analysis and the sensitivity analysis by numerical and graphical simulation of the model parameters are presented and discussed. The graphical solutions of the model equations were developed using MATLAB as well as computer simulations, and we calculated the model prediction and then compared it to the reported data. Finally, we discuss the accuracy and reliability of the model’s results.

Keywords

Component, Tuberculosis (TB), SEIR Model, Infectious Disease, Mycobacterium, Basic Reproduction Number, Epidemic Model, Stability Analysis

Share and Cite:

1. Introduction

Throughout human history, numerous epidemics, including tuberculosis, influenza, SARS, MERS, and Ebola, have had substantial impacts on many aspects of society, encompassing health, politics, and economics. Consequently, both the scientific community and world agencies have attempted to mitigate the effects of such outbreaks. Since the progress of today’s modern world cannot be separated from mathematics. According to [1], almost all human activities are related to mathematics. The use of mathematical models is a valuable method for addressing real-life problems. It can offer valuable insights into the effectiveness of different TB control measures and is a crucial instrument for investigating the dynamics of tuberculosis (TB) [2]. Several mathematicians, statisticians, and biologists have developed several transmission dynamic models of tuberculosis in the past few decades. The ongoing COVID-19 pandemic serves as a timely reminder of the importance of epidemiological research and the development of mathematical models to comprehensively explain the dynamics governing the behavior of these epidemics. A prominent application involves studying the spread of infectious diseases within specific regions. To understand and predict the spread of such diseases, researchers employ well-known disease spread models, which can be either deterministic or stochastic. The models have their own characteristics based on the type and size of the observed infectious disease spread, for example, the SIR, SEIR, and BSEIR models of tuberculosis [3]. The coronavirus infectious disease has recently brought to our attention the significance of epidemic research and the creation of mathematical models to understand the dynamics of these epidemic diseases [1].

1.1. Mycobacterium Tuberculosis (TB)

Tuberculosis has left its mark on history, with evidence of spinal tubercular decay found in Egyptian mummies dating back to 3000-2400 BCE. However, its recognition as a distinct disease was delayed until the 1820s, largely due to its diverse symptoms. The term tuberculosis was coined in 1834 by Johann Lukas Schönlein. The causative agent, Mycobacterium tuberculosis, was later identified by Nobel Laureate Robert Koch in March 1882, marking a milestone in medical microbiology. It is one of the most contagious diseases and poses a significant global health challenge, with millions of new cases being reported annually. It is caused by the Mycobacterium tuberculosis bacteria and which can damage the kidneys, bones, brain, gastrointestinal tract, and other organs in addition to the lungs. It spreads through the air like the common cold, and the major symptom of tuberculosis is a long-term cough with bloody sputum and fever. The most common symptoms of pulmonary tuberculosis (TB affecting the lungs) include persistent cough, chest pain, fatigue, weight loss, night sweats, fever, and loss of appetite. According to World Health Organization (WHO) [4], the number of people affected with TB each year is declining. Tuberculosis is an infectious disease and a serious public health problem in Bangladesh, with the number of patients increasing at an alarming rate. One person is infected with tuberculosis every minute in the country, and one person dies every 12 minutes. On the occasion of World Tuberculosis Day, this information was given in an event jointly organized by National Tuberculosis Control Program of Bangladesh Health Department and International Centre for Diarrheal Disease Research, Bangladesh (icddr,b).

1.2. Current Status of TB in Bangladesh and Globally

According to the WHO’s 2022 Global TB Report, an estimated 10.6 million people contracted tuberculosis (TB) in 2021, a 4.5% increase from 2020. TB remains the 13th leading cause of death and the second deadliest infectious disease, surpassing HIV/AIDS but trailing COVID-19.

Despite being curable and preventable, TB persists worldwide. In Bangladesh, high population density, poverty, overcrowding and poor living conditions fuel its spread. The report estimates that 3.75 million people fell ill with TB, with around 42,000 deaths, posing a serious threat to the country’s health and economy [4].

In Figure 1, Figure 2(a), and Figure 2(b), it is evident that the number of infected people is increasing in Bangladesh while decreasing globally, except for the year 2021. Additionally, the number of deaths is decreasing in Bangladesh but fluctuating on a global scale.

Figure 1. Bar diagram for TB infected people of Bangladesh and Global from 2016-2021.

(a)

(b)

Figure 2. (a) Bar diagram of TB mortality people of Bangladesh and Global from 2016-2021. (b) Bar diagram of TB mortality people percentage in Bangladesh through globally.

2. Background Study

Mathematical models serve as crucial tools in comprehending and managing TB disease. In light of this, numerous researchers endeavor to construct mathematical frameworks aimed at deciphering the intricate architecture of TB, formulating effective control strategies, innovating novel diagnostic techniques and treatments, and comprehending the disease’s dynamics. Furthermore, these models enable us to anticipate and plan for future disease management and epidemic scenarios. The foundation of utilizing mathematics to grasp disease propagation can be traced back to 1766 when Bernoulli introduced the initial mathematical model for disease spread [5]. Since the 20th century, there has been an escalating interest in employing mathematical models for the study of infectious diseases, with earlier works documented in sources such as [5].

In 1927, Kermack and McKendrick introduced the deterministic Susceptible-Infected-Recovered (SIR) model to describe the dynamics of epidemic propagation [6]. While this foundational model has been widely used to understand disease behavior, its simplicity—omitting additional compartments and control strategies such as vaccination, treatment, quarantine, isolation, and demographic factors like age and sex—limits its realism. Consequently, researchers have developed more sophisticated models that incorporate these complexities to better capture the multifaceted nature of disease spread [7].

The exposed group’s inclusion in the SIR model is an important improvement, representing those infected yet not yet capable of transmitting the disease. This extended construct is known as the Susceptible-Exposed-Infectious-Recovered (SEIR) model [8]. Aron and Schwartz were pioneers in investigating the influence of seasonal variations on epidemic transmission using the SEIR model [8]. In a context involving variable total population size, Li et al. scrutinized the global dynamics of the SEIR model [8]. Newton and Reiter then applied an SEIR model to study the behavior of dengue fever. As studies progressed, the SEIR framework was adopted for analyzing TB dynamics. Chavez and Feng focused on four SEIR models to unravel the intricacies of TB transmission dynamics. Röst and Wu proposed a novel SEIR model that incorporated age-dependent infectivity. Additionally, numerous authors delved into the global stability aspects of TB models.

Epidemic models often use nonlinear differential equations with parameters and state variables. Accurately calibrating these parameters to reflect real-world dynamics is the main challenge, typically tackled via nonlinear optimization. The precise estimation of these parameters is crucial as it allows us to calculate the basic reproduction number, denoted as $R_0$ . This critical values represents the expected number of new infectious cases generated by a single infected individual. The magnitude of $R_0$ holds the key to understanding how rapidly the disease will spread and the potential impact of various control measures. If $R_0>1$ , the disease will lead to epidemics, indicating sustained transmission within the population. Conversely, if $R_0<1$ , the disease will eventually die out as transmission becomes limited, possibly leading to containment and control of the outbreak.

Our goal is to analyze the dynamics of Tuberculosis (TB) transmission under various conditions and identify key factors influencing its spread through mathematical modeling. Specifically, we assess its future impact and the effectiveness of vaccine therapy using the SEIR model. Model parameters were estimated by fitting the original WHO-reported data from 2017 to 2021.

In Section 3, the mathematical model formulation of tuberculosis, the SEIR model. In Section 4, the numerical simulation and the sensitivity analysis of the parameters of the model are investigated, presented, and discussed. In Sections 5 and 6, TB prediction in Bangladesh by using the SEIR model, the solution of the model, and results analysis have been discussed.

3. Mathematical Model Formulation of Tuberculosis

3.1. The SEIR Model

We implement the standard SEIR model to predict the prevalence of tuberculosis, shown in Figure 3. Where S = susceptibles, E = exposeds, I = infectives, and R = recovered populations. The first compartment is the susceptible (S) group, consisting of healthy individuals who have not yet been exposed to the TB bacterium. The second compartment is the exposed (E) group, comprising individuals who have come into contact with the disease but are not yet infectious. The third compartment is the infective (I) group, which includes individuals who have contracted TB and are capable of transmitting the disease. Finally, the fourth compartment is the recovered (R) group, consisting of individuals who have successfully recovered from TB. At any given time t, we denote the proportions of individuals in each compartment of the population as $S\left(t\right),E\left(t\right),I\left(t\right)$ , and $R\left(t\right)$ respectively. The total population, which is homogeneous and isolated, is denoted by $N=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)$ for all $t$ .

Figure 3. Compartmental flow pattern for SEIR model.

Model Assumptions

Our model incorporates the following assumptions:

1) Susceptibles are infected only through direct contact with infectives and standard infection rate is considered.

2) The likelihood of an individual being infected is not influenced by age, sex, societal position, race, or climatic conditions.

3) Equal Interactions: All individuals in the population have an equal probability of interacting with one another.

4) Closed Population: The disease transmission occurs within a confined area without any emigration or immigration.

5) The birth ($b$ ) and death rates ($\mu$ ) are constants and unequal ($b\ne \mu$ ).

3.2. The Model Equations

Undetected or late-detected infections play a significant role in the spread of disease. Therefore, it is crucial to identify such cases promptly and ensure they are isolated in healthcare facilities for immediate treatment and education to prevent further transmission. The time at which the susceptible class changes are equal to the rate at which the infections occur. This occurs when the disease is passed from an infective individual to a susceptible one. The number of susceptible-infective contacts is relative to the product of $S\left(t\right)$ and $I\left(t\right)$ . Hence

$\frac{\text{d}S\left(t\right)}{\text{d}t}=bN-\frac{\beta S\left(t\right)I\left(t\right)}{N}-\mu S\left(t\right)$ .

where $\beta SI$ rate of infection. The term $\frac{\beta SI}{N}$ is negative because the number of susceptible persons decreases.

Let $ϵE$ be the rate at which an exposed individual transitions to an infectious state. Then we have

$\frac{\text{d}E\left(t\right)}{\text{d}t}=\frac{\beta S\left(t\right)I\left(t\right)}{N}-\left(ϵ+\mu \right)E\left(t\right)$

Also consider $\gamma I$ be represent an individuals which are potentially recovers. Therefore, the change in the infected and recovered populations are

$\frac{\text{d}I\left(t\right)}{\text{d}t}=ϵE\left(t\right)-\left(\gamma +\mu \right)\left(t\right)$ and $\frac{\text{d}R\left(t\right)}{\text{d}t}=\gamma I\left(t\right)-\mu R\left(t\right)$ respectively.

Hence, the nonlinear system of differential equations including the tuberculosis patient with unequal death and birth rates has the following form:

$\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=bN-\frac{\beta S\left(t\right)I\left(t\right)}{N}-\mu S\left(t\right)\\ \frac{\text{d}E\left(t\right)}{\text{d}t}=\frac{\beta S\left(t\right)I\left(t\right)}{N}-\left(ϵ+\mu \right)E\left(t\right)\\ \frac{\text{d}I\left(t\right)}{\text{d}t}=ϵE\left(t\right)-\left(\gamma +\mu \right)\left(t\right)\\ \frac{\text{d}R\left(t\right)}{\text{d}t}=\gamma I\left(t\right)-\mu R\left(t\right)\end{array}\}$ (1)

with subject to $S\left(0\right)\ge 0,E\left(0\right)\ge 0,I\left(0\right)\ge 0,R\left(0\right)\ge 0$ .

This model makes the assumption that the birth rate and death rate within the population are unequal. The total population $N\left(t\right)$ , which is homogeneous and isolated for all time t. Moreover, the transmission rate, recovery rate, death rate, and birth rate are denoted by $\beta ,\gamma ,\mu$ and $b$ respectively, and are all non-negative as ($\beta >0$ , $\gamma >0$ , $\mu >0$ , $b>0$ ).

3.3. Non-Dimensionalization form of the SEIR Model

First we consider the following transformatio

$\tilde{S}=\frac{S}{N}$ , $\tilde{I}=\frac{I}{N}$ , $\tilde{R}=\frac{R}{N}$ , and $\tilde{E}=\frac{E}{N}$

where $\tilde{S},\tilde{I},\tilde{R}$ and $\tilde{E}$ represents the fractions of the class susceptible, exposed, and infective and recovery in the population respectively.

$\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=b-\beta S\left(t\right)I\left(t\right)-\mu S\left(t\right)\\ \frac{\text{d}E\left(t\right)}{\text{d}t}=\beta S\left(t\right)I\left(t\right)-\left(ϵ+\mu \right)E\left(t\right)\\ \frac{\text{d}I\left(t\right)}{\text{d}t}=ϵE\left(t\right)-\left(\gamma +\mu \right)I\left(t\right)\\ \frac{\text{d}R\left(t\right)}{\text{d}t}=\gamma I\left(t\right)-\mu R\left(t\right)\end{array}\}$ (2)

We obtain the following non-dimensional form after neglecting tilde (~) sign with subject to $S\left(0\right)\ge 0,E\left(0\right)\ge 0,I\left(0\right)\ge 0,R\left(0\right)\ge 0$ .

3.4. Basic Reproductive Number, $R_0$ for TB Case

In epidemiology field, understanding the transmission dynamics of infectious diseases is crucial for effective disease control and prevention strategies. One of the fundamental concepts used to know epidemic behavior is the basic reproductive number, which is denoted by $R_0$ .

The basic reproductive number is the average number of secondary infections caused by one infectious individual in a fully susceptible population. In other words, it quantifies the contagiousness or transmissibility of a pathogen in a specific setting. By estimating$R_0$ , researchers and public health officials can evaluate the severity of an outbreak, assess the effectiveness of interventions, and predict the future course of an epidemic. Therefore when $R_0<1$ , the infection eventually dies out if each infected individual, on average, produces less than one new infected individual (i.e., $R_0<1$ ). On the other hand, if $R_0>1$ , the infection can spread within the population, as each infected individual generates more than one new case on average.In our model, we employ the next generation matrix to calculate the basic reproductive number, which serves as an important measure for this research.

Once the disease-free equilibrium (DFE) has been established, it becomes necessary to construct a sub-model that focuses solely on the “disease” compartments, which are a subset of the equations in the SEIR model. Consequently, the model (2) transforms into a sub-model that exclusively includes the E and I equations, namely:

$\{\begin{array}{l}\frac{\text{d}E\left(t\right)}{\text{d}t}=\beta S\left(t\right)I\left(t\right)-\left(ϵ+\mu \right)E\left(t\right)\\ \frac{\text{d}I\left(t\right)}{\text{d}t}=ϵE\left(t\right)-\left(\gamma +\mu \right)I\left(t\right) \end{array}$

Next, we need to linearize around the DFE by computing the Jacobian, a matrix of partial derivatives, evaluated at the DFE. We consider:

$F=\beta S\left(t\right)I\left(t\right)-\left(ϵ+\mu \right)E\left(t\right)$ and $G=ϵE\left(t\right)-\left(\gamma +\mu \right)I\left(t\right)$

The sub-model’s Jacobian matrix at the DFE is:

$\begin{array}{c}J\left(S_1^{*},E_1^{*},I_1^{*},R_1^{*}\right)={\left[\begin{array}{cc}\frac{\partial F}{\partial E}& \frac{\partial F}{\partial I}\\ \frac{\partial G}{\partial E}& \frac{\partial G}{\partial I}\end{array}\right]}_{\left(\frac{b}{\mu },0,0,0\right)}={\left[\begin{array}{cc}-\left(ϵ+\mu \right)& \beta S_1^{*}\\ ϵ& -\left(\gamma +\mu \right)\end{array}\right]}_{\left(\frac{b}{\mu },0,0,0\right)}\\ =\left[\begin{array}{cc}-\left(ϵ+\mu \right)& \frac{\beta b}{\mu }\\ ϵ& -\left(\gamma +\mu \right)\end{array}\right]\end{array}$

Finally, we have

$\frac{\text{d}\overrightarrow{x}}{\text{d}t}=J\overrightarrow{X}$, Where $\frac{\text{d}\overrightarrow{x}}{\text{d}t}=\left[\begin{array}{c}\frac{\text{d}E}{\text{d}t}\\ \frac{\text{d}I}{\text{d}t}\end{array}\right]$, $\overrightarrow{X}=\left[\begin{array}{c}E\\ I\end{array}\right]$ , $J=\left[\begin{array}{cc}-\left(ϵ+\mu \right)& \frac{\beta b}{\mu }\\ ϵ& -\left(\gamma +\mu \right)\end{array}\right]$

Then we can write in the form

$\frac{\text{d}\overrightarrow{x}}{\text{d}t}=\left(F-V\right)\overrightarrow{X}$, Where $F=\left[\begin{array}{cc}0& \frac{\beta b}{\mu }\\ 0& 0\end{array}\right]$ , $V=\left[\begin{array}{cc}ϵ+\mu & 0\\ -ϵ& \gamma +\mu \end{array}\right]$

So that $V^{-1}=\frac{1}{\left(ϵ+\mu \right)\left(\gamma +\mu \right)}\left[\begin{array}{cc}\gamma +\mu & 0\\ ϵ& ϵ+\mu \end{array}\right]=\left[\begin{array}{cc}\frac{1}{ϵ+\mu }& 0\\ \frac{ϵ}{\left(ϵ+\mu \right)\left(\gamma +\mu \right)}& \frac{1}{\gamma +\mu }\end{array}\right]$

Therefore, the next generation matrix is given by

$F{V}^{-1}=\left[\begin{array}{cc}0& \frac{\beta b}{\mu }\\ 0& 0\end{array}\right]\left[\begin{array}{cc}\frac{1}{ϵ+\mu }& 0\\ \frac{ϵ}{\left(ϵ+\mu \right)\left(\gamma +\mu \right)}& \frac{1}{\gamma +\mu }\end{array}\right]=\left[\begin{array}{cc}\frac{\beta bϵ}{\mu \left(ϵ+\mu \right)\left(\gamma +\mu \right)}& \frac{\beta b}{\mu \left(\gamma +\mu \right)}\\ 0& 0\end{array}\right]$

The basic reproduction number is $R_0=\rho \left(F{V}^{-1}\right)$ .

Where the spectral radius ($\rho$ ) is the eigenvalue of the given next generation matrix with the largest magnitude. That is:

$R_0=\frac{\beta bϵ}{\mu \left(ϵ+\mu \right)\left(\gamma +\mu \right)}$ (3)

Hence, $R_0$ is called the basic reproduction number and it is the threshold quantity ($R_0$ ) considered to be the result of the contact rate ($\beta$ ).

The number of contacts between susceptible and infective individuals is:

$\sigma =\frac{\beta }{\gamma }$ (4)

3.5. The Herd Immunity Threshold, $H_1$

The herd immunity threshold ($H_1$ ) is the level of immunization needed in a population to stop disease transmission. Therefore the disease spreads or not depends on this threshold. The herd immunity threshold is estimated by the equation:

$H_1=1-\frac{1}{R_0}$ (5)

3.6. The Equilibrium Points of the Model

The fourth equation in system (1) can be excluded as it is independent of the first three variable or equation. To find equilibrium points, the system (2) set equal to zero:

$\begin{array}{lll}b-\beta S\left(t\right)I\left(t\right)-\mu S\left(t\right)=0\hfill & \hfill & (\text{a})\hfill \\ \beta S\left(t\right)I\left(t\right)-\left(ϵ+\mu \right)E\left(t\right)=0\hfill & \hfill & (\text{b})\hfill \\ ϵE\left(t\right)-\left(\gamma +\mu \right)I\left(t\right)=0\hfill & \hfill & (\text{c})\hfill \end{array}\}$ (6)

3.6.1. Disease or Infection Free Equilibrium Point

At infection free equilibrium point, the absence of disease in the system is not specified, so substituting $I=0$ in the system (6), so that, the DFE point is

$D_1=\left(S_1^{*},E_1^{*},I_1^{*}\right)=\left(\frac{b}{\mu },0,0\right)$

3.6.2. Endemic Equilibrium Point

The endemic equilibrium point of (6) is determined when $I\ne 0$ . Then we have the following endemic equilibrium point is

$D_2=\left(S_2^{*},E_2^{*},I_2^{*}\right)=\left(\frac{b}{\mu R_0},\frac{b\left(R_0-1\right)}{R_0\left(ϵ+\mu \right)},\frac{\mu \left(R_0-1\right)}{\beta }\right)$

3.7. Stability Analysis

In this section, we examine the stability of the system (2) in the vicinity of two equilibrium points: the infection-free stability point $\left(S_1^{*},E_1^{*},I_1^{*}\right)$ and the endemic stability point $\left(S_2^{*},E_2^{*},I_2^{*}\right)$ .

After linearization of the system (2) using taylor series expansion about the equlibrium points $\left(S^{*},E^{*},I^{*}\right)$ , we have the following Jacobian matrix:

$J=\left[\begin{array}{ccc}-\beta I^{*}-\mu & 0& -\beta S\\ \beta I^{*}& -\left(ϵ+\mu \right)& \beta S\\ 0& ϵ& -\left(\gamma +\mu \right)\end{array}\right]$ (7)

For system (2) to have stable equilibrium points, all eigenvalues must be negative or have negative real parts otherwise to have unstable points and which is obtained from the Jacobian’s characteristic Equation of (6).

3.7.1. Stability Analysis of DFE Point

At the point $D_1=\left(S_1^{*},E_1^{*},I_1^{*}\right)=\left(\frac{b}{\mu },0,0\right)$ Corresponding Jacobian matrix is

$J=\left[\begin{array}{ccc}-\mu & 0& \frac{-\beta b}{\mu }\\ 0& -\left(ϵ+\mu \right)& \frac{\beta b}{\mu }\\ 0& ϵ& -\left(\gamma +\mu \right)\end{array}\right]$

And the corresponding characteristics equations are:

$\left|J-\lambda I_3\right|=\left|\begin{array}{ccc}-\mu -\lambda & 0& -\beta \\ 0& -\left(ϵ+\mu \right)-\lambda & \beta \\ 0& ϵ& -\left(\gamma +\mu \right)-\lambda \end{array}\right|=0$

$\Rightarrow \left(-\mu -\lambda \right)\left\{\left(\mu +ϵ+\lambda \right)\left(\gamma +\mu +\lambda \right)-\frac{\beta bϵ}{\mu }\right\}=0$

or,

$\lambda =-\mu$ and ${\lambda }^2+\left(\gamma +2\mu +ϵ\right)\lambda +\left(\mu +ϵ\right)\left(\mu +\gamma \right)\left(1-R_0\right)=0$

where $R_0=\frac{\beta bϵ}{\mu \left(\mu +ϵ\right)\left(\mu +\gamma \right)}$

From the quadratic equation formula we have:

$\lambda =\frac{-\left(\gamma +2\mu +ϵ\right)\pm \sqrt{{\left(\gamma +2\mu +ϵ\right)}^2-4\left(\mu +ϵ\right)\left(\mu +\gamma \right)\left(1-R_0\right)}}{2}$

Therefore, the eigenvalues are: ${\lambda }_1=-\mu$

${\lambda }_2=\frac{-\left(\gamma +2\mu +ϵ\right)+\sqrt{{\left(\gamma +2\mu +ϵ\right)}^2-4\left(\mu +ϵ\right)\left(\mu +\gamma \right)\left(1-R_0\right)}}{2}$

${\lambda }_3=\frac{-\left(\gamma +2\mu +ϵ\right)-\sqrt{{\left(\gamma +2\mu +ϵ\right)}^2-4\left(\mu +ϵ\right)\left(\mu +\gamma \right)\left(1-R_0\right)}}{2}$

From these three eigenvalue we see that ${\lambda }_1$ is clearly negative. Another eigenvalue ${\lambda }_2$ and ${\lambda }_3$ also negative when $R_0<1$ .

Hence, $D_1$ is locally asymptotically stable if and only if $R_0<1$ .

3.7.2. Stability Analysis of EEE Point

At the point $D_2=\left(S_2^{*},E_2^{*},I_2^{*}\right)=\left(\frac{b}{\mu R_0},\frac{b\left(R_0-1\right)}{R_0\left(ϵ+\mu \right)},\frac{\mu \left(R_0-1\right)}{\beta }\right)$ Corresponding Jacobian

Matrix is:

$\begin{array}{c}J=\left[\begin{array}{ccc}-\beta \frac{\mu \left(R_0-1\right)}{\beta }-\mu & 0& \frac{-\beta b}{\mu R_0}\\ \beta \frac{\mu \left(R_0-1\right)}{\beta }& -\left(ϵ+\mu \right)& \frac{\beta b}{\mu R_0}\\ 0& ϵ& -\left(\gamma +\mu \right)\end{array}\right]\\ =\left[\begin{array}{ccc}-\mu R_0& 0& \frac{-\beta b}{\mu R_0}\\ \mu \left(R_0-1\right)& -\left(ϵ+\mu \right)& \frac{\beta b}{\mu R_0}\\ 0& ϵ& -\left(\gamma +\mu \right)\end{array}\right]\end{array}$

And the corresponding characteristics equations are:

$\left|J-\lambda I_3\right|=\left|\begin{array}{ccc}-\mu R_0-\lambda & 0& \frac{-\beta }{R_0}\\ \mu \left(R_0-1\right)& -\left(ϵ+\mu \right)-\lambda & \frac{\beta }{R_0}\\ 0& ϵ& -\left(\gamma +\mu \right)-\lambda \end{array}\right|=0$

Or,

${\lambda }^3+p{\lambda }^2+\mu R_0\left(ϵ+\gamma +2\mu \right)\lambda +\mu \left(\mu +ϵ\right)\left(\gamma +\mu \right)\left(R_0-1\right)=0$ (8)

where $p=\left(ϵ+\gamma +\mu \left(2+R_0\right)\right)$

Methods that assess system stability by determining if all roots have negative real parts, without solving the characteristic equation, are crucial. The Routh-Hurwitz criterion offers necessary and sufficient conditions for stability.

3.7.3. Routh-Hurwitz Stability Criterion [9]

Let us consider given the polynomial,

$P\left(\lambda \right)=a_0{\lambda }^n+a_1{\lambda }^{n-1}+a_2{\lambda }^{n-2}+a_3{\lambda }^{n-3}+\cdots +a_{n-1}\lambda +a_n$

where $a_i$ are the real constants, $i=1,2,3,\cdots ,n$ , define the $n$ Hurwitz matrices using the coefficient $a_i$ of the characteristics polynomial:

$H_1=\left[a_1\right],H_2=\left[\begin{array}{cc}{a}_1& 1\\ a_3& a_2\end{array}\right],H_3=\left[\begin{array}{ccc}{a}_1& 1& 0\\ a_3& a_2& a_1\\ a_5& a_4& a_3\end{array}\right]$ and

$H_n=\left[\begin{array}{ccccc}{a}_1& 1& 0& \cdots & 0\\ a_3& a_2& a_1& \cdots & 0\\ a_5& a_4& a_3& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ 0& 0& 0& \cdots & a_n\end{array}\right]$

where $a_i=0$ if $j>n$ .

All the roots of $P\left(\lambda \right)$ are negative or have a negative real part if and only if all Hurwitz determinants are positive: $det{H}_j>0$ , for $j=1,2,3,\cdots ,n$ .

The Routh-Hurwitz criterion for $n=3$ simplifies to:

$det{H}_1=a_1>0$ ; $det{H}_2=a_1{a}_2-a_3>0$ ; $det{H}_3=a_3\left(a_1{a}_2-a_3\right)>0$ .

From the characteristics of Equation (8), we can rewrite it in the following form:

1) $det{H}_1=ϵ+\gamma +\mu \left(2+R_0\right)>0$

2) $\begin{array}{l}det{H}_2=\mu R_0\left\{2R_0\left(ϵ+\gamma +2\mu \right)+3\mu \left(ϵ+\gamma +\mu \right)+{ϵ}^2+{\gamma }^2+\gamma ϵ\right\}\\ +\mu \left(\mu +ϵ\right)\left(\gamma +\mu \right)>0\end{array}$

3) $det{H}_3=a_3\left(a_1{a}_2-a_3\right)>0$ ; if $a_3=\mu \left(\mu +ϵ\right)\left(\gamma +\mu \right)\left(R_0-1\right)>0$

And hence, $a_3>0$ if and only if $R_0>1$ . Hence, EE $D_2$ is locally asymptotically stable if and only if $R_0>1$ , otherwise unstable.

4. Numerical Simulation in Case of Bangladesh

The numerical simulations of the model systems in this part are shown under different parameters. World Health Organization (WHO) and Bangladesh’s National Tuberculosis Control Program are the sources we use to get all the data for these model parameters. TB incidence statistics from the WHO Global Tuberculosis Report for the years 2017 to 2021 were used to modify the model’s parameters. As shown in Table 1, some of the parameters were found in the literature, while others had to be predicted or fitted using data. The table below presents the model’s parameter estimates.

Table 1. Parameters and their estimated values used in the model.

Parameter	Symbol	Value	Source
Birth rate	b	0.0180	WHO
Death rate	μ	0.0055	WHO
Transmission rate of infected population	β	0.6650	Fitted
Recovery rate	γ	0.5840	Fitted
The rate of exposed people become infective	$ϵ$	0.8510	Fitted
Initial number of susceptible	S(0)	139,181,445	Assumed
Initial number of exposed	E(0)	250,000	Assumed
Initial number of infected	I(0)	353,123	Assumed
Initial number of recovered	R(0)	20,000,000	Assumed
Total population in 2016 (N) = 159,784,568

4.1. The SEIR Model

The values of the parameter estimate from Table 1 are substituted into Equation (2), and then we obtain the following system:

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=0.0180-0.665S\left(t\right)I\left(t\right)-0.0055S\left(t\right)\\ \frac{\text{d}E\left(t\right)}{\text{d}t}=0.665S\left(t\right)I\left(t\right)-0.8565E\left(t\right)\\ \frac{\text{d}I\left(t\right)}{\text{d}t}=0.8510E\left(t\right)-0.5895I\left(t\right)\\ \frac{\text{d}R\left(t\right)}{\text{d}t}=0.584I\left(t\right)-0.0055R\left(t\right)\end{array}$ (9)

Here, in the following figure, we show the solution curve of the system (9) for different parameters value and also draw the phase portrait. In Figure 4, we see that the infected population curve passes above the point of intersection of the exposed and recovered population. In this graph, use the parameter value from Table 1. Also, Figure 5 shows the phase portrait of the system (9).

Figure 4. The SEIR model solution curve with the corresponding parameters.

(a)

(b)

(c)

(d)

Figure 5. (a) Phase Portrait for susceptible, infected and exposed population. (b) Phase Portrait for susceptible, infected and recovered population. (c) Phase Portrait for exposed, infected and recovered population. (d) Phase Portrait for exposed, susceptible and recovered population.

4.2. Estimation of the Basic Reproductive Number, $R_0$

From Equation (3), the SEIR model’s basic reproduction number $R_0$ is

$R_0=\frac{0.0180\times 0.665\times 0.851}{0.0055\left(0.5840+0.0055\right)\left(0.851+0.0055\right)}=3.6681>1$

Since $R_0>1$ , the prevalence of tuberculosis indicates the potential for an epidemic. This occurs because the transmission rate exceeds the recovery rate, allowing the infection to spread within the population.

The number of contacts between susceptible and infectious individuals is given by Equation (4), expressed as:

$\sigma =\frac{0.6650}{0.5840}=1.1386$ .

This indicates that, on average, one tuberculosis patient contacts 1.14 susceptible individuals during their infectious period.

4.3. Estimation of the Herd Immunity Threshold, $H_1$

From Equation (5), the herd immunity threshold is given as

$H_1=1-\frac{1}{3.6681}=0.7274$

To effectively control an epidemic during an outbreak, approximately 72.74% of the population must be immunized.

4.4. Estimation of the Equilibrium Points

The equilibrium points from system (2) are the disease-free equilibrium (DFE) and endemic equilibrium (EE) as follows:

$\begin{array}{l}{D}_1=\left(S_1^{*},E_1^{*},I_1^{*},R_1^{*}\right)=\left(3.2727,0,0,0\right),\\ D_2=\left(S_2^{*},E_2^{*},I_2^{*},R_2^{*}\right)=\left(0.8922,0.0153,0.022,2.336\right)\end{array}$

4.5. Stability Analysis of the DFE and EE Points

4.5.1. Stability Analysis at DFE $D_1$

The Jacobian matrices for the system (2) using (7) evaluated at the DFE point is given by,

$J=\left[\begin{array}{ccc}-0.0055& 0& -2.1764\\ 0& -0.8565& 2.1764\\ 0& 0.851& -0.5895\end{array}\right]$

Therefore, the eigenvalues are

${\lambda }_1=-0.0055$ ; ${\lambda }_2=0.6444$ ; ${\lambda }_3=-2.090$

From these three eigenvalue we see that ${\lambda }_1$ and ${\lambda }_3$ clearly negative. Another eigenvalue ${\lambda }_2$ is positive. Also we have $R_0=3.6681>1$ . The disease-free equilibrium is an unstable steady state, meaning that introducing an individual infected with Mycobacterium tuberculosis into a susceptible population will likely cause an outbreak.

4.5.2. Stability Analysis at EE $D_2$

The Jacobian matrices for the system (2) using (7) evaluated at the EE point are given by

$J=\left[\begin{array}{ccc}-0.0202& 0& -0.5933\\ 0.0147& -0.8565& 0.5933\\ 0& 0.851& -0.5895\end{array}\right]$ (10)

The characteristic polynomial equation of the above Jacobian matrix is given by

$\Rightarrow {\lambda }^3+1.4662{\lambda }^2+0.0292\lambda +0.0074=0$

And according to Routh-Hurwitz Stability Criterion, we have

$det{H}_1=1.4662>0$ ; $det{H}_2=0.0354>0$ ; $det{H}_3=0.0003>0$

The characteristic polynomial Equation of (10) satisfies the Routh-Hurwitz Stability Criterion. Therefore, all the eigenvalues are negative and the endemic equilibrium points of the system (2) is asymptotically stable.

4.6. Sensitivity Analysis

The partial derivatives of the basic reproduction number $R_0$ with respect to model parameters $\beta ,\gamma$ and $ϵ$ have been calculated. Here $b$ and $\mu$ constant. So we have the following

1) $\frac{\partial R_0}{\partial \beta }=\frac{bϵ}{\mu \left(ϵ+\mu \right)\left(\gamma +\mu \right)}>0$ , it means that, the value of $R_0$ will increase when the value of $\beta$ also increases, that is there is a proportional relation between them and also shows that the number of infected will increase faster.

2) $\frac{\partial R_0}{\partial \gamma }=-\frac{b\beta ϵ}{\mu \left(ϵ+\mu \right){\left(\gamma +\mu \right)}^2}<0$ , it can be said that, the value of $R_0$ will decrease when the value of parameter $\gamma$ also increases.

3) $\frac{\partial R_0}{\partial ϵ}=\frac{b\beta \mu }{\mu \left(\gamma +\mu \right){\left(ϵ+\mu \right)}^2}>0$ , it can be said that, the value of $R_0$ will increases when the value of parameter $ϵ$ also increases.

4.6.1. Sensitivity Analysis by Mathematically

For DFE point:

1) If $\beta$ is increased, $\gamma$ and $ϵ$ are remain constant or same: Let $\beta =0.958$ , ${\lambda }_1=-0.0055<0$ ; ${\lambda }_2=0.9158>0$ ; ${\lambda }_3=-2.3619<0$ , $R_0=\frac{b\beta ϵ}{\mu \left(ϵ+\mu \right)\left(\gamma +\mu \right)}=5.2844>1$

2) If $\beta$ is reduced, $\gamma$ and $ϵ$ are remain constant or same: Let $\beta =0.175$

${\lambda }_1=-0.0055<0$ ; ${\lambda }_2=-0.0122<0$ ; ${\lambda }_3=-1.434<0$ and $R_0=0.9653<1$

3) If $\beta$ and $ϵ$ is remain same and $\gamma$ is increase: Let $\gamma =2.75$

${\lambda }_1=-0.0055<0$ ; ${\lambda }_2=-0.1462<0$ ; ${\lambda }_3=-3.4658<0$ and $R_0=0.7853<1$

4) If $\beta$ and $ϵ$ is remain same and $\gamma$ is increase: Let, $\gamma =0.384$

${\lambda }_1=-0.0055<0$ ; ${\lambda }_2=0.7578>0$ ; ${\lambda }_3=-2.0032<0$ and $R_0=5.5517>1$

5) If $\beta$ and $\gamma$ is remain same and $ϵ$ is increase: Let $ϵ=1.50$

${\lambda }_1=-0.0055<0$ ; ${\lambda }_2=0.8164>0$ ; ${\lambda }_3=-2.9114<0$ and $R_0=3.6784>1$

6) If $\beta$ and $\gamma$ is remain same and $ϵ$ is decreases: Let $ϵ=0.002$

${\lambda }_1=-0.0055<0$ ; ${\lambda }_2=-0.0001<0$ ; ${\lambda }_3=-0.5969<0$ and $R_0=0.9845<1$ .

That is, the disease-free equilibrium point is stable when $\gamma$ is increased, $\beta$ and $ϵ$ is decreased. Also, unstable when vice-versa.

For Endemic Equilibrium (EE) point:

According to Routh-Hurwitz Stability Criterion, we have

1) If $\beta$ is increased, $\gamma$ and $ϵ$ are remaining constant or same: Let $\beta =0.958$ .

$det{H}_1=1.29>0$ ; $det{H}_2=0.0579>0$ ; $det{H}_3=0.0007>0$ ; $R_0=\frac{b\beta ϵ}{\mu \left(ϵ+\mu \right)\left(\gamma +\mu \right)}=7.998>1$ .

2) If $\beta$ is reduced, $\gamma$ and $ϵ$ are remain constant or same: Let $\beta =0.105$ , we have

$det{H}_1=1.2508>0$ ; $det{H}_2=0.0077>0$ ; $det{H}_3=-0.000002<0$ ; $R_0=0.8766<1$ .

3) If $\beta$ and $ϵ$ is remain same and $\gamma$ is increase: Let $\gamma =2.75$

$det{H}_1=3.6163>0$ ; $det{H}_2=0.0592>0$ ; $det{H}_3=-0.0001<0$ ; $R_0=\mathrm{0.0.7848}<1$ .

4) If $\beta$ and $ϵ$ is remain same and $\gamma$ is decreases: Let $\gamma =0.384$

$det{H}_1=1.2765>0$ ; $det{H}_2=0.0402>0$ ; $det{H}_3=0.0003>0$ ; $R_0=5.5517>1$ .

5) If $\beta$ and $\gamma$ is remain same and $ϵ$ is increase: Let $ϵ=1.50$

$det{H}_1=2.1152>0$ ; $det{H}_2=0.0766>0$ ; $det{H}_3=0.0010>0$ ; $R_0=3.6784>1$ .

6) If $\beta$ and $\gamma$ is remain same and $ϵ$ is decreases: Let $ϵ$ is very small, that is $ϵ=0.001$

$det{H}_1=0.5991>0$ ; $det{H}_2=0.0011>0$ ; $det{H}_3=-1\times {10}^{-8}<0$ ; $R_0=0.5680<1$ .

That is, the endemic equilibrium point is unstable when $\gamma$ is increased, $\beta$ and $ϵ$ is decreased. Also, unstable when vice versa.

4.6.2. Sensitivity Analysis by Simulation

The data were used to numerically simulate the SEIR model for TB. The parameter values are listed in Table 1, and it was produced using MATLAB R2023a. The MATLAB codes are in the appendix. Below are the parameter settings for Susceptible (S), Expose Infected (I), and Recovered (R), along with the effects and modifications that the model would undergo as a result of changing the values of each compartment. Time was measured in months throughout the course of a year, depending on when TB prevalence develops. Some graphing and assumptions were made to show the effects of changes in each compartment of the model.

In Figure 6: There were 500 individuals staying there at first, and they were all deemed vulnerable. This may indicate that the community was free of any exposed, infected, or recovered individuals. The graphs in the subsequent Figure 6 were generated by the simulation. As seen in the graph, there were 500 susceptible people throughout the study period, whereas the number of exposed, infected, and recovered people stayed at zero. This means that there won’t be any influence on any compartment when every member of the population falls into the vulnerable group.

Figure 6. Solution Curve of SEIR model for DFE point corresponding to Table 1 parameter value.

From Figure 7(a), it can be seen that when the value of transmission rate β is changes then there is a greater effect on susceptible compartment people and the population of susceptible decrease throughout time with increasing β and stabilized at the end of the research period. In Figure 7(b) we observed that β have a significant effect on the exposed compartment and the exposed population increases and acquired a pick point, then slowly decrease when β is increases. From Figure 7(c), it can be observed that the value of transmission rate (β) when changed, then there is a momentous impact on the infected compartment people. If the value of β increased, we see that the number of infected population increases and attained a pick point then decrease slowly. According to Figure 7(d), we have seen that the population of recovered increase graphically with the time as the value of β increases.

According to Figure 8(a), it can be seen that the rate of exposed people become infective ($ϵ$ ) then there is less effect on susceptible compartment people and the population of susceptible decreases when $ϵ$ increases. From Figure 8(b), we observed that when $ϵ$ increases then the exposed population is decreases. In Figure 8(c), we found that the value of $ϵ$ when changes, then there is a significant change on the infected compartment people. If the value of $ϵ$ increases we see that the number of infected population increases. According to Figure 8(d), we have seen that the population of recovered increases with the value of $ϵ$ increases.

(a)

(b)

(c)

(d)

Figure 7. (a) Susceptible Population decreases as β increases. (b) Exposed population increase as β increases. (c) Infected Population increases as β increases. (d) Recovered Population increases as β increases.

According to Figure 9(a), it can be seen that when the value of recovery rate increases then there is less effect on susceptible compartment people and the population of susceptible increase throughout time with increasing γ. From Figure 9(b), we observed that when γ is increases then the exposed population decreases. In Figure 9(c), we found that the value of recovery rate (γ) when changed, then there is a significant change on the infected compartment people. If the value of γ increases, we see that the number of infected population decreases. According to Figure 9(d), we have seen that the population of recovered increase significantly with the value of γ increases.

(a)

(b)

(c)

(d)

Figure 8. For the different value of the parameter $ϵ$ respective figures shows (a) Susceptible Population; (b) Exposed Population; (c) Infected Population; (d) Recovered Population.

From the above discussion, we have the following findings:

For DFE point

(a)

(b)

(c)

(d)

Figure 9. (a) Susceptible Population for different value of γ; (b) Exposed Population for different value of γ; (c) Infected Population for different value of γ; (d) Recovered Population for different value of γ.

1) The disease-free equilibrium point is more unstable when β and $ϵ$ increase, γ decreases.

2) The disease-free equilibrium point is more stable when β and $ϵ$ decreases, γ increases.

For EE point

1) The endemic equilibrium point is more stable when β and $ϵ$ increases, γ decreases.

2) The endemic equilibrium point is more unstable when β and $ϵ$ decreases, γ increases.

5. TB Prediction in Bangladesh by Using Model

In this section, Initial data have been taken as the following values in order to present the numerical results for all three models. Since the epidemic data from 2017 to 2021 have been taken into account, now we consider the total initial population as $N\left(0\right)=159784568$ and this initial population is the same as the reported population data of Bangladesh in 2016. Also, we consider the initial infected population $I\left(0\right)=353123$ , which is the same in the TB reported data 2017 was given by World Health Organization (WHO), and the initial exposed and recovered population has been assumed by comparing the initial infected data. The initial susceptible population is calculated by the following instruction:

$S\left(0\right)=N-E\left(0\right)-I\left(0\right)-R\left(0\right)$ .

5.1. Solution of the SEIR Model

The SEIR epidemic model, represented by a system of ODEs, was solved using MATLAB R2023a with the “ode15s” command. This command employs a backward differentiation formula with a quasi-constant step technique and a maximum order of k = 5. Since real incidence data is in years, the numerical solution timelines were also expressed in years.

Figure 10 shows a bar graph comparing annual total infections from actual cases and SEIR model predictions. The results indicate that the SEIR model underestimates actual cases, with the real data showing a linear increase in infections, while the model’s predictions fluctuate and increase annually.

Figure 10. The bar diagram for the SEIR model prediction and reported data.

In Figure 11, we observed in the plotted graph the behaviours of the individual population for the SEIR model. As seen in the plotted graph, the susceptible population increases with time because of the unequal birth and death rates and the birth rate is greater than the death rate. The exposed population fluctuates with time, and after a particular time, the exposed population increases. The infected population depends on the exposed population, so the infected population also fluctuates and after a time increases over time. Hence, the population of recovered is also increased. Therefore, it can be concluded that the TB disease is not under control in Bangladesh according to the SEIR model.

Figure 11. The behaviors of the compartments population over time for the SEIR model.

5.2. The SEIR Model Suggests Some Recommendations for Control of This Epidemic Disease, Which Are as Follows

Figure 12. The total number of infected individuals for different values of β.

Figure 13. The total number of infected individuals for different values of γ.

Figure 14. The total number of infected individuals for different values of $ϵ$ .

Figures 12-14 demonstrate how the model’s dynamics are influenced by the basic reproduction number, $R_0$ , which reflects the average number of new infections caused by one infected individual. As $R_0$ depends on model parameters, sensitivity analysis was used to assess their impact on disease transmission. The figures show how variations in parameters from Table 1 affect disease dynamics, with curves indicating a 20% - 30% increase in specific parameters. It is clear that increasing the recovery rate $\gamma$ or reducing the transmission rate $\beta$ and exposed rate $ϵ$ can reduce the total number of infections. Controlling these parameters could lead to the disease being brought under control after some time.

6. Results Analysis

The article studied highlighted the alarming situation of tuberculosis (TB) in Bangladesh, as depicted by the SEIR model. The model reveals that the disease is currently out of control, with the number of infected individuals, susceptible and exposed population steadily increasing over time. This worrisome trend emphasizes the urgent need for effective measures to curb the spread of TB in Bangladesh.

Based on the obtained parameters, the basic reproductive number $R_0=3.6681$ , this value indicates that the disease will spread in the event of an outbreak and one tuberculosis patient comes into contact with 1.1386 susceptible individuals. The herd immunity threshold was determined to be 0.7274. This indicates that approximately 72.74% of the population needs to be immunized to effectively control the transmission of the disease.

The disease-free equilibrium $\left(S_1^{*},E_1^{*},I_1^{*},R_1^{*}\right)=\left(3.2727,0,0,0\right)$ was identified as an unstable steady state, whereas the endemic equilibrium $\left(S_2^{*},E_2^{*},I_2^{*},R_2^{*}\right)=\left(0.8922,0.0153,0.0220,2.336\right)$ was found to be asymptotically stable.

According to the sensitivity analysis, the results illustrate that an increase in the transmission and infected from exposed rate or a decrease in the recovery rate or leads to $R_0>1$ , resulting in an unstable disease-free equilibrium. This implies that the disease will spread during an outbreak. Conversely, decreasing the transmission and infected from exposed rate or increasing the recovery rate results in $R_0<1$ , making the disease-free equilibrium asymptotically stable. Consequently, the disease will not spread in such scenarios.

Moreover, when the transmission infected from exposed rate is increased or the recovery rate is decreased, $R_0>1$ , indicating the stability of the endemic equilibrium. This suggests that the disease will spread during an outbreak. Conversely, as the transmission and infected from exposed rate decreases or the recovery rate increases, $R_0<1$ , and the endemic equilibrium becomes unstable. This implies that the disease will not spread.

7. Conclusions

Bangladesh, a resource-limited country with a high tuberculosis burden, faces challenges in understanding the disease’s transmission dynamics and epidemiology. This study aims to improve understanding of tuberculosis and other infectious diseases in Bangladesh and their future impact. The SEIR model was used to investigate epidemic spread within the population, with model parameters estimated using the least squares method to fit reported infection data. Stability analysis and the significance of the basic reproduction number were also conducted, showing its critical role in determining the model’s dynamics.

Results indicate that tuberculosis is an epidemic in Bangladesh, and control strategies must be revised. Without effective interventions, TB will pose a serious threat to the health system. Sensitivity analysis of the model parameters was performed, revealing that models without vaccination fail to capture real-world complexities, likely due to the linear nature of the models. Key limitations include the absence of age-structure, immigration, and seasonality factors, and assuming uniform immune responses across the population. Future research should incorporate non-linear terms, improve immigration modeling, and develop sub-models for quarantine measures to better reflect their impacts.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ucakan, Y., Gulen, S. and Koklu, K. (2021) Analysing of Tuberculosis in Turkey through SIR, SEIR and BSEIR Mathematical Models. Mathematical and Computer Modelling of Dynamical Systems, 27, 179-202.[CrossRef]
[2]	Kuddus, M.A., Meehan, M.T., Sayem, M.A. and McBryde, E.S. (2021) Scenario Analysis for Programmatic Tuberculosis Control in Bangladesh: A Mathematical Modelling Study. Scientific Reports, 11, Article No. 4354.[CrossRef] [PubMed]
[3]	Rangkuti, Y.M., Sinaga, M.S., Marpaung, F. and Side, S. (2014) A VSEIR Model for Transmission of Tuberculosis (TB) Disease in North Sumatera, Indonesia. AIP Conference Proceedings, 1635, 201-208.[CrossRef]
[4]	World Health Organization (2010) Tuberculosis, Fact Sheet No. 104. https://www.who.int/en/news-room/fact-sheets/detail/tuberculosis
[5]	Gutierrez, M.C., Brisse, S., Brosch, R., Fabre, M., Omaïs, B., Marmiesse, M., et al. (2005) Ancient Origin and Gene Mosaicism of the Progenitor of Mycobacterium Tuberculosis. PLOS Pathogens, 1, e5.[CrossRef] [PubMed]
[6]	Kermack, W.O. and McKendrick, A.G. (1927) A Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Society of London. Series A, 115, 700-721.
[7]	Capasso, V. and Serio, G. (1978) A Generalization of the Kermack-Mckendrick Deterministic Epidemic Model. Mathematical Biosciences, 42, 43-61.[CrossRef]
[8]	Li, M.Y., Graef, J.R., Wang, L. and Karsai, J. (1999) Global Dynamics of a SEIR Model with Varying Total Population Size. Mathematical Biosciences, 160, 191-213.[CrossRef] [PubMed]
[9]	Gantmacher, F.R. (1964) Matrix Theory. Interscience. Vol. II. Chelsea Pub. Co.