Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control

Erick Mutwiri Kirimi; Jeconiah Okelo; Mark Kimathi; Kenneth Ngure

doi:doi:10.11648/j.ajmcm.20251004.12

Research Article |

| Peer-Reviewed

Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control

Erick Mutwiri Kirimi^*

, Jeconiah Okelo

, Mark Kimathi

, Kenneth Ngure

Published in American Journal of Mathematical and Computer Modelling (Volume 10, Issue 4)

Received: 9 August 2025 Accepted: 21 August 2025 Published: 14 October 2025

Views: Downloads:

Download PDF

Share This Article

Twitter
Linked In
Facebook

Abstract

In this paper, a deterministic mathematical model illustrating the transmission dynamics of pulmonary tuberculosis and pneumonia co-infection is formulated, incorporating a drug-resistant strain. The model employs a Holling-type functional response to capture the impact of natural immunity on the progression from latent tuberculosis infection to active disease, as well as its role in controlling drug-resistant pulmonary tuberculosis-pneumonia co-infections. The model is extended to include optimal control theory, aimed at identifying strategies to minimize co-infections using prevention, screening of latently infected individuals, and treatment as control variables. Pontryagin’s Maximum Principle is applied to characterize the optimal control system. The resulting optimality system is then solved numerically using the Runge-Kutta-based forward-backward sweep method. Numerical simulations demonstrate that enhancing natural immunity among latently infected individuals significantly reduces the number of co-infected cases. The optimal control analysis indicates that the most effective strategy for controlling or reducing co-infections of drug-resistant tuberculosis and pneumonia is the combined optimization of infection prevention and screening of latently infected individuals. These findings underscore the importance of scaling up preventive measures against pulmonary tuberculosis and opportunistic pneumonia, alongside screening efforts, to effectively control co-infections. Additionally, the study recommends strengthening immunity among latently infected populations to further reduce the prevalence of co-infections.

Published in	American Journal of Mathematical and Computer Modelling (Volume 10, Issue 4)
DOI	10.11648/j.ajmcm.20251004.12
Page(s)	121-144
Creative Commons	This is an Open Access article, distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution and reproduction in any medium or format, provided the original work is properly cited.
Copyright	Copyright © The Author(s), 2025. Published by Science Publishing Group

Keywords

Pulmonary Tuberculosis, Co-infection, Latently Infected, Natural Immunity, Drug-resistant Strain, Optimal Control

1. Introduction

Pulmonary tuberculosis (TB), which primarily affects the lungs, remains one of the most persistent and deadly contagious diseases globally

[1]

. For instance, in 2023, 10.8 million people worldwide fell ill with pulmonary TB, and approximately 1.2 million succumbed to the disease

[2]

. The impact of pulmonary tuberculosis is exacerbated by its frequent co-occurrence with other infections, notably pneumonia

[3]

. The coexistence of pulmonary tuberculosis and pneumonia limits treatment options and hinders effective case management, thereby contributing to the emergence of drug-resistant strains. This presents a significant challenge in controlling disease transmission and leads to increased mortality and adverse outcomes associated with pulmonary tuberculosis

[4]

. Preventing the occurrence of pneumonia in individuals with drug-sensitive pulmonary TB is therefore crucial. This underscores the urgent need for integrated approaches, including mathematical modeling, to address both the interaction between pulmonary tuberculosis and opportunistic pneumonia, and the challenges posed by resistant strains.

Mathematical modeling serves as a valuable tool for analyzing the transmission dynamics of co-infections and for designing and optimizing intervention strategies

[5]

. When integrated with optimal control theory, such models can evaluate the timing, intensity, and cost-effectiveness of strategies aimed at minimizing disease burden

[6]

. In this study, a mathematical model is formulated to capture the co-dynamics of pulmonary tuberculosis and opportunistic pneumonia in the presence of a drug-resistant TB strain. The model incorporates time-dependent control strategies designed to reduce the burden of co-infections and curb the spread of the resistant strain. By applying Pontryagin’s Maximum Principle, the necessary conditions for optimal control are derived, and numerical simulations are conducted to assess the impact of combined interventions.

Numerous mathematical models have been formulated to study pulmonary tuberculosis as an independent disease, with a focus on strategies to control the transmission of infection

[7-11]

. For example, Ochieng

[12]

conducted a study on the mathematical modeling of pulmonary tuberculosis transmission, incorporating environmental factors and optimal control. The study compared the use of vaccination and public awareness campaigns as control strategies through optimal control analysis. The results indicated that vaccination is the more cost-effective strategy for combating pulmonary tuberculosis transmission in Kenya. Kang et al.

[13]

, in their study, considered the effects of vaccination, awareness of tuberculosis prevention, screening of exposed individuals, and clearance of pathogens in a contaminated environment on the transmission of pulmonary tuberculosis. The optimal control analysis revealed that combining all four strategies is the most cost-effective approach to reducing the prevalence of pulmonary tuberculosis in the community. A deterministic mathematical model of pulmonary tuberculosis, taking into consideration vaccination and a drug-resistant strain, was formulated and analyzed by

[14]

. The study considered prevention, treatment, and vaccination as control measures to reduce the transmission of pulmonary tuberculosis in the population. The optimal control analysis revealed that the combination of prevention and vaccination is the most cost-effective strategy for combating disease transmission and reducing drug-resistant cases. Mengistu et al.

[15]

conducted a study on the mathematical modeling of pulmonary tuberculosis incorporating drug-resistant strains. Their study analyzed the effects of vaccination, distancing, and treatment for both drug-susceptible and drug-resistant strains to eliminate drug-resistant pulmonary tuberculosis in the population using an optimal control approach. The findings revealed that successful treatment of the drug-susceptible strain, as a single strategy, is the most effective approach for eliminating drug-resistant pulmonary tuberculosis. The most cost-effective combination for combating drug-resistant pulmonary tuberculosis was found to be distancing and the successful treatment of the drug-susceptible strain.

Recent mathematical models have been developed to explore various aspects of pneumonia as an independent disease, contributing to efforts aimed at its containment. For example, Aldila et al.

[16]

formulated a mathematical model to investigate the potential impact of pneumonia treatment on disease transmission dynamics. The study included an optimal control analysis of the model, and the results showed that implementing treatment during the early stages of infection is critical to averting outbreaks in the community. A mathematical model of pneumonia in the presence of drug resistance threats was formulated by

[17]

. An optimal control analysis was performed to determine the effects of various control strategies against antibiotic resistance. The results indicated that a combination of prevention, treatment, and immunity-enhancing efforts is the most effective strategy for minimizing the prevalence of resistant strains in the population. Swai et al.

[18]

formulated a mathematical model for the transmission of pneumonia in the presence of drug resistance. The model considered control interventions, including vaccination, public health education, and treatment. The study applied optimal control analysis to determine the interventions that minimize infections from both drug-sensitive and drug-resistant strains. The results indicated that a vaccination program is the most cost-effective strategy for controlling infections in resource-limited settings. A mathematical model incorporating optimal control analysis was formulated by

[19]

to describe bimodal pneumonia transmission behavior in a vulnerable compartment. The model included prevention and treatment as control measures. Simulation results showed that a combination of prevention and treatment is an effective strategy for controlling the spread of pneumonia during epidemics.

Mathematical modeling of pulmonary tuberculosis and pneumonia co-infection has received relatively little attention, despite clinical evidence indicating a lethal synergism between the two diseases

[20]

. The only recent study on the mathematical modeling of pulmonary tuberculosis and pneumonia is presented in

[21]

. That study employed an SEIR (Susceptible, Exposed, Infectious, Recovered) transmission framework to formulate the model, and its findings indicated that treating tuberculosis is the most effective strategy for containing the co-infection. The present study builds on the work of Gweryina et al.

[21]

by formulating a novel mathematical model for pulmonary tuberculosis in the presence of opportunistic pneumonia, with the aim of accurately representing the natural progression of the co-infection and proposing strategies to prevent or mitigate their coexistence.

The novelty of the formulated model is presented in this section. Previous mathematical models describing the transmission dynamics of pulmonary tuberculosis and pneumonia co-infection have not incorporated several key aspects addressed in the current study. First, the model includes a drug-resistant strain of tuberculosis, acknowledging that the coexistence of pulmonary tuberculosis and pneumonia often results in limited treatment options and suboptimal case management, both of which contribute to the emergence of drug-resistant strains

[22]

. Accordingly, this study identifies strategies to reduce the emergence and spread of drug-resistant tuberculosis within the population. Second, many existing models assume that the progression rate from latent infection to active pulmonary tuberculosis is directly proportional to the number of latently infected individuals. However, in reality, this progression is significantly influenced by an individual’s natural immunity

[23]

. Through mathematical modeling, this study rigorously explores the role of natural immunity in the reactivation of latent infections and its impact on controlling or reducing the co-existence of drug-resistant pulmonary tuberculosis and opportunistic pneumonia. Third, the model is extended to incorporate optimal control theory, enabling the identification of strategies that should be optimized to minimize the co-existence of both drug-sensitive and drug-resistant tuberculosis with opportunistic pneumonia in the community. The primary objective of this research is to determine the optimal strategies for reducing the co-infection burden of pulmonary tuberculosis and opportunistic pneumonia in the presence of a drug-resistant strain. Based on an extensive review of the literature, no previous study has addressed the specific research gaps identified in the present work.

The paper is organized as follows: Section 2 presents the model formulation, while Section 3 provides the mathematical analysis. Section 4 extends the model to incorporate optimal control theory. Section 5 presents numerical simulations to determine optimal strategies for reducing the prevalence of pulmonary tuberculosis and opportunistic pneumonia co-infections in the presence of a drug-resistant strain. Finally, Section 6 concludes the study.

2. Model Formulation

The model partitions the total human population at time

t

, denoted by

N (t)

, into thirteen compartments based on their epidemiological status. The model variables used in the formulation are

S (t)

V (t)

E_{w} (t)

E_{r} (t)

I_{w} (t)

I_{r} (t)

I_{wP} (t)

I_{rP} (t)

T_{w} (t)

T_{r} (t)

T_{wP} (t)

T_{rP} (t)

, and

R (t)

, as defined in Table 1.

Table 1. Definition of the model variables.

Variable	Definition
$S (t)$	Individuals susceptible to pulmonary tuberculosis
$V (t)$	Individuals vaccinated against pulmonary tuberculosis
$E_{w} (t)$	Individuals latently infected with the drug-sensitive strain of pulmonary TB
$E_{r} (t)$	Individuals latently infected with the drug-resistant strain of pulmonary TB
$I_{w} (t)$	Individuals infected with the drug-sensitive strain of pulmonary tuberculosis
$I_{r} (t)$	Individuals infected with the drug-resistant strain of pulmonary tuberculosis
$I_{wP} (t)$	Individuals co-infected with the drug-sensitive strain of TB and pneumonia
$I_{rP} (t)$	Individuals co-infected with the drug-resistant strain of TB and pneumonia
$T_{w} (t)$	Individuals undergoing treatment for the drug-sensitive strain of pulmonary TB
$T_{r} (t)$	Individuals undergoing treatment for the drug-resistant strain of pulmonary TB
$T_{wP} (t)$	Individuals undergoing treatment for co-infection with the drug-sensitive strain of pulmonary TB and pneumonia
$T_{rP} (t)$	Individuals undergoing treatment for co-infection with the drug-resistant strain of pulmonary TB and pneumonia
$R (t)$	Individuals who have recovered

Thus,

N (t) = S (t) + V (t) + E_{w} (t) + E_{r} (t) + I_{w} (t) + I_{r} (t)

+ I_{wP} (t) + I_{rP} (t) + T_{w} (t) + T_{r} (t) + T_{wP} (t) + T_{rP} (t) + R (t)

The model introduced in this study incorporates a saturated progression rate from latent infection to active disease, accounting for natural immunity through Holling-type saturation functions:

\frac{ψ_{1} E_{w} (t)}{1 + m E_{w} (t)}

and

\frac{ψ_{2} E_{r} (t)}{1 + m E_{r} (t)}

for drug-sensitive tuberculosis and drug-resistant tuberculosis, respectively, as described in

[24]

. The number of latently infected individuals effectively saturates when immunity levels are high. In these functions,

ψ_{1}

and

ψ_{2}

represent the reactivation rates of latent infections, while

m

quantifies the delay in reactivation due to enhanced immunity. When

m

approaches zero, corresponding to cases of extremely low immunity, the functions approximate

ψ_{1} E_{w} (t)

and

ψ_{2} E_{w} (t)

, respectively.

The formulation of the model involves transitions between compartments, representing changes from one disease state to another, and utilizes the parameters described in Table 2.

Table 2. Description of the model parameters.

Parameter	Description
$Γ$	Recruitment rate of individuals into the population.
$Q$	Proportion of the population that is vaccinated.
$ρ$	Efficacy of the pulmonary tuberculosis vaccine.
$λ_{w}$ , $λ_{r}$ , $λ_{P}$ .	Forces of infection for drug-sensitive tuberculosis, drug-resistant tuberculosis, and pneumonia, respectively.
$μ$	Natural death rate.
$φ_{1}$ , $φ_{2} .$	Rates at which individuals latently infected with the drug-sensitive and drug-resistant strains of pulmonary TB are screened.
$θ_{1}$ , $θ_{2} .$	Treatment rates of individuals with drug-sensitive and drug-resistant tuberculosis.
$θ_{3}$ , $θ_{4}$ .	Rates at which individuals with drug-sensitive and drug-resistant TB-pneumonia co-infection are screened.
$ψ_{1}$ , $ψ_{2}$ .	Rates of progression from latent drug-sensitive and drug-resistant strain infections to TB disease.
$b_{1}$ , $b_{2}$ .	Rates at which individuals undergoing treatment for the drug-sensitive strain of TB, and for co-infection with the drug-sensitive strain of TB and pneumonia, respectively, develop a resistant strain.
$ξ_{1}$ , $ξ_{2}$ .	Rates at which individuals co-infected with the drug-sensitive and drug-resistant strains of TB and pneumonia recover from acute pneumonia.
$υ_{1}$ , $υ_{2}$ .	Rates at which individuals co-infected with the drug-sensitive and drug-resistant strains of TB and pneumonia recover from acute pneumonia while undergoing treatment.
$δ_{1}$ , $δ_{2}$ , $δ_{3}$ , $δ_{4}$ , $δ_{5}$ , $δ_{6}$ , $δ_{7}$ , $δ_{8}$ .	Mortality rates due to the drug-sensitive strain of TB; drug-resistant strain of TB; co-infection with the drug-sensitive strain of TB and pneumonia; co-infection with the drug-resistant strain of TB and pneumonia; drug-sensitive TB in individuals undergoing treatment; drug-resistant TB in individuals undergoing treatment; co-infection with the drug-sensitive strain of TB and pneumonia in individuals receiving treatment; and co-infection with the drug-resistant strain of TB and pneumonia in individuals undergoing treatment, respectively.
$φ_{3}, φ_{3} .$	Rates at which individuals with the drug-sensitive and drug-resistant strains of TB recover while undergoing treatment, respectively.
$m$	An immunity parameter that delays the progression to severe tuberculosis.
$κ_{1}$	Rate at which individuals who have recovered from TB become susceptible again
$β_{w}$ , $β_{r}$ , $β_{P}$ .	Transmission rates of drug-sensitive TB, drug-resistant TB, and pneumonia, respectively.
$ϵ_{1}$ , $ϵ_{2}$ , $ϵ_{3}$ , $ϵ_{4}$ , $ϵ_{5}$ , $ϵ_{6}$ .	Transmission coefficients of tuberculosis from: the drug-sensitive TB strain only; co-infection with the drug-sensitive TB strain and pneumonia during treatment; the drug-sensitive TB strain during treatment; the drug-resistant TB strain; co-infection with the drug-resistant TB strain and pneumonia during treatment; and the drug-resistant TB strain only during treatment, respectively.
$ϵ_{7}$ , $ϵ_{8}$ , $ϵ_{9}$ .	Transmission coefficients of pneumonia from: individuals co-infected with the drug-resistant strain of TB and pneumonia; individuals co-infected with the drug-sensitive strain of TB and pneumonia who are undergoing treatment; and individuals co-infected with the drug-resistant strain of TB and pneumonia during treatment, respectively.

The formulation of the new model is guided by the following assumptions:

1) Susceptible individuals and those who have been vaccinated may become infected with either the drug-sensitive or drug-resistant strain of pulmonary TB upon exposure.

2) Defaulting on treatment for drug-sensitive pulmonary TB can lead to the development of a drug-resistant strain.

3) Successful treatment is available for both drug-sensitive and drug-resistant TB strains.

4) Individuals undergoing treatment for co-infection with pulmonary TB and pneumonia recover from pneumonia first, as it is an acute disease.

5) Progression from the latently infected state (both drug-sensitive and drug-resistant) to active pulmonary TB disease depends on natural immunity.

Based on the variables, parameters, and assumptions, the schematic diagram illustrating the transmission dynamics of pulmonary tuberculosis in the presence of opportunistic pneumonia is presented in Figure 1.

Download: Download full-size image

Figure 1. Schematic diagram illustrating the transmission dynamics of pulmonary tuberculosis and pneumonia co-infection incorporating a drug-resistant strain.

The transmission model is formulated as a system of ordinary differential equations, as follows:

\frac{dS}{dt} = (1 - Q) Γ + κ_{1} R - (μ + λ_{r} + λ_{w}) S

\frac{dV}{dt} = Q Γ - [μ + (1 - ρ) λ_{r} + (1 - ρ) λ_{w}] V

\frac{d E_{w}}{dt} = λ_{w} S + (1 - ρ) λ_{w} V - (μ + \frac{ψ_{1}}{1 + m E_{w}} + φ_{1}) E_{w}

\frac{d E_{r}}{dt} = λ_{r} S + (1 - ρ) λ_{r} V - (μ + \frac{ψ_{2}}{1 + m E_{r}} + φ_{2}) E_{r}

\frac{d I_{w}}{dt} = \frac{ψ_{1} E_{w}}{1 + m E_{w}} + ξ_{1} I_{wP} - (μ + θ_{1} + δ_{1} + λ_{P}) I_{w}

\frac{d I_{r}}{dt} = \frac{ψ_{2} E_{r}}{1 + m E_{r}} + ξ_{2} I_{rP} - (μ + δ_{2} + θ_{2} + λ_{P}) I_{r}

,(1)

\frac{d I_{wP}}{dt} = λ_{P} I_{w} - (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) I_{wP},

\frac{d I_{rP}}{dt} = λ_{P} I_{r} + b_{3} I_{wP} - (θ_{4} + μ + δ_{4} + ξ_{2}) I_{rP}

\frac{{dT}_{w}}{dt} = θ_{1} I_{w} + υ_{1} T_{wP} - (φ_{3} + b_{1} + δ_{5} + μ) T_{w}

\frac{{dT}_{r}}{dt} = θ_{2} I_{r} + υ_{2} T_{rP} + b_{1} T_{w} - (φ_{4} + δ_{6} + μ) T_{r}

\frac{{dT}_{wP}}{dt} = θ_{3} I_{wP} - (υ_{1} + b_{2} + δ_{7} + μ) T_{wP}

\frac{{dT}_{rP}}{dt} = θ_{4} I_{rP} + b_{2} T_{wP} - (υ_{2} + δ_{8} + μ) T_{rP}

\frac{dR}{dt} = φ_{1} E_{w} + φ_{2} E_{r} + φ_{3} T_{w} + φ_{4} T_{r} - (μ + κ_{1}) R

where:

λ_{w}

λ_{r}

, and

λ_{P}

are defined as follows:

λ_{w} = β_{w} (I_{wP} + ϵ_{1} I_{w} + ϵ_{2} T_{wP} + ϵ_{3} T_{w})

λ_{r} = β_{r} (I_{rP} + ϵ_{4} I_{r} + ϵ_{5} T_{rP} + ϵ_{6} T_{r})

λ_{P} = β_{P} (I_{wP} + ϵ_{7} I_{rP} + ϵ_{8} T_{wP} + ϵ_{9} T_{rP})

It is noted that

ϵ_{3} < ϵ_{2} < ϵ_{1}

, and

0 < ϵ_{1}

ϵ_{2}

ϵ_{3} < 1

Similarly,

ϵ_{6} < ϵ_{5} < ϵ_{4}

and

0 < ϵ_{4}

ϵ_{5}

ϵ_{6} < 1

Moreover,

ϵ_{9} < ϵ_{8} < ϵ_{7}

and

0 < ϵ_{7}

ϵ_{8}

ϵ_{9} < 1

All parameters are assumed to be positive constants, and the initial conditions of system (1) are given as follows:

S (0) > 0

V (0) \geq 0

E_{w} (0) \geq 0

E_{r} (0)

I_{w} (0) \geq 0

I_{r} (0) \geq 0

I_{wP} (0) \geq 0

I_{rP} (0) \geq 0

T_{w} (0) \geq 0

T_{r} (0) \geq 0

T_{wP} (0) \geq 0

T_{rP} (0) \geq 0

,and

R (0) \geq 0

3. Model Analysis

In this section, the key properties of the model are analyzed to identify the conditions that influence the transmission of infections in the population. The analysis addresses the positivity and boundedness of solutions, equilibrium points, reproduction numbers, and stability analysis.

3.1. Positivity of Solutions

The model system (1) describes the dynamics of the human population, and therefore its solutions must remain positive for all

t > 0

, as demonstrated in Theorem 1.

Theorem 1: Given that the initial conditions

S (0) > 0

V (0) \geq 0

E_{w} (0) \geq 0

E_{r} (0)

I_{w} (0) \geq 0

I_{r} (0) \geq 0

I_{wP} (0) \geq 0

I_{rP} (0) \geq 0

T_{w} (0) \geq 0

T_{r} (0) \geq 0

T_{wP} (0) \geq 0

T_{rP} (0) \geq 0

R (0) \geq 0

, then the solution set

\{S (t), V (t), E_{w} (t), E_{r} (0), I_{w} (t), I_{r} (t), I_{wP} (t), I_{rP} (t), T_{w} (t), T_{r} (t), T_{wP} (t), T_{rP} (t), R (t)\}

remains non-negative for all

t > 0 .

Proof: The first equation of system (1) is considered to demonstrate the positivity of the model solutions. The equation is given by:

\frac{dS}{dt} = (1 - Q) Γ + κ_{1} R - (μ + λ_{r} + λ_{w}) S

(2)

By the comparison theorem, equation (2) can be rewritten as:

\frac{dS}{dt} \geq - (μ + λ_{r} + λ_{w}) S

(3)

Separating variables and integrating from

0

t

yields:

S (t) \geq S (0) e^{- \int_{0}^{t} (μ + λ_{r} + λ_{w}) dt} \geq 0

(4)

Equation (4) proves that the solution of equation (1) remains positive, provided that

S (0) > 0

Similarly, it can be shown that the solutions for the other model equations in system (1) are also non-negative.

Thus, the solution set

\{S (t), V (t), E_{w} (t), E_{r} (t), I_{w} (t), I_{r} (t), I_{wP} (t), I_{rP} (t), T_{w} (t), T_{r} (t), T_{wP} (t), T_{rP} (t), R (t)\} \geq 0 \forall t > 0

(5)

This proves that the solutions of model system (1) remain positive for all

t > 0 .

3.2. Boundedness of Solutions

The model system (1) pertains to the human population; therefore, it is necessary to demonstrate that its solutions remain bounded for all

t > 0

. The following theorem is used to prove the boundedness of the solutions:

Theorem 2: Given positive initial conditions, the feasible region is defined as:

Ω = \{[S (t), V (t), E_{w} (t), E_{r} (t), I_{w} (t), I_{r} (t), I_{wP} (t), I_{rP} (t), T_{w} (t), T_{r} (t), T_{wP} (t), T_{rP} (t), R (t)] \in R_{+}^{13} : N \leq \frac{Γ}{μ}\}

(6)

Proof: The sum of all the equations in system (1) represents the total human population in the model and satisfies the following differential equation:

\frac{dN}{dt} = Γ - μN - δ_{1} I_{w} - δ_{2} I_{r} - δ_{3} I_{wP} - δ_{4} I_{rP} - δ_{5} T_{w} - δ_{6} T_{r} - δ_{7} T_{wP} - δ_{8} T_{rP}

(7)

In the absence of disease-induced mortality (i.e.,

δ_{1} = δ_{2} = δ_{3} = δ_{4} = δ_{5} = δ_{6} = δ_{7} = δ_{8} = 0

), equation (7) reduces to:

\frac{dN}{dt} = Γ - μN

(8)

Integrating equation (8) and applying the initial condition yields:

N (t) \leq \frac{Γ}{μ} + (N (0) - \frac{Γ}{μ}) e^{- μt}

(9)

t \to \infty

N (t) \to \frac{Γ}{μ}

, implying that

0 \leq N (t) \leq \frac{Γ}{μ}

Thus, the feasible solution of system (1) enters and remains in the region:

Ω = \{[S (t), V (t), E_{w} (t), E_{r} (t), I_{w} (t), I_{r} (t), I_{wP} (t), I_{rP} (t), T_{w} (t), T_{r} (t), T_{wP} (t), T_{rP} (t), R (t)] \in R_{+}^{13} : N \leq \frac{Γ}{μ}\}

(10)

Therefore, the basic model is well-posed both epidemiologically and mathematically, and it is sufficient to study its dynamics within the region

Ω

3.3. Disease Free Equilibrium Point and the Control Reproduction Number

The disease-free equilibrium point of system (1) is obtained by setting all latently infected classes, infectious classes, treatment classes, and the recovered class to zero, i.e.

E_{w} = E_{r} = I_{w} = I_{r} = I_{wP} = I_{rP} = T_{w} = T_{r} = T_{wP} = I_{rP} = R = 0

Therefore, system (1) has a disease-free equilibrium given by:

B^{0} (S^{0}, V^{0}, E_{w}^{0}, E_{r}^{0}, I_{w}^{0}, I_{r}^{0}, I_{wP}^{0}, I_{rP}^{0}, T_{w}^{0}, T_{r}^{0}, T_{wP}^{0}, T_{rP}^{0}, R^{0}) = (\frac{(1 - Q) Γ}{μ}, \frac{Q Γ}{μ}, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

(11)

3.4. Control Reproduction Number

The control reproduction number is defined as the expected number of secondary infections produced by a single infected individual during their entire infectious period in a population that is not entirely susceptible due to the presence of control measures

[25]

. The Next-Generation Matrix method, as described by

[26]

, is used to compute the control reproduction number, which is determined as a spectral radius of the matrix

F V^{- 1}

Let

X = {(E_{w}, E_{r}, I_{w}, I_{r}, I_{wP}, I_{rP}, T_{w}, T_{r}, T_{wP}, T_{rP})}^{T}

Then, from system (1), it follows that:

\frac{dX}{dt} = f - v

,(12)

where

f

and

v

are the matrices representing the new infection and transition terms, respectively, defined as:

f = [\begin{matrix} λ_{w} S + (1 - ρ) λ_{w} V \\ λ_{r} S + (1 - ρ) λ_{r} V \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}]

v = [\begin{matrix} (\frac{ψ_{1}}{1 + m E_{w}} + φ_{1} + μ) E_{w} \\ (\frac{ψ_{2}}{1 + m E_{r}} + φ_{2} + μ) E_{r} \\ (λ_{P} + δ_{1} + μ + θ_{1}) I_{w} - \frac{ψ_{1} E_{w}}{1 + m E_{w}} - ξ_{1} I_{wP} \\ (λ_{P} + δ_{2} + μ + θ_{2}) I_{r} - \frac{ψ_{2} E_{r}}{1 + m E_{w}} - ξ_{2} I_{rP} \\ (ξ_{1} + δ_{3} + μ + θ_{3} + b_{3}) I_{wP} - λ_{P} I_{w} - λ_{P} I_{r} \\ (ξ_{2} + δ_{4} + μ + θ_{4}) I_{rP} - b_{3} I_{wP} \\ (φ_{3} + b_{1} + δ_{5} + μ) T_{w} - θ_{1} I_{w} - υ_{1} T_{wP} \\ (φ_{4} + δ_{6} + μ) T_{r} - θ_{2} I_{r} - υ_{2} T_{rP} - b_{1} T_{w} \\ (υ_{1} + b_{2} + δ_{7} + μ) T_{wP} - θ_{3} I_{wP} \\ (υ_{2} + δ_{8} + μ) T_{rP} - θ_{4} I_{rP} - b_{2} T_{wP} \end{matrix}]

.(13)

The Jacobian matrices associated with new infections

(F)

and transitions

(V)

at the disease-free equilibrium are given, respectively, by:

F = [\begin{matrix} 0 & 0 & {ϵ_{1} g}_{1} & 0 & g_{1} & 0 & {ϵ_{3} g}_{1} & 0 & ϵ_{2} g_{1} & 0 \\ 0 & 0 & 0 & {ϵ_{4} g}_{2} & 0 & g_{2} & 0 & {ϵ_{6} g}_{2} & 0 & {ϵ_{5} g}_{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}]

,(14)

where:

g_{1} = β_{w} [S^{0} + (1 - ρ) V^{0}]

g_{2} = β_{r} [S^{0} + (1 - ρ) V^{0}]

The matrix

V

is given by:

V = [\begin{matrix} r_{1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & r_{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ {- ψ}_{1} & 0 & r_{3} & 0 & - ξ_{1} & 0 & 0 & 0 & 0 & 0 \\ 0 & - ψ_{2} & 0 & r_{4} & 0 & - ξ_{2} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & r_{5} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - b_{3} & r_{6} & 0 & 0 & 0 & 0 \\ 0 & 0 & - θ_{1} & 0 & 0 & 0 & r_{7} & 0 & - υ_{1} & 0 \\ 0 & 0 & 0 & - θ_{2} & 0 & 0 & - b_{1} & r_{8} & 0 & - υ_{2} \\ 0 & 0 & 0 & 0 & - θ_{3} & 0 & 0 & 0 & r_{9} & 0 \\ 0 & 0 & 0 & 0 & 0 & - θ_{4} & 0 & 0 & - b_{2} & r_{10} \end{matrix}]

.(15)

where:

r_{1} = (ψ_{1} + φ_{1} + μ)

r_{2} = (ψ_{2} + φ_{2} + μ)

r_{3} = (δ_{1} + θ_{1} + μ)

r_{4} = (δ_{2} + θ_{2} + μ)

r_{5} = (ξ_{1} + δ_{3} + θ_{3} + μ + b_{3})

r_{6} = (ξ_{2} + δ_{4} + θ_{4} + μ)

r_{7} = (φ_{3} + δ_{5} + b_{1} + μ)

r_{8} = (φ_{4} + δ_{6} + μ)

r_{9} = (υ_{1} + δ_{7} + b_{2} + μ)

r_{10} = (υ_{2} + δ_{8} + μ)

The dominant eigenvalues, corresponding to the spectral radius

ρ (F V^{- 1})

, give the control reproduction numbers, expressed as:

R_{cw} = \frac{β_{w} [ϵ_{1} ψ_{1} (φ_{3} + δ_{5} + b_{1} + μ) + ϵ_{3} θ_{1} ψ_{1}] [\frac{(1 - Q) Γ}{μ} + (1 - ρ) \frac{Q Γ}{μ}]}{(ψ_{1} + φ_{1} + μ) (δ_{1} + θ_{1} + μ) (φ_{3} + δ_{5} + b_{1} + μ)}

(16)

R_{cr} = \frac{β_{r} [ϵ_{4} ψ_{2} (φ_{4} + δ_{6} + μ) + ϵ_{6} θ_{2} ψ_{2}] [\frac{(1 - Q) Γ}{μ} + (1 - ρ) \frac{Q Γ}{μ}]}{(δ_{2} + θ_{2} + μ) (ψ_{2} + φ_{2} + μ) (φ_{4} + δ_{6} + μ)}

(17)

Here,

R_{cw}

is the number of secondary drug-sensitive tuberculosis infections in the presence of opportunistic pneumonia produced by a single infected individual, while

R_{cr}

represents the number of secondary drug-resistant tuberculosis infections in the presence of opportunistic pneumonia produced by an infected individual during their entire infectious period in a population that is not fully susceptible due to existing control efforts.

3.5. Local Stability of the Disease-free Equilibrium

In this section, a small perturbation around the disease-free equilibrium point is investigated to determine whether it leads to disease prevalence for either of the two strains of tuberculosis in the presence of opportunistic pneumonia, or if the system returns to the disease-free state.

The local stability of the disease-free equilibrium, in the absence of the drug-resistant strain, is examined by employing Theorem 3 below.

Theorem 3: The disease-free equilibrium point is locally asymptotically stable if

R_{cw} < 1

and unstable if

R_{cw} > 1

Proof: To prove the local stability of the disease-free equilibrium, the Jacobian matrix of system (1), evaluated at the disease-free equilibrium point

B_{w}^{0}

, in the absence of the drug-resistant strain, is derived as follows:

J (B_{w}^{0}) = [\begin{matrix} - μ & 0 & 0 & {- y}_{1} & - y_{2} & - y_{3} & - y_{4} & κ_{1} \\ 0 & - μ & 0 & {- y}_{5} & {- y}_{6} & {- y}_{7} & - y_{8} & 0 \\ 0 & 0 & - r_{1} & y_{9} & y_{10} & y_{11} & y_{12} & 0 \\ 0 & 0 & ψ_{1} & {- r}_{3} & ξ_{1} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & {- r}_{5} & 0 & 0 & 0 \\ 0 & 0 & 0 & θ_{1} & 0 & - r_{7} & υ_{1} & 0 \\ 0 & 0 & 0 & 0 & θ_{3} & 0 & - r_{9} & 0 \\ 0 & 0 & φ_{1} & 0 & 0 & φ_{3} & 0 & {- r}_{11} \end{matrix}]

(18)

where:

y_{1} = β_{w} ϵ_{1} S^{0}

y_{2} = β_{w} S^{0}

y_{3} = β_{w} ϵ_{3} S^{0}

y_{4} = β_{w} ϵ_{2} S^{0}

y_{5} = β_{w} (1 - ρ) ϵ_{1} V^{0}

y_{6} = β_{w} (1 - ρ) V^{0}

y_{7} = β_{w} (1 - ρ) ϵ_{3} V^{0}

y_{8} = β_{w} (1 - ρ) ϵ_{2} V^{0}

y_{9} = β_{w} ϵ_{1} [S^{0} + (1 - ρ) V^{0}]

y_{10} = β_{w} [S^{0} + (1 - ρ) V^{0}]

y_{11} = β_{w} ϵ_{3} [S^{0} + (1 - ρ) V^{0}]

y_{12} = β_{w} ϵ_{2} [S^{0} + (1 - ρ) V^{0}]

The characteristic polynomial of equation (18) is given by:

(μ + λ) (μ + λ) (r_{5} + λ) (r_{9} + λ) (r_{11} + λ) (M_{1} λ^{3} + {M_{2} λ}^{2} + M_{3} λ + M_{4}) = 0

(19)

By the Routh-Hurwitz criterion, equation (19) has strictly negative roots, including:

λ_{1} = λ_{2} = - μ

λ_{3} = {- r}_{5}

λ_{4} = - r_{9}

λ_{5} = {- r}_{11}

The remaining part of the characteristic polynomial in equation (19) is:

M_{1} λ^{3} + M_{2} λ^{2} + M_{3} λ + M_{4} = 0

(20)

where

M_{1}, M_{2}, M_{3}

and

M_{4}

are given by:

M_{1} = 1 > 0

M_{2} = r_{1} + r_{3} + r_{7}

M_{3} = r_{1} r_{3} + r_{1} r_{7} + r_{3} r_{7} - ψ_{1} y_{9}

M_{4} = r_{1} r_{3} r_{7} - (r_{7} ψ_{1} y_{9} + y_{11} ψ_{1} θ_{1})

By the Routh-Hurwitz criterion, the following conditions must be satisfied for local asymptotic stability:

M_{1} > 0

M_{2} > 0

M_{3} > 0

,and

M_{4} > 0

From the condition

M_{4} > 0

, we obtain:

r_{1} r_{3} r_{7} - (r_{7} ψ_{1} y_{9} + y_{11} ψ_{1} θ_{1}) > 0

(21)

Inequality (21) can be rewritten as:

(r_{7} ψ_{1} y_{9} + y_{11} ψ_{1} θ_{1}) < r_{1} r_{3} r_{7}

(22)

which implies:

\frac{(r_{7} ψ_{1} y_{9} + y_{11} ψ_{1} θ_{1})}{r_{1} r_{3} r_{7}} < 1

(23)

Substituting the expressions:

r_{1} = (ψ_{1} + φ_{1} + μ)

r_{3} = (δ_{1} + θ_{1} + μ)

r_{7} = (φ_{3} + δ_{5} + b_{1} + μ)

y_{9} = β_{w} ϵ_{1} [S^{0} + (1 - ρ) V^{0}]

y_{11} = β_{w} ϵ_{3} [S^{0} + (1 - ρ) V^{0}]

S^{0} = \frac{(1 - Q) Γ}{μ}

, and

V^{0} = \frac{Q Γ}{μ}

into inequality (23), and rearranging yields:

\frac{β_{w} [ϵ_{1} ψ_{1} (φ_{3} + δ_{5} + b_{1} + μ) + ϵ_{3} θ_{1} ψ_{1}] [\frac{(1 - Q) Γ}{μ} + (1 - ρ) \frac{Q Γ}{μ}]}{(ψ_{1} + φ_{1} + μ) (δ_{1} + θ_{1} + μ) (φ_{3} + δ_{5} + b_{1} + μ)} < 1

(24)

The left-hand side of inequality (24) defines the control reproduction number:

R_{cw} = \frac{β_{w} [ϵ_{1} ψ_{1} (φ_{3} + δ_{5} + b_{1} + μ) + ϵ_{3} θ_{1} ψ_{1}] [\frac{(1 - Q) Γ}{μ} + (1 - ρ) \frac{Q Γ}{μ}]}{(ψ_{1} + φ_{1} + μ) (δ_{1} + θ_{1} + μ) (φ_{3} + δ_{5} + b_{1} + μ)}

(25)

Therefore, inequality (24) can be written as:

R_{cw} < 1

(26)

The disease-free equilibrium point is thus locally asymptotically stable if the control reproduction number

R_{cw}

is less than 1. This implies that if

R_{cw} < 1

, the drug-sensitive tuberculosis strain, in the presence of opportunistic pneumonia, will be eliminated from the population, provided the initial conditions of system (1) lie within the basin of attraction of the drug-sensitive tuberculosis-free equilibrium.

Similarly, it can be shown that the disease-free equilibrium of the system in the presence of the drug-resistant tuberculosis strain only is locally asymptotically stable if

R_{cr} < 1

, and unstable if

R_{cr} > 1

3.6. The Endemic Equilibrium

The endemic equilibrium refers to a stable state in which a disease persists in the population at a constant level over time. It implies that the disease remains present but does not lead to an exponential increase in the number of cases. To find the endemic equilibrium, we set the derivatives of the system (1) to zero, indicating a steady state where the number of new infections is balanced by the number of recoveries and other transitions out of the infected states. Hence, the model system (1) can be represented as:

0 = (1 - Q) Γ + κ_{1} R^{*} - (μ + λ_{r} + λ_{w}) S^{*}

0 = Q Γ - [μ + (1 - ρ) λ_{r} + (1 - ρ) λ_{w}] V^{*}

0 = λ_{w} S^{*} + (1 - ρ) λ_{w} V^{*} - (μ + \frac{ψ_{1}}{1 + m E_{w}^{*}} + φ_{1}) E_{w}^{*}

0 = λ_{r} S^{*} + (1 - ρ) λ_{r} V^{*} - (μ + \frac{ψ_{2}}{1 + m E_{r}^{*}} + φ_{2}) E_{r}^{*}

0 = \frac{ψ_{1} E_{w}^{*}}{1 + m E_{w}^{*}} + ξ_{1} I_{wP}^{*} - (μ + θ_{1} + δ_{1} + λ_{P}) I_{w}^{*}

0 = \frac{ψ_{2} E_{r}^{*}}{1 + m E_{r}^{*}} + ξ_{2} I_{rP}^{*} - (μ + δ_{2} + θ_{2} + λ_{P}) I_{r}^{*}

0 = λ_{P} I_{w}^{*} - (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) I_{wP}^{*}

,(27)

0 = λ_{P} I_{r}^{*} + b_{3} I_{wP}^{*} - (θ_{4} + μ + δ_{4} + ξ_{2}) I_{rP}^{*}

0 = θ_{1} I_{w}^{*} + υ_{1} T_{wP}^{*} - (φ_{3} + b_{1} + δ_{5} + μ) T_{w}^{*}

0 = θ_{2} I_{r}^{*} + υ_{2} T_{rP}^{*} + b_{1} T_{w}^{*} - (φ_{4} + δ_{6} + μ) T_{r}^{*}

0 = θ_{3} I_{wP}^{*} - (υ_{1} + b_{2} + δ_{7} + μ) T_{wP}^{*}

0 = θ_{4} I_{rP}^{*} + b_{2} T_{wP}^{*} - (υ_{2} + δ_{8} + μ) T_{rP}^{*}

0 = φ_{1} E_{w}^{*} + φ_{2} E_{r}^{*} + φ_{3} T_{w}^{*} + φ_{4} T_{r}^{*} - (μ + κ_{1}) R^{*}

The steady-state solution for system (27) is given by:

B^{*} = (S^{*}, V^{*}, E_{w}^{*}, E_{r}^{*}, I_{w}^{*}, I_{r}^{*}, I_{wP}^{*}, I_{rP}^{*}, T_{w}^{*}, T_{r}^{*}, T_{wP}^{*} {, T_{rP}^{*}, R}^{*})

where:

S^{*} = \frac{(1 - Q) Γ + κ_{1} R^{*}}{λ_{w}^{* *} + λ_{r}^{* *} + μ}

V^{*} = \frac{Q Γ}{(1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ}

E_{w}^{*} = \frac{λ_{w}^{* *} [((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ]}{((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) (μ + χ_{1} + φ_{1})}

E_{r}^{*} = \frac{λ_{r}^{* *} [((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ]}{((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) (μ + χ_{2} + φ_{2})}

I_{w}^{*} = \frac{(μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) {χ_{1} λ}_{w}^{* *} [((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ]}{(μ + χ_{1} + φ_{1}) {A_{5} D}_{10}}

I_{r}^{*} = \frac{(θ_{4} + μ + δ_{4} + ξ_{2}) {χ_{2} λ}_{r}^{* *} [((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ]}{(μ + χ_{2} + φ_{2}) A_{5} D_{11}}

I_{wP}^{*} = \frac{λ_{P}^{* *} (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) {χ_{1} λ}_{w}^{* *} [((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ]}{(μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) (μ + χ_{1} + φ_{1}) A_{5} D_{10}}

I_{rP}^{*} = \frac{λ_{P}^{* *} (θ_{4} + μ + δ_{4} + ξ_{2}) {χ_{2} λ}_{r}^{* *} [((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ]}{(θ_{4} + μ + δ_{4} + ξ_{2}) (μ + χ_{2} + φ_{2}) A_{5} D_{11}}

T_{w}^{*} = \frac{(μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) {χ_{1} λ}_{w}^{* *} [((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ] D_{13}}{(φ_{3} + b_{1} + δ_{5} + μ) A_{5} D_{10} D_{12}}

T_{r}^{*} = \frac{(θ_{4} + μ + δ_{4} + ξ_{2}) [θ_{2} D_{15} + υ_{2} (φ_{3} + b_{1} + δ_{5} + μ) λ_{P}^{* *} D_{14} + b_{1} D_{16}] D_{9}}{(φ_{3} + b_{1} + δ_{5} + μ) (μ + χ_{2} + φ_{2}) A_{5} D_{10} D_{11} D_{12} D_{17}}

T_{wP}^{*} = \frac{θ_{3} λ_{P}^{* *} (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) {χ_{1} λ}_{w}^{* *} [((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ]}{(μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) (μ + χ_{1} + φ_{1}) (υ_{1} + b_{2} + δ_{7} + μ) A_{5} D_{10}}

T_{rP}^{*} = \frac{λ_{P}^{* *} (θ_{4} + μ + δ_{4} + ξ_{2}) D_{14}}{(θ_{4} + μ + δ_{4} + ξ_{2}) (μ + χ_{2} + φ_{2}) A_{5} D_{10} D_{11} D_{12} (υ_{2} + δ_{8} + μ)}

R^{*} = \frac{D_{19} [{λ_{w}^{* *} (μ + χ_{2} + φ_{2}) φ}_{1} A_{6} + {λ_{r}^{* *} (μ + χ_{1} + φ_{1}) φ}_{2} A_{6} + λ_{w}^{* *} φ_{3} A_{7} + φ_{4} A_{8}]}{(κ_{1} + μ) (μ + χ_{1} + φ_{1}) (μ + χ_{2} + φ_{2}) {(φ_{3} + b_{1} + δ_{5} + μ) A}_{5} D_{10} D_{11} D_{17} D_{18}}

where:

A_{5} = ((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) (λ_{w}^{* *} + λ_{r}^{* *} + μ)

A_{6} = (φ_{3} + b_{1} + δ_{5} + μ) D_{10} D_{11} D_{17} D_{18}

A_{7} = (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) χ_{1} (μ + χ_{2} + φ_{2}) D_{11} {D_{13} D}_{17}

A_{8} = (θ_{4} + μ + δ_{4} + ξ_{2}) [θ_{2} D_{15} + υ_{2} (φ_{3} + b_{1} + δ_{5} + μ) λ_{P}^{* *} D_{14} + b_{1} D_{16}]

D_{9} = ((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) ((1 - Q) Γ + κ_{1} R^{*}) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ

D_{10} = (μ + θ_{1} + δ_{1} + λ_{P}^{* *}) (μ + θ_{3} + δ_{3} + ξ_{1}) - λ_{P}^{* *} ξ_{1}

D_{11} = (μ + θ_{2} + δ_{2} + λ_{P}^{* *}) (θ_{4} + μ + δ_{4} + ξ_{2}) - λ_{P}^{* *} ξ_{2}

D_{12} = (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) (μ + χ_{1} + φ_{1}) (υ_{1} + b_{2} + δ_{7} + μ)

D_{13} = [θ_{1} (μ + θ_{3} + δ_{3} + ξ_{1}) (υ_{1} + b_{2} + δ_{7} + μ) + υ_{1} θ_{3} λ_{P}^{* *}]

D_{14} = θ_{4} {χ_{2} λ}_{r}^{* *} D_{10} D_{12} + b_{2} (μ + χ_{2} + φ_{2}) D_{11} θ_{3} (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) {χ_{1} λ}_{w}^{* *}

D_{15} = (φ_{3} + b_{1} + δ_{5} + μ) (υ_{2} + δ_{8} + μ) D_{10} D_{12} (θ_{4} + μ + δ_{4} + ξ_{2}) {χ_{2} λ}_{r}^{* *}

D_{16} = (μ + χ_{2} + φ_{2}) {(υ_{2} + δ_{8} + μ) (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) {χ_{1} λ}_{w}^{* *} D_{13} D}_{11}

D_{17} = (θ_{4} + μ + δ_{4} + ξ_{2}) (υ_{2} + δ_{8} + μ) (φ_{4} + δ_{6} + μ)

D_{18} = (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) (υ_{1} + b_{2} + δ_{7} + μ)

D_{19} = ((1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *}) + μ) (S^{*} (λ_{w}^{* *} + λ_{r}^{* *} + μ)) + (1 - ρ) (λ_{w}^{* *} + λ_{r}^{* *} + μ) Q Γ

χ_{1} = \frac{ψ_{1}}{1 + m E_{w}^{*}}

χ_{2} = \frac{ψ_{2}}{1 + m E_{r}^{*}}

3.7. Loal Stability of the Endemic Equilibrium

The endemic equilibrium of system (1), in the absence of drug-resistant tuberculosis, is given by:

B_{w}^{*} = (S^{*}, V^{*}, E_{w}^{*}, I_{w}^{*}, I_{wP}^{*}, T_{w}^{*}, T_{wP}^{*} {, R}^{*})

The local stability of this endemic equilibrium is analyzed by employing Theorem 5.

Theorem 5: A positive endemic equilibrium exists and is locally asymptotically stable whenever

R_{cw} > 1

Proof: The drug-sensitive tuberculosis strain becomes endemic if:

\frac{d E_{w}}{dt} > 0

\frac{d I_{w}}{dt} > 0

\frac{d I_{wP}}{dt} > 0

\frac{d T_{w}}{dt} > 0

\frac{d T_{wP}}{dt} > 0

Considering

\frac{d E_{w}}{dt} > 0

, and letting the natural immunity parameter

m

approach zero at the endemic equilibrium, gives:

\frac{d E_{w}}{dt} = λ_{w} S + (1 - ρ) λ_{w} V - (μ + ψ_{1} + φ_{1}) E_{w} > 0

(28)

Inequality (28) can be rewritten as:

E_{w} < \frac{λ_{w} S + (1 - ρ) λ_{w} V}{(μ + ψ_{1} + φ_{1})}

(29)

The endemicity of the coexistence of drug-sensitive tuberculosis and pneumonia is primarily driven by the transmission of drug-sensitive TB, since pneumonia is considered an opportunistic infection. Therefore, the force of infection for drug-sensitive TB from tuberculosis patients only is considered, and is given by:

λ_{w} = β_{w} (ϵ_{1} I_{w} + ϵ_{3} T_{w})

(30)

Substituting equation (30) into inequality (29) gives:

E_{w} < \frac{β_{w} ({ϵ_{1} I}_{w} + ϵ_{3} T_{w}) [S + (1 - ρ) V]}{(μ + ψ_{1} + φ_{1})}

(31)

From the steady-state system (27), the expressions for

I_{w}^{*}

and

T_{w}^{*}

, representing the population spreading of drug-sensitive TB only, are respectively:

I_{w}^{*} = \frac{ψ_{1} E_{w}^{*}}{(θ_{1} + μ + δ_{1})},

(32)

T_{w}^{*} = \frac{θ_{1} I_{w}^{*}}{(φ_{3} + μ + δ_{5} + b_{1})}

.(33)

Substituting equations (32) and (33) into inequality (31) yields:

E_{w}^{*} < \frac{β_{w} [\frac{ϵ_{1} ψ_{1} E_{w}^{*}}{θ_{1} + μ + δ_{1}} + \frac{ϵ_{3} θ_{1} ψ_{1} E_{w}^{*}}{(φ_{3} + μ + δ_{5} + b_{1}) (θ_{1} + μ + δ_{1})}] [S + (1 - ρ) V]}{μ + ψ_{1} + φ_{1}}

(34)

Simplifying inequality (34) leads to:

1 < \frac{β_{w} [(ϵ_{1} ψ_{1} (φ_{3} + μ + δ_{5} + b_{1}) + ϵ_{3} θ_{1} ψ_{1}) \begin{matrix}  \end{matrix}] [\frac{(1 - Q) Γ}{μ} + (1 - ρ) \frac{QΓ}{μ}]}{(φ_{3} + μ + δ_{5} + b_{1}) (θ_{1} + μ + δ_{1}) (μ + ψ_{1} + φ_{1})}

(35)

The right-hand side of inequality (35) defines the control reproduction number

R_{cw}

, given by:

R_{cw} = \frac{β_{w} [(ϵ_{1} ψ_{1} (φ_{3} + μ + δ_{5} + b_{1}) + ϵ_{3} θ_{1} ψ_{1}) \begin{matrix}  \end{matrix}] [\frac{(1 - Q) Γ}{μ} + (1 - ρ) \frac{Q Γ}{μ}] E_{w}^{*}}{(φ_{3} + μ + δ_{5} + b_{1}) (θ_{1} + μ + δ_{1}) (μ + ψ_{1} + φ_{1})}

(36)

Therefore, inequality (35) can be written as follows:

1 < R_{cw}

\Rightarrow R_{cw} > 1

(37)

Thus, when the reproduction number

R_{cw} > 1

, the system reaches an endemic equilibrium. This equilibrium occurs when each infected individual causes, on average, more than one new infection in a partially susceptible population, leading to the sustained transmission of the disease.

Similarly, it can be shown that the local stability of the endemic equilibrium in the system with only the drug-resistant strain of tuberculosis, given by

B_{r}^{*} = (S^{*}, V^{*}, E_{r}^{*}, I_{r}^{*} {, I}_{rP}^{*}, T_{r}^{*} {, T_{rP}^{*}, R}^{*})

is locally asymptotically stable if

R_{cr} > 1

and unstable if

R_{cr} < 1

3.8. Sensitivity Analysis of the Control Reproduction Numbers

This section aims to determine the parameters that fuel the transmission and prevalence of pulmonary tuberculosis, considering the presence of opportunistic pneumonia in the population. This is done by performing sensitivity analysis of the reproduction numbers relative to the model parameters. The normalized forward sensitivity index method, as described in

[27]

, has been used to determine the sensitivity analysis of the control reproduction numbers as follows:

Υ_{K}^{R_{c}} = \frac{\partial R_{c}}{\partial K} \times \frac{K}{R_{c}}

(38)

Here,

K

is the parameter of interest, while

R_{c}

is the control reproduction number.

The parameter values in Table 3 are used to calculate the sensitivity indices of

R_{Cw}

with respect to the parameters

β_{w}

ϵ_{1}

ϵ_{3}

ψ_{1}

θ_{1}

φ_{1}

φ_{3}

ρ

b_{1}

and sensitivity indices of

R_{Cr}

with respect to the parameters

β_{r}

ϵ_{4}

ϵ_{6}

ψ_{2}

θ_{2}

φ_{2}

φ_{4}

ρ

The computed sensitivity indices are presented in Table 4. A positive sensitivity index indicates that the reproduction number increases with an increase in the corresponding parameter, while a negative sensitivity index indicates that the reproduction number decreases as the parameter increases

[28]

. From Table 4, it is observed that the parameters

β_{w}

ϵ_{1}

ψ_{1}

, and

ϵ_{3}

have positive sensitivity index values, indicating that an increase in these parameters leads to a corresponding increase in the number of individuals co-infected with drug-sensitive tuberculosis and pneumonia. Moreover, the parameters

θ_{1}

φ_{1}

φ_{3}

ρ

, and

b_{1}

have negative sensitivity index values, meaning that increasing these parameters results in a decrease in the number of individuals co-infected with drug-sensitive tuberculosis and pneumonia.

On the other hand, increasing the parameters

β_{r}

ϵ_{4}

ψ_{2}

, and

ϵ_{6}

leads to a corresponding increase in the number of individuals co-infected with drug-resistant tuberculosis and pneumonia, while increasing the parameters

φ_{2}

θ_{2}

φ_{4}

, and

ρ

leads to a corresponding decrease in the number of individuals co-infected with drug-resistant tuberculosis and pneumonia.

Table 3. Baseline parameter values used in the tuberculosis-pneumonia co-infection model.

Parameter

Value

Reference

$Γ$

1508563 ${year}^{- 1}$

[

29 ]

$Q$

0.8 ${year}^{- 1}$

[

29 ]

$ρ$

0.6

[1]

$μ$

0.0057 ${year}^{- 1}$

[

29 ]

$ψ_{1}$ , $ψ_{2}$

0.7, 0.3 ${year}^{- 1}$

0 ]

$θ_{1}$ , $θ_{2}$

0.6, 0.8 ${year}^{- 1}$

0 ]

$θ_{3}$ , $θ_{4}$

0.023, 0.053 ${year}^{- 1}$

Data fitted

$m$

0.0003 ${year}^{- 1}$

Data fitted

$δ_{1}$ , $δ_{2}$

0.006, 0.06 ${year}^{- 1}$

1 ]

$δ_{3}$

0.005 ${year}^{- 1}$

Data fitted

$δ_{4}$

0.03 ${year}^{- 1}$

Data fitted

$δ_{5}$ , $δ_{6}$

0.06, 0.007 ${year}^{- 1}$

[31]

$δ_{7}$ , $δ_{8}$

0.00031, 0.00011 ${year}^{- 1}$

Data fitted

$φ_{1}$ , $φ_{2}$

0.34,0.05 ${year}^{- 1}$

1 ]

$φ_{3}$ , $φ_{4}$

0.75,0.7 ${year}^{- 1}$

1 ]

$b_{1}$

0.08 ${year}^{- 1}$

Data fitted

$b_{2}$

0.065 ${year}^{- 1}$

Data fitted

$b_{3}$

0.035 ${year}^{- 1}$

Data fitted

$ν_{1}$ , $ν_{2}$

0.001, 0.00219 ${year}^{- 1}$

Data fitted

$κ_{1}$

0.003 ${year}^{- 1}$

2 ]

$β_{w}$

$2.2258 \times 10^{- 7} {year}^{- 1}$

1 ]

$β_{r}$

$1.78 \times 10^{- 8} {year}^{- 1}$

1 ]

$β_{P}$

$2.79684 \times 10^{- 7} {year}^{- 1}$

1 ]

$ξ_{1}$ , $ξ_{2}$

0.002, 0.003 ${year}^{- 1}$

Data fitted

$ϵ_{1}$ , $ϵ_{2}$ , $ϵ_{3}$

0.3704, 0.1655, 0.1527

Data fitted

$ϵ_{4}$ , $ϵ_{5}$ , $ϵ_{6}$

0.1487, 0.1389, 0.1054

Data fitted

$ϵ_{7}$ , $ϵ_{8}$ , $ϵ_{9}$

0.0315, 0.0576, 0.6012

Data fitted

The sensitivity indices from Table 4 further elucidate that the parameters with the greatest effects on the spread of both drug-sensitive and drug-resistant tuberculosis-pneumonia co-infections are the transmission rate and vaccine efficacy. For instance, a 10% increase in transmission rates results in a 10% increase in the number of new co-infections from index cases, while a 10% increase in vaccine efficacy leads to a 9.23077% decrease in the number of new co-infections from index cases.

Table 4. Sensitivity indices of the reproduction numbers with respect to selected model parameters.

Parameter	Sensitivity index
Sensitivity indices of $R_{cw}$
$ρ$	$- 0.923077$
$θ_{1}$	$- 0.742107$
$φ_{1}$	$- 0.325141$
$φ_{3}$	$- 0.212753$
$b_{1}$	$- 0.022693$
$ϵ_{3}$	$+ 0.238766$
$ψ_{1}$	$+ 0.330592$
$ϵ_{1}$	$+ 0.761234$
$β_{w}$	+1
Sensitivity indices of $R_{cr}$
$ρ$	$- 0.923077$
$θ_{2}$	$- 0.804797$
$φ_{2}$	$- 0.3211300$
$φ_{4}$	$- 0.0888458$
$ϵ_{6}$	$+ 0.0904577$
$ψ_{2}$	$+ 0.357739$
$ϵ_{4}$	$+ 0.909542$
$β_{r}$	+1

4. Extension of the Model to Optimal Control

In this study, mathematical optimal control theory is applied to determine the most effective strategies for combating the spread of pulmonary tuberculosis. A model of pulmonary tuberculosis in the presence of opportunistic pneumonia, incorporating a drug-resistant strain, is extended to include time-dependent control strategies. These strategies target the prevention of the drug-sensitive TB strain, the drug-resistant TB strain, and opportunistic pneumonia, corresponding to the control functions

u_{1} (t)

u_{2} (t)

, and

u_{3} (t)

, respectively. Additionally, the time-dependent controls

u_{4} (t)

and

u_{5} (t)

represent screening interventions for drug-sensitive latent TB and drug-resistant latent TB, respectively, while

u_{6} (t)

and

u_{7} (t)

correspond to treatment efforts for the drug-sensitive and drug-resistant TB strains, respectively.

The prevention of the drug-sensitive TB strain, represented by control

u_{1} (t)

, involves efforts to reduce the number of susceptible individuals becoming infected through community outreach activities, including awareness campaigns, contact tracing, and targeted active TB case detection. The prevention of the drug-resistant TB strain, represented by control

u_{2} (t)

, involves providing social support for drug-resistant patients to optimize treatment outcomes, tracing contacts of drug-resistant TB cases, conducting drug susceptibility testing for TB patients, and monitoring individuals to ensure the completion of the TB treatment course. The prevention of opportunistic pneumonia among TB patients, represented by control

u_{3} (t)

, entails providing TB patient-centered care to diagnose co-morbidities and raising awareness about the co-occurrence of TB with other lung diseases. Control

u_{4} (t)

refers to the case-finding strategy for individuals with drug-sensitive latent TB, while

u_{5} (t)

refers to the strategy for identifying individuals with drug-resistant latent TB. Controls

u_{6} (t)

and

u_{7} (t)

represent the medical treatment coverage strategies for active drug-sensitive and drug-resistant TB strains, respectively. Therefore, incorporating the controls

u_{i} (t)

, where

i =

1,2,3,…,7, into the model yields the following optimal control system:

\frac{dS}{dt} = (1 - Q) Γ + κ_{1} R - ((1 - u_{1}) λ_{w} + {(1 - u_{2}) λ}_{r} + μ) S

\frac{dV}{dt} = Q Γ - [(1 - u_{1}) (1 - ρ) λ_{w} + (1 - u_{2}) (1 - ρ) λ_{r} + μ] V

\frac{d E_{w}}{dt} = {(1 - u_{1}) λ}_{w} S + (1 - u_{1}) (1 - ρ) λ_{w} V - (μ + \frac{ψ_{1}}{1 + m E_{w}} + (φ_{1} + u_{4})) E_{w}

\frac{d E_{r}}{dt} = {(1 - u_{2}) λ}_{r} S + (1 - u_{2}) (1 - ρ) λ_{r} V - (μ + \frac{ψ_{2}}{1 + m E_{r}} + (φ_{2} + u_{5})) E_{r}

\frac{d I_{w}}{dt} = \frac{ψ_{1} E_{w}}{1 + m E_{w}} + ξ_{1} I_{wP} - (μ + (θ_{1} + u_{6}) + δ_{1} + {(1 - u_{3}) λ}_{P}) I_{w}

\frac{d I_{r}}{dt} = \frac{ψ_{2} E_{r}}{1 + m E_{r}} + ξ_{2} I_{rP} - (μ + δ_{2} + (θ_{2} + u_{7}) + {(1 - u_{3}) λ}_{P}) I_{r}

,(39)

\frac{d I_{wP}}{dt} = {(1 - u_{3}) λ}_{P} I_{w} - (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3}) I_{wP}

\frac{d I_{rP}}{dt} = {(1 - u_{3}) λ}_{P} I_{r} + b_{3} I_{wP} - (θ_{4} + μ + δ_{4} + ξ_{2}) I_{rP}

\frac{{dT}_{w}}{dt} = (θ_{1} + u_{6}) I_{w} + υ_{1} T_{wP} - (φ_{3} + b_{1} + δ_{5} + μ) T_{w}

\frac{{dT}_{r}}{dt} = (θ_{2} + u_{7}) I_{r} + υ_{2} T_{rP} + b_{1} T_{w} - (φ_{4} + δ_{6} + μ) T_{r}

\frac{{dT}_{wP}}{dt} = θ_{3} I_{wP} - (υ_{1} + b_{2} + δ_{7} + μ) T_{wP}

\frac{{dT}_{rP}}{dt} = θ_{4} I_{rP} + b_{2} T_{wP} - (υ_{2} + δ_{8} + μ) T_{rP}

\frac{dR}{dt} = (φ_{1} + u_{4}) E_{w} + (φ_{2} + u_{5}) E_{r} + φ_{3} T_{w} + φ_{4} T_{r} - (μ + κ_{1}) R

The study aims to minimize the number of infected individuals while keeping the cost of intervention measures as low as possible. To optimize the control strategies, a bounded Lebesgue measurable control set is introduced and defined as follows:

U = \{u_{1} (t), u_{2} (t), \dots, u_{7} (t) : u_{i \min} \leq u_{i} \leq u_{i \max}, t \in [0, T], i = 1, 2, \dots, 7\}

(40)

The lower bounds

u_{i \min}

(for

i =

1, 2,…,7) represent suboptimal existing control strategies, while the upper bounds

u_{i \max}

represent the maximum feasible efforts to minimize the infected population and associated costs. The optimal levels of effort for controlling the infections are determined by minimizing the following objective functional:

J (u_{1}, u_{2}, u_{3}, u_{4}, u_{5}, u_{6}, u_{7}) = \int_{0}^{T} M dt

(41)

where:

{M = B}_{1} E_{w} + B_{2} E_{r} + B_{3} I_{w} + B_{4} I_{r} + B_{5} I_{wP}

+ B_{6} I_{rP} + B_{7} T_{w} + B_{8} T_{r} + B_{9} T_{wP} + B_{10} T_{rP} + \frac{1}{2} \sum_{i = 1}^{7} w_{i} u_{i}^{2}

The coefficients

B_{1}

B_{2}

,…,

B_{10}

represent the costs associated with minimizing the infected populations, while the term

\frac{1}{2} \sum_{i = 1}^{7} w_{i} u_{i}^{2}

represents the costs associated with implementing the controls

u_{i}

. Additionally, the fixed constant

T

denotes the final time of intervention. The cost expression is quadratic because the cost of control interventions increases non-linearly with the intensity of the intervention. Quadratic functions effectively capture such real-world cost structures, where expenses escalate faster than the applied control effort.

Thus, we determining the optimal control

u^{*} = (u_{1}^{*}, u_{2}^{*}, \dots, u_{7}^{*})

such that:

J (u^{*}) =

min

\{J (u) : u = (u_{1}, u_{2}, u_{3}, u_{4}, u_{5}, u_{6}, u_{7}) \in U\}

.(42)

The necessary condition for the existence of an optimal control is determined by applying Pontryagin’s Maximum Principle, which transforms the model system of equations (39) and the minimization functional (41) into a problem of minimizing the pointwise Hamiltonian with respect to the control functions

u_{1}

u_{2}

u_{3}

u_{4}

u_{5}

u_{6}

u_{7}

[33]

. Accordingly, the Hamiltonian

H

is defined using Pontryagin’s Maximum Principle as follows:

H = M + P_{1} (t) \{(1 - Q) Γ + κ_{1} R - ((1 - u_{1}) λ_{w} + {(1 - u_{2}) λ}_{r} + μ) S\}

{+ P}_{2} (t) \{Q Γ - ((1 - u_{1}) (1 - ρ) λ_{w} + (1 - u_{2}) (1 - ρ) λ_{r} + μ) V\}

+ P_{3} (t) \{{(1 - u_{1}) λ}_{w} S + (1 - u_{1}) (1 - ρ) λ_{w} V - (μ + \frac{ψ_{1}}{1 + m E_{w}} + (φ_{1} + u_{4})) E_{w}\}

+ P_{4} (t) \{{(1 - u_{2}) λ}_{r} S + (1 - u_{2}) (1 - ρ) λ_{r} V - (μ + \frac{ψ_{2}}{1 + m E_{r}} + (φ_{2} + u_{5})) E_{r}\}

+ P_{5} (t) \{\frac{ψ_{1} E_{w}}{1 + m E_{w}} + ξ_{1} I_{wP} - (μ + (θ_{1} + u_{6}) + δ_{1} + {(1 - u_{3}) λ}_{P}) I_{w}\}

+ P_{6} (t) \{\frac{ψ_{2} E_{r}}{1 + m E_{r}} + ξ_{2} I_{rP} - (μ + δ_{2} + (θ_{2} + u_{7}) + {(1 - u_{3}) λ}_{P}) I_{r}\}

+ P_{7} (t) \{{(1 - u_{3}) λ}_{P} I_{w} - (μ + θ_{3} + δ_{3} + ξ_{1} {+ b}_{3}) I_{wP}\}

(43)

+ P_{8} (t) \{{(1 - u_{3}) λ}_{P} I_{r} + b_{3} I_{wP} - (θ_{4} + μ + δ_{4} + ξ_{2}) I_{rP}\}

+ P_{9} (t) \{(θ_{1} + u_{6}) I_{w} + υ_{1} T_{wP} - (φ_{3} + b_{1} + δ_{5} + μ) T_{w}\}

+ P_{10} (t) \{(θ_{2} + u_{7}) I_{r} + υ_{2} T_{rP} + b_{1} T_{w} - (φ_{4} + δ_{6} + μ) T_{r}\}

+ P_{11} (t) \{θ_{3} I_{wP} - (υ_{1} + b_{2} + δ_{7} + μ) T_{wP}\}

+ P_{12} (t) \{θ_{4} I_{rP} + b_{2} T_{wP} - (υ_{2} + δ_{8} + μ) T_{rP}\}

+ P_{13} (t) \{(φ_{1} + u_{4}) E_{w} + (φ_{2} + u_{5}) E_{r} + φ_{3} T_{w} + φ_{4} T_{r} - (μ + κ_{1}) R\}

where:

P_{i} (t)

i = 1

2

,…,13, are the corresponding adjoint or costate variables, determined by differentiating the Hamiltonian

H

, given by equation (43), with respect to state variables

S, V, E_{w}, E_{r}, I_{w}, I_{r}, I_{wP}, I_{rP}, T_{w}, T_{r}, T_{wP}, T_{rP}, R .

For an optimal control

u^{*} = (u_{1}^{*}, u_{2}^{*}, \dots, u_{7}^{*})

that minimizes

J (u_{1}, u_{2}, u_{3}, u_{4}, u_{5}, u_{6}, u_{7})

over

U

, there exist adjoint variables

P_{i}

, for

i = 1

,…,13, such that:

\frac{d P_{1}}{dt} = - \frac{\partial H}{\partial S}

\frac{d P_{2}}{dt} = - \frac{\partial H}{\partial V}

\frac{d P_{3}}{dt} = - \frac{\partial H}{\partial E_{w}}

\frac{d P_{4}}{dt} = - \frac{\partial H}{\partial E_{r}}

\frac{d P_{5}}{dt} = - \frac{\partial H}{\partial I_{w}}

\frac{d P_{6}}{dt} = - \frac{\partial H}{\partial I_{r}}

\frac{d P_{7}}{dt} = - \frac{\partial H}{\partial I_{wP}}

\frac{d P_{8}}{dt} = - \frac{\partial H}{\partial I_{rP}}

\frac{d P_{9}}{dt} = - \frac{\partial H}{\partial T_{w}}

\frac{d P_{10}}{dt} = - \frac{\partial H}{\partial T_{r}}

\frac{d P_{11}}{dt} = - \frac{\partial H}{\partial T_{wP}}

\frac{d P_{12}}{dt} = - \frac{\partial H}{\partial T_{rP}}

\frac{d P_{13}}{dt} = - \frac{\partial H}{\partial R}

.(44)

Evaluating system (44) yields the following adjoint system:

\frac{d P_{1}}{dt} = P_{1} ((1 - u_{1}) λ_{w} + (1 - u_{2}) λ_{r} + μ) - P_{3} (1 - u_{1}) λ_{w} - P_{4} (1 - u_{2}) λ_{r}

\frac{d P_{2}}{dt} = P_{2} ((1 - u_{1}) {(1 - ρ) λ}_{w} + (1 - u_{2}) (1 - ρ) λ_{r} + μ) - P_{3} (1 - u_{1}) (1 - ρ) λ_{w} {- P}_{4} (1 - u_{2}) (1 - ρ) λ_{r}

\frac{d P_{3}}{dt} = - B_{1} + P_{3} (\frac{ψ_{1}}{{(1 + m E_{w})}^{2}} + φ_{1} + u_{4} + μ) - P_{5} \frac{ψ_{1}}{{(1 + m E_{w})}^{2}} - P_{13} (φ_{1} + u_{4})

\frac{d P_{4}}{dt} = - B_{2} + P_{4} (\frac{ψ_{2}}{{(1 + m E_{r})}^{2}} + φ_{2} + u_{5} + μ) - P_{6} \frac{ψ_{2}}{{(1 + m E_{r})}^{2}} - P_{13} (φ_{2} + u_{5})

,(45)

\frac{d P_{5}}{dt} = - B_{3} + P_{1} ({1 - u}_{1}) β_{w} ϵ_{1} S + P_{2} ({1 - u}_{1}) (1 - ρ) β_{w} ϵ_{1} V - P_{9} (θ_{1} + u_{6})

- P_{3} (({1 - u}_{1}) β_{w} ϵ_{1} S + ({1 - u}_{1}) (1 - ρ) β_{w} ϵ_{1} V)

{+ P}_{5} (θ_{1} + u_{6} + δ_{1} + {(1 - u_{3}) λ}_{P} + μ) - P_{7} (1 - u_{3}) λ_{P}

\frac{d P_{6}}{dt} = - B_{4} + P_{1} ({1 - u}_{2}) β_{r} ϵ_{4} S + P_{2} ({1 - u}_{2}) (1 - ρ) β_{r} ϵ_{4} V

- P_{10} (θ_{2} + u_{7}) - P_{4} (({1 - u}_{2}) β_{r} ϵ_{4} S + ({1 - u}_{2}) (1 - ρ) β_{r} ϵ_{4} V)

{+ P}_{6} (θ_{2} + u_{7} + δ_{2} + {(1 - u_{3}) λ}_{P} + μ) - P_{8} (1 - u_{3}) λ_{P}

\frac{d P_{7}}{dt} = - B_{5} + P_{1} ({1 - u}_{1}) β_{w} S + P_{2} ({1 - u}_{1}) (1 - ρ) β_{w} V - P_{11} θ_{3}

+ P_{6} (1 - u_{3}) β_{P} I_{r} {- P}_{3} (({1 - u}_{1}) β_{w} S + ({1 - u}_{1}) (1 - ρ) β_{w} V) +

P_{5} (({1 - u}_{3}) β_{P} I_{w} - ξ_{1}) {- P}_{7} ({(1 - u_{3}) β}_{P} I_{w} - (μ + θ_{3} + δ_{3} + ξ_{1} + b_{3})) - P_{8} ((1 - u_{3}) β_{P} I_{r} - b_{3})

\frac{d P_{8}}{dt} = - B_{6} + P_{1} ({1 - u}_{2}) β_{r} S + P_{2} ({1 - u}_{2}) (1 - ρ) β_{r} V - P_{12} θ_{4} - P_{7} (1 - u_{3}) β_{P} ϵ_{7} I_{w}

- P_{4} (({1 - u}_{2}) β_{r} S + ({1 - u}_{2}) (1 - ρ) β_{r} V) + P_{6} (({1 - u}_{3}) β_{P} ϵ_{7} I_{r} - ξ_{2})

- P_{8} ({(1 - u_{3}) β}_{P} ϵ_{7} I_{r} - (μ + θ_{4} + δ_{4} + ξ_{2})) + P_{5} (1 - u_{3}) β_{P} ϵ_{7} I_{w}

\frac{d P_{9}}{dt} = - B_{7} + P_{1} ({1 - u}_{1}) β_{w} ϵ_{3} S + P_{2} ({1 - u}_{1}) (1 - ρ) β_{w} ϵ_{3} V - P_{10} b_{1} - P_{13} φ_{3}

- P_{3} (({1 - u}_{1}) β_{w} ϵ_{3} S + ({1 - u}_{1}) (1 - ρ) β_{w} ϵ_{3} V) {+ P}_{9} (φ_{3} + b_{1} + δ_{5} + μ)

\frac{d P_{10}}{dt} = - B_{8} + P_{1} ({1 - u}_{2}) β_{r} ϵ_{6} S + P_{2} ({1 - u}_{2}) (1 - ρ) β_{r} ϵ_{6} V - P_{13} φ_{4}

- P_{4} (({1 - u}_{2}) β_{r} ϵ_{6} S + ({1 - u}_{2}) (1 - ρ) β_{r} ϵ_{6} V) {+ P}_{10} (φ_{4} + δ_{6} + μ)

\frac{d P_{11}}{dt} = - B_{9} + P_{1} ({1 - u}_{1}) β_{w} ϵ_{2} S + P_{2} ({1 - u}_{1}) (1 - ρ) β_{w} ϵ_{2} V {- P}_{3} (({1 - u}_{1}) β_{w} ϵ_{2} S + ({1 - u}_{1}) (1 - ρ) β_{w} ϵ_{2} V)

+ P_{5} ({1 - u}_{3}) β_{P} {ϵ_{8} I}_{w} {+ P}_{6} {(1 - u_{3}) β}_{P} ϵ_{8} I_{r} - P_{7} (1 - u_{3}) β_{P} ϵ_{8} I_{w} {- P}_{8} {(1 - u_{3}) β}_{P} ϵ_{8} I_{r}

- P_{9} υ_{1} + P_{11} (υ_{1} + b_{2} + δ_{7} + μ) - P_{12} b_{2}

\frac{d P_{12}}{dt} = - B_{10} + P_{1} ({1 - u}_{2}) β_{r} ϵ_{5} S + P_{2} ({1 - u}_{2}) (1 - ρ) β_{r} ϵ_{5} V {- P}_{4} (({1 - u}_{2}) β_{r} ϵ_{5} S + ({1 - u}_{2}) (1 - ρ) β_{r} ϵ_{5} V)

+ P_{5} ({1 - u}_{3}) β_{P} {ϵ_{9} I}_{w} {+ P}_{6} {(1 - u_{3}) β}_{P} ϵ_{9} I_{r} - P_{7} (1 - u_{3}) β_{P} ϵ_{9} I_{w} {- P}_{8} {(1 - u_{3}) β}_{P} ϵ_{9} I_{r} - P_{10} υ_{2} + P_{12} (υ_{2} + δ_{8} + μ)

\frac{d P_{13}}{dt} = - P_{1} κ_{1} + P_{13} (μ + κ_{1})

with the transversality conditions

P_{i} (T) = 0

, for

i = 1

,…,13.

The optimal controls

u_{1}^{*}, u_{2}^{*}, \dots, u_{7}^{*}

are obtained by minimizing the Hamiltonian

H

with respect to the control variables

u = (u_{1}, u_{2}, u_{3}, u_{4}, u_{5}, u_{6}, u_{7}) \in U

, such that

\frac{\partial H}{\partial u_{i}} = 0

. Thus, the following set of optimality conditions is obtained:

\frac{\partial H}{\partial u_{1}} = P_{1} λ_{w} S + P_{2} (1 - ρ) λ_{w} V - P_{3} λ_{w} S - P_{3} (1 - ρ) λ_{w} V + w_{1} u_{1} = 0

\frac{\partial H}{\partial u_{2}} = P_{1} λ_{r} S + P_{2} (1 - ρ) λ_{r} V - P_{4} λ_{r} S - P_{4} {(1 - ρ) λ}_{r} V + w_{2} u_{2} = 0

\frac{\partial H}{\partial u_{3}} = P_{5} λ_{P} I_{w} + P_{6} λ_{P} I_{r} - P_{7} λ_{P} I_{w} - P_{8} λ_{P} I_{r} + w_{3} u_{3} = 0

\frac{\partial H}{\partial u_{4}} = {- P}_{3} E_{w} + P_{13} E_{w} + w_{4} u_{4} = 0

,(46)

\frac{\partial H}{\partial u_{5}} = {- P}_{4} E_{r} + P_{13} E_{r} + w_{5} u_{5} = 0

\frac{\partial H}{\partial u_{6}} = {- P}_{5} I_{w} + P_{9} I_{w} + w_{6} u_{6} = 0

\frac{\partial H}{\partial u_{7}} = {- P}_{6} I_{r} + P_{10} I_{r} + w_{7} u_{7} = 0

u_{1} = {\overset{̅}{u}}_{1}

u_{2} = {\overset{̅}{u}}_{2}

u_{3} = {\overset{̅}{u}}_{3}

u_{4} = {\overset{̅}{u}}_{4}

u_{5} = {\overset{̅}{u}}_{5}

u_{6} = {\overset{̅}{u}}_{6}

, and

u_{7} = {\overset{̅}{u}}_{7}

, respectively.

By solving system (46), the controls

{\overset{̅}{u}}_{1}, {\overset{̅}{u}}_{2}, \dots, {\overset{̅}{u}}_{7}

are obtained as follows:

\frac{\partial H}{\partial u_{1}} = P_{1} λ_{w} S + P_{2} (1 - ρ) λ_{w} V - P_{3} λ_{w} S - P_{3} (1 - ρ) λ_{w} V + w_{1} u_{1} = 0

\frac{\partial H}{\partial u_{2}} = P_{1} λ_{r} S + P_{2} (1 - ρ) λ_{r} V - P_{4} λ_{r} S - P_{4} {(1 - ρ) λ}_{r} V + w_{2} u_{2} = 0

\frac{\partial H}{\partial u_{3}} = P_{5} λ_{P} I_{w} + P_{6} λ_{P} I_{r} - P_{7} λ_{P} I_{w} - P_{8} λ_{P} I_{r} + w_{3} u_{3} = 0

\frac{\partial H}{\partial u_{4}} = {- P}_{3} E_{w} + P_{13} E_{w} + w_{4} u_{4} = 0

,(47)

\frac{\partial H}{\partial u_{5}} = {- P}_{4} E_{r} + P_{13} E_{r} + w_{5} u_{5} = 0

\frac{\partial H}{\partial u_{6}} = {- P}_{5} I_{w} + P_{9} I_{w} + w_{6} u_{6} = 0

\frac{\partial H}{\partial u_{7}} = {- P}_{6} I_{r} + P_{10} I_{r} + w_{7} u_{7} = 0

Using the bounds for the controls in

U

, that is,

u_{i \min} \leq {\overset{̅}{u}}_{i} \leq u_{i \max}

for all

i = 1

,2,…,7, the optimal controls

u_{1}^{*}, u_{2}^{*}, \dots, u_{7}^{*}

are obtained as follows:

u_{1}^{*} =

min

\{\max (u_{1 \min}, \frac{β_{w} (I_{wP} + ϵ_{1} I_{w} + ϵ_{2} T_{wP} + ϵ_{3} T_{w}) \{(P_{3} - P_{1}) S + (P_{3} - P_{2}) (1 - ρ) V\}}{w_{1}}), u_{1 \max}\}

u_{2}^{*} =

min

\{\max (u_{2 \min}, \frac{β_{r} (I_{rP} + ϵ_{4} I_{r} + ϵ_{5} T_{rP} + ϵ_{6} T_{r}) \{(P_{4} - P_{1}) S + (P_{4} - P_{2}) (1 - ρ) V\}}{w_{2}}), u_{2 \max}\}

u_{3}^{*} =

min

\{\max (u_{3 \min}, \frac{β_{P} (I_{wP} + ϵ_{7} I_{rP} + ϵ_{8} T_{wP} + ϵ_{9} T_{rP}) \{(P_{7} - P_{5}) I_{w} + (P_{8} - P_{6}) I_{r}\}}{w_{3}}), u_{3 \max}\}

u_{4}^{*} =

min

\{\max (u_{4 \min}, \frac{(P_{3} - P_{13}) E_{w}}{w_{4}}), u_{4 \max}\}

,(48)

u_{5}^{*} =

min

\{\max (u_{5 \min}, \frac{(P_{4} - P_{13}) E_{r}}{w_{5}}), u_{5 \max}\}

u_{6}^{*} =

min

\{\max (u_{6 \min}, \frac{(P_{5} - P_{9}) I_{w}}{w_{6}}), u_{6 \max}\}

u_{7}^{*} =

min

\{\max (u_{7 \min}, \frac{(P_{6} - P_{10}) I_{r}}{w_{7}}), u_{7 \max}\}

The parameters

u_{i min​}

represent the suboptimal existing control strategies derived from

[33]

and are given as follows:

u_{1 \min} = 0.47

u_{2 \min} = 0.43

u_{3 \min} = 0

u_{4 \min} = 0.32

u_{5 \min} = 0

u_{6 \min} = 0.63

u_{7 \min} = 0.52

(49)

Additionally, we set

u_{i \max} = 1

for all

i = 1

,2,…,7.

It is worth noting that the optimality system consists of the state system (39) with its initial conditions, the adjoint system (45) with transversality conditions, and the control characterization (48).

5. Model Simulations

Numerical simulations of the system of model equations (1) were conducted to predict the impact of natural immunity on the reactivation of latent infections and its role in controlling or reducing the co-existence of both drug-sensitive and drug-resistant strains of pulmonary tuberculosis and pneumonia. These simulations were performed utilizing MATLAB’s built-in ordinary differential equation solver, ode45.

Additionally, the optimal strategies are determined by solving the optimality system using a fourth-order Runge-Kutta iterative scheme. The forward-backward sweep method is employed because the state variables have initial conditions, while the adjoint system has final conditions. The state variables are solved using a forward difference scheme, whereas the adjoint variables are solved using a backward difference scheme. The iterative solution scheme begins by solving the state system using the suboptimal initial values of the control variables, as provided in equation (49). The solution of the state system, along with the initial control values, is then used to solve the adjoint system. Subsequently, the control variables are updated by combining the previous controls with values obtained from the characterizations. These updated controls are then used to re-solve both the state and adjoint systems. The iterations continue until convergence is achieved. The final approximations of the control, state, and adjoint systems are considered the solution to the optimal control problem. The aim is to determine the optimal control strategies that minimize the populations co-infected with drug-sensitive TB and pneumonia, as well as those co-infected with drug-resistant TB and pneumonia.

The parameter values used are provided in Table 3, and the initial state values are as follows:

S (0 = 5, 062, 291

;

V (0) = 22, 780, 311

;

E_{w} (0) = 1, 687, 405

;

E_{r} (0) = 30364

;

I_{w} (0) = 41, 733

;

I_{r} (0) = 5000

;

I_{wP} (0) = 41, 253

;

I_{rP} =

10,167;

T_{w} (0) =

21,787;

T_{r} (0) =

3978;

T_{wP} =

27,427;

T_{rP} =

6,588; and

R (0) =

5,705,328. These values were obtained from published Kenyan data

[29-32]

. The simulations are conducted over a time span of 0 to 20 years, taking into account that pulmonary tuberculosis infections often take a long time to reactivate into chronic disease and eventually be cured. The simulation results are presented graphically in Figures 2-11.

5.1. Effects of Varying the Immunity Parameter on Co-infected Populations

Download: Download full-size image

Figure 2. Effects of varying the immunity parameter (m) on drug-resistant tuberculosis-pneumonia co-infections.

Download: Download full-size image

Figure 3. Effects of varying the immunity parameter (m) on drug-sensitive tuberculosis-pneumonia co-infections.

In the model flow chart shown in Figure 1,

m

represents the natural immunity parameter. Figures 2 and 3 illustrate the effects of varying this parameter on the co-existence of drug-resistant and drug-sensitive strains of pulmonary tuberculosis and pneumonia, respectively. An increase in the natural immunity parameter is observed to reduce the number of individuals co-infected with pulmonary tuberculosis and pneumonia. This reduction occurs because enhanced natural immunity diminishes the transmission of both tuberculosis strains by slowing the progression of latent tuberculosis infections to active pulmonary disease. Furthermore, reduced progression of latent infections results in fewer individuals requiring treatment for drug-sensitive pulmonary tuberculosis, thereby lowering the risk of developing drug-resistant strains due to improper, incomplete, or excessive use of antibiotics. Consequently, the decline in infections with both drug-resistant and drug-sensitive strains significantly reduces the number of individuals co-infected with pulmonary tuberculosis and opportunistic pneumonia.

5.2. Effects of Optimizing Prevention on Co-infected Populations

Download: Download full-size image

Figure 4. Effects of optimizing prevention on drug-resistant tuberculosis-pneumonia co-infections.

Download: Download full-size image

Figure 5. Effects of optimizing prevention on drug-sensitive tuberculosis-pneumonia co-infections. Effects of optimizing prevention on drug-sensitive tuberculosis-pneumonia co-infections.

Figures 4 and 5 depict the effects of optimizing prevention strategies on drug-resistant and drug-sensitive tuberculosis and pneumonia co-infections, respectively. It is observed that optimizing prevention strategies, while keeping other strategies at suboptimal levels, reduces drug-resistant tuberculosis and pneumonia co-infections by 18.70% and drug-sensitive tuberculosis and pneumonia co-infections by 19.10%. These reductions are attributed to the disruption of transmission chains for drug-resistant and drug-sensitive tuberculosis strains, as well as pneumonia infections, through preventive measures such as isolation, contact tracing, targeted detection of active pulmonary TB cases, social support for drug-resistant TB patients to complete treatment, monitoring individuals to ensure full treatment adherence, raising awareness about pulmonary TB and its co-morbidity with other lung infections, and patient-centered care for diagnosing co-morbidities.

5.3. Effects of Optimizing Treatment on Co-infected Populations

Download: Download full-size image

Figure 6. Effects of optimizing treatment on drug-resistant tuberculosis-pneumonia co-infections.

Download: Download full-size image

Figure 7. Effects of optimizing prevention on drug sensitive tuberculosis-pneumonia co-infections.

Figures 6 and 7 show the effects of optimizing treatment for both drug-resistant and drug-sensitive tuberculosis strains on their co-infections with pneumonia, while keeping other strategies at suboptimal levels. The results indicate that optimizing treatment alone reduces drug-resistant tuberculosis-pneumonia co-infections by 4.4%, but it has no significant effect on the population co-infected with drug-sensitive tuberculosis and pneumonia. Optimizing treatment alone is relatively ineffective because timely administration of drugs is challenging, and false negatives may occur during diagnosis, limiting its impact on drug-sensitive tuberculosis-pneumonia co-infections. However, the optimization of treatment helps manage drug-resistant tuberculosis patients more effectively, as they require longer and more complex treatment regimens. This reduces their contribution to ongoing transmission and, consequently, lowers the number of individuals co-infected with opportunistic pneumonia. Additionally, optimizing treatment improves adherence to the full drug regimen, which helps prevent the emergence of drug-resistant tuberculosis-pneumonia co-infections.

5.4. Effects of Optimizing Prevention and Treatment on Co-infected Populations

Download: Download full-size image

Figure 8. Effects of optimizing prevention and treatment on drug-resistant tuberculosis-pneumonia co-infections.

The impact of optimizing the prevention and treatment of both drug-resistant and drug-sensitive tuberculosis, coupled with the optimization of opportunistic pneumonia prevention, on populations co-infected with drug-resistant tuberculosis-pneumonia and drug-sensitive tuberculosis-pneumonia is illustrated in Figures 8 and 9, respectively. The results indicate that optimizing prevention and treatment reduces drug-resistant tuberculosis-pneumonia co-infections by 19.47% and drug-sensitive tuberculosis-pneumonia co-infections by 19.26%. Simultaneous optimization of prevention and treatment has a synergistic effect in reducing both types of co-infections. Prevention strategies interrupt transmission chains, thereby reducing new infections from both tuberculosis strains, while optimized treatment ensures adherence to full therapy and proper case management, which are critical for controlling existing drug-resistant tuberculosis cases. This reduces the infectious period and lowers the likelihood of spreading the drug-resistant strain or developing further resistance. As a result, the combined approach significantly decreases the prevalence of both drug-resistant and drug-sensitive tuberculosis-pneumonia co-infections.

Download: Download full-size image

Figure 9. Effects of optimizing prevention and treatment on drug-sensitive tuberculosis-pneumonia co-infections.

5.5. Effects of Optimizing Prevention and Screening on Co-infected Populations

Figures 10 and 11 illustrate the effects of optimizing the prevention of drug-resistant tuberculosis, drug-sensitive tuberculosis, and opportunistic pneumonia infections, coupled with the optimization of screening for both latent drug-resistant and drug-sensitive infections, on populations co-infected with drug-resistant tuberculosis-pneumonia and drug-sensitive tuberculosis-pneumonia, respectively. The results indicate that optimizing prevention and screening, while keeping the other strategies at suboptimal levels, reduces the population co-infected with drug-resistant tuberculosis and pneumonia by 25.0%, and the population co-infected with drug-sensitive tuberculosis and pneumonia by 35.09%. These reductions are attributed to the screening of latent tuberculosis infections, which allows for timely intervention before the infections progress to active disease. Additionally, preventive measures interrupt transmission chains, reducing the number of individuals exposed to both tuberculosis and pneumonia. Consequently, fewer cases of active tuberculosis results to a lower risk of immunosuppression and pulmonary complications, thereby decreasing susceptibility to opportunistic pneumonia. Screening for latent tuberculosis also lowers the progression to active drug-resistant tuberculosis, resulting in fewer patients who require prolonged treatment and, in turn, reducing community transmission and cases of co-infection with pneumonia. The reduction in drug-sensitive tuberculosis-pneumonia co-infections is greater than that of drug-resistant tuberculosis-pneumonia co-infections because latent drug-sensitive tuberculosis is easier to manage effectively with standard antibiotics.

Download: Download full-size image

Figure 10. Effects of optimizing prevention and screening on drug-resistant tuberculosis-pneumonia co-infections.

Download: Download full-size image

Figure 11. Effects of optimizing prevention and screening on drug-sensitive tuberculosis-pneumonia co-infections.

6. Conclusion

In this research study, a mathematical model of pulmonary tuberculosis in the presence of opportunistic pneumonia, incorporating the effects of a drug-resistant tuberculosis strain, was formulated. The model accounted for the role of natural immunity in controlling or reducing the transmission of opportunistic pneumonia among both drug-resistant and drug-sensitive tuberculosis patients by incorporating a saturated progression from latent infection to active pulmonary tuberculosis using Holling’s function. Additionally, the model integrated optimal control theory to identify strategies that should be optimized to minimize the co-existence of both drug-sensitive and drug-resistant tuberculosis with opportunistic pneumonia in the community.

Numerical simulations were conducted to evaluate the impact of natural immunity on the progression from latent tuberculosis infection to active disease and its role in controlling drug-resistant pulmonary tuberculosis-pneumonia co-infections. Furthermore, Pontryagin’s Maximum Principle was used to determine the optimal strategies for reducing the co-existence of pulmonary tuberculosis and opportunistic pneumonia in the presence of a drug-resistant strain. The numerical simulations revealed that enhancing immunity among latently infected individuals significantly reduces the number of co-infected cases. The optimal control analysis indicated that optimizing a combination of infection prevention and screening of latently infected individuals is the most effective strategy to control or reduce co-infections, followed by a combination of prevention and treatment. Optimizing prevention alone ranked third, while treatment alone was the least effective strategy.

Therefore, this study recommends placing greater emphasis on preventive measures against both strains of pulmonary tuberculosis and opportunistic pneumonia, coupled with the screening of latently infected individuals, to reduce the prevalence of pulmonary tuberculosis and pneumonia co-infections. It also recommends enhancing natural immunity among latently infected individuals to further reduce co-infections with both strains of pulmonary tuberculosis and opportunistic pneumonia.

The limitations of this study, which could be addressed in future research, include applying a stochastic approach to investigate the impact of randomness on the transmission dynamics of the diseases considered; utilizing a fractional-order model to enhance the accuracy of the results; and employing artificial intelligence techniques to improve the detection of infected individuals within the population.

Abbreviations

TB	Tuberculosis
S	Susceptible
E	Exposed
I	Infectious
R	Recovered
N	Total Human Population

Author Contributions

Erick Mutwiri Kirimi: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Resources, Software, Validation, Visualization, Writing - original draft, Writing - review & editing

Jeconiah Okelo: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review & editing

Mark Kimathi: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review & editing

Kenneth Ngure: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review & editing

Data Availability

The data used to support the findings of this study is included in the article.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding publication of this paper.

References

[1]	World Health Organization. (2023). Global tuberculosis report 2023. Retrieved May 10, 2025, from https://www.who.int/publications/i/item/9789240083851
[2]	World Health Organization. (2024). Tuberculosis fact sheet. Retrieved May 14, 2025, from https://www.who.int/news-room/fact-sheets/detail/tuberculosis
[3]	Garcia, R. (2019). Community‐acquired pneumonia due to Streptococcus pneumoniae: when to consider coinfection with active pulmonary tuberculosis. Case Reports in Infectious Diseases, 2019(1), 4618413. https://doi.org/10.1155/2019/4618413
[4]	Chavarría, A. G. C., & Quirós, L. B. (2018). Coinfection with Staphylococcus aureus community-acquired pneumonia and mycobacterium tuberculosis in a diabetic patient: a case report. Int Phys Med Rehab J, 3(3), 262-263.
[5]	Yan, P., & Chowell, G. (2019). Quantitative methods for investigating infectious disease outbreaks (Vol. 70). Cham, Switzerland: Springer. https://link.springer.com/book/10.1007/978-3-030-21923-9
[6]	Baguelin, M., Medley, G. F., Nightingale, E. S., O’Reilly, K. M., Rees, E. M., Waterlow, N. R., & Wagner, M. (2020). Tooling-up for infectious disease transmission modelling. Epidemics, 32, 100395. https://doi.org/10.1016/j.epidem.2020.100395
[7]	Kirimi, E. M., Muthuri, G. G., Ngari, C. G., & Karanja, S. (2024). Modeling the effects of vaccine efficacy and rate of vaccination on the transmission of pulmonary tuberculosis. Informatics in Medicine Unlocked, 46, 101470. https://doi.org/10.1016/j.imu.2024.101470
[8]	Kirimi, E. M., Muthuri, G. G., Ngari, C. G., & Karanja, S. (2024). A Model for the Propagation and Control of Pulmonary Tuberculosis Disease in Kenya. Discrete Dynamics in Nature and Society, 2024(1), 5883142. https://doi.org/10.1155/2024/5883142
[9]	Li, Q., & Wang, F. (2023). An epidemiological model for tuberculosis considering environmental transmission and reinfection. Mathematics, 11(11), 2423. https://doi.org/10.3390/math11112423
[10]	Alfiniyah, C., Soetjianto, W. S. P. A., Aziz, M. H. N., & Ghadzi, S. M. B. S. (2024). Mathematical modeling and optimal control of tuberculosis spread among smokers with case detection. AIMS Mathematics, 9(11), 30472-30492. https://www.aimspress.com/journal/Math
[11]	Xu, C., Cheng, K., Guo, S., Yuan, D., & Zhao, X. (2024). A dynamic model to study the potential TB infections and assessment of control strategies in China. arXiv preprint arXiv: 2401.12462. https://doi.org/10.48550/arXiv.2401.12462
[12]	Ochieng, F. O. (2025). SEIRS model for TB transmission dynamics incorporating the environment and optimal control. BMC Infectious Diseases, 25(1), 490. https://doi.org/10.1002/eng2.13068
[13]	Kang, T. L., Huo, H. F., & Xiang, H. (2024). Dynamics and optimal control of tuberculosis model with the combined effects of vaccination, treatment and contaminated environments. Mathematical Biosciences and Engineering, 21(4), 5308-5334.
[14]	Otoo, D., Osman, S., Poku, S. A., & Donkoh, E. K. (2021). Dynamics of Tuberculosis (TB) with Drug Resistance to First‐Line Treatment and Leaky Vaccination: A Deterministic Modelling Perspective. Computational and Mathematical Methods in Medicine, 2021(1), 5593864.
[15]	Mengistu, A. K., & Witbooi, P. J. (2023). Cost‐Effectiveness Analysis of the Optimal Control Strategies for Multidrug‐Resistant Tuberculosis Transmission in Ethiopia. International Journal of Differential Equations, 2023(1), 8822433.
[16]	Aldila, D., Awdinda, N., Herdicho, F. F., Ndii, M. Z., & Chukwu, C. W. (2023). Optimal control of pneumonia transmission model with seasonal factor: learning from Jakarta incidence data. Heliyon, 9(7).
[17]	Yano, T. K., & Bitok, J. (2022). Computational Modelling of Pneumonia Disease Transmission Dynamics with Optimal Control Analysis. Appl. Comput. Math, 11(5), 130-139. https://doi.org/10.11648/j.acm.20221105.13
[18]	Swai, M. C., Shaban, N., & Marijani, T. (2021). Optimal control in two strain pneumonia transmission dynamics. Journal of Applied Mathematics, 2021(1), 8835918.
[19]	Legesse, F. M., Rao, K. P., & Keno, T. D. (2023). Cost effectiveness and optimal control analysis for bimodal pneumonia dynamics with the effect of children's breastfeeding. Frontiers in Applied Mathematics and Statistics, 9, 1224891.
[20]	World Health Organization. (2025). Integrated approach to tuberculosis and lung health: Policy brief. Geneva: World Health Organization. Licence: CC BY-NC-SA 3.0 IGO. Retrieved May 14, 2025, from https://iris.who.int/bitstream/handle/10665/380840/9789240107526-eng.pdf%20?Sequence =1
[21]	Gweryina, R. I., Madubueze, C. E., Bajiya, V. P., & Esla, F. E. (2023). Modeling and analysis of tuberculosis and pneumonia co-infection dynamics with cost-effective strategies. Results in Control and Optimization, 10, 100210.
[22]	Xiao- Ying, L., Hong- fei, D., Li-mei, Y., & Wei-min, L. (2016). Pulmonary tuberculosis and bacterial pneumonia co-infection; a clinical etiology investigation in China. Chinese Journal of Zoonoses, 30 (4), 364-372.
[23]	Chandra, P., Grigsby, S. J., & Philips, J. A. (2022). Immune evasion and provocation by Mycobacterium tuberculosis. Nature Reviews Microbiology, 20(12), 750-766. https://www.nature.com/articles/s41579-022-00763-4
[24]	Banerjee, S. (2021). Mathematical modeling: models, analysis and applications. Chapman and Hall/CRC. https://doi.org/10.1201/9781351022941
[25]	Khajanchi, S., Das, D. K., & Kar, T. K. (2018). Dynamics of tuberculosis transmission with exogenous reinfections and endogenous reactivation. Physica A: Statistical Mechanics and its Applications, 497, 52-71. https://doi.org/10.1016/j.physa.2018.01.014
[26]	Driessche, K., Khajanchi, S., & Kar, T. K. (2020). Transmission dynamics of tuberculosis with multiple re-infections. Chaos, Solitons & Fractals, 130, 109450. https://doi.org/10.1016/j.chaos.2019.109450
[27]	Khajanchi, S., Bera, S., & Roy, T. K. (2021). Mathematical analysis of the global dynamics of a HTLV-I infection model, considering the role of cytotoxic T-lymphocytes. Mathematics and Computers in Simulation, 180, 354-378. https://doi.org/10.1016/j.matcom.2020.09.009
[28]	Chowell, G., & Kiskowski, M. (2016). Modeling ring-vaccination strategies to control Ebola virus disease epidemics. Mathematical and statistical modeling for emerging and re-emerging infectious diseases, 71-87. https://doi.org/10.1007/978-3-319-40413-4
[29]	KNBS 2023. Statistical Abstract 2023. Nairobi, Kenya. Retrieved May 21, 2025, from https://www.knbs.or.ke/wp-content/uploads/2023/11/2023-Statistical-Abstract.pdf
[30]	National Tuberculosis, Leprosy and Lung Disease Program. National strategic plan for tuberculosis, leprosy and lung health 2019-2023. Ministry of Health, Kenya. Retrieved May 23, 2025, from https://nltp.co.ke/wp-content/uploads/2020/12/National-Strategic-Plan-2019-2023.pdf
[31]	National Tuberculosis, Leprosy and Lung Disease Program. 2022 Annual Report. Ministry of Health, Kenya. Retrieved May 21, 2025, from https://nltp.co.ke/reports/annual-reports/
[32]	National Tuberculosis, Leprosy and Lung Disease Program. 2023 Annual Report. Ministry of Health, Kenya. Retrieved May 21, 2025, from https://nltp.co.ke/wp-content/uploads /2025/ 03/2023-Annual-Report.pdf
[33]	Obsu, L. L. (2022). Optimal control analysis of a tuberculosis model. Journal of Biological Systems, 30(04), 837-855.

Cite This Article

Plain Text BibTeX RIS

APA Style

Kirimi, E. M., Okelo, J., Kimathi, M., Ngure, K. (2025). Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control. American Journal of Mathematical and Computer Modelling, 10(4), 121-144. https://doi.org/10.11648/j.ajmcm.20251004.12

Copy | Download

ACS Style

Kirimi, E. M.; Okelo, J.; Kimathi, M.; Ngure, K. Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control. Am. J. Math. Comput. Model. 2025, 10(4), 121-144. doi: 10.11648/j.ajmcm.20251004.12

Copy | Download

AMA Style

Kirimi EM, Okelo J, Kimathi M, Ngure K. Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control. Am J Math Comput Model. 2025;10(4):121-144. doi: 10.11648/j.ajmcm.20251004.12

Copy | Download

@article{10.11648/j.ajmcm.20251004.12,
  author = {Erick Mutwiri Kirimi and Jeconiah Okelo and Mark Kimathi and Kenneth Ngure},
  title = {Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control
},
  journal = {American Journal of Mathematical and Computer Modelling},
  volume = {10},
  number = {4},
  pages = {121-144},
  doi = {10.11648/j.ajmcm.20251004.12},
  url = {https://doi.org/10.11648/j.ajmcm.20251004.12},
  eprint = {https://article.sciencepublishinggroup.com/pdf/10.11648.j.ajmcm.20251004.12},
  abstract = {In this paper, a deterministic mathematical model illustrating the transmission dynamics of pulmonary tuberculosis and pneumonia co-infection is formulated, incorporating a drug-resistant strain. The model employs a Holling-type functional response to capture the impact of natural immunity on the progression from latent tuberculosis infection to active disease, as well as its role in controlling drug-resistant pulmonary tuberculosis-pneumonia co-infections. The model is extended to include optimal control theory, aimed at identifying strategies to minimize co-infections using prevention, screening of latently infected individuals, and treatment as control variables. Pontryagin’s Maximum Principle is applied to characterize the optimal control system. The resulting optimality system is then solved numerically using the Runge-Kutta-based forward-backward sweep method. Numerical simulations demonstrate that enhancing natural immunity among latently infected individuals significantly reduces the number of co-infected cases. The optimal control analysis indicates that the most effective strategy for controlling or reducing co-infections of drug-resistant tuberculosis and pneumonia is the combined optimization of infection prevention and screening of latently infected individuals. These findings underscore the importance of scaling up preventive measures against pulmonary tuberculosis and opportunistic pneumonia, alongside screening efforts, to effectively control co-infections. Additionally, the study recommends strengthening immunity among latently infected populations to further reduce the prevalence of co-infections.
},
 year = {2025}
}

Copy | Download

TY  - JOUR
T1  - Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control

AU  - Erick Mutwiri Kirimi
AU  - Jeconiah Okelo
AU  - Mark Kimathi
AU  - Kenneth Ngure
Y1  - 2025/10/14
PY  - 2025
N1  - https://doi.org/10.11648/j.ajmcm.20251004.12
DO  - 10.11648/j.ajmcm.20251004.12
T2  - American Journal of Mathematical and Computer Modelling
JF  - American Journal of Mathematical and Computer Modelling
JO  - American Journal of Mathematical and Computer Modelling
SP  - 121
EP  - 144
PB  - Science Publishing Group
SN  - 2578-8280
UR  - https://doi.org/10.11648/j.ajmcm.20251004.12
AB  - In this paper, a deterministic mathematical model illustrating the transmission dynamics of pulmonary tuberculosis and pneumonia co-infection is formulated, incorporating a drug-resistant strain. The model employs a Holling-type functional response to capture the impact of natural immunity on the progression from latent tuberculosis infection to active disease, as well as its role in controlling drug-resistant pulmonary tuberculosis-pneumonia co-infections. The model is extended to include optimal control theory, aimed at identifying strategies to minimize co-infections using prevention, screening of latently infected individuals, and treatment as control variables. Pontryagin’s Maximum Principle is applied to characterize the optimal control system. The resulting optimality system is then solved numerically using the Runge-Kutta-based forward-backward sweep method. Numerical simulations demonstrate that enhancing natural immunity among latently infected individuals significantly reduces the number of co-infected cases. The optimal control analysis indicates that the most effective strategy for controlling or reducing co-infections of drug-resistant tuberculosis and pneumonia is the combined optimization of infection prevention and screening of latently infected individuals. These findings underscore the importance of scaling up preventive measures against pulmonary tuberculosis and opportunistic pneumonia, alongside screening efforts, to effectively control co-infections. Additionally, the study recommends strengthening immunity among latently infected populations to further reduce the prevalence of co-infections.

VL  - 10
IS  - 4
ER  -

Copy | Download

Author Information

Erick Mutwiri Kirimi

Department of Mathematics, Pan African University Institute for Basic Sciences, Technology, and Innovation, Nairobi, Kenya

Contact Email

http://orcid.org/0009-0009-0496-9592
Jeconiah Okelo

Department of Pure and Applied Mathematics, Jomo Kenyatta University of Agriculture and Technology, Nairobi, Kenya

Contact Email

http://orcid.org/0009-0002-9439-4110
Mark Kimathi

Department of Mathematics and Statistics, Machakos University, Machakos, Kenya

Contact Email

http://orcid.org/0000-0002-5822-8687
Kenneth Ngure

Department of Public and Community Health, Jomo Kenyatta University of Agriculture and Technology, Nairobi, Kenya

Contact Email

http://orcid.org/0000-0002-8062-0933

Download PDF

Submit an Article

Plain Text BibTeX RIS

APA Style

Kirimi, E. M., Okelo, J., Kimathi, M., Ngure, K. (2025). Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control. American Journal of Mathematical and Computer Modelling, 10(4), 121-144. https://doi.org/10.11648/j.ajmcm.20251004.12

Copy | Download

ACS Style

Kirimi, E. M.; Okelo, J.; Kimathi, M.; Ngure, K. Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control. Am. J. Math. Comput. Model. 2025, 10(4), 121-144. doi: 10.11648/j.ajmcm.20251004.12

Copy | Download

AMA Style

Kirimi EM, Okelo J, Kimathi M, Ngure K. Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control. Am J Math Comput Model. 2025;10(4):121-144. doi: 10.11648/j.ajmcm.20251004.12

Copy | Download

@article{10.11648/j.ajmcm.20251004.12,
  author = {Erick Mutwiri Kirimi and Jeconiah Okelo and Mark Kimathi and Kenneth Ngure},
  title = {Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control
},
  journal = {American Journal of Mathematical and Computer Modelling},
  volume = {10},
  number = {4},
  pages = {121-144},
  doi = {10.11648/j.ajmcm.20251004.12},
  url = {https://doi.org/10.11648/j.ajmcm.20251004.12},
  eprint = {https://article.sciencepublishinggroup.com/pdf/10.11648.j.ajmcm.20251004.12},
  abstract = {In this paper, a deterministic mathematical model illustrating the transmission dynamics of pulmonary tuberculosis and pneumonia co-infection is formulated, incorporating a drug-resistant strain. The model employs a Holling-type functional response to capture the impact of natural immunity on the progression from latent tuberculosis infection to active disease, as well as its role in controlling drug-resistant pulmonary tuberculosis-pneumonia co-infections. The model is extended to include optimal control theory, aimed at identifying strategies to minimize co-infections using prevention, screening of latently infected individuals, and treatment as control variables. Pontryagin’s Maximum Principle is applied to characterize the optimal control system. The resulting optimality system is then solved numerically using the Runge-Kutta-based forward-backward sweep method. Numerical simulations demonstrate that enhancing natural immunity among latently infected individuals significantly reduces the number of co-infected cases. The optimal control analysis indicates that the most effective strategy for controlling or reducing co-infections of drug-resistant tuberculosis and pneumonia is the combined optimization of infection prevention and screening of latently infected individuals. These findings underscore the importance of scaling up preventive measures against pulmonary tuberculosis and opportunistic pneumonia, alongside screening efforts, to effectively control co-infections. Additionally, the study recommends strengthening immunity among latently infected populations to further reduce the prevalence of co-infections.
},
 year = {2025}
}

Copy | Download

TY  - JOUR
T1  - Modeling Pulmonary Tuberculosis-Pneumonia Co-dynamics Incorporating Drug Resistance with Optimal Control

AU  - Erick Mutwiri Kirimi
AU  - Jeconiah Okelo
AU  - Mark Kimathi
AU  - Kenneth Ngure
Y1  - 2025/10/14
PY  - 2025
N1  - https://doi.org/10.11648/j.ajmcm.20251004.12
DO  - 10.11648/j.ajmcm.20251004.12
T2  - American Journal of Mathematical and Computer Modelling
JF  - American Journal of Mathematical and Computer Modelling
JO  - American Journal of Mathematical and Computer Modelling
SP  - 121
EP  - 144
PB  - Science Publishing Group
SN  - 2578-8280
UR  - https://doi.org/10.11648/j.ajmcm.20251004.12
AB  - In this paper, a deterministic mathematical model illustrating the transmission dynamics of pulmonary tuberculosis and pneumonia co-infection is formulated, incorporating a drug-resistant strain. The model employs a Holling-type functional response to capture the impact of natural immunity on the progression from latent tuberculosis infection to active disease, as well as its role in controlling drug-resistant pulmonary tuberculosis-pneumonia co-infections. The model is extended to include optimal control theory, aimed at identifying strategies to minimize co-infections using prevention, screening of latently infected individuals, and treatment as control variables. Pontryagin’s Maximum Principle is applied to characterize the optimal control system. The resulting optimality system is then solved numerically using the Runge-Kutta-based forward-backward sweep method. Numerical simulations demonstrate that enhancing natural immunity among latently infected individuals significantly reduces the number of co-infected cases. The optimal control analysis indicates that the most effective strategy for controlling or reducing co-infections of drug-resistant tuberculosis and pneumonia is the combined optimization of infection prevention and screening of latently infected individuals. These findings underscore the importance of scaling up preventive measures against pulmonary tuberculosis and opportunistic pneumonia, alongside screening efforts, to effectively control co-infections. Additionally, the study recommends strengthening immunity among latently infected populations to further reduce the prevalence of co-infections.

VL  - 10
IS  - 4
ER  -

Copy | Download