Optimal Control Approach for Modeling the Introduction of Sterilizing Treatment for Tuberculosis in Low- and Middle-Income Countries

1. Introduction

Tuberculosis has been a significant public health concern for centuries, and it remains a burden in many middle- and low-income countries. Tuberculosis (TB) is a bacterial infectious disease caused by Mycobacterium tuberculosis (Mtb), a rod-shaped (bacillus) bacterium also known as Koch’s bacilli. It primarily affects the lungs (pulmonary TB) but can also spread to other organs (extrapulmonary TB). TB is transmitted through airborne droplets when an infected person coughs, sneezes, or speaks, releasing tubercle bacilli into the air.

Notably, TB infections can be latent or active. Latent TB infection leaves the bacilli dormant in the body without causing symptoms, but they can reactivate, mainly due to malnutrition, comorbidities, or immunosuppressive treatment. Latently infected individuals are not infectious; they do not transmit the bacilli to others. People with active bacilli that multiply in their bodies are said to undergo active TB infection. They are contagious and usually show symptoms, such as persistent cough, weight loss, fever, and night sweats. A simple X-ray test may confirm the presence of active bacilli at an earlier non-symptomatic disease stage. However, health insurance will only cover such an X-ray test if a person exhibits some alarming symptoms. Thus, only symptomatic patient may receive treatment in many middle- and low-income countries. If untreated, TB can be fatal.

Understanding the distinction between sterilizing and non-sterilizing activities of antituberculous drugs is crucial in treating active TB infections. Sterilizing activity refers to a drug’s ability to eliminate all persistent bacilli (both active and dormant) that can cause relapse if not eradicated. On the other hand, non-sterilizing (or bactericidal) activity targets actively replicating bacteria, reducing the bacterial load but not necessarily eliminating all dormant forms [1,2,3]. It is also important to underline that there is no long-lasting immunity to TB after recovering from an active disease.

In simpler terms, sterilizing treatment leaves the patients bacilli-free by eradicating all active and dormant bacilli. In contrast, after the completion of non-sterilizing treatment and the disappearance of all the disease symptoms, some dormant bacilli may remain in the patient’s body and become active when the patient’s immune system is weakened. In other words, patients who received non-sterilizing treatment turn latently infected after completing the treatment.

The standard strict-scheme course of sterilizing TB treatment approved by the World Health Organization [4], is based on a combination of rifampicin (replaceable by rifapentine or rifabutin) and pyrazinamide, both featuring a strong sterilizing effect, and twp bactericidal drugs (usually isoniazid and ethambutol). Among these drugs, rifampicin and its analogs are the most expensive, and healthcare entities in middle- and low-income countries may not have sufficient funds for their acquisition. As a result, many patients in such countries receive non-sterilizing TB treatment with cheaper bactericidal drugs only. In many middle- and low-income countries, there are no programs aimed at detection of latent TB infections, and such infections may only be occasionally discovered, e.g., before an immunosuppressive treatment for some comorbidities is to be assigned. In general, latent TB is not considered “sickness” because of no symptoms it brings forward for justifying the medical test to run. Thus, we actually have no estimation of how many latent TB infections there are, and the outcomes of non-sterilizing treatment may also leads to the accumulation of latently infected people who may become actively infected later in life, especially when their immune status changes.

Several studies (see [5,6] and references therein) highlight that COVID-19 provoked active TB infections in latently infected individuals through immune suppression, lung damage, immunosuppressive treatments, and socioeconomic stressors during the pandemic attributed to the limited access to healthcare facilities, lack of medical personnel and laboratory space, as well as poverty increase due to economic problems caused by quarantine. This effect was more visible in high-TB burden regions composed of middle- and low-income countries, where latent TB is common. As a result, TB morbidity notably increased in such regions [4]. Therefore, reduction of the latent TB infections is crucial for global TB control and eventual eradication, and healthcare entities should try to provide sterilizing treatment to actively TB-infected patients whenever possible.

This paper addresses the challenges related to the introduction of sterilizing treatment from a mathematical standpoint while also considering the limited funds available for TB treatment in middle- and low-income countries. Our goal is to showcase the future benefits of not only switching entirely to sterilizing therapy (which is an ideal case) but also to assess the effect of implementing a variable (non-constant in time) distribution between the two treatment types within the limits of an assigned budget. Such a distribution has been designed using the optimal control approach.

The paper is organized as follows. Section 2 presents the new TB transmission model featuring two classes of treated individuals, where a constant share defines the distribution between the two treatment types. In Section 3, we conduct the local analysis of the introduced model, including the computation of its basic reproductive number and possible equilibria, and analyze the long-term behavior of the model with a constant share of new patients to undergo each type of treatment.

In Section 4, we introduce a time-variable share that defines the distribution between the two treatment types and formulate an optimal control problem of finding an optimal share that minimizes the total number of new latent and active TB infections under a budget constraint. We also justify the existence of the optimal solution to the underlying optimal control problem and anticipate its structure.

The main practical results are gathered in Section 5. In SubSection 5.1, we estimate and compare the outcomes of both treatment types (solely non-sterilizing and solely sterilizing) from the short-, medium-, and long-term perspectives. Such estimations are then contrasted with those obtained for the optimal time-dependent distributions between the two types of treatment, which are designed based on the optimal control approach. The latter is performed in SubSection 5.2, where we showcase and discuss the numerical solutions to the optimal control problem under budget constraints while assuming the relative costs of sterilizing drugs are either constant or decreasing in time. A practical outlook of the numerical results is given in SubSection 5.3.

Finally, the key outcomes of the present work are outlined in Section 6.

2. The TB Transmission Model with Two Treatment Types

Taking as the basis a simple model developed in [7] that only features non-sterilizing treatment, we now propose its extension that helps model the introduction of sterilizing treatment while still accounting for patients undergoing the habitual non-sterilizing treatment. Our model includes five time-dependent variables or classes of human individuals, each describing a particular group of people detailed below:

• $S\left(t\right)$ – the number of susceptible people (free of bacilli).
• $L\left(t\right)$ – the number of latently infected people (non-infectious carriers of inactive bacilli) with no disease symptoms.
• $I\left(t\right)$ – the number of infected people (infectious carriers of active bacilli) who have disease symptoms but are not being treated yet.
• $T_1\left(t\right)$ – the number of individuals receiving non-sterilizing treatment.
• $T_2\left(t\right)$ – the number of individuals receiving sterilizing treatment.

It is worth pointing out that during the initial phase of treatment (either sterilizing or not), patients may still be infectious, but as the treatment progresses, their infectiousness decreases. Therefore, we may assume that people from the classes $T_1$ and $T_2$ exhibit reduced infectivity, i.e., they are partially infectious.

We also suppose that all five classes are well-mixed and that the disease is transmitted only through encounters between infectious and non-infectious individuals. The total human population is denoted by

$$N\left(t\right):=S\left(t\right)+L\left(t\right)+I\left(t\right)+T_1\left(t\right)+T_2\left(t\right).$$

(1)

for all $t\ge 0$. The dynamic interflow between five disjoint classes $S,L,I,T_1,$ and $T_2$ is expressed through by the dynamical system

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& =\Lambda -\beta S\left(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2\right)+{\delta }_2{T}_2-\mu S\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dL}{dt}}& =(1-p)\beta S\left(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2\right)+{\delta }_1{T}_1-ϵL-\mu L\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dI}{dt}}& =p\beta S\left(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2\right)+ϵL+{\kappa }_1{T}_1+{\kappa }_2{T}_2-\gamma I-\alpha I-\mu I\hfill \end{array}$$

(2c)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{T}_1}{dt}}& =(1-q)\gamma I-{\delta }_1{T}_1-{\kappa }_1{T}_1-\mu T_1\hfill \end{array}$$

(2d)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{T}_2}{dt}}& =q\gamma I-{\delta }_2{T}_2-{\kappa }_2{T}_2-\mu T_2\hfill \end{array}$$

(2e)

with nonnegative initial conditions

$$S\left(0\right)=S_0,L\left(0\right)=L_0,I\left(0\right)=I_0,T_1\left(0\right)=T_0^1,T_2\left(0\right)=0.$$

(3)

The last initial condition in (3) implies that no patients are undergoing sterilizing treatment at the moment $t=0$, and this type of treatment is just being introduced at $t=0$.

Figure 1 displays the flow diagram of the TB transmission model (2), and the detailed descriptions of the model’s parameters are given in Table 1.

Let us briefly explain the model (2), paying attention to its novelty compared to the one developed in [7]. First, in the model (2), we have a share $q\in [0,1]$ of infectious people who are assigned to undergo the sterilizing treatment (per unit time), while the rest of contagious individuals, $(1-q)$, will be sent for habitual non-sterilizing treatment (see Eqs. (2d) and (2e) above). Thus, when $q=0$, we have $T_2\left(t\right)\equiv 0$ for all $t\ge 0$, and the system (2) is reduced to the one thoroughly studied in [7].

Second, we also assume that the rates of treatment completion $({\delta }_i,i=1,2)$ and abandonment $({\kappa }_i,i=1,2)$ are different for the patients from class $T_1$ undergoing non-sterilizing treatment ($i=1$) and those from class $T_2$ assigned for sterilizing treatment ($i=2$). It should also be noted that, after completing the treatment, people from the class $T_1$ become non-infectious but still latently infected (see Eqs. (2d) and (2b)) while people from the class $T_2$ become bacillus-free (see Eqs. (2e) and (2a)). However, by abandoning either treatment type, the patients become actively infected (see Eqs. (2c)-(2e)).

Third, as indicated in [8], patients may remain contagious for some time under both types of treatment, but their infectiousness gradually decreases as the treatment progresses. Let us also recall that non-sterilizing treatment makes the bacilli dormant or inactive while sterilizing treatment kills them. Therefore, it is fair to assume that the infectiousness of people undergoing treatment and belonging to the classes $T_1$ and $T_2$ is reduced, respectively, by factors ${\varphi }_1$ and ${\varphi }_2$ satisfying the relationship

$$0<{\varphi }_2\le {\varphi }_1<1,$$

which stays in line with experimental observations [3,9] and is reflected by Eqs. (2a)-(2c), where the mass-action incidence models the disease transmission through the term $\beta \left(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2\right)$ with $\beta$ denoting, as in [7], the “per capita rate of effective respiratory contacts with infectious people”. Thus, using this definition of $\beta$ introduced in [7] we have

$$\beta \Rrightarrow {\displaystyle \frac{\left(\mathrm{No}.\mathrm{of}\mathrm{contacts}\mathrm{per}\mathrm{person}\mathrm{per}\mathrm{unit}\mathrm{time}\right)\times \left(\mathrm{contagion}\mathrm{probability}\right)}{\mathrm{Total}\mathrm{size}\mathrm{of}\mathrm{the}\mathrm{population}\left(\approx \Lambda /\mu \right)}}.$$

Since there is no evidence of vertical TB transmission [10], we assume that all newborns, recruited at the constant rate $\Lambda >0$, are bacillus-free. Here, we assume all individuals have the same natural mortality rate denoted by $\mu >0$, while individuals from class I additionally exhibit the disease-induced mortality rate denoted by $\alpha >0$.

After inhaling active bacilli, a susceptible person from the class S, may either enter the I-class of contagious people with a probability $p\in (0,1)$ by acquiring an active TB form or pass to the L-class with a probability $(1-p)$ by getting a latent TB infection. Here, the probability p of becoming actively infected and infectious depends on the immune response of the person receiving active bacilli. If a person’s immune system is weakened by comorbidities, malnutrition, or immunosuppressive treatments (for HIV, cancer, or transplants), this person has a higher risk of developing an active TB form [11,12].

People belonging to the class L of latently infected may develop an active TB infection if their immune status gets altered. In such a case, they will progress to the infectious class I at a rate $ϵ>0$ (see Eqs. (2b)-(2c)). In other words, latently TB infected individuals may remain non-infectious during the time $1/ϵ$ before developing an active TB infection.

Contagious people (class I) are then assigned for either sterilizing or non-sterilizing treatment at a rate $\gamma >0$. Thus, a share q of them passes to the class $T_2$, and a share $(1-q)$ passes to the class $T_1$ (see Eqs. (2c)-(2e)). It is worth noting that people carrying active TB bacilli may remain untreated for a time $1/\gamma$ during which medical test are run to confirm the disease and the treatment is assigned.

The well-posedness of the model (2) is guaranteed by the following statement.

Proposition 1. 

Solution to the dynamical system (2), corresponding to any set of nonnegative initial conditions (3), is unique, nonnegative, and bounded for all $t\ge 0$.

The proof of Proposition 1 is almost identical to the one presented in [7] and therefore is left to the readers. Furthermore, Proposition 1 also entails the existence of a biologically meaningful region

$$\Omega :=\left\{\left(S,L,I,T_1,T_2\right)\in {\mathbb{R}}_{+}^5:S+L+I+T_1+T_2\le {\displaystyle \frac{\Lambda }{\mu }}\right\},$$

(4)

which is positively invariant and constitutes the absorbing set of the system (2) because it contains all possible solutions of (2) engendered by the initial conditions satisfying the inequality $S\left(0\right)+L\left(0\right)+I\left(0\right)+T_1\left(0\right)+T_2\left(0\right)\le {\displaystyle \frac{\Lambda }{\mu }}$. Therefore, $\Omega$ must contain all possible equilibria of the system (2).

3. Local Analysis of the Model

The steady states or equilibria of the dynamical system (2) are its constant solutions and can be obtained by solving the following algebraic system:

$$\begin{array}{cc}\hfill 0& =\Lambda -\beta S(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2)+{\delta }_2{T}_2-\mu S\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill 0& =(1-p)\beta S(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2)+{\delta }_1{T}_1-ϵL-\mu L\hfill \end{array}$$

(5b)

$$\begin{array}{cc}\hfill 0& =p\beta S(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2)+ϵL+{\kappa }_1{T}_1+{\kappa }_2{T}_2-\gamma I-\alpha I-\mu I\hfill \end{array}$$

(5c)

$$\begin{array}{cc}\hfill 0& =(1-q)\gamma I-{\delta }_1{T}_1-{\kappa }_1{T}_1-\mu T_1\hfill \end{array}$$

(5d)

$$\begin{array}{cc}\hfill 0& =q\gamma I-{\delta }_2{T}_2-{\kappa }_2{T}_2-\mu T_2\hfill \end{array}$$

(5e)

In this section, we will identify the steady states of the dynamical system (2) as nonnegative solutions of the nonlinear system (Section 3) belonging to the absorbing set $\Omega \subset {\mathbb{R}}_{+}^5$, which was defined by (4) in the previous section.

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

For the system (2), the disease-free equilibrium (DFE), denoted by ${\mathbf{E}}^0$, can be obtained by setting the coordinates of all infective classes to zero, that is, by substituting $L=I=T_1=T_2=0$ in (Section 3). As a result, we have

$${\mathbf{E}}^0=\left({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0\right),$$

which is nonnegative and belongs to the absorbing set $\Omega$.

One of the fundamental concepts in the epidemiology of infectious diseases is the basic reproductive number, $\mathcal{R}>0$, which quantifies the transmission potential of an infectious disease within a susceptible population. This metric represents the average number of secondary infections generated by a single infected individual in a completely susceptible population [13]. From a mathematical standpoint, the basic reproductive number $\mathcal{R}$ is usually calculated as the spectral radius of the next-generation matrix at the DFE [14]. According to this approach, we should first define a sub-vector $\mathbf{Y}\in {\mathbb{R}}_{+}^4$ containing the human compartments that carry the infection

$$\mathbf{Y}:=\left(L,I,T_1,T_2\right).$$

Next, we extract four differential equations from the dynamical system (2) that correspond to the components of $\mathbf{Y}$ and write them in the form

$${\displaystyle \frac{d\mathbf{Y}}{dt}}=\mathcal{F}\left(\mathbf{Y}\right)-\mathcal{V}\left(\mathbf{Y}\right),$$

where

$$\mathcal{F}\left(\mathbf{Y}\right):=\left[\begin{array}{c}(1-p)\beta S(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2)\\ p(\beta S(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2)\\ 0\\ 0\end{array}\right],\mathcal{V}\left(\mathbf{Y}\right):=\left[\begin{array}{c}(ϵ+\mu )L-{\delta }_1{T}_1\\ (\gamma +\alpha +\mu )I-({\kappa }_1{T}_1+{\kappa }_2{T}_2+ϵL)\\ ({\delta }_1+{\kappa }_1+\mu )T_1-(1-q)\gamma I\\ ({\delta }_2+{\kappa }_2+\mu )T_2-q\gamma I\end{array}\right].$$

It is worth noting that $\mathcal{F}\left(\mathbf{Y}\right)\ge \mathbf{0}$ stands for the rate of appearance of new infections (i.e., the disease transmission) and $\mathcal{V}\left(\mathbf{Y}\right)\ge \mathbf{0}$ accounts for the disease transition rate. Both indicated inequalities are understood in the “component-by-component” sense.

According to [14], the next-generation matrix is defined by $\mathcal{F}{\mathcal{V}}^{-1}$ where $\mathcal{F}$ and $\mathcal{V}$ are the Jacobian matrices of $\mathcal{F}\left(\mathbf{Y}\right)$ and $\mathcal{V}\left(\mathbf{Y}\right)$ evaluated at the disease-free equilibrium ${\mathbf{E}}^0=\left({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0\right)$, that is,

$$\mathcal{F}:={\displaystyle \frac{\partial \mathcal{F}}{\partial \mathbf{Y}}}{|}_{{\mathbf{E}}^0}=\left[\begin{array}{cccc}0& (1-p)\beta {\displaystyle \frac{\Lambda }{\mu }}& (1-p)\beta {\varphi }_1{\displaystyle \frac{\Lambda }{\mu }}& (1-p)\beta {\varphi }_2{\displaystyle \frac{\Lambda }{\mu }}\\ 0& p\beta {\displaystyle \frac{\Lambda }{\mu }}& p\beta {\varphi }_1{\displaystyle \frac{\Lambda }{\mu }}& p\beta {\varphi }_2{\displaystyle \frac{\Lambda }{\mu }}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],$$

$$\mathcal{V}:={\displaystyle \frac{\partial \mathcal{V}}{\partial \mathbf{Y}}}{|}_{{\mathbf{E}}^0}=\left[\begin{array}{cccc}ϵ+\mu & 0& -{\delta }_1& 0\\ -ϵ& \gamma +\alpha +\mu & -{\kappa }_1& -{\kappa }_2\\ 0& -(1-q)\gamma & {\delta }_1+{\kappa }_1+\mu & 0\\ 0& -q\gamma & 0& {\delta }_2+{\kappa }_2+\mu \end{array}\right].$$

To simplify our calculations, we introduce the following notations for the four positive elements located on the main diagonal of $\mathcal{V}$

$$A:=ϵ+\mu ,B:=\gamma +\alpha +\mu ,C:={\delta }_1+{\kappa }_1+\mu ,D:={\delta }_2+{\kappa }_2+\mu$$

(6)

so the determinant of $\mathcal{V}$ becomes

$$\begin{array}{cc}\hfill det\mathcal{V}& =Adet\left|\begin{array}{ccc}B& -{\kappa }_1& -{\kappa }_2\\ -(1-q)\gamma & C& 0\\ -q\gamma & 0& D\end{array}\right|+ϵdet\left|\begin{array}{ccc}0& -{\delta }_1& 0\\ -(1-q)\gamma & C& 0\\ -q\gamma & 0& D\end{array}\right|\hfill \\ \hfill & =ABCD-(1-q)\gamma \left({\kappa }_{1}A+{\delta }_1ϵ\right)D-q\gamma {\kappa }_{2}AC\hfill \end{array}$$

(7)

The above form of $det\mathcal{V}$ will be used in Section 3.3 to simplify some computations in the sequel. However, one can also obtain an equivalent form of (7) that reveals that $det\mathcal{V}>0$, and, therefore, is helpful for further analysis. Indeed,

$$\begin{array}{cc}\hfill det\mathcal{V}& =ABCD-(1-q)\gamma {\kappa }_{1}AD-q\gamma {\kappa }_{2}AC-(1-q)\gamma ϵ{\delta }_{1}D-(1-q)\gamma {\delta }_1\mu D+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =ABCD-(1-q)\gamma {\kappa }_{1}AD-q\gamma {\kappa }_{2}AC-(1-q)\gamma {\delta }_{1}AD+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =A\left[BCD-(1-q)\gamma {\kappa }_{1}D-q\gamma {\kappa }_{2}C-(1-q)\gamma {\delta }_{1}D\right]+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =A\left[CD(\gamma +\alpha +\mu )-(1-q)\gamma ({\kappa }_1+{\delta }_1)D-q\gamma {\kappa }_{2}C\right]+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =A\left[CD(\alpha +\mu )+\gamma CD-(1-q)\gamma ({\kappa }_1+{\delta }_1)D-q\gamma {\kappa }_{2}C\right]+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =A\left[CD(\alpha +\mu )+\gamma \left[CD-(1-q)({\kappa }_1+{\delta }_1)D-q{\kappa }_{2}C\right]\right]+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =A\left[CD(\alpha +\mu )+\gamma \left[CD-(1-q)({\kappa }_1+{\delta }_1+\mu )D+(1-q)\mu D-q{\kappa }_{2}C\right]\right]+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =A\left[CD(\alpha +\mu )+\gamma \left[CD-CD+qCD+(1-q)\mu D-q{\kappa }_{2}C\right]\right]+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =A\left[CD(\alpha +\mu )+\gamma \left[q{\delta }_{2}C+q{\kappa }_{2}C+q\mu C+(1-q)\mu D-q{\kappa }_{2}C\right]\right]+(1-q)\gamma {\delta }_1\mu D\hfill \\ \hfill & =(ϵ+\mu )\left[CD(\alpha +\mu )+(1-q)\gamma \mu D+q\gamma ({\delta }_2+\mu )C\right]+(1-q)\gamma {\delta }_1\mu D>0.\hfill \end{array}$$

Thus, let us introduce the following notation

$$E:=det\mathcal{V}>0$$

(8)

where $det\mathcal{V}$ admits either the form (7) or

$$det\mathcal{V}=(ϵ+\mu )\left[CD(\alpha +\mu )+(1-q)\gamma \mu D+q\gamma ({\delta }_2+\mu )C\right]+(1-q)\gamma {\delta }_1\mu D,$$

(9)

as derived above. To write the next-generation matrix in a more compact form, let us introduce the following auxiliary quantities:

$$\begin{array}{cccc}\hfill \mathcal{A}:=& C(D+q{\varphi }_2\gamma )+(1-q){\varphi }_1\gamma D,\hfill & \hfill \mathcal{B}:=& D(AB{\varphi }_1+A{\kappa }_1+{\delta }_1ϵ),\hfill \\ \hfill \mathcal{C}:=& A{\varphi }_2{\kappa }_1-A{\varphi }_1{\kappa }_2+{\delta }_1ϵ{\varphi }_2,\hfill & \hfill \mathcal{D}:=& AC(B{\varphi }_2+{\kappa }_2).\hfill \end{array}$$

(10)

Then, we can write the next-generation matrix $\mathcal{F}{\mathcal{V}}^{-1}$ as

$$\mathcal{F}{\mathcal{V}}^{-1}={\displaystyle \frac{\beta {\displaystyle \frac{\Lambda }{\mu }}}{E}}\left[\begin{array}{cccc}(1-p)ϵ\mathcal{A}& (1-p)A\mathcal{A}& (1-p)[\mathcal{B}+q\gamma \mathcal{C}]& (1-p)[\mathcal{D}-(1-q\left)\gamma \mathcal{C}\right]\\ pϵ\mathcal{A}& pA\mathcal{A}& p[\mathcal{B}+q\gamma \mathcal{C}]& p[\mathcal{D}-(1-q\left)\gamma \mathcal{C}\right]\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].$$

The characteristic polynomial $\mathcal{P}\left(\lambda \right)$ corresponding to $\mathcal{F}{\mathcal{V}}^{-1}$ is

$$\mathcal{P}\left(\lambda \right)={\lambda }^4-{\displaystyle \frac{\beta }{E}}{\displaystyle \frac{\Lambda }{\mu }}\mathcal{A}\left[pA+(1-p)ϵ\right]{\lambda }^3$$

and using the left-hand relationship in (10), we have

$$\mathcal{P}\left(\lambda \right)={\lambda }^4-{\displaystyle \frac{\beta }{E}}{\displaystyle \frac{\Lambda }{\mu }}\left[CD+(1-q){\varphi }_1\gamma D+q{\varphi }_2\gamma C\right]\left[pA+(1-p)ϵ\right]{\lambda }^3.$$

According to the structure of $\mathcal{P}\left(\lambda \right)$, three eigenvalues of $\mathcal{F}{\mathcal{V}}^{-1}$ are zero, and the positive one is precisely equal to the spectral radius of $\mathcal{F}{\mathcal{V}}^{-1}$, that is

$$\varrho \left(\mathcal{F}{\mathcal{V}}^{-1}\right)={\displaystyle \frac{\beta }{E}}{\displaystyle \frac{\Lambda }{\mu }}\left[CD+(1-q){\varphi }_1\gamma D+q{\varphi }_2\gamma C\right]\left[pA+(1-p)ϵ\right].$$

(11)

Thus, we arrive at the basic reproductive number, ${\mathcal{R}}_0>0,$ for the TB transmission model (2) that expresses a mean number of secondary TB infections produced by one infectious individual when introduced into a totally susceptible population:

$${\mathcal{R}}_0={\displaystyle \frac{\beta {\displaystyle \frac{\Lambda }{\mu }}\left[CD+(1-q){\varphi }_1\gamma D+q{\varphi }_2\gamma C\right]\left[p(ϵ+\mu )+(1-p)ϵ\right]}{(ϵ+\mu )\left[CD(\alpha +\mu )+(1-q)\gamma \mu D+q\gamma ({\delta }_2+\mu )C\right]+(1-q){\delta }_1\gamma \mu D}}$$

(12)

Note that, using the notations of $C,D,$ and E defined in (6), (8), the full form of ${\mathcal{R}}_0$ is

$${\mathcal{R}}_0={\displaystyle \frac{\beta {\displaystyle \frac{\Lambda }{\mu }}\left[\left({\delta }_1+{\kappa }_1+\mu \right)\left({\delta }_2+{\kappa }_2+\mu \right)+(1-q){\varphi }_1\gamma \left({\delta }_2+{\kappa }_2+\mu \right)+q{\varphi }_2\gamma \left({\delta }_1+{\kappa }_1+\mu \right)\right]\left[p(ϵ+\mu )+(1-p)ϵ\right]}{(ϵ+\mu )\left[\left({\delta }_1+{\kappa }_1+\mu \right)\left({\delta }_2+{\kappa }_2+\mu \right)(\alpha +\mu )+(1-q)\gamma \mu \left({\delta }_2+{\kappa }_2+\mu \right)+q\gamma ({\delta }_2+\mu )\left({\delta }_1+{\kappa }_1+\mu \right)\right]+(1-q){\delta }_1\gamma \mu \left({\delta }_2+{\kappa }_2+\mu \right)}}$$

3.2. Existence of the Endemic Equilibrium

To obtain another feasible solution $(S^{*},L^{*},I^{*},T_1^{*},T_2^{*})$ of the nonlinear system (Section 3), we start from the equations (5a), (5d), and (5e) and express

$$\begin{array}{cc}\hfill T_1^{*}& =\left({\displaystyle \frac{(1-q)\gamma }{{\delta }_1+k_1+\mu }}\right)I^{*},\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill T_2^{*}& =\left({\displaystyle \frac{q\gamma }{{\delta }_2+k_2+\mu }}\right)I^{*},\hfill \end{array}$$

(13b)

$$\begin{array}{cc}\hfill S^{*}& ={\displaystyle \frac{\Lambda +{\delta }_2{T}_2^{*}}{\beta (I^{*}+{\varphi }_1{T}_1^{*}+{\varphi }_2{T}_2^{*})+\mu }}.\hfill \end{array}$$

(13c)

Then, by substituting (13a) and (13b) into (5b) we can also express

$$L^{*}=\left[{\displaystyle \frac{(1-p)\beta S^{*}}{ϵ+\mu }}+{\displaystyle \frac{(1-p)\beta {\varphi }_1(1-q)\gamma S^{*}}{(ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )}}+{\displaystyle \frac{(1-p)\beta {\varphi }_{2}q\gamma S^{*}}{(ϵ+\mu )({\delta }_2+{\kappa }_2+\mu )}}+{\displaystyle \frac{(1-q)\gamma {\delta }_1}{(ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )}}\right]I^{*}.$$

(13d)

Now, we substitute the quantities $T_1^{*}$, $T_2^{*}$, and $L^{*}$ given by (13a), (13b), and (13d) into equation (5c) and get

$$\begin{array}{cc}\hfill & p\beta S^{*}\left[1+{\displaystyle \frac{(1-q)\gamma {\varphi }_1}{{\delta }_1+{\kappa }_1+\mu }}+{\displaystyle \frac{q\gamma {\varphi }_2}{{\delta }_2+{\kappa }_2+\mu }}\right]I^{*}+{\displaystyle \frac{(1-q)\gamma {\kappa }_1}{{\delta }_1+{\kappa }_1+\mu }}{I}^{*}+{\displaystyle \frac{q\gamma {\kappa }_2}{{\delta }_2+{\kappa }_2+\mu }}{I}^{*}\hfill \\ \hfill & +ϵ\left[{\displaystyle \frac{(1-p)\beta S^{*}}{ϵ+\mu }}+{\displaystyle \frac{(1-p)\beta {\varphi }_1(1-q)\gamma S^{*}}{(ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )}}+{\displaystyle \frac{(1-p)\beta {\varphi }_{2}q\gamma S^{*}}{(ϵ+\mu )({\delta }_2+{\kappa }_2+\mu )}}+{\displaystyle \frac{(1-q)\gamma {\delta }_1}{(ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )}}\right]I^{*}\hfill \\ \hfill & -(\gamma +\alpha +\mu )I^{*}=0.\hfill \end{array}$$

An obvious solution of the above equation $I^{*}=0$ leads us to the disease-free equilibrium ${\mathbf{E}}^0=\left({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0\right)$, while its other possible solution(s) must satisfy the equality

$$\left[p\beta +{\displaystyle \frac{p\beta (1-q)\gamma {\varphi }_1}{{\delta }_1+{\kappa }_1+\mu }}+{\displaystyle \frac{p\beta q\gamma {\varphi }_2}{{\delta }_2+{\kappa }_2+\mu }}+{\displaystyle \frac{(1-p)\beta ϵ}{ϵ+\mu }}+{\displaystyle \frac{(1-p)\beta (1-q)\gamma {\varphi }_1ϵ}{(ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )}}+{\displaystyle \frac{(1-p)\beta q\gamma {\varphi }_2ϵ}{(ϵ+\mu )({\delta }_2+{\kappa }_2+\mu )}}\right]S^{*}$$

$$+\left[{\displaystyle \frac{(1-q)\gamma {\kappa }_1}{{\delta }_1+{\kappa }_1+\mu }}+{\displaystyle \frac{q\gamma {\kappa }_2}{{\delta }_2+{\kappa }_2+\mu }}+{\displaystyle \frac{(1-q)\gamma {\delta }_1ϵ}{(ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )}}-(\gamma +\alpha +\mu )\right]=0.$$

(14)

Note that the above equation is linear with respect to $S^{*}$. Using the definitions of ${\mathcal{R}}_0$ and E given in (12) and (8), respectively, we can express its independent term as

$$-{\displaystyle \frac{E}{(ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )({\delta }_2+{\kappa }_2+\mu )}}$$

and its leading coefficient as

$${\displaystyle \frac{E{\mathcal{R}}_0\mu }{\Lambda (ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )({\delta }_2+{\kappa }_2+\mu )}}.$$

Thus, (14) can be written in the following form:

$${\displaystyle \frac{E{\mathcal{R}}_0\mu }{\Lambda (ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )({\delta }_2+{\kappa }_2+\mu )}}{S}^{*}-{\displaystyle \frac{E}{(ϵ+\mu )({\delta }_1+{\kappa }_1+\mu )({\delta }_2+{\kappa }_2+\mu )}}=0$$

meaning that a unique solution of the linear equation (14) is

$$S^{*}={\displaystyle \frac{\Lambda }{\mu {\mathcal{R}}_0}}>0.$$

(15)

Notably, when ${\mathcal{R}}_0>1$, we have that $S^{*}<{\displaystyle \frac{\Lambda }{\mu }}$. Comparing the expression for $S^{*}$ given by (15) with the one in formula (13c), we can write

$${\displaystyle \frac{\Lambda }{\mu {\mathcal{R}}_0}}={\displaystyle \frac{\Lambda +{\delta }_2{T}_2^{*}}{\beta \left(I^{*}+{\varphi }_1{T}_1^{*}+{\varphi }_2{T}_2^{*}\right)+\mu }}\iff \Lambda \mu {\mathcal{R}}_0+{\delta }_2{T}_2^{*}\mu {\mathcal{R}}_0=\beta \Lambda I^{*}+\beta \Lambda {\varphi }_1{T}_1^{*}+\beta \Lambda {\varphi }_2{T}_2^{*}+\mu \Lambda .$$

Rearranging some terms in the last expression and substituting $T_1^{*}$ and $T_2^{*}$ from (13a) and (13b), we then obtain

$$\Lambda \mu ({\mathcal{R}}_0-1)=\beta \Lambda I^{*}+{\displaystyle \frac{\beta \Lambda {\varphi }_1(1-q)\gamma I^{*}}{{\delta }_1+{\kappa }_1+\mu }}+{\displaystyle \frac{\beta \Lambda {\varphi }_{2}q\gamma I^{*}}{{\delta }_2+{\kappa }_2+\mu }}-{\displaystyle \frac{{\delta }_{2}q\gamma \mu {\mathcal{R}}_0{I}^{*}}{{\delta }_2+{\kappa }_2+\mu }}$$

Again, using the definitions of ${\mathcal{R}}_0$ and E given in (12) and (8), we arrive at

$$\Lambda \mu ({\mathcal{R}}_0-1)={\displaystyle \frac{\mu {\mathcal{R}}_0\left[E-q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )\right]I^{*}}{(ϵ+p\mu )\left({\delta }_1+{\kappa }_1+\mu \right)\left({\delta }_2+{\kappa }_2+\mu \right)}},$$

from where the coordinate $I^{*}$ is finally expressed as

$$I^{*}={\displaystyle \frac{\Lambda \left({\delta }_1+{\kappa }_1+\mu \right)\left({\delta }_2+{\kappa }_2+\mu \right)(ϵ+p\mu )({\mathcal{R}}_0-1)}{{\mathcal{R}}_0\left[E-q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )\right]}}.$$

(16)

To establish the sign of the denominator in (16), we recall that $p\in [0,1]$ and $\mu >0$. Hence,

$$p\mu \le \mu \Rightarrow ϵ+p\mu \le (ϵ+\mu )\Rightarrow q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )\le q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+\mu )$$

because the parameters $ϵ,q,\gamma ,{\delta }_i,{\kappa }_i,i=1,2$ are nonnegative. Furthermore, using the definition of E given by (8) we have

$$\begin{array}{cc}\hfill q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )& <q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+\mu )+(1-q)\gamma {\delta }_1\mu \left({\delta }_2+{\kappa }_2+\mu \right)\hfill \\ \hfill & +(ϵ+\mu )\left[\left({\delta }_1+{\kappa }_1+\mu \right)\left({\delta }_2+{\kappa }_2+\mu \right)(\alpha +\mu )+(1-q)\gamma \mu \left({\delta }_2+{\kappa }_2+\mu \right)+q\gamma \mu \left({\delta }_1+{\kappa }_1+\mu \right)\right]=E.\hfill \end{array}$$

Thus, it holds that

$$E-q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )>0,$$

meaning that $I^{*}$ defined by formula (16) is positive only if ${\mathcal{R}}_0>1.$

From the preceding rationale and to culminate this section, we formulate the following statement concerned with the existence of equilibria of the TB transmission system given by (2).

Proposition 2. 

If ${\mathcal{R}}_0<1$, the system (2) has only the disease-free equilibrium ${\mathbf{E}}^0=\left({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0\right)$. If ${\mathcal{R}}_0>1$, the system (2) has two equilibria: the disease-free equilibrium $E^0$ and the endemic (strictly positive) equilibrium ${\mathbf{E}}^{*}:=\left(S^{*},L^{*},I^{*},T_1^{*},T_2^{*}\right)$ whose coordinates are defined by

$$\begin{array}{cc}\hfill S^{*}& ={\displaystyle \frac{\Lambda }{\mu {\mathcal{R}}_0}},\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill L^{*}& ={\displaystyle \frac{\Lambda \left[(1-p)E+(1-q)\gamma {\delta }_1(ϵ+p\mu )\left({\delta }_2+{\kappa }_2+\mu \right)\right]({\mathcal{R}}_0-1)}{{\mathcal{R}}_0(ϵ+\mu )\left[E-q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )\right]}},\hfill \end{array}$$

(17b)

$$\begin{array}{cc}\hfill I^{*}& ={\displaystyle \frac{\Lambda \left({\delta }_1+{\kappa }_1+\mu \right)\left({\delta }_2+{\kappa }_2+\mu \right)(ϵ+p\mu )({\mathcal{R}}_0-1)}{{\mathcal{R}}_0\left[E-q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )\right]}},\hfill \end{array}$$

(17c)

$$\begin{array}{cc}\hfill T_1^{*}& ={\displaystyle \frac{\Lambda (1-q)\gamma \left({\delta }_2+{\kappa }_2+\mu \right)(ϵ+p\mu )({\mathcal{R}}_0-1)}{{\mathcal{R}}_0\left[E-q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )\right]}},\hfill \end{array}$$

(17d)

$$\begin{array}{cc}\hfill T_2^{*}& ={\displaystyle \frac{\Lambda q\gamma \left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )({\mathcal{R}}_0-1)}{{\mathcal{R}}_0\left[E-q\gamma {\delta }_2\left({\delta }_1+{\kappa }_1+\mu \right)(ϵ+p\mu )\right]}},\hfill \end{array}$$

(17e)

where E has the form given in (7)-(8). Furthermore, the DFE is a boundary equilibrium (that is, $E^0\in \partial \Omega$), while the endemic equilibrium is interior (${\mathbf{E}}^{*}\in \mathrm{int}\Omega$) whenever $\alpha >0$, and becomes a boundary equilibrium only if $\alpha =0$.

A formal proof of Proposition 2 is provided in Appendix A.

3.3. Long-Term Behavior of the TB Transmission Dynamics

According to Proposition 2, the TB transmission system may possess either one or two equilibria, whose existence is determined by the threshold value of basic reproductive number ${\mathcal{R}}_0$ defined by (12), which also plays a crucial role in the long-term behavior of the solutions of the system (2).

It is worth recalling that if ${\mathcal{R}}_0<1$, then each contagious individual induces less than one new infection, on average, during their infectious period, meaning that the disease will disappear sooner or later regardless of its current state. From the mathematical perspective, it can be expected that the trajectories of the system (2) would eventually converge to ${\mathbf{E}}^0=\left({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0\right)$ when $t\to \infty$ as ${\mathbf{E}}^0$ is the only equilibrium that exists when ${\mathcal{R}}_0<1$.

On the other hand, if ${\mathcal{R}}_0>1$, then each contagious individual induces more than one new infection, on average, during their infectious period, meaning that the disease will persist because new infections will consequently arise. In other words, the trajectories of (2) will converge to the endemic equilibrium ${\mathbf{E}}^{*}$ with strictly positive coordinates defined by (2), thus diverging from the DFE, ${\mathbf{E}}^0$.

To formalize this intuitive rationale, we formulate the following result.

Theorem 1. 

The long-term evolution of the system (2) depends on the threshold value of the basic reproductive number ${\mathcal{R}}_0$ defined by (12) in the following sense:

(i)

When ${\mathcal{R}}_0<1$, the disease-free equilibrium ${\mathbf{E}}^0=\left({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0\right)$ is globally asymptotically stable (or GAS) on the absorbing set $\Omega \subset {\mathbb{R}}_{+}^5$ given by (4).

(ii)

When ${\mathcal{R}}_0>1$, the disease-free equilibrium ${\mathbf{E}}^0$ is unstable, and the endemic equilibrium ${\mathbf{E}}^{*}$ with strictly positive coordinates defined by (2) is GAS on the open subset ${\Omega }_0\subset \Omega$ defined by

$${\Omega }_0:=\left\{\left(S,L,I,T_1,T_2\right)\in \Omega :L+I+T_1+T_2>0\right\}$$

(18)

belonging to the absorbing set Ω.

A formal proof of Theorem 1 is provided in Appendix A.

From the common-sense standpoint, allowing for any equilibrium in which the disease will always persist is undesirable. Therefore, instead of paying too much attention to the stability properties of ${\mathbf{E}}^{*}$, we would like to focus more on achieving the disease-free state ${\mathbf{E}}^0$ by adequately changing the parameter $q\in [0,1]$ that mimics the assignment of actively infected people for sterilizing treatment.

In the next section, we propose employing the optimal control framework to obtain an optimal dynamic assignment for sterilizing treatment. We seek to minimize active and latent human TB infections while keeping the total costs of both treatment types (non-sterilizing and sterilizing) within a determined budget.

4. Optimal Control Approach

The share $q\in [0,1]$ of infectious patients assigned for sterilizing treatment was assumed constant in the previous sections. However, this share can also be modeled by a piecewise continuous real function $u\left(t\right):{\mathbb{R}}_{+}\to [0,1]$ dependent on time $t\ge 0$. Thus, we introduce a control function $u\left(t\right)$ that mimics the dynamic (time-variable) share of infectious patients who will further undergo the sterilizing treatment. In contrast, $\left(1-u\left(t\right)\right)$ will model the time-variable share of infectious patients who will further undergo non-sterilizing treatment. As both $u\left(t\right)$ and $\left(1-u\left(t\right)\right)$ are piecewise continuous and bounded to the range $[0,1]$, replacement of q by $u\left(t\right)$ in the dynamical system (2) will not alter the existence and uniqueness of its solution, and $\Omega \subset {\mathbb{R}}_{+}^5$ will remain being the absorbing set of the controlled system

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& =\Lambda -\beta S\left(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2\right)+{\delta }_2{T}_2-\mu S,\hfill \end{array}$$

(19a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dL}{dt}}& =(1-p)\beta S\left(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2\right)+{\delta }_1{T}_1-ϵL-\mu L,\hfill \end{array}$$

(19b)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dI}{dt}}& =p\beta S\left(I+{\varphi }_1{T}_1+{\varphi }_2{T}_2\right)+ϵL+{\kappa }_1{T}_1+{\kappa }_2{T}_2-\gamma I-\alpha I-\mu I,\hfill \end{array}$$

(19c)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{T}_1}{dt}}& =\left(1-u\left(t\right)\right)\gamma I-{\delta }_1{T}_1-{\kappa }_1{T}_1-\mu T_1,\hfill \end{array}$$

(19d)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{T}_2}{dt}}& =u\left(t\right)\gamma I-{\delta }_2{T}_2-{\kappa }_2{T}_2-\mu T_2,\hfill \end{array}$$

(19e)

for all admissible piecewise continuous real functions $u\left(t\right)\in [0,1]$ defined for $t\ge 0.$

The choice of the control function $u\left(t\right)$ can pursue one or more desirable (but also realistic) goals that will benefit public health. In particular, the costs of sterilizing and non-sterilizing treatment may be essential in choosing $u\left(t\right)$. In this context, it is worth noting that sterilizing TB treatment is usually more expensive than non-sterilizing treatment, even though the total duration of both treatment types is basically the same [15,16].

To model the introduction of sterilizing TB treatment, we assume that the observation period starts at $t=0$ and ends at $t=T$ where T is fixed. We also suppose that before the observation period $[0,T]$, all actively infected people have been receiving only non-sterilizing treatment. This leads us to assign the following initial conditions to the controlled system (Section 4)

$$S\left(0\right)=S_0,L\left(0\right)=L_0,I\left(0\right)=I_0,T_1\left(0\right)=T_0^1,T_2\left(0\right)=0,$$

(20a)

such that

$$N\left(0\right)=S\left(0\right)+L\left(0\right)+I\left(0\right)+T_1\left(0\right)+T_2\left(0\right)\le \Lambda /\mu$$

(20b)

Here, it is essential to point out that all people undergoing treatment at $t=0$ should continue receiving non-sterilizing drugs, meaning that they may either become latently infected (L-class) after completing the treatment or become actively infected (I-class) in the case of treatment abandonment.

In this section, we propose modeling the introduction of sterilizing TB treatment, while assuming that the additional costs related to the use of sterilizing drugs must be covered by the budget allocated for the habitual use of non-sterilizing drugs (budget constraints). This approach appears more realistic and practical in application, even though it requires solving an optimal control problem with a budget constraint.

Let $u\left(t\right):[0,T]\to [0,1]$ denote a piecewise continuous control function and $[0,T]$ be a fixed observation interval. The purpose of control consists in minimizing the number of actively and latently infected people during the observation period, expressed by the objective functional:

$$\underset{\begin{array}{c}0\le u\left(t\right)\le 1\end{array}}{min}J\left(u\right)=\underset{\begin{array}{c}0\le u\left(t\right)\le 1\end{array}}{min}{\displaystyle \underset{0}{\overset{T}{\int }}}\left[C_{1}I\left(t\right)+L\left(t\right)\right]dt,$$

(21)

where $C_1>1$ defines the priority of minimizing the current number of actively infected individuals over the latently infected ones because the former ones do transmit the disease while the latter ones do not.

To introduce the budget constraints, we have to model the costs related to each treatment type. Here, we propose using relative costs and thus avoid dealing with absolute values of non-sterilizing and sterilizing drugs that differ from country to country.

First, we note that the treatment cost depends on the dynamic compartments $T_i\left(t\right),i=1,2$, which contain all people undergoing treatment. Suppose we have the time interval $[0,T]$ as a planning horizon and, for simplicity, assume that the cost (per individual) related to the non-sterilizing treatment of one patient entering the class $T_1\left(t\right)$ is normalized to unity. If only non-sterilizing treatment is assigned to all identified infectious individuals (meaning that $u\left(t\right)=0$ and $T_2\left(t\right)\equiv 0$ for all $t\in [0,T]$), then instantaneous cost is equal to $1·\gamma I\left(t\right)$, and the cumulative cost related solely to non-sterilizing treatment can be assessed by the integral

$$B:={\displaystyle \underset{0}{\overset{T}{\int }}}\gamma I\left(t\right)dt,$$

(22)

which will serve as the “baseline cumulative cost”. Notably, $I\left(t\right)$ in (22) is the solution to the system (Section 4) with $u\left(t\right)=0,t\in [0,T]$.

As sterilizing drugs are more expensive compared to non-sterilizing ones [15], we may also assume that the cost (per unit time) related to the sterilizing treatment is $C_2$-times more expensive with $C_2>1$ compared to non-sterilizing treatment. For example, $C_2=2$ implies that sterilizing drugs are twice as expensive as non-sterilizing drugs.

On the other hand, a recent WHO report [4] highlights the ongoing efforts to reduce the cost of sterilizing TB drugs, emphasizing partnerships like the Global Drug Facility to negotiate lower prices and improve affordability in low- and middle-income countries. Currently, the key sterilizing TB drugs (such as pyrazinamide and rifampicin) remain rather expensive, primarily due to the patent protections and limited generic production. However, advocacy campaigns are pressuring manufacturers to reduce prices and increase the accessibility of sterilizing TB drugs by increasing their generic production [17]. All these actions raise hope that the prices of sterilizing TB drugs may gradually go down in the future. To model such a situation, we can define the relative cost (per unit time) of sterilizing treatment using a decreasing function depending on time (see Figure 2), that is,

$$c\left(t\right):=C_2-{\displaystyle \frac{C_3}{T}}·t\mathrm{where}{C}_2>1,C_3\ge 0,\mathrm{and}{C}_2-C_3\ge 1.$$

(23)

This function $c\left(t\right)$ can be used to define two scenarios related to the relative cost of sterilizing treatment, namely:

Constant cost. 

In this case, we assume $C_3=0$ and $c\left(t\right)=C_2>1$, meaning that the relative (per unit time) cost of sterilizing treatment remains constant during the entire observation period $[0,T]$.

Decreasing cost. 

In this case, we assume $C_3>0$ and $C_2-C_3\ge 1$ in (23), meaning that the relative (per unit time) cost of sterilizing treatment $c\left(t\right)\ge 1$ is a monotone decreasing function of time $t\in [0,T]$.

At the beginning of this section, we introduced a piecewise continuous control function $u\left(t\right)$ that mimics the dynamic (time-variable) share of infectious patients who will further undergo the sterilizing treatment, while $\left(1-u\left(t\right)\right)$ models the time-variable share of infectious patients who will further undergo non-sterilizing treatment. Therefore, the instantaneous cost related to both types of treatment is

$$\cos \mathrm{t}\left(t\right)=\left(1-u\left(t\right)\right)\gamma I\left(t\right)+c\left(t\right)u\left(t\right)\gamma I\left(t\right)=\left[\left(c\left(t\right)-1\right)u\left(t\right)+1\right]\gamma I\left(t\right)$$

To assess the cumulative total cost related to both treatments for a given $u\left(t\right):[0,T]\to [0,1]$, we employ the following integral

$${\displaystyle \underset{0}{\overset{T}{\int }}}\left[\left(c\left(t\right)-1\right)u\left(t\right)+1\right]\gamma I\left(t\right)dt,$$

(24)

where $I\left(t\right)$ is the solution to the system (Section 4)-(Section 4) for a given $u\left(t\right)$, while $c\left(t\right)\ge 1,t\in [0,T]$ is defined by (23). Here, we also recall that the cost (per individual) related to the non-sterilizing treatment is normalized to unity.

Using the baseline cumulative cost B defined by (22), we can now formulate the budget constraint expressing that the additional expenses related to the use of sterilizing drugs must come from the budget calculated for the habitual use of non-sterilizing drugs

$${\displaystyle \underset{0}{\overset{T}{\int }}}\left[\left(C_2-{\displaystyle \frac{C_3}{T}}·t-1\right)u\left(t\right)+1\right]\gamma I\left(t\right)dt\le B,$$

(25)

where the constant $C_3$ can be either zero or positive and thus be adjusted to either constant or decreasing costs described above. The integral constraint defined by (25) (known as an isoperimetric constraint in the optimal control theory) must be satisfied alongside the dynamics of the controlled system (Section 4)-(Section 4). Notably, the constraint (25), when added to the optimal control problem, will narrow down the class of admissible piecewise continuous control functions $u\left(t\right):[0,T]\to [0,1]$ and their corresponding states $\mathbf{Z}\left(t\right)=\mathbf{Z}\left(t;u\left(t\right)\right)$, which are solutions to the dynamical system (Section 4)-(Section 4). However, this constraint is feasible in the sense that there exists at least one control function (namely, $u\left(t\right)=0$ for all $t\in [0,T]$) that fulfills it (cf. relationship (22)).

Taking into account the constraint (25), the underlying optimal control problem can be formulated as follows. For each $t\in [0,T]$, find an optimal share $u^{*}\left(t\right)\in [0,1]$ of individuals from the class $I\left(t\right)$ to undergo a sterilizing TB treatment to minimize the objective functional (21) subject to the dynamical system (Section 4)-(Section 4) together with the constraint (25).

Dealing with the budget constraint (25) is challenging. However, as stipulated in [18], an inequality-type constraint can be “merged” into the objective using a Lagrange multiplier $\eta \ge 0$ thus forming the “augmented objective” (or Lagrangian functional), that is,

$$\begin{array}{cc}\hfill \mathcal{L}\left(u\right)& =J\left(u\right)+\eta \left({\displaystyle \underset{0}{\overset{T}{\int }}}\gamma I\left(t\right)\left[\left(c\left(t\right)-1\right)u\left(t\right)+1\right]dt-B\right)\hfill \\ \hfill & ={\displaystyle \underset{0}{\overset{T}{\int }}}\left(C_{1}I\left(t\right)+L\left(t\right)+\eta \gamma I\left(t\right)\left[\left(C_2-{\displaystyle \frac{C_3}{T}}·t-1\right)u\left(t\right)+1\right]\right)dt-\eta B,\hfill \end{array}$$

(26)

where the constant term $-\eta B$ is of no concern for minimization. The multiplier $\eta \ge 0$ is chosen to satisfy the complementary slackness condition

$$\eta ·\left({\displaystyle \underset{0}{\overset{T}{\int }}}\gamma I\left(t\right)\left[\left(c\left(t\right)-1\right)u^{*}\left(t\right)+1\right]dt-B\right)=0.$$

(27)

at the optimal control $u^{*}\left(t\right)$. If the budget is inactive at the optimum (meaning the inequality (25) is strict in $u^{*}$), then $\eta =0$ and the problem reduces to the unconstrained minimization. However, if the budget is active, then $\eta >0$ is chosen so that the budget inequality (25) holds as equality at $u^{*}$.

Remark 1. 

In more intuitive terms, the multiplier $\eta >$ can also be viewed as a “penalty weight” in the augmented objective $\mathcal{L}\left(u\right)$ as its larger value will put more importance on satisfying the budget constraint (25) compared to minimizing the integrand of the original objective $J\left(u\right)$ [18].

Remark 2. 

If no budget constraint is imposed (meaning that (25) is ignored), we can deliberately set $\eta =0$ in the sequel to obtain the unconstrained optimal control problem where $\mathcal{L}\left(u\right)=J\left(u\right)$ is minimized on the solutions to the ODE system (Section 4).

To account for the budget constraint (25), we formulate the optimal control problem in a slightly different form. For each $t\in [0,T]$, find an optimal share $u^{*}\left(t\right)\in [0,1]$ of the actively infected and symptomatic individuals from the class $I\left(t\right)$ to undergo a sterilizing TB treatment to minimize the augmented objective functional (26) with $\eta \ge 0$ satisfying the complementary slackness condition (27), on the solutions to the ODE system (Section 4) with initial conditions (Section 4).

This problem can be formally solved by applying the Pontryagin minimum principle (see, e.g., [18,19] or other similar textbooks). However, one has to recall [18] that the necessary condition for the solvability of the optimal control problem (Section 4)-(21) is the solvability of the initial value problem (Section 4)-(Section 4) for any piecewise continuous control function $u\left(t\right)\in [0,1],t\in [0,T],$ together with the boundedness of all its possible trajectories $\mathbf{Z}\left(t\right):=\left(S\left(t\right),L\left(t\right),I\left(t\right),T_1\left(t\right),T_2\left(t\right)\right)$. To validate this necessary condition, we formulate the following result.

Proposition 3. 

For any measurable function $u\left(t\right):[0,T]\to [0,1]$, the initial value problem (Section 4)-(Section 4) has a unique continuous solution $\mathbf{Z}\left(t\right)=\mathbf{Z}\left(t;u\left(t\right)\right)$ bounded for all $t\in [0,T]$ in the sense that $\mathbf{Z}\left(t\right)\in \Omega .$

Proof. 

The initial value problem (Section 4)-(Section 4) can be written as

$${\displaystyle \frac{d\mathbf{Z}}{dt}}=\Phi \left(\mathbf{Z}\left(t\right),u\left(t\right)\right),\mathbf{Z}\left(0\right)=\left(S_0,L_0,I_0,T_0^1,0\right),$$

where $\Phi \left(\mathbf{Z}\left(t\right),u\left(t\right)\right)$ is the vector function of the right-hand sides of the controlled dynamical system (Section 4). This vector field is locally Lipschitz in $\mathbf{Z}$ and for any measurable $u\left(t\right)\in [0,1]$ defined for $t\in [0,T]$, and the Carathéodory existence and uniqueness theorem applies [20,21]. Therefore, a unique solution $\mathbf{Z}\left(t\right)=\mathbf{Z}\left(t;u\left(t\right)\right)$ exists for any initial condition $\mathbf{Z}\left(0\right)$ and for all measurable $u\left(t\right)\in [0,1],t\in [0,T]$. Furthermore, in light of Proposition 1 and the property (20b), all the solutions of the form $\mathbf{Z}\left(t\right)=\mathbf{Z}\left(t;u\left(t\right)\right)$ are bounded and belong to the closed set $\Omega \subset {\mathbb{R}}_{+}^5$ defined by (4). □

After establishing the solvability of the controlled dynamical system (Section 4)-(Section 4) for all measurable $u\left(t\right)\in [0,1]$ defined for $t\in [0,T]$, including piecewise continuous functions, we formulate the following statement to ensure the existence of an optimal control $u^{*}$ that minimizes the augmented objective functional (26) subject to the dynamical system (Section 4)-(Section 4).

Proposition 4. 

There exists a measurable optimal control function $u^{*}\left(t\right):[0,T]\to [0,1]$, and the corresponding solution ${\mathbf{Z}}^{*}\left(t\right)=\mathbf{Z}\left(t;u^{*}\left(t\right)\right)=\left(S^{*}\left(t\right),L^{*}\left(t\right),I^{*}\left(t\right),T_1^{*}\left(t\right),T_2^{*}\left(t\right)\right)$ to the initial value problem (Section 4)-(Section 4) that minimize the augmented objective functional $\mathcal{L}\left(u\right)$ defined by (26).

Proof. 

The existence of an optimal control $u^{*}\left(t\right)$ belonging to the class of measurable (that is, integrable) functions, can be proven by applying the standard result [22,23], according to which the fulfillment of certain conditions will be sufficient for the existence of an optimal control. In our case, the validity of the following conditions will suffice to guarantee the existence of an optimal control:

(i)

Boundedness of the state space. It was already proved by Proposition 3 that all the state trajectories $\mathbf{Z}\left(t;u\left(t\right)\right)$ of the controlled system (Section 4) remain bounded for any $u\left(t\right)\in [0,1],t\in [0,T]$. Moreover, the sets of initial condition $\mathbf{Z}\left(0\right)$ defined by (Section 4) and the set of the terminal conditions $\mathbf{Z}\left(T\right):=\left(S\left(T\right),L\left(T\right),I\left(T\right),T_1\left(T\right),T_2\left(T\right)\right)$ are bounded by the compact set $\Omega$.

(ii)

Compactness and convexity of the control set. Since $u\left(t\right)\in [0,1]$, the control set is bounded and convex.

(iii)

Linear dependence on control. The right-hand sides of the ODE system (Section 4) are linear in u.

(iv)

Continuity of the objective function. The integrand in (21) is continuous with respect to $L,I,u$, and t.

(v)

Lower semicontinuity of the objective functional. The value of $\mathcal{L}\left(u\right)$ is bounded by $-\eta B$ from below.

Finally, the absorbing set $\Omega$ defined by (4) is convex for any $u\left(t\right)\in [0,1]$ and for all $t\in [0,T]$, and the existence result from [22] (see Theorem 9.3.i) is applicable. The latter, together with the validity of conditions (i)-(v) ensures the existence of $u^{*}:[0,T]\to [0,1]$ that minimizes the augmented functional (26) over any finite time interval $t\in [0,T]$ subject to the dynamical system (Section 4) with initial conditions (Section 4).

It is also important to point out that the formulated optimal control problem with budget constraint is feasible in the sense that there is at least one control function ($u\left(t\right)=0,t\in [0,T]$) that satisfies the budget constraint (25). □

Once the existence of optimal control is proved by Proposition 4, we proceed to formulate the necessary condition of optimality in the form of the Pontryagin minimum principle (see, e.g., [18,19] or similar textbooks). Here, the Hamiltonian corresponding to the optimal control problem (Section 4)-(21) is

$$H\left(\mathbf{Z},u,\mathbf{\lambda },t\right)=C_{1}I+L+\eta \gamma I\left[\left(c\left(t\right)-1\right)u+1\right]+〈\mathbf{\lambda },\Phi (\mathbf{Z},u)〉,$$

(28)

where $\langle ·,·\rangle$ stands for scalar product, $\Phi (\mathbf{Z},u)$ is the vector field of the system (Section 4), and $\mathbf{\lambda }=\left({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5\right)\in {\mathbb{R}}^5$ denotes the vector of co-states (or adjoint variables) corresponding to the state variables $\mathbf{Z}=\left(S,L,I,T_1,T_2\right)=\left(Z_1,Z_2,Z_3,Z_4,Z_5\right)$. In the above expression, the co-state vector function $\mathbf{\lambda }=\mathbf{\lambda }\left(t\right)$ is absolutely continuous and satisfies the adjoint system with transversality condition

$${\displaystyle \frac{d\mathbf{\lambda }}{dt}}=-{\displaystyle \frac{\partial H\left(\mathbf{Z},u,\mathbf{\lambda },t\right)}{\partial \mathbf{Z}}},\mathbf{\lambda }\left(T\right)=\mathbf{0}.$$

(29)

Note that the detailed form of (29) is provided in Appendix B.

The optimal control $u^{*}\left(t\right)$ whose existence was proved by Proposition 4 must satisfy the Pontryagin minimum principle, which is a necessary condition of optimality. Therefore, almost at each $t\in [0,T]$ along the optimal path ${\mathbf{Z}}^{*}\left(t\right)=\mathbf{Z}\left(t;u^{*}\left(t\right)\right)$, the Hamiltonian (28) must fulfill the condition

$$H\left({\mathbf{Z}}^{*},u^{*},\mathbf{\lambda },t\right)\le H\left({\mathbf{Z}}^{*},u,\mathbf{\lambda },t\right)\mathrm{for}\mathrm{all}u\in [0,1].$$

Here, it is worth noting that Hamiltonian is linear in u, and its partial derivative

$${\displaystyle \frac{\partial H}{\partial u}}=\gamma I\left[\eta \left(c\left(t\right)-1\right)+{\lambda }_5-{\lambda }_4\right]$$

(30)

is independent of u, meaning that $H\left(\mathbf{Z},u,\mathbf{\lambda },t\right)$ attains its minimum with respect to u at the extreme points of the closed interval $[0,1]$. From this rationale, it stems that

$$u^{*}\left(t\right)=\left\{\begin{array}{ccc}0\hfill & \mathrm{if}& \phi \left(t\right)>0\\ 1\hfill & \mathrm{if}& \phi \left(t\right)<0\end{array}\right.$$

(31)

where

$$\phi \left(t\right):=\gamma I\left(t\right)\left[\eta \left(c\left(t\right)-1\right)+{\lambda }_5-{\lambda }_4\right]$$

(32)

denotes the switching function that represents the right-hand side of the relationship (30) evaluated at each $t\in [0,T].$ From the expression (32), we can see that the sign of the switching function $\phi \left(t\right)$ is defined by the relationship between the extra cost, $\left(c\left(t\right)-1\right)\ge 0,$ weighted by $\eta \ge 0$, and the adjoint functions ${\lambda }_4\left(t\right)$ and ${\lambda }_5\left(t\right)$ corresponding to both types of treatment.

Let us recall the meaning of the adjoint functions ${\lambda }_4\left(t\right)$ and ${\lambda }_5\left(t\right)$ associated with the state variables $T_1\left(t\right)$ and $T_2\left(t\right)$. In general terms, adjoint variables encode how changes in the final state ($t=T$) propagate backward in time to influence the optimal control trajectory and the value of the objective functional. As the adjoint variables are determined by backward integration (see Eq. (29) above), their values reflect how the future goals (e.g., minimizing the cumulative number of all active and latent TB infections and the compliance with the budget, both expressed via $\mathcal{L}\left(u\right)$ in (26)) influence the present-day decisions. Here, the decision is to distribute all infectious TB patients to undergo either sterilizing or non-sterilizing treatment. Thus, ${\lambda }_4\left(t\right)$ and ${\lambda }_5\left(t\right)$ quantify how increasing the number of individuals assigned for non-sterilizing and sterilizing treatment ($u\left(t\right)=0$ and $u\left(t\right)=1$, respectively, at each $t\in [0,T]$) will affect the cumulative number of infections while complying with the assigned budget.

Furthermore, each adjoint variable can be interpreted as a shadow price representing the opportunity value of increasing the corresponding state variable at a given time t by one unit (ceteris paribus). From this angle, the quantities ${\lambda }_4\left(t\right)$ and ${\lambda }_5\left(t\right)+\eta \left(c\left(t\right)-1\right)$ express the opportunity values of adding one individual to the class $T_1$ and $T_2$ on the overall objective. Revising the formulas (31) and (32), one may easily guess that, for each single $t\in [0,T]$, the relationship ${\lambda }_4\left(t\right)<{\lambda }_5\left(t\right)+\eta \left(c\left(t\right)-1\right)$ will lead to employing the non-sterilizing treatment, while ${\lambda }_4\left(t\right)>{\lambda }_5\left(t\right)+\eta \left(c\left(t\right)-1\right)$ will motivate for the sterilizing treatment. Here, the multiplier $\eta \ge 0$ also acts like a shadow price per unit of extra sterilizing-treatment cost, because it may shift the switching threshold.

Thus, the optimal control expressed by (31) is piecewise constant and switches between two regimes: $u^{*}\left(t\right)=1$ (only sterilizing treatment) and $u^{*}\left(t\right)=0$ (only non-sterilizing treatment). The number of switchings is defined by the number of roots of the switching function (32) and depends on the relative magnitude of ${\lambda }_4\left(t\right)$ and ${\lambda }_5\left(t\right)+\eta \left(c\left(t\right)-1\right)$ at each $t\in [0,T]$. This type of optimal control is referred to as “bang-bang” with the switchings occurring at the isolated points $t_s\in [0,T]$ when $\phi \left(t_s\right)=0$, or equivalently when ${\lambda }_4\left(t_s\right)={\lambda }_5\left(t_s\right)+\eta \left(c\left(t_s\right)-1\right)$. The latter implies that both treatment types render the same opportunity value at $t_s$ while these values change at $t=t_s\pm \Delta t$ for arbitrary small $\Delta t$.

Remark 3. 

From the theoretical standpoint, there may exist an interval $[t_1,t_2]\subseteq [0,T]$ with $t_1<t_2$ where it holds ${\lambda }_4\left(t\right)={\lambda }_5\left(t\right)+\eta \left(c\left(t_s\right)-1\right)$ for all $t\in [t_1,t_2]$ (meaning that $\phi \left(t\right)$ vanish within $[t_1,t_2]$). In such a case, the optimal control may admit interior values when $t\in [t_1,t_2]$, that is, $0<u^{*}\left(t\right)<1$, meaning that it is optimal to assign both types of TB treatment during the period $[t_1,t_2]$. Situations like this are less probable than the straightforward bang-bang control, especially in our case, because of the form of the switching function (32) that explicitly includes the non-increasing cost function $c\left(t\right)$ and the constant multiplier η. Still, they may eventually take place and are known in the literature as singular arcs $u^{*}\left(t\right)\in (0,1),t\in [t_1,t_2]$ (see more details in [24] and similar textbooks).

Given the nonlinearity of the optimal control problem (Section 4)-(26) and the complexity of the adjoint system (29) (see its detailed form in Appendix B), the optimal control $u^{*}\left(t\right)$ and corresponding optimal states ${\mathbf{Z}}^{*}\left(t\right)=\mathbf{Z}\left(t;u^{*}\left(t\right)\right)$ can only be found numerically. Nonetheless, the general structure of $u^{*}\left(t\right)$ formalized by (31) and predicted by the Pontryagin minimum principle will serve as a reference to validate its numerical solution, which will be provided in the next section.

5. Numerical Solutions and Discussion

Let us recall that, according to [4], a combination of four drugs (rifampicin, pyrazinamide, isoniazid, and ethambutol) is necessary for the standard strict-scheme course of sterilizing treatment. If one or more drugs become (temporarily or permanently) unavailable to a patient, the treatment may lose its sterilizing effect and become non-sterilizing. This is precisely what is observed in many low- and middle-income countries due to the deficit of rifampicin, the most powerful first-line anti-TB sterilizing drug. Its intermittent presence or absence in the proper treatment scheme usually leads to non-sterilizing treatment based on three other drugs (isoniazid, pyrazinamide, and ethambutol) that are more accessible and cheap.

From this rationale, we can safely assume that both treatment types (sterilizing and non-sterilizing) have about the same duration, usually about 6 months [4], and the treatment escape rate can also be assumed the same for both treatments as the patients usually are not aware of the treatment type they receive. These two assumptions will remain valid for all numerical simulations presented in this section, along with the values of the parameters of the dynamical systems (2) and (Section 4), which are displayed in Table 2.

In [8], it was shown that both treatment types reduce the infectiousness of patients with active TB forms with solely sterilizing treatment contributing a bit more to such reduction [3], we focus on the “pessimistic” scenario and assume ${\varphi }_1={\varphi }_2$ in numerical simulations. We also suppose that treatment is assigned immediately after clinical tests confirm the presence of active bacilli, so the time between the beginning of the infectiousness and the commencement of therapy, $1/\gamma$, is around two months.

Numerous studies have revealed that tuberculosis tends to have a notably higher incidence in bigger cities (1 million inhabitants or more), featuring social disparities. Using the values of $\Lambda$ and $\mu$ given in Table 2, we can model the total human population of about $2.4-2.5$ million inhabitants. Thus, we take

$$N_h:=S\left(0\right)+L\left(0\right)+I\left(0\right)+T_1\left(0\right)+T_2\left(0\right)=2.4\mathrm{million}\mathrm{people}$$

and ${\displaystyle \frac{\Lambda }{\mu }}$ close to $2.5$ million inhabitants. At the same time, the TB transmission rate, $\beta$, is chosen to fit the annual TB incidence rate of around 200 cases per $100000$ inhabitants. We have also assumed that initial conditions (Section 4) assigned to the ODE system (Section 4) are close to the endemic equilibrium (2) of the system (2) with $q=0$, namely:

$$S\left(0\right)=N_h-L\left(0\right)-I\left(0\right)-T_1\left(0\right)-T_2\left(0\right),$$

(33a)

$$L\left(0\right)=1.15L^{*},I\left(0\right)=1.15I{L}^{*},T_1\left(0\right)=1.15T_1^{*},T_2\left(0\right)=0.$$

(33b)

In the sequel, we provide numerical simulations when the share of actively TB-infected patients to undergo the sterilizing treatment is either constant, $q\in [0,1]$ or variable in time, $u\left(t\right)\in [0,1],t\in [0,T]$.

5.1. Useful Insights from ODE System (2) with Constant $q\in [0,1]$

To show off the expected effect of sterilizing TB treatment, we present the evolution of the basic reproductive number ${\mathcal{R}}_0$ as a function of the parameter $q\in [0,1]$ (the share of actively infected TB patients assigned for the sterilizing treatment), that is, ${\mathcal{R}}_0={\mathcal{R}}_0\left(q\right)$ using the formula (12). Here, we employed the values of parameters of the model (2) defined in Table 2. The function ${\mathcal{R}}_0\left(q\right),q\in [0,1]$ is displayed in Figure 3, where ${\mathcal{R}}_0\left(0\right)=1.085>1$ meaning that if only the non-sterilizing treatment ($q=0$) is given to all actively-infected TB patients, the disease will persist. However, if $10\%$ or more ($0.1\le q\le 1$) of such patients will undergo the sterilizing treatment, there is a perspective to eradicate the disease (${\mathcal{R}}_0\left(u\right)<1,q\in [0.1,1]$).

Notably, ${\mathcal{R}}_0$ is a decreasing function of $q\in [0,1]$ (see Figure 3). Therefore, the greater the share of actively TB-infected patients assigned for the sterilizing treatment, the sooner the disease eradication can be achieved. The latter highlights the benefits of sterilizing TB treatment.

To quantify the effect of variations in model parameters on disease transmission potential, we conducted a sensitivity analysis of ${\mathcal{R}}_0$ using the parameter ranges listed in Table 2. The method described in [33,34] was applied to estimate the first-order and total-order sensitivity indices of the basic reproduction number. The analysis was based on $34000$ simulations, accounting for the variability of all model parameters. The first-order indices represent the direct contribution of each parameter to the variance of ${\mathcal{R}}_0$, excluding interaction effects with other parameters. The results presented on the left-hand chart of Figure 4 reveal that only a few parameters exert a substantial direct influence on ${\mathcal{R}}_0$, as summarized below:

• $\Lambda$ (recruitment rate of individuals): $0.32$ (95% CI: $0.21$–$0.42$) exerts the strongest direct influence on ${\mathcal{R}}_0$.
• $\beta$ (transmission rate): $0.24$ (95% CI: $0.15$–$0.30$) plays a significant role in determining the transmission dynamics.
• q (fraction of patients receiving sterilizing treatment): $0.19$ (95% CI: $0.09$–$0.27$) has a notable effect, although with considerable uncertainty.
• $\gamma$ (rate of screening/recruitment for treatment): $0.06$ (95% CI: $0.03$–$0.09$) shows a moderate but non-negligible impact.
• $ϵ$ (rate of slow infection progression): $0.03$ (95% CI: $-0.002$–$0.06$) contributes weakly, with higher uncertainty.
• Other parameters ($\mu$, ${\varphi }_1$, ${\varphi }_2$, p, ${\delta }_1$, ${\delta }_2$, $\alpha$, ${\kappa }_1$, ${\kappa }_2$, $N_h$) exhibit negligible direct contributions.

The total-order indices reflect the contribution of each parameter, including its interactions with all other parameters. Unlike the first-order indices, which measure the effect of varying one parameter at a time while keeping the others fixed, the total-order indices capture the combined influence that emerges when all parameters vary simultaneously. This provides a more realistic assessment of parameter importance under joint uncertainty. As shown in Figure 4 (right chart), the parameters with the highest total-order sensitivity indices are detailed below:

• q (fraction of patients receiving sterilizing treatment): 0.40 (95% CI: 0.32 – 0.48) indicates that, beyond its direct effect, it also interacts significantly with other parameters.
• $\Lambda$: 0.37 (95% CI: 0.25 – 0.46) keeps its robust effect.
• $\beta$: 0.26 (95% CI: 0.16 – 0.36) maintains its key role in the sensitivity of ${\mathcal{R}}_0$.
• $ϵ$: 0.12 (95% CI: 0.05 – 0.20) gains importance when interactions are considered.
• The remaining parameters ($\gamma$, ${\delta }_1$, ${\delta }_2$, ${\varphi }_1$, ${\varphi }_2$, $\mu$, $\alpha$, ${\kappa }_1$, ${\kappa }_2$,$N_h$, p) show values close to zero in their total indices, indicating a marginal influence.

To summarize, the basic reproduction number ${\mathcal{R}}_0$ exhibits the highest sensitivity to the fraction q of infectious individuals assigned to sterilizing treatment when all the model’s parameters vary within their respective ranges. Namely, by increasing q, the value of ${\mathcal{R}}_0$ will decrease even if other parameters of the model are perturbed, and the situation shown in Figure 3 is expected. It implies that, from an epidemiological perspective, interventions aimed at increasing the sterilizing treatment coverage will have a significant impact on reducing ${\mathcal{R}}_0$ and mitigating the disease propagation.

It is also relevant to consider two borderline scenarios, namely:

1.

All actively TB-infected patients will undergo the non-sterilizing treatment ($q=0$ in system (2)).

2.

All actively TB-infected patients will undergo the sterilizing treatment ($q=1$ in system (2)).

These two borderline scenarios may give us the worst- and best-expected outcomes regarding disease propagation and persistence in the short-, medium-, and long-term perspectives.

Numerical simulations of the dynamical system (2) with constant $q=0$ (solely non-sterilizing treatment) and $q=1$ (solely sterilizing treatment) are presented in Figure 5, using blue and red curves, respectively. The left- and right-hand charts in Figure 5 correspond to the solutions $L\left(t\right)$ and $I\left(t\right)$ of the ODE system (2), respectively. Besides fully agreeing with the expected value of ${\mathcal{R}}_0$ (cf. Figure 3), the graphs in Figure 5 clearly reflect the benefits of sterilizing treatment ($q=1$) compared to a non-sterilizing one ($q=0$). In effect, the use of solely non-sterilizing treatment leads to the endemic persistence of the disease (see the blue-colored curves in Figure 5), while solely sterilizing treatment is capable of reducing rapidly the current number of latent and active TB infections even in a short-term perspective (cf. the difference between the blue- and red-colored curves in Figure 5 at $t=5$ years, marked by the dashed vertical lines). Moreover, the benefits of solely sterilizing treatment ($q=1$) become more visible in the medium- and long-term perspective if we compare the blue- and red-colored curves in Figure 5 at $t=15$ and $t=30$.

It is worthwhile to highlight that the numerical outcomes displayed in Figure 5 fully agree with the theoretical results presented in Section 3 and Section 3.3. When $q=0$, we have ${\mathcal{R}}_0>1$ (see Figure 3) and the trajectories of the ODE system (2) converge to the endemic equilibrium ${\mathbf{E}}^{*}$ defined by (2). In contrast, when $q=1$, we have ${\mathcal{R}}_0<1$ (see Figure 3) meaning that ${\mathbf{E}}^{*}$ does not exist (as affirmed in Proposition 2) and the trajectories of the ODE system (2) should converge to ${\mathbf{E}}^0=\left({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0\right)$, its unique equilibrium, when $t\to \infty .$ This conclusion is perfectly consistent with Theorem 1.

There are other important dynamic indicators that help visualize the potential outcomes of both TB treatment types. Among them, we choose for illustration the following:

• Cumulative number of latent TB infections, $C_L\left(t\right),t\in [0,T]$, arising during the period $[0,T]$ as a result of contagion or non-sterilizing treatment, is defined by

  $$C_L\left(t\right):={\displaystyle \underset{0}{\overset{t}{\int }}}\left[(1-p)\beta S\left(\tau \right)\left(I\left(\tau \right)+{\varphi }_1{T}_1\left(\tau \right)+{\varphi }_2{T}_2\left(\tau \right)\right)+{\delta }_1{T}_1\left(\tau \right)\right]d\tau ,$$

  (34)

  where $S,I,T_1,T_2$ are solutions to the system (2) with either $q=0$ or $q=1.$
• Cumulative active TB incidence, $C_I\left(t\right),t\in [0,T]$, which sums up all new active TB infections that occur during the period $[0,T]$, is defined by

  $$C_I\left(t\right):={\displaystyle \underset{0}{\overset{t}{\int }}}\left[p\beta S\left(\tau \right)\left(I\left(\tau \right)+{\varphi }_1{T}_1\left(\tau \right)+{\varphi }_2{T}_2\left(\tau \right)\right)+ϵL\left(\tau \right)\right]d\tau ,$$

  (35)

  where $S,L,I,T_1,T_2$ are solutions to the system (2) with either $q=0$ or $q=1.$
• Cumulative TB-induced mortality, $C_M\left(t\right)$ that sums up all TB-induced death cases that appear during the period $[0,T]$ is defined by

  $$C_M\left(t\right):={\displaystyle \underset{0}{\overset{t}{\int }}}\alpha I\left(\tau \right)d\tau$$

  (36)

  where I the component of the solution to the dynamical system (2) with either $q=0$ or $q=1$.

Figure 6 displays the graphs for $C_L\left(t\right),C_I\left(t\right),$ and $C_M\left(t\right),t\in [0,T]$. Even though the reduction of cumulative active TB incidence, $C_I\left(t\right)$, and the disease-induced mortality, $C_M\left(t\right)$, is relatively moderate during the short-term perspective (as marked by the vertical line $t=5$ in all three charts), such a reduction becomes more visible in the medium- and long-term perspectives.

Moreover, the cumulative latent TB incidence, $C_L\left(t\right)$, exhibits a noticeable decrease when $q=1$ (see the left-hand chart in Figure 6). Numerical estimations of the key quantities $C_L\left(T\right),C_I\left(T\right),$ and $C_M\left(T\right)$ that can be contrasted by comparing the end-point at $t=T$ of the blue and red curves in the left and right columns of Figure 6, are provided in Table 3 (three upper rows). Furthermore, using a rule similar to (Section 5.1), we can also assess other important quantities characterizing the sole use of either non-sterilizing or sterilizing treatment for active TB infections. Namely, we can estimate:

• The total number of treatments started during the period $[0,T]$ is defined by ${\displaystyle \underset{0}{\overset{T}{\int }}}\gamma I\left(\tau \right)d\tau .$
• The total number of treatment dropouts during the period $[0,T]$ is defined by ${\displaystyle \underset{0}{\overset{T}{\int }}}{\kappa }_i{T}_i\left(\tau \right)d\tau ,i=1,2.$
• The total number of treatment completions during the period $[0,T]$ is defined by ${\displaystyle \underset{0}{\overset{T}{\int }}}{\delta }_i{T}_i\left(\tau \right)d\tau ,i=1,2.$

Numerical estimations of the three additional features defined above are also provided in Table 3 (three lower rows) for both treatment types and in the short, medium, and long-term perspectives. Thus, Figure 5, Figure 6, and Table 3 jointly showcase the benefits of the solely sterilizing treatment over the non-sterilizing one, when no extra cost for sterilizing drugs is involved. In other words, if the budget constraint (22) is ignored, meaning that the functional $J\left(u\right)$ defined by (21) is minimized subject to the ODE system (Section 4), the optimal solution will always be $u^{*}=1,t\in [0,T]$ for any length of the observation period ($T=5,15,$ or 30).

In the following subsection, we provide numerical solutions to the optimal control problem with a budget constraint and discuss how the additional cost of sterilizing drugs may affect the type of treatment assigned to new patients under short-, medium-, and long-term planning.

5.2. Numerical Solutions for Optimal Control

In Section 4, we proposed an optimal control approach for introducing sterilizing TB treatment under the budget constraint (22), which naturally depends on the planning horizon T. Similar to SubSection 5.1, let us focus on short-term ($T=5$ years), medium-term ($T=15$ years), and long-term ($T=30$ years) perspectives and estimate the corresponding values of the budget using the formula (22), that is,

$$B_5:={\displaystyle \underset{0}{\overset{5}{\int }}}\gamma I\left(t\right)dt=52268,B_{15}:={\displaystyle \underset{0}{\overset{15}{\int }}}\gamma I\left(t\right)dt=156758,B_{30}:={\displaystyle \underset{0}{\overset{30}{\int }}}\gamma I\left(t\right)dt=313331,$$

where $I\left(t\right)$ is the solution to the ODE system (2) with $q=0$ (solely non-sterilizing treatment), while assuming that the cost (per unit time) related to the non-sterilizing treatment of one patient entering the class $T_1\left(t\right)$ is normalized to unity. The baseline cumulative costs $B_5,B_{15},$ and $B_{30}$ defined above are obtained using the parameter values from Table 2 for $T=5,T=15,$ and $T=30$, respectively. They are then replaced in the constraint (25) and also in the augmented functional (26) together with the complementary slackness condition (27).

It is worthwhile to point out that we have set $C_1=5$ in the objective functional (26) for all numerical simulations to emphasize that minimizing active TB infections is 5 times more important than reducing latent TB infections because actively infected patients contribute a lot more to the infection spread compared to latently infected ones.

As a reference, we can also assess the cumulative costs corresponding to solely sterilizing treatment using the formula (24) with $u\left(t\right)=1,t\in [0,T]$. In Section 4, we considered two cases (constant and decreasing costs, see p. Section 4) corresponding to the constant and decreasing relative costs of sterilizing treatment and expressed through the values of coefficients $C_2,C_3$ of the cost function $c\left(t\right)$ defined by (23) (see also Figure 2). In the sequel, we determine the cost function $c\left(t\right)$ for these two cases.

Constant cost. 

In this case, we assume $C_3=0$ in (23) and $C_2\in \{1.2,1.5\}$, meaning that the relative (per unit time) cost of sterilizing treatment remains is 20% or 50% more expensive compared to non-sterilizing treatment during the entire observation period $[0,T]$, and for all the three perspectives: short-term ($T=5$ years), medium-term ($T=15$ years), and long-term ($T=30$ years). Notably, such an increment of the relative cost of sterilizing treatment (20% to 50%) was estimated from a recent study by [16] that has taken into account the costs of two sterilizing drugs (Rifampicin and Pyrazinamide). These drugs are manufactured only in India, China, Europe, and the US, and their supply to other countries is somewhat intermittent. The prices of the four TB medications expressed in US dollars (data borrowed from [16]) are summarized in Table 4 below:

Thus, we have two constant cost functions to test,

$$c_1\left(t\right):=1.2\mathrm{and}{c}_2\left(t\right):=1.5,t\in [0,T].$$

Decreasing cost. 

In this case, it is prudent to assume that the longer the observation period, the cheaper the relative cost of sterilizing treatment will be by the end of the observation period. Thus, for both values $C_2\in \{1.2,1.5\}$ we can project that the relative costs of sterilizing and non-sterilizing become equal in $T=30$ years, that is $c\left(30\right)=1$. Under such an assumption, we then define two decreasing cost functions as follows:

$$c_1^d\left(t\right):=1.2-{\displaystyle \frac{t}{150}}\mathrm{and}{c}_2^d\left(t\right):=1.5-{\displaystyle \frac{t}{60}}$$

These functions are illustrated in Figure 7, where the corresponding values of $C_3$ for $T=5,15,30$ are also provided to have consistency with the formulas (23), (25).

Figure 8 displays the cost accumulations under solely non-sterilizing treatment (blue-colored curves) and solely sterilizing treatment with smaller constant and decreasing relative costs on the left ($c_1\left(t\right),c_1^d\left(t\right)$, red and yellow curves, respectively) and the larger ones on the right ($c_2\left(t\right),c_2^d\left(t\right)$, red and yellow curves, respectively). While in the short-term perspective ($T=5$ years) the spending on sterilizing treatment always exceeds the baseline cost (blue curves in both charts), this tendency changes in the medium-term planning ($T=15$) if a smaller relative cost is assumed (meaning that either $c_1\left(t\right)$ or $c_1^d\left(t\right)$ is employed). However, if a larger relative cost is assumed (expressed by the cost functions $c_2\left(t\right)$ or $c_2^d\left(t\right)$), the spending on sterilizing treatment drops below the baseline cost only in the long-term perspective ($T=30$). Moreover, in the long-term perspective ($T=30$ years), the cumulative spending exhibits a notable reduction in all scenarios compared to the habitual practices of non-sterilizing treatment. Table 5 provides the details of cumulative spending given in Figure 8.

From Figure 8 and Table 5 we can easily guess that $u^{*}\left(t\right)=1,t\in [O,T]$ is the optimal solution in the long-term perspective ($T=30$ years) and also in the medium-term perspective ($T=15$ years) if either constant or decreasing relative cost of sterilizing drugs is assumed, that is, if the cost functions $c_1\left(t\right)$ or $c_1^d\left(t\right)$ are employed. In such scenarios, the budget constraint (25) becomes inactive, leading to $\eta =0$ in the complementary slackness condition (27).

Now, we focus on the numerical solving the optimal control problem with budget constraint (25) in the short-term perspective ($T=5$) with the relative costs defined by $c_i\left(t\right)$ and $c_i^d\left(t\right),i=1,2$ and in the medium-term perspective ($T=15$) when the relative costs are higher, meaning that the cost functions $c_2\left(t\right)$ and $c_2^d\left(t\right)$ are employed.

All numerical solutions of the optimal control problem of minimizing the functional (26) subject to the dynamical system (Section 4) with initial conditions (Section 5) were obtained using following scheme, consisting of two steps.

In the first step, we apply direct time discretization to obtain a finite-dimensional problem. The solution to this problem provides an initial guess regarding the structure of the solution. Specifically, for our optimal control problem, we obtain a bang-bang control with two arcs. In the second step, we use the arc-parametrization method, as described in [35], to compute a refined optimal solution. Although sufficient conditions for the optimality of the bang-bang control are not formally proven here, such verification can be performed using the NUDOCCCS software developed in [36] and further employed in[37].

For the numerical implementation of both methods, we use the Applied Modeling Programming Language (AMPL), [38] in combination with the IPOPT solver [39]. The key parameters chosen are $N=2000$ discretization points and an IPOPT tolerance of ${10}^{-8}$.

Table 6 summarizes the types of optimal solutions returned by the AMPL-NUDOCCCS solver. Even though all types of optimal control shown in Table 6 formally comply with the predicted form (31), there are only four scenarios where the optimal control exhibits the verily bang-bang form. Notably, for the two short-term scenarios ($T=5$ years), the optimal bang-bang controls are almost identical, so the difference between them cannot be visually perceived. Therefore, we only show the graph of optimal control for the decreasing cost function $c_1^d\left(t\right)$ in Figure 9 because $u^{*}\left(t\right)$ for the constant cost function $c_1\left(t\right)$ is almost identical. The only difference is the switching point $t_s\in [0,5]$, namely,

$$t_s=0.05\approx 18\mathrm{days}\mathrm{for}{c}_1\left(t\right)=1.2\mathrm{and}{t}_s=0.06\approx 22\mathrm{days}\mathrm{for}{c}_1^d\left(t\right)=1.2-{\displaystyle \frac{t}{150}}.$$

Figure 9 exhibits the bang-bang optimal controls $u^{*}\left(t\right)$ (left column) and the corresponding switching functions $\phi$ (right column), obtained numerically using the formula (32).

Let us now discuss the results exhibited in all charts of Figure 9, starting from the upper row. Because the planning horizon $T=5$ is too short to reflect the long-term advantages of the costly sterilizing TB treatment, and the budget is constrained, this treatment is applied only during the first 18 or 22 days to stay within budget limits, while the shorter period corresponds to the constant lower cost $c_1\left(t\right)$ and a bit longer period applies to the decreasing lower cost $c_1^d\left(t\right)$.

However, there is no way to match the budget constraint in the short-term perspective when the relative cost is assumed higher (that is, cost functions $c_2\left(t\right)$ or $c_2^d\left(t\right)$ are employed). In such situation, only non-sterilizing treatment has to be assigned to all patients with active TB to stay within the limits of the budget.

Nonetheless, from the outcomes of the short-term scenarios one can deduce that the smaller is the relative cost, the longer can be the period for applying the sterilizing treatment while satisfying the budget constraint. This tendency is also illustrated by the medium-term scenarios ($T=15$ years), as shown in the middle and lower rows on Figure 9.

Indeed, when the relative cost of sterilizing treatment is higher and expressed by $c_2\left(t\right)=1.5$ (see the charts (c) and (d) in Figure 9), this better kind of treatment can be assigned to all patients with active TB forms up to the middle of the observation period ($t_s=7.46$ years $\approx 0.5T$), while for a bit lower decreasing cost $c_2^d\left(t\right)$ the application of sterilizing treatment may last far more without violating the budget as shown in the the charts (e) and (f) in Figure 9 where $t_s=12.17\approx 0.8T$. Moreover, if a lower relative cost expressed by the function $c_1\left(t\right)$ or $c_1^d\left(t\right)$ is assumed, the budget constraint (25) becomes inactive leading to $\eta =0$ in the complementary slackness condition (27) and constant optimal control $u^{*}\left(t\right)$.

Notably, in the long-term perspective ($T=30$ years) the optimal control becomes constant, $u^{*}\left(t\right)=1,t\in [0,30]$ regardless of the function expressing the relative cost of sterilizing treatment, meaning that all infectious patients can receive sterilizing treatment without surpassing the assigned budget. These outcomes plainly agree with two charts presented in Figure 8 that exhibit the cumulative relative costs by blue curves (baseline non-sterilizing treatment), yellow curves (lower-cost sterilizing treatment), and red curves (higher-cost sterilizing treatment) when $u\left(t\right)=q=1$ for all $t\in [0,30]$. Since the yellow and red curves are below the blue ones on both charts of Figure 8 at $t=30$ years, the budget constraint is inactive in the case of the long-term planning. However, when $t=15$ years, the yellow curves are already below the blue ones while the red curves are still above them. Therefore, in the case of higher relative cost and medium-term planning (that is, with $c_2\left(t\right)$ or $c_2^d\left(t\right)$ and $T=15$ years), the bang-bang optimal controls help match the budget.

Furthermore, when $t=5$ years, the red curves in Figure 8 are well above the blue ones, while the yellow curves are slightly above the blue ones. Thus, in the case of lower relative cost and short-term planning (that is, with $c_1\left(t\right)$ or $c_1^d\left(t\right)$ and $T=5$ years), the bang-bang optimal controls also help match the budget.

Table 7 summarizes the core indicators of the four optimal controls $u^{*}\left(t\right)$ exhibiting the bang-bang structure, while similar indicators for $u^{*}\left(t\right)=0$ (solely non-sterilizing treatment with $c_2\left(t\right)$ or $c_2^d\left(t\right)$ and $T=5$ years) and for $u^{*}\left(t\right)=1$ (solely non-sterilizing treatment with $c_1\left(t\right)$ or $c_1^d\left(t\right)$ and $T=15$ years) are presented in Table 3. Note that, in the long-term perspective ($T=30$ years), these indicators are the same as presented in the fourth column of Table 5 as they do not depend on the cost function $c\left(t\right)$.

Keeping in mind a minimal difference between the optimal bang-bang controls corresponding to the lower relative cost functions $c_1\left(t\right)$ and $c_1^d\left(t\right)$, we now compare the values of indicators presented in the second and third columns of Table 7 and observe that the decreasing nature of the cost function $c_1^d\left(t\right)$ has a positive effect on these indicators even in a short-term perspective ($T=5$ years). By contrasting the values of indicators in the second and third columns of Table 7 with their baseline quantities given in the fifth column of Table 3 corresponding to the solely non-sterilizing treatment, we can observe that the application of sterilizing treatment during only the first 18-22 day (in accordance with the optimal bang-bang controls presented in Figure 9(a)) may improve the epidemiological situation without incurring into extra-costs.

Under the medium-term planning ($T=15$ years), the benefits of optimal bang-bang controls shown in the charts (c) and (e) of Figure 9 become more evident compared to the baseline case of solely non-sterilizing treatment, while total costs remain unchanged. Such benefits can be perceived by comparing the indicator values in the last two columns of Table 7 with their baseline quantities in the sixth column of Table 3. In other words, bang-bang optimal controls enhance the cost-effectiveness of treatment as they improve the indicators predicted in the baseline case at no extra cost.

Notably, under short-term planning ($T=5$ years) and under medium-term perspective ($T=15$ years) with a higher relative cost (bang-bang optimal control), the budget constraint (25) is active, meaning that the total budget is spent. However, if the relative cost of sterilizing treatment is moderate (that is, cost functions $c_1\left(t\right)$ or $c_1^d\left(t\right)$ are employed) or if the planning is done for the long-term perspective ($T=30$ years), the optimal control is $u^{*}\left(t\right)=1$ and with this control the budget constraint (25) becomes inactive, meaning that there will be surplus compared to the baseline budget. The latter emphasizes even more the cost-effectiveness of sterilizing treatment in the medium- and long-term perspective.

5.3. Practical Outlook

To conclude this section, we summarize the outcomes obtained in the previous subsections and outline general recommendations stemming from our computational experiments.

Using the estimated quantities displayed in Table 3, Table 5, and Table 7, we can assess the potential benefits and extra costs of switching entirely to the sterilizing treatment ($q=1$ in system (2) or $u\left(t\right)=1,t\in [0,T]$ in system (Section 4)) and employing the optimal control strategy $u^{*}\left(t\right),t\in [0,T]$ while assuming either constant or decreasing relative costs of the sterilizing drugs. Such assessments will be made in comparison to the habitual or “baseline” scenario, which encompasses assigning the non-sterilizing treatment to all new patients with active TB infections. In mathematical terms, this baseline scenario is expressed by setting $q=0$ in the ODE system (2) or taking $u\left(t\right)=0$ for all $t\in [0,T]$ in the system (Section 4).

In the sequel, we consider three key indicators characterizing the benefits of either solely sterilizing treatment ($u\left(t\right)=1$ for all $t\in [0,T]$) or the bang-bang optimal controls $u^{*}\left(t\right),t\in [0,T]$ displayed in the charts (a), (c), and (e) in Figure 9, compared to the baseline scenario ($u\left(t\right)=0$ for all $t\in [0,T]$), namely:

• The total number of avoided active TB infections is assessed by comparing the cumulative number of new active TB infections under the baseline and other scenarios during the observation period $[0,T]$.
• The total number of avoided latent TB infections is assessed by comparing the cumulative number of new latent TB infections under the baseline and other scenarios during the observation period $[0,T]$.
• The total number of avoided TB-induced deaths is assessed by comparing the cumulative mortality under the baseline and other scenarios during the observation period $[0,T]$.

To assess the extra costs or surpluses, we contrast the data from Table 5 for solely non-sterilizing and solely sterilizing treatments, and recall that the optimal bang-bang control strategies $u^{*}\left(t\right)$ do not incur any extra costs.

First, we address the short-term perspective where the observation period equals 5 years ($T=5$). The bar diagram in Figure 10 displays the three key indicators defined above and the extra costs, all in the form of percentages relative to the baseline scenario (non-sterilizing treatment only), where the latter corresponds to the horizontal axis ($0\%$). This diagram clearly shows how solely sterilizing treatment (to a much greater extent) and optimal bang-bang control strategies (just moderately) contribute to avoiding future latent TB infections, even in the short term. It is also visible that solely sterilizing treatment guarantees higher benefits compared to the optimal control strategies with either constant or decreasing relative cost of sterilizing drugs.

Nonetheless, such higher benefits come with a considerable extra cost (see the last column of two bars expressing over budget in Figure 10). Indeed, we can see that an over-budget between $33.61\%$ and $37.28\%$ corresponding to the sole use of sterilizing treatment (compared to the optimal control strategies) actually can cover about $80\%$ of extra-avoided active latent TB infections together with about $9\%$ of additionally averted active TB infections and $8\%$ of the TB-induced deaths. Alternatively, bang-bang optimal control strategies do not seem very promising under short-term planning, as their benefits relative to the baseline non-sterilizing treatment are rather moderate.

Let us now examine the same indicators in the medium-term perspective with the observation period of 15 years ($T=15$). The bar diagram in Figure 11 clearly shows that the outcomes of solely sterilizing treatment and bang-bang optimal control strategies become more aligned. Moreover, when a lower relative cost of sterilizing drugs is assumed (expressed by the cost functions $c_1\left(t\right)$ and $c_1^d\left(t\right)$), the healthcare system may even save between $11.63\%$ and $14.87\%$ of an initially projected budget by assigning only sterilizing treatment to all new patients with active TB infections. Therefore, from a medium-term perspective, the use of more expensive sterilizing drugs at least during the first half of the observation period will eventually pay off by drastically reducing the pool of latently TB-infected individuals and averting more than $20\%$ of active TB infections and TB-induced deaths.

Finally, Figure 12 shows the expected benefits and potential savings for the healthcare system when planning is done at a long-term perspective. Here, we highlight that, through long-term planning and providing solely more expensive sterilizing treatment to all patients, over $45\%$ of active TB infections and TB-related deaths, together with more than $90\%$ of latent TB infections, can be averted at a lower total cost than offering a cheaper non-sterilizing treatment. This outcome has a logical explanation: sterilizing treatment drastically reduces the number of latent TB infections and slows down the disease transition even if the disease transmission and the treatment abandonment rates remain unaltered.

6. Conclusions

In this work, we developed and analyzed a mathematical model for tuberculosis (TB) transmission that incorporates two types of treatments, sterilizing and non-sterilizing, where a fraction of actively TB-infected patients may receive a sterilizing treatment while the rest will undergo a non-sterilizing treatment due to resource-limited settings. Our model extends the prior TB-transmission model presented in [7] by explicitly considering the effects of non-sterilizing therapy, which reduces bacterial load but allows latent infections to persist, in contrast to sterilizing treatment, which eliminates both active and dormant bacilli of Mycobacterium tuberculosis.

We first performed a theoretical analysis of the model, identifying the disease-free and endemic equilibria and deriving the basic reproductive number, ${\mathcal{R}}_0$, as a function of the proportion of infectious individuals assigned for sterilizing treatment. Our results confirm that TB eradication is only possible when ${\mathcal{R}}_0<1$, and the latter is achievable only if a sufficiently high fraction of patients undergo sterilizing treatment. However, in settings where budget constraints limit access to sterilizing drugs, TB may still persist at endemic levels due to the increase of latently TB-infected individuals.

The model’s long-term behavior was further explored through stability analysis, demonstrating that if ${\mathcal{R}}_0<1$, the disease-free equilibrium is globally asymptotically stable (GAS), meaning TB will eventually be eradicated. Conversely, if ${\mathcal{R}}_0>1$, the system stabilizes at an endemic equilibrium, whose coordinates are influenced by the proportion of patients receiving sterilizing treatment. In other words, by increasing the share of patients who receive sterilizing treatment, the endemic equilibrium may vanish and the epidemiological system may then evolve toward the disease-free equilibrium. This underscores the importance of ensuring adequate access to sterilizing drugs to maximize their epidemiological impact.

To determine an optimal distribution between the two TB treatment types, we employed an optimal control approach where the proportion of infectious individuals receiving sterilizing treatment may vary in time. The results reveal that, in the absence of budget constraints, the immediate and total change to sterilizing treatment for all new patients leads to the fastest and most sustainable decline in TB incidence. However, when a budget constraint is introduced, it is optimal to start immediately with a more expensive sterilizing treatment, at least at the beginning of the observation period, and then continue with a more affordable non-sterilizing treatment to match the budget, not the other way around. The latter was mathematically expressed by the bang-bang optimal controls with only one switch between the two treatment types.

Also, it has been shown that switching time is linked to the planning horizon, which plays an essential role in achieving significant reductions in both active and latent TB cases and in disease-induced mortality. Namely, as the planning horizon becomes longer, more active and latent TB infections and disease-related deaths can be averted under the bang-bang optimal control strategies.

Numerical simulations further demonstrated that bang-bang optimal controls are also sensitive to the relative cost of the sterilizing treatment. In scenarios where the cost decreases over time, an exclusively sterilizing treatment becomes increasingly more affordable in the medium- and long-term perspectives and leads to more visible epidemiological outcomes. If the planning horizon is long enough and even if the cost of sterilizing drugs remains relatively high and non-decreasing over time, the total cost of administering sterilizing drugs to all patients over the entire observation period will stay below the projected baseline budget.

Moreover, in quantitative terms, with a planning horizon of 30 years, the healthcare system may save up to $38\%$ of the baseline budget if only a more expensive sterilizing treatment is provided to all patients with active TB infections. Such an outcome can be intuitively explained by the fact that non-sterilizing treatment increases the pool of latently TB-infected people, while sterilizing treatment reduces it.

From a healthcare policy perspective, our findings highlight the need for long-term planning of the TB treatment programs in low- and middle-income countries. While complete and immediate replacement of non-sterilizing treatment with sterilizing drugs remains the best scenario, a well-planned, a budget-sensitive longer-term planning can still provide substantial epidemiological benefits. Future work may focus on incorporating additional real-world constraints, such as patient motivation for treatment adherence and completion, to enhance decision-making in TB treatment policies.

In conclusion, this study underscores the critical role of longer-term planning and optimizing the distribution between two types of TB treatments, particularly in resource-limited settings. By balancing epidemiological effectiveness with financial constraints, policymakers can design TB control programs that maximize health benefits while remaining economically sustainable.