Mathematical Analysis of the Transmission Dynamics of Tuberculosis

Jannatun Nayeem^*, Israt Sultana

Department of Arts and Sciences, Ahsanullah University of Science and Technology, Dhaka, Bangladesh.

DOI: 10.4236/ajcm.2019.93012                

Abstract

We develop a dynamical model to understand the underlying dynamics of TUBERCULOSIS infection at population level. The model, which integrates the treatment of individuals, the infections of latent and recovery individuals, is rigorously analyzed to acquire insight into its dynamical features. The phenomenon resulted due to the exogenous infection of TUBERCULOSIS disease. The mathematical analysis reveals that the model exhibits a backward bifurcation when TB treatment remains of infected class. It is shown that, in the absence of treatment, the model has a disease-free equilibrium (DEF) which is globally asymptotically stable (GAS) and the associated reproduction threshold is less than unity. Further, the model has a unique endemic equilibrium (EEP), for a special case, whenever the associated reproduction threshold quantity exceeds unity. For a special case, the EEP is GAS using the central manifold theorem of Castillo-Chavez.

Keywords

TUBERCULOSIS Model, SIR Model, Equilibria, Stability, Castillo-Chavez Theorem, Disease Dynamics, Disease Endemic Equilibrium

1. Introduction

A differential equation which describes some physical process is often called a mathematical model of the process. Again a differential equation is a mathematical equation that rates some functions of one or more variables with its derivatives, differential equation arises whenever a deterministic relation involving some continuously varying quantities and their rates of change in space and time. These equations occupy the place at center stage of both pure and applied mathematics. For the mathematicians, mathematical modeling offers an important tool in the study of the evolution of diseases such as Tuberculosis, HIV, Hepatitis, Ebola, etc. Various epidemiology models such as SIR, SEIR, SIRS, SEIS, MSEIR, etc. can be built to analyze these types of diseases. Among them the SIR model is widely used in epidemiology and public health to compute number of individuals in each category of the population and to explain the change in the number of people needing medical attention during an epidemic as well as evaluate policies effectively during the endemic Tuberculosis [1]. The susceptible infected removal (SIR) is a system of ordinary differential equation in three dimensions. SIR model is analyzed by building a mathematical theorem which guarantees the existence of a case of Tuberculosis, the disease free equilibrium phase and stage of disease endemic Tuberculosis [2].

Tuberculosis is an infectious disease usually caused by the Mycobacterium tuberculosis (MTB). Tuberculosis is spread through air, just like a cold or flu when people have active TB in their lungs, they are suffering from cough, spit, speak or sneeze. Tuberculosis generally affects the lungs but it can also be other parts of the body like brain and spine. Tuberculosis is contagious, but it is not easy to catch. It has slow intrinsic dynamics, the incubation period and the infectious period spam long term intervals in the order of years on average. Therefore, a mathematical model is needed to have a better insight on the dynamics of the disease [3] [4]. Mathematical models are a simplified representation of how an infection spreads across a population over time. Most epidemic models are based on dividing the population into a small number of compartments.

Tuberculosis is not only a health problem, but also an economic problem of mankind as outbreaks usually lead to enormous expenditure on healthcare. Over 80% of all TB patients live in 22 countries, mostly in Sub-Saharan Africa and Asia [5]. Owing to the rising TB mortality and infection rates (especially in developing countries), the World Health Organization (WHO) declared TB as a global public health emergency in 1993 [6] [7]. Over the years, a number of global initiatives, spear headed by WHO, were embarked upon the hope of minimizing the burden of TB worldwide. These include the Stop TB Partnership, International Standards of Tuberculosis Care and Patient’s Care and the Global Plan to Stop TB [8]. A notably medical contribution in TB control was the introduction of antibiotics that resulted in significant decrease in mortality (for instance, a 70% reduction in TB related mortality was recorded USA between 1945 to 1955) [9] [10]. At present, global per capita rate of TB is increasing at approximately 1.1% per year and the number of cases at 2.4% per year [11]. According to the 2004 WHO report “Global Tuberculosis Control” [12], there were 8.8 million new cases of TB worldwide in 2002, with close to 2 million TB-related deaths, more than any other infectious diseases. In 2011, the treatment success rate continued to be high at 87% among all new TB cases [13]. On the other hand, treatment interruptions are frequent in active TB cases during the intensive phase and the continuation phase because of a wide range of reasons [14]. It may be recognized that treatment interruptions and the missed TB cases are the key factors to cause the more drug-resistant TB cases and the high TB mortality [15]. In 2012, there was estimation that 450 thousand individuals developed multidrug-resistant TB (MDR-TB) and an estimation of 170 thousand deaths from MDR-TB [13], which is currently a main threat to tuberculosis control programs and community health [16]. Moreover, an estimated 1.3 million died from this disease and 8.6 million people developed TB under treatment. From 2013 to 2015 it has seen that 6.1 million people developed TB where infected rate was 34%. But TB treatment averted 49 million deaths globally between 2000 and 2015 [16]. It has seen that the rate of infected Mycobacterium tuberculosis (MTB) will reduce day by day, influencing on Bacilli Calmette-Guerin (BCG) vaccines which were first used in 1921 medically in USA and some TB treatment therapies. Tuberculosis infection can be transmitted through primary progression after a recent infection, re-activation of a latent infection and re-infection of a previously infected individual [17]. A small proportion of those infected will develop primary disease within several years of their first infection. Those who escape primary disease may eventually re-activate this latent infection decade after an initial transmission event. Lastly, latently infected individual can be re-infected by a process known as exogenous re-infection and develop the disease as a result of this new exposure [18]. However, if a person got the tuberculosis, then it couldn’t only be seen by the main caution of factors which is causing tuberculosis disease, but also the possibility of exogenous re-infection happening. Tuberculosis is a vaccine prevented disease. The current practice in most part of world especially in developed countries is that children aged 12 to 15 months are given a dose of the BCG vaccine. In practice, even vaccinated individual may still be susceptible if vaccination failure occurred or their vaccine induced immunity waned.

Recent years have seen an increasing trend in the representation of mathematical models in publications in the epidemiological literature, from specialist journals of medicine, biology and mathematics to the highest impact generalist journals [19], showing the importance of inter disciplinary.

In this paper, we have formulated the transmission dynamics of Tuberculosis in the presence of treatment and investigated its role in the dynamics of the disease.

2. Formulation of Model

Following the classical assumptions, we formulate a deterministic, compact mental, mathematical model to describe the transmission dynamics of measles. The population is homogeneously mixing and reflecting the demography of a typical developing country, as it experiments an exponentially increasing dynamics. In other to describe the model equations, the total population (N) is divided into three classes: Susceptible (X), infected (Y) and Recovered (Z). Here we shall detail the transitions among these four classes as depicted in Figure 1.

The class X of susceptible is increased by birth or immigration at a rate $\Lambda$. It is decreased by infection following contact with infected individuals at a rate β, and diminished by natural death at a rate μ. The class Y is decreased by testing and therapy at a rate r, breakthrough into infected class at a rate β and diminished by natural death at a rate μ. The class Y of infected individuals is generated by breakthrough of individuals at a rate k. The class is decreased by recovery from infection at a rate β and diminished by natural death at a rate μ. The model assumes that both recovered exposed individuals and recovered infected individuals become permanently immune to the disease. This generates a class R of individuals who have complete protection against the disease.

The transitions between model classes can now be expressed by the following system of first order differential equations:

The description of Variables of the TB Model is shown in Table 1.

$\frac{\text{d}X}{\text{d}t}=\Lambda -\beta \left(Y+\eta Z\right)X-\mu X+r_{1}Y+r_{2}Z$ (1)

$\frac{\text{d}Y}{\text{d}t}=\beta \left(Y+\eta Z\right)X-k_{1}Y$ (2)

$\frac{\text{d}Z}{\text{d}t}=kY-k_{2}Z$ (3)

Since the model monitors human population, all the associated parameters and state variables are non negative t ≥ 0. It is easy to show that the state variables of the model remain non-negative for all non-negative initial conditions. Consider the biological feasible region.

$\Omega =\left\{\left(S,I,R\right)\in R_{+}^3:N\to \frac{\Lambda }{\mu }\right\}$

From the model Equation (1) to (3) it will be shown that the region is positive. The total population of individuals is given by

Table 1. Description of variables of the TB model.

$N=S+I+R$

3. Analysis of Model

Disease Free Equilibrium (DFE): The equilibrium points of the system can be obtained by equating the rate of changes of zero, given by ${\epsilon }_0$ ,

$\therefore {\epsilon }_0=\left\{X=\frac{\Lambda }{\mu },Y=0,Z=0\right\}$ (4)

The stability of the DFE will be analyzed using the next generation method [22]. The non-negative matrix F (of the new infection terms) and the non-singular M-matrix V (of the remaining transfer terms) are given, respectively by,

$F=\left(\begin{array}{cc}\frac{\beta \Lambda }{\mu }& \frac{\beta \eta \Lambda }{\mu }\\ 0& 0\end{array}\right)$

and $V=\left(\begin{array}{cc}{k}_1& 0\\ -k& k_2\end{array}\right)$

where, $k_1=r_1+\mu +k,k_2=r_2+\mu +d$.

The associated reproduction number, denoted by $R_0$ , is given by $R_0=\rho \left(F{V}^{-1}\right)$ , where $\rho$ denotes the spectral radius (dominant eigenvalue in magnitude) of the next generation matrix $F{V}^{\prime }$. It follows that

$V^{-1}=\left(\begin{array}{cc}\frac{1}{k_1}& 0\\ \frac{k}{k_1{k}_2}& \frac{1}{k_2}\end{array}\right)$

$\therefore F{V}^{-1}=\left(\begin{array}{c}\frac{\beta \Lambda \left(k_2+\eta k\right)}{\mu k_1{k}_2}\\ 0\end{array}\right)$

$\therefore R_0=\frac{\beta \Lambda \left(k_2+\eta k\right)}{\mu k_1{k}_2}$ (5)

Lemma: The disease free equilibrium ${\epsilon }_0$ of the model (1), (2) and (3), is locally asymptotically stable (LAS) if $R_0<1$ and unstable if $R_0>1$.

The threshold quantity, $R_0$ , is the reproduction number for the model. The epidemiological implication of Lemma 1 is that Tuberculosis spread can be effectively controlled in the community (when $R_0<1$ ) if the initial sizes of the populations of the model are in the basin of attraction of the disease free equilibrium ${\epsilon }_0$.

Since we have considered Tuberculosis model in some stages, are shown the backward bifurcation, where the stable DFE co-exists with a stable endemic equilibrium when the associated reproduction threshold ( $R_0$ ) is less than unity, it is instructive to determine whether or not the model also exhibits this dynamical feature. This is investigated below.

Theorem 1: The model (1), (2) and (3) undergoes a backward bifurcation at $R_0=1$ if the inequality holds.

Proof: The proof of theorem, which is based on the use of center manifold theory. The backward bifurcation phenomenon of the model is numerically illustrated below. It is convenient to let $X=x_1,Y=x_2,Z=x_3$ , so that $N=x_1+x_2+x_3$. Further, by introducing the vector notation $X={\left(x_1,x_2,x_3\right)}^{\text{T}}$ , the model can be written in the form,

$\frac{\text{d}x}{\text{d}t}=F\left(x\right)$ ,

where $F={\left(f_1,f_2,f_3\right)}^{\text{T}}$ , as follows

$\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=f_1=\Lambda -\beta \left(x_2+\eta x_3\right)x_1-\mu x_1+r_1{x}_2+r_2{x}_3\\ \frac{\text{d}{x}_2}{\text{d}t}=f_2=\beta \left(x_2+\eta x_3\right)x_1-k_1{x}_2\\ \frac{\text{d}{x}_3}{\text{d}t}=f_3=k{x}_2-k_2{x}_3\end{array}$ (6)

where, $\lambda =\beta \left(Y+\eta Z\right)$. The jacobian of the system at the DFE $\left({\epsilon }_0\right)$ is given by,

$J\left({\epsilon }_0\right)=\left[\begin{array}{ccc}-\mu & -\frac{\beta \Lambda }{\mu }& -\frac{\beta \eta \Lambda }{\mu }\\ 0& \frac{\beta \Lambda }{\mu }-k-\mu & \frac{\beta \eta \Lambda }{\mu }\\ 0& k& -d-\mu \end{array}\right]$

To analyze the dynamics of the model and we compute the eigenvalues of the jacobian of the equations at the disease free equilibrium (DEF). It can be shown that this jacobian has a left eigenvector is given by $V={\left(v_1,v_2,v_3\right)}^{\text{T}}$ where,

$v_1=0$

$v_2=\text{free}$

and $v_3=\frac{\beta \eta \Lambda }{\mu k_2}{v}_2$

The right eigenvector is given by, $W={\left(w_1,w_2,w_3\right)}^{\text{T}}$

$w_1=-\frac{\beta \Lambda \left(k_2+\eta k\right)}{{\mu }^2{k}_2}{w}_2$

$w_2=\text{free}$ and $w_3=\frac{k{w}_2}{k_2}$

Theorem 2: (Castillo-Chavez and Song)

Consider the following general system of ordinary differential equations with a parameter $\phi$.

$\frac{\text{d}x}{\text{d}t}=f\left(x,\phi \right)$ , $f:{\Re }^n\times \Re \to \Re$ and $f\in C\left({\Re }^n\times \Re \right)$

where 0 is an equilibrium of the system (i.e. $f\left(0,\phi \right)=0$ for all $\phi$ and assume

A1: $A=D_{x}f\left(0,0\right)=\left(\frac{\partial f_i}{\partial x_j}\left(0,0\right)\right)$ is the linearization matrix of the system (6) around the equilibrium 0 with $\phi$ evaluated at 0. Zero is a simple eigenvalue of A and other eigenvalues of A have negative real parts.

A2: Matrix A has a right eigenvector w and a left eigenvector v (each corresponding to the zero eigenvalue).

Let, $f_k$ be the kth component of f and

$\begin{array}{l}a={\displaystyle \underset{k,i,j=1}{\overset{n}{\sum }}{v}_k{w}_i{w}_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}\left(0,0\right)}\\ b={\displaystyle \underset{k,i=1}{\overset{n}{\sum }}{v}_k{w}_i\frac{{\partial }^2{f}_k}{\partial x_i\partial \phi }\left(0,0\right)}\end{array}$

The local dynamics of the system around the equilibrium point 0 is totally determined by the sings of a and b. Particularly, $a<0,b>0$ , the system does not show backward bifurcation at $R_0=1$. In these cases, $0<\phi \ll 1$ , 0 becomes unstable and there exists a positive locally asymptotically stable equilibrium.

Computations of a and b:

$a={\displaystyle \underset{k,i,j=1}{\overset{n}{\sum }}{v}_k{w}_i{w}_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}\left(0,0\right)}=\beta v_2{w}_1\left(w_2+\eta w_3\right)<0$ (7)

$b={\displaystyle \underset{k,i=1}{\overset{n}{\sum }}{v}_k{w}_i\frac{{\partial }^2{f}_k}{\partial x_i\partial \phi }\left(0,0\right)}=\frac{\pi v_2\left(\eta w_3+w_2\right)}{\mu }>0$ (8)

This result is summarized below.

Theorem 3: The model (1), (2) and (3) exhibit backward bifurcation at $R_0=1$ whenever $a<0,b>0$. It should be noted that, in the absence of recovery exposed stage and infected stage the backward bifurcation co-efficient, a is given in below,

$a=\beta v_2{w}_1\left(\eta w_3+w_2\right)<0$

since all the model parameters and the eigenvectors $w_i\left(i=2,3,\cdots \right)$ and $v_i\left(i=1,2,\cdots \right)$ are non-negative and $0<\epsilon <1$. Thus, since the inequality does not hold in this case, the model (1), (2) and (3) will not undergo backward bifurcation in the absence of recovery exposed stage and infected stage. This result is summarized below.

Lemma: The model (1), (2) and (3) does not undergo backward bifurcation at $R_0=1$ in the absence of treatment ( $r_1=r_2=0$ ). If we consider $r_1$ and $r_2$ exist then the coefficient of a maybe positive and b is also positive.

The backward bifurcation phenomenon of the model is numerically illustrated in below:

Simulations of the model shows that Figure 1 is backward bifurcation diagram for the population of susceptible individuals, Figure 2 is backward bifurcation diagram for the population of latent indiviuals and Figure 3 is backward bifurcation diagram for the population of infected indiviuals. Parameter values used are as given in Table 2, with $\pi =2000$ [20], $\mu =0.02$ [20], $d=0.1$ [20], $\eta =0.08$ , $r_1=0.85$ , $r_2=0.9$ , $k=0.7$ [21].

4. Global Stability of DFE of the TB Model

Let,

Figure 1. Backward bifurcation diagram of susceptible class.

Figure 2. Backward bifurcation diagram of latent class.

Figure 3. Backward bifurcation diagram for population of Infected class.

Table 2. Data summary of Parameters of the TB Model.

$\begin{array}{l}\frac{\text{d}X}{\text{d}t}=H\left(X,Z\right)\\ \frac{\text{d}Z}{\text{d}t}=G\left(X,Z\right),G\left(X,0\right)=0.\end{array}$ (9)

where,

$X=\left(X,0\right)$ and $Z=\left(Y,Z\right)$ with the components of $X\in R^1$ denoting the uninfected population and the components of $Z\in R^2$ denoting the infected population.

The disease free equilibrium is now denoted as,

$E_0=\left(X^{*},0,0\right),X^{*}=\left\{\frac{\Lambda }{\mu },0,0\right\}$

Now, $\frac{\text{d}X}{\text{d}t}=H\left(X,0\right)$ , $X^{*}$ is globally asymptotically stable (GAS)

$G\left(X,Z\right)=Pz-\hat{G}\left(X,Z\right),\hat{G}\left(X,Z\right)\ge 0$ for $\left(X,Z\right)\in \Omega$. (10)

where, $P=D_{Z}G\left(X^{*},0\right)$ is an M-matrix (the off diagonal elements of P are non negative) and $\Omega$ is the region where the model makes biological sense. If the system (9) satisfies the conditions of (10) then the theorem below holds.

Theorem 4: The fixed point ${\epsilon }_0\left(X^{*},0\right)$ is a globally asymptotically stable equilibrium of system (9) provided that $R_0<1$ and the assumptions in (10) are satisfied.

Proof:

From the system (1) and (2) we have,

$H\left(X,0\right)=\left(\begin{array}{c}\Lambda -\mu X+r_{1}Y+r_{2}Z\\ 0\end{array}\right)$

$\begin{array}{l}G\left(X,Z\right)=P\left(Z\right)-\hat{G}\left(X,Z\right)\\ \Rightarrow \hat{G}\left(X,Z\right)=P\left(Z\right)-G\left(X,Z\right)\end{array}$ (11)

where,

$P\left(Z\right)=\left(\begin{array}{cc}-k_1& 0\\ k& -k_2\end{array}\right)$

and $G\left(X,Z\right)=\left(\begin{array}{c}\beta YX+\beta \eta ZX-k_{1}Y\\ kY-k_{2}Z\end{array}\right)$

Putting values $P\left(Z\right)$ and $G\left(X,Z\right)$ in (11) no equation and we obtain,

$\hat{G}\left(X,Z\right)=\left(\begin{array}{c}{\hat{G}}_1\left(X,Z\right)\\ {\hat{G}}_2\left(X,Z\right)\end{array}\right)=0$ (12)

It is clear that $\hat{G}\left(X,Z\right)=0$ for all $\left(X,Z\right)\in \Omega$ we also note that matrix P is an M-matrix since its off diagonal elements are non-negative.

5. Endemic Equilibrium Point (EEP) of the TB Model

Let, ${\epsilon }_1=\left(X^{**},Y^{**},Z^{**}\right)$ represents any arbitrary endemic equilibrium of the Tuberculosis model. Solving the Equations (1)-(3), the model has the following endemic equilibrium points (EEP),

$\begin{array}{l}{X}^{**}=\frac{\Lambda }{\mu R_0}\\ Y^{**}=\frac{\Lambda k_2\left(R_0-1\right)}{R_0\left\{\mu r_2+\left(\mu +d\right)\left(\mu +k\right)\right\}}\\ Z^{**}=\frac{\Lambda k\left(R_0-1\right)}{R_0\left\{\mu r_2+\left(\mu +d\right)\left(\mu +k\right)\right\}}\end{array}$

Existence of Endemic Equilibrium Point (EEP): special case

In this section, the possible existence and stability of endemic (positive) equilibria of the model (1), (2) and (3) (i.e. equilibria where at least one of the infected of the model is non-zero) will be consider for the special case where recovery rate from exposed stage and infected stage does not occur (i.e. $r_1=r_2=0$ ).

Let, ${\epsilon }_1=\left(X^{**},Y^{**},Z^{**}\right)$ represents any arbitrary endemic equilibrium of the model (1), (2) and (3) with $r_1=r_2=0$. Solving the equations of the system at steady-stage gives,

$\begin{array}{l}{X}^{**}=-\left\{\beta \left(y+\eta z\right)\right\}-\mu \\ Y^{**}=\beta \left(y+\eta z\right)\\ Z^{**}=0\end{array}$ (13)

The expression for $\lambda$ , defined in (1), (2) and (3) at the endemic steady-state is given by

${\lambda }^{**}=\beta \left(y+\eta z\right)$ (14)

For mathematical convenience, the expression in (14) is re-written,

${\lambda }^{**}=\beta \left(\frac{\Lambda k_2\left(R_0-1\right)}{R_0\left\{\mu r_2+\left(\mu +d\right)\left(\mu +k\right)\right\}}+\frac{\eta k\Lambda \left(R_0-1\right)}{R_0\left\{\mu r_2+\left(\mu +d\right)\left(\mu +k\right)\right\}}\right)$

And we get,

${\lambda }^{**}=\frac{\beta \Lambda \left(R_0-1\right)\left(\eta k+k_2\right)}{R_0\left\{\mu r_2+\left(\mu +d\right)\left(\mu +k\right)\right\}}>0$ (15)

The components of the unique endemic equilibrium ${\epsilon }_1$ can be obtained by substituting the unique value of ${\lambda }^{**}$ , given into the expression in (14). Thus, the following has been established.

Lemma: The model with recovery rate from exposed stage and infected stage $r_1=r_2=0$ has a unique endemic equilibrium, given by ${\epsilon }_1$ , whenever $R_0>1,{\lambda }^{**}>0$.

6. Global Stability of EEP by Non-Linear Lynapunov Function

Theorem 5: The unique EEP

$\begin{array}{c}{\epsilon }_1=\left\{X^{**},Y^{**},Z^{**}\right\}\\ =\left\{\frac{\Lambda }{\mu R_0},\frac{\Lambda k_2\left(R_0-1\right)}{R_0\left\{\mu r_2+\left(\mu +d\right)\left(\mu +k_1\right)\right\}},\frac{\Lambda k\left(R_0-1\right)}{R_0\left\{\mu r_2+\left(\mu +d\right)\left(\mu +k_1\right)\right\}}\right\}\end{array}$ of the model with $r_1=r_2=0$ is globally asymptotically stable (GAS), whenever $R_0>1$.

Proof:

Let, $r_1=r_2=0$ and $R_0>1$ , so the EEP, ${\epsilon }_1$ exists. Consider the following non-linear Lyapunov function,

$L=\left(X-X^{**}-X^{**}\ln \frac{X}{X^{**}}\right)+\left(Y-Y^{**}-Y^{**}\ln \frac{Y}{Y^{**}}\right)+a_1\left(Z-Z^{**}-Z^{**}\ln \frac{Z}{Z^{**}}\right)$ (16)

With Lyapunov derivative is given by,

$\stackrel{\dot{}}{L}=\left(1-\frac{X^{**}}{X}\right)\stackrel{\dot{}}{X}+\left(1-\frac{Y^{**}}{Y}\right)\stackrel{\dot{}}{Y}+\frac{\beta \eta X^{**}}{k_2}\left(1-\frac{Z^{**}}{Z}\right)\stackrel{\dot{}}{Z}$

Now,

$\begin{array}{c}\stackrel{\dot{}}{L}=\beta Y^{**}{X}^{**}+\beta \eta Z^{**}{X}^{**}+2\mu X^{**}-\mu X-\frac{\beta Y^{**}{\left(X^{**}\right)}^2}{X}\\ -\frac{\beta \eta Z^{**}{\left(X^{**}\right)}^2}{X}-\frac{\mu {\left(X^{**}\right)}^2}{X}+\beta Y{X}^{**}-\mu Y-kY\\ -\beta X{Y}^{**}-\frac{\beta \eta Y^{**}Z}{Y}+\beta Y^{**}{X}^{**}+\beta \eta X^{**}{Z}^{**}\end{array}$ (17)

And

$\begin{array}{l}\frac{\beta \eta X^{**}}{k_2}\left(1-\frac{Z^{**}}{Z}\right)\left(kY-k_{2}Z\right)\\ =\beta \eta X^{**}\frac{Y}{Y^{**}}{Z}^{**}-\beta \eta X^{**}Z-\beta \eta X^{**}{Z}^{**}\frac{Z^{**}}{Z}\frac{Y}{Y^{**}}+\beta \eta X^{**}{Z}^{**}\end{array}$ (18)

Adding (17) and (18) and we get,

$\begin{array}{l}=\mu X^{**}\left[2-\frac{X}{X^{**}}-\frac{X^{**}}{X}\right]+\beta X^{**}{Y}^{**}\left[2-\frac{X}{X^{**}}-\frac{X^{**}}{X}\right]\\ +\beta \eta X^{**}{Z}^{**}\left[3-\frac{X^{**}}{X}-\frac{Y}{Y^{**}}\frac{Z^{**}}{Z}-\frac{X}{X^{**}}\frac{Y^{**}}{Y}\frac{Z}{Z^{**}}\right]\end{array}$

Since the arithmetic mean exceeds the geometric mean, it follows then that

$\begin{array}{l}2-\frac{X}{X^{**}}-\frac{X^{**}}{X}\le 0\\ 3-\frac{X^{**}}{X}-\frac{Y}{Y^{**}}\frac{Z^{**}}{Z}-\frac{X}{X^{**}}\frac{Y^{**}}{Y}\frac{Z}{Z^{**}}\le 0\end{array}$ (19)

So that $\stackrel{\dot{}}{L}\le 0$ for $R_0>1$. Hence, L is a Lyapunov function of the system with $r_1=r_2=0$ on $\Omega$. In other words, $\underset{t\to \infty }{\lim }\left(X,Y,Z\right)=\left(X^{**},Y^{**},Z^{**}\right)$.

Thus, by the Lyapunov function L and LaSalle’s Invariance Principle every solution to the equation in the model, with $r_1=r_2=0$ approaches ${\epsilon }_1$ as $t\to \infty$ for $R_0>1$.

7. Numerical Simulation and Discussions

The effect of the TB transmission dynamics is monitored by simulating the model with the parameters from Table 2. In Figure 4, TB transmission decreases gradually with the time where $R_0<1$ (with $\beta =0.9$ , $\beta =0.09$ , $\beta =0.09$ ). This transmission has to observe as a latent stage with a certain time. In Figure 5, one infected TB patient takes medicine at the rate $\left(r_1=0.5,r_2=0.5\right)$ in a proper time, his infection is decreases rapidly and the infection can be eradicated from the community within very short time. If he takes treatment at the moderate rate $\left(r_1=0.05,r_2=0.05\right)$ , he cures from the disease gradually day by day with time. On the other hand, if he takes treatment at the low rate $\left(r_1=0.005,r_2=0.005\right)$ , he will be cure from the disease at very slow rate but he gets remove from TB disease at a certain period.

In Figure 6, it is shown that when $R_0=0.8544<1$ , TB infected person gets remove from tuberculosis rapidly after a certain period with proper treatment and TB disease is stable for the value of $R_0$. On the other hand, in Figure 7, it is shown that when $R_0=4.7466>1$ , TB infected person spreads there infection at a stable stage with a certain time and TB disease is unstable for the value of $R_0$ and the disease can’t eradicate from the community permanantly.

Figure 4. Diagram of TB transmission rate.

Figure 5. Diagram of TB treatment.

Figure 6. Simulation of the model showing the total number of infected individuals as a function of time, using the parameters in Table 2 where $R_0=0.8544<1$.

Figure 7. Simulation of the model showing the total number of infected individuals as a function of time, using the parameters in Table 2 where $R_0=4.7466>1$.

8. Result

We rigorously analyzed (mathematically and numerically) the dynamics of TB in the model. Some mathematical and epidemiological findings of the study are given below:

1) The model has a disease free equilibrium (DFE) which is asymptotically stable $R_0<1$ and unstable if $R_0>1$. The model is also globally asymptotically stable for a special case when $r_1=r_2=0$.

2) When $R_0=0.8544<1$ , the rate of infected individuals increases and after a certain time it smoothly decreases.

3) And lastly, the prevalence is very high when $R_0=4.7466>1$.

9. Conclusions

A deterministic model for the transmission dynamics of TB in the population level is designed and rigorously analyzed. Some of the main findings of the study include the following:

1) The model exhibits a phenomenon of backward bifurcation, when DFE is locally asymptotically stable.

2) The model has an EEP which is globally asymptotically stable for special case (i.e. $r_1=r_2=0$ ).

The model does not undergo backward bifurcation in the absence of treatment stage.

Conflicts of Interest

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

References

[1]	Kuddus, A., Rahman, A., Talukder, M.R. and Haque, A. (2014) A Modified Sir Model to Study on Physical Behavior among Smallpox Infective Population in Bangladesh. American Journal of Mathematics and Statistics, 4, 231-239.
[2]	Side, S., Sanusi, W., Aidid, M.K. and Sidjara, S. (2016) Global Stability of SIR and SEIR Model for TB Disease Transmission with Lyapunov Function Method. Asian Journal of Applied Sciences, 9, 87-96. https://doi.org/10.3923/ajaps.2016.87.96
[3]	Dago, M.M., Ibrahim, M.O. and Tosin, A.S. (2015) Stability Analysis of a Deterministic Mathematical Model for Transmission Dynamics of TB. Proceedings of 32nd The IIER International Conference, Dubai, UAE, 8 August 2015, 42-45.
[4]	Adewale, S.O., Podder, C.N. and Gumble, A.B. (2009) Mathematical Analysis of a TB Transmission Model with Dots. Canadian Applied Mathematics Quarterly, 17.
[5]	Magombedze, G., Garira, W. and Mwenje, E. (2006) Modelling the Human Immune Response Mechanics to Mycobacterium Tuberculosis Infection in the Lungs. Mathematical Biosciences & Engineering, 3, 661-682. https://doi.org/10.3934/mbe.2006.3.661
[6]	Murray, G.J. and Salomon, J.A. (1998) Expanding the WHO Tuberculosis Control Strategy Re-Thinking the Role of Active Case Finding. International Journal of Tuberculosis and Lung Disease, 2, S9-S15.
[7]	Guo, H.B. and Li, M.Y. (2006) Global Stability in a Mathematical Model of Tuberculosis. Canadian Applied Mathematics Quarterly, 14, 185-197.
[8]	Gumel, A.B. and Song, B. (2008) Existence of Multiple-Stable Equilibria for a Multi-Drug-Resistant Model of Mycobacterium Tuberculosis. Mathematical Biosciences & Engineering, 5, 437-455. https://doi.org/10.3934/mbe.2008.5.437
[9]	Dubos, R. and Dubos, J. (1952) The White Plogue: Tuberculosis Man and Society. Little and Brown Press, Boston, MA.
[10]	Loweell, A., et al. (1996) Tuberculosis. Harvard University Press, Cambridge, MA.
[11]	Blower, S., Mclean, A.R., Porco, T., Moss, A.R., et al. (1995) The Intrinsic Transmission Dynamics of Tuberculosis Epidemics. Nature Medicine, 1, 815-822. https://doi.org/10.1038/nm0895-815
[12]	WHO (2004) Global TB Control WHO Reports.
[13]	World Health Organization (2013)
[14]	Jakubowiak, W., Bogorods Kaya, E., Borisov, S., Danilova, I. and Kourbatova, E. (2009) Treatment Interruptions and Duration Associated with Default among New Patients with Tuberculosis in Six Regions of Russia. International Journal of Infectious Diseases, 13, 362-368. https://doi.org/10.1016/j.ijid.2008.07.015
[15]	Feng, Z., Iannelli, M. and Milner, F.A. (2001) A 2 Strain TB Model with Age of Infection. SIAM Journal on Applied Mathematics, 62, 1634-1656. https://doi.org/10.1137/S003613990038205X
[16]	World Health Organization (2015)
[17]	Cohen, T., Colijin, T.C., Finklea, B. and Murray, M. (2007) Exogenous Re-Infection and the Dynamics of TB Epidemics: Local Effects in a Network Model Transmission. Journal of the Royal Society Interface, 4, 523-531. https://doi.org/10.1098/rsif.2006.0193
[18]	Castillo Chavez, C., et al. (2005) A Model of TB with Exogenous Re-Infection. Theoretical Population Biology, 132, 235-239.
[19]	Ferguson, N.M., Cummings, D.A.T., Fraser, C., Cajka, J.C., Cooley, P.C. and Burke, D.S. (2006) Strategies for Mitigating an Influenza Pandemic. Nature, 442, 448-452. https://doi.org/10.1038/nature04795
[20]	Okuonghae, D. and Korobeinikov, A. (2007) Dynamics of Tuberculosis: The Effect of Direct Observation Therapy Strategy (DOTS) in Nigeria. Mathematical Modelling of Natural Phenomena, 2, 101-113. https://doi.org/10.1051/mmnp:2008013
[21]	Sharomi, O., Podder, C.N., Gumel, A.B. and Song, B. (2008) Mathematical Analysis of the Transmission Dynamics of HIV/TB Coinfection in the Presence of Treatment. Mathematical Biosciences & Engineering, 5, 145-174. https://doi.org/10.3934/mbe.2008.5.145
[22]	Nayeem, J. and Podder, C.N. (2014) A Mathematical Study on the Vaccination Impact on the Disease Dynamics of HBV. IOSR Journal of Mathematics, 10, 26-44. https://doi.org/10.9790/5728-10612644