1. Background

Tuberculosis (TB) is an infectious disease caused by the bacterium Mycobacterium Tuberculosis, which typically affects lungs (known as pulmonary TB of PTB) and that mostly occurs in adults, i.e. in individuals above 14 years of age, but could also affect other parts of body (known as extra-pulmonary TB or EPTB). Tuberculosis is one of the leading causes of death worldwide after HIV and remains a major public threat in many countries. An estimate shows that globally about one-third of the population is infected with TB bacteria. The global incidence of all forms of TB cases during 2008 was estimated to be 9.4 million, at the rate of 139/100,000 population and in 2014 there were an estimated 9.6 million incident cases of TB, which is equivalent to 130/100,000 population.

In this paper we present a model for the disease situation in India, one of the countries in which the disease is endemic. This work is an extension of previous results, [3,15,18] and in particular of the work in progress [19], in that we account for the diagnostic system to be partioned into the public and private sectors, contrary to what was assumed in the earlier works.

---

2. Motivation

In India, the tuberculosis situation is characterized by high prevalence (total number of TB cases over a period of one year) and incidence (new TB cases over a period of one year) of disease or active TB (when the individual is infectious and transmits bacteria to others) and high rate of transmission of infection (the latent TB-individual is infected with TB bacteria but cannot transmit them to others; also, in this situation, bacteria remain in the dormant state). Primary surveys show that about 30% -50% of India population is latently infected. This means that people have been infected by TB bacteria but are not ill with the disease and do not transmit it but may become infectious in the future.

The countrywide National Tuberculosis Program (NTP) to control TB was originally undertaken in 1962, but it did not achieve the goal of disease burden reduction. The Government of India has intensified anti-tuberculosis activities by implementing the DOTS strategy under the Revised National Tuberculosis Control Program (RNTCP) since 1998. DOTS is the WHO recommended treatment strategy to cure TB. The 2015 TB statistics show that the incidence rate is 217 [112-315] per thousand individuals in the population [16] and the prevalence rate estimate for the year 2014 is 195 [131-271] per 100,000 individuals in the population, [17]. The implementation of DOTS across the world has shown a decline of 1.5% in incidence over the past decade. In the new 2016-2020 Global Plan to end TB, "The Paradigm Shift", one of the targets is to achieve a 10% annual decline in TB incidence. The question of resistance to treatment and reactivation has been addressed in [13,14].

By formulating a model and analysing it, the objective of this paper is to give some estimates of the crucial parameters so that acting on these parameters may help in achieving the eradication of TB. These results could help in defining polices to bring down the TB burden in India. To achieve these goals the need for modelling is evident.

As it will become clear in the following sections, there is a discrepancy between the rates at which TB is diagnosed between the public and the private sectors. While one would expect the latter to be more efficient, and therefore to achieve a higher diagnosis rate, in fact exactly the opposite occurs. One highly improbable possibility is that the equipment or the doctors in the private sector is less efficient than in the public one. Another alternative considers instead psychologic reasons: doctors would be reluctant to notify paying patients that they have been exposed to the disease or maybe are even asymptomatic disease-carriers. At this point, independently of the reason, the question arises whether this diagnosis failure rate contributes in a substantial way to the endemicity of the disease.

On the basis of the above discussion, we would like then to address two main questions: namely whether the somewhat surprising difference in the diagnosis rates of the public and private sector makes a relevant difference for the disease eradication, and whether this eradication is at all possible.

---

3. The model assumptions

The flow diagram in Figure 1 below corresponds to the understanding among epidemiologists in India of how a susceptible person in a population may become infected and infectious, move through treatment, recover and then possibly become infectious again. This provides the road map for writing the mathematical model for TB. In India, the healthcare sector is segregated into two sectors: the public one (i.e. government run hospitals) and the private one (clinics and/or hospitals run by private practitioners). Research shows that diagnosis and treatment in the public sector has a larger success rate in comparison to seeking care in the private sector. But, unfortunately, due to different reasons TB patients tend to seek care in private sector clinics/hospitals. TB spread is determined largely by the nature of interaction of patients with active TB with the rest of the population.

We assume that the total human population is $N(t)$. This population is divided into six classes consisting of susceptible $S(t)$, latently-infected $L_1(t)$, latent but recovered $L_2(t)$, the infectious or diseased population $D(t)$, population treated in public sector $T_1(t)$ and lastly the population treated in private sector $T_2(t)$. Note that the $L_2$ individuals are cured (recovered) after completing six-months of anti-TB treatment (ATT). These cured individuals remain latent as TB bacteria may remain dormant in the host body. Thus

$N(t) = S(t) + L_1(t) + L_2(t) + D(t) + T_1(t) + T_2(t).$

(1)

For convenience we omit the explicit dependence on time in these population classes.

The model, given by equations (2) below, assumes standard incidence between the diseased and susceptible population and mass action interaction between the latent but recovered $L_2$ individuals and diseased populations. A fraction $\sigma$ of the latent persons exposed to infection $L_1$ rapidly moves to the infectious/diseased state $D$, because the disease progresses spontaneously to the virulent form. The rest ($1- \sigma$) remain latently infected ($L_1$) for a long period of time. In the model we have considered only endogenous reactivation which means that the TB bacteria which remain dormant have now become active and thus the individual becomes infectious. Thus an $L_1$ individual becomes diseased only by reactivation. Also, it is known that the relapse to the infectious state ($D$) from the latent state ($L_2$) is enabled by contact with an infectious individual. Thus recovered individuals, $L_2$, become diseased by re-infection. This can only happen after interacting with the diseased population, $D$. Therefore, we have considered mass-action of $L_2$ and $D$ individuals to model these interactions and thus to represent the re-infection process of the $L_2$ class.

Further, note that the individuals under treatment cannot infect the suceptibles and recovered populations, since the TB bacterium becomes inactive already in the first few weeks of the treatment. Also, other safety measures, like wearing masks, reduce the chances of transmitting disease early in the treatment. Therefore, in the model we have not assumed transmission of disease from individuals undergoing treatment i.e. $T_1$ and $T_2$.

The parameter $\beta$ is reported by epidemiologists as the number of people infected by one infected person per year. Usually people who are in class $D$ show symptoms and seek treatment. The rate at which treatment is sought from public clinics is $\nu_1$ and from private clinics is $\nu_2$. Not all diseased individuals are successfully treated. Those seeking treatment from the public and the private sectors remain diseased at rates $\mu_{12}$ and $\mu_{22}$ respectively. This happens for several reasons including noncompliance of drug regimen, drug resistance or other health complications. Recovery after treatment occurs at rates $\mu_{11}$ and $\mu_{21}$ from the public and the private sector, respectively. Recovered individuals are considered latently infected as the TB bacilli remain in their lungs. So these fall into the class $L_2$ and are more likely to relapse into the diseased ($D$) class than those in class $L_1$. The progression to diseased state from the $L_1$ class occurs at rate $\phi_1$. The progression rate to the diseased state from the $L_2$ class is $\phi_2$ with $\phi_2>\phi_1$. In India 40% of the population lies in classes $L_1$ and $L_2$, with 1% of the population being added every year.

---

4. The mathematical model

Recalling (1) and using standard incidence for modeling the disease spread among susceptibles and diseased, as explained in the previous section, we have

$\begin {eqnarray}\label{eq:Model} \frac{dS}{dt}&=& A-\beta \frac{SD}{N}-\alpha_0 S \\ \nonumber \frac{dL_1}{dt}&=&(1-\sigma)\beta \frac{SD}{N} - \alpha_0 L_1 - \phi_1 L_1 \\ \nonumber \frac{dD}{dt}&=& \sigma \beta \frac{SD}{N} -\alpha_2 D + \phi_1 L_1 +\phi_2 L_2 D + (\mu_{12} T_1+\mu_{22} T_2) -(\nu_1 + \nu_2) D \\ \nonumber \frac{dT_1}{dt}&=& \nu_1 D - (\alpha_3 +\mu_{11} +\mu_{12})T_1 \\ \nonumber \frac{dT_2}{dt}&=& \nu_2 D - (\alpha_4 + \mu_{21} +\mu_{22})T_2 \\ \nonumber \frac{dL_2}{dt}&=&\mu_{11} T_1 + \mu_{21} T_2 -\alpha_1 L_2 - \phi_2 L_2 D \nonumber\end{eqnarray}$

(2)

For later purposes, we give here the Jacobian of (2)

$J= \left( \begin{array}{cccccc} J_{11}&\beta \dfrac {SD}{N^2}&\beta \dfrac SN \left(\dfrac DN-1\right) & \beta \dfrac {SD}{N^2}&\beta \dfrac {SD}{N^2}&\beta \dfrac {SD}{N^2} \\ J_{21}&J_{22}&J_{23} &J_{24} & J_{25} &J_{26} \\ J_{31}&J_{32}&J_{33}&J_{34} & J_{35}& J_{36}\\ 0&0&\nu_1&J_{44}&0&0 \\ 0&0&\nu_2&0&J_{55}&0 \\ 0&0&-\phi_2 D&\mu_{11}&\mu_{21}&J_{66} \end{array} \right)$

(3)

with

$\begin{eqnarray*} J_{11}=\beta \dfrac DN \left(\dfrac SN-1\right)-\alpha_0, J_{21}=(1-\sigma) \beta \dfrac DN \left(1-\dfrac SN\right), \\ J_{22}=-(1-\sigma) \beta \dfrac {SD}{N^2} -\left(\alpha_0+\phi_1\right), J_{23}=J_{24}=J_{25} =J_{26}=-(1-\sigma) \beta \dfrac {SD}{N^2}, \\ J_{31}=\sigma \beta \dfrac DN \left(1- \dfrac SN\right), J_{32}=\phi_1-\sigma \beta \dfrac {SD}{N^2},\\ J_{33}=\sigma \beta \dfrac SN \left(1- \dfrac DN\right)+\phi_2 L_2 -(\nu_1+\nu_2+\alpha_2), J_{34}=\mu_{12}-\sigma \beta \dfrac {SD}{N^2},\\ J_{35}=\mu_{22}-\sigma \beta \dfrac {SD}{N^2}, J_{36}=\phi_2 D-\sigma \beta \dfrac {SD}{N^2}, J_{44}=-(\alpha_3+\mu_{11}+\mu_{12}),\\ J_{55}= -(\alpha_4+\mu_{21}+\mu_{22}), J_{66}= -(\alpha_1+\phi_2 D).\end{eqnarray*}$

---

5. Basic reproduction number

It is easily seen that there are only two possible equilibria, as all other combinations of population values lead to some inconsistency in the solution of the equilibrium system.

The disease-free equilibrium (DFE) $E_0=(S_0,0,0,0,0,0)$ is easily assessed, the value of susceptibles coming from the first equilibrium equation of (2), $S_0=A \alpha_0^{-1}$, with all other populations vanishing. Then there is possibly the coexistence equilibrium $E^*=(S^*,L_1^*,D^*,T_1^*,T_2^*,L_2^*)$, which will be discussed later.

Note that presently there are no available data on the value of $R_0$ for TB in India. We thus calculate the basic reproduction number $R_0$ as the spectral radius of next generation matrix [30].

There are two infected classes $L_1$ and $D$, so evaluating the gains and losses of each such compartment, we have:

$\left[ \begin{array}{c} \textrm{New infections, i.e. Gains to } L_1\\ \textrm{New infections, i.e. Gains to } D \\ \textrm{Losses from } L_1\\ \textrm{Losses from } D\\ \end{array} \right] = \left[ \begin{array}{c} (1-\sigma)\beta \displaystyle \frac{SD}{N}\\ \sigma\beta \displaystyle \frac{SD}{N} + \phi_2 L_2D + \phi_1L_1 \\ (\alpha_0+\phi_1)L_1 \\ (\alpha_2 + \nu_1 + \nu_2)D - (\mu_{12} T_1 + \mu_{22} T_2) \end{array} \right]$

By suitably taking partial derivatives, we find

$\begin{align} &\ \ \ \ \ \ \ \ \ \ \ F=\left[ \begin{matrix} \frac{\partial }{\partial {{L}_{1}}}\left[ (1-\sigma )\beta \frac{SD}{N} \right]&\frac{\partial }{\partial D}\left[ (1-\sigma )\beta \frac{SD}{N} \right] \\ \frac{\partial }{\partial {{L}_{1}}}\left( \sigma \beta \frac{SD}{N}+{{\phi }_{2}}{{L}_{2}}D \right)&\frac{\partial }{\partial D}\left( \sigma \beta \frac{SD}{N}+{{\phi }_{2}}{{L}_{2}}D \right) \\ \end{matrix} \right] \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\left[ \begin{matrix} -(1-\sigma )\beta \frac{SD}{{{N}^{2}}}&(1-\sigma )\beta \frac{S}{N}-(1-\sigma )\beta \frac{SD}{{{N}^{2}}} \\ {{\phi }_{1}}-\sigma \beta \frac{SD}{{{N}^{2}}}&\sigma \beta \frac{S}{N}-\sigma \beta \frac{SD}{{{N}^{2}}}+{{\phi }_{2}}{{L}_{2}} \\ \end{matrix} \right], \\ &V=\left[ \begin{matrix} \frac{\partial }{\partial {{L}_{1}}}(({{\alpha }_{0}}+{{\phi }_{1}}){{L}_{1}})&\frac{\partial }{\partial D}(({{\alpha }_{0}}+{{\phi }_{1}}){{L}_{1}}) \\ {{V}_{21}}&{{V}_{22}} \\ \end{matrix} \right]=\left[ \begin{matrix} {{\alpha }_{0}}+{{\alpha }_{1}}&0 \\ 0&{{\alpha }_{2}}+{{\nu }_{1}}+{{\nu }_{2}} \\ \end{matrix} \right], \\ \end{align}$

with

$\begin {eqnarray*} V_{21}= \dfrac{\partial}{\partial L_1}( (\alpha_2 + \nu_1 + \nu_2)D - (\mu_{12} T_1 + \mu_{22} T_2))\\ V_{22}= \dfrac{\partial}{\partial D}( (\alpha_2 + \nu_1 + \nu_2)D - (\mu_{12} T_1 + \mu_{22} T_2)). \end{eqnarray*}$

Evaluation at $E_0$ then gives

$F_{E_0} = \left[ \begin{array}{cc} 0&(1-\sigma) \beta \\ \phi_1 &\sigma \beta \end{array} \right], V_{E_0} = \left[ \begin{array}{cc} \alpha_0+\alpha_1 &0 \\ 0&\alpha_2+\nu_1 + \nu_2 \end{array} \right]$

so that

$\begin{eqnarray*} V_{E_0}^{-1} = \left[ \begin{array}{cc} \displaystyle\frac{1}{\alpha_0 + \phi_1} & 0 \\ 0 &\displaystyle \frac{1}{\alpha_2 + \nu_1 + \nu_2} \end{array} \right], \\ G = F_{E_0} V_{E_0}^{-1} = \left[ \begin{array}{cc} 0 &\displaystyle\frac{(1-\sigma)\beta}{\alpha_2 + \nu_1 + \nu_2} \\ \ \\ \displaystyle\frac{\phi_1}{\alpha_0 + \phi_1} & \displaystyle\frac{\sigma\beta }{\alpha_2 + \nu_1 + \nu_2} \\ \end{array} \right].\end{eqnarray*}$

It follows that the dominant eigenvalue of $G$ gives the value of the basic reproduction number

$R_0 = \frac 12 \left\{ \displaystyle\frac{\sigma\beta }{\alpha_2 + \nu_1 + \nu_2} + \left[ \displaystyle\frac{\sigma^2 \beta^2 }{(\alpha_2 + \nu_1 + \nu_2)^2} + 4 \displaystyle\frac{(1-\sigma) \beta \phi_1 }{(\alpha_0 + \phi_1)(\alpha_2 + \nu_1 + \nu_2)} \right]^{\frac 12} \right\}.$

With the values reported in Table 1, it appears that $R_0\approx 1.8414$, a value that indicates that the disease is endemic, and this value of the basic reproduction number is also not too close to the eradication threshold.

To gain more insight, we provide a sensitivity analysis of $R_0$ in terms of the most important parameters of the model. Since the basic reproduction number depends on six them, we choose to fix $\beta$ as a reference. In each parameter space given by $\beta$ coupled with one of the remaining parameters, namely $\sigma$, $\phi_1$, $\alpha_0$, $\alpha_2$, $\nu_1$ and above all $\nu_2$, we plot the surface $R_0$ and the contour line corresponding to the threshold $R_0=1$. The resulting surfaces and contour curves are plotted in Figures 2-7.

Recall that to achieve disease eradication, a value of $R_0<1$ must necessarily be attained. We thus concentrate now on possible means to achieve this goal. The following discussion aims at ascertaining whether acting on each parameter that appears in the expression for $R_0$, a sufficient reduction of this threshold can at all be obtained. This could perhaps lead to indications for suitable policies to be pursued by the authorities in order to curb the burden of this endemic and pernicious disease.

From Figure 2 it appears that from the current situation, variations in $\sigma$ will not lead to any improvement, in fact if this parameter is enlarged, the endemicity of the disease will be more pronounced, while even if it is reduced to values very close to zero, the disease is not going to disappear, without a simultaneous reduction of the transmission rate $\beta$. A similar situation appears from Figure 3. The system behavior is sensitive to the parameter $\phi_1$, in that a relatively small reduction of the latter induces disease eradication. But this parameter is an intrinsic parameter of the disease, the spontaneous progression to the active form from an unknown latently infected individual. It therefore appears that acting on this parameter in practice is extremely difficult if at all possible.

Sensitivity with respect to the diagnosis and treatment rates is shown in Figures 4 and 5. The same result is obtained in Figure 7, by an increased diseased-related mortality, which however goes in the opposite direction of the goal of saving lives. Incidentally, this is an example of what in ecology and epidemiology is sometimes found, that something harmful at the individual level is beneficial for the community, and vice versa. Such an observation is found for instance in beehives, [2]. It appears that an increase in the identification of the cases and their cure will help in reducing the endemicity range of the disease. Clearly Figure 6 shows that an increase in the natural death rate will contribute also to disease eradication, but this is certainly not a practical recommendation to follow.

In Figure 8 we observe that changes in $\phi_2$ are irrelevant on the value of the basic reproduction number. This however is to be expected, in that $\phi_2$ does not appear in the definition of $R_0$.

Finally, in Figure 9 we plot the value of the basic reproduction number as function of the treatment rates in both public and private sectors. The findings indicate that the difference in the public and private sectors diagnosis rate is not essential for eradicating the disease. In fact, even if they were 100%, the plot indicates that the surface for $R_0$ would still be above the critical threshold 1, so that the disease could not be eradicated.

---

6. Nonlinear stability of the disease-free equilibrium

We can assess the asymptotic stability of the disease-free equilibrium $E_0$ by means of a suitable nonlinear Lyapunov function. As candidate, choose:

${\mathcal{L}}_1= a S_0 \left[ \frac {S-S_0}{S_0} - \ln \left( 1+ \frac {S-S_0}{S_0} \right) \right] + bL_1 + c D + eT_1 + fT_2 + gL_2.$

Note that ${\mathcal{L}}_1$ is nonnegative and ${\mathcal{L}}_1(E_0)=0$. Differentiation with respect to time and use of (2), with $N_0=S_0$, leads to

$\begin{eqnarray*} {\mathcal{L}}_1'=-aA \frac {(S-S_0)^2}{SS_0} + \beta \frac {(S-S_0)D}{N} [b(1-\sigma) -a + c \sigma]\\ + \beta \frac DN (N-N_0) S_0 [b(1-\sigma) + c \sigma] + L_1 [c\phi_1 -b(\alpha_0 + \phi_1)] + D L_2 (c-g)\phi_2\\ + D \left[ \beta \frac {S_0}{N_0} [b(1-\sigma) + c\sigma ] + e \nu_1 + f \nu_2 - c (\alpha_2+\nu_1+\nu_2)\right]- \alpha_1 g L_2 \\ + T_1 [c\mu_{12} - e(\alpha_3+\mu_{11}+\mu_{12}) + g \mu_{11}] + T_2 [c\mu_{22} - f(\alpha_4+\mu_{21}+\mu_{12}) + g \mu_{21}].\end{eqnarray*}$

Taking then $c=g$ and $a=b(1-\sigma) + c \sigma$, since $N-N_0 \le N$, we find the upper bound

$\ \ \ \ \ \ \ \ {\mathcal{L}}_1' \le -aA \frac {(S-S_0)^2}{SS_0} + L_1 [c\phi_1 -b(\alpha_0 + \phi_1)]\\ \ \ \ \ \ \ \ \ + D \left[ a \beta S_0 + a \beta \frac {S_0}{N_0} + e \nu_1 + f \nu_2 - c (\alpha_2+\nu_1+\nu_2)\right] - \alpha_1 c L_2 \\ + T_1 [c\mu_{12} - e(\alpha_3+\mu_{11}+\mu_{12}) + c \mu_{11}] + T_2 [c\mu_{22} - f(\alpha_4+\mu_{21}+\mu_{12}) + c \mu_{21}].$

Finally imposing the following inequalities

$\begin{eqnarray*} c\phi_1 < b(\alpha_0 + \phi_1), a \beta S_0 \left( 1+\frac 1{N_0}\right) + e \nu_1 + f \nu_2 < c (\alpha_2+\nu_1+\nu_2), \\ c(\mu_{12}+\mu_{11}) < e(\alpha_3+\mu_{11}+\mu_{12}), c(\mu_{22}+\mu_{21}) < f(\alpha_4+\mu_{21}+\mu_{22}),\end{eqnarray*}$

it follows that

${\mathcal{L}}_1'\le -a A \frac {(S-S_0)^2}{SS_0} - \alpha_1 c L_2 \le 0$

showing that it is nonpositive definite. This is not enough to ensure global stability. Rather, it could be the starting point to assess the largest domain of attraction of the DFE. Following [20], one can search for the largest invariant subset in $\mathbf R^6$ for which ${\mathcal{L}}_1'=0$, i.e. $S=\widetilde S$ and $L_2=0$. In principle it can be obtained by making such substitutions into the right hand side of (2) and solving the resulting system. This appears to be rather complicated, so that we will only investigate the situation via numerical methods below.

---

7. Local stability of the disease-free equilibrium and backward bifurcation

At the disease-free equilibrium the Jacobian (3) has immediately three easy eigenvalues $-\alpha_0<0$, $-(\alpha_0+\phi_1)<0$, $-\alpha_1<0$. The remaining ones, that characterize the stability of this equilibrium, are those of the submatrix

$J_0=J\left( E_0\right)= \left( \begin{array}{ccc} \sigma \beta -(\nu_1+\nu_2+\alpha_2)&\mu_{12}&\mu_{22} \\ \nu_1&-(\alpha_3+\mu_{11}+\mu_{12})&0 \\ \nu_2&0&-(\alpha_4+\mu_{21}+\mu_{22}) \\ \end{array} \right)$

(4)

We now investigate numerically the disease-free equilibrium and the endemic one. Departing from the Table values only for the new individuals recruitment rate, setting it to the value $A=100$, for which the susceptibles at disease-free equilibrium would attain the value $S_0=14085.50704$ we find that $R_0=1.8414>1$, so that the disease is endemic. However, the computation of the eigenvalues shows that they are negative, $-1.5390$, $-0.7265$, $-1.0528$ and this is further convalidated by the Routh-Hurwitz conditions at the DFE. Denoting by $M_2$ the sum of the principal minors of order 2 of $J_0$, they are

$\begin{eqnarray*} -{\textrm{tr}}\left( J_0\right) = 3.3184>0, \ \ \ -\det \left( J_0\right) = 1.1772>0,\\ -{\textrm{tr}}\left( J_0\right) M_2 \left( J_0\right) -\det \left( J_0\right) = 10.4481>0.\end{eqnarray*}$

This shows that the DFE is stable, therefore giving rise to a bistability situation and to a backward, or subcritical, bifurcation, see [29] and [12] page 28. We investigated numerically the possibility of attaining the DFE, searching the initial values space for values that upon integration would lead to the DFE. In fact, this occurs, but for unrealistically small values of the diseased classes, namely $L_1(0)=D(0)=10^{-15}$, but the value for the susceptibles can be large, reasonably far away from the equilibrium value, $S(0)=S_0+9950$, while the remaining initial conditions have always been chosen as

$T_1(0) = 0.49D(0), \ \ \ \ T_2(0) = 0.41 D(0), \ \ \ \ L_2(0) = 0.85 T_1(0) + 0.51 T_2(0).$

Any increase in the values of $L_1(0)$, $D(0)$ or $S(0)$ given above would instead make the system tend to the endemic equilibrium.

We then tried to study this situation as $\beta$ decreases. For $\beta=2.3$, we find $R_0=1.2044$ but the eigenvalues at DFE are still negative, $-1.5623$, $-0.7476$, $-1.0535$, so that this equilibrium is stable. The initial conditions leading to the DFE can now be chosen a bit larger, namely $L_1(0)=D(0)=10^{-12}$ while the choice for the susceptibles is now much wider, up to $S(0)=S_0+179000$.

A further decrease, $\beta=1.5$, gives now $R_0=0.97045$, for which apparently the disease would be eradicated. But again, in view of the backward bifurcation, the endemic equilibrium persists and is attained for the following choice of initial conditions

$\begin{eqnarray*} S(0)=S_0+180000=194084.5070, \ \ \ D(0)=0.3 S(0)=58225.3521, \\ \nonumber T_1(0) = 28530.4225,\ \ \ T_2(0) = 23872.3944, \\ \nonumber L_2(0) = 36425.7803, \ \ \ L_1(0) = .8 S(0)=155267.6056.\end{eqnarray*}$

(5)

It is plotted in Figure 10.

A similar situation arises for $\beta=0.9$, giving $R_0=0.7501$, but the endemic equilibrium is found for the same initial conditions (5), see Figure 11. Here susceptibles more than double the value of Fig 10, latently infected do not change sensibly, infected are found at a smaller steady level.

For $\beta=0.6$, we find $R_0=0.6116$, in this case using the initial conditions (5), the endemic equilibrium seems now to have disappeared, and the same holds for $\beta=0.4$, for which $R_0=0.4988$.

---

8. Coexistence equilibrium analysis

Although the result of the previous section indicates that the disease could in principle be eradicated, in addition, we could now pursue an alternative road for trying to curb it.

First of all, we investigate whether the the coexistence equilibrium can be assessed analytically. We solve the fourth and the fifth equilibrium equations in terms of $D$ and solve the sixth one for $L_2$, to get $T_1= D \eta_1$, $T_2= D \eta_2$ with

$\begin {eqnarray}\label{eq:Ts} \eta_1= \frac {\nu_1 }{\alpha_3 + \mu_{11} + \mu_{12}}, \ \ \ \eta_2= \frac {\nu_2 }{\alpha_4 + \mu_{21} + \mu_{22}}, \ \ \ L_2= \frac {\mu_{11} T_1 + \mu_{21} T_2 }{\alpha_1 + \phi_2 D}. \end{eqnarray}$

(6)

Using the first two above equations in the last one, we find $L_2$ in terms of $D$,

$\begin {eqnarray}\label{eq:L2} L_2= \frac D{\alpha_1 + \phi_2 D} \left[ {\frac {\mu_{11} \nu_1 }{\alpha_3 + \mu_{11} + \mu_{12}} + \frac {\mu_{21} \nu_2 }{\alpha_4 + \mu_{21} + \mu_{22}} }\right] = \frac {D(\mu_{11} \eta_1 + \mu_{21} \eta_2)}{\alpha_1 + \phi_2 D} . \end{eqnarray}$

(7)

Taking the linear combination of the first two equilibrium equations with weights $1-\sigma$ and $1$ respectively, we find

$\begin {eqnarray}\label{eq:L1} L_1= (1-\sigma) \frac {A-\alpha_0 S}{\alpha_0 + \phi_1}. \end{eqnarray}$

(8)

Adding the first, with weight $\sigma$, and the third equilibrium equations, setting $\Omega=\mu_{12}\eta_1+\mu_{22}\eta_2-\alpha_2-\nu_1-\nu_2 \in \mathbf R$ and using (7), we obtain

$\begin {eqnarray}\label{eq:10} \sigma A -\alpha_0 \sigma S + \phi_1 L_1 + \phi_2 D^2 \frac {\mu_{11}\eta_1+\mu_{21}\eta_2}{\alpha_1 + \phi_2 D} + D \Omega=0. \end{eqnarray}$

(9)

Rearranging (9) and letting

$W=\Omega +\mu_{11}\eta_1+\mu_{21}\eta_2\in \mathbf R, \ \ \ \theta=\frac {\sigma \alpha_0+\phi_1}{\alpha_0+\phi_1}>0,$

leads to

$\begin {eqnarray}\label{eq:11} \sigma A \alpha_1 + \phi_1 \alpha_1 L_1 + D (\sigma A \phi_2 + \Omega \alpha_1) + \phi_1 \phi_2 L_1 D + \phi_2 D^2 W= \alpha_0 \sigma S (\alpha_1 + \phi_2 D). \end{eqnarray}$

(10)

Use of (8) into (10) gives

$\begin {eqnarray}\label{eq:16} \Phi(S,D):= \theta A \alpha_1 -\alpha_0 \alpha_1 \theta S + D (\theta A \phi_2 + \Omega \alpha_1) - \alpha_0 \phi_2 \theta SD + \phi_2 D^2 W=0. \end{eqnarray}$

(11)

Finally, the first equilibrium equation can be rewritten as $(A -\alpha_0 S) [S+L_1+D+T_1+T_2+L_2] - \beta SD =0$ from which, substituting into it (6), (8) and (7), we obtain

$\begin {eqnarray}\label{eq:first} &\Psi(S,D):=- \beta SD \\ \nonumber &+(A -\alpha_0 S) \left[ S+(1-\sigma) \dfrac {A-\alpha_0 S}{\alpha_0+\phi_1} +D(\eta_1+\eta_2+1) + D \dfrac {\mu_{11}\eta_1+\mu_{21}\eta_2}{\alpha_1+\phi_2 D}\right] =0. \end{eqnarray}$

(12)

The curve obtained by taking the common denominator in (12) and setting the numerator to zero is a third order implicit function and therefore very difficult to study, even numerically. The values of the diseased and susceptible populations at the coexistence equilibrium would be obtained by the intersection of $\Phi$ and $\Psi$ in the first quadrant of the $S-D$ plane, while the remaining populations would come from (6), (7) and (8).

The mathematical problem appears to be a hard task. For this reason, in order to gain anyway some insight into the actual situation, we make a very strong simplifying assumption. The mathematical difficulty arises from the denominator in the last fraction of (12). To simplify it, we assume $\phi_2=0$, which implies that there are no relapses after treatment. As said, this is very unlikely, but we use it as a probe into the problem.

In view of the simplification, the equations for $\Phi$ and $\Psi$ become easier to handle, in fact the first one is a straight line and the second one a conic section, namely:

$\Phi_s(S,D):= \theta A -\alpha_0 \theta S + D \Omega =0,$

(13)

$\label{eq:psi_s} \Psi_s(S,D):= \rho A^2 + \nu A S + \pi A D - (\alpha_0 \pi + \beta) SD - \alpha_0 \theta S^2=0,$

(14)

where

$\rho = \frac {1-\sigma}{\alpha_0 + \phi_1}>0 , \ \ \ \nu = \frac {2\alpha_0 \sigma- \alpha_0 + \phi_1}{\alpha_0 + \phi_1} \in {\mathbf{R}}, \ \ \ \pi = \eta_1+\eta_2 + 1 + \frac {\mu_{11}\eta_1+\mu_{21}\eta_2}{\alpha_1}>0.$

Note that the straight line $\Phi_s(S,D)$ in the $S-D$ plane has slope and height at the origin of uncertain signs, $\alpha_0 \theta \Omega ^{-1} \in \mathbf R$, $A \theta \Omega^{-1} \in \mathbf R$. But its intercept with the $D$ axis is positive, $S_0=A \alpha_0^{-1}>0$. Now $\Psi_s(S,D)$ is a nondegenerate curve if its first invariant does not vanish, a condition that in fact we assume:

$\frac 14 A^2 \left[ \alpha_0 \theta \pi^2 -(\alpha_0 \pi + \beta) \pi \nu - (\alpha_0 \pi + \beta)^2 \rho \right] \ne 0.$

In particular it is a hyperbola, since its second invariant is negative, $- 4^{-1} (\alpha_0 \pi + \beta)^2<0$. We study it by intersecting it with vertical lines, $S=k$. Its intersections with these lines are the points

$D= \frac {\alpha_0 \theta k^2- k \nu A -\rho A^2}{A \pi - (\alpha_0 \pi + \beta) k}.$

The latter are feasible when positive, which shows that there is a feasible branch of $\Psi_s(S,D)$ in between the values $S_*$ and $\widehat S=\widehat k=A\pi (\alpha_0 \pi + \beta)^{-1}$. Here $S_*=S_{\pm}=k_{\pm}=A(2\alpha_0 \theta)^{-1} [\nu \pm (\nu^2 + 4 \alpha_0 \rho \theta)^{\frac 12}]$, where the plus sign is taken when $\nu>0$ since then $S_*=S_+>0>S_-$ and the minus sign whenever $\nu<0$, since in such case $S_*=S_->0>S_+$. This situation leads to two possible configurations. Thus, depending on whether $S_*>\widehat S$ or $\widehat S> S_*$, $\Psi_s(S,D)$ has a branch raising up to $+\infty$ at $S=\widehat S$ from the zero at $S=S_*$, in the former case, or conversely decaying from $+\infty$ at $S=\widehat S$ to zero at $S=S_*$ in the second one. An intersection with $\Phi_s(S,D)$ cannot occur for sure if the zero $S_0$ of this straight line, whenever its slope is positive, lies beyond the zero of $\Psi_s(S,D)$, or conversely if the slope is negative, namely either for $S_0> \max \{ S_*, \widehat S\}$, $\Omega >0$ or for $S_0 < \min \{ S_*, \widehat S\}$, $\Omega<0$. In our case, taking again $A=30$, we find $\Omega = -1.0754<0$, $S_*=4225.352112676057$ $\widehat S=213.0447$, $S_0=4225.352112676056<S_*$. We would need $S_0<213.0447$ which is not true, indicating that the intersection exists. In fact, in Figure 12 we plot the situation and discover that in fact there are 2 intersections. In order to eradicate the disease, one could try to make them vanish, through a saddle-node bifurcation, by influencing the slope of the straight line $m=\alpha_0 \theta \Omega^{-1} <0$, but that cannot certainly be achieved by reducing the natural mortality rate $\alpha_0$. An increase of $m$ is then necessary so that the intersections do not occur. The same result can be achieved by lowering the height at the origin, $A\theta (-\Omega)^{-1}>0$. A decrease in the recruitment rate $A$ could be viable, but highly improbable. One can try then to reduce $\theta$ or to increase $-\Omega$. For the former, one could try to reduce $\sigma$, but that at most gives $\theta|_{\sigma=0}=\phi_1(\alpha_0+\phi_1)^{-2}$, which might not be enough; on the other hand, observe that

$\dfrac {d \theta}{d \phi_1} = \alpha_0 \dfrac {1-\sigma}{(\alpha_0+\phi_1)^2}>0,$

implying that we must reduce $\phi_1$. The second alternative it to increase $-\Omega=\alpha_2+\nu_1+\nu_2-\mu_{12}\eta_1-\mu_{22}\eta_2$, the maximum being of course $-\Omega=2+\alpha_2$, if the diagnosis rates achieve 100% of precision and the failure rates after treatment vanish altogether.

Disease eradication can in fact be achieved for an almost extreme case, taking $\phi_1=0.001$, $\phi_2=0$, $\sigma=0$, $\mu_{12}=0$, $\mu_{22}=0$, $\nu_1=1$, $\nu_2=1$ and in this case $\Omega = -2.320$, see left column of Figure 13. In the right column, a similar situation occurs, for the following parameters that do not achieve the extreme values, but are really quite close to them: $\phi_1=0.001$, $\phi_2=0$, $\sigma=0.01$, $\mu_{12}=0.001$, $\mu_{22}=0.001$, $\nu_1=.999$, $\nu_2=.999$ for which $\Omega = -2.316$. In any case these considerations would hold only for $\phi_2=0$, which is a very restrictive situation.

---

9. Discussion

The model that has been introduced here was meant to compare the TB treatments performed in the public and private sectors of health care in India. The bottom line of the results of our investigation are the considerations that can be inferred from Fig. 9, namely that the low rate of diagnosis actually found in the private sector, whether it be due to malpractice, poorer diagnostic means or simply ascribed to the fact that doctors may be reluctant to let paying patients know that they are infected, seems not to constitute the main problem in the disease endemicity. It appears thus that even achieving $100\%$ correct diagnosis in both sectors would not help in eradicating the disease.

From the extensive simulations that we ran, see Figures 2-8, it appears that apart from the disease contact rate, the effect of the other parameters affecting $R_0$ is rather thin and essentially negligible in the disease eradication issue. At most, changing some of them might render the set in the parameter space where $R_0<1$ smaller, but without a corresponding substantial change in the disease contact rate $\beta$ the critical threshold cannot be crossed. On the other hand, by drastically reducing this parameter, the basic reproduction number can be reduced to values less than one. In this case, however, another unexpected and unpleasant phenomenon occurs, the onset of a backward, or subcritical, bifurcation, for which, in spite of the theoretical result assessing the possibility of achieving disease eradication, because the system is coming from an endemic equilibrium, it will continue on this manifold, even if the disease-free equilibrium exists. This happens until really low values of $\beta$ are achieved, which might be unrealistic in practice.

Alternatively, for values of $\beta$ that bring $R_0$ below the critical threshold, or even larger values of $\beta$, one could try to achieve the disease-free equilibrium by acting on the initial conditions of the system, drastically and suddenly reducing them, so that this change brings the system into the domain of attraction of the disease-free equilibrium. But our experiments indicate that in order to fall into the domain of attraction of this equilibrium, the number of infected must be so small to be essentially unreacheable in ordinary life.

A different approach has also been attempted, namely to try to render the endemic equilibrium unfeasible. This approach appears to be analytically untractable, except for the unrealistic case of no disease relapses, i.e. $\phi_2=0$, from the class $L_2$ of latent and cured individuals, and also difficult to address numerically. In principle however this approach indicates an alternative way of fighting the disease, namely to render the coexistence equilibrium unfeasible.

For the particular case $\phi_2=0$, a highly improbable situation to achieve, since in practice it is nearly impossible to prevent relapses, the system is reduced to the intersections of the curves (11) and (12), i.e. $\Phi_s$ and $\Psi_s$. The point at which they meet provides the population values for susceptibles and diseased at the endemic equilibrium. To attain this situation some means of rendering disease relapses after cure impossible or negligible enough should be devised, to keep $\phi_2$ at zero or at a very low level. Although we did not perform on this an exhaustive investigation, nevertheless some information can be gathered. Our simulations in the previous section indicate that disease eradication is possible but for values of the some of the parameters that are almost extreme. This entails for instance that $\phi_1$ must be reduced. However, this is an almost impossible task, as this parameter models the intrinsic progress from the latently infected individuals to active disease outbreak. In other words, even if a drug would be discovered to slow down this progression, it would need to be administered to a set of unknown individuals, the latently infected, or it could be given to the same set of people after their identification, following a suitable screening, which is probably a measure that has enourmous economical costs. The other relevant parameters are the disease diagnosis rates $\nu_i$, in both public and private sectors. These should achieve almost 100% certainty. Further, $\mu_{1i}$, the relapse rates after treatment, should drop almost to zero and the proportion of the primary latently infected progressing to the active disease stage should be lowered too, but not that dramatically. However this last task is rather difficult, as $\sigma$ is an intrinsic disease parameter and it is difficult to act on it, even if some drug to this effect were discovered, as this parameter pertains also to the unknown population of the latently infected as mentioned above. Hence, also this approach unfortunately does not appear to be viable in practice.

The rater sad conclusion that we must draw from all these considerations is therefore that eradicating the disease in the present state of affairs is rather difficult if not at all impossible.

---

Acknowledgments

This work was undertaken when the first author visited the University of Torino, with a WWS2 grant, which is thankfully acknowledged. The research has also been partially supported by the project "Metodi numerici nelle scienze applicate" of the Dipartimento di Matematica "Giuseppe Peano". EV gratefully acknowledges very useful discussions with Antoine Perasso and Rafael Bravo de la Parra and the referees for their constructive comments.

---