Exploring a Mathematical Model with Saturated Treatment for the Co-Dynamics of Tuberculosis and Diabetes

Abstract

In this research, we present a deterministic epidemiological mathematical model that delves into the intricate dynamics of the coexistence of tuberculosis and diabetes. Our comprehensive analysis explores the interplay and the influence of diabetes on tuberculosis incidence within a human population segregated into diabetic and non-diabetic groups. The model incorporates a saturated incidence rate and treatment regimen for latent tuberculosis infections, offering insights into their impact on tuberculosis control. The theoretical findings reveal the emergence of a phenomenon known as backward bifurcation, attributed to exogenous reinfection and saturated treatment. Additionally, our study employs both local and global sensitivity analyses to identify pivotal parameters crucial to the spread of tuberculosis within the population. This investigation contributes valuable insights to the understanding of the complex relationship between tuberculosis and diabetes, offering a foundation for more effective disease control strategies.

1. Introduction

Tuberculosis (TB) and diabetes mellitus (DM) are both global health issues and are responsible for numerous deaths annually. TB is a life-threatening contagious disease caused by Bacillus Mycobacterium tuberculosis bacteria. People can become infected with TB when they inhale pathogens released into the air by individuals with active TB who cough or sneeze. Latent tuberculosis infection (LTBI) is a subclinical disease in which the pathogenic organism causing TB remains in the body without causing symptoms, see [1,2,3].

Individuals with LTBIs do not experience symptoms, and the condition is non-infectious [1]. It is estimated that nearly 25% of the world’s population has an LTBI [3,4]. Many people with LTBIs will never develop active TB infections in their lifetime [1]. Even though individuals with LTBIs cannot transmit TB infections, they are at a 10% risk of developing active TB later [5,6]. Certain medical conditions, such as alcoholism, malnutrition, substance abuse, silicosis, renal failure, smoking, steroid therapy, indoor air pollution, and diabetes mellitus, can trigger the reactivation of LTBI. Immunocompromised individuals, such as those with HIV infection, and those with immunosuppressive conditions, such as autoimmune diseases, allergic diseases, and post-organ transplant patients, are at increased risk of developing active tuberculosis disease [7,8]. As part of the United Nations’ efforts to combat tuberculosis [9,10], it is recommended that diabetic and non-diabetic individuals who test positive for LTBI be treated.

Diabetes mellitus is a non-communicable disease that is caused by the pancreas’ inability to produce adequate levels of insulin or the body’s inability to effectively use the insulin it produces. There are a significant number of individuals globally who have diabetes, with an estimated 540 million adults living with the condition (types 1 and 2) [8,11]. Type 1 diabetes is considered an autoimmune disease, while type 2 is the most common form among adults [11,12]. People living with diabetes who are also latently infected with tuberculosis are at an increased risk of developing active tuberculosis infection [1,7,9].

The saturated incidence rate has been utilized in several disease models due to its superiority over the bilinear incidence rate [13,14,15,16,17,18,19]. This rate is considered more realistic because it takes into account the psychological and crowding effects in epidemic models. The concept of the saturated incidence rate was introduced by Capasso and Serio [15] in 1978, following their study on the transmission of cholera disease in Bari, Italy. Bentalab et al. [14] applied the saturated incidence rate in their research on the multi-strain SEIR disease model, using the force of infection as follows:

$${\displaystyle \frac{{\eta }_1{I}_1}{1+b{I}_1}}+{\displaystyle \frac{{\eta }_2{I}_2}{1+b{I}_2}},$$

where $I_1$ and $I_2$ are the infective classes, ${\eta }_1$ and ${\eta }_2$ are the transmission rates of the different strains, and b is the crowding, psychological, or inhibitory effect of the population.

There are treatments for both latent tuberculosis infections (LTBIs) and active tuberculosis. It is widely acknowledged that tuberculosis is both preventable and curable, whereas diabetes can be managed effectively. Due to the limited availability of medical resources and the large number of people with LTBIs in the population, it is not feasible to provide simultaneous treatment to all patients. As a result, saturated treatment is needed. In public health interventions, saturated treatment involves reaching and providing medical services to as many eligible patients as possible to maximize the effects of intervention strategies [4,19,20].

In the case of saturated treatment for LTBIs in both the DM and non-DM populations, health officials would first identify key groups, such as people living with HIV and their relatives, colleagues, and roommates with active TB, and test for LTBIs using the tuberculin skin test (TST) or TB blood test. Individuals with positive results are then identified and treated [4,20]. Researchers have proposed various compartmental and fractional mathematical models to better understand the transmission dynamics of infectious diseases [21,22,23,24,25,26,27,28,29].

Baba et al. [13] conducted a study on an SEI tuberculosis model that incorporates the saturated incidence rate and some control variables, which has been shown to effectively reduce the prevalence of TB in the population. Cao et al. [4] compared various LTBIs treatment strategies and examined their impact on effective TB control in China. The SEIS model, which incorporates both different treatment functions and exogenous reinfections, provides policymakers with the best approach of controlling the spread of tuberculosis (TB) disease. Egonmwan and Okuonghae [30] developed a TB mathematical model that incorporated the effect of early diagnosis, which they found to be a critical factor in reducing the TB burden in Nigeria. Sulayman et al. [6] explored the influence of imperfect vaccines and other exogenous factors in a TB model. Das et al. [31] conducted a mathematical analysis of a deterministic SEIR tuberculosis epidemiological model. Mushayabasa and Bhunu [32] modeled the consequences of early therapy for individuals with LTBI and performed an optimal control analysis. Sulayman and Abdullah [33] investigated the extent of the saturated treatment of active TB infection in the transmission mechanism of TB. Faniran et al. [34] examined the asymptotic behavior of TB in a population divided into smokers and non-smokers, with both home and public health education as control measures.

Several mathematical models that analyze the comorbidity of TB and other illnesses, including both infectious and non-infectious diseases, have been reported in the literature. For instance, Agwu et al. [35] used a deterministic mathematical model to study the co-dynamics of tuberculosis and diabetes. Similarly, Moualeu et al. [5] investigated the influence of diabetes mellitus on the transmission of TB among both diabetic and non-diabetic groups. Additionally, Ojo et al. [36] developed and examined a deterministic model for the co-dynamics of COVID-19 and TB, and demonstrated the necessary conditions for the coexistence or elimination of both diseases in the population.

Inspired by prior research on the transmission dynamics of infectious diseases [4,5,35,37], we aim to explore the effect of the DM burden on the incidence and prevalence of tuberculosis in society, the impact of saturated treatment of LTBIs among diabetic and non-diabetic populations in controlling TB disease, and the influence of exogenous reinfection on TB burden in the community. In this paper, we propose and analyze a deterministic compartmental SEIS model that naturally captures the transmission mechanism of TB among DM and non-DM individuals while incorporating the saturated treatment regimen of LTBIs in the population. To the best of our knowledge, this is the first SEIS compartmental model of the co-dynamics of tuberculosis and diabetes mellitus. The rest of the paper is structured as follows.

In Section 2, we formulate the TB-DM co-dynamics model. Section 3 describes the analysis of the proposed TB-DM model. In Section 4, we present a sensitivity analysis of the proposed model. In Section 5, we conduct uncertainty and global sensitivity analyses. In Section 6, numerical simulations are conducted. Lastly, Section 7 concludes the study.

2. Mathematical Model Formulation

Our deterministic autonomous model system is composed of nine mutually exclusive compartments where the human populations are subdivided into one of the following several groups: the state variable S describes the populations of persons that are susceptible to tuberculosis and also potentially diabetic, L describes the population of potentially diabetic individuals who have been exposed to the latent form of TB, A is the population of those non-diabetic persons infected with active tuberculosis, D describes the population of the diabetic population who are not suffering from any complication as a result of the disease, $L_d$ describes the population of those who are exposed to TB in the latent form and also have noncomplicated diabetes, $A_d$ is the population of individuals who are living with the comorbidity of active TB and diabetes in the absence of complications, those in the C class have diabetes with complications and still susceptible to TB, the people in the $L_c$ group are latently infected with TB and still suffering from diabetes with complications, and finally, members in the $A_c$ group have diabetes with complications and also living with full-blown tuberculosis disease. The total human population at time t represented by $N\left(t\right)$ is defined below as

$$N\left(t\right)=S\left(t\right)+L\left(t\right)+A\left(t\right)+D\left(t\right)+L_d\left(t\right)+A_d\left(t\right)+C\left(t\right)+L_c\left(t\right)+A_c\left(t\right).$$

2.1. Underlying Biological Assumptions of the Model

The following biological assumptions are used in the model:

• There is a random mixing of the individuals in the population and interaction between individuals in different compartments.
• We incorporate the saturated treatment of LTBIs, saturated incidence rate, endogenous reactivation, and exogenous reinfections into the nonlinear system of autonomous ordinary differential equations that govern the co-interaction dynamics of tuberculosis and diabetes.
• Tuberculosis (TB) disease progresses only from asymptomatic infection (LTBI) to the active TB infection stage.
• Individuals first develop diabetes without complications.
• The per capita natural death rate $\mu$ is used in all the compartments.
• The populations are also reduced by either TB or diabetes (with complication)-related deaths.
• The force of infection denoted by $\xi$ is defined by

  $$\xi =\omega \left({\displaystyle \frac{A}{1+bA}}+{\displaystyle \frac{{\beta }_2{A}_d}{1+b{A}_d}}+{\displaystyle \frac{{\beta }_1{A}_c}{1+b{A}_c}}\right).$$
  It is assumed that ${\beta }_1>{\beta }_2>1$ takes into consideration the high infectiousness of individuals living with diabetes with complications and active tuberculosis infection [5,35,38]. Individuals living with diabetes, even without complications, and who have acquired the TB disease potentially could indirectly increase the risk of community transmission of the Mycobacterium because of prolonged infectiousness and higher bacterial loads [5,38,39].
• People die due to complications as a result of diabetes triggered by TB disease.
• The per capita saturated treatment regimen is used for the treatment of latent TB infections at the rate ${\displaystyle \frac{q_1}{1+q_{a}L}}$ (for example, in the $L\left(t\right)$ compartment). We use the same treatment rate for diabetic individuals with latent TB infections. For the definitions of $q_1$ and $q_a$, refer to

  Table 1. Definitions of the model parameters used and their numerical values.

  Parameter	Definition	Baseline Values per Year	Range	Source
  $\Lambda$	Recruitment rate	667,685	[600,000, 700,000]	[5]
  $\mu$	Per capita natural death rate	0.02041	[0.0202, 0.02189]	[30]
  ${\sigma }_1$	Rate of acquiring diabetes mellitus free of complications by the susceptible population	0.009	[0.00466, 0.0133]	[5]
  ${\sigma }_2$	Rate of developing diabetes with complications	0.01	[0.00413, 0.0159]	[35]
  $r_1$	Per capita progression rate from LTBIs to full-blown TB disease among non-diabetic individuals	0.023	[0.00565, 0.0404]	[18]
  $r_2$	Per capita progression rate from LTBIs to full-blown TB among diabetic individuals without complications	$2r_1$	[0.0169, 0.0751]	[35]
  $r_3$	Per capita progression rate from LTBIs to full-blown TB among diabetic individuals with complications	$2r_2$	[0.0827, 0.101]	[35]
  $q_1$	Per capita treatment rate for LTBIs	1.5	[1.3, 1.8]	[19]
  ${\alpha }_1$	Transmission rate of exogenous reinfection for latently infected non-diabetic individuals	0.05	[0.0221, 0.0779]	[4]
  ${\alpha }_2$	Transmission rate of exogenous reinfection for the population with diabetes without complications but with latent TB	$1.01{\alpha }_1$	[0.00897, 0.0920]	[35]
  ${\alpha }_3$	Transmission rate of exogenous reinfection for the population of diabetic individuals with complications but with latent TB	$1.01{\alpha }_2$	[0.0349, 0.0671]	[35]
  ${\delta }_1$	The disease-induced mortality rate as a result of active TB in the population of non-diabetic individuals	0.0025	[0.00219, 0.00281]	[4]
  ${\delta }_2$	The active tuberculosis infection-induced mortality rate in the $A_d\left(t\right)$ compartment	$1.25{\delta }_1$	[0.00158, 0.00467]	[35]
  ${\delta }_3$	The active tuberculosis infection-induced mortality rate in the $A_c\left(t\right)$ compartment	$1.25{\delta }_2$	[0.000807, 0.00701]	[35]
  $d_1$	The mortality rate because of severe diabetes complications in the $L_c$ compartment	0.005	[0.00376, 0.00624]	[35]
  $d_2$	The mortality rate because of severe diabetes complications in the $A_c\left(t\right)$ compartment	0.1 $d_1$	[0.000398, 0.000602]	[35]
  $d_c$	Deaths related to diabetes complications	0.1 $d_1$	[0.000398, 0.000602]	[35]
  $\varphi$	The modification parameter for the increased rate of developing TB in latent form after being diabetic without complications	1.001	[0.459, 1.543]	Assumed
  $\phi$	The modification parameter for the increased rate of developing TB in latent form after being diabetic with complications	2.851	[2.247, 3.455]	Assumed
  ${\gamma }_1$	The modification parameter for the increased rate at which people in the $A\left(t\right)$ compartment acquired diabetes without complications	0.0001	[0.0000396, 0.000160]	Assumed
  ${\gamma }_2$	The modification parameter for the increased rate at which people in the $A_d\left(t\right)$ compartment acquired diabetes with complications	0.0001	[0.0000396, 0.000160]	Assumed
  ${\mathsf{\tau }}_1$	The rate of developing diabetes without complications among non-diabetic individuals having TB in the latent form	0.025	[0.0191, 0.0309]	Assumed
  ${\mathsf{\tau }}_2$	The rate of developing diabetes with complications among those having TB in the latent form in the $L_d$ compartment	1.01	[0.709, 1.311]	[35]
  $c_1$	The cure rate of active tuberculosis infections among non-diabetic individuals	0.546	[0.214, 0.878]	[18]
  $c_2$	The cure rate of active tuberculosis infections among diabetic individuals without complications	0.546	[0.214, 0.878]	[18]
  $c_3$	The cure rate of active tuberculosis infections among non-diabetic individuals with complications	0.546	[0.214, 0.878]	[18]
  $\omega$	Effective contact rate	$1.65\times {10}^{-8}$	[0.0000000143, 0.0000000187]	[18]
  $q_a$	Saturating factor	$6.7\times {10}^{-10}$	[0.000000000453, 0.000000000887]	[4]
  ${\beta }_1$	Modification parameter	5.597	[5.168, 6.026]	Assumed
  ${\beta }_2$	Modification parameter	5.14	[4.723, 5.557]	Assumed
  b	Saturation	0.7	[0, 1]	Assumed

  and [4,40].

The proposed SEIS model compartmentalizes humans into diabetics and nondiabetics who are susceptible to tuberculosis infections. The nondiabetic susceptible class ($S\left(t\right)$) is increased when individuals are recruited into the population either by birth or immigration at a constant rate $\Lambda$. The individuals in the susceptible compartment can develop diabetes in its uncomplicated form at the rate of ${\sigma }_1$ and also contract TB in its dormant form (latency stage) at the rate of $\xi .$ Those nondiabetic population treated for latent TB infections in the $L\left(t\right)$ compartment add to the population of the susceptible nondiabetics. People who are cured of active TB infections and are nondiabetics in the $A\left(t\right)$ compartment also increase the population of $S\left(t\right).$ The diabetic susceptible compartments $D\left(t\right)$ and $C\left(t\right)$ are also populated and depopulated similarly to the $S\left(t\right)$ compartment but we have used $\varphi$ and $\phi$ as modification parameters for the increased rates at which individuals contract latent TB infections in the $D\left(t\right)$ and $C\left(t\right)$ compartments, respectively. Diabetic individuals in the $D\left(t\right)$ compartment developed diabetes-associated complications at the rate ${\sigma }_2$. Compartments $A\left(t\right)$ and $A_d\left(t\right)$ are depopulated whenever people become diabetic without complications at the rate ${\gamma }_1{\sigma }_1$ and with complications at the rate ${\gamma }_2{\sigma }_2$, respectively. Every compartment is depleted at the rate $\mu$ when people die naturally. The population in the $L\left(t\right)$, $L_d\left(t\right)$ and $L_c\left(t\right)$ compartments are also reduced when individuals latently infected with TB develop full-blown TB disease at the rates $r_i$ and ${\alpha }_i\xi$ for $i=1,2,3$, respectively, to account for endogenous reactivation and exogenous reinfection. Individuals in the $L\left(t\right)$ and $L_d\left(t\right)$ compartments leave when they acquire diabetes at the rate of ${\mathsf{\tau }}_1$ and ${\mathsf{\tau }}_2$, respectively. The per capita saturated treatment rate of LTBIs is used in all those compartments with latently infected persons. The populations of the compartments with active TB infected persons are reduced either by TB disease-related deaths at the rate ${\delta }_1,{\delta }_2,$ and ${\delta }_3$ or when they are cured at the rates $c_1,c_2,$ and $c_3.$ The populations in the $C\left(t\right),L_c\left(t\right),$ and $A_c\left(t\right)$ compartments are also decreased when people die due to complications as a result of diabetes or comorbidity of diabetes and TB at the rate $d_c,d_1,$ and $d_2$, respectively. For detailed explanations of the parameters used in the epidemiological model and the movement of individuals between different compartments, see Table 1 and the schematic diagram of the model in Figure 1.

The proposed autonomous model is stated as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& =\Lambda -(\xi +\mu +{\sigma }_1)S+c_{1}A+{\displaystyle \frac{q_{1}L}{1+q_{a}L}}\hfill \\ \hfill {\displaystyle \frac{dL}{dt}}& =\xi S-({\mathsf{\tau }}_1+{\alpha }_1\xi +r_1+\mu )L-{\displaystyle \frac{q_{1}L}{1+q_{a}L}}\hfill \\ \hfill {\displaystyle \frac{dA}{dt}}& =({\alpha }_1\xi +r_1)L-(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)A\hfill \\ \hfill {\displaystyle \frac{dD}{dt}}& ={\sigma }_{1}S-(\varphi \xi +\mu +{\sigma }_2)D+c_2{A}_d+{\displaystyle \frac{q_1{L}_d}{1+q_a{L}_d}}\hfill \\ \hfill {\displaystyle \frac{d{L}_d}{dt}}& ={\mathsf{\tau }}_{1}L+\varphi \xi D-({\alpha }_2\xi +r_2+{\mathsf{\tau }}_2+\mu )L_d-{\displaystyle \frac{q_1{L}_d}{1+q_a{L}_d}}\hfill \\ \hfill {\displaystyle \frac{d{A}_d}{dt}}& =({\alpha }_2\xi +r_2)L_d+{\gamma }_1{\sigma }_{1}A-(\mu +{\delta }_2+c_2+{\gamma }_2{\sigma }_2)A_d\hfill \\ \hfill {\displaystyle \frac{dC}{dt}}& ={\sigma }_{2}D-(d_c+\mu +\phi \xi )C+c_3{A}_c+{\displaystyle \frac{q_1{L}_c}{1+q_a{L}_c}}\hfill \\ \hfill {\displaystyle \frac{d{L}_c}{dt}}& ={\mathsf{\tau }}_2{L}_d+\phi \xi C-(\mu +{\alpha }_3\xi +r_3+d_1)L_c-{\displaystyle \frac{q_1{L}_c}{1+q_a{L}_c}}\hfill \\ \hfill {\displaystyle \frac{d{A}_c}{dt}}& =({\alpha }_3\xi +r_3)L_c+{\gamma }_2{\sigma }_2{A}_d-(\mu +c_3+{\delta }_3+d_2)A_c\hfill \end{array}$$

(1)

where $S\left(0\right)>0,L\left(0\right)\ge 0,A\left(0\right)\ge 0,D\left(0\right)\ge 0,L_d\left(0\right)\ge 0,A_d\left(0\right)\ge 0,C\left(0\right)\ge 0,$$L_c\left(0\right)\ge 0,A_c\left(0\right)\ge 0.$

2.2. Positivity of Solutions

Theorem 1. 

The solutions of model system (1) with non-negative initial conditions, $S\left(0\right)$; $L\left(0\right)$; $A\left(0\right)$; $D\left(0\right)$; $L_d\left(0\right)$; $A_d\left(0\right)$; $C\left(0\right)$; $L_c\left(0\right)$; $A_c\left(0\right)$, will remain non-negative for all time $t>0$.

Proof. 

Let ${\tilde{t}}_f=sup\{t>0:S>0,L>0,A>0,D>0,L_d>0,A_d>0,C>0,L_c>0,A_c>0\}.$ We see that the first equation of model system (1) can be written as

$${\displaystyle \frac{dS}{dt}}+(\xi +\mu +{\sigma }_1)S\ge \Lambda ,$$

(2)

and by the integrating factor method on (2) we have that

$${\displaystyle \frac{d}{dt}}\left[S\left(t\right)exp\left({\int }_0^t\xi \left(s\right)ds+(\mu +{\sigma }_1)t\right)\right]\ge \Lambda exp\left({\int }_0^t\xi \left(s\right)ds+(\mu +{\sigma }_1)t\right).$$

Thus,

$$\begin{array}{c}\hfill S\left({\tilde{t}}_f\right)exp\left[(\mu +{\sigma }_1){\tilde{t}}_f+{\int }_0^{{\tilde{t}}_f}\xi \left(s\right)ds\right]-S\left(0\right)\ge {\int }_0^{{\tilde{t}}_f}\Lambda exp\left[(\mu +{\sigma }_1)\eta +{\int }_0^{\eta }\xi \left(s\right)ds\right]d\eta .\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill S\left({\tilde{t}}_f\right)& \ge exp\left[-(\mu +{\sigma }_1){\tilde{t}}_f-{\int }_0^{{\tilde{t}}_f}\xi \left(s\right)ds\right]S\left(0\right)\hfill \\ \hfill & +exp\left[-(\mu +{\sigma }_1){\tilde{t}}_f-{\int }_0^{{\tilde{t}}_f}\xi \left(s\right)ds\right]\left({\int }_0^{{\tilde{t}}_f}\Lambda exp\left[(\mu +{\sigma }_1)\eta +{\int }_0^{\eta }\xi \left(s\right)ds\right]d\eta \right)>0.\hfill \end{array}\end{array}$$

(3)

Hence, we see from the above inequality (3) that $S\left({\tilde{t}}_f\right)$ is non-negative. Therefore, by using a similar procedure, we can show that the remaining state variables, $L\left(t\right)$; $A\left(t\right)$; $D\left(t\right)$; $L_d\left(t\right)$; $A_d\left(t\right)$; $C\left(t\right)$; $L_c\left(t\right)$; and $A_c\left(t\right)$, are non-negative for all time $t>0.$ □

2.3. Domain Invariant

Let $\mathrm{\Omega }=\{(S\left(t\right),L\left(t\right),A\left(t\right),D\left(t\right),L_d\left(t\right),A_d\left(t\right),C\left(t\right),L_c\left(t\right),A_c\left(t\right))\in {\mathbb{R}}_{+}^9:N\left(t\right)\le {\displaystyle \frac{\Lambda }{\mu }}\}.$

Theorem 2. 

The biologically meaningful domain Ω of the TB-DM model (1) is positively invariant.

Proof. 

Since

$$\begin{array}{cc}\hfill {\displaystyle \frac{dN\left(t\right)}{dt}}& =\Lambda -\mu N\left(t\right)-({\delta }_{1}A+{\delta }_2{A}_d+{\delta }_3{A}_c+d_1{L}_c+d_2{A}_c+d_{c}C)\hfill \\ \hfill & \le \Lambda -\mu N\left(t\right)-\zeta (A+A_d+A_c+L_c+A_c+C)\hfill \\ \hfill & \le \Lambda -\mu N\left(t\right),\hfill \end{array}$$

where $\zeta =\min \{{\delta }_1,{\delta }_2,{\delta }_3,d_1,d_2,d_c\},$ we have that

$$0\le N\left(t\right)\le {\displaystyle \frac{\Lambda }{\mu }}+\left(N\left(0\right)-{\displaystyle \frac{\Lambda }{\mu }}\right)e^{-\mu t}.$$

Note that $N\left(0\right)$ is the evaluation of $N\left(t\right)$ at the initial values of each state variable. Thus, $0\le N\left(t\right)\le {\displaystyle \frac{\Lambda }{\mu }}$ as $t\to \infty .$ Thus, the domain $\mathrm{\Omega }$ is positively invariant. □

Therefore, it is sufficient to consider the dynamics of the TB-DM co-infection transmitted generated by $\left(\right)$ in $\mathrm{\Omega }$, the region where the model is mathematically and epidemiologically well-posed.

3. Equilibria and Bifurcation Analysis

The aim of this section is to find the equilibrium points and determine the type of bifurcation exhibited by the autonomous model system (1).

3.1. Tuberculosis-Free Equilibrium

The tuberculosis-free equilibrium point of the model system (1), denoted by ${\mathcal{T}}_0$, is stated as

$$\begin{array}{cc}\hfill {\mathcal{T}}_0& =\left({\mathcal{S}}_0,0,0,{\mathcal{D}}_0,0,0,{\mathcal{C}}_0,0,0\right)\hfill \end{array}$$

where ${\mathcal{S}}_0={\displaystyle \frac{\Lambda }{\mu +{\sigma }_1}},{\mathcal{D}}_0={\displaystyle \frac{{\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}$, and ${\mathcal{C}}_0={\displaystyle \frac{{\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}.$ By applying the approach in [22,41,42], we will establish the local stability of ${\mathcal{T}}_0$. From (1), it can be shown that

$${\displaystyle \frac{d}{dt}}\left(\begin{array}{c}L\left(t\right)\\ A\left(t\right)\\ L_d\left(t\right)\\ A_d\left(t\right)\\ L_c\left(t\right)\\ A_c\left(t\right)\end{array}\right)=\left(\begin{array}{c}\xi S\\ 0\\ \varphi \xi D\\ 0\\ \phi \xi C\\ 0\end{array}\right)-\left(\begin{array}{c}{\displaystyle ({\mathsf{\tau }}_1+{\alpha }_1\xi +r_1+\mu )L+\frac{q_{1}L}{1+q_{a}L}}\\ -({\alpha }_1\xi +r_1)L+(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)A\\ {\displaystyle -{\mathsf{\tau }}_{1}L+({\alpha }_2\xi +r_2+{\mathsf{\tau }}_2+\mu )L_d+\frac{q_1{L}_d}{1+q_a{L}_d}}\\ -({\alpha }_2\xi +r_2)L_d-{\gamma }_1{\sigma }_{1}A+(\mu +{\delta }_2+c_2+{\gamma }_2{\sigma }_2)A_d\\ {\displaystyle (\mu +{\alpha }_3\xi +r_3+d_1)L_c+\frac{q_1{L}_c}{1+q_a{L}_c}-{\mathsf{\tau }}_2{L}_d}\\ -({\alpha }_3\xi +r_3)L_c-{\gamma }_2{\sigma }_2{A}_d+(\mu +c_3+{\delta }_3+d_2)A_c\end{array}\right)$$

(4)

Then, let F and V be the matrices of new infections and transfers between compartments, respectively, given by

$$F=\left(\begin{array}{cccccc}0& {\displaystyle \frac{\omega \Lambda }{\mu +{\sigma }_1}}& 0& {\displaystyle \frac{\omega {\beta }_2\Lambda }{\mu +{\sigma }_1}}& 0& {\displaystyle \frac{\omega {\beta }_1\Lambda }{\mu +{\sigma }_1}}\\ 0& 0& 0& 0& 0& 0\\ 0& {\displaystyle \frac{\varphi \omega {\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}& 0& {\displaystyle \frac{\varphi \omega {\beta }_2{\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}& 0& {\displaystyle \frac{\varphi \omega {\beta }_1{\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}\\ 0& 0& 0& 0& 0& 0\\ 0& {\displaystyle \frac{\phi \omega {\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}& 0& {\displaystyle \frac{\phi \omega {\beta }_2{\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}& 0& {\displaystyle \frac{\phi \omega {\beta }_1{\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}\\ 0& 0& 0& 0& 0& 0\end{array}\right)$$

(5)

and

$$V=\left(\begin{array}{cccccc}{B}_1& 0& 0& 0& 0& 0\\ -r_1& B_2& 0& 0& 0& 0\\ -{\mathsf{\tau }}_1& 0& B_3& 0& 0& 0\\ 0& -{\gamma }_1{\sigma }_1& -r_2& B_4& 0& 0\\ 0& 0& -{\mathsf{\tau }}_2& 0& B_5& 0\\ 0& 0& 0& -{\gamma }_2{\sigma }_2& -r_3& B_6\end{array}\right)$$

(6)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill B_1& ={\mathsf{\tau }}_1+r_1+\mu +q_1,\hfill \\ \hfill B_2& =\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1,\hfill \\ \hfill B_3& =r_2+{\mathsf{\tau }}_2+\mu +q_1,\hfill \\ \hfill B_4& =\mu +{\delta }_2+c_2+{\gamma }_2{\sigma }_2,\hfill \\ \hfill B_5& =\mu +r_3+d_1+q_1,\mathrm{and}\hfill \\ \hfill B_6& =\mu +c_3+{\delta }_3+d_2.\hfill \end{array}\end{array}$$

(7)

Hence, the spectral radius of the matrix $F{V}^{-1}$ i.e., $\rho \left(F{V}^{-1}\right)$ defined the basic reproduction number, denoted by ${\mathcal{R}}_0$, which is obtained as

$${\mathcal{R}}_0=\omega ({\mathcal{R}}_{s_1}+{\mathcal{R}}_{s_2}+{\mathcal{R}}_d+{\mathcal{R}}_c)$$

(8)

where

$$\begin{array}{cc}\hfill {\mathcal{R}}_{s_1}& ={\displaystyle \frac{\Lambda }{(\mu +{\sigma }_1)}}\left({\displaystyle \frac{r_1}{B_1{B}_2}}+{\displaystyle \frac{r_1{\beta }_2{\gamma }_1{\sigma }_1}{B_1{B}_2{B}_4}}+{\displaystyle \frac{r_1{\beta }_1{\gamma }_1{\gamma }_2{\sigma }_1{\sigma }_2}{B_1{B}_2{B}_4{B}_6}}\right)\hfill \\ \hfill {\mathcal{R}}_{s_2}& ={\displaystyle \frac{\Lambda }{(\mu +{\sigma }_1)}}\left({\displaystyle \frac{{\beta }_1{r}_3{\mathsf{\tau }}_1{\mathsf{\tau }}_2}{B_1{B}_3{B}_5{B}_6}}+{\displaystyle \frac{{\beta }_2{r}_2{\mathsf{\tau }}_1}{B_1{B}_3{B}_4}}+{\displaystyle \frac{{\beta }_1{\gamma }_2{\sigma }_2{r}_2{\mathsf{\tau }}_1}{B_1{B}_3{B}_4{B}_6}}\right),\hfill \\ \hfill {\mathcal{R}}_d& ={\displaystyle \frac{\varphi {\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}\left({\displaystyle \frac{{\beta }_1{r}_3{\mathsf{\tau }}_2}{B_3{B}_5{B}_6}}+{\displaystyle \frac{r_2{\beta }_2}{B_3{B}_4}}+{\displaystyle \frac{{\beta }_1{r}_2{\gamma }_2{\sigma }_2}{B_3{B}_4{B}_6}}\right),\mathrm{and}\hfill \\ \hfill {\mathcal{R}}_c& ={\displaystyle \frac{{\sigma }_1{\sigma }_2\Lambda r_3{\beta }_1\phi }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)B_5{B}_6}}.\hfill \end{array}$$

The result that follows is established by a theorem in [42].

Theorem 3. 

The tuberculosis-free equilibrium, ${\mathcal{T}}_0$, of the model (1) is locally asymptotically stable in Ω if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

The basic reproduction number, ${\mathcal{R}}_0$, is a critical threshold in epidemiology and it is the expected average number of secondary cases of tuberculosis infections generated when an infected person with tuberculosis disease is introduced to the population subdivided into diabetic and non-diabetic groups where nobody is immune to TB disease. The biological implication of ${\mathcal{R}}_0$ is that it tells whether the tuberculosis disease is declining in the population or still spreading. In the following theorem, we investigate the impact of some model parameters on ${\mathcal{R}}_0$, the basic reproduction number.

Theorem 4. 

The partial derivative of the threshold dynamic, ${\mathcal{R}}_0$, stated in (8), with respect to the model parameter $q_1$ gives the following property:

$${\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial q_1}}<0.$$

Proof. 

We know that

$${\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial q_1}}={\displaystyle \frac{\partial {\mathcal{R}}_{s_1}}{\partial q_1}}+{\displaystyle \frac{\partial {\mathcal{R}}_{s_2}}{\partial q_1}}+{\displaystyle \frac{\partial {\mathcal{R}}_d}{\partial q_1}}+{\displaystyle \frac{\partial {\mathcal{R}}_c}{\partial q_1}},$$

and that

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathcal{R}}_{s_1}}{\partial q_1}}& =-{\displaystyle \frac{\Lambda }{(\mu +{\sigma }_1)B_1^2}}\left({\displaystyle \frac{r_1}{B_2}}+{\displaystyle \frac{r_1{\beta }_2{\gamma }_1{\sigma }_1}{B_2{B}_4}}+{\displaystyle \frac{r_1{\beta }_1{\gamma }_1{\gamma }_2{\sigma }_1{\sigma }_2}{B_2{B}_4{B}_6}}\right),\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_{s_2}}{\partial q_1}}& =-{\displaystyle \frac{\Lambda (B_1+B_3)}{(\mu +{\sigma }_1)B_1^2{B}_3^2}}\left({\displaystyle \frac{{\beta }_1{r}_3{\mathsf{\tau }}_1{\mathsf{\tau }}_2}{B_5{B}_6}}+{\displaystyle \frac{{\beta }_2{r}_2{\mathsf{\tau }}_1}{B_4}}+{\displaystyle \frac{{\beta }_1{\gamma }_2{\sigma }_2{r}_2{\mathsf{\tau }}_1}{B_4{B}_6}}+{\displaystyle \frac{{\beta }_1{r}_3{\mathsf{\tau }}_1{\mathsf{\tau }}_2{B}_1{B}_3}{B_5^2{B}_6(B_1+B_3)}}\right),\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_d}{\partial q_1}}& =-{\displaystyle \frac{\varphi {\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}\left({\displaystyle \frac{{\beta }_1{r}_3(B_3+B_5){\mathsf{\tau }}_2}{B_3^2{B}_5^2{B}_6}}+{\displaystyle \frac{r_2{\beta }_2}{B_3^2{B}_4}}+{\displaystyle \frac{{\beta }_1{r}_2{\gamma }_2{\sigma }_2}{B_3^2{B}_4{B}_6}}\right),\mathrm{and}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_c}{\partial q_1}}& =-{\displaystyle \frac{{\sigma }_1{\sigma }_2\Lambda {\beta }_1\phi r_3}{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)B_5^2{B}_6}}\hfill \end{array}$$

are clearly negative. Hence, the result follows. □

Remark 1. 

Theorem 4 implies that the saturated treatment of latent tuberculosis infections among both diabetic and non-diabetic individuals will decrease the value of ${\mathcal{R}}_0$, and will eventually ease the burden of tuberculosis disease in the population. Using the threshold quantity, ${\mathcal{R}}_0$, in (8) to investigate the outcome of saturated treatment of individuals harboring the latent forms of tuberculosis infections and healing of those with active TB diseases on the dynamics of the tuberculosis disease in a population, we see that

$$\underset{\begin{array}{c}{c}_1,c_2,c_3\to 1\\ q_1\to \infty \end{array}}{lim}{\mathcal{R}}_0=0.$$

(9)

The limit in (9) indicates that a tuberculosis control strategy that prioritizes providing large treatment of active tuberculosis infections and a high frequency of treatment to those with LTBIs in the population will inevitably result in efficient TB management because we have zero on the right-hand side of (9). Also,

$$\underset{\begin{array}{c}{c}_1,c_2,c_3\to \infty \\ q_1\to 1\end{array}}{lim}{\mathcal{R}}_0=0.$$

(10)

The limit in (10) signifies that a tuberculosis control strategy that gives much attention to providing a rapid rate of treatment of active tuberculosis infections and large treatment of latent tuberculosis infections in the population will eventually lead to effective TB management because we obtain zero on the right-hand side of (10).

3.2. Existence of Tuberculosis-Present Equilibria

We investigate the existence of the endemic equilibria for a special case of the model system (1) when the saturated treatment of the latent tuberculosis (TB) infections is assumed to be negligible (i.e., $q_1=0$). The possible endemic equilibria are those for the diabetes-free individuals and also for the coexistence of diabetics and non-diabetic populations. We neglect the endemic equilibrium for only diabetic individuals (with and without complications) because it is infeasible and impractical. Assume that ${\mathcal{T}}_e=(S^{**},L^{**},A^{**},0,0,0,0,0,0)$ is the endemic equilibrium point of the special case of model (1) examined in this section.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Lambda -(\xi +\mu +{\sigma }_1)S+c_{1}A& =0\hfill \\ \hfill \xi S-({\mathsf{\tau }}_1+{\alpha }_1\xi +r_1+\mu )L& =0\hfill \\ \hfill ({\alpha }_1\xi +r_1)L-(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)A& =0\hfill \end{array}\end{array}$$

(11)

Solving the above system of Equation (11) for $S^{**},L^{**},$ and $A^{**}$ in terms of the force of infection given by

$${\xi }^{**}=\omega {\displaystyle \frac{A^{**}}{1+b{A}^{**}}},$$

(12)

we have

$$\begin{array}{cc}\hfill S^{**}& ={\displaystyle \frac{({\mathsf{\tau }}_1+{\alpha }_1{\xi }^{**}+r_1+\mu )B_2\Lambda }{B_2({\xi }^{**}+\mu +{\sigma }_1)({\mathsf{\tau }}_1+{\alpha }_1{\xi }^{**}+r_1+\mu )-c_1{\xi }^{**}({\alpha }_1{\xi }^{**}+r_1)}},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill L^{**}& ={\displaystyle \frac{{\xi }^{**}{B}_2\Lambda }{B_2({\xi }^{**}+\mu +{\sigma }_1)({\mathsf{\tau }}_1+{\alpha }_1{\xi }^{**}+r_1+\mu )-c_1{\xi }^{**}({\alpha }_1{\xi }^{**}+r_1)}},\mathrm{and}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill A^{**}& ={\displaystyle \frac{{\xi }^{**}({\alpha }_1{\xi }^{**}+r_1)\Lambda }{B_2({\xi }^{**}+\mu +{\sigma }_1)({\mathsf{\tau }}_1+{\alpha }_1{\xi }^{**}+r_1+\mu )-c_1{\xi }^{**}({\alpha }_1{\xi }^{**}+r_1)}}.\hfill \end{array}$$

(15)

and can define

$$\begin{array}{c}\hfill {\mathcal{R}}_0^{nd}={\displaystyle \frac{\omega r_1\Lambda }{(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)({\mathsf{\tau }}_1+r_1+\mu )(\mu +{\sigma }_1)}}.\end{array}$$

By substituting $A^{**}$ into (12) and simplifying, we obtain the polynomial equation

$$M_1{\left({\xi }^{**}\right)}^2+M_2\left({\xi }^{**}\right)+M_3=0,$$

(16)

where

$$\begin{array}{cc}\hfill M_1& =(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1){\alpha }_1+(b\Lambda -c_1){\alpha }_1,\hfill \\ \hfill M_2& =(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)[({\mathsf{\tau }}_1+r_1+\mu )+{\alpha }_1(\mu +{\sigma }_1)]+b{r}_1\Lambda -(\omega {\alpha }_1\Lambda +c_1{r}_1),\mathrm{and}\hfill \\ \hfill M_3& =({\mathsf{\tau }}_1+r_1+\mu )(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)(\mu +{\sigma }_1)(1-{\mathcal{R}}_0^{nd}).\hfill \end{array}$$

If we solve (16) for ${\xi }^{**}$ and substitute it into (13), we obtain the positive equilibrium point for the nondiabetic groups. We observe that $M_3$ is always positive if and only if the critical threshold ${\mathcal{R}}_0^{nd}<1$. $M_3$ will always be negative if and only if ${\mathcal{R}}_0^{nd}$ exceeds unity. Since the polynomial Equation (16) is quadratic in ${\xi }^{**}$, its discriminant is computed as follows:

$${\left(M_2\right)}^2-4M_1{M}_3.$$

(17)

By setting (17) to zero and solving it in terms of ${\mathcal{R}}_0^{nd}$, we have that ${\mathcal{R}}_0^{nd}={\mathcal{R}}_c^{nd}$ where

$$\begin{array}{c}\hfill {\mathcal{R}}_c^{nd}=1-{\displaystyle \frac{{\left(M_2\right)}^2}{4M_1(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)(\mu +{\sigma }_1)}}.\end{array}$$

We see that the discriminant (i.e., ${\left(M_2\right)}^2-4M_1{M}_3$) is positive if and only if ${\mathcal{R}}_0^{nd}>{\mathcal{R}}_c^{nd}$, negative if and only if ${\mathcal{R}}_0^{nd}<{\mathcal{R}}_c^{nd}$, and zero if and only if ${\mathcal{R}}_0^{nd}={\mathcal{R}}_c^{nd}.$ There can be zero, one, or two tuberculosis-present (endemic) equilibrium point(s) depending on the signs of $M_1$, $M_2$ and $M_3$. We can obtain the following results for the existence of the endemic equilibrium.

Theorem 5. 

The special case of the model system (1) has

1.

a unique tuberculosis-present (endemic) equilibrium if ${\alpha }_1$ is positive, $M_1>0$, and ${\mathcal{R}}_0^{nd}$ is greater than unity;

2.

a unique tuberculosis-present (endemic) equilibrium if ${\mathcal{R}}_0^{nd}$ is unity, $M_1>0$, $M_2<0$ and ${\alpha }_1>0$;

3.

a unique tuberculosis-present (endemic) equilibrium with multiplicity of 2 if ${\mathcal{R}}_0^{nd}={\mathcal{R}}_c^{nd}$ and $M_1<0$, $M_2>0$ and ${\alpha }_1>0$;

4.

two tuberculosis-present (endemic) equilibriums if ${\mathcal{R}}_0^{nd}$ is strictly between ${\mathcal{R}}_c^{nd}$ and $M_1>0$, $M_2<0$ where ${\alpha }_1$ is positive;

5.

no tuberculosis-present (endemic) equilibrium otherwise whenever ${\mathcal{R}}_0^{nd}<{\mathcal{R}}_c^{nd}$, $M_1>0$ and $M_2<0$ or if $M_2>0$ and ${\mathcal{R}}_0^{nd}$ is less than unity where ${\alpha }_1$ is positive;

6.

a unique tuberculosis-present (endemic) equilibrium if ${\alpha }_1=0,$ $M_2>0$ and ${\mathcal{R}}_0^{nd}>1$, and zero tuberculosis-present (endemic) equilibrium if ${\mathcal{R}}_0^{nd}$ is less than or equal to unity.

Remark 2. 

Note that items (1)–(5) in Theorem 5 suggest the possibility of backward bifurcation at ${\mathcal{R}}_0^{nd}=1.$ The backward bifurcation is a phenomenon that occurs in several disease models. This is when a stable tuberculosis-free equilibrium (TBFE) point exists together with a stable endemic equilibrium point even when ${\mathcal{R}}_0^{nd}<1$. The epidemiological explanation of this phenomenon is that while ${\mathcal{R}}_0<1$ is necessary for disease elimination, it is not sufficient for the effective control of the tuberculosis disease. Hence, there must be a new critical threshold to wipe out the disease [43,44,45].

3.3. Existence of Backward Bifurcation

In this part, the backward bifurcation analysis of the system (1) will be performed. Let

$$E{E}^{**}=(S^{**},L^{**},A^{**},D^{**},L_d^{**},A_d^{**},C^{**},L_c^{**},A_c^{**})$$

be an arbitrary tuberculosis-present (endemic) equilibrium point. For the purpose of analyzing the bifurcation of the system (1) at ${\mathcal{R}}_0=1$, we suppose

$$S=y_1,L=y_2,A=y_3,D=y_4,L_d=y_5,A_d=y_6,C=y_7,L_c=y_8,A_c=y_9.$$

Then,

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{y}_1}{dt}}& =\Lambda -(\xi +\mu +{\sigma }_1)y_1+c_1{y}_3+{\displaystyle \frac{q_1{y}_2}{1+q_a{y}_2}}:=f_1\hfill \\ \hfill {\displaystyle \frac{d{y}_2}{dt}}& =\xi y_1-({\mathsf{\tau }}_1+{\alpha }_1\xi +r_1+\mu )y_2-{\displaystyle \frac{q_1{y}_2}{1+q_a{y}_2}}:=f_2\hfill \\ \hfill {\displaystyle \frac{d{y}_3}{dt}}& =({\alpha }_1\xi +r_1)y_2-(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)y_3:=f_3\hfill \\ \hfill {\displaystyle \frac{d{y}_4}{dt}}& ={\sigma }_1{y}_1-(\varphi \xi +\mu +{\sigma }_2)y_4+c_2{y}_6+{\displaystyle \frac{q_1{y}_5}{1+q_a{y}_5}}:=f_4\hfill \\ \hfill {\displaystyle \frac{d{y}_5}{dt}}& ={\mathsf{\tau }}_1{y}_2+\varphi \xi y_4-({\alpha }_2\xi +r_2+{\mathsf{\tau }}_2+\mu )y_5-{\displaystyle \frac{q_1{y}_5}{1+q_a{y}_5}}:=f_5\hfill \\ \hfill {\displaystyle \frac{d{y}_6}{dt}}& =({\alpha }_2\xi +r_2)y_5+{\gamma }_1{\sigma }_1{y}_3-(\mu +{\delta }_2+c_2+{\gamma }_2{\sigma }_2)y_6:=f_6\hfill \\ \hfill {\displaystyle \frac{d{y}_7}{dt}}& ={\sigma }_2{y}_4-(d_c+\mu +\phi \xi )y_7+c_3{y}_9+{\displaystyle \frac{q_1{y}_8}{1+q_a{y}_8}}:=f_7\hfill \\ \hfill {\displaystyle \frac{d{y}_8}{dt}}& =\phi \xi y_7+{\mathsf{\tau }}_2{y}_5-(\mu +{\alpha }_3\xi +r_3+d_1)y_8-{\displaystyle \frac{q_1{y}_8}{1+q_a{y}_8}}:=f_8\hfill \\ \hfill {\displaystyle \frac{d{y}_9}{dt}}& =({\alpha }_3\xi +r_3)y_8+{\gamma }_2{\sigma }_2{y}_6-(\mu +c_3+{\delta }_3+d_2)y_9:=f_9\hfill \end{array}$$

(18)

where

$$\xi =\omega \left({\displaystyle \frac{y_3}{1+b{y}_3}}+{\displaystyle \frac{{\beta }_2{y}_6}{1+b{y}_6}}+{\displaystyle \frac{{\beta }_1{y}_9}{1+b{y}_9}}\right).$$

Suppose that ${\mathcal{R}}_0=\omega \psi$ where $\psi ={\mathcal{R}}_{s_1}+{\mathcal{R}}_{s_2}+{\mathcal{R}}_d+{\mathcal{R}}_c$. Assume that $\omega$ is the bifurcation parameter. Then, when ${\mathcal{R}}_0=1$, we have that

$$\omega ={\omega }^{*}={\displaystyle \frac{1}{\psi }}.$$

Now, we evaluate the Jacobian matrix of the system (18) at ${\mathcal{T}}_0$ to have

$$J=\left(\begin{array}{ccccccccc}-(\mu +{\sigma }_1)& q_1& -\omega y_1^{*}+c_1& 0& 0& -\omega {\beta }_2{y}_1^{*}& 0& 0& -\omega {\beta }_1{y}_1^{*}\\ 0& -B_1& \omega y_1^{*}& 0& 0& \omega {\beta }_2{y}_1^{*}& 0& 0& \omega {\beta }_1{y}_1^{*}\\ 0& r_1& -B_2& 0& 0& 0& 0& 0& 0\\ {\sigma }_1& 0& -\varphi \omega y_4^{*}& -(\mu +{\sigma }_2)& q_1& -\varphi \omega {\beta }_2{y}_4^{*}+c_2& 0& 0& -\varphi \omega {\beta }_1{y}_4^{*}\\ 0& {\mathsf{\tau }}_1& \varphi \omega y_4^{*}& 0& -B_3& \varphi \omega {\beta }_2{y}_4^{*}& 0& 0& \varphi \omega {\beta }_1{y}_4^{*}\\ 0& 0& {\gamma }_1{\sigma }_1& 0& r_2& -B_4& 0& 0& 0\\ 0& 0& -\phi \omega y_7^{*}& {\sigma }_2& 0& -\phi \omega {\beta }_2{y}_7^{*}& -(d_c+\mu )& q_1& c_3-\phi \omega {\beta }_1{y}_7^{*}\\ 0& 0& \phi \omega y_7^{*}& 0& {\mathsf{\tau }}_2& \phi \omega {\beta }_2{y}_7^{*}& 0& -B_5& \phi \omega {\beta }_1{y}_7^{*}\\ 0& 0& 0& 0& 0& {\gamma }_2{\sigma }_2& 0& r_3& -B_6\end{array}\right),$$

(19)

where $B_i$ for $i=1,2,\cdots ,6$ are as defined in (7). Let ${\Delta }_1=u_3+{\beta }_2{u}_6+{\beta }_1{u}_9$, and ${\Delta }_2=u_3^2+{\beta }_2{u}_6^2+{\beta }_1{u}_9^2.$ Now, we use the approach based on the center manifold developed in [46] to determine the direction of the bifurcation at ${\mathcal{R}}_0=1$. Note that zero is a simple eigenvalue of the matrix J. Let $u={[u_1,u_2,u_3,u_4,u_5,u_6,u_7,u_8,u_9]}^{\top }$ be the right eigenvector that corresponds to the zero eigenvalue of J. So,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{u}_1\hfill & ={\displaystyle \frac{(q_1{B}_2+r_1(c_1-\omega y_1^{*}))u_3-\omega {\beta }_2{r}_1{y}_1^{*}{u}_6-\omega {\beta }_1{r}_1{y}_1^{*}{u}_9}{(\mu +{\sigma }_1)r_1}},\hfill \\ u_2\hfill & ={\displaystyle \frac{\omega {\Delta }_1}{B_1}}{y}_1^{*},\hfill \\ u_3\hfill & =u_3>0,\hfill \\ u_4\hfill & ={\displaystyle \frac{{\sigma }_1{u}_1+q_1{u}_5+\left(c_2-\varphi \omega {\beta }_2{y}_4^{*}\right)u_6-\varphi \omega y_4^{*}{u}_3-\varphi \omega {\beta }_1{y}_4^{*}{u}_9}{\mu +{\sigma }_2}},\hfill \\ u_5\hfill & ={\displaystyle \frac{B_4{u}_6-{\gamma }_1{\sigma }_1{u}_3}{r_2}},\hfill \\ u_6\hfill & =u_6>0,\hfill \\ u_7\hfill & ={\displaystyle \frac{{\sigma }_2{u}_4+\left(c_3-\phi \omega {\beta }_1{y}_7^{*}\right)u_9+q_1{u}_8-\phi \omega y_7^{*}{u}_3-\phi \omega {\beta }_2{y}_7^{*}{u}_6}{(d_c+\mu )}},\hfill \\ u_8\hfill & ={\displaystyle \frac{B_6{u}_9-{\gamma }_2{\sigma }_2{u}_6}{r_3}},\mathrm{and}\hfill \\ u_9\hfill & =u_9>0.\hfill \end{array}\right.\end{array}$$

The left eigenvector represented by v that corresponds to the zero eigenvalue of J is $v={[v_1,v_2,v_3,v_4,v_5,v_6,v_7,v_8,v_9]}^{\top }$ where,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & v_1=v_4=v_7=0,\hfill \\ \hfill & v_2=v_2>0,\hfill \\ \hfill & v_3={\displaystyle \frac{B_1{v}_2-{\mathsf{\tau }}_1{v}_5}{r_1}},\hfill \\ \hfill & v_5={\displaystyle \frac{r_2{B}_5{v}_6+{\mathsf{\tau }}_2{r}_3{v}_9}{B_3{B}_5}},\hfill \\ \hfill & v_6={\displaystyle \frac{\omega {\beta }_2{y}_1^{*}{v}_2+\varphi \omega {\beta }_2{y}_4^{*}{v}_5+\phi \omega {\beta }_2{y}_7^{*}{v}_8+{\gamma }_2{\sigma }_2{v}_9}{B_4}},\hfill \\ \hfill & v_8={\displaystyle \frac{r_3}{B_5}}{v}_9,\mathrm{and}\hfill \\ \hfill & v_9=v_9>0.\hfill \end{array}\right.\end{array}$$

Then, we compute the bifurcation coefficients given by

$$\begin{array}{cc}\hfill a& =\sum _{k,i,j=1}^n{v}_k{u}_i{u}_j{\displaystyle \frac{{\partial }^2{f}_k}{\partial y_i\partial y_j}}(0,0)\mathrm{and}\hfill \\ \hfill b& =\sum _{k,i=1}^n{v}_k{u}_i{\displaystyle \frac{{\partial }^2{f}_k}{\partial y_i\partial {\omega }^{*}}}(0,0)\hfill \end{array}$$

Hence,

$$b=(u_3+{\beta }_2{u}_6+{\beta }_1{u}_9)\left(v_2{y}_1^{*}+\varphi v_5{y}_4^{*}+\phi v_8{y}_7^{*}\right)>0$$

and

$$\begin{array}{cc}\hfill a& =v_2\omega \left((u_1-{\alpha }_1{u}_2){\Delta }_1-2b{\Delta }_2{y}_1^{*}\right)+v_5\varphi \omega \left((u_4-{\alpha }_2{u}_5){\Delta }_1-2b{\Delta }_2{y}_4^{*})\right)\\ & +({\alpha }_1{u}_2+{\alpha }_2{u}_5+{\alpha }_3{u}_8)\omega {\Delta }_1+v_8\phi \omega \left((u_7-{\alpha }_3{u}_8){\Delta }_1-2b{\Delta }_2{y}_7^{*})\right).\end{array}$$

We draw the following conclusion to classify the type of bifurcation at ${\mathcal{R}}_0=1$ based on Theorem 4.1 in Castillo-Chavez and Song [46].

Theorem 6. 

If the coefficient a is positive, then the tuberculosis-diabetes co-dynamics governed by the model given in (1) will undergo a backward bifurcation at ${\mathcal{R}}_0=1.$ Conversely, if a is negative, then the model system (1) will undergo a forward bifurcation at ${\mathcal{R}}_0=1.$

In Figure 2, we have used the following parameters: $\mu =0.020$, ${\delta }_1=0.025$, $c_1=0.005$, ${\gamma }_1=1.01$, ${\sigma }_1=0.0009$, ${\alpha }_1=4.28$, $b=5\times {10}^{-8}$, $\Lambda =667685$, ${\mathsf{\tau }}_1=0.25$, $r_1=0.023$, and $\omega =1.65\times {10}^{-8}$. We choose ${\alpha }_1$ such that ${\mathcal{R}}_0^{nd}=0.8128$ and ${\mathcal{R}}_c^{nd}=0.4942.$ The backward bifurcation is induced by the exogenous reinfection ${\alpha }_1$ when ${\mathcal{R}}_c^{nd}<{\mathcal{R}}_0^{nd}<1.$

3.4. Global Asymptotic Stability of the Tuberculosis-Free Equilibrium Point

In this study, we will employ the method outlined in [47] to investigate the global asymptotic stability of the tuberculosis-free equilibrium point (TBFE) for the model system (1). To confirm the global asymptotic stability of the tuberculosis-free equilibrium point, it is essential to fulfill two specific conditions. One, we rewrite the system (1) in the form

$$\begin{array}{cc}\hfill {\displaystyle \frac{dX}{dt}}& =G(X,I),\hfill \\ \hfill {\displaystyle \frac{dI}{dt}}& =H(X,I),H(X,0)=0\hfill \end{array}$$

(20)

where $X\in {\mathbb{R}}^m$ represents (its components) the number of persons free of infections and $I\in {\mathbb{R}}^n$ represents (its components) the number of infected persons, which includes latent, infectious, etc. $U_0=(X^{*},0)$ is the tuberculosis-free equilibrium point of this system.

The subsequent conditions $P1$ and $P2$ below must be satisfied to guarantee the global asymptotic stability of the TBFE.

(P1)

For ${\displaystyle \frac{dX}{dt}=G(X,0)}$, $X^{*}$ is globally asymptotically stable (GAS).

(P2)

$H(X,I)=BI-\widehat{H}(X,I)$, $\widehat{H}(X,I)\ge 0$ for $(X,I)\in \mathrm{\Omega }$,

where

$$B=D_{I}H(X^{*},0)$$

is an M-matrix, i.e., the off diagonal elements of B are non-negative and $\mathrm{\Omega }$ is the region where the model is biologically meaningful.

If the system (20) satisfies the conditions $\left(P1\right)$ and $\left(P2\right)$ above, then the theorem below holds.

Theorem 7. 

The fixed point $U_0=(X^{*},0)$ of the model system (20) is globally asymptotic stable (GAS) provided that ${\mathcal{R}}_0<1$ (l.a.s) and that assumptions $\left(P1\right)$ and $\left(P2\right)$ are met.

Proof. 

Since ${\mathcal{R}}_0<1$, we will now verify the conditions $\left(P1\right)$ and $\left(P2\right)$. For the model system (1), we have that

$$\begin{array}{c}\hfill {\displaystyle \frac{dX}{dt}}=G(X,I)=\left(\begin{array}{c}{\displaystyle \Lambda -(\xi +\mu +{\sigma }_1)S+c_{1}A+\frac{q_{1}L}{1+q_{a}L}}\\ {\displaystyle {\sigma }_{1}S-(\varphi \xi +\mu +{\sigma }_2)D+c_2{A}_d+\frac{q_1{L}_d}{1+q_a{L}_d}}\\ {\displaystyle {\sigma }_{2}D-(d_c+\mu +\phi \xi )C+c_3{A}_c+\frac{q_1{L}_c}{1+q_a{L}_c}}\end{array}\right)\end{array}$$

(21)

and

$$\begin{array}{c}\hfill {\displaystyle \frac{dX}{dt}}=G(X,0)=\left(\begin{array}{c}\Lambda -(\mu +{\sigma }_1)S\\ {\sigma }_{1}S-(\mu +{\sigma }_2)D\\ {\sigma }_{2}D-\mu C\end{array}\right).\end{array}$$

(22)

This implies that

$${\displaystyle X^{*}=\left(\frac{\Lambda }{\mu +{\sigma }_1},0,0,\frac{{\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)},0,0,\frac{{\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)},0,0\right)}$$

is globally asymptotic stable.

Let ${\displaystyle S^{*}=\frac{\Lambda }{\mu +{\sigma }_1}}$, ${\displaystyle D^{*}=\frac{{\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)},}$ and ${\displaystyle C^{*}=\frac{{\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}.}$ We also have that

$$\begin{array}{c}\hfill {\displaystyle \frac{dZ}{dt}}=H(X,I)=\left(\begin{array}{c}{\displaystyle \xi S-({\mathsf{\tau }}_1+{\alpha }_1\xi +r_1+\mu )L-\frac{q_{1}L}{1+q_{a}L}}\\ ({\alpha }_1\xi +r_1)L-(\mu +{\delta }_1+c_1+{\gamma }_1{\sigma }_1)A\\ {\displaystyle {\mathsf{\tau }}_{1}L+\varphi \xi D-({\alpha }_2\xi +r_2+{\mathsf{\tau }}_2+\mu )L_d-\frac{q_1{L}_d}{1+q_a{L}_d}}\\ ({\alpha }_2\xi +r_2)L_d+{\gamma }_1{\sigma }_{1}A-(\mu +{\delta }_2+c_2+{\gamma }_2{\sigma }_2)A_d\\ {\displaystyle \phi \xi C+{\mathsf{\tau }}_2{L}_d-(\mu +{\alpha }_3\xi +r_3+d_1)L_c-\frac{q_1{L}_c}{1+q_a{L}_c}}\\ ({\alpha }_3\xi +r_3)L_c+{\gamma }_2{\sigma }_2{A}_d-(\mu +c_3+{\delta }_3+d_2)A_c\end{array}\right)\end{array}$$

(23)

and

$$B=\left(\begin{array}{cccccc}-B_1& {\displaystyle \frac{\omega \Lambda }{\mu +{\sigma }_1}}& 0& {\displaystyle \frac{\omega {\beta }_2\Lambda }{\mu +{\sigma }_1}}& 0& {\displaystyle \frac{\omega {\beta }_1\Lambda }{\mu +{\sigma }_1}}\\ r_1& -B_2& 0& 0& 0& 0\\ {\mathsf{\tau }}_1& {\displaystyle \frac{\varphi \omega {\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}& -B_3& {\displaystyle \frac{\varphi \omega {\beta }_2{\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}& 0& {\displaystyle \frac{\varphi \omega {\beta }_1{\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}\\ 0& {\gamma }_1{\sigma }_1& r_2& -B_4& 0& 0\\ 0& {\displaystyle \frac{\phi \omega {\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}& {\mathsf{\tau }}_2& {\displaystyle \frac{\phi \omega {\beta }_2{\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}& -B_5& {\displaystyle \frac{\phi \omega {\beta }_1{\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}}\\ 0& 0& 0& {\gamma }_2{\sigma }_2& r_3& -B_6\end{array}\right).$$

Therefore,

$$\begin{array}{cc}\hfill BI& =\left(\begin{array}{c}{\displaystyle -B_{1}L+\frac{\omega \Lambda }{\mu +{\sigma }_1}A+\frac{\omega {\beta }_2\Lambda }{\mu +{\sigma }_1}{A}_d+\frac{\omega {\beta }_1\Lambda }{\mu +{\sigma }_1}{A}_c}\\ r_{1}L-B_{2}A\\ {\displaystyle {\mathsf{\tau }}_{1}L+\frac{\varphi \omega {\sigma }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}A-B_3{L}_d+\frac{\varphi \omega {\sigma }_1{\beta }_2\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}{A}_d+\frac{\varphi \omega {\sigma }_1{\beta }_1\Lambda }{(\mu +{\sigma }_1)(\mu +{\sigma }_2)}{A}_c}\\ {\gamma }_1{\sigma }_{1}A+r_2{L}_d-B_4{A}_d\\ {\displaystyle \frac{\phi \omega {\sigma }_1{\sigma }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}A+{\mathsf{\tau }}_2{L}_d+\frac{\phi \omega {\sigma }_1{\sigma }_2{\beta }_2\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}{A}_d-B_5{L}_c+\frac{\phi \omega {\sigma }_1{\sigma }_2{\beta }_1\Lambda }{(\mu +d_c)(\mu +{\sigma }_1)(\mu +{\sigma }_2)}{A}_c}\\ {\gamma }_2{\sigma }_2{A}_d+r_3{L}_c-B_6{A}_c\end{array}\right).\hfill \end{array}$$

Then,

$$H(X,I)=BI-\widehat{H}(X,I)$$

where

$$\widehat{H}=\left(\begin{array}{c}{\displaystyle {\alpha }_1\xi L-\frac{q_1{q}_a{L}^2}{1+q_{a}L}+\omega A\left(\frac{(S^{*}-S)+bA{S}^{*}}{1+bA}\right)+\omega {\beta }_2{A}_d\left(\frac{(S^{*}-S)+b{A}_d{S}^{*}}{1+b{A}_d}\right)+\omega {\beta }_1{A}_c\left(\frac{(S^{*}-S)+b{A}_c{S}^{*}}{1+b{A}_c}\right)}\\ -{\alpha }_1\xi L\\ {\displaystyle -\frac{q_1{q}_a{L}_d^2}{1+q_{a}Ld}+{\alpha }_2\xi L_d+\varphi \omega A\left(\frac{(D^{*}-D)+bA{D}^{*}}{1+bA}\right)+\varphi \omega {\beta }_2\left(\frac{(D^{*}-D)+b{A}_d{D}^{*}}{1+b{A}_d}\right)+\varphi \omega {\beta }_1\left(\frac{(D^{*}-D)+b{A}_c{D}^{*}}{1+b{A}_c}\right)}\\ -{\alpha }_2\xi L_d\\ {\displaystyle -\frac{q_1{q}_a{L}_c^2}{1+q_a{L}_c}+\phi \omega A\left(\frac{(C^{*}-C)+bA{C}^{*}}{1+bA}\right)+\phi \omega {\beta }_2\left(\frac{(C^{*}-C)+b{A}_d{C}^{*}}{1+b{A}_d}\right)+\phi \omega {\beta }_1\left(\frac{(C^{*}-C)+b{A}_c{C}^{*}}{1+b{A}_c}\right)}\\ -{\alpha }_3\xi L_c\end{array}\right).$$

(24)

□

We see that the condition $\left(P2\right)$ is not satisfied because $\widehat{H}(X,I)$ is not positive semi-definite. As a consequence, the disease-free equilibrium $U_0$ may not be globally asymptotically stable. If we let $q_1={\alpha }_1={\alpha }_2={\alpha }_3=0$ in (24), then we have that $\widehat{H}(X,I)\ge 0.$ Thus, the treatment coefficient of latent tuberculosis infections and the exogenous re-infection rate cause the phenomenon known as the backward bifurcation.

4. Sensitivity Analysis

The primary goal of sensitivity analysis is to identify the critical model parameters and their significance in the occurrence and prevalence of tuberculosis. The normalized forward sensitivity index of a variable in relation to a parameter is considered as the proportion of the relative change in the variable ${\mathcal{R}}_0$ to the relative change in the parameter [5,34,48]. In cases where the parameters are known with minimal uncertainty, we can calculate the partial derivative of the variable with respect to the parameters.

Given that ${\mathcal{R}}_0$ is a differentiable function of the parameters and serves as a crucial threshold for understanding the disease’s future progression, we conducted a sensitivity analysis to determine the model parameters that have a relative impact on the value of ${\mathcal{R}}_0$. As an illustration, we compute the sensitivity index of ${\mathcal{R}}_0$ with respect to $\omega$ below:

$$\left({\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial \omega }}\right)\times \left({\displaystyle \frac{\omega }{{\mathcal{R}}_0}}\right)=1>0.$$

It should be noted that we have used the values of the parameters given in Table 1 in the above computation.

Table 2. Sensitivity index of parameters.

Parameter	Index	Parameter	Index
$r_1$	$0.0685$	$r_2$	$0.1304$
$r_3$	$0.7489$	$c_1$	$-0.0670$
$c_2$	$-0.1306$	$c_3$	$-0.7595$
${\delta }_1$	$-3.0664\times {10}^{-4}$	${\delta }_2$	$-7.4721\times {10}^{-4}$
${\delta }_3$	$-0.0054$	${\sigma }_1$	$0.6078$
${\sigma }_2$	$0.3088$	$\omega$	1
${\beta }_1$	0.7940	${\beta }_2$	0.1362
${\mathsf{\tau }}_1$	0.0150	${\mathsf{\tau }}_2$	0.0590
$\varphi$	0.3046	$\phi$	0.6093
${\gamma }_1$	$4.5639\times {10}^{-7}$	${\gamma }_2$	$2.0675\times {10}^{-8}$
$q_1$	$-1.0057$	$\mu$	$-1.9515$
$d_1$	$-0.0025$	$d_2$	$-6.9553\times {10}^{-4}$
$d_c$	$-0.0146$	$\Lambda$	1

shows the sensitivity indices of the remaining parameters.

From Table 2, parameters with negative sensitivity indices will lower the value of ${\mathcal{R}}_0$, while those with positive sensitivity indices have strong influences on the transmission of tuberculosis.

5. Uncertainty and Global Sensitivity Analysis

Given that those single values assigned as estimates of the parameters used in mathematical models are just approximations that come with uncertainties, and also due to the unreliability of the local sensitivity analysis when it comes to assessing uncertainty and sensitivity precisely in any biological system [48,49], there is a need to perform the global uncertainty and sensitivity analyses to identify those influential parameters on the dynamics of tuberculosis in the TB-DM model (1).

The model of the comorbid interaction of tuberculosis and diabetes (1) contains 31 parameters whose baseline values given in Table 1 are expected to have uncertainties either due to measurement errors, variability, random sampling error, etc. Following the same methodology in [49], we employ a stratified sampling technique known as Latin hypercube sampling (LHS) to measure the degree of confidence in the estimation of the model parameters, while the partial rank correlation coefficient (PRCC) is used to perform the global sensitivity analysis and ${\mathcal{R}}_0$, the basic reproduction number, is used as response function. The analysis is performed with the baseline values and range of the parameters stated in Table 1. A total of 1000 simulations were executed with R software (version Rstudio2023.12.0+369).

Figure 3 shows that the parameters with the greatest impact on the transmission dynamics of tuberculosis in a population divided into diabetic and non-diabetic populations are the rate of developing diabetes without complications by the susceptible individuals (${\sigma }_1)$, the rate of developing diabetes with complications (${\sigma }_2)$, the treatment rate of tuberculosis among diabetes with complication patients ($c_3$), the saturated treatment rate of LTBIs ($q_1$), the effective contact rate ($\omega$), and the modification parameter that takes into consideration the high rate of contracting TB in its latent form in those sick with diabetes with or without complications ($\varphi ,\phi$). With the basic reproduction number (${\mathcal{R}}_0$) as response function, the PRCC-ranked parameter values of the TB-DM model (1) are listed in

Table 3. PRCC values for the parameters of the TB-DM Model (1) using the ${\mathcal{R}}_0$ as response function. the top influential model parameters are highlighted in bold font.

Parameter	${\mathcal{R}}_0$	Parameter	${\mathcal{R}}_0$	Parameter	${\mathcal{R}}_0$
${\beta }_1$	0.28317105	${\beta }_2$	0.06460933	$r_3$	0.31721437
$\mu$	−0.35729983	${\gamma }_1$	0.04053560	$d_1$	0.01963876
${\gamma }_2$	−0.01760227	${\delta }_1$	0.03087338	$d_2$	−0.07371650
${\delta }_2$	−0.01503126	${\delta }_3$	−0.04763940	$\varphi$	0.64292017
${\mathsf{\tau }}_1$	0.02483570	${\sigma }_1$	0.82610899	$\phi$	0.51019187
${\sigma }_2$	0.68273337	$c_1$	−0.25271347	$c_2$	−0.44545779
$r_1$	0.34626172	$r_2$	0.39143764	$c_3$	−0.92444959
$q_1$	−0.58589007	$\omega$	0.51280401	$\Lambda$	0.36138587
${\mathsf{\tau }}_2$	0.09809109	$d_c$	−0.05637532

.

6. Numerical Simulation

In this segment, we validate some theoretical findings by implementing the numerical simulations of the model system (1). We solve the model system (1) numerically with ode45 routine of MATLAB R2022a [50] and the simulations were performed using various biological parameters given in Table 1. We choose, arbitrarily, the following initial conditions: $S\left(0\right)=49,336,500$, $L\left(0\right)=4,148,950$, $A\left(0\right)=32,180$, $D\left(0\right)=1,128,800$, $L_d\left(0\right)=12,300$, $A_d\left(0\right)=10,100$, $C\left(0\right)=54,400$, $L_c\left(0\right)=27,700$ and $A_c\left(0\right)=1050$.

In Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, the simulation results show that the absence of the saturated treatment of LTBIs (setting $q_1=0$) will lead to the persistence of the tuberculosis disease in all the compartments with actively and latently infected individuals. The basic reproduction number ${\mathcal{R}}_0$ is $3.2653$ when $q_1=0$ and this established our theoretical finding about the endemicity of the disease when ${\mathcal{R}}_0>1$. With the saturated treatment of LTBIs, the basic reproduction number ${\mathcal{R}}_0$ is $0.1384$. It is observed that tuberculosis disease dies out eventually in both the diabetic and non-diabetic populations within 20 years.

In a case when $q_a=0$, we have the unsaturated treatment regimen of latent tuberculosis infections (LTBIs), i.e., the treatment is now assumed to be proportional to the number of individuals exposed to tuberculosis infections. Also, by setting $b=0$, we have the bilinear incidence rate and the per capita treatment rate of latent tuberculosis infections $q_1$ is decreased by about 87% to perform the simulations in Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15, i.e., $q_1=0.2$. We compare model (1) with saturated incidence rate and treatment of LTBIs with the unsaturated treatment intervention of LTBIs and bilinear incidence rate. We observe that tuberculosis infections in both latent and active forms can be eliminated from the population very fast when individuals exposed to TB are treated with the saturated treatment regimen and saturated incidence rate utilized. Both active TB infections and latent TB infections stay around for a long with the bilinear incidence rate and unsaturated treatment regimen but eventually die out in the long run. Due to the concerning rise in the incidence of diabetes over the past decades [12], we quadruple the parameter values ${\sigma }_1,{\sigma }_2,{\mathsf{\tau }}_1,$ and ${\mathsf{\tau }}_2$ in the simulation for the bilinear incidence rate and saturated treatment with the intention of investigating the influence of the increased incidence of diabetes on the spread of TB infections. The numerical simulations in Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15 confirmed the persistence of tuberculosis disease in the population with the increased incidence of diabetes.

7. Conclusions and Recommendations

In this study, we considered a nine-compartment, mutually exclusive compartmental model of the comorbid infection of tuberculosis and diabetes. Our focus was on the impact of saturated treatment with LTBIs on tuberculosis prevention and control in a population that was divided into diabetic and non-diabetic subgroups. Additionally, we analyzed the consequences of diabetes and exogenous reinfections on the prevalence and incidence of tuberculosis infections in the community.

Our theoretical and epidemiological results indicate that Model (1) possesses a tuberculosis-free equilibrium point (TBFE) that is locally asymptotically stable when ${\mathcal{R}}_0<1$ and globally asymptotically stable when there is no exogenous reinfection or saturated treatment of those who harbor the latent form of TB infections. Model (1) also exhibits a phenomenon known as backward bifurcation due to the saturated treatment of LTBIs and exogenous reinfection. In this scenario, maintaining the crucial threshold ${\mathcal{R}}_0<1$ is no longer sufficient for controlling tuberculosis infections in the population. We determined the endemic equilibrium of a specific case of Model (1) and established the existence of the tuberculosis-present (endemic) equilibrium points.

The model parameters that influence the basic reproduction number (${\mathcal{R}}_0$) significantly are the rates of developing diabetes without or with complications (${\sigma }_1,{\sigma }_2)$, the effective contact rate ($\omega$), the treatment rate of tuberculosis among diabetes with complication patients ($c_3$), the saturated treatment rate of LTBIs ($q_1$), and the modification parameter that takes into consideration the high rate of getting TB in its latent form in those sick with complicated and uncomplicated diabetes $(\varphi ,\phi )$.

The findings from the numerical simulations and sensitivity analysis of Model (1) revealed that a higher incidence of diabetes mellitus is a risk factor for the elevated incidence and prevalence of TB in the population. As the number of diabetic individuals in society increases, it becomes more challenging to effectively control the spread of tuberculosis. Additionally, we discovered that the saturated treatment of LTBIs in the population is a successful end TB strategy to adopt for the purpose of achieving the WHO target of ending the global tuberculosis (TB) epidemic by the year 2035. Furthermore, we established that without a saturated treatment regimen for LTBIs, tuberculosis infections will persist. Therefore, providing the saturated treatment of LTBIs will also alleviate the double burden of the comorbidity of TB and DM in the community.

It is strongly recommended that diabetic patients be screened and treated for latent tuberculosis infection (LTBI) to prevent its progression to active tuberculosis disease. This approach will help reduce the spread of tuberculosis, morbidity, and mortality [8,51]. The converging epidemic of infectious diseases like tuberculosis and non-infectious diseases (e.g., diabetes mellitus) also places a huge burden on the concerted efforts to effectively control the TB disease in the population. The direct impact of diabetes on active tuberculosis and tuberculosis relapses is a major concern. Furthermore, poor glycemic control among people with diabetes mellitus can result in a high prevalence of tuberculosis in the population [2,9,51]. According to Lee et al. [9], good glycemic control may decrease the risk of tuberculosis among people with diabetes mellitus and contribute to tuberculosis control in populations with a high prevalence of both tuberculosis and diabetes.

In future research endeavors, additional treatment approaches can be integrated into our model system (1) and in-depth analyses conducted. The aim will be to evaluate and compare the most effective treatment strategies for LTBIs in diabetic and non-diabetic populations. This information will be valuable in guiding policymakers and public health officials in the development of effective policies to combat the spread of TB. We also intend to perform optimal control and cost-effectiveness analyses of the model. Furthermore, other mathematical and epidemiological tools can be employed to investigate the interactions between drug-resistant TB and diabetes mellitus, while taking into account important factors such as sex and the influence of age on the development of diabetes [52].