A synthesized model of tuberculosis transmission featuring treatment abandonment

1.   Introduction

Pulmonary tuberculosis (TB) is a public health problem in many countries worldwide and constitutes one of the leading causes of morbidity and mortality. This infectious disease puts a large number of people with weakened immune systems at risk, for example, due to the HIV/AIDS or diabetes comorbidities, chemo/radio-therapies, transplants, or malnutrition. Also, tuberculosis is a social disease that mainly affects poor people living in crowded and precarious conditions of big cities or penitentiary institutions, especially in low- and middle-income countries ^{[1, 2]}.

Tuberculosis is caused by a species of pathogenic bacteria called mycobacterium tuberculosis or Koch's bacilli. The bacteria can attack any part of the body, but usually attack the lungs. This disease is transmitted from person to person through microscopic airborne droplets from the cough or sneeze of a person with an active form of the disease. Thus, a susceptible person can become infected by inhaling airborne bacilli. Common symptoms of active pulmonary tuberculosis are severe cough (often with blood in the sputum), chest pain, weakness, fever, poor appetite, and weight loss.

However, people who are infected with TB do not always feel sick. The World Health Organization ^[3] estimates that a quarter of the world's population has tuberculosis in a "latent" or inactive form; that is, these people are carriers of inactive bacilli, they have no symptoms of the disease, and they cannot transmit the infection to others. However, inactive bacilli can regain their vitality together with replication capacity and cause an active form of the disease when the immune response of a latent person is disturbed or altered.

A noticeable epidemiological feature of tuberculosis is its long period of latency that may last from several weeks to a lifetime. On the one hand, people with a latent form of tuberculosis can last in this state for many years without actively developing the disease, as long as their immune system does not present alterations. On the other hand, latent people whose immune systems are weakened due to malnutrition or comorbidities (influenza, pneumonia, COVID19, kidney disease, diabetes, etc.) or those receiving immunosuppressive treatments for HIV, cancer, or transplants are at much higher risk of developing an active form of tuberculosis.

The active form of tuberculosis can be diagnosed with the sputum smear test (also known as acid fast bacillus testing—AFB+) that detects the presence of active bacilli, whereas the latent form of tuberculosis can be detected by the tuberculin skin test or the measurement of interferon-$\gamma$ in whole blood (also known as interferon-gamma release assays—IGRAs) ^[4]. Both forms of tuberculosis are treatable by antibiotics, although the treatment periods are rather long (6–12 months) and require prolonged medical monitoring ^[5].

Most patients may recover from a primary active TB infection without further manifestation of the disease if they strictly follow all the guidelines of medical treatment. Even though the drugs that are widely used nowadays for standard TB treatment (namely, isoniazid and rifampicin) are nominally characterized as "sterilizing", several studies question this characterization and manifest that a patient does not always achieve true sterilization after completing treatment ^{[6, 7]}. Besides, total bacilli extermination cannot be reliably confirmed by existing test methods, such as blood and sputum markers ^[8].

Moreover, low- and middle-income countries lack the infrastructure, funding, or capacity to diagnose and treat patients with the latent form of tuberculosis. Therefore, TB control in these countries is centered on preventing latent TB infections from developing into active disease and on diagnosing and treating patients with the active form of the disease ^{[9, 10]}.

Several studies (see e.g., ^{[11, 12]} and references therein) pointed out that treatment non-adherence and abandonment mostly occurs due to the patient-related factors (such as low income, unemployment, illiteracy), to the health services (reduced availability, intermittent access, and bureaucratic structure), and to the characteristics of the treatment itself (long duration and possible adverse effects). To prevent the spread of the contagious (or active) form of the disease it is important to ensure that all patients take their anti-tuberculosis medications exactly as prescribed and do not abandon the treatment before its full completion. Notably, the patients who do not follow the strict guidelines for treatment or those who abandon treatment before completion usually remain infectious and develop the so-called "drug-resistant" form of active tuberculosis. Moreover, people with active TB not receiving proper medical treatment may die from this infectious disease.

This paper aims to propose and analyze a lower-dimensional model of TB transmission that accounts for treatment abandonment. Among a wide variety of mathematical models describing the dynamics of tuberculosis transmission, there are only a few that are deliberately kept mathematically tractable and apt for theoretical analysis^*. The prime feature of such compartmental models is the relatively small number of state variables (four or fewer) that denote the disjoint classes compounding the target human population. In particular, there are several models where the human population is defined by the sum of four classes or compartments ^{[14, 15, 16, 17, 18]}: Susceptible and latently infected people (they are non-infectious), actively infected people (they are fully infectious), and people undergoing treatment (they are partially infectious). However, these models do not account for patients who abandon the treatment before its completion while the treatment abandonment enhances the disease spread.

^*A more detailed survey of simplest models describing the transmission of tuberculosis has been performed in ^[13].

On the other hand, more sophisticated models of higher dimension where all individuals who abandon treatment are gathered in an additional class or compartment have been proposed and studied during the last decade ^{[19, 20, 21, 22, 23]}. These and, possibly, other studies have disclosed through mathematical modeling the adverse effects of treatment abandonment on the disease persistence and spread.

It is worth mentioning that the four-dimensional model presented in ^[24] includes two classes of infectious people, those who receive treatment (partially infectious) and those who refuse ("lost to follow up" cases or fully infectious), along with susceptible and latent classes (non-infectious). This is an interesting formulation though it is based on a strong assumption that does not seem realistic. Namely, all individuals are timely identified as latently or actively infected, then they are immediately offered to start treatment, and then they either accept or refuse that offer. In other words, people with identified TB infections either start being treated or become "lost to follow up" cases without starting treatment. In contrast to the model presented in ^[24], the model proposed in this work assumes that all identified carriers of the active bacilli start being treated. Afterward, they either complete the treatment and become non-infectious or abandon the treatment before completion and become infectious.

The proposed model is presented in Section 2, while its basic reproductive number and its possible equilibria are derived in Section 3. Section 4 is centered on establishing the properties of global stability for two possible equilibria of the proposed model. The role of the treatment abandonment rate is then illustrated in Section 5 by examining two scenarios, with and without treatment abandonment. Finally, Section 6 provides conclusions of our work.

2.   Model for the TB transmission

In this section, we propose a four-dimensional model that incorporates a group of people with an active form of the disease, who start the anti-tuberculosis treatment and then quit without completing it. The distinctive feature of this model is that such people are not separated into an additional class. Thus, our model has only four state variables, namely:

● $S(t)$ – the number of susceptible individuals at the moment $t \geq 0$ who are bacillus-free.

● $L(t)$ – the number of latently infected individuals at the moment $t \geq 0$ who are carriers of inactive bacilli but they have no symptoms of the disease and cannot transmit it to the others (i.e., they are non-infectious).

● $I(t)$ – the number of actively infected individuals at the moment $t \geq 0$ who have not commenced the treatment yet; they are carriers of active bacilli, have symptoms of the disease, and are fully capable of transmitting it to the others (i.e., they are fully infectious).

● $T(t)$ – the number of individuals undergoing treatment at the moment $t \geq 0$; they are carriers of active bacilli but their infectiousness is reduced by the treatment (i.e., they are partially infectious).

In this model, we assume that all four classes are homogeneously mixed and the disease transmission depends on the frequency of respiratory contacts between the bacillus-free individuals and those carrying active bacilli.

Thus, the total population at each moment $t \geq 0$ denoted by

The interaction between four disjoint classes $S(t), L(t), I(t),$ and $T(t)$ is described by the following ODE system

with nonnegative initial conditions

It is assumed that all newborn individuals are bacillus-free, they belong to the $S$-class, and they are recruited at the constant rate $\Lambda > 0$ (see Eq (2.2a) above). Notably, even though the immunological changes during pregnancy may activate the TB infection in pregnant women with the latent form of the disease, there is no evidence of TB vertical transmission ^[25]. The average rate of natural mortality of all human individuals is denoted by $\mu > 0$.

The disease transmission is modeled by the mass action incidence through the term $\beta (I + \phi T)$ in Eqs (2.2a)–(2.2c). Here, $\beta > 0$ denotes the so-called per capita rate of effective respiratory contacts with infectious people (i.e., those leading to the contagion) and can be viewed as

When a susceptible person inhales active bacilli, he/she may either develop an active form of the disease with probability $p \in (0, 1)$ or progress to the class of latently infected people $L(t)$ with probability $(1-p)$ (cf. the first terms in Eqs (2.2b)–(2.2c)). The probability $p$ of developing an active form of the disease is related to the immune status of the human host that receives active bacilli. Thus, people with a weakened immune response due to comorbidities (influenza, pneumonia, COVID19, kidney disease, diabetes, etc.), malnutrition, or immunosuppressive treatments for HIV, cancer, or transplants are at much higher risk of developing an active form of tuberculosis.

Susceptible people may become infected through respiratory contacts with either fully infectious people (from the class $I(t)$) or partially infectious people (from the class $T(t)$). It has been shown that treatment reduces, but does not fully eliminate, the infectiousness of patients with active forms of tuberculosis ^[26]. Therefore, our model (2.2) assumes that people from the $T$-class are capable of transmitting the disease at a reduced rate $\phi \beta,$ where $\phi \in (0, 1)$ can be viewed as a coefficient of the infectivity reduction. The latter stays in line with other models describing the transmission of tuberculosis ^{[27, 28]}.

Non-infectious people from the class $L(t)$ (who are the carriers of inactive bacilli) may also develop the active form of the disease when their immune status is changed or altered. Such people progress to the class $I(t)$ at the rate $\epsilon > 0$ (see Eqs (2.2b)–(2.2c)). Fully infectious people from the class $I(t)$ are either recruited for treatment at the rate $\gamma > 0$ and then progress to the class $T(t)$ or they are removed at the disease-induced mortality rate $\alpha > 0$ (see Eqs (2.2c)–(2.2d)).

Thus, people commencing the treatment become partially infectious and belong to the class $T(t)$. They remain in this class before occurs one of the following events:

● they complete treatment and then return to the class $L(t)$ at the rate $\delta > 0$, meaning they may still undergo a TB reactivation or relapse towards an active form of the disease;

● they abandon treatment before completion, progress to the class $I(t)$ at the rate $k > 0$, and become fully infectious.

The latter is expressed by Eqs (2.2b)–(2.2d). It is worthwhile to emphasize that here we assume, on the grounds of arguments presented in ^{[6, 7, 8]}, the non-sterilizing treatment, which is typical for low- and middle-income countries^†. Finally, Figure 1 exhibits the flow diagram of the TB transmission model (2.2) and Table 1 recaps a detailed description of the model's parameters.

^†It should be recalled that TB control in low- and middle-income countries is mainly centered on diagnosing and treating active infections with pharm products bearing low sterilizing effects that are more accessible due to their moderate costs ^{[9, 10]}.

To showcase the well-posedness of the model (2.2), we formulate the following result.

Proposition 1. For any set of positive initial conditions (2.3), the ODE system (2.2) has a unique nonnegative solution which is ultimately uniformly bounded for all $t \geq 0$.

Proof. First, we observe that the right-hand sides of the ODE system (2.2) are continuous with respect to the state variables $S, L, I,$ and $T.$ Therefore, this system has a unique solution for any initial condition.

To prove that all trajectories of the system (2.2) engendered by nonnegative initial conditions (2.3) remain nonnegative for all $t \geq 0$, let us rewrite the Eq (2.2a) in the following form:

Then using the integrating factor, we obtain

This yields

and finally leads to

A similar rationale can be used to show that $L(t) \geq 0, I(t) \geq 0,$ and $T(t) \geq 0$ for all $t \geq 0.$ Thus, the lower bound for all solutions of the system (2.2) is zero.

To establish the upper bounds, let us consider the equation for the total population (2.1) which is obtained by summing up the four ODEs of the system (2.2):

From the above equation, we obtain

Thus, if $N_0 \leq \dfrac{\Lambda}{\mu}$, we have

Otherwise, if $N_0 > \dfrac{\Lambda}{\mu}$, then $\dfrac{d N}{dt} < 0$ for some $t \in [0, \bar{t}]$ meaning that $N(t)$ is decreasing towards $\dfrac{\Lambda}{\mu}$ and its upper bound is $N_0$. Thus, we conclude that

3.   Local analysis of the model

A direct consequence of Proposition 1 is the existence of the biologically feasible region

that contains all possible solutions of the system (2.2) engendered by the initial conditions that satisfy $S_0 + L_0 + I_0 + T_0 \leq \dfrac{\Lambda}{\mu}$. Therefore, $\Omega$ is positively invariant and attracting. The latter implies that $\Omega$ contains all possible equilibria of the system (2.2).

In the present section, we derive the basic reproductive number of the TB transmission model (2.2) and establish the conditions for existence of possible equilibria of the system (2.2) that are nonnegative solutions of the following algebraic system:

3.1.   Disease-free equilibrium and basic reproductive number $\mathscr{R}_0$

The disease-free equilibrium (DFE) of the system (2.2) is the solution of (3.2) with $L = I = T = 0$, that is,

It is easy to see that ${\bf{E}}^0 \in \Omega$ and this equilibrium always exists.

In epidemiology, one of the core metrics related to the speed of the disease propagation in the human population is the so-called basic reproductive number $\mathscr{R}_0 > 0$. This quantity expresses the expected number of secondary infections produced by one infectious individual in a completely susceptible population during his/her entire period of the infectiousness ^[29].

For compartmental epidemiological models, including the proposed model (2.2), the basic reproductive number $\mathscr{R}_0$ can be calculated as the largest eigenvalue (or spectral radius) of the next-generation matrix evaluated at the disease-free steady state ^[30]. Following this approach, let us first define the state sub-vector ${\bf{Y}} \in \mathbb{R}^3_+$ that contains the classes (or compartments) of human hosts carrying the infection, that is,

Second, we extract three differential equations corresponding to the components of ${\bf{Y}}$ from the ODE system (2.2) and write them in the following form

where $\mathcal{F}({\bf{Y}}) \geq 0$ represents the rate of appearance of new infections (or rate of the disease transmission) and $\mathcal{V}({\bf{Y}}) \geq 0$ stands for the rate of the disease transition:

Third, we evaluate the Jacobian matrices of $\mathcal{F}({\bf{Y}})$ and $\mathcal{V}({\bf{Y}})$ at the disease-free equilibrium ${\bf{E}}^0 = \left(\dfrac{\Lambda}{\mu}, 0, 0, 0 \right)$ and obtain

Then the next-generation matrix $\mathscr{F} \mathscr{V}^{-1}$ can be written as

where

denote the elements of $\mathscr{V}$ located on its main diagonal, and

Notably, the characteristic polynomial $\mathcal{P}(\lambda)$ of next-generation matrix $\mathscr{F} \mathscr{V}^{-1}$ admits the following form

Given the structure of $\mathcal{P}(\lambda)$, there are two eigenvalues of $\mathscr{F} \mathscr{V}^{-1}$ which are equal to zero, and the other eigenvalue is exactly the spectral radius of $\mathscr{F} \mathscr{V}^{-1}$:

Thus we have

The above formula expresses the average number of all secondary infections produced by one infectious individual and, therefore, accounts for the fast and slow development of the disease. Let us recall that the parameter $p$ denotes the fraction of effective contacts leading to the fast TB infections, while $(1-p)$ denotes the fraction of effective contacts leading to the slow TB infections. Thus, an alternative form of the formula (3.4) is

where

denotes the basic reproductive number of fast tuberculosis, and

denotes the basic reproductive number of slow tuberculosis. In the expressions (3.6), (3.7), the positive quantity $D$ is defined by (3.3b).

The stylized definitions of $\mathscr{R}_{0f}$ and $\mathscr{R}_{0f}$ given by (3.6) and (3.7) will play a notable role in the process of finding other possible solution(s) to (3.2) as illustrated in the next subsection.

3.2.   Existence of the endemic equilibrium

Our goal is to find out whether the algebraic system (3.2) possesses another feasible solution $(S^*, L^*, I^*, T^*) \in \Omega$ besides the disease-free equilibrium ${\bf{E}}^0 = \left(\dfrac{\Lambda}{\mu}, 0, 0, 0 \right)$.

First, from the Eqs (3.2d) and (3.2a) we can express $T^*$ and $S^*$, respectively:

Then, by plugging (3.8a) into (3.2b) we can also express

Now, in Eq (3.2c) we replace $T^*$ and $L^*$ by their underlying expressions (3.8a) and (3.8c) to obtain

The above equation has an obvious solution $I^{*} = 0$ which leads us to the disease-free equilibrium ${\bf{E}}^0 = \left(\dfrac{\Lambda}{\mu}, 0, 0, 0\right)$, and another possible solution is the one satisfying

Using the definitions of $T^{*}$ and $S^*$ given by Eqs (3.8a) and (3.8b), we have

Thus, the expression (3.9) may also be written as

where

Applying the definitions of $\mathscr{R}_{0s}, \mathscr{R}_{0f}, \mathscr{R}_0,$ and $D$ (see Eqs (3.3b) and (3.5)–(3.7), respectively), let us try to simplify the expressions for $P$ and $Q$ defined by (3.11). Notably, it follows from (3.11a) that

The first and second terms in (3.12) can be rewritten as

and

respectively. Hence, their sum renders

Similarly, the third and fourth terms in (3.12) can be written as

and

respectively. Hence, their sum renders

Thus, we obtain the final simplified expression for $P$ by summing up (3.13) and (3.14), that is,

To simplify the expression for $Q$ defined by (3.11b), we rewrite it as follows:

Now we have

instead of (3.11), and their substitution into (3.10) leads to

or, equivalently,

that can be rewritten as

Solving this last expression for $I^*$, we then obtain

Finally, we proceed to formulate the following result regarding the existence of a strictly positive endemic equilibrium ${\bf{E}}^*: = \big(S^*, L^*, I^*, T^* \big)$ of the ODE system (2.2).

Proposition 2. When $\mathscr{R}_0 > 1$, the ODE system (2.2) has a strictly positive endemic equilibrium ${\bf{E}}^* = \big(S^*, L^*, I^*, T^* \big) \in \Omega$ defined by

Furthermore, ${\bf{E}}^* \in \mathit{{\rm{int}}} \Omega$ if $\alpha > 0$, and ${\bf{E}}^* \in \partial \Omega$ if $\alpha = 0$.

Proof. In the first place, we note that Eq (3.16c) for $I^*$ has been already derived (see Eq (3.15) above). The coordinates $S^*, L^*,$ and $T^*$ of ${\bf{E}}^*$ can be obtained by plugging $I^*$ into Eq (3.8) (the underlying computations are omitted here). From the form of Eq (3.16), it is easy to conclude that all coordinates of ${\bf{E}}^*$ are positive whenever $\mathscr{R}_0 > 1$. Moreover, $S^* < \dfrac{\Lambda}{\mu}$ only if $\mathscr{R}_0 > 1$.

To prove that ${\bf{E}}^* \in \Omega$, let us define $N^* : = S^* + L^* + I^* + T^* > 0$ and show that $N^* \leq \dfrac{\Lambda}{\mu}$. Notably, $S^*$ admits the following form:

Using this form for $S^*$ together with the formulas (3.16b)–(3.16d), we have

Applying the definition of $\mathscr{R}_0$ given by (3.4) together with the identity $p \mu + \epsilon = p(\epsilon + \mu) + (1-p) \epsilon$, the term inside the square brackets in the first row of (3.17) can be rewritten as

Thus, Eq (3.17) becomes

The above relationship clearly indicates that $N^* < \dfrac{\Lambda}{\mu}$ when $\mathscr{R}_0 > 1$ and $\alpha > 0$, meaning that ${\bf{E}}^*$ is an interior equilibrium with respect to $\Omega$. On the other hand, if $\mathscr{R}_0 > 1$ and $\alpha = 0$, then $N^* = \dfrac{\Lambda}{\mu}$ meaning that ${\bf{E}}^*$ is a boundary equilibrium with respect to $\Omega$. Indeed, for $\alpha = 0$ we have that $N(t)$ defined by (2.1) solves the differential equation $\dfrac{d N}{d t} = \Lambda - \mu N.$

4.   Stability properties of the TB transmission model

In the previous section, it has been shown that the ODE system (2.2) possesses only one equilibrium ${\bf{E}}^0$ contained in $\Omega$ if $\mathscr{R}_0 \leq 1$ and has two equilibria, ${\bf{E}}^0$ and ${\bf{E}}^*$ (both contained in $\Omega$) if $\mathscr{R}_0 > 1$.

In the present section, we analyze the long-term behavior of the TB transmission system (2.2) and establish the stability properties of its equilibria. Notably, the basic reproductive number $\mathscr{R}_0$ defined by (3.4) will play the core role in this analysis.

4.1.   Stability for the disease-free equilibrium

Let us recall that $\mathscr{R}_0 < 1$ implies that one infectious individual produces, on average, less than one new infection during his/her period of contagiousness. Furthermore, when $\mathscr{R}_0 < 1$, the dynamical system (2.2) has only the disease-free equilibrium ${\bf{E}}^0$. Thus, it is expected that $\big(S(t), L(t), I(t), T(t) \big) \to {\bf{E}}^0$ as $t \to \infty$ in the "component-by-component" sense, meaning that the disease will eventually die out regardless of its current state. This intuitive rationale is formalized by the following theorem.

Theorem 1 (Stability properties of the disease-free equilibrium ${\bf{E}}^0$). The disease-free equilibrium ${\bf{E}}^0 = \left(\dfrac{\Lambda}{\mu}, 0, 0, 0 \right)$ is globally asymptotically stable when $\mathscr{R}_0 < 1$ and is unstable when $\mathscr{R}_0 > 1$.

Proof. Let us suppose that $\mathscr{R}_0 < 1$. Under this condition, the disease-free equilibrium ${\bf{E}}^0$ is the unique steady state of (2.2) and ${\bf{E}}^0 \in \Omega$. To establish global stability of ${\bf{E}}^0$ we employ the result proposed in ^[31]. First, we write the vector of states as $(S, L, I, T): = ({\bf{X}}, {\bf{Y}})$ where ${\bf{X}}: = S$ and ${\bf{Y}}: = (L, I, T)$ denote the noninfected and infected classes or compartments, respectively. Using this notation, the disease-free equilibrium becomes ${\bf{E}}^0 = ({\bf{X}}^*, {\bf{0}})$ where ${\bf{X}}^* = \dfrac{\Lambda}{\mu}$. Second, the dynamical system (2.2) should be put into the form

The latter is accomplished by defining

According to ^[31], the following conditions must be met to guarantee both local and global stability of ${\bf{E}}^0 = ({\bf{X}}^*, {\bf{0}})$:

(H1) ${\bf{X}}^*$ is globally asymptotically stable equilibrium of the system $\dfrac{d {\bf{X}}}{dt} = F({\bf{X}}, {\bf{0}})$.

(H2) If $\Omega$ is the set where the model (4.1) makes biological sense, then for all $({\bf{X}}, {\bf{Y}}) \in \Omega$ the function $G({\bf{X}}, {\bf{Y}})$ in (4.1b) admits the following form:

where $\mathcal{A}: = \dfrac{\partial G}{\partial {\bf{Y}}} \bigg|_{ {\bf{E}}^0}$ is a Metzler matrix (all its off-diagonal elements are nonnegative).

Third, we verify the conditions (H1)–(H2) stated above bearing in mind that $\Omega \subset \mathbb{R}^4_+$ is defined by (3.1). Notably, for $F({\bf{X}}, {\bf{Y}})$ defined by (4.2a), we have that $\dfrac{d {\bf{X}}}{dt} = F({\bf{X}}, {\bf{0}})$ reduces to a single ODE

This equation possesses a unique equilibrium ${\bf{X}}^* = \dfrac{\Lambda}{\mu},$ which is globally asymptotically stable for all ${\bf{X}} \geq 0$. Thus, condition (H1) holds. Regarding the condition (H2), we observe that (4.2b) can be written as $G({\bf{X}}, {\bf{Y}}) = \mathcal{A} {\bf{Y}} - \widehat{G}({\bf{X}}, {\bf{Y}})$ with $\mathcal{A}$ and $\widehat{G}({\bf{X}}, {\bf{Y}})$ given by

and

respectively. It is immediate to check via (4.3) that $\mathcal{A}_{ij} \geq 0$ when $i \neq j$ and $i, j = 1, 2, 3$. Therefore, $\mathcal{A}$ is a Metzler matrix. From (4.4) it also follows that $\widehat{G}({\bf{X}}, {\bf{Y}}) \geq {\bf{0}}$ for all $({\bf{X}}, {\bf{Y}}) \in \Omega$ since it holds that $0 \leq S \leq \dfrac{\Lambda}{\mu}.$ Thus, both conditions (H1)–(H2) are verified, and this completes the proof that ${\bf{E}}^0$ is a globally asymptotically stable equilibrium of the system (2.2) whenever $\mathscr{R}_0 < 1.$

Now, let us suppose that $\mathscr{R}_0 > 1.$ To prove that ${\bf{E}}^0$ is unstable under this condition, it is sufficient to show that the Jacobian matrix $\mathcal{J} ({\bf{E}}^0)$ of the system (2.2) evaluated at ${\bf{E}}^0$ has at least one eigenvalue with positive real part whenever $\mathscr{R}_0 > 1$. Using the positive quantities $A, B,$ and $C$ introduced by the relationships (3.3a) in Section 3, it is easy to deduce that

Let us also recall that

where $\lambda_i, i = 1, 2, 3, 4$ denote four eigenvalues of $\mathcal{J} ({\bf{E}}^0)$. Therefore, it suffices to show that $\det \mathcal{J} ({\bf{E}}^0) < 0$ whenever $\mathscr{R}_0 > 1$. Indeed,

In the above expression, we observe that

in virtue of the relationship (3.3b). Furthermore, $\mathscr{R}_0$ admits the form

Thus, we have

Therefore, $\det \mathcal{J} ({\bf{E}}^0) < 0$ whenever $\mathscr{R}_0 > 1$ and ${\bf{E}}^0$ is unstable if $\mathscr{R}_0 > 1$. This completes the proof of Theorem 1.

It is worthwhile to note that ${\bf{E}}^0$ cannot be a repeller when $\mathscr{R}_0 > 1$ since $\mathcal{J} ({\bf{E}}^0)$ possesses at least one strictly negative eigenvalue ($\lambda_1 = - \mu < 0$). Thus, ${\bf{E}}^0$ is a saddle point. In fact, when $\mathscr{R}_0 > 1$ the solution of the system (2.2) engendered by the initial conditions $(S_0, 0, 0, 0) \neq {\bf{E}}^0$ with $S_0 \neq \dfrac{\Lambda}{\mu}$ converges to ${\bf{E}}^0.$

4.2.   Stability of the endemic equilibrium

Let us recall that $\mathscr{R}_0 > 1$ implies that one infectious individual produces, on average, more than one new infection during his/her period of contagiousness. It was also shown in Section 3 that, when $\mathscr{R}_0 > 1$, the dynamical system (2.2) has two equilibria: the disease-free equilibrium ${\bf{E}}^0$ and the endemic equilibrium ${\bf{E}}^* \in \Omega$ whose coordinates are strictly positive and defined by (3.16). Moreover, Theorem 4.1 has established that the disease-free equilibrium ${\bf{E}}^0$ is unstable when $\mathscr{R}_0 > 1$ meaning that ${\bf{E}}^0$ is hardly reachable when $\mathscr{R}_0 > 1$. Thus it is expected that $\big(S(t), L(t), I(t), T(t) \big) \to {\bf{E}}^*$ as $t \to \infty$ in the "component-by-component" sense, the endemic equilibrium ${\bf{E}}^*$ can be eventually reached, and the disease may persist regardless of its current state. This intuitive rationale is formalized by the following theorem.

Theorem 2 (Stability properties of the endemic equilibrium ${\bf{E}}^*$). When $\mathscr{R}_0 > 1$, the endemic equilibrium ${\bf{E}}^* = \big(S^{*}, L^{*}, I^{*}, T^{*} \big)$ is globally asymptotically stable in the open set $\Omega_* \subset \Omega$ where

Proof. As stated by Proposition 2, the endemic equilibrium ${\bf{E}}^* = \big(S^{*}, L^{*}, I^{*}, T^{*} \big)$ defined by (3.16) exists whenever $\mathscr{R}_0 > 1$ and ${\bf{E}}^* \in \Omega_*$. To prove that ${\bf{E}}^*$ is globally stable in $\Omega_*$, we will use the so-called Lyapunov-LaSalle theorem (see, e.g., Theorem 6.2 in ^[32]). The first step is to construct a Lyapunov function for the ODE system (2.2). Let us recall that a continuously differentiable scalar function $V: \mathbb{R}^4_+ \mapsto \mathbb{R}$ is a Lyapunov function of model (2.2) if it has the following properties ^[33]:

(i) $V({\bf{E}}^*) = 0$;

(ii) $V(S, L, I, T)$ is radially unbounded in $R_+^4$;

(iii) $\dfrac{d}{dt} V(S, L, I, T) \leq 0$ for all $(S, L, I, T) \in \mathbb{R}_+^4$ and $t \geq 0$.

A good candidate for $V(S, L, I, T)$ is the separable scalar function proposed by B. Goh [33, p10] which is radially unbounded and has been used by many scholar in the context of TB transmission models (see, e.g., ^{[20, 21, 27, 34]} and references therein). Therefore, let us consider the following form of $V$:

where $B_1$, $B_2,$ and $B_3$ are positive constants to be determined. Notably, function $V$ defined by (4.5) fulfills items (i) and (ii) mentioned above, while item (iii) requires further considerations. Namely, we should find the proper positive values for $B_1$, $B_2,$ and $B_3$ so that the orbital derivative of $V$ (i.e., the time derivative of $V$ along all solutions of the system (2.2)) be nonpositive. For further analysis, it will be convenient to consider an alternative form of (3.2) which is satisfied by the coordinates $(S^*, L^*, I^*, T^*)$ of $E^*$:

These relationships will play an important role in the sequel. Using the chain rule together with (4.6), the orbital derivative of $V$ can be written as

where the first term is

The second and the third terms in the right-hand side of (4.7) can be written, respectively, as

and

while the last term of (4.7) gives

Using the expressions from the four previous formulas and making some laborious rearrangements, the orbital derivative of $V$ can be written as

Now the positive constants $B_1$, $B_2,$ and $B_3$ should be chosen in a way that the coefficients of $S(I+\phi T)$, $L$, $I,$ and $T$ in the above expression become equal to zero, that is, the following four equations must be satisfied:

This is a tedious task, so we present here the final solution of (4.8), while the meticulous details are moved to Appendix A:

With $B_1$, $B_2$ and $B_3$ defined by (4.9), the orbital derivative of $V$ becomes

For convenience, we introduce new variables

to eliminate $S, L, I$ and $T$ in the expression for $\dfrac{d V}{dt}$. Here, we observe that for $B_1$ and $B_2$ satisfying the Eq (4.8a) it holds that

Keeping in mind this relationship, the orbital derivative of $V$ can be written as

Following the idea proposed in ^[27], we multiply the Eq (4.6b) by $B_1$ and the Eq (4.8d) by $L^*$ in order to obtain

Hence, it follows that

Similarly, we multiply the Eq (4.6d) by $B_3$ and the Eq (4.8c) by $T^*$ and obtain

Then we have

Let $F_1({\bf{Z}})$ and $F_2({\bf{Z}}), {\bf{Z}}: = (x, y, z, u)$ be two scalar functions which will be defined later. Then it is also fulfilled that

in virtue of (4.11) and (4.12). By summing the above expressions (both equal to zero) to the right-hand side of formula (4.10) and rearranging some terms, the orbital derivative of $V$ becomes

To make vanish the last two summands in the above formula, we choose the functions $F_1({\bf{Z}})$ and $F_2({\bf{Z}})$ in the following way:

Using this selection of $F_1({\bf{Z}})$ and $F_2({\bf{Z}})$, the orbital derivative of $V$ is now written as

Finally, to demonstrate that $\dfrac{d V}{dt} \leq 0$, we observe that the first term in the right-hand side of (4.14) is strictly negative whenever $x \neq 1$ or, equivalently, whenever $S \neq S^*$. However, this term vanishes when $S = S^*$. Further, we apply the so-called "arithmetic mean – geometric mean inequality" (see, e.g., ^[35] or similar textbooks) to all other terms on the right-hand side of (4.14) and obtain that they are nonpositive for any positive $(x, y, z, u)$ and vanish only if

Thus $V = V(S, L, I, T)$ defined by (4.5) with positive constants $B_1, B_2$ and $B_3$ satisfying (4.9) is a Lyapunov function for the system (2.2) and we also have

where

Notably, the largest invariant subset of $\Omega_0$ is the singleton ${\bf{E}}^*$ (endemic equilibrium). Then in virtue of the Lyapunov-LaSalle theorem [32, Theorem 6.2, p. 138], every solution of (2.2) engendered by an initial condition $\big(S_0, L_0, I_0, T_0 \big) \in \Omega_*$ converges to ${\bf{E}}^*$ when $t \to \infty.$ In other words, ${\bf{E}}^*$ is globally asymptotically stable in $\Omega_*$. This completes the proof of Theorem 2.

Summarizing the results presented so far, we conclude that the number of equilibria of the system (2.2) located in $\Omega$, together with their stability, is determined by the value of the basic reproductive number $\mathscr{R}_0$. It is also clear that the dynamical system (2.2) undergoes a forward (or transcritical) bifurcation at $\mathscr{R}_0 = 1$ which is depicted in Figure 2. In the bifurcation diagram, the values of $\mathscr{R}_0$ are located on the horizontal axis, and the vertical axis corresponds to the equilibrium size of population groups contributing to the disease spread (that is, $L^* + I^* + T^*$). In Figure 2, an asymptotically stable equilibrium (${\bf{E}}^0$ if $\mathscr{R}_0 < 1$ or ${\bf{E}}^*$ if $\mathscr{R}_0 > 1$) is depicted by a red solid line, while a black dashed line displays the unstable equilibrium ${\bf{E}}^0$ only if $\mathscr{R}_0 > 1$.

Thus, when $\mathscr{R}_0$ changes its value due to perturbation of some parameter(s) included in its expression (3.4), either from less than unity to greater than unity or vice versa, it alters the long-term behavior of the system (2.2).

5.   Numerical solutions of the model and the role of treatment abandonment

In the preceding sections, it has been shown that the four-dimensional TB transmission model (2.2) is biologically meaningful and well-posed. Furthermore, this model exhibits rather explicit stability properties tightly related to the basic reproductive number $\mathscr{R}_0$ that are highly desirable in epidemiological modeling.

Let us briefly illustrate that the proposed model (2.2) exhibits the behavior predicted by Theorems 4.1 and 4.2 given in Section 4. Having assigned plausible values to all constant parameters of the model (2.2) (see Table 2), we proceed to solve numerically the ODE system (2.2) using Mathematica 12 software tool under two illustrative scenarios:

Scenario 1: TB transmission without treatment abandonment ($k = 0$);

Scenario 2: TB transmission with treatment abandonment ($k = 0.3 > 0$).

We first recall that tuberculosis is a social disease that is easily spread under crowded conditions of big cities ^[1]. Therefore, we take as an example a city with an average population of $2.5$ million people and scale the parameter $\Lambda$ accordingly. The average life expectancy varies between $66$ and $88$ depending on the country, and we assume it is equal to $77$ years^‡. Furthermore, in 2019, the TB incidence per year varied between $0$ and $654$ cases per $100 \:000$ inhabitants in different countries^§, and we fit the value of $\beta$ to match the TB incidence of approximately $250$ cases per $100 \:000$ inhabitants per year. Numerical values of other parameters and their underlying ranges are borrowed from the existing literature.

^‡For more information on life expectancy, please refer to the WORLDOMETER website (https://www.worldometers.info/demographics/life-expectancy/)

^§More information is available at The World Bank Data (https://data.worldbank.org/indicator/SH.TBS.INCD)

For both scenarios, we assign the same set of initial conditions that are relatively close to the endemic equilibrium ${\bf{E}}^*$. The coordinates of ${\bf{E}}^{*}$ are calculated using four formulas of (3.16):

The initial conditions should fulfill the relationship $N(0) = S(0)+L(0)+I(0)+T(0) = 2.5 \times 10^6$, and we assume that

We are interested in the time evolution of two infectious groups: $I$-class (fully capable of transmitting the disease) and $T$-class (partially infectious). Their profiles are presented in Figure 3 for short-term (left column) and long-term evolution (right column).

Notably, Figure 3 reveals different short- and long-term tendencies of the disease propagation for Scenario 1 (dashed-line curves) and Scenario 2 (solid-line curves), and this has a straightforward explanation from the standpoint of Theorems 1 and 2. Namely, since

the trajectories of the system (2.2) converge either to the disease-free equilibrium ${\bf{E}}^0$ (under Scenario 1) or to the endemic equilibrium ${\bf{E}}^*$ (under Scenario 2). Thus, theoretical findings presented in Sections 3 and 4 are confirmed by numerical simulations of the proposed TB transmission model.

Let us now recall the primary reason why this model has been proposed. Namely, we intended to include the group of patients who abandon treatment before its completion, without assigning them to an additional class and, thus, retaining only the four original compartments of the model ($S, L, I$ and $T$). Therefore, we proceed now to analyze the effect of treatment abandonment on disease transmission.

The model should reflect that TB-positive patients who abandon treatment before completion contribute to the disease spread by passing from "partially infectious" to "fully infectious". Let us verify if the proposed model (2.2) fulfills this expectation by analyzing the partial derivative of $\mathscr{R}_0$ with respect to $k$, the parameter expressing the rate of treatment abandonment:

Thus $\mathscr{R}_0$ is an increasing function of $k$ whenever

that is, for all $\alpha$ such that

Notably, the relationship (5.3) holds trivially if $\alpha = 0$, i.e. if there is no disease-induced mortality. Furthermore, for any positive value of $\alpha$ below the threshold $\alpha^*$ (defined by formula (5.4) above), an increase in the rate of treatment abandonment $k$ would increase the value of $\mathscr{R}_0$. In other words, by increasing $k$, a greater average number of secondary TB infections could be expected from one person carrying active bacilli.

However, if $\alpha > \alpha^*$, the basic reproductive number $\mathscr{R}_0$ is a decreasing function of the treatment abandonment rate $k$. In this case, a smaller average number of secondary TB infections could be expected from one actively infected person if $k$ is increased. The rationale behind this statement is rather simple and intuitive. When the disease-induced mortality is high ($\alpha > \alpha^*$), more people are removed from the $I$-class per unit time, meaning that they stop spreading the infection. As $k$ increases, more people pass from the "partially infectious" $T$-class to the "fully infectious" $I$-class, from which they are fastly removed (because $\alpha$ is high) instead of remaining in the system and bearing a reduced capacity of infecting the others. As a consequence, when $k$ increases, one infectious individual produces, on average, a smaller number of secondary infections, meaning that $\mathscr{R}_0$ is a decreasing function of $k$ when $\alpha > \alpha^*$. Furthermore, it is easy to check that $\mathscr{R}_0$ is a decreasing function of $\alpha$:

Let us now examine the behavior of $\mathscr{R}_0$ as a function of $k$ (treatment abandonment rate) and $\alpha$ (disease-induced mortality rate), which is denoted by $R_0: = R_0(k, \alpha)$ and displayed by a yellow-colored surface in Figure 4. In this figure, the green-colored horizontal plane indicates the threshold surface $\mathscr{R}_0 = 1$. The intersection between $R_0(k, \alpha)$ and $\mathscr{R}_0 = 1$ is marked by a red-colored curve.

When both arguments of the function $R_0(k, \alpha)$ take values corresponding to Scenario 1, that is, $k = 0$ and $\alpha = 0.03$, the respective point $R_0(0, 0.03)$ lies below the threshold plane $\mathscr{R}_0 = 1$. On the other hand, the point $R_0(0.3, 0.03)$ corresponding to Scenario 2 is located above the threshold plane $\mathscr{R}_0 = 1$. It is also displayed in Figure 4 that if $k = 0.3$ but $\alpha$ is enhanced to $\alpha = 0.15$, the respective point $R_0(0.3, 0.15)$ of the yellow-colored surface $R_0(k, \alpha)$ lies below the green-colored threshold plane $\mathscr{R}_0 = 1$.

Thus, from a purely mathematical standpoint, the current value of $\mathscr{R}_0$ can be driven below $1$ either by reducing the treatment abandonment rate $k$ or by augmenting the disease-induced mortality rate $\alpha$. Naturally, the former seems a lot more reasonable and compelling than the latter, meaning that even a small reduction in $k$ (ceteris paribus) can certainly abate the disease propagation. Two charts displayed in Figure 5 depict these tendencies for variation of each parameter within its range, while other parameters retain their values determined in Table 2. Namely, the left chart shows the change in the value of $\mathscr{R}_0$ as a function $R_0^k (k)$ that depends only on $k \in [0, 0.5]$. Similarly, the right chart shows the change in the value of $\mathscr{R}_0$ as a function $R_0^{\alpha} (\alpha)$ that depends only on $\alpha \in [0, 0.3]$. The black points on both charts of Figure 5 mark the critical values $k^c = 0.1771$ and $\alpha^c = 0.0493$ such that $R_0^k (k^c) = 1$ and $R_0^{\alpha} (\alpha^c) = 1$, respectively.

To conclude this section, let us assess the negative effect of the treatment abandonment in terms of additional active TB infections it may cause and additional deaths it may induce. For that purpose, let us introduce two auxiliary variables:

The first one, $C(t)$, denotes the cumulative incidence and expresses the overall number of active TB infections that appear during the period $[0, t]$. The second one, $M(t)$, denotes the cumulative disease-induced mortality and expresses the overall number of TB-infected people who died of the disease during the period $[0, t]$. Notably, in the above expressions, the variables $S, L, I$ and $T$ are solutions of the system (2.2) engendered by the initial conditions (5.1) when other parameters take values from Table 2. Hence $C(t)$ and $M(t)$ can be viewed as solutions of the following differential equations with underlying initial conditions

Thus, we can compose a six-dimensional ODE system by complementing the original model (2.2) with two auxiliary Eqs (5.5) and solve it numerically using the Mathematica software package for two scenarios established earlier. Figure 6 exhibits the cumulative incidence $C(t)$ (left chart) and cumulative mortality $M(t)$ (right chart) under Scenario 1 (dashed-line curves) and Scenario 2 (solid-line curves) when $t \in [0, 3].$

Our numerical simulations also allow to estimate, through the model, how many new infections and fatalities may occur by the end of the observation period ($3$ years) under both scenarios. In effect, we have that

Thus, if the treatment abandonment rate is reduced from $k = 0.3$ to $k = 0$ and then for three years no patient abandons treatment, it may help to avoid about $3725$ active TB infections together with $44$ TB-induced fatalities in a population of $2.5$ million people.

6.   Conclusions

In the present paper, we have proposed and justified a synthesized version of the tuberculosis transmission model that accounts for treatment abandonment. A distinctive feature of this model is that it contains only four standard state variables or compartments expressing the classes of susceptible, latently infected, actively infected, and treated individuals, while other models accounting for treatment abandonment contain more state variables ^{[19, 20, 21, 22, 23]}.

It was also shown that the proposed model is biologically meaningful and well-posed from a mathematical standpoint. We have also rigorously proved that the proposed model exhibits the properties of global stability that are highly desirable in epidemiological modeling, and this constitutes another characteristic attribute of our model. The qualitative analysis of the system behavior with treatment desertion is helpful for a better understanding of the endemic persistence of tuberculosis. It also explains why this infectious disease is difficult to exterminate.

Furthermore, the proposed model enables us to visualize that the treatment abandonment enhances the incidence of the disease and disease-induced mortality. This is an important feature that may motivate the patients with active TB not to abandon their ongoing treatment and thus to avoid new infections and fatalities. On the other hand, the model also revealed that a reduction of the treatment abandonment rate has a positive effect on the disease incidence and results in avoiding disease-related fatalities.

Once our lower-dimensional model with treatment desertion is rigorously justified from the theoretical standpoint, this model can be further used for fitting the realistic data of active TB infections detected, patients being treated, and treatment abandonment cases. Since tuberculosis is a social disease, its spread depends on social strata that may feature lower/higher transmission and treatment abandonment rates. In this context, adding social heterogeneity through the metapopulation or network-based modeling approaches ^{[39, 40]} in combination with parameter estimations ^[41] brings forward a germane outlook for our future research.

Acknowledgements

This research was funded by the National Fund for Science, Technology, and Innovation (Autonomous Heritage Fund Francisco José de Caldas) by way of the Research Program No. 1106-852-69523 (Principal Investigator: Hector J. Martinez-Romero), Contract: CT FP 80740-439-2020 (Colombian Ministry of Science, Technology, and Innovation-Minciencias), Grant ID: CI-71243 (Universidad del Valle, Colombia). Hanner E. Bastidas-Santacruz acknowledges support from Universidad del Valle through Project CI-71309 for institutional support of Master course students.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix A: Solution of the system (4.8)

For convenience, the algebraic System (4.8) is replicated here

and we immediately observe that system (A-1) is over-determined for it has four equations in only three unknowns $B_1, B_2,$ and $B_3.$ Nonetheless, we can prove that the positive constants $B_1, B_2$ and $B_3$ defined by (4.9) are solutions of all four equations of (A-1).

First, from Eq (A-1d) we obtain

and plugging this expression into Eq (A-1a) we get

so that

coincide with (4.9a) and (4.9b).

Further, Eq (A-1b) multiplied by $- \phi$ is added to Eq (A-1c) and yields

This formula agrees with (4.9c).

Let us now verify that $B_1, B_2,$ and $B_3$ defined by (4.9) effectively satisfy all four equations of the System (A-1) (which are identical to (4.8)). Verification of Eqs (A-1a) and (A-1d) seems rather trivial and we omit it here. To verify (A-1b) and (A-1c), we recall the definition of $S^{*}, \mathscr{R}_0$ and $D$ provided by the formulas (3.16a), (3.4) and (3.3b), respectively.

From the above relationships, we can also deduce an alternative form of $B_2$ that will be used in the sequel:

We begin by proving that Eq (A-1b) is satisfied:

The latter implies that Eq (A-1b) is fulfilled by $B_1, B_2$ and $B_3$ defined by (4.9). To complete our verification, we check the Eq (A-1c):

Thus, we have proven that positive quantities $B_1, B_2$ and $B_3$ defined by (4.9) are solutions of all four Eqs (A-1).