1. Introduction

Tuberculosis (TB) is an infectious disease whose etiological agent is the Mycobacterium tuberculosis (Mtb), in 2015 was one of top 10 causes of death worldwide and it is the second leading cause of death due to communicable diseases, preceded only by the human immunodeficiency virus (HIV). From 2014 to 2015 the rate of decline in TB incidence was $1.5\%$, although the battle against the tuberculosis epidemic is gaining WHO established that is necessary to accelerate to a $4-5\%$ annual decline by 2020 to reach the first milestones of the End TB Strategy [37].

In 2015, there were an estimated 10.4 million new (incident) TB cases worldwide, 1.4 million TB deaths and an additional 0.4 million deaths resulting from TB disease among people living with HIV. Without treatment, the death rate from TB is high. Studies of the natural history of TB disease in the absence of treatment with anti-TB drugs (that were conducted before drug treatments became available) found that about $70\%$ of people with sputum smear-positive pulmonary TB died within 10 years, as did about $20\%$ of people with culture-positive (but smear-negative) pulmonary TB [37].

Most of the infected individuals with Mtb are capable to control the infection and remain in a latent stage in which they cannot transmit the disease. It is estimated that about one third of the world population has latent TB, and around 10$\%$ of the infected population develop the active form of the illness, whether in short-term (primary infection) or long-term (reactivation) [36]. Factors associated to reactivation are malnutrition, diabetes, tobacco, HIV and immune compromised situations.

The Mtb bacteria may affect different tissues of the organism, but the most common form of the disease is pulmonary TB. In the lung, Mtb is restricted to discrete sites of infection called granulomas which are well-organized, dynamical structures formed at the site of the bacteria and mediated by specific immune responses during the infection process. The granuloma formation process starts shortly after infection, when the inhaled Mtb is ingested and transported across the alveolar epithelium into the lung tissue and adjacent lymph nodes.

A granuloma is composed of immune cells at various stages of differentiation with the infected macrophages forming the centre of the cellular accumulation. The recruited T cells secrete cytokines that activate infected cells to control their mycobacterial load and activate cytotoxic T cells. The cellular composition of TB granulomatous lesions includes blood-derived infected and uninfected macrophages, foamy macrophages, epithelioid cells (uniquely differentiated macrophages), and multinucleated giant cells (Langerhans cells), B and T lymphocytes, and fibroblasts [10,25,26].

A characteristic of the granulomas is the formation of caseous centre containing necrotic tissue, cell debris and killed mycobacteria. Bacteria are found within macrophages (intracellular bacteria) and within the zone between the necrotic centre and the cellular wall of the granuloma (extracellular bacteria) [7].

The value of an experimental model of mycobacterial persistence has at least two-fold: to uncover fundamental processes associated with clinical latency, and to guide new interventions, diagnostics, antibiotics, and vaccines, to detect, manage, and prevent disease [3,32]. Unfortunately, despite extensive studies on the interactions between Mtb and macrophages, and the granuloma formation, the mechanisms by which pathogen evades anti-microbial responses and establishes persistence within the host cell is not well understood [10].

Numerous theoretical studies have been done to understand the information available on Mtb infection, both from the point of view epidemiological as immune. To this respect, diverse mathematical models have been proposed to assess the impact on the infection progression of factors like Mtb population dynamics, immune system, treatment, and bacterial resistance. See for example [1,2,4,5,6,8,18,35,22,21,24,29,30].

In particular, in this work we propose a model to evaluate the impact of Mtb growth on the outcome of infection. For this end, we formulate a model that takes into account two ways of bacterial growth, the first one results of the average number of bacteria produced in the interior of an infected macrophage, and the other of logistic type that takes into account the competition to infect a macrophage between the outer bacteria. This model is a continuation of previous studies given in [13,14,15,16], and its content is organized in the following way: in the second section we formulate the mathematical model. In the third and fourth sections we do the qualitative analysis of the model. Finally, in the fifth, sixth and seventh sections we present the sensitivity analysis, numerical results and the discussion.

---

2. Model formulation

Following [13], we denote by $\bar{M}_U(t)$, $\bar{M}_I(t)$, $\bar B(t)$, and $\bar{T}(t)$ the populations densities at time $t$ of non infected macrophages, infected macrophages, bacilli Mtb, $T$ cells, respectively.

We assume that non infected macrophages are recruited at a constant rate $\Lambda_U$, become infected at a rate $\beta \bar M_U\bar B$, and are removed at a per capita constant rate $\mu_U$.

$\frac{d\bar{M}_U}{dt} = \Lambda_U- \mu_U\bar{M}_U-\bar \beta \bar B\bar{M}_U.$

(1)

Infected macrophages grow at a rate $\beta \bar M_U\bar B$, die at a per capita constant rate $\mu_I\geq\mu_U$, and are eliminated by $T$ cells at rate proportional to the product of $\bar{M_I}$ and $\bar{T},$ with proportionality constant $\bar{\alpha}_T$.

$\frac{d\bar{M}_I}{dt} = \bar \beta \bar B\bar{M}_U-\bar{\alpha}_T\bar{M}_I\bar{T}- \mu_I\bar{M}_I$

(2)

Mtb population grows inside of an infected macrophage up to a limit where the macrophage dies and releases bacteria. Accordingly to this, we assume that growth rate of Mtb inside macrophages is given by $\bar r\mu_I\bar M_I$, where $\bar r$ is the average number of bacilli produced by an infected macrophage. The released bacteria begin to spread outside the macrophage and start to compete among themselves by infecting new macrophages. Therefore, we assume that outside the macrophage, the Mtb population has a logistic growth with intrinsic reproduction rate $\nu$, and carrying capacity $K$. Finally, we assume that bacteria die at a per capita constant rate $\mu_B \geq\nu$, and that non infected macrophages eliminate Mtb at a rate $\bar{\gamma_U}\bar{M}_U\bar B$, with proportiona constant $\bar{\gamma_U}$.

$\frac{d \bar B}{dt} = \bar r\mu_I\bar M_I+ \nu\left(1-\frac{\bar B}{K}\right)\bar B-\bar\gamma_U\bar{M}_U \bar B- \mu_{B}\bar B$

(3)

$T$ cells are recruited up to a maximum number $T_{max}$ at a proportional rate to the number of infected macrophages with proportionality constant $k_I$, and they die at a per capita constant rate $\mu_T$.

$\frac{d\bar{T}}{dt} = \bar{k}_I\left(1-\dfrac{\bar{T}}{T_{max}}\right)\bar{M}_I-\mu_T\bar{T}.$

(4)

Figure 1 shows the flow diagram of macrophages, T cells and bacteria described in the differential equations (1)-(4). In order to reduce the number of parameters we introduce the following change of variables:

$M_U = \frac{\bar M_U}{\Lambda_U/\mu_U}, \, \, M_I = \frac{\bar M_I}{\Lambda_U/\mu_U}, \, \, B = \frac{\bar B}{K}, \, \, T = \frac{\bar T}{T_{max}}.$

(5)

Replacing the new variables the system (1)-(4) becomes

$\frac{dM_U}{dt} = \mu_U-\mu_UM_U-\beta BM_U \\ \frac{dM_I}{dt} = \beta BM_U-\alpha_TM_IT-\mu_IM_I$

$\frac{dB}{dt} = rM_I+ \nu \left(1-B\right)B-\gamma_UM_UB-\mu_{B}B \\ \frac{dT}{dt} = k_I\left(1-T\right)M_I-\mu_TT,$

(6)

where

$\alpha_T = \bar\alpha_TT_{max}, \;\;\; \beta = \bar \beta K, \;\;\; \gamma_U = \bar\gamma_U\dfrac{\Lambda_U}{\mu_U}, \;\;\; r = \dfrac{\bar r}{K}\mu_I\dfrac{\Lambda_U}{\mu_U}, \;\;\; k_I = \frac{\bar k_I\Lambda_U}{\mu_U}.$

(7)

It is a simple matter to verify that system (6) satisfies the existence and uniqueness conditions. Moreover, the region of biological interest is given by

$\label{omega2} \Omega = \left\{\left( \begin{array}{c} M_U \\ M_I \\ B \\ T \\ \end{array} \right) \in\mathbb{R}^4:0\leq M_U, M_I\leq 1, 0\leq M_U+M_I\leq 1, 0 \leq B\leq B_M, 0\leq T\leq T_c\right\},$

(8)

where $T_c = \dfrac{k_I}{k_I+\mu_T}$, and $B_M = \dfrac{1+\sqrt{1+4r/\nu}}{2}.$.

The following lemma establishes that system (6) is well posed in the sense that solutions with initial conditions in $\Omega$ remain there for all $t\ge 0$.

Lemma 2.1. The set $\Omega$ defined by (8) is positively invariant for system (6).

The proof is similar to the one given in Lemma 1 of [14].

---

3. Equilibrium solutions

The equilibria of system (6) are given by the solutions of the following algebraic system

$\mu_U-\mu_UM_U-\beta BM_U = 0 \\ \beta BM_U-\alpha_TM_IT-\mu_IM_I = 0 \\ rM_I+ \nu\left(1-B\right)B-\gamma_UM_UB-\mu_{B}B = 0 \\ k_I\left(1-T\right)M_I-\mu_TT = 0.$

(9)

It is clear that $P_0 = (1, 0, 0, 0)$ is an equilibrium of system (6) which represents the state of non infection. The solutions of (9) with $B\neq 0$ are called the bacteria-present equilibria, and correspond to the chronic infection state. In order to find these equilibria, we start to solving $M_U$ and $T$ from the first and fourth equations of (9):

$M_U = \frac{\mu_U}{\mu_U+\beta B}\, \, \, \text{and}\, \, \, T = \frac{k_IM_I}{k_IM_I+\mu_T}.$

(10)

Replacing $M_U$ defined by (10) in the second equation of (9) we get

$M_I = \frac{\beta B\mu_U}{(\mu_U+\beta B)(\alpha_TT+\mu_I)},$

(11)

which is equivalent to

$M_I = \left(\frac{\beta B}{\mu_U+\beta B}\right)\left(\frac{\mu_I}{\alpha_TT+\mu_I}\right)\frac{\mu_U}{\mu_I}.$

We observe that $0\leq M_U\leq1$ and $0\leq T\leq1$, and since $\mu_I\geq\mu_U$, then $0\leq M_I\leq1$. Replacing the expression for $T$ given by (10) in (11) we get

$M_I = \frac{\beta B\mu_U(k_IM_I+\mu_T)}{(\mu_U+\beta B)\left[(\alpha_T+\mu_I)k_IM_I+\mu_I\mu_T\right]}.$

(12)

From (12) we obtain the following quadratic equation

$M_I^2+b(B)M_I-c(B) = 0,$

(13)

where

$b(B) = \frac{\mu_I\mu_T}{(\alpha_T+\mu_I)k_I}-\frac{\beta B\mu_U}{(\alpha_T+\mu_I)(\mu_U+\beta B)} \\ c(B) = \frac{\mu_T\beta B\mu_U}{k_I(\alpha_T+\mu_I)(\mu_U+\beta B)}.$

(14)

Since $c(B)>0$, the only positive solution of (13) is

$M_I = \bar g_1(B) = \frac{-b(B)+\sqrt{[b(B)]^2+4c(B)}}{2}.$

(15)

Now, replacing $M_U$ defined in (10) in the third equation of (9), and solving for $M_I$ we obtain

$M_I = \frac{\left[\beta\nu B^2+\left(\mu_U\nu-\beta\nu+\beta\mu_B\right)B+\mu_U(\gamma_U+\mu_B)(1-R_0)\right]B}{r(\mu_U+\beta B)},$

(16)

with

$R_0 = \frac{\nu}{\gamma_U+\mu_B}.$

(17)

Equation (16) can be written as

$M_I = \bar g_2(B) = \frac{\beta\nu \left[ B^2+\dfrac{\gamma_U}{\nu}\left(\sigma-\sigma_c\right)B+\mu_U \left(\beta R_0\right)^{-1}(1-R_0)\right]B}{r(\mu_U+\beta B)},$

(18)

with

$\sigma = \frac{\nu\mu_U}{\gamma_U\beta} \, \, \text{and}\, \, \sigma_c = \frac{\nu-\mu_{B}}{\gamma_U}.$

(19)

Let us observe that $\bar g_1$, and $\bar g_2$ are two different expression for $M_I$, to study the intersections of above functions is important to determine the intersections of the functions:

$\label{g1g2} g_1(B) = r(\mu_U + \beta B)\bar g_1(B) \\ g_2(B) = \left[\beta\nu B^2+\left(\mu_U\nu-\beta\nu+\beta\mu_B\right)B+\mu_U(\gamma_U+\mu_B)(1-R_0)\right]B.\\ $

(20)

From (20) we obtain that $g_1(0) = g_2(0)$, in addition since $\bar g_1$ defined in (15) is positive and strictly increasing, then $g_1$ is positive, strictly increasing, and concave. The following proposition establishes some obvious properties of $g_1$ and $g_2$.

Proposition 1. The functions $g_1$ and $g_2$ intersect in $B = 0$. Also, $g_1$ is positive, strictly increasing, and concave in the first quadrant.

The roots of the cubic polynomial $g_2$ are $B = 0$, and

$B^{\pm} = \frac{-\dfrac{\gamma_U}{\nu}\left(\sigma-\sigma_c\right) \pm\sqrt{\left[\frac{\gamma_U}{\nu}\left(\sigma-\sigma_c\right)\right]^2-\dfrac{4\mu_UR_0^{-1}(1-R_0)}{\beta}}}{2}.$

(21)

Furthermore, we see that the derivatives of $g_1$ y $g_2$ are given by

$g_1'(B) = r\beta g_1(B) + \frac{r(\mu_U + \beta B)}{2}\left[-b'(B)+\frac{b(B)b'(B)+2c'(B)}{\sqrt{[b(B)]^2+4c(B)}}\right] \\ g_2'(B) = \beta \nu \left[3B^2 + 2\frac{\gamma_U}{\nu}(\sigma-\sigma_c)B + \frac{\mu_U}{\beta}R_0^{-1}(1-R_0)\right].$

(22)

From (22) we obtain

$g_1'(0) = \frac{r\mu_U c'(0)}{b(0)} = \frac{r\mu_U\beta}{\mu_I} \\ g_2'(0) = \nu\mu_U\left(\frac{1}{R_0}-1 \right).$

(23)

Observe that

$g'_1(0)-g'_2(0) = \mu_U (\gamma_U + \mu_B)(R_0+R_1-1),$

where

$R_1 = \frac{r\beta}{\mu_I(\gamma_U+\mu_B)}.$

(24)

In order to have a biological interpretation of the existence results for the bacteria-present equilibria in terms of dimensionless variables, in addition to $R_0$, and $R_1$, we introduce the following parameters (the significance of each of these parameters is given in section 3.1).

$R_B = \frac{\nu}{\mu_B}, ~~~~ R_{\beta} = \frac{\beta}{\mu_I}, ~~~~ R_{\gamma_U} = \frac{\gamma_U}{\mu_B}.$

(25)

In terms of the above parameters, $\sigma$ and $\sigma_c$ defined in (19) can be rewritten as

$\sigma = \frac{\mu_U}{\mu_I}\frac{R_B}{R_{\beta}R_{\gamma_U}}, \, \, \text{and}\, \, \sigma_c = \frac{R_B-1}{R_{\gamma_U}}.$

(26)

Furthermore, when $R_B>1$ we have the following results:

1. $\sigma < \sigma_c$ is equivalent to $R_{\beta}>\rho$

2. $\sigma = \sigma_c$ is equivalent to $R_{\beta} = \rho$

3. $\sigma>\sigma_c$ is equivalent to $R_{\beta} < \rho$,

where

$\rho = \frac{\mu_U}{\mu_I}\frac{R_B}{R_B-1}.$

To analyze the existence of bacteria-present equilibria, we consider two cases, $R_0\ge 1$, and $R_0 < 1$. For the first case we have the following result:

Proposition 2. If $R_0\geq 1$, there exists a unique bacteria-present equilibrium.

Proof. Assume first $R_0 >1$. From (21) we see that in this case the non zero roots of $g_2$ satisfy $B^+>0$ and $B^- < 0$. Also, $g_1$ is concave and $g_2$ is convex in the first quadrant, $g'_1(0)>0$, $g'_2(0) < 0$, $\lim_{B\to\infty}g_1(B) = \infty,$ and $\lim_{B\to\infty}g_2(B) = \infty.$ All these conditions imply that $g_1$ and $g_2$ intersect only once in the positive quadrant. For $R_0 = 1$ we have

$B^{\pm} = \frac{-\dfrac{\gamma_U}{\nu}\left(\sigma-\sigma_c\right) \pm\left|\dfrac{\gamma_U}{\nu}\left(\sigma-\sigma_c\right)\right|}{2}.$

Then, for $\sigma < \sigma_c$, $B^{+} = 0$ and $B^{-}>0$, when $\sigma = \sigma_c$, $B^{\pm} = 0$ and for $\sigma>\sigma_c$, $B^{+} < 0$ and $B^{-} = 0$. In all cases it is easy to verify form the qualitative behavior of $g_1$ y $g_2$ that these functions have only one positive intersection for $R_0 = 1$.

In the case $R_0 < 1$ we have the following results.

Proposition 3. Assume $R_0 < 1$ and $R_{\beta} \le \rho$. If $R_0+R_1 \ge 1,$ there is a unique bacteria-present equilibrium, an none if $R_0+R_1 < 1$.

Proof. The assumptions of the proposition are equivalent to $1- R_0 > 0$, and $\sigma-\sigma_c \ge 0$, which in turn imply that $g_2$ is positive and strictly increasing in the first quadrant. It follows that both $g_1$ and $g_2$ are strictly increasing positive functions with $g_1(0) = g_2(0) = 0$, and the existence of a unique intersection in the first quadrant will depend only on their derivatives evaluated at $B = 0$. Since $g_1$ and $g_2$ are respectively convex and concave functions, we conclude that $g_1$ and $g_2$ have a unique positive intersection if and only if $g_1'(0)\geq g_2'(0)$. Since this inequality is equivalent to $R_0 + R_1\geq1$, we have proved the proposition.

In the following we will assume $R_B >1$, $R_{\beta} >\rho$ (or equivalent $\sigma < \sigma_c)$, and $R_0 < 1$. Now, $B^{\pm}$ defined in (21) are positive real zeros if and only if

$\left[\dfrac{\gamma_U}{\nu}\left(\sigma-\sigma_c\right)\right]^2-\dfrac{4\mu_UR_0^{-1}(1-R_0)}{\beta}\geq 0,$

or equivalently, $R_0\geq R_0^*$, where

$R_0^* = \frac{1}{1+\frac{\beta}{4\mu_U}\left[\frac{\gamma_U}{\nu}(\sigma-\sigma_c)\right]^2}.$

When $\sigma < \sigma_c,$ and $R_0 < 1$, then $g'_2(0)>0$ and $g''_2(0) = 2\beta \gamma_U(\sigma-\sigma_c) < 0,$ which implies $g_2$ increasing and concave in $B = 0$. Also since $g_2$ intersects the B-axis in the positive quadrant, it gets its maximum and minimum in two positive values of $B$ denoted by $B^{\max}$ and $B^{\min}$ respectively. This implies the following result:

Proposition 4. Assume $R_0^* \le R_0 < 1$ and $R_{\beta}> \rho$ then

1. If $R_0+R_1\leq 1,$ $g_1$ and $g_2$ intersect twice in the positive quadrant, and therefore there exist two bacteria-present equilibria.

2. If $R_0+R_1>1$ and $g_1(B^{\max})>g_2(B^{\max})$, then $g_1$ and $g_2$ have only one positive intersection, and therefore only one bacteria-present equilibrium.

3. If $R_0+R_1>1$ and $g_1(B^{\max}) < g_2(B^{\max})$, then $g_1$ y $g_2$ have three positive intersections, and therefore three bacteria-present equilibria.

Now, if $R_0 < R_0^*$ the roots $B^{\pm}$ are complex, however since $g_2$ is increasing and concave in $B = 0$, then it gets its local maximum and minimum $B^{\max}$ and $B^{\min}$ respectively. This implies similar results to the ones in Proposition 4.

Proposition 5. Assume $R_0 < R_0^* < 1$, $R_B>1$, and $R_{\beta}> \rho$ then

1. If $R_0+R_1\leq1$ and $g_1(B^{\min}) < g_2(B^{\min})$, then $g_1$ and $g_2$ have no positive intersection, and therefore there is no bacteria-present equilibrium.

2. If $R_0+R_1\leq1$, $g_1(B^{\min})>g_2(B^{\min})$ and $g_1(B^{\max}) < g_2(B^{\max})$, $g_1$ y $g_2$ intersect twice in the positive quadrant, and there are two bacteria-present equilibria.

3. If $R_0+R_1>1$ and $g_1(B^{\min}) < g_2(B^{\min})$, then $g_1$ and $g_2$ have only one positive intersection, and therefore only one bacteria-present equilibrium.

4. If $R_0+R_1>1$, $g_1(B^{\min})>g_2(B^{\min})$ and $g_1(B^{\max}) < g_2(B^{\max})$, then $g_1$ y $g_2$ have three positive intersections and, therefore three bacteria-present equilibria.

5. If $R_0+R_1>1$ and $g_1(B^{\max})>g_2(B^{\max})$, then $g_1$ and $g_2$ intersect in one positive point, and there is only one bacteria-present equilibrium.

---

3.1. Interpretation of the bifurcation parameters

In this section we will give a biological interpretation of the parameters $R_B$, $R_{\beta}$, $R_{\gamma_U}$, $R_0$ and $R_1$ defined in the section above.

In the formulation of the model is not consider the explicit distinction between internal and external bacteria. However, for the purposes of interpretation we will denote interior (exterior) bacteria the ones in the interior (exterior) of the infected macrophages. In this sense, the product between the average number of bacteria produced by an infected macrophage, $\bar r$, times the rate of infection $\bar\beta$ is interpreted as the growth rate of the interior or intracellular bacteria, while the logistic growth rate $\nu$ is the the growth rate of external or extracellular bacteria. Under these considerations, and knowing that $1/\mu_B$ is the average life expectancy of bacteria,

$R_B = \frac{\nu}{\mu_B}$

(27)

represents the average number of bacteria generated by an exterior Mtb.

On the other hand, in a healthy organism, the population of uninfected macrophages is given by $\Lambda_U/\mu_U,$ and this population eliminates invasive bacteria at a rate $\gamma_U = \bar\gamma_U\Lambda_U/\mu_U$, thus

$R_{\gamma_U} = \frac{\gamma_U}{\mu_B} = \frac{\bar\gamma_U\frac{\Lambda_U}{\mu_U}}{\mu_B},$

(28)

represents the average number of invasive bacteria eliminated during their lifetime. Therefore, this number is a measure of the effectiveness of macrophages in controlling bacteria.

Now, once infected, a macrophage on average generates $K/\mu_I$ bacteria, which in turn infect $\bar \beta K/\mu_I$ macrophages. Therefore, the parameter

$R_{\beta} = \frac{\bar\beta K}{\mu_I} = \frac{\beta}{\mu_I},$

(29)

is named the Basic Reproductive Number of the infection since it represents the average number of infected macrophages derived from one infected macrophage when bacteria is introduced for the first time into the organism.

We notice that the parameter $R_0$ defined in (19) can be rewritten in terms of $R_B$ as

$R_0 = \frac{\nu}{\gamma_U+\mu_B} = F(\gamma_U)R_B,$

(30)

where

$F(\gamma_U) = \frac{\mu_B}{\gamma_U+\mu_B} = 1-\frac{\gamma_U}{\gamma_U+\mu_B}.$

(31)

Since $\gamma_U = \bar\gamma_U\frac{\Lambda_U}{\mu_U}$ represents the rate at which non infected macrophages eliminate bacteria, then can be interpreted as the fraction of invasive bacteria eliminated by macrophages. From (31) we conclude that $F$ is equal to the bacteria fraction that survive macrophage attack. Therefore, $R_0$ represents the bacteria produced by the fraction of external bacteria that survive to macrophages attack, and it is called the Associate Reproductive Number.

Finally, we see that

$R_1 = \frac{r\beta}{\mu_I(\gamma_U+\mu_B)} = \frac{\bar r\bar\beta K }{\gamma_U+\mu_B}\frac{\Lambda_U}{\mu_U},$

can be interpreted as the bacteria produced by the fraction of internal bacteria that survive to the control of the population of infected macrophages at equilibrium.

---

4. Stability of equilibrium solutions

In this section we analyze conditions for stability of the equilibrium points. For this, we calculate the eigenvalues relative to the Jacobian of system (6) evaluated at the equilibrium points, given by

$J\left( \begin{array}{c} M_U \\ M_I \\ B \\ T \\ \end{array} \right) = \left( \begin{array}{cccc} -(\mu_U+\beta B) &0 & -\beta M_U & 0 \\ \beta B & -(\alpha_TT+\mu_I) & \beta M_U & -\alpha_TM_I \\ -\gamma_UB & r & a & 0\\ 0& (1-T)k_I & 0 & -(k_IM_I+\mu_T) \\ \end{array} \right),$

(32)

where

$a = \nu\left(1-\frac{2B}{K}\right)-\gamma_UM_U-\mu_B.$

(33)

For the infection free equilibrium $P_0 = (1, 0, 0, 0),$ the Jacobian is given by

$J\left(P_0 \right) = \left( \begin{array}{ccccc} -\mu_U &0 & -\beta & 0 \\ 0 & -\mu_I & \beta & 0 \\ 0 & r & \nu-(\gamma_U+\mu_B) & 0\\ 0& k_I & 0 & -\mu_T \\ \end{array} \right).$

(34)

Simple calculations show that the eigenvalues are given by $\lambda_1 = -\mu_U$, $\lambda_2 = -\mu_T,$ and the roots of the quadratic equation

$\lambda^2+[\mu_I+\gamma_U+\mu_B-\nu]]\lambda+\mu_I(\gamma_U+\mu_B)[1-(R_0+R_1)] = 0$

or equivalently

$\lambda^2+[\mu_I+(\gamma_U+\mu_B)(1-R_0)]]\lambda+\mu_I(\gamma_U+\mu_B)[1-(R_0+R_1)] = 0.$

(35)

From Routh-Hurwitz criterion we conclude that all the eigenvalues of the equation (35) have negative real part if and only if $R_0+R_1 < 1$. Therefore we have the following

Proposition 6. The infection free equilibrium $P_0 = (1, 0, 0, 0)$ is locally asymptotically stable if $R_0+R_1 < 1$, and unstable when $R_0+R_1>1$.

Now, we analyze the stability of bacteria-present equilibria which reflects the infection persistence. From the equations at equilibrium (9) we get the following equalities

$\frac{\mu_U}{M_U} = \mu_U+\beta B \\ \frac{\beta BM_U}{M_I} = \alpha_TT+\mu_I \\ \nu\left(1-2B\right)-\gamma_UM_U-\mu_B = -\left(\frac{rM_I}{B}+\nu B\right)\\ \frac{k_IM_I}{T} = k_IM_I+\mu_T.$

(36)

Replacing (36) in (34) we obtain

$J\left(P_i \right) = \left( \begin{array}{cccc} -\dfrac{\mu_U}{M_U} &0 & -\beta M_U & 0 \\ \beta B & -\dfrac{\beta BM_U}{M_I} & \beta M_U & -\alpha_TM_I \\ -\gamma_UB & r & -\left(\dfrac{rM_I}{B}+\nu B\right) & 0\\ 0& (1-T)k_I & 0 & -\dfrac{k_IM_I}{T} \\ \end{array} \right).$

(37)

To get the conditions for negative eigenvalues of $J(P_i)$, $i = 1, 2, 3$, we obtain after tedious calculations its characteristic polynomial given by

$p_1(\lambda) = \left(\lambda+\frac{\mu_U}{M_U}\right)\left(\lambda+\frac{\beta BM_U}{M_I}\right) \left(\lambda+\frac{rM_I}{B}+\nu B\right)\left(\lambda+\frac{k_IM_I}{T}\right)\\ +r\beta M_U\left(\lambda+\frac{k_IM_I}{T}\right)\left[\beta B-\left(\lambda+\frac{\mu_U}{M_U}\right)\right] \\ +\beta M_U\alpha_TM_I(1-T)k_I\left(\lambda+\frac{rM_I}{B}+\nu B\right)\left(\lambda+\frac{\mu_U}{M_U}\right)\\ -\beta M_U\gamma_UB\left[\left(\lambda+\frac{\beta B M_U}{M_I}\right) \left(\lambda+\frac{k_IM_I}{T}\right)+\alpha_TM_I(1-T)k_I\right]\\ = \lambda^4+s_1\lambda^3+s_2\lambda^2+s_3\lambda+s_4,$

(38)

where

$s_1 = \frac{\mu_U}{M_U}+\frac{\beta B M_U}{M_I}+\frac{rM_I}{B}+\nu B+\frac{k_IM_I}{T}\\ s_2 = \left(\frac{rM_I}{B}+\nu B\right)\frac{k_IM_I}{T}+\frac{\mu_U}{M_U}\frac{\beta BM_U}{M_I} +\left(\frac{\mu_U}{M_U}+\frac{\beta BM_U}{M_I}\right) \left(\frac{rM_I}{B}+\nu B+\frac{k_IM_I}{T}\right) \\ +\beta M_U\alpha_TM_I(1-T)k_I-r\beta M_U-\beta M_U\gamma_UB\\ s_3 = \frac{\mu_U}{M_U}\frac{\beta BM_U}{M_I}\left(\frac{rM_I}{B}+\nu B+\frac{k_IM_I}{T}\right) + \left(\frac{rM_I}{B}+\nu B\right)\frac{k_IM_I}{T}\left(\frac{\mu_U}{M_U}+\frac{\beta BM_U}{M_I}\right)\\ +\beta M_U\alpha_TM_I(1-T)k_I\left(\frac{rM_I}{B}+\nu B+\frac{\mu_U}{M_U}\right)+r\beta M_U\gamma_UB \\ -r\beta M_U\left(\frac{\mu_U}{M_U}+\frac{k_IM_I}{T}\right)-\beta B\gamma_UM_U\left(\frac{\beta B M_U}{M_I}+\frac{k_IM_I}{T}\right)\\s_4 = \frac{\mu_U}{M_U}\frac{\beta B M_U}{M_I}\left(\frac{rM_I}{B}+\nu B\right) \frac{k_IM_I}{T} +\alpha_TM_I(1-T)k_I\left(\frac{rM_I}{B}+\nu B \right) \frac{\mu_U}{M_U}\\ +r\beta M_U\beta B\frac{k_IM_I}{T}-\beta M_U\gamma_UB\frac{\beta B M_U}{M_I}\frac{k_IM_I}{T}-\beta M_U\gamma_UB\alpha_TM_I(1-T)k_I\\ -r\beta M_U\frac{\mu_U}{M_U}\frac{k_IM_I}{T}.$

(39)

The coefficient $s_1 >0$ since all the parameters are positive. We rewrite the other coefficients as:

$s_2 = \frac{\beta BM_U}{M_I}\left(\frac{\mu_U}{M_U}+\nu B +\frac{k_IM_I}{T}\right) + \frac{k_IM_I}{T}\left(\frac{\mu_U}{M_U}+\frac{rM_I}{B}+\nu B\right) \\ +\frac{\mu_U}{M_U}\frac{rM_I}{B}+\beta M_U\alpha_TM_I(1-T)k_I+\frac{\beta\gamma_UB}{M_U}\left(\sigma-M_U^2\right)\\ s_3 = \left(\frac{\beta BM_U}{M_I}+ \frac{k_IM_I}{T}\right) \frac{\beta\gamma_UB}{M_U}\left(\sigma-M_U^2\right)+r\beta M_U\beta B\\ +\frac{rM_I}{B}\left(\frac{\mu_U}{M_U}\frac{k_IM_I}{T}+\beta B\alpha_TM_I(1-T)k_I\right)\\ +\left[\frac{\beta BM_U}{M_I}\frac{k_IM_I}{T}+\beta B\alpha_TM_I(1-T)k_I\right] \left(\frac{\mu_U}{M_U}+\nu B\right)\\ s_4 = \left[\frac{\beta B M_U}{M_I}\frac{k_IM_I}{T}+\alpha_Tk_IM_I(1-T)\right] \frac{\beta\gamma_UB}{M_U}\left(\sigma-M_U^2\right) \\ +\alpha_TM_I(1-T)k_I\frac{\mu_U}{M_U}\frac{rM_I}{B}+r\beta M_U\beta B\frac{k_IM_I}{T},$

(40)

where $\sigma$ is given in (19). Since $s_1 >0$, the Routh-Hurwitz criteria sets that the roots of a polynomial $p_1(\lambda)$ of order four have negative real part if and only if its coefficients satisfy

$s_4 > 0 \\ D_1 = s_1>0 \\ D_2 = s_1 s_2-s_3 > 0 \\ D_3 = (s_1s_2-s_3)s_3-s_1^2s_4 >0.$

(41)

See [9], in order to determine the conditions for which the previous inequalities are satisfied; we define the following constants:

$A = \frac{\mu_U}{M_U}, ~~N = \frac{\beta BM_U}{M_I},~~ C = \frac{rM_I}{B}, ~~D = \nu B, E = \frac{k_IM_I}{T}, \\ F = \beta B\alpha_Tk_IM_I(1-T), \, \bar F = \alpha_Tk_IM_I(1-T), \, G = r\beta M_U\beta B \\ X(M_U) = AD-\beta M_U\gamma_UB = \frac{\beta\gamma_UB}{M_U}\left(\sigma-M_U^2\right).\\ $

(42)

Replacing $A$, $N$, $C$, $D$, $E$, $F$, $\bar F$, $G$ and $X(M_U)$ in $s_1, \dots, s_4$ we obtain

$s_1 = A+N+C+D+E \\ s_2 = N(A+D+E)+E(A+C+D)+AC+F+X(M_U) \\ s_3 = (N+E)X(M_U)+G+C(AE+F)+(NE+F)(A+D) \\ s_4 = (NE+\bar F)X(M_U) +\bar F A C+GE$

(43)

Furthermore, after some simplifications, $D_2$, and $D_3$ can be written as

$D_2 = (N+E)[(A+D)^2+C(C+D)+(N+E)(A+C+D)+AC+F]\\ +(A+D)[C(A+E)+X(M_U)]+C[AC+X(M_U)]\\ D_3 = r_0+r_1(EN-\bar F)+r_2(ANC-G)+r_3(ADEN-\bar FX(M_U)) \\ +AN^2(DEN-C\bar F)+(AC+2AN)(D^2EN-C^2\bar F)\\ +GN(AD-X(M_U)),$

(44)

with

$r_0 = \sum\limits_{n = 1}^{172}a_nA^{\alpha_1}N^{\alpha_2}C^{\alpha_3}D^{\alpha_4}E^{\alpha_5}F^{\alpha_6}\bar F^{\alpha_7}G^{\alpha_8}[X(M_U)]^{\alpha_9}\\ r_1 = \left[(A+D)^2+AC+2CD\right](AC+X(M_U))\\ +AC(AC+2AE+2CD+2DE+E^2+2AN+2DN+EN) \\ r_2 = s_1E^2+(AC+AN+DN+X(M_U))E+(A+C+D)F+G \\r_3 = (2A+2C+2D+E+N)E+(2A+2C+D+N),$

and $a_n\in\{1, 2, 3, 4, 5, 6\}$ for $n = 1, \ldots, 172$ and $\alpha_k\in\{0, 1, 2, 3\}$ for $k = 1, \ldots, 9$.

The following theorem summarizes the stability results of the unique bacteria-present equilibrium when $R_0>1.$

Theorem 4.1. If $R_0 >1$ there exists a unique bacteria-present equilibrium $P_1$ which is locally asymptotically stable.

Proof. The existence of a unique bacteria-present equilibrium under the hypothesis of the theorem is proved in Proposition 2. It can easily verify that the Routh-Hurwitz conditions for the coefficients given in (43) are satisfied if $X(M_U) = \sigma-M^2_U$, and $D-C = \nu-\gamma M_U- \mu_B$ are both bigger or equal to zero, therefore it is enough to show that those expressions are positive when $R_0 >1$.

Indeed, if $R_0>1$ then $\nu-\gamma_U-\mu_B >0$, and since $M_U\leq 1$ it follows that $D-C >0$. On the other hand, from $R_0>1$ we get

$\sigma = \frac{\nu-\mu_B}{\gamma_U} >\frac{\gamma_U}{\gamma_U} = 1,$

which implies $\sigma \ge M^2_U$, or equivalently $X(M_U) \ge 0$. These results prove the local stability of $P_1$.

In the following we will assume $R_0 \le 1$, and $R_0 + R_1 >1$. Propositions 3, 4, and 5 assure the existence of a unique bacteria-present equilibrium when $g_1(B^{min}) < g_1(B^{min})$, or $g_1(B^{max})>g_2(B^{max})$. Under these condition we have the following stability results:

Theorem 4.2. Assume $R_0 < 1$, $R_0 + R_1 >1$, and $g_1(B^{min}) < g_1(B^{min})$, or $g_1(B^{max})>g_2(B^{max})$, then

a) if $R_{\beta}\leq \rho$, and $R_B \ge R_{\gamma_U} +1$,

or

b) $R_{\beta} > \rho$ and $\sigma >1$,

the bacteria-present equilibrium $P_1$ is locally asymptotically stable.

Proof. As in the above theorem it is enough to show that $D-C$ and $X(M_U)$ are bigger than zero. To prove $a)$ we observe that $D-C$ can be written as

$D-C = \nu - \gamma_U M_U - \mu_B = \gamma_U(\sigma_c - M_U).$

Since $R_B \ge R_{\gamma_U} +1$ is equivalent to $\sigma_c\ge 1$ then $D-C\ge0$. On the other hand, as $R_{\beta}\leq\rho$ is equivalent to $\sigma_c\leq\sigma$, and $M_U\leq\sigma_c$ then $X(M_U)>0$. Therefore, we conclude that $P_1$ is l.a.s.

The stability of $P_1$ in the case b) follows from the assumption $\sigma >1$.

In the following we assume $R_\beta >\rho$. The next theorem summarizes the stability results when there are more than one bacteria-present equilibrium.

Theorem 4.3.1 Assume $R_{\beta}>\rho.$

Ⅰ. If $R_0^* < R_0 < 1,$ then

  a. if $R_0+R_1\leq 1$ there are two bacteria-present equilibria, one l.a.s and the other one unstable,

  b. if $R_0+R_1>1$ and $g_1(B^{\max})>g_2(B^{\max})$, there are three equilibria, two l.a.s and the third one unstable.

Ⅱ. If $R_0 < R_0^* < 1,$ then

  a. if $R_0+R_1\leq 1,$ $g_1(B^{\min})>g_2(B^{\min})$, and $g_1(B^{\max}) < g_2(B^{\max}),$ there are two bacteria-present equilibria, one l.a.s, and the other one unstable

  b. if $R_0+R_1>1,$ $g_1(B^{\min})>g_2(B^{\min})$, and $g_1(B^{\max}) < g_2(B^{\max}),$ there are three bacteria-present equilibria, two l.a.s, and the third one unstable.

Proof. We will prove the stability properties in the cases where there are three equilibria finding intervals in the $B$ axis were $X(M_U)$ takes positive and negative values. A similar argument works in the cases with two bacteria-present equilibria. For this end, let $\tilde M_U = \sqrt{\sigma}$, then $X(\tilde M_U) = 0$. Since

$M_U = f_1(B) = \frac{\mu_U}{\mu_U +\beta B}$

is continuous and strictly decreasing function of $B$, then

$\tilde{B} = \frac{\mu_U}{ \beta }\left(\frac{1}{\sqrt{\sigma}}-1\right),$

is the unique value of $B$ that satisfies $\tilde M_U = f_1(\tilde B) = \sqrt{\sigma}$. Substituting $R_0 = \nu/(\gamma_U + \mu_B)$ in the function $g_2$, and evaluating in $\tilde B$ we obtain

$g_2(\tilde B) = \beta \nu\left\lbrace\left[\displaystyle{\mu_U \over \beta} \left( \displaystyle{1 \over \sqrt{\sigma}}-1\right) \right]^2 + \displaystyle{\gamma_U \over \nu} (\sigma - \sigma_c)\displaystyle{\mu_U \over \beta}\left(\displaystyle{1 \over \sqrt{\sigma}}-1\right) + \displaystyle{\mu_U \over \nu\beta}(\gamma_U + \mu_B - \nu) \right\rbrace\tilde B$

(45)

After some simplifications (45) becomes

$g_2(\tilde B) = \displaystyle{2\mu_U \gamma_U \over \sqrt{\sigma}}\left(\sqrt{\sigma}- \displaystyle{\sigma + \sigma_c \over 2}\right) \tilde B.$

(46)

We observe that for $\sigma < < 1$, the function $g_2$ satisfies $g_2(\tilde B) < -\mu_U \gamma_U \sigma_c < 0$. In addition, $g_2$ has one root in $B = 0$ and two positive roots $B^-$ and $B^+$, that satisfy $0 < B^- < B^+$. Since $g_2$ is positive in $(0, B^-)$, and negative in $( B^-, B^+)$ then $g_2$ reaches a maximum $B^{max} \in (0, B^-)$, and a minimum $B^{min} \in ( B^-, B^+)$. It follows that $B^- < \tilde B < B^+$. Now, since $f_1$ is continuous and strictly decreasing, there are positive $M_U^-$ and $M_U^+$ such that $M_U^- = f(B^-)$ and $M_U^+ = f(B^+)$, with $M_U^+ < \tilde M_U < M_U^-$. Furthermore, it can be verified that for all $B\in (0, \tilde B)$ there exists $M_U < \tilde M_U$ such that $X(M_U) < 0$, and for all $B\in(\tilde B, \infty),$ $M_U > \tilde M_U$ such that $X(M_U) < 0$. Now, let $B_1 < B_2 < B_3$ the points where $g_1$ and $g_2$ intersect, then from the properties of these functions, it is verified that $B_1 < \tilde B, B_2 < \tilde B$ and $B_3 > \tilde B$, and it follows that there are $M_{1U} = M_U(B_1) > \tilde M_U$, $M_{2U} = M_U(B_2) > \tilde M_U$, and $M_{3U} = M_U(B_3) > \tilde M_U$, which implies that $X(M_{1U})>0$, $X(M_{2U})>0$, $X(M_{3U})>0$, (see Figure 2). This concludes the proof of the Theorem.

---

5. Sensitivity analysis

The results of model (6) depend of several parameters, hence it is expected that uncertainities arise in the numerical estimates of those parameters which affect the model results. In this section we are interesting to perform global sensitivity analysis of the parameters related to the bacterial growth, infection rate, and elimination by macrophages ($\nu$, $\mu_{B}$, $\bar \beta$, $\Lambda_U/\mu_U$, $\bar \gamma_U,$ $\bar r$) in order to quantify the impact of their variations in the model outcome. For this end we follow [22,27] to carry out a global sensitivity analysis using latin hypercube sampling (LHS) to account for the effect of the uncertainties using $R_0$ and $R_1$ given by (17), and (24) as the response functions, and the parameter ranges of Table 1.

We sampled the space of the input values using LHS with a uniform probability distribution. In LHS, each parameter probability distribution is divided into $N$ equal intervals, and sampling from each interval exactly once guarantees that the entire parameter space is explored. Furthermore, a Monte Carlo simulation was done by drawing $N = 10000$ independent parameters set with $i = 1, ..., N$, and evaluating $R_0$ and $R_1$ for the corresponding parameter set. Assuming that the relation between output and input is linear, a linear regression model is used to assess the $R_i$, $i = 0, 1$ sensitivity to each parameter.

Figure 3 shows the standard regression coefficients (SCR) for $R_0 = \displaystyle{\nu \over \gamma_U + \mu_B}$ with $\gamma_U = \bar \gamma_U\displaystyle{\Lambda_U \over \mu}$ assuming the parameter values and ranges given in Table 1. In all simulations, the sensitivity based on the SCR can capture 96% of the variation on $R_0$. As it is expected $R_0$ increases when Mtb growth rate $\nu$ increases, and decreases when Mtb mortality, $\mu_{B},$ and bacteria elimination rate by macrophage population $\gamma_U$ increases. The sensitivity analysis reveals that the top parameters that play the more dominant role on the dynamics of the Mtb are the growth rate $\nu$, and death rate $\mu_B$ of Mtb bacteria. Thus, variations of the Mtb growth rate accounts for almost $60\%$ of the $R_1$ positive variation, while the death rate for almost the $80\%$ of the negative one, which implies that the bacteria dynamics have more influence on the initial evolution of the Mtb infection than the immune response due to macrophages.

Next, we quantify the impact of the variations or sensitivity of the parameters $\bar r$, $\bar \beta$, $\Lambda_U/\mu$, $\bar \gamma_U$, and $K$ on the bacteria produced and releases by the internal bacteria. For this, we perform sensitivity analysis using $R_1$ given in (24) as the response function. The results are illustrated in Figure 4. We observe that the dominant parameters are the bacteria generated by a macrophage with almost $80\%$ of the positive variation, followed by the infection rate of Mtb with almost $50\%$ of the positive variation, respectively. According to these results, $R_1$ is very sensitive to the generation and release of the Mtb bacteria by the macrophages. These findings for $R_0$, and $R_1$ could explain why the macrophage response is not enough to control an initial invasion of bacteria, and the need of the immune system to carry out a more complex defensive mechanisms to contain infections by Mtb such as the recruitment of different elements of the immune system, and the formation of granulomas.

---

6. Numerical solutions

In this section we present numerical simulations of system (6). It is important to make clear that the parameters variability depends of the immunological conditions of each patient. However, we will present some estimations based in a bibliographic revision.

In the following we verify numerically the existence of three equilibria for conditions according to the results given in last sections. Taking the parameters $\Lambda_U = 1000$, $\beta = 0.000025$, $\mu_U = 0.033$, $\mu_I = 0.11$, $\mu_B = 0.12$, $\mu_T = 0.33$, $\bar\alpha_T = 0.00003$, $\bar r = 2$, $\nu = 0.4$, $\bar\gamma_U = 0.000029$, $\bar k_I = 0.00015$, $K = 25000$ and $T_{max} = 50000$ we obtain $\sigma < \sigma_c$, $R_0 < 1$, $R_0+R_1>1$, $R_0^* < R_0$ and $g_1\left(B^{\max}\right) < g_2\left(B^{\max}\right)$, which implies by Proposition 4 that $g_1$ and $g_2$ intersect in three points, $B_1 = 0.0081$, $B_2 = 0.0765,$ and $B_3 = 0.5566$.

Therefore, the equilibrium solutions are

$P_1 = \left( \begin{array}{c} 0.8298 \\ 0.1847 \\ 0.0081\\ 0.2028 \\ \end{array} \right), P_2 = \left( \begin{array}{c} 0.4575 \\ 0.36 \\ 0.0765\\ 0.3315 \\ \end{array} \right), P_3 = \left( \begin{array}{c} 0.1099 \\ 0.4875 \\ 0.5566\\ 0.4018 \\ \end{array} \right).$

(47)

Numerical simulations confirm that $P_1$ and $P_2$ are l.a.s, and $P_3$ unstable. Even more, they suggest that the unstable branch of $P_2$ divides the stability regions of $P_1$ and $P_3$. In fact, we verified numerically that the initial condition

$P(0) = (0.457461, 0.360034, B(0), 0.331513),$

with $B(0) < 0.0765$ ($B(0)>0.0765$) is in the attraction region of $P_1$ ($P_3$), as can be seen in Figure 5 which shows the temporal course of $B(t)$ with ten initial conditions. In this case we have two stable equilibria ($P_1$ and $P_2$), and two unstable ones ($P_0$ and $P_3$). In the following simulations we show a bi-stability region in the case where there are only two bacteria-present equilibria (one stable and one unstable), and the trivial equilibrium is stable. The parameters in the simulations are $\Lambda_U = 1000$, $\bar\beta = 0.00025$, $\mu_U = 0.033$, $\mu_I = 0.1$, $\mu_B = 0.12$, $\mu_T = 0.15$, $\bar\alpha_T = 0.003$, $\bar r = 5$, $\nu = 0.4$, $\bar\gamma_U = 0.0029$, $\bar k_I = 0.01$, $K = 250000$ and $T_{max} = 50000$. It can be verified that the conditions $R_{\beta}>\rho$, $R_0^* < R_0 < 1$ and $R_0+R_1 < 1$ are satisfied, then again by Proposition 4 it follows the existence of two bacteria-present equilibria $P_1$ and $P_2$ with

$P_1 = \left( \begin{array}{c} 0.0019 \\ 0.0004 \\ 0.2679 \\ 0.4802 \\ \end{array} \right), \, P_2 = \left( \begin{array}{c} 0.0009 \\ 0.0004 \\ 0.4315 \\ 0.1 \\ \end{array} \right).$

(48)

Theorem 4.3 implies that $P_0$ and $P_2$ are l.a.s. and $P_1$ unstable. We observe this behavior in Figure 6 which illustrate the temporal course of bacteria for ten initial conditions. From the existence and stability results it follows that for $R_{\beta}\leq \rho$ the bifurcation diagram corresponds to a translation of a forward bifurcation in which the stable disease-free equilibrium $P_0$ bifurcates to the stable bacteria-present equilibrium $P_0$ in the value $R_0 = 1-R_1$ (see Figure 7). In contrast, from the existence and stability analysis of equilibria when $R_{\beta} >\rho$, the bifurcation diagrams are the ones shown in Figure 8. It is worth to notice that the bacteria-present equilibrium is stable when $R_0\geq1$, that is, all the bifurcation diagrams have the same behavior for $R_0\geq1.$

---

7. Discussion

In this work we explore the effect on the progression to TB disease due to the population growth of tuberculosis bacterium and the patient immune system control. For this end we formulated a non linear system of ordinary differential equations to describe in a simple way the interaction of Mtb with T cells and macrophages. The immune response to TB infection is a very complex phenomena that involve process of cellular differentiation and activation that have been described in several works [30,35], but due to the difficulties that involve modelling all of these process, we only considered the most important cells in the activation of the immune system against Mtb.

The model formulated in this work arises as a necessity to complement previous works given in [13,14,15]. Here we assume two forms of bacterial growth, the first one is the growth in the interior of the infected macrophages considered in previous works, and the second one is a logistic growth of external bacteria competing for the resources. Both assumptions are justified [17,26]

As it was expected, the complexity of the results increased with the assumption of logistic growth. The qualitative analysis of the model revealed different scenarios in which there is always the infection-free state, while depending on certain conditions there may be one, two or even three bacteria-present equilibria. An interesting fact is that for certain values of the parameters there are two kinds of bi-stability regions. In the first one the disease-free equilibrium and the bacteria-present equilibrium coexist, which means that depending on the initial conditions of the host and bacteria, the infection will be cleared out or will progress to TB, either in a latent or active form. In the second case the introduction of bacteria always will progress to infection, and depending on the initial conditions, the population will approach to a state with low or with high number of Mtb. The first state could be associated to latent TB, and the second one to active TB.

The above results were obtained in terms of the following parameters: ⅰ) number of bacteria generated by an external bacteria, $R_B$; number of bacteria generated by external bacteria that survive to macrophages control, $R_0$; ⅲ) number of bacteria generated by internal bacteria that survive to macrophages control at equilibrium, $R_1$; ⅳ) number of infected macrophages derived of one infected macrophage, $R_{\beta}$; and finally, ⅴ) threshold parameters $\rho$ and $R^{*}_0$ that do not have biological interpretation but are involved in the bifurcation of equilibrium solutions.

The qualitative analysis and numerical results suggest that for $R_{\beta} \le \rho$ there is a forward bifurcation in which the infection free equilibrium $P_0$ bifurcates to the bacteria-present equilibrium. Unlike the classical bifurcation which appears for $R_0 = 1$, in this case it occurs for $R_0 = 1- R_1$. This indicates the existence of a range in which the external bacteria population can not grow, but with the participation of the internal bacteria, reactivation of TB may occur in a patient with latent TB.

When $R_{\beta} > \rho$ could be occur three kind of bifurcations, one forward in the case of a unique bacteria-present equilibrium, a second one backwards [14] with an horizontal translation in which two bacteria-present equilibria coexist forming two branches of a tangential bifurcation when $0 < R_0 < 1-R_1$, and a S-shaped bifurcation when there are three bacteria-present solutions. In terms of the infection, these bifurcations indicate that depending the relation between the growth rates of external and internal bacteria as well as the bacteria elimination rate by $T$ cells and macrophages the infection can dye out, progress to latent or active TB, or to both types of TB.

On the other hand, sensitive analysis of the model parameters indicates why macrophagues are not enough to control an initial invasion by Mtb and the need of the immune system to carry out a more complex defensive mechanisms to contain infection by Mtb such as the recruitment of different elements of the immune system, and the formation of granulomas.

Concluding, in this work we proved that including competition between bacteria it is possible to obtain a greater variety of scenarios observed in the development of pulmonary tuberculosis.

---

Acknowledgments

We want to thank anonymous referees for their valuable comments that helped us to improve the paper. E. Ibarguen-Mondragón and E. M. Burbano-Rosero acknowledge support from project approved by ACUERDO No 182-01/11/2016 (VIPRI-UDENAR). Lourdes Esteva acknowledges support from project IN-112713, PAPIIT-UNAM.

---

Grant No 182-01/11/201, Vicerrectoría de Investigaciones, Posgrados y Relaciones Internacionales de la Universidad de Nariño.