<figure><img src="image_1.png" /></figure>

HAL
open science

# Analysis of the Dynamics of Tuberculosis in Algeria Using a Compartmental VSEIT Model with Evaluation of the Vaccination and Treatment Effects

Bouchra Chennaf, Mohammed Salah Abdelouahab, René Lozi

► To cite this version:

Bouchra Chennaf, Mohammed Salah Abdelouahab, René Lozi. Analysis of the Dynamics of Tuberculosis in Algeria Using a Compartmental VSEIT Model with Evaluation of the Vaccination and Treatment Effects. Computation, 2023, 11 (7), pp.146. 10.3390/computation11070146. hal-04549426

HAL Id: hal-04549426

https://hal.science/hal-04549426

Submitted on 17 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Analysis of the Dynamics of Tuberculosis in Algeria Using a Compartmental VSEIT Model with Evaluation of the Vaccination and Treatment Effects

Bouchra Chennaf <fnref n="1" /> ID, Mohammed Salah Abdelouahab <fnref n="1" />,* ID and René Lozi <fnref n="2" /> ID

**Abstract:** Despite low tuberculosis (TB) mortality rates in China, Europe, and the United States, many countries are still struggling to control the epidemic, including India, South Africa, and Algeria. This study aims to contribute to the body of knowledge on this topic and provide a valuable tool and evidence-based guidance for the Algerian healthcare managers in understanding the spread of TB and implementing control strategies. For this purpose, a compartmental mathematical model is proposed to analyze TB dynamics in Algeria and investigate the vaccination and treatment effects on disease breaks. A qualitative study is conducted to discuss the stability property of both disease-free equilibrium and endemic equilibrium. In order to adopt the proposed model for the Algerian case, we estimate the model parameters using Algerian TB-reported data from 1990 to 2020. The obtained results using the proposed mathematical compartmental model show that the reproduction number ($R_0$) of TB in Algeria is less than one, suggesting that the disease can be eradicated or effectively controlled through a combination of interventions, including vaccination, high-quality treatment, and isolation measures.

**Keywords:** tuberculosis model; epidemic; vaccination; parameter estimation

<figure><img src="image_1.png" /></figure>

check for
updates

**Citation:** Chennaf, B.; Abdelouahab, M.S.; Lozi, R. Analysis of the Dynamics of Tuberculosis in Algeria Using a Compartmental VSEIT Model with Evaluation of the Vaccination and Treatment Effects. *Computation* **2023**, *11*, 146. https://doi.org/10.3390/computation11070146

Academic Editors: Demos T. Tsahalis and Alexander Pchelintsev

Received: 9 June 2023

Revised: 5 July 2023

Accepted: 16 July 2023

Published: 21 July 2023

<figure><img src="image_2.png" /></figure>

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The recent outbreak of the coronavirus disease, known as COVID-19 has indeed highlighted the critical role of epidemic research, particularly mathematical modeling, in understanding and combating infectious diseases, since it provides a powerful tool to analyze the dynamic transmission of diseases and assess the potential impact of various interventions and control measures.

Tuberculosis is a contagious infection caused by bacteria called Mycobacterium tuberculosis that primarily affects the lungs. It can also spread to other body parts, including the brain and spine. Importantly, it can be contracted not only through direct contact with an infected individual but also through the inhalation of airborne droplets containing the bacteria.

TB is one of the top 10 killers worldwide and causes 1.8 million deaths each year. Of all new TB cases recorded in 2020, 86% occurred in the 30 countries with the highest disease burden. Two-thirds of the cases are concentrated in eight countries, with India leading, followed by China, Indonesia, Philippines, Pakistan, Nigeria, Bangladesh and South Africa (According to the World Health Organization (WHO)) [1]. This demonstrates that TB poses a threat to human health and has a detrimental impact on social and economic life.

Although Algeria may not be among the top eight countries with the highest concentration of TB cases globally, it is still a significant concern in Algeria. Thus, it is imperative that government agencies and scientists work together to manage and combat the spread of TB epidemics.

<footnote n="1">¹ Laboratory of Mathematics and Their Interactions, Abdelhafid Boussouf University Center, Mila 43000, Algeria; bouchra.chennaf@centre-univ-mila.dz</footnote>

<footnote n="2">² Laboratory J.A. Dieudonné, CNRS, Université Côte d'Azur, 06108 Nice, France; rene.lozi@univ-cotedazur.fr</footnote>

<footnote>* Correspondence: m.abdelouahab@centre-univ-mila.dz</footnote>
<figure><img src="image_3.png" /></figure>

Article

# Analysis of the Dynamics of Tuberculosis in Algeria Using a Compartmental VSEIT Model with Evaluation of the Vaccination and Treatment Effects

Bouchra Chennaf <fnref n="1" /> ID, Mohammed Salah Abdelouahab <fnref n="1" />,* ID and René Lozi <fnref n="2" /> ID

**Abstract:** Despite low tuberculosis (TB) mortality rates in China, Europe, and the United States, many countries are still struggling to control the epidemic, including India, South Africa, and Algeria. This study aims to contribute to the body of knowledge on this topic and provide a valuable tool and evidence-based guidance for the Algerian healthcare managers in understanding the spread of TB and implementing control strategies. For this purpose, a compartmental mathematical model is proposed to analyze TB dynamics in Algeria and investigate the vaccination and treatment effects on disease breaks. A qualitative study is conducted to discuss the stability property of both disease-free equilibrium and endemic equilibrium. In order to adopt the proposed model for the Algerian case, we estimate the model parameters using Algerian TB-reported data from 1990 to 2020. The obtained results using the proposed mathematical compartmental model show that the reproduction number ($R_0$) of TB in Algeria is less than one, suggesting that the disease can be eradicated or effectively controlled through a combination of interventions, including vaccination, high-quality treatment, and isolation measures.

**Keywords:** tuberculosis model; epidemic; vaccination; parameter estimation

<figure><img src="image_4.png" /></figure>

check for
updates

**Citation:** Chennaf, B.; Abdelouahab, M.S.; Lozi, R. Analysis of the Dynamics of Tuberculosis in Algeria Using a Compartmental VSEIT Model with Evaluation of the Vaccination and Treatment Effects. *Computation* **2023**, *11*, 146. https://doi.org/10.3390/computation11070146

Academic Editors: Demos T. Tsahalis and Alexander Pchelintsev

Received: 9 June 2023

Revised: 5 July 2023

Accepted: 16 July 2023

Published: 21 July 2023

<figure><img src="image_5.png" /></figure>

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The recent outbreak of the coronavirus disease, known as COVID-19 has indeed highlighted the critical role of epidemic research, particularly mathematical modeling, in understanding and combating infectious diseases, since it provides a powerful tool to analyze the dynamic transmission of diseases and assess the potential impact of various interventions and control measures.

Tuberculosis is a contagious infection caused by bacteria called Mycobacterium tuberculosis that primarily affects the lungs. It can also spread to other body parts, including the brain and spine. Importantly, it can be contracted not only through direct contact with an infected individual but also through the inhalation of airborne droplets containing the bacteria.

TB is one of the top 10 killers worldwide and causes 1.8 million deaths each year. Of all new TB cases recorded in 2020, 86% occurred in the 30 countries with the highest disease burden. Two-thirds of the cases are concentrated in eight countries, with India leading, followed by China, Indonesia, Philippines, Pakistan, Nigeria, Bangladesh and South Africa (According to the World Health Organization (WHO)) [1]. This demonstrates that TB poses a threat to human health and has a detrimental impact on social and economic life.

Although Algeria may not be among the top eight countries with the highest concentration of TB cases globally, it is still a significant concern in Algeria. Thus, it is imperative that government agencies and scientists work together to manage and combat the spread of TB epidemics.

<footnote n="1">¹ Laboratory of Mathematics and Their Interactions, Abdelhafid Boussouf University Center, Mila 43000, Algeria; bouchra.chennaf@centre-univ-mila.dz</footnote>

<footnote n="2">² Laboratory J.A. Dieudonné, CNRS, Université Côte d'Azur, 06108 Nice, France; rene.lozi@univ-cotedazur.fr</footnote>

<footnote>* Correspondence: m.abdelouahab@centre-univ-mila.dz</footnote>
Mathematical modeling plays a critical role in the planning and implementation of TB control programs. Although Bernoulli used mathematical models for smallpox in 1760 [2], the research on infectious diseases using deterministic mathematical models actually started in the 20th century. Other major works in mathematical epidemiology are due to P.D. En'ko between 1873 and 1894. However, it can be said that the foundations of mathematical epidemiology based on compartmental models are due to Sir Ronald Ross, who gave the first mathematical model of malaria transmission in 1911 [3]. The early research in this field are available in [4–6].

The main difference between compartmental models and other models of diseases is that compartmental models explicitly consider the different stages of disease progression and the transitions between them. This allows for a more detailed understanding of how diseases spread and how interventions can be implemented to control their spread [7]. Other models, such as statistical models or network models, may not explicitly include this level of detail.

Susceptible-Infected-Recovered (SIR) is a deterministic model that Kermack and McKendrick proposed in 1927 to characterize the behavior of epidemic spread [8]. Despite the fact that this model has been used successfully to represent the behavior of disease, it is unrealistic by ignoring other compartments and control techniques, such as vaccination, treatment, isolation, and the impact of age and sex.

Epidemiological compartmental models can be broadly classified into two categories: differential equation-based models that describe the dynamics of infectious diseases using continuous functions which can capture the continuous changes in the state variables over time, typically represented by systems of ordinary differential equations and difference equation-based models, often used when data are collected at discrete time points or when the population dynamics are better captured in a discrete manner.

The first mathematical model of TB was developed in 1962 by Waaler and Anderson, who divided the entire population into different groups [9]. Since then, numerous academics have created various mathematical models to investigate and control TB in countries heavily impacted by the disease; see for instance [10,11]. These models have been instrumental in guiding public health policies and interventions aimed at reducing the burden of TB.

Vaccination is one of the most vital factors in stopping and controlling the spread of TB. The tuberculosis vaccination against Bacillus Calmette–Guérin (BCG) was first given to a human in 1921. The World Health Organization (WHO) currently advises immunizing newborns with a single intradermal injection of BCG as soon as possible after the birth [12]. To address the prevailing epidemiological situation, the Algerian Health Care Administration implemented a dedicated vaccination schedule and mandatory vaccination campaigns for children. As a result, the BCG vaccination coverage reached a remarkable rate of over 98% across newborns. However, BCG vaccination is not typically recommended for adults, as its effectiveness in this age group is limited. Therefore, this paper neglects it.

The mathematical modeling of tuberculosis relies on vaccination for its importance in giving predictions to eradicate the disease. There are many previous studies concerned with this topic; for example, in [13] the goal was to determine the dynamics of tuberculosis in Turkey, and the impact of vaccine therapy on the disease. Yang et al. [7] formulated a mathematical model to investigate the effects of immunization and treatment on the dynamics of tuberculosis transmission. Egonmwan et al. [14] developed a mathematical model that includes immunization of newborn children and older susceptible people in the dynamics of TB transmission in a population, with the goal of providing protection to older susceptible people. Revelle et al. [15] formulated models for the economic allocation of activities to control tuberculosis in developing countries.

To the best of our knowledge, the proposed model is not considered elsewhere in its present form and there is no research on modeling the dynamic transmission of tuberculosis in Algeria using a compartmental model while simultaneously estimating the relevant
Mathematical modeling plays a critical role in the planning and implementation of TB control programs. Although Bernoulli used mathematical models for smallpox in 1760 [2], the research on infectious diseases using deterministic mathematical models actually started in the 20th century. Other major works in mathematical epidemiology are due to P.D. En'ko between 1873 and 1894. However, it can be said that the foundations of mathematical epidemiology based on compartmental models are due to Sir Ronald Ross, who gave the first mathematical model of malaria transmission in 1911 [3]. The early research in this field are available in [4–6].

The main difference between compartmental models and other models of diseases is that compartmental models explicitly consider the different stages of disease progression and the transitions between them. This allows for a more detailed understanding of how diseases spread and how interventions can be implemented to control their spread [7]. Other models, such as statistical models or network models, may not explicitly include this level of detail.

Susceptible-Infected-Recovered (SIR) is a deterministic model that Kermack and McKendrick proposed in 1927 to characterize the behavior of epidemic spread [8]. Despite the fact that this model has been used successfully to represent the behavior of disease, it is unrealistic by ignoring other compartments and control techniques, such as vaccination, treatment, isolation, and the impact of age and sex.

Epidemiological compartmental models can be broadly classified into two categories: differential equation-based models that describe the dynamics of infectious diseases using continuous functions which can capture the continuous changes in the state variables over time, typically represented by systems of ordinary differential equations and difference equation-based models, often used when data are collected at discrete time points or when the population dynamics are better captured in a discrete manner.

The first mathematical model of TB was developed in 1962 by Waaler and Anderson, who divided the entire population into different groups [9]. Since then, numerous academics have created various mathematical models to investigate and control TB in countries heavily impacted by the disease; see for instance [10,11]. These models have been instrumental in guiding public health policies and interventions aimed at reducing the burden of TB.

Vaccination is one of the most vital factors in stopping and controlling the spread of TB. The tuberculosis vaccination against Bacillus Calmette–Guérin (BCG) was first given to a human in 1921. The World Health Organization (WHO) currently advises immunizing newborns with a single intradermal injection of BCG as soon as possible after the birth [12]. To address the prevailing epidemiological situation, the Algerian Health Care Administration implemented a dedicated vaccination schedule and mandatory vaccination campaigns for children. As a result, the BCG vaccination coverage reached a remarkable rate of over 98% across newborns. However, BCG vaccination is not typically recommended for adults, as its effectiveness in this age group is limited. Therefore, this paper neglects it.

The mathematical modeling of tuberculosis relies on vaccination for its importance in giving predictions to eradicate the disease. There are many previous studies concerned with this topic; for example, in [13] the goal was to determine the dynamics of tuberculosis in Turkey, and the impact of vaccine therapy on the disease. Yang et al. [7] formulated a mathematical model to investigate the effects of immunization and treatment on the dynamics of tuberculosis transmission. Egonmwan et al. [14] developed a mathematical model that includes immunization of newborn children and older susceptible people in the dynamics of TB transmission in a population, with the goal of providing protection to older susceptible people. Revelle et al. [15] formulated models for the economic allocation of activities to control tuberculosis in developing countries.

To the best of our knowledge, the proposed model is not considered elsewhere in its present form and there is no research on modeling the dynamic transmission of tuberculosis in Algeria using a compartmental model while simultaneously estimating the relevant
biological parameters specific to the country. Therefore, conducting research in this domain would make a valuable contribution to the field of TB modeling and control in Algeria.

In this study, we propose a VSEIT epidemiological model to investigate the dynamics of TB disease in Algeria. To confirm its performance we estimated the biological model’s pa- rameters using specific TB data, including disease incidence, prevalence, and other relevant epidemiological information from 1990 to 2020 from the WHO Global TB Report [1].

This paper is organized as follows: Section 2 presents the formulation of the VSEIT TB model and analyzes its dynamic properties. In Section 3, the estimation of model parameters is conducted, along with their sensitivity analysis. Section 4 is focused on the discussion of the obtained results. Finally, the paper is concluded in Section 5.

## 2. Mathematical Model and Dynamic Analysis

In this section, we present the proposed mathematical model for TB infection, and we examine its dynamic.

### 2.1. Model Formulation

The population is divided into five distinct subgroups: vaccinated individuals (V), susceptible individuals (S), exposed or exposed individuals (E), infected individuals with active TB (I), and individuals currently receiving treatment (T). Hence, the total population is

$$N(t) = V(t) + S(t) + E(t) + I(t) + T(t).$$

This model aims to provide a comprehensive understanding of the spread and pro- gression of TB within a population, allowing for more effective prevention and treatment strategies to be developed.

The number of people that have received vaccination (V) is increased through a small proportion of immunized newborns, $p\Lambda$. The vaccinated population decreases as vaccinated individuals become susceptible at a rate $k$ (the vaccine’s efficacy wanes over time), during the protection period, they will not become infected even if they contact infected individuals as long as the vaccination provides immunity to all of them. The natural death rate in the class V is $\mu$. Hence, the population of vaccinated individuals is given by the first equation in system (1).

The population of susceptible individuals (S) is increased by the small proportion of newborns who are not immunized from TB, $(1-p)\Lambda$, and also increases as vaccinated individuals become susceptible, at a rate $k$. This population decreases when there is contact with infected people, at a rate $\beta$. As older people die naturally at a rate $\mu$, the population of susceptible individuals decreases. Hence, the population of susceptible individuals is given by the second equation in system (1).

We assume that the population of exposed individuals E increases when the suscepti- ble population makes effective contact with infected individuals and decreases as latently infected peoples progress from exposed to active TB, at a rate $\epsilon$, and die naturally, at a rate $\mu$. The population of exposed individuals E is also increased by an inflow of a fraction, $\delta(1 - \alpha)$ of individuals under treatment, where the parameter $\alpha$ represents the treatment failure rate. Particularly, $\alpha = 0$ means that all treated individuals will move to an exposed state, whereas $\alpha = 1$, means that the treatment has failed, and all treated individuals will remain infectious. Hence, the population of exposed individuals is given by the third equation in system (1).

The population of infected individuals (I) increases as latently infected individuals progress from exposed to active TB, and the effectively treated patients return to active TB at a rate $a\delta$, which significantly increases the population I. As infected people receive treatment, the population I decreases at a rate $\gamma$. Both natural death and TB disease kill people at rates $\sigma$ and $\mu$, respectively. Hence, the population of infected individuals is given by the fourth equation in system (1).

Finally, as infected people are treated, the population of treated individuals (T) grows at a rate $\gamma$. As individuals who have been successfully treated become reinfected, the
biological parameters specific to the country. Therefore, conducting research in this domain would make a valuable contribution to the field of TB modeling and control in Algeria.

In this study, we propose a VSEIT epidemiological model to investigate the dynamics of TB disease in Algeria. To confirm its performance we estimated the biological model’s parameters using specific TB data, including disease incidence, prevalence, and other relevant epidemiological information from 1990 to 2020 from the WHO Global TB Report [1].

This paper is organized as follows: Section 2 presents the formulation of the VSEIT TB model and analyzes its dynamic properties. In Section 3, the estimation of model parameters is conducted, along with their sensitivity analysis. Section 4 is focused on the discussion of the obtained results. Finally, the paper is concluded in Section 5.

## 2. Mathematical Model and Dynamic Analysis

In this section, we present the proposed mathematical model for TB infection, and we examine its dynamic.

### 2.1. Model Formulation

The population is divided into five distinct subgroups: vaccinated individuals (V), susceptible individuals (S), exposed or exposed individuals (E), infected individuals with active TB (I), and individuals currently receiving treatment (T). Hence, the total population is

$$N(t) = V(t) + S(t) + E(t) + I(t) + T(t).$$

This model aims to provide a comprehensive understanding of the spread and progression of TB within a population, allowing for more effective prevention and treatment strategies to be developed.

The number of people that have received vaccination (V) is increased through a small proportion of immunized newborns, $p\Lambda$. The vaccinated population decreases as vaccinated individuals become susceptible at a rate $k$ (the vaccine’s efficacy wanes over time), during the protection period, they will not become infected even if they contact infected individuals as long as the vaccination provides immunity to all of them. The natural death rate in the class V is $\mu$. Hence, the population of vaccinated individuals is given by the first equation in system (1).

The population of susceptible individuals (S) is increased by the small proportion of newborns who are not immunized from TB, $(1-p)\Lambda$, and also increases as vaccinated individuals become susceptible, at a rate $k$. This population decreases when there is contact with infected people, at a rate $\beta$. As older people die naturally at a rate $\mu$, the population of susceptible individuals decreases. Hence, the population of susceptible individuals is given by the second equation in system (1).

We assume that the population of exposed individuals E increases when the susceptible population makes effective contact with infected individuals and decreases as latently infected peoples progress from exposed to active TB, at a rate $\epsilon$, and die naturally, at a rate $\mu$. The population of exposed individuals E is also increased by an inflow of a fraction, $\delta(1-\alpha)$ of individuals under treatment, where the parameter $\alpha$ represents the treatment failure rate. Particularly, $\alpha = 0$ means that all treated individuals will move to an exposed state, whereas $\alpha = 1$, means that the treatment has failed, and all treated individuals will remain infectious. Hence, the population of exposed individuals is given by the third equation in system (1).

The population of infected individuals (I) increases as latently infected individuals progress from exposed to active TB, and the effectively treated patients return to active TB at a rate $a\delta$, which significantly increases the population I. As infected people receive treatment, the population I decreases at a rate $\gamma$. Both natural death and TB disease kill people at rates $\sigma$ and $\mu$, respectively. Hence, the population of infected individuals is given by the fourth equation in system (1).

Finally, as infected people are treated, the population of treated individuals (T) grows at a rate $\gamma$. As individuals who have been successfully treated become reinfected, the
population *T* decreases at a rate *δ*. The population continues to decline due to natural mortality (*μ*) and deaths from TB (*η*). Hence, the population of treated individuals is given by the fifth equation in system (1).

The dynamic of TB infection is described by the following system of differential equations:

$$
\left\{
\begin{aligned}
\frac{dV(t)}{dt} &= p\Lambda - (k + \mu)V(t), \\
\frac{dS(t)}{dt} &= (1 - p)\Lambda + kV(t) - \beta S(t)I(t) - \mu S(t), \\
\frac{dE(t)}{dt} &= \beta S(t)I(t) - (\epsilon + \mu)E(t) + (1 - \alpha)\delta T(t), \\
\frac{dI(t)}{dt} &= \epsilon E(t) + \alpha\delta T(t) - (\gamma + \mu + \sigma)I(t), \\
\frac{dT(t)}{dt} &= \gamma I(t) - (\mu + \delta + \eta)T(t).
\end{aligned}
\right.
\qquad (1)
$$

With: $V(0) \ge 0$, $S(0) \ge 0$, $E(0) \ge 0$, $I(0) \ge 0$ and $T(0) \ge 0$ with $N(0) > 0$.

The flowchart of the model is shown in Figure 1. The model variables are presented in Table 1 and the model parameters are presented in Table 2.

<figure><img src="image_3.png" /><figcaption>Figure 1. Flowchart of VSEIT model.</figcaption></figure>

<table><caption>**Table 1.** Description of variables of the model (1).</caption><thead><tr><td>Variable</td><td>Description</td></tr></thead><tbody><tr><td>V(t)</td><td>The vaccinated population at time t.</td></tr><tr><td>S(t)</td><td>The susceptible population which is able to be infected at any time t.</td></tr><tr><td>E(t)</td><td>The exposed population which is not yet infectious.</td></tr><tr><td>I(t)</td><td>The infected population at time t.</td></tr><tr><td>T(t)</td><td>The treated population at time t.</td></tr></tbody></table>

The values of parameters in Table 2 will be given in Section 3.1.
population *T* decreases at a rate *δ*. The population continues to decline due to natural mortality (*μ*) and deaths from TB (*η*). Hence, the population of treated individuals is given by the fifth equation in system (1).

The dynamic of TB infection is described by the following system of differential equations:

$$
\left\{
\begin{aligned}
\frac{dV(t)}{dt} &= p\Lambda - (k + \mu)V(t), \\
\frac{dS(t)}{dt} &= (1 - p)\Lambda + kV(t) - \beta S(t)I(t) - \mu S(t), \\
\frac{dE(t)}{dt} &= \beta S(t)I(t) - (\epsilon + \mu)E(t) + (1 - \alpha)\delta T(t), \\
\frac{dI(t)}{dt} &= \epsilon E(t) + \alpha\delta T(t) - (\gamma + \mu + \sigma)I(t), \\
\frac{dT(t)}{dt} &= \gamma I(t) - (\mu + \delta + \eta)T(t).
\end{aligned}
\right.
\qquad (1)
$$

With: $V(0) \ge 0$, $S(0) \ge 0$, $E(0) \ge 0$, $I(0) \ge 0$ and $T(0) \ge 0$ with $N(0) > 0$.

The flowchart of the model is shown in Figure 1. The model variables are presented in Table 1 and the model parameters are presented in Table 2.

<figure><img src="image_6.png" /><figcaption>Figure 1. Flowchart of VSEIT model.</figcaption></figure>

<table><caption>**Table 1.** Description of variables of the model (1).</caption><thead><tr><td>Variable</td><td>Description</td></tr></thead><tbody><tr><td>V(t)</td><td>The vaccinated population at time t.</td></tr><tr><td>S(t)</td><td>The susceptible population which is able to be infected at any time t.</td></tr><tr><td>E(t)</td><td>The exposed population which is not yet infectious.</td></tr><tr><td>I(t)</td><td>The infected population at time t.</td></tr><tr><td>T(t)</td><td>The treated population at time t.</td></tr></tbody></table>

The values of parameters in Table 2 will be given in Section 3.1.
<table><caption>**Table 2.** Parameters of model (1).</caption><thead><tr><th>Model Parameter</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>Λ</td><td>Recruitment rate</td><td>year<sup>-1</sup></td></tr><tr><td>μ</td><td>Natural death rate</td><td>year<sup>-1</sup></td></tr><tr><td>k</td><td>Rate of moving from V to S</td><td>year<sup>-1</sup></td></tr><tr><td>β</td><td>Transmission rate</td><td>year<sup>-1</sup></td></tr><tr><td>γ</td><td>Treatment rate</td><td>year<sup>-1</sup></td></tr><tr><td>ε</td><td>Progression rate</td><td>year<sup>-1</sup></td></tr><tr><td>α</td><td>Treatment failure rate</td><td>year<sup>-1</sup></td></tr><tr><td>δ</td><td>Rate at which the treated population leave the class T</td><td>year<sup>-1</sup></td></tr><tr><td>σ</td><td>Disease death rate in I</td><td>year<sup>-1</sup></td></tr><tr><td>η</td><td>Disease death rate in T</td><td>year<sup>-1</sup></td></tr><tr><td>p</td><td>Vaccination rate</td><td>year<sup>-1</sup></td></tr></tbody></table>

## 2.2. Model Analysis

In this subsection the proposed model (1), will be qualitatively analyzed.

### 2.2.1. Invariance of the Feasible Region

The TB model (1) will be studied in a biologically feasible region $\Omega \subset \mathbb{R}_+^5$ given by

$$ \Omega = \{ (V(t), S(t), E(t), I(t), T(t)) \in \mathbb{R}_+^5 : N(t) \le \frac{\Lambda}{\mu} \}. \quad (2) $$

**Lemma 1.** For all $t > 0$ and non-negative initial conditions, the solution of TB model (1) is positive whenever it exists. Furthermore, if $0 \le N(0) \le \frac{\Lambda}{\mu}$, then

$$ 0 \le N(t) \le \frac{\Lambda}{\mu}, \text{ for all } t > 0. $$

**Proof.** For positive values of $V(t), S(t), E(t), I(t)$ and $T(t)$ we have

$$ \left\{ \begin{array}{l} V'|_{V=0} = p\Lambda \ge 0, \\ S'|_{S=0} = (1-p)\Lambda + kV \ge 0, \\ E'|_{E=0} = \beta SI + (1-\alpha)\delta T \ge 0, \\ I'|_{I=0} = \epsilon E + \alpha \delta T \ge 0, \\ T'|_{T=0} = \gamma I \ge 0. \end{array} \right. \qquad (3) $$

Hence, for non-negative initial conditions, the solution remains positive $\forall t \ge 0$.

It follows from the addition of the VSEIT model Equations (1) that

$$ \frac{dN(t)}{dt} = \Lambda - \mu(V(t) + S(t) + E(t) + I(t) + T(t)) - (\sigma I(t) + \eta T(t)), \quad (4) $$

$$ = \Lambda - \mu N(t) - (\sigma I(t) + \eta T(t)) \le \Lambda - \mu N(t). \quad (5) $$

For $N(t) \le \frac{\Lambda}{\mu}$, we have $\frac{dN(t)}{dt} \le 0$.

Thus, for $0 \le N(0) \le \frac{\Lambda}{\mu}$, we obtain $0 \le N(t) \le \frac{\Lambda}{\mu}$ for all $t \ge 0$.

It follows that the feasible region $\Omega$ is positively invariant. □
<table><caption>**Table 2.** Parameters of model (1).</caption><thead><tr><th>Model Parameter</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>Λ</td><td>Recruitment rate</td><td>year<sup>-1</sup></td></tr><tr><td>μ</td><td>Natural death rate</td><td>year<sup>-1</sup></td></tr><tr><td>k</td><td>Rate of moving from V to S</td><td>year<sup>-1</sup></td></tr><tr><td>β</td><td>Transmission rate</td><td>year<sup>-1</sup></td></tr><tr><td>γ</td><td>Treatment rate</td><td>year<sup>-1</sup></td></tr><tr><td>ε</td><td>Progression rate</td><td>year<sup>-1</sup></td></tr><tr><td>α</td><td>Treatment failure rate</td><td>year<sup>-1</sup></td></tr><tr><td>δ</td><td>Rate at which the treated population leave the class T</td><td>year<sup>-1</sup></td></tr><tr><td>σ</td><td>Disease death rate in I</td><td>year<sup>-1</sup></td></tr><tr><td>η</td><td>Disease death rate in T</td><td>year<sup>-1</sup></td></tr><tr><td>p</td><td>Vaccination rate</td><td>year<sup>-1</sup></td></tr></tbody></table>

## 2.2. Model Analysis

In this subsection the proposed model (1), will be qualitatively analyzed.

### 2.2.1. Invariance of the Feasible Region

The TB model (1) will be studied in a biologically feasible region $\Omega \subset \mathbb{R}_+^5$ given by

$$ \Omega = \{ (V(t), S(t), E(t), I(t), T(t)) \in \mathbb{R}_+^5 : N(t) \le \frac{\Lambda}{\mu} \}. \quad (2) $$

**Lemma 1.** For all $t > 0$ and non-negative initial conditions, the solution of TB model (1) is positive whenever it exists. Furthermore, if $0 \le N(0) \le \frac{\Lambda}{\mu}$, then

$$ 0 \le N(t) \le \frac{\Lambda}{\mu}, \text{ for all } t > 0. $$

**Proof.** For positive values of $V(t), S(t), E(t), I(t)$ and $T(t)$ we have

$$ \left\{ \begin{array}{l} V'|_{V=0} = p\Lambda \ge 0, \\ S'|_{S=0} = (1-p)\Lambda + kV \ge 0, \\ E'|_{E=0} = \beta SI + (1-\alpha)\delta T \ge 0, \\ I'|_{I=0} = \epsilon E + \alpha \delta T \ge 0, \\ T'|_{T=0} = \gamma I \ge 0. \end{array} \right. \qquad (3) $$

Hence, for non-negative initial conditions, the solution remains positive $\forall t \ge 0$.

It follows from the addition of the VSEIT model Equations (1) that

$$ \frac{dN(t)}{dt} = \Lambda - \mu(V(t) + S(t) + E(t) + I(t) + T(t)) - (\sigma I(t) + \eta T(t)), \quad (4) $$

$$ = \Lambda - \mu N(t) - (\sigma I(t) + \eta T(t)) \le \Lambda - \mu N(t). \quad (5) $$

For $N(t) \le \frac{\Lambda}{\mu}$, we have $\frac{dN(t)}{dt} \le 0$.

Thus, for $0 \le N(0) \le \frac{\Lambda}{\mu}$, we obtain $0 \le N(t) \le \frac{\Lambda}{\mu}$ for all $t \ge 0$.

It follows that the feasible region $\Omega$ is positively invariant. □
### 2.2.2. Equilibrium Points and Their Stability

To find equilibrium points of the model (1), one solves the equations:

$$ \frac{dV}{dt} = \frac{dS}{dt} = \frac{dE}{dt} = \frac{dI}{dt} = \frac{dT}{dt} = 0. $$

One gets two equilibrium points:
The disease-free equilibrium point “DFE”

$$ E_1 = (V_1^*, S_1^*, E_1^*, I_1^*, T_1^*) = \left( \frac{p\Lambda}{k+\mu}, \frac{(k+\mu-\mu p)\Lambda}{\mu(k+\mu)}, 0, 0, 0 \right), \text{ at which } N = V_1^* + S_1^* = \frac{\Lambda}{\mu}. $$

and the endemic equilibrium point “EE”

$$ E_2 = (V_2^*, S_2^*, E_2^*, I_2^*, T_2^*) = \left( \frac{p\Lambda}{k+\mu}, \frac{(k+\mu-\mu p)\Lambda}{(k+\mu)(\beta I_2^* + \mu)}, \frac{(\gamma+\mu+\sigma)(\mu+\delta+\eta) - \alpha\delta\gamma}{\epsilon(\mu+\delta+\eta)} I_2^*, I_2^*, \frac{\gamma}{\mu+\delta+\eta} I_2^* \right), $$

where:

$$ I_2^* = \frac{(k+\mu-\mu p)\epsilon\Lambda(\mu+\delta+\eta)}{(k+\mu)((\epsilon+\mu)(\gamma+\mu+\sigma)(\mu+\delta+\eta) - (\epsilon+\mu)\alpha\delta\gamma - (1-\alpha)\gamma\delta\epsilon)} - \frac{\mu}{\beta}. $$ $$ = \frac{\mu}{\beta}(\mathcal{R}_0 - 1). $$

Then, the endemic equilibrium point $E_2$, exists for $\mathcal{R}_0 > 1$, and at this equilibrium, one has $N = \frac{\Lambda - (\sigma I_2^* + \eta T_2^*)}{\mu} < \frac{\Lambda}{\mu}$.

### 2.2.3. The Basic Reproduction Number $\mathcal{R}_0$

The basic reproduction number, denoted $\mathcal{R}_0$, is the expected number of secondary cases produced in a completely susceptible population, by a typical infectious individual during its infective period [16].

If $\mathcal{R}_0 < 1$, the disease will not be able to spread among the population, whereas, if $\mathcal{R}_0 > 1$, the disease has the potential to spread among the population and become endemic.

Using the next generation matrix one gets the basic reproduction number $\mathcal{R}_0$ for the proposed model (1) as

$$ \mathcal{R}_0 = \frac{\epsilon(k + \mu - \mu p)\Lambda\beta k_3}{\mu(k + \mu)(k_1 k_2 k_3 - \alpha\gamma\delta k_1 - (1 - \alpha)\delta\gamma\epsilon)}. $$

The details of the calculation of $\mathcal{R}_0$ are presented in the Appendix A and B.

### 2.2.4. Local Stability Analysis of DFE

**Theorem 1.** The disease-free equilibrium point $E_1$ (DFE) of model (1) is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

**Proof.** The Jacobian matrix $J_{E_1}$ of system (1) at the DFE $E_1$ is given by

$$ J_{E_1} = \begin{bmatrix} -(k+\mu) & 0 & 0 & 0 & 0 \\ k & -\mu & 0 & -\beta S_1^* & 0 \\ 0 & 0 & -k_1 & \beta S_1^* & (1-\alpha)\delta \\ 0 & 0 & \epsilon & -k_2 & \alpha\delta \\ 0 & 0 & 0 & \gamma & -k_3 \end{bmatrix} $$

The characteristic equation of $J_{E_1}$ is

$$ -(k+\mu)-\lambda)(-\mu-\lambda)[\lambda^3 + a_1\lambda^2 + a_2\lambda + a_3] = 0, \quad (6) $$
### 2.2.2. Equilibrium Points and Their Stability

To find equilibrium points of the model (1), one solves the equations:

$$ \frac{dV}{dt} = \frac{dS}{dt} = \frac{dE}{dt} = \frac{dI}{dt} = \frac{dT}{dt} = 0. $$

One gets two equilibrium points:
The disease-free equilibrium point “DFE”

$$ E_1 = (V_1^*, S_1^*, E_1^*, I_1^*, T_1^*) = \left( \frac{p\Lambda}{k+\mu}, \frac{(k+\mu-\mu p)\Lambda}{\mu(k+\mu)}, 0, 0, 0 \right), \text{ at which } N = V_1^* + S_1^* = \frac{\Lambda}{\mu}. $$

and the endemic equilibrium point “EE”

$$ \begin{aligned} E_2 &= (V_2^*, S_2^*, E_2^*, I_2^*, T_2^*) \\ &= \left( \frac{p\Lambda}{k+\mu'}, \frac{(k+\mu-\mu p)\Lambda}{(k+\mu)(\beta I_2^* + \mu')}, \frac{(\gamma+\mu+\sigma)(\mu+\delta+\eta) - \alpha\delta\gamma}{\epsilon(\mu+\delta+\eta)} I_2^*, I_2^*, \frac{\gamma}{\mu+\delta+\eta} I_2^* \right), \end{aligned} $$

where:

$$ \begin{aligned} I_2^* &= \frac{(k+\mu-\mu p)\epsilon\Lambda(\mu+\delta+\eta)}{(k+\mu)((\epsilon+\mu)(\gamma+\mu+\sigma)(\mu+\delta+\eta) - (\epsilon+\mu)\alpha\delta\gamma - (1-\alpha)\gamma\delta\epsilon)} - \frac{\mu}{\beta}. \\ &= \frac{\mu}{\beta}(\mathcal{R}_0 - 1). \end{aligned} $$

Then, the endemic equilibrium point $E_2$, exists for $\mathcal{R}_0 > 1$, and at this equilibrium, one has
$$ N = \frac{\Lambda - (\sigma I_2^* + \eta T_2^*)}{\mu} < \frac{\Lambda}{\mu}. $$

### 2.2.3. The Basic Reproduction Number $\mathcal{R}_0$

The basic reproduction number, denoted $\mathcal{R}_0$, is the expected number of secondary cases produced in a completely susceptible population, by a typical infectious individual during its infective period [16].

If $\mathcal{R}_0 < 1$, the disease will not be able to spread among the population, whereas, if $\mathcal{R}_0 > 1$, the disease has the potential to spread among the population and become endemic.

Using the next generation matrix one gets the basic reproduction number $\mathcal{R}_0$ for the proposed model (1) as

$$ \mathcal{R}_0 = \frac{\epsilon(k+\mu-\mu p)\Lambda\beta k_3}{\mu(k+\mu)(k_1k_2k_3 - \alpha\gamma\delta k_1 - (1-\alpha)\delta\gamma\epsilon)}. $$

The details of the calculation of $\mathcal{R}_0$ are presented in the Appendix A and B.

### 2.2.4. Local Stability Analysis of DFE

**Theorem 1.** The disease-free equilibrium point $E_1$ (DFE) of model (1) is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

**Proof.** The Jacobian matrix $J_{E_1}$ of system (1) at the DFE $E_1$ is given by

$$ J_{E_1} = \begin{bmatrix} -(k+\mu) & 0 & 0 & 0 & 0 \\ k & -\mu & 0 & -\beta S_1^* & 0 \\ 0 & 0 & -k_1 & \beta S_1^* & (1-\alpha)\delta \\ 0 & 0 & \epsilon & -k_2 & \alpha\delta \\ 0 & 0 & 0 & \gamma & -k_3 \end{bmatrix} $$

The characteristic equation of $J_{E_1}$ is

$$ (-k+\mu-\lambda)(-\mu-\lambda)[\lambda^3 + a_1\lambda^2 + a_2\lambda + a_3] = 0, \quad (6) $$
where

$$
\begin{aligned}
a_1 &= [k_1 + k_2 + k_3], \\
a_2 &= [k_1 k_2 + k_1 k_3 + k_2 k_3 + \alpha \delta \gamma + \epsilon \beta S_1^*], \\
a_3 &= [-\epsilon \beta S_1^* k_3 + k_1 k_2 k_3 - k_1 \alpha \delta \gamma - \epsilon (1 - \alpha) \delta \gamma] \\
&= [-\epsilon \beta S_1^* k_3 + \frac{\epsilon \beta S_1^* k_3}{\mathcal{R}_0}] \\
&= \epsilon \beta S_1^* k_3 (\frac{1}{\mathcal{R}_0} - 1).
\end{aligned}
$$

Then, all the eigenvalues of the characteristic Equation (6) have a negative real part if the coefficients $a_i$, $i = 1, 2, 3$ fulfill the Routh-Hurwitz conditions, which are $a_1 > 0$, $a_3 > 0$ and $a_1a_2 - a_3 > 0$.

Hence, the disease-free equilibrium of model (1) is locally asymptotically stable, providing that $\mathcal{R}_0 < 1$. $\square$

### 2.2.5. Global Stability Analysis of DFE

**Theorem 2.** The disease-free equilibrium point $E_1$ (DFE) of model (1) is globally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

**Proof.** To prove the theorem, we consider the following Lyapunov function:

$$W(V, S, E, I, T) = b_1E + b_2I + b_3T,$$

where $b_i$, for $i = 1, 2, 3$, are positive constants to be chosen later. Calculating the derivative of $W$ with respect to time along the solutions of system (1), we obtain:

$$
\begin{aligned}
\frac{dW}{dt} &= b_1 \frac{dE}{dt} + b_2 \frac{dI}{dt} + b_3 \frac{dT}{dt} \\
&= b_1[\beta SI - k_1 E + (1-\alpha)\delta T] + b_2[\epsilon E + \alpha \delta T - k_2 I] + b_3[\gamma I - k_3 T] \\
&\le b_1 \left[ \frac{\Lambda \beta}{\mu} I - k_1 E + (1-\alpha)\delta T \right] + b_2[\epsilon E + \alpha \delta T - k_2 I] + b_3[\gamma I - k_3 T], \text{ because } S \le \frac{\Lambda}{\mu} \\
&= \left[ b_1 \frac{\Lambda \beta}{\mu} + b_3 \gamma - b_2 k_2 \right] I + [b_2 \epsilon - b_1 k_1] E + [b_1(1-\alpha)\delta + b_2 \alpha \delta - b_3 k_3] T \\
&\le \frac{(b_2 k_2 - b_3 \gamma)(k+\mu-\mu p)}{(k+\mu)} \left[ \frac{b_1 \Lambda \beta}{\mu(b_2 k_2 - b_3 \gamma)} - 1 \right] I + [b_2 \epsilon - b_1 k_1] E \\
&\quad + [b_1(1-\alpha)\delta + b_2 \alpha \delta - b_3 k_3] T.
\end{aligned}
$$

Choosing $b_1 = \frac{\epsilon k_3}{(k+\mu-\mu p)}(k+\mu)$, $b_2 = \frac{k_1 k_3}{(k+\mu-\mu p)}(k+\mu)$, and $b_3 = \frac{(1-\alpha)\epsilon+\alpha\delta k_1}{(k+\mu-\mu p)}(k+\mu)$, one obtains

$$\frac{dW}{dt} \le \frac{b_1 \Lambda \beta}{\mu \mathcal{R}_0} (\mathcal{R}_0 - 1) I.$$

Hence, if $\mathcal{R}_0 < 1$, then $\frac{dW}{dt}$ is negative. By LaSalle's invariant principle (A1), this implies that $E_1$ is globally asymptotically stable. $\square$
where

$$
\begin{aligned}
a_1 &= [k_1 + k_2 + k_3], \\
a_2 &= [k_1 k_2 + k_1 k_3 + k_2 k_3 + \alpha \delta \gamma + \epsilon \beta S_1^*], \\
a_3 &= [-\epsilon \beta S_1^* k_3 + k_1 k_2 k_3 - k_1 \alpha \delta \gamma - \epsilon (1 - \alpha) \delta \gamma] \\
&= [-\epsilon \beta S_1^* k_3 + \frac{\epsilon \beta S_1^* k_3}{\mathcal{R}_0}] \\
&= \epsilon \beta S_1^* k_3 (\frac{1}{\mathcal{R}_0} - 1).
\end{aligned}
$$

Then, all the eigenvalues of the characteristic Equation (6) have a negative real part if the coefficients $a_i$, $i = 1, 2, 3$ fulfill the Routh-Hurwitz conditions, which are $a_1 > 0$, $a_3 > 0$ and $a_1a_2 - a_3 > 0$.

Hence, the disease-free equilibrium of model (1) is locally asymptotically stable, providing that $\mathcal{R}_0 < 1$. $\square$

### 2.2.5. Global Stability Analysis of DFE

**Theorem 2.** The disease-free equilibrium point $E_1$ (DFE) of model (1) is globally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

**Proof.** To prove the theorem, we consider the following Lyapunov function:

$$W(V, S, E, I, T) = b_1E + b_2I + b_3T,$$

where $b_i$, for $i = 1, 2, 3$, are positive constants to be chosen later. Calculating the derivative of $W$ with respect to time along the solutions of system (1), we obtain:

$$
\begin{aligned}
\frac{dW}{dt} &= b_1 \frac{dE}{dt} + b_2 \frac{dI}{dt} + b_3 \frac{dT}{dt} \\
&= b_1[\beta SI - k_1 E + (1-\alpha)\delta T] + b_2[\epsilon E + \alpha \delta T - k_2 I] + b_3[\gamma I - k_3 T] \\
&\le b_1 \left[ \frac{\Lambda \beta}{\mu} I - k_1 E + (1-\alpha)\delta T \right] + b_2[\epsilon E + \alpha \delta T - k_2 I] + b_3[\gamma I - k_3 T], \text{ because } S \le \frac{\Lambda}{\mu} \\
&= \left[ b_1 \frac{\Lambda \beta}{\mu} + b_3 \gamma - b_2 k_2 \right] I + [b_2 \epsilon - b_1 k_1] E + [b_1(1-\alpha)\delta + b_2 \alpha \delta - b_3 k_3] T \\
&\le \frac{(b_2 k_2 - b_3 \gamma)(k+\mu-\mu p)}{(k+\mu)} \left[ \frac{b_1 \Lambda \beta}{\mu(b_2 k_2 - b_3 \gamma)} - 1 \right] I + [b_2 \epsilon - b_1 k_1] E \\
&\quad + [b_1(1-\alpha)\delta + b_2 \alpha \delta - b_3 k_3] T.
\end{aligned}
$$

Choosing $b_1 = \frac{\epsilon k_3}{(k+\mu-\mu p)}(k+\mu)$, $b_2 = \frac{k_1 k_3}{(k+\mu-\mu p)}(k+\mu)$, and $b_3 = \frac{(1-\alpha)\epsilon+\alpha\delta k_1}{(k+\mu-\mu p)}(k+\mu)$, one obtains

$$\frac{dW}{dt} \le \frac{b_1 \Lambda \beta}{\mu \mathcal{R}_0} (\mathcal{R}_0 - 1) I.$$

Hence, if $\mathcal{R}_0 < 1$, then $\frac{dW}{dt}$ is negative. By LaSalle's invariant principle (A1), this implies that $E_1$ is globally asymptotically stable. $\square$
### 2.2.6. Local Stability Analysis of EE

The Jacobian matrix $J_{E_2}$ of system (1) at the EE $E_2$ is given by

$$J_{E_2} = \begin{bmatrix} -(k+\mu) & 0 & 0 & 0 & 0 \\ k & -\beta I_2^* - \mu & 0 & -\beta S_2^* & 0 \\ 0 & \beta I_2^* & -k_1 & \beta S_2^* & (1-\alpha)\delta \\ 0 & 0 & \epsilon & -k_2 & \alpha\delta \\ 0 & 0 & 0 & \gamma & -k_3 \end{bmatrix}$$

The characteristic equation of $J_{E_2}$ is

$$-(k+\mu)-\lambda)[\lambda^4 + a_1\lambda^3 + a_2\lambda^2 + a_3\lambda + a_4] = 0, \quad (7)$$

where

$$a_1 = (\mu + k_1 + k_2 + k_3 + \beta I^*),$$

$$a_2 = [k_1(\mu + \beta I^*) + k_3(\mu + k_1 + k_2 + \beta I^*) + k_2(\mu + k_1 + \beta I^*) - \gamma\alpha\delta - \epsilon\beta S^*],$$

$$a_3 = [\epsilon(\beta^2 I^* S^* + k_1\beta S^*) + k_3(k_2(\mu + k_1 + \beta I^*) + k_1(\mu + \beta I^*) - \beta\epsilon S^*) + k_1 k_2(\mu + \beta I^*) \\ + \gamma(\alpha\delta k_2 + \epsilon\delta(\alpha - 1)) - \gamma\alpha\delta(\mu + k_1 + k_2 + \beta I^*) - \beta\epsilon S^*(\mu + k_1 + \beta I^*)],$$

$$a_4 = [-\gamma\alpha\delta(k_2(\mu + k_1 + \beta I^*) + k_1(\mu + \beta I^*) - \beta\epsilon S^*) + k_3(k_1 k_2(\mu + \beta I^*) + \epsilon(\beta^2 I^* S^* + k_1\beta S^*) \\ + \beta\epsilon S^*(\mu + k_1 + \beta I^*)) - \gamma(\alpha\delta k_2 - \delta\epsilon(\alpha - 1))(\mu + k_1 + k_2 + \beta I^*) - \gamma(\epsilon(\delta k_1(\alpha - 1) + \beta\alpha\delta S^*)) \\ + k_2(\alpha\delta k_2 + \delta\epsilon(\alpha - 1))].$$

Then, all the eigenvalues of the characteristic Equation (7) have negative real parts if the coefficients $a_i$, $i = 1, 2, 3, 4$ fulfill the Routh-Hurwitz conditions, which are $a_i > 0$ for $i = 1, 2, 3, 4$ and $a_1a_2a_3 > a_3^2 + a_1^2a_4$. Hence, the endemic equilibrium of model (1) is locally asymptotically stable if $R_0 > 1$.

### 2.2.7. Global Stability Analysis of EE

In this section, we prove the global asymptotic stability of the EE of model (1). Using the method described in [17], at the EE, from system (1) we obtain

$$\left\{ \begin{array}{l} p\Lambda = (k + \mu)V_2^*, \\ kV_2^* = -(1-p)\Lambda + \beta S_2^* I_2^* + \mu S_2^*, \\ k_1E_2^* = \beta S_2^* I_2^* + (1-\alpha)\delta T_2^*, \\ k_2I_2^* = \epsilon E_2^* + \alpha\delta T_2^*, \\ \gamma I_2^* = k_3T_2^*. \end{array} \right.$$

**Theorem 3.** The endemic equilibrium point $E_2$ (EE) of model (1) is globally asymptotically stable if $R_0 > 1$.

**Proof.** To prove the theorem, we consider the following Lyapunov function:

$$W(V, S, E, I, T) = k \left[ V(t) - V^* - V^* \ln \frac{V(t)}{V^*} \right] + \epsilon \left[ S(t) - S^* - S^* \ln \frac{S(t)}{S^*} \right] \\ + \epsilon \left[ E(t) - E^* - E^* \ln \frac{E(t)}{E^*} \right] + k_1 \left[ I(t) - I^* - I^* \ln \frac{I(t)}{I^*} \right] \\ + \frac{\delta T^*(k_1\alpha + \epsilon(1-\alpha))}{\gamma I^*} \left[ T(t) - T^* - T^* \ln \frac{T(t)}{T^*} \right].$$

Calculating the derivative of $W$ with respect to time along the solutions of system (1), we obtain:
### 2.2.6. Local Stability Analysis of EE

The Jacobian matrix $J_{E_2}$ of system (1) at the EE $E_2$ is given by

$$J_{E_2} = \begin{bmatrix} -(k+\mu) & 0 & 0 & 0 & 0 \\ k & -\beta I_2^* - \mu & 0 & -\beta S_2^* & 0 \\ 0 & \beta I_2^* & -k_1 & \beta S_2^* & (1-\alpha)\delta \\ 0 & 0 & \epsilon & -k_2 & \alpha\delta \\ 0 & 0 & 0 & \gamma & -k_3 \end{bmatrix}$$

The characteristic equation of $J_{E_2}$ is

$$-(k+\mu)-\lambda)[\lambda^4 + a_1\lambda^3 + a_2\lambda^2 + a_3\lambda + a_4] = 0, \quad (7)$$

where

$$a_1 = (\mu + k_1 + k_2 + k_3 + \beta I^*),$$

$$a_2 = [k_1(\mu + \beta I^*) + k_3(\mu + k_1 + k_2 + \beta I^*) + k_2(\mu + k_1 + \beta I^*) - \gamma\alpha\delta - \epsilon\beta S^*],$$

$$a_3 = [\epsilon(\beta^2 I^* S^* + k_1\beta S^*) + k_3(k_2(\mu + k_1 + \beta I^*) + k_1(\mu + \beta I^*) - \beta\epsilon S^*) + k_1k_2(\mu + \beta I^*) \\ + \gamma(\alpha\delta k_2 + \epsilon\delta(\alpha - 1)) - \gamma\alpha\delta(\mu + k_1 + k_2 + \beta I^*) - \beta\epsilon S^*(\mu + k_1 + \beta I^*)],$$

$$a_4 = [-\gamma\alpha\delta(k_2(\mu + k_1 + \beta I^*) + k_1(\mu + \beta I^*) - \beta\epsilon S^*) + k_3(k_1k_2(\mu + \beta I^*) + \epsilon(\beta^2 I^* S^* + k_1\beta S^*) \\ + \beta\epsilon S^*(\mu + k_1 + \beta I^*)) - \gamma(\alpha\delta k_2 - \delta\epsilon(\alpha - 1))(\mu + k_1 + k_2 + \beta I^*) - \gamma(\epsilon(\delta k_1(\alpha - 1) + \beta\alpha\delta S^*)) \\ + k_2(\alpha\delta k_2 + \delta\epsilon(\alpha - 1))].$$

Then, all the eigenvalues of the characteristic Equation (7) have negative real parts if the coefficients $a_i$, $i = 1, 2, 3, 4$ fulfill the Routh-Hurwitz conditions, which are $a_i > 0$ for $i = 1, 2, 3, 4$ and $a_1a_2a_3 > a_3^2 + a_1^2a_4$. Hence, the endemic equilibrium of model (1) is locally asymptotically stable if $R_0 > 1$.

### 2.2.7. Global Stability Analysis of EE

In this section, we prove the global asymptotic stability of the EE of model (1). Using the method described in [17], at the EE, from system (1) we obtain

$$\begin{cases} p\Lambda = (k + \mu)V_2^*, \\ kV_2^* = -(1-p)\Lambda + \beta S_2^* I_2^* + \mu S_2^*, \\ k_1E_2^* = \beta S_2^* I_2^* + (1-\alpha)\delta T_2^*, \\ k_2I_2^* = \epsilon E_2^* + \alpha\delta T_2^*, \\ \gamma I_2^* = k_3T_2^*. \end{cases}$$

**Theorem 3.** The endemic equilibrium point $E_2$ (EE) of model (1) is globally asymptotically stable if $R_0 > 1$.

**Proof.** To prove the theorem, we consider the following Lyapunov function:

$$W(V, S, E, I, T) = k \left[ V(t) - V^* - V^* \ln \frac{V(t)}{V^*} \right] + \epsilon \left[ S(t) - S^* - S^* \ln \frac{S(t)}{S^*} \right] \\ + \epsilon \left[ E(t) - E^* - E^* \ln \frac{E(t)}{E^*} \right] + k_1 \left[ I(t) - I^* - I^* \ln \frac{I(t)}{I^*} \right] \\ + \frac{\delta T^*(k_1\alpha + \epsilon(1-\alpha))}{\gamma I^*} \left[ T(t) - T^* - T^* \ln \frac{T(t)}{T^*} \right].$$

Calculating the derivative of $W$ with respect to time along the solutions of system (1), we obtain:
$$ \frac{dW}{dt} = k \left[ \left(1 - \frac{V^*}{V}\right) V' \right] + \epsilon \left[ \left(1 - \frac{S^*}{S}\right) S' + \left(1 - \frac{E^*}{E}\right) E' \right] + k_1 \left[ \left(1 - \frac{I^*}{I}\right) I' \right] + \frac{\delta T^*(k_1\alpha + \epsilon(1-\alpha))}{\gamma I^*} \left[ \left(1 - \frac{T^*}{T}\right) T' \right]. $$

A straightforward simple calculation gives

$$ \begin{aligned} k \left(1 - \frac{V^*}{V}\right) V' &= k \left(1 - \frac{V^*}{V}\right) \left[p\Lambda - (k + \mu)V\right] \\ &= k \left(1 - \frac{V^*}{V}\right) \left[(k + \mu)V^* - (k + \mu)V\right] \\ &= k(k + \mu)V^* \left(1 - \frac{V^*}{V}\right) \left(1 - \frac{V}{V^*}\right) \\ &= k(k + \mu)V^* \left(2 - \frac{V^*}{V} - \frac{V}{V^*}\right). \end{aligned} \quad (8) $$

$$ \begin{aligned} \epsilon \left(1 - \frac{S^*}{S}\right) S' &= \epsilon \left(1 - \frac{S^*}{S}\right) \left[(1 - p)\Lambda + kV - \beta SI - \mu S\right] \\ &= \epsilon\mu S^* \left(2 - \frac{S^*}{S} - \frac{S}{S^*}\right) + \epsilon\beta S^* I^* \left(1 - \frac{S^*}{S} + \frac{I}{I^*}\right) \\ &\quad - \epsilon\beta SI. \end{aligned} \quad (9) $$

$$ \begin{aligned} \epsilon \left(1 - \frac{E^*}{E}\right) E' &= \epsilon \left(1 - \frac{E^*}{E}\right) \left[\beta SI - k_1 E + (1 - \alpha)\delta T\right] \\ &= \epsilon\beta SI - \epsilon\beta SI \frac{E^*}{E} - k_1\epsilon E + k_1\epsilon E^* + (1 - \alpha)\epsilon\delta T - (1 - \alpha)\delta\epsilon \frac{TE^*}{E} \quad (10) \\ &= \epsilon\beta SI - \epsilon\beta SI \frac{E^*}{E} - k_1\epsilon E + \epsilon\beta S^* I^* + (1 - \alpha)\delta\epsilon T^* + (1 - \alpha)\epsilon\delta T \\ &\quad - (1 - \alpha)\delta\epsilon \frac{TE^*}{E}. \end{aligned} $$

$$ \begin{aligned} k_1 \left(1 - \frac{I^*}{I}\right) I' &= k_1 \left(1 - \frac{I^*}{I}\right) \left[\epsilon E + \alpha \delta T - k_2 I\right] \\ &= k_1 \epsilon E - k_1 \epsilon E \frac{I^*}{I} + k_1 \alpha \delta T - k_1 \alpha \delta T \frac{I^*}{I} - k_1 k_2 I + k_1 k_2 I^* \\ &= k_1 \epsilon E - k_1 \epsilon E^* \frac{I^* E}{I E^*} + k_1 \alpha \delta T - k_1 \alpha \delta T \frac{I^*}{I} + \epsilon \beta S^* I^* + k_1 \alpha \delta T^* + \epsilon (1 - \alpha) \delta T^* \\ &\quad - \epsilon \beta S^* I^* \frac{I}{I^*} - (1 - \alpha) \epsilon \delta T^* \frac{I}{I^*} - k_1 \alpha \delta T^* \frac{I}{I^*} \\ &= k_1 \epsilon E - \epsilon \beta S^* I^* \frac{E I^*}{E^* I} - (1 - \alpha) \epsilon \delta T^* \frac{E I^*}{E^* I} + k_1 \alpha \delta T - k_1 \alpha \delta T \frac{I^*}{I} + \epsilon \beta S^* I^* + k_1 \alpha \delta T^* \\ &\quad + \epsilon (1 - \alpha) \delta T^* - \epsilon \beta S^* I^* \frac{I}{I^*} - (1 - \alpha) \epsilon \delta T^* \frac{I}{I^*} - k_1 \alpha \delta T^* \frac{I}{I^*}. \end{aligned} \quad (11) $$

$$ \begin{aligned} \frac{\delta T^*(k_1\alpha + \epsilon(1-\alpha))}{\gamma I^*} \left(1 - \frac{T^*}{T}\right) \left[\gamma I - k_3 T\right] &= \frac{\delta T^*(k_1\alpha + \epsilon(1-\alpha))}{I^*} \left(1 - \frac{T^*}{T}\right) \left[I - \frac{I^*}{T^*} T\right] \\ &= \delta\alpha k_1 T^* \frac{I}{I^*} + \delta(1-\alpha)\epsilon T^* \frac{I}{I^*} - \delta\alpha k_1 T^* \frac{I^*}{I^* T} \\ &\quad - \delta(1-\alpha)\epsilon T^* \frac{I^*}{I^* T} - \delta\alpha k_1 T - \delta(1-\alpha)\epsilon T \\ &\quad + \delta\alpha k_1 T^* + \delta(1-\alpha)\epsilon T^*. \end{aligned} \quad (12) $$
$$ \frac{dW}{dt} = k \left[ \left(1 - \frac{V^*}{V}\right) V' \right] + \epsilon \left[ \left(1 - \frac{S^*}{S}\right) S' + \left(1 - \frac{E^*}{E}\right) E' \right] + k_1 \left[ \left(1 - \frac{I^*}{I}\right) I' \right] \\ + \frac{\delta T^*(k_1\alpha + \epsilon(1-\alpha))}{\gamma I^*} \left[ \left(1 - \frac{T^*}{T}\right) T' \right]. $$

A straightforward simple calculation gives

$$ \begin{aligned} k \left(1 - \frac{V^*}{V}\right) V' &= k \left(1 - \frac{V^*}{V}\right) \left[p\Lambda - (k + \mu)V\right] \\ &= k \left(1 - \frac{V^*}{V}\right) \left[(k + \mu)V^* - (k + \mu)V\right] \\ &= k(k + \mu)V^* \left(1 - \frac{V^*}{V}\right) \left(1 - \frac{V}{V^*}\right) \\ &= k(k + \mu)V^* \left(2 - \frac{V^*}{V} - \frac{V}{V^*}\right). \end{aligned} \tag{8} $$

$$ \begin{aligned} \epsilon \left(1 - \frac{S^*}{S}\right) S' &= \epsilon \left(1 - \frac{S^*}{S}\right) \left[(1 - p)\Lambda + kV - \beta SI - \mu S\right] \\ &= \epsilon\mu S^* \left(2 - \frac{S^*}{S} - \frac{S}{S^*}\right) + \epsilon\beta S^* I^* \left(1 - \frac{S^*}{S} + \frac{I}{I^*}\right) \\ &\quad - \epsilon\beta SI. \end{aligned} \tag{9} $$

$$ \begin{aligned} \epsilon \left(1 - \frac{E^*}{E}\right) E' &= \epsilon \left(1 - \frac{E^*}{E}\right) \left[\beta SI - k_1 E + (1 - \alpha)\delta T\right] \\ &= \epsilon\beta SI - \epsilon\beta SI \frac{E^*}{E} - k_1\epsilon E + k_1\epsilon E^* + (1 - \alpha)\epsilon\delta T - (1 - \alpha)\delta\epsilon \frac{TE^*}{E} \tag{10} \\ &= \epsilon\beta SI - \epsilon\beta SI \frac{E^*}{E} - k_1\epsilon E + \epsilon\beta S^* I^* + (1 - \alpha)\delta\epsilon T^* + (1 - \alpha)\epsilon\delta T \\ &\quad - (1 - \alpha)\delta\epsilon \frac{TE^*}{E}. \end{aligned} $$

$$ \begin{aligned} k_1 \left(1 - \frac{I^*}{I}\right) I' &= k_1 \left(1 - \frac{I^*}{I}\right) \left[\epsilon E + \alpha \delta T - k_2 I\right] \\ &= k_1 \epsilon E - k_1 \epsilon E \frac{I^*}{I} + k_1 \alpha \delta T - k_1 \alpha \delta T \frac{I^*}{I} - k_1 k_2 I + k_1 k_2 I^* \\ &= k_1 \epsilon E - k_1 \epsilon E^* \frac{I^* E}{I E^*} + k_1 \alpha \delta T - k_1 \alpha \delta T \frac{I^*}{I} + \epsilon \beta S^* I^* + k_1 \alpha \delta T^* + \epsilon (1 - \alpha) \delta T^* \\ &\quad - \epsilon \beta S^* I^* \frac{I}{I^*} - (1 - \alpha) \epsilon \delta T^* \frac{I}{I^*} - k_1 \alpha \delta T^* \frac{I}{I^*} \\ &= k_1 \epsilon E - \epsilon \beta S^* I^* \frac{E I^*}{E^* I} - (1 - \alpha) \epsilon \delta T^* \frac{E I^*}{E^* I} + k_1 \alpha \delta T - k_1 \alpha \delta T \frac{I^*}{I} + \epsilon \beta S^* I^* + k_1 \alpha \delta T^* \\ &\quad + \epsilon (1 - \alpha) \delta T^* - \epsilon \beta S^* I^* \frac{I}{I^*} - (1 - \alpha) \epsilon \delta T^* \frac{I}{I^*} - k_1 \alpha \delta T^* \frac{I}{I^*}. \end{aligned} \tag{11} $$

$$ \begin{aligned} \frac{\delta T^*(k_1\alpha + \epsilon(1-\alpha))}{\gamma I^*} \left(1 - \frac{T^*}{T}\right) \left[\gamma I - k_3 T\right] &= \frac{\delta T^*(k_1\alpha + \epsilon(1-\alpha))}{I^*} \left(1 - \frac{T^*}{T}\right) \left[I - \frac{I^*}{T^*} T\right] \\ &= \delta\alpha k_1 T^* \frac{I}{I^*} + \delta(1-\alpha)\epsilon T^* \frac{I}{I^*} - \delta\alpha k_1 T^* \frac{I^*}{I^* T} \\ &\quad - \delta(1-\alpha)\epsilon T^* \frac{I^*}{I^* T} - \delta\alpha k_1 T - \delta(1-\alpha)\epsilon T \\ &\quad + \delta\alpha k_1 T^* + \delta(1-\alpha)\epsilon T^*. \end{aligned} \tag{12} $$
Using Equations (8)–(12), one gets

$$
\begin{equation}
\begin{aligned}
\frac{dW}{dt} = & k(k + \mu)V^*\left(2 - \frac{V^*}{V} - \frac{V}{V^*}\right) + \epsilon\mu S^*\left(2 - \frac{S^*}{S} - \frac{S}{S^*}\right) \\
& + \epsilon\beta S^* I^*\left(3 - \frac{S^*}{S} - \frac{I}{I^*} - \frac{EI^*}{E^*I} - \frac{SIE^*}{S^*I^*E}\left(1 - \frac{ES^*}{E^*S}\right)\right) \\
& + (1 - \alpha)\epsilon\delta T^*\left(3 - \frac{I^*E}{IE^*} - \frac{E^*T}{ET^*} - \frac{T^*I}{TI^*}\right) + \delta\alpha k_1 T^*\left(2 - \frac{I^*T}{IT^*} - \frac{IT^*}{I^*T}\right).
\end{aligned}
\tag{13}
\end{equation}
$$

Using the properties of the geometric and arithmetic means in Equation (13), one obtains

$$
\left\{
\begin{array}{l}
2 - \frac{V}{V^*} - \frac{V^*}{V} \le 0, \\
2 - \frac{S}{S^*} - \frac{S^*}{S} \le 0, \\
3 - \frac{S^*}{S} - \frac{I}{I^*} - \frac{EI^*}{E^*I} - \frac{SIE^*}{S^*I^*E} \left(1 - \frac{ES^*}{E^*S}\right) \le 0, \\
3 - \frac{I^*E}{IE^*} - \frac{E^*T}{FT^*} - \frac{T^*I}{TI^*} \le 0, \\
2 - \frac{I^*T}{IT^*} - \frac{I^*T}{I^*T} \le 0.
\end{array}
\right.
$$

Since none of the parameters are negative, it follows that $\frac{dW}{dt} \le 0$ when $R_0 > 1$. As a result, according to LaSalle's Invariance Principle (A1), $(V, S, E, I, T) \to (V^*, S^*, E^*, I^*, T^*)$ as $t \to \infty$. $\square$

2.2.8. Transcritical Bifurcation Analysis

Here, we discuss the existence of transcritical bifurcation of system (1). At $R_0 = 1$, if we take $\beta$ as a bifurcation parameter, we obtain

$$
\beta^* = \beta = \frac{\mu(k + \mu)(k_1 k_2 k_3 - \alpha \gamma \delta k_1 - (1 - \alpha) \delta \gamma \epsilon)}{\epsilon(k + \mu - \mu p) \Lambda k_3}.
$$

The following modification are made in the variables ofsystem (1) so that $V = x_1, S = x_2, L = x_3, I = x_4$, and $T = x_5$. Further using vector notation $x = (x_1, x_2, x_3, x_4, x_5)^T$, model (1)
can then be written in the form $\frac{dx}{dt} = F$, with $F = (f_1, f_2, f_3, f_4, f_5)^T$ as shown below

$$
\left\{
\begin{aligned}
\dot{x}_1 &= p\Lambda - (k + \mu)x_1, \\
\dot{x}_2 &= (1 - p)\Lambda + kx_1 - \beta x_2 x_4 - \mu x_2, \\
\dot{x}_3 &= \beta x_2 x_4 - (\epsilon + \mu)x_3 + (1 - \alpha)\delta x_5, \\
\dot{x}_4 &= \epsilon x_3 + \alpha\delta x_5 - (\gamma + \mu + \sigma)x_4, \\
\dot{x}_5 &= \gamma x_4 - (\mu + \delta + \eta)x_5.
\end{aligned}
\right.
\tag{14}
$$

with $N = \sum_{i=1}^{5} x_i$.

The Jacobian matrix evaluated at the disease-free equilibrium $E_1$ (DFE) for $\beta = \beta^*$ is

$$
J_{E_1} = \begin{bmatrix}
-(k + \mu) & 0 & 0 & 0 & 0 \\
k & -\mu & 0 & -\beta^* S^* & 0 \\
0 & 0 & -k_1 & \beta^* S^* & (1 - \alpha)\delta \\
0 & 0 & \epsilon & -k_2 & \alpha\delta \\
0 & 0 & 0 & \gamma & -k_3
\end{bmatrix}
$$

The Jacobian matrix $J_{E_1}$ has a simple zero eigenvalue calculated at $\beta^*$.
Using Equations (8)–(12), one gets

$$ \frac{dW}{dt} = k(k + \mu)V^*\left(2 - \frac{V^*}{V} - \frac{V}{V^*}\right) + \epsilon\mu S^*\left(2 - \frac{S^*}{S} - \frac{S}{S^*}\right) \\ + \epsilon\beta S^* I^*\left(3 - \frac{S^*}{S} - \frac{I}{I^*} - \frac{EI^*}{E^*I} - \frac{SIE^*}{S^*I^*E}\left(1 - \frac{ES^*}{E^*S}\right)\right) \\ + (1 - \alpha)\epsilon\delta T^*\left(3 - \frac{I^*E}{IE^*} - \frac{E^*T}{ET^*} - \frac{T^*I}{TI^*}\right) + \delta\alpha k_1 T^*\left(2 - \frac{I^*T}{IT^*} - \frac{IT^*}{I^*T}\right). \quad (13) $$

Using the properties of the geometric and arithmetic means in Equation (13), one obtains

$$ \left\{ \begin{array}{l} 2 - \frac{V}{V^*} - \frac{V^*}{V} \le 0, \\ 2 - \frac{S}{S^*} - \frac{S^*}{S} \le 0, \\ 3 - \frac{S^*}{S} - \frac{I}{I^*} - \frac{EI^*}{E^*I} - \frac{SIE^*}{S^*I^*E}\left(1 - \frac{ES^*}{E^*S}\right) \le 0, \\ 3 - \frac{I^*E}{IE^*} - \frac{E^*T}{FT^*} - \frac{T^*I}{TI^*} \le 0, \\ 2 - \frac{I^*T}{IT^*} - \frac{I^*T}{I^*T} \le 0. \end{array} \right. $$

Since none of the parameters are negative, it follows that $\frac{dW}{dt} \le 0$ when $R_0 > 1$. As a result, according to LaSalle's Invariance Principle (A1), $(V, S, E, I, T) \to (V^*, S^*, E^*, I^*, T^*)$ as $t \to \infty$. $\square$

### 2.2.8. Transcritical Bifurcation Analysis

Here, we discuss the existence of transcritical bifurcation of system (1). At $R_0 = 1$, if we take $\beta$ as a bifurcation parameter, we obtain

$$ \beta^* = \beta = \frac{\mu(k + \mu)(k_1 k_2 k_3 - \alpha \gamma \delta k_1 - (1 - \alpha) \delta \gamma \epsilon)}{\epsilon(k + \mu - \mu p) \Lambda k_3}. $$

The following modification are made in the variables ofsystem (1) so that $V = x_1$, $S = x_2$, $L = x_3$, $I = x_4$, and $T = x_5$. Further using vector notation $x = (x_1, x_2, x_3, x_4, x_5)^T$, model (1) can then be written in the form $\frac{dx}{dt} = F$, with $F = (f_1, f_2, f_3, f_4, f_5)^T$ as shown below

$$ \left\{ \begin{aligned} \dot{x}_1 &= p\Lambda - (k + \mu)x_1, \\ \dot{x}_2 &= (1 - p)\Lambda + kx_1 - \beta x_2 x_4 - \mu x_2, \\ \dot{x}_3 &= \beta x_2 x_4 - (\epsilon + \mu)x_3 + (1 - \alpha)\delta x_5, \\ \dot{x}_4 &= \epsilon x_3 + \alpha\delta x_5 - (\gamma + \mu + \sigma)x_4, \\ \dot{x}_5 &= \gamma x_4 - (\mu + \delta + \eta)x_5. \end{aligned} \right. \quad (14) $$

with $N = \sum_{i=1}^{5} x_i$.

The Jacobian matrix evaluated at the disease-free equilibrium $E_1$ (DFE) for $\beta = \beta^*$ is

$$ J_{E_1} = \begin{bmatrix} -(k + \mu) & 0 & 0 & 0 & 0 \\ k & -\mu & 0 & -\beta^* S^* & 0 \\ 0 & 0 & -k_1 & \beta^* S^* & (1 - \alpha)\delta \\ 0 & 0 & \epsilon & -k_2 & \alpha\delta \\ 0 & 0 & 0 & \gamma & -k_3 \end{bmatrix} $$

The Jacobian matrix $J_{E_1}$ has a simple zero eigenvalue calculated at $\beta^*$.
The right and left eigenvectors denoted by: $Y = (y_1, y_2, y_3, y_4, y_5)$ and $Z = (z_1, z_2, z_3, z_4, z_5)$, respectively, are obtained as follows

$$y_1 = 0, y_2 = \frac{(k_1 k_2 k_3 - \alpha \gamma \delta k_1 - (1 - \alpha) \delta \gamma \epsilon)}{\mu \epsilon k_3}, y_3 = \frac{k_1 k_2 k_3 - \alpha \gamma \delta k_1}{\epsilon k_3}, y_4, y_4 > 0, y_5 = \frac{\gamma}{k_3} y_4.$$

and

$$z_1 = 0, z_2 = 0, z_3 > 0, z_4 = \frac{k_1}{\epsilon} z_3, z_5 = \frac{(1 - \alpha) \delta \epsilon + \alpha \delta k_1}{\epsilon k_3} z_3.$$

We have

$$\left\{ \begin{array}{l} Y^t D_\beta F(E_1, \beta^*) = 0, \\ Y^t D_x D_\beta F(E_1, \beta^*) Z = -S^* y_2 z_4 + S^* y_3 z_4, \\ Y^t D_x^2 F(E_1, \beta^*) (Z, Z) = (-\beta^* y_2 z_4, -\beta^* y_2 z_4). \end{array} \right. \qquad (15)$$

The following conditions are satisfied:

$$\left\{ \begin{array}{l} Y^t D_\beta F(E_1, \beta^*) = 0, \\ Y^t D_x D_\beta F(E_1, \beta^*) Z \neq 0, \\ Y^t D_x^2 F(E_1, \beta^*) (Z, Z) \neq 0. \end{array} \right. \qquad (16)$$

Hence, there is a transcritical bifurcation at $\beta = \beta^*$ as illustrated in Figure 2.

<figure><img src="image_4.png" /><figcaption>Figure 2. Transcritical bifurcation of the model (1).</figcaption></figure>

### 3. Parameters Estimation and Numerical Simulation

In this section, six model parameters will be estimated based on TB incidence data from the WHO Global TB Report [1] between 1990 and 2020 (see Table 3) and the other parameters will be inspired by the statistical data in the literature.

#### 3.1. Parameters Estimation

Using the data of Algeria's population from [18], one takes the death rate $\mu$ as the average death rate per year from 1990 to 2020, $\mu = 0.00498$, and the recruitment rate $\Lambda$, as the average birth per year from 1990 to 2020, $\Lambda = 811,085$.
The right and left eigenvectors denoted by: $Y = (y_1, y_2, y_3, y_4, y_5)$ and $Z = (z_1, z_2, z_3, z_4, z_5)$, respectively, are obtained as follows

$$y_1 = 0, y_2 = \frac{(k_1 k_2 k_3 - \alpha \gamma \delta k_1 - (1 - \alpha) \delta \gamma \epsilon)}{\mu \epsilon k_3}, y_3 = \frac{k_1 k_2 k_3 - \alpha \gamma \delta k_1}{\epsilon k_3}, y_4, y_4 > 0, y_5 = \frac{\gamma}{k_3} y_4.$$

and

$$z_1 = 0, z_2 = 0, z_3 > 0, z_4 = \frac{k_1}{\epsilon} z_3, z_5 = \frac{(1 - \alpha) \delta \epsilon + \alpha \delta k_1}{\epsilon k_3} z_3.$$

We have

$$\left\{ \begin{array}{l} Y^t D_\beta F(E_1, \beta^*) = 0, \\ Y^t D_x D_\beta F(E_1, \beta^*) Z = -S^* y_2 z_4 + S^* y_3 z_4, \\ Y^t D_x^2 F(E_1, \beta^*) (Z, Z) = (-\beta^* y_2 z_4, -\beta^* y_2 z_4). \end{array} \right. \qquad (15)$$

The following conditions are satisfied:

$$\left\{ \begin{array}{l} Y^t D_\beta F(E_1, \beta^*) = 0, \\ Y^t D_x D_\beta F(E_1, \beta^*) Z \neq 0, \\ Y^t D_x^2 F(E_1, \beta^*) (Z, Z) \neq 0. \end{array} \right. \qquad (16)$$

Hence, there is a transcritical bifurcation at $\beta = \beta^*$ as illustrated in Figure 2

<figure><img src="image_7.png" /><figcaption>Figure 2. Transcritical bifurcation of the model (1).</figcaption></figure>

### 3. Parameters Estimation and Numerical Simulation

In this section, six model parameters will be estimated based on TB incidence data from the WHO Global TB Report [1] between 1990 and 2020 (see Table 3) and the other parameters will be inspired by the statistical data in the literature.

#### 3.1. Parameters Estimation

Using the data of Algeria's population from [18], one takes the death rate $\mu$ as the average death rate per year from 1990 to 2020, $\mu = 0.00498$, and the recruitment rate $\Lambda$, as the average birth per year from 1990 to 2020, $\Lambda = 811, 085$.
The child immunization rate, BCG, is the ratio of children aged 12–23 months who have received BCG vaccination. Figure 3 displays the percentage of one-year-old children who have received the BCG immunization in Algeria between 1990 and 2020, according to data from officially recognized sources compiled by the World Bank [19]. Hence, one gets the average vaccination rate $p = 0.977$. The BCG has shown an overall efficacy of 70% to 80% against childhood TB, namely meningitis and miliary TB [20]. Hence, in this paper one takes the average rate of moving from $V$ to $S$ as the BCG immunization failure $k = 1 - 0.75 = 0.25$. The treatment success from 2000 to 2020 [21] is used to calculate the treatment failure rate $\alpha = 1 - 0.8905 = 0.1095$.

<figure><img src="image_5.png" /><figcaption>Figure 3. Percentage of one-year-old children in Algeria receiving BCG vaccination during the time period 1990–2020.</figcaption></figure>

The initial conditions were carefully selected as follows: The total initial population, $N(0)$, was set to 25,518,074, which corresponds to the population of Algeria in 1990, as reported in [18]. The initial infected population, $I(0) = 11,607$, was obtained from the WHO Global TB Report [1]. The initial exposed population, $E(0)$, was assumed to be 8852. Additionally, the initial treated population, $T(0)$, was set to 20,000, whereas the number of vaccinated individuals, $V(0)$, was determined to be 8,109,389. As a result of these values, the initial susceptible population can be calculated as $S(0) = N(0) - E(0) - I(0) - T(0) - V(0) = 17,368,226$. These initial conditions were carefully chosen to ensure accurate and reliable numerical simulations of the system under study.

One estimates the parameters $\beta, \gamma, \epsilon, \sigma, \alpha, \delta, k$, and $\eta$ by minimizing the error between actual TB incidence data and the solution of the proposed model (1). The objective function used in this parameter estimation is given by

$$ \psi = \sum_{i=1}^{n} (I_{t_i} - I_{t_i}^*)^2, \qquad (17) $$

where $I_{t_i}^*$ denotes the actual TB-infected case, $I_{t_i}$ is the corresponding model solution at time $t_i$, and $n$ is the number of available actual data. The MATLAB function 'fitnlm', which solves nonlinear regression problems based on the Levenberg–Marquardt algorithm in MATLAB R2020b, was employed to minimize the function (17).

In Figure 4 and Table 4, the incidence data are shown along with the model-fitted values, obtained using the values in Table 3.
The child immunization rate, BCG, is the ratio of children aged 12–23 months who have received BCG vaccination. Figure 3 displays the percentage of one-year-old children who have received the BCG immunization in Algeria between 1990 and 2020, according to data from officially recognized sources compiled by the World Bank [19]. Hence, one gets the average vaccination rate $p = 0.977$. The BCG has shown an overall efficacy of 70% to 80% against childhood TB, namely meningitis and miliary TB [20]. Hence, in this paper one takes the average rate of moving from $V$ to $S$ as the BCG immunization failure $k = 1 - 0.75 = 0.25$. The treatment success from 2000 to 2020 [21] is used to calculate the treatment failure rate $\alpha = 1 - 0.8905 = 0.1095$.

<figure><img src="image_8.png" /><figcaption>Figure 3. Percentage of one-year-old children in Algeria receiving BCG vaccination during the time period 1990–2020.</figcaption></figure>

The initial conditions were carefully selected as follows: The total initial population, $N(0)$, was set to 25,518,074, which corresponds to the population of Algeria in 1990, as reported in [18]. The initial infected population, $I(0) = 11,607$, was obtained from the WHO Global TB Report [1]. The initial exposed population, $E(0)$, was assumed to be 8852. Additionally, the initial treated population, $T(0)$, was set to 20,000, whereas the number of vaccinated individuals, $V(0)$, was determined to be 8,109,389. As a result of these values, the initial susceptible population can be calculated as $S(0) = N(0) - E(0) - I(0) - T(0) - V(0) = 17,368,226$. These initial conditions were carefully chosen to ensure accurate and reliable numerical simulations of the system under study.

One estimates the parameters $\beta, \gamma, \epsilon, \sigma, \alpha, \delta, k$, and $\eta$ by minimizing the error between actual TB incidence data and the solution of the proposed model (1). The objective function used in this parameter estimation is given by

$$ \psi = \sum_{i=1}^{n} (I_{t_i} - I_{t_i}^*)^2, \qquad (17) $$

where $I_{t_i}^*$ denotes the actual TB-infected case, $I_{t_i}$ is the corresponding model solution at time $t_i$, and $n$ is the number of available actual data. The MATLAB function 'fitnlm', which solves nonlinear regression problems based on the Levenberg–Marquardt algorithm in MATLAB R2020b, was employed to minimize the function (17).

In Figure 4 and Table 4, the incidence data are shown along with the model-fitted values, obtained using the values in Table 3.
<table><caption>**Table 3.** Parameters and initial data of the model (1).</caption><thead><tr><th>Parameters</th><th>Description</th><th>Algeria's Parameters</th><th>References</th></tr></thead><tbody><tr><td>V(0)</td><td>Initial number of vaccinated</td><td>8,109,389</td><td>Assumed</td></tr><tr><td>S(0)</td><td>Initial number of susceptible</td><td>17,368,226</td><td>Calculated</td></tr><tr><td>E(0)</td><td>Initial number of exposed</td><td>8852</td><td>Assumed</td></tr><tr><td>I(0)</td><td>Initial number of infected</td><td>11,607</td><td>[1]</td></tr><tr><td>T(0)</td><td>Initial number of treated</td><td>20,000</td><td>Assumed</td></tr><tr><td>Λ</td><td>Recruitment rate</td><td>811,085</td><td>[18]</td></tr><tr><td>μ</td><td>Natural death rate</td><td>0.00498</td><td>[18]</td></tr><tr><td>k</td><td>Rate of moving from V to S</td><td>0.25</td><td>[19]</td></tr><tr><td>β</td><td>Transmission rate</td><td>6.6752 × 10<sup>-11</sup></td><td>Fitted</td></tr><tr><td>γ</td><td>Treatment rate</td><td>0.0043</td><td>Fitted</td></tr><tr><td>ε</td><td>Progression rate</td><td>0.0656</td><td>Fitted</td></tr><tr><td>α</td><td>Treatment failure rate</td><td>0.1095</td><td>[21]</td></tr><tr><td>δ</td><td>Rate at which the treated population leaves the class T</td><td>0.1325</td><td>Fitted</td></tr><tr><td>σ</td><td>Disease death rate in I</td><td>0.0136</td><td>Fitted</td></tr><tr><td>η</td><td>Disease death rate in T</td><td>4.2327 × 10<sup>-6</sup></td><td>Fitted</td></tr><tr><td>p</td><td>Vaccination rate</td><td>0.977</td><td>[19]</td></tr></tbody></table>

<table><caption>**Table 4.** The reported data and the model fitted values of TB cases in Algeria.</caption><thead><tr><th>Year</th><th>Reported Data</th><th>Numerical Value</th><th>Year</th><th>Reported Data</th><th>Numerical Value</th></tr></thead><tbody><tr><td>1990</td><td>11,607</td><td>11,607</td><td>2006</td><td>21,143</td><td>20,613</td></tr><tr><td>1991</td><td>11,332</td><td>12,162</td><td>2007</td><td>21,369</td><td>20,884</td></tr><tr><td>1992</td><td>11,428</td><td>12,793</td><td>2008</td><td>20,588</td><td>21,116</td></tr><tr><td>1993</td><td>13,345</td><td>13,471</td><td>2009</td><td>21,701</td><td>21,313</td></tr><tr><td>1994</td><td>13,345</td><td>14,137</td><td>2010</td><td>22,336</td><td>21,474</td></tr><tr><td>1995</td><td>13,507</td><td>14,880</td><td>2011</td><td>21,429</td><td>21,604</td></tr><tr><td>1996</td><td>15,329</td><td>15,578</td><td>2012</td><td>21,880</td><td>21,705</td></tr><tr><td>1997</td><td>16,522</td><td>16,255</td><td>2013</td><td>20,701</td><td>21,778</td></tr><tr><td>1998</td><td>15,324</td><td>16,255</td><td>2014</td><td>22,517</td><td>21,825</td></tr><tr><td>1999</td><td>16,647</td><td>17,517</td><td>2015</td><td>23,705</td><td>21,850</td></tr><tr><td>2000</td><td>18,572</td><td>18,090</td><td>2016</td><td>22,801</td><td>21,854</td></tr><tr><td>2001</td><td>18,250</td><td>18,621</td><td>2017</td><td>23,077</td><td>21,838</td></tr><tr><td>2002</td><td>18,934</td><td>19,108</td><td>2018</td><td>23,465</td><td>21,805</td></tr><tr><td>2003</td><td>19,730</td><td>19,550</td><td>2019</td><td>20,879</td><td>21,757</td></tr><tr><td>2004</td><td>19,929</td><td>19,947</td><td>2020</td><td>17,212</td><td>21,694</td></tr><tr><td>2005</td><td>21,336</td><td>20,301</td><td>2021</td><td>-</td><td>21,619</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Data fitting of the number of TB cases in Algeria.</figcaption></figure>
<table><caption>**Table 3.** Parameters and initial data of the model (1).</caption><thead><tr><th>Parameters</th><th>Description</th><th>Algeria's Parameters</th><th>References</th></tr></thead><tbody><tr><td>V(0)</td><td>Initial number of vaccinated</td><td>8,109,389</td><td>Assumed</td></tr><tr><td>S(0)</td><td>Initial number of susceptible</td><td>17,368,226</td><td>Calculated</td></tr><tr><td>E(0)</td><td>Initial number of exposed</td><td>8852</td><td>Assumed</td></tr><tr><td>I(0)</td><td>Initial number of infected</td><td>11,607</td><td>[1]</td></tr><tr><td>T(0)</td><td>Initial number of treated</td><td>20,000</td><td>Assumed</td></tr><tr><td>Λ</td><td>Recruitment rate</td><td>811,085</td><td>[18]</td></tr><tr><td>μ</td><td>Natural death rate</td><td>0.00498</td><td>[18]</td></tr><tr><td>k</td><td>Rate of moving from V to S</td><td>0.25</td><td>[19]</td></tr><tr><td>β</td><td>Transmission rate</td><td>6.6752 × 10<sup>-11</sup></td><td>Fitted</td></tr><tr><td>γ</td><td>Treatment rate</td><td>0.0043</td><td>Fitted</td></tr><tr><td>ε</td><td>Progression rate</td><td>0.0656</td><td>Fitted</td></tr><tr><td>α</td><td>Treatment failure rate</td><td>0.1095</td><td>[21]</td></tr><tr><td>δ</td><td>Rate at which the treated population leaves the class T</td><td>0.1325</td><td>Fitted</td></tr><tr><td>σ</td><td>Disease death rate in I</td><td>0.0136</td><td>Fitted</td></tr><tr><td>η</td><td>Disease death rate in T</td><td>4.2327 × 10<sup>-6</sup></td><td>Fitted</td></tr><tr><td>p</td><td>Vaccination rate</td><td>0.977</td><td>[19]</td></tr></tbody></table>

<table><caption>**Table 4.** The reported data and the model fitted values of TB cases in Algeria.</caption><thead><tr><th>Year</th><th>Reported Data</th><th>Numerical Value</th><th>Year</th><th>Reported Data</th><th>Numerical Value</th></tr></thead><tbody><tr><td>1990</td><td>11,607</td><td>11,607</td><td>2006</td><td>21,143</td><td>20,613</td></tr><tr><td>1991</td><td>11,332</td><td>12,162</td><td>2007</td><td>21,369</td><td>20,884</td></tr><tr><td>1992</td><td>11,428</td><td>12,793</td><td>2008</td><td>20,588</td><td>21,116</td></tr><tr><td>1993</td><td>13,345</td><td>13,471</td><td>2009</td><td>21,701</td><td>21,313</td></tr><tr><td>1994</td><td>13,345</td><td>14,137</td><td>2010</td><td>22,336</td><td>21,474</td></tr><tr><td>1995</td><td>13,507</td><td>14,880</td><td>2011</td><td>21,429</td><td>21,604</td></tr><tr><td>1996</td><td>15,329</td><td>15,578</td><td>2012</td><td>21,880</td><td>21,705</td></tr><tr><td>1997</td><td>16,522</td><td>16,255</td><td>2013</td><td>20,701</td><td>21,778</td></tr><tr><td>1998</td><td>15,324</td><td>16,255</td><td>2014</td><td>22,517</td><td>21,825</td></tr><tr><td>1999</td><td>16,647</td><td>17,517</td><td>2015</td><td>23,705</td><td>21,850</td></tr><tr><td>2000</td><td>18,572</td><td>18,090</td><td>2016</td><td>22,801</td><td>21,854</td></tr><tr><td>2001</td><td>18,250</td><td>18,621</td><td>2017</td><td>23,077</td><td>21,838</td></tr><tr><td>2002</td><td>18,934</td><td>19,108</td><td>2018</td><td>23,465</td><td>21,805</td></tr><tr><td>2003</td><td>19,730</td><td>19,550</td><td>2019</td><td>20,879</td><td>21,757</td></tr><tr><td>2004</td><td>19,929</td><td>19,947</td><td>2020</td><td>17,212</td><td>21,694</td></tr><tr><td>2005</td><td>21,336</td><td>20,301</td><td>2021</td><td>-</td><td>21,619</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Data fitting of the number of TB cases in Algeria.</figcaption></figure>
## 3.2. Sensitivity Analysis

By using sensitivity, the spread and prevalence of diseases can be analyzed for each parameter. As a result of errors in data collection and assumed parameters, it is commonly used to measure the robustness of model predictions. To determine the relative significance of these parameters on disease transmission, we examined the impact of various model parameters. It has been determined how the model parameters $\beta$, $\gamma$, $k$, and $\epsilon$ affect the partial derivatives of the basic reproduction number $R_0$. Since $\frac{\partial R_0}{\partial \beta} > 0$, it follows that the transmission rate can be lowered to lessen the infection. However, given the partial derivatives $\frac{\partial R_0}{\partial \gamma} < 0$, it is implied that TB infection can be controlled by raising the parameter $\gamma$. Increasing $k$ results in an increase in $R_0$ because $\frac{\partial R_0}{\partial k} > 0$. It demonstrates that the diseased population will grow more quickly. Parameter $\alpha$ is the failure rate of treatment and $\frac{\partial R_0}{\partial \alpha} > 0$. Therefore, by lowering the treatment failure rate $\alpha$, the cumulative number of infected people can be minimized.

Sensitivity indices should be computed to estimate the relative change in a variable when parameters change. These indications were calculated using the following definition.

**Definition 1.** For a certain value $\sigma$, the normalized forward sensitivity index of $R_0$ is determined by

$$S_{\sigma}^{R_0} = \frac{\sigma}{R_0} \frac{\partial R_0}{\partial \sigma} \qquad (18)$$

For the baseline model parameters are determined by Formula (18), the derived sensitivity indices of the basic reproduction number $R_0$ are shown in Table 5.

<table><caption>**Table 5.** Parameters and sensitivity index.</caption><thead><tr><th>Parameter</th><th>Sensitivity Index</th></tr></thead><tbody><tr><td>Λ</td><td>+1</td></tr><tr><td>μ</td><td>-1.6502</td></tr><tr><td>k</td><td>+0.0012</td></tr><tr><td>β</td><td>+1</td></tr><tr><td>γ</td><td>-0.1671</td></tr><tr><td>ε</td><td>+0.0005</td></tr><tr><td>α</td><td>+2.1364 × 10<sup>-10</sup></td></tr><tr><td>δ</td><td>+1.4003 × 10<sup>-09</sup></td></tr><tr><td>σ</td><td>-0.4043</td></tr><tr><td>η</td><td>-2.5311 × 10<sup>-11</sup></td></tr><tr><td>p</td><td>-0.0194</td></tr></tbody></table>

We observe from Table 5 that the values of $S_{\beta}^{R_0}$ and $S_{\Lambda}^{R_0}$ are exactly +1. This indicates that a rise in $\beta$, and $\Lambda$ will result in an increase in $R_0$ that is proportionate to both parameters. Additionally, we demonstrate that the parameters $k$, $\epsilon$, $\alpha$, and $\delta$ are exactly proportional to $R_0$ because $S_{k}^{R_0} > 0$, $S_{\epsilon}^{R_0} > 0$, $S_{\alpha}^{R_0} > 0$, and $S_{\delta}^{R_0} > 0$. Moreover, the terms $S_{\mu}^{R_0} < 0$, $S_{\gamma}^{R_0} < 0$, $S_{\sigma}^{R_0} < 0$, $S_{p}^{R_0} < 0$, and $S_{\eta}^{R_0} < 0$ denote that the parameters $\mu$, $\gamma$, $\sigma$, $p$, and $\eta$ are inversely proportional to $R_0$.

## 4. Results and Discussion

The results of parameters estimation are reported in Table 3, and Figure 4 illustrates the incidence data along with the model-fitted curve, obtained using the values in Table 3. The goodness fit of our model is supported by a height value coefficient of determination, namely $R^2 = 0.9016$; this indicates that the model fits the reported real data well.
## 3.2. Sensitivity Analysis

By using sensitivity, the spread and prevalence of diseases can be analyzed for each parameter. As a result of errors in data collection and assumed parameters, it is commonly used to measure the robustness of model predictions. To determine the relative significance of these parameters on disease transmission, we examined the impact of various model parameters. It has been determined how the model parameters $\beta$, $\gamma$, $k$, and $\epsilon$ affect the partial derivatives of the basic reproduction number $R_0$. Since $\frac{\partial R_0}{\partial \beta} > 0$, it follows that the transmission rate can be lowered to lessen the infection. However, given the partial derivatives $\frac{\partial R_0}{\partial \gamma} < 0$, it is implied that TB infection can be controlled by raising the parameter $\gamma$. Increasing $k$ results in an increase in $R_0$ because $\frac{\partial R_0}{\partial k} > 0$. It demonstrates that the diseased population will grow more quickly. Parameter $\alpha$ is the failure rate of treatment and $\frac{\partial R_0}{\partial \alpha} > 0$. Therefore, by lowering the treatment failure rate $\alpha$, the cumulative number of infected people can be minimized.

Sensitivity indices should be computed to estimate the relative change in a variable when parameters change. These indications were calculated using the following definition.

**Definition 1.** For a certain value $\sigma$, the normalized forward sensitivity index of $R_0$ is determined by

$$S_{\sigma}^{R_0} = \frac{\sigma}{R_0} \frac{\partial R_0}{\partial \sigma} \qquad (18)$$

For the baseline model parameters are determined by Formula (18), the derived sensitivity indices of the basic reproduction number $R_0$ are shown in Table 5.

<table><caption>**Table 5.** Parameters and sensitivity index.</caption><thead><tr><td>Parameter</td><td>Sensitivity Index</td></tr></thead><tbody><tr><td>Λ</td><td>+1</td></tr><tr><td>μ</td><td>-1.6502</td></tr><tr><td>k</td><td>+0.0012</td></tr><tr><td>β</td><td>+1</td></tr><tr><td>γ</td><td>-0.1671</td></tr><tr><td>ε</td><td>+0.0005</td></tr><tr><td>α</td><td>+2.1364 × 10<sup>-10</sup></td></tr><tr><td>δ</td><td>+1.4003 × 10<sup>-09</sup></td></tr><tr><td>σ</td><td>-0.4043</td></tr><tr><td>η</td><td>-2.5311 × 10<sup>-11</sup></td></tr><tr><td>p</td><td>-0.0194</td></tr></tbody></table>

We observe from Table 5 that the values of $S_{\beta}^{R_0}$ and $S_{\Lambda}^{R_0}$ are exactly +1. This indicates that a rise in $\beta$, and $\Lambda$ will result in an increase in $R_0$ that is proportionate to both parameters. Additionally, we demonstrate that the parameters $k$, $\epsilon$, $\alpha$, and $\delta$ are exactly proportional to $R_0$ because $S_{k}^{R_0} > 0$, $S_{\epsilon}^{R_0} > 0$, $S_{\alpha}^{R_0} > 0$, and $S_{\delta}^{R_0} > 0$. Moreover, the terms $S_{\mu}^{R_0} < 0$, $S_{\gamma}^{R_0} < 0$, $S_{\sigma}^{R_0} < 0$, $S_{p}^{R_0} < 0$, and $S_{\eta}^{R_0} < 0$ denote that the parameters $\mu$, $\gamma$, $\sigma$, $p$, and $\eta$ are inversely proportional to $R_0$.

## 4. Results and Discussion

The results of parameters estimation are reported in Table 3, and Figure 4 illustrates the incidence data along with the model-fitted curve, obtained using the values in Table 3. The goodness fit of our model is supported by a height value coefficient of determination, namely $R^2 = 0.9016$; this indicates that the model fits the reported real data well.
Using the estimate parameters value, one obtains $R_0 = 0.5228$, which is less than 1. This suggests that there is a possibility of decreasing or eliminating the disease by maintaining effective treatment and isolation measures in the future as illustrated by the model fitted curve for the period times 2020–2050 in Figure 4.

On the contrary, assuming that the government stops enforcing strict health measures against tuberculosis for children, such as vaccination, effective treatment strategies, and isolation of infected individuals after the year 2020, let us consider a hypothetical scenario. For instance, let us assume that $p = 10^{-2}$, $\gamma = 10^{-3}$, and $\beta = 4 \times 10^{-10}$. The other parameters are taken from Table 3. With this parameter set, the basic reproduction number is found to be greater than one $R_0 = 3.248 > 1$, indicating that the non-endemic equilibrium $E_1$ is unstable, whereas the endemic equilibrium $E_2$ is asymptotically stable. Obviously, the solutions of model (1) converge to $E_2$ as shown in Figure 5.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The expected situation if the government decides to abandon strict health measures after 2020, where the used parameters are $\beta = 4 \times 10^{-10}$, $\gamma = 10^{-3}$, $\sigma = 0.136$, $\eta = 4.2326 \times 10^{-6}$, $\delta = 0.1325$, $\alpha = 0.1095$, $\epsilon = 0.0656$, $k = 0.25$, and $p = 10^{-2}$.</figcaption></figure>

In order to gain a deeper understanding of how certain parameters affect the spread of the disease, one plots $R_0$ versus six parameters as shown in Figure 6. Obviously, there is a proportional relationship between the basic reproduction number $R_0$ and the three parameters $\beta$, $\epsilon$, and $\alpha$. This suggests that an increase in any of these parameters will result in an increase in the basic reproduction number, which in turn will lead to a greater spread of the disease.

On the other hand, this study found that there is an inverse relationship between the basic reproduction number $R_0$ and the other three parameters $\gamma$, $p$ and $\mu$. This implies that an increase in any of these parameters will result in a decrease in the basic reproduction number, which will lead to a slower spread of the disease. It is important to note that these results are in good agreement with reality.

The findings indicate that any plan aimed to prevent the spread of TB must consider four factors:

* Enhancing the precision and quality of TB diagnosis to facilitate appropriate actions towards affected individuals.

* Imposing isolation measures on infected individuals and monitoring their families medically to minimize contact with contagious patients.

* Sustaining a high rate of vaccination for children to provide immunity.

* Increasing the treatment rate by training specialized doctors, acquiring the most potent medicines, and establishing dedicated facilities to combat this disease.
Using the estimate parameters value, one obtains $R_0 = 0.5228$, which is less than 1. This suggests that there is a possibility of decreasing or eliminating the disease by maintaining effective treatment and isolation measures in the future as illustrated by the model fitted curve for the period times 2020–2050 in Figure 4.

On the contrary, assuming that the government stops enforcing strict health measures against tuberculosis for children, such as vaccination, effective treatment strategies, and isolation of infected individuals after the year 2020, let us consider a hypothetical scenario. For instance, let us assume that $p = 10^{-2}$, $\gamma = 10^{-3}$, and $\beta = 4 \times 10^{-10}$. The other parameters are taken from Table 3. With this parameter set, the basic reproduction number is found to be greater than one $R_0 = 3.248 > 1$, indicating that the non-endemic equilibrium $E_1$ is unstable, whereas the endemic equilibrium $E_2$ is asymptotically stable. Obviously, the solutions of model (1) converge to $E_2$ as shown in Figure 5.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** The expected situation if the government decides to abandon strict health measures after 2020, where the used parameters are $\beta = 4 \times 10^{-10}$, $\gamma = 10^{-3}$, $\sigma = 0.136$, $\eta = 4.2326 \times 10^{-6}$, $\delta = 0.1325$, $\alpha = 0.1095$, $\epsilon = 0.0656$, $k = 0.25$, and $p = 10^{-2}$.</figcaption></figure>

In order to gain a deeper understanding of how certain parameters affect the spread of the disease, one plots $R_0$ versus six parameters as shown in Figure 6. Obviously, there is a proportional relationship between the basic reproduction number $R_0$ and the three parameters $\beta$, $\epsilon$, and $\alpha$. This suggests that an increase in any of these parameters will result in an increase in the basic reproduction number, which in turn will lead to a greater spread of the disease.

On the other hand, this study found that there is an inverse relationship between the basic reproduction number $R_0$ and the other three parameters $\gamma$, $p$ and $\mu$. This implies that an increase in any of these parameters will result in a decrease in the basic reproduction number, which will lead to a slower spread of the disease. It is important to note that these results are in good agreement with reality.

The findings indicate that any plan aimed to prevent the spread of TB must consider four factors:

* Enhancing the precision and quality of TB diagnosis to facilitate appropriate actions towards affected individuals.

* Imposing isolation measures on infected individuals and monitoring their families medically to minimize contact with contagious patients.

* Sustaining a high rate of vaccination for children to provide immunity.

* Increasing the treatment rate by training specialized doctors, acquiring the most potent medicines, and establishing dedicated facilities to combat this disease.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Influence of the parameters $\beta, \gamma, p, \epsilon, \alpha$, and $\mu$, respectively, on the basic reproduction number $R_0$.</figcaption></figure>

**5. Conclusions**

In this study, we have developed a mathematical VSEIT model, which takes into account the biological factors of TB and certain realistic assumptions to analyze the transmission dynamics of this disease in Algeria. By using the reported infection data, we have estimated the model parameters using the least squares method. Our research has revealed that controlling the spread of TB is heavily dependent on some key factors, especially the contact parameter, $\beta$, the treatment parameter $\gamma$, and the vaccination parameter $p$. By identifying these crucial elements, we can better understand how to prevent and treat TB in Algeria. Using the estimated model parameters for Algeria we found that the reproduction number of the disease is less than one, meaning that TB can be eradicated by maintaining effective vaccination, high-quality treatment, and isolation measures. This research will have important implications for public health policymakers and healthcare professionals who are working to combat the spread of TB in Algeria and other regions where the disease is prevalent. The insights gained from this study can be used to develop more effective prevention and treatment strategies, which can ultimately help to reduce the burden of TB on affected communities.

**Author Contributions:** Conceptualization, M.S.A.; Methodology, M.S.A. and R.L.; Software, B.C.; Validation, R.L.; Investigation, B.C. and M.S.A.; Data curation, B.C.; Writing—original draft, B.C. and M.S.A.; Writing—review & editing, M.S.A. and R.L.; Supervision, M.S.A. and R.L.; Project administration, M.S.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the Algerian Directorate General for Scientific Research and Technological Development OF FUNDER grant number C00L03CU430120230003.

**Data Availability Statement:** The data presented in this study are available in [1,18,19,21].

**Conflicts of Interest:** The authors declare no conflict of interest.

**Appendix A. The Basic Reproduction Number $R_0$**

The basic reproduction number $R_0$ can be obtained using the next generation method [22], as $R_0 = \rho(FV^{-1})$.

The associated matrices $\mathcal{F}$ for new infection in the infected compartments and $\mathcal{V}$ for the remaining transfer terms are given respectively by
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Influence of the parameters $\beta, \gamma, p, \epsilon, \alpha$, and $\mu$, respectively, on the basic reproduction number $R_0$.</figcaption></figure>

**5. Conclusions**

In this study, we have developed a mathematical VSEIT model, which takes into account the biological factors of TB and certain realistic assumptions to analyze the transmission dynamics of this disease in Algeria. By using the reported infection data, we have estimated the model parameters using the least squares method. Our research has revealed that controlling the spread of TB is heavily dependent on some key factors, especially the contact parameter, $\beta$, the treatment parameter $\gamma$, and the vaccination parameter $p$. By identifying these crucial elements, we can better understand how to prevent and treat TB in Algeria. Using the estimated model parameters for Algeria we found that the reproduction number of the disease is less than one, meaning that TB can be eradicated by maintaining effective vaccination, high-quality treatment, and isolation measures. This research will have important implications for public health policymakers and healthcare professionals who are working to combat the spread of TB in Algeria and other regions where the disease is prevalent. The insights gained from this study can be used to develop more effective prevention and treatment strategies, which can ultimately help to reduce the burden of TB on affected communities.

**Author Contributions:** Conceptualization, M.S.A.; Methodology, M.S.A. and R.L.; Software, B.C.; Validation, R.L.; Investigation, B.C. and M.S.A.; Data curation, B.C.; Writing—original draft, B.C. and M.S.A.; Writing—review & editing, M.S.A. and R.L.; Supervision, M.S.A. and R.L.; Project administration, M.S.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the Algerian Directorate General for Scientific Research and Technological Development OF FUNDER grant number C00L03CU430120230003.

**Data Availability Statement:** The data presented in this study are available in [1,18,19,21].

**Conflicts of Interest:** The authors declare no conflict of interest.

**Appendix A. The Basic Reproduction Number $R_0$**

The basic reproduction number $R_0$ can be obtained using the next generation method [22], as $R_0 = \rho(FV^{-1})$.

The associated matrices $\mathcal{F}$ for new infection in the infected compartments and $\mathcal{V}$ for the remaining transfer terms are given respectively by
$$ \mathcal{F} = \begin{bmatrix} \beta SI \\ 0 \\ 0 \end{bmatrix}, \quad \mathcal{V} = \begin{bmatrix} (\epsilon + \mu)E - (1 - \alpha)\delta T \\ -\epsilon E - \alpha\delta T + (\gamma + \mu + \sigma)I \\ -\gamma I + (\mu + \delta + \eta)T \end{bmatrix}. $$

Let $k_1 = (\epsilon + \mu)$, $k_2 = (\gamma + \mu + \sigma)$, $k_3 = (\mu + \delta + \eta)$, then

$$ \mathcal{V} = \begin{bmatrix} k_1 E - (1 - \alpha)\delta T \\ -\epsilon E - \alpha\delta T + k_2 I \\ -\gamma I + k_3 T \end{bmatrix}. $$

Next, we evaluate F and V which are the Jacobian of $\mathcal{F}$ and $\mathcal{V}$ respectively at $E_1$ such that F is non-negative and V is a non-singular matrix.

We denote the Jacobian of $\mathcal{F}$ and $\mathcal{V}$ by $J(\mathcal{F})$ and $J(\mathcal{V})$ respectively. Thus,

$$ J(\mathcal{F}) = \begin{pmatrix} 0 & \beta S & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad J(\mathcal{V}) = \begin{bmatrix} k_1 & 0 & -(1-\alpha)\delta \\ -\epsilon & k_2 & -\alpha\delta \\ 0 & -\gamma & k_3 \end{bmatrix} $$

$F = J(\mathcal{F})$ and $V = J(\mathcal{V})$ at $E_1$, Thus,

$$ F = \begin{pmatrix} 0 & \frac{(k + \mu - \mu p) \Lambda \beta}{\mu(k + \mu)} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \text{and} \quad V = \begin{bmatrix} k_1 & 0 & -(1-\alpha)\delta \\ -\epsilon & k_2 & -\alpha\delta \\ 0 & -\gamma & k_3 \end{bmatrix} $$

It follows that

$$ V^{-1} = \frac{1}{(k_1 k_2 k_3 - \alpha \gamma \delta k_1 - (1-\alpha) \delta \gamma \epsilon)} \begin{bmatrix} k_2 k_3 - \alpha \delta \gamma & \gamma(1-\alpha) & -k_2(1-\alpha) \\ \epsilon k_3 & k_1 k_2 & k_1 \alpha \delta - (1-\alpha) \epsilon \\ \epsilon \gamma & \gamma k_1 & k_1 k_2 \end{bmatrix}, $$

then

$$ FV^{-1} = \frac{1}{(k_1 k_2 k_3 - \alpha \gamma \delta k_1 - (1 - \alpha) \delta \gamma \epsilon)} \begin{bmatrix} \frac{\epsilon(k+\mu-\mu p)\Lambda\beta}{\mu(k+\mu)} k_3 & \frac{\epsilon(k+\mu-\mu p)\Lambda\beta}{\mu(k+\mu)} k_1 k_2 & \frac{\epsilon(k+\mu-\mu p)\Lambda\beta(k_1\alpha\delta-\epsilon(1-\alpha))}{\mu(k+\mu)} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

Thus,

$$ \mathcal{R}_0 = \rho(FV^{-1}) = \frac{\epsilon(k + \mu - \mu p)\Lambda\beta k_3}{\mu(k + \mu)(k_1 k_2 k_3 - \alpha\gamma\delta k_1 - (1 - \alpha)\delta\gamma\epsilon)}, $$

## Appendix B. (LaSalle's Invariance Principle)

Instead of solely focusing on the stability of a specific equilibrium point, this principle offers insights into the overall behavior and trajectories of the system's solutions. Consider the autonomous system

$$ \dot{x} = f(x), \quad x \in \mathbb{R}^n, \tag{A1} $$

where f is of class $C^1$. The LaSalle's invariance principle [23] is stated in the following theorem.

**Theorem A1.** Let $\Omega \subset \mathbb{R}^n$ be a bounded closed (compact) set with the property that every solution of (A1) which begins in $\Omega$ remains for all future time in $\Omega$.

Suppose there is a scalar function $V(x)$ which has continuous first partials in $\Omega$ and is such that $\dot{V}(x) \le 0$ in $\Omega$. Let E be the set of all points in $\Omega$ where $\dot{V}(x) = 0$.

Let M be the largest invariant set in E. Then, every solution starting in $\Omega$ approaches M as $t \to +\infty$.
$$ \mathcal{F} = \begin{bmatrix} \beta SI \\ 0 \\ 0 \end{bmatrix}, \quad \mathcal{V} = \begin{bmatrix} (\epsilon + \mu)E - (1 - \alpha)\delta T \\ -\epsilon E - \alpha\delta T + (\gamma + \mu + \sigma)I \\ -\gamma I + (\mu + \delta + \eta)T \end{bmatrix}. $$

Let $k_1 = (\epsilon + \mu)$, $k_2 = (\gamma + \mu + \sigma)$, $k_3 = (\mu + \delta + \eta)$, then

$$ \mathcal{V} = \begin{bmatrix} k_1 E - (1 - \alpha)\delta T \\ -\epsilon E - \alpha\delta T + k_2 I \\ -\gamma I + k_3 T \end{bmatrix}. $$

Next, we evaluate F and V which are the Jacobian of $\mathcal{F}$ and $\mathcal{V}$ respectively at $E_1$ such that F is non-negative and V is a non-singular matrix.

We denote the Jacobian of $\mathcal{F}$ and $\mathcal{V}$ by $J(\mathcal{F})$ and $J(\mathcal{V})$ respectively. Thus,

$$ J(\mathcal{F}) = \begin{pmatrix} 0 & \beta S & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad J(\mathcal{V}) = \begin{bmatrix} k_1 & 0 & -(1-\alpha)\delta \\ -\epsilon & k_2 & -\alpha\delta \\ 0 & -\gamma & k_3 \end{bmatrix} $$

$F = J(\mathcal{F})$ and $V = J(\mathcal{V})$ at $E_1$, Thus,

$$ F = \begin{pmatrix} 0 & \frac{(k + \mu - \mu p)\Lambda\beta}{\mu(k + \mu)} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \text{and} \quad V = \begin{bmatrix} k_1 & 0 & -(1-\alpha)\delta \\ -\epsilon & k_2 & -\alpha\delta \\ 0 & -\gamma & k_3 \end{bmatrix} $$

It follows that

$$ V^{-1} = \frac{1}{(k_1 k_2 k_3 - \alpha \gamma \delta k_1 - (1-\alpha) \delta \gamma \epsilon)} \begin{bmatrix} k_2 k_3 - \alpha \delta \gamma & \gamma(1-\alpha) & -k_2(1-\alpha) \\ \epsilon k_3 & k_1 k_2 & k_1 \alpha \delta - (1-\alpha) \epsilon \\ \epsilon \gamma & \gamma k_1 & k_1 k_2 \end{bmatrix}, $$

then

$$ FV^{-1} = \frac{1}{(k_1 k_2 k_3 - \alpha \gamma \delta k_1 - (1 - \alpha) \delta \gamma \epsilon)} \begin{bmatrix} \frac{\epsilon(k+\mu-\mu p)\Lambda\beta}{\mu(k+\mu)} k_3 & \frac{\epsilon(k+\mu-\mu p)\Lambda\beta}{\mu(k+\mu)} k_1 k_2 & \frac{\epsilon(k+\mu-\mu p)\Lambda\beta(k_1\alpha\delta-\epsilon(1-\alpha))}{\mu(k+\mu)} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

Thus,

$$ \mathcal{R}_0 = \rho(FV^{-1}) = \frac{\epsilon(k + \mu - \mu p)\Lambda\beta k_3}{\mu(k + \mu)(k_1 k_2 k_3 - \alpha\gamma\delta k_1 - (1 - \alpha)\delta\gamma\epsilon)}, $$

## Appendix B. (LaSalle's Invariance Principle)

Instead of solely focusing on the stability of a specific equilibrium point, this principle offers insights into the overall behavior and trajectories of the system's solutions. Consider the autonomous system

$$ \dot{x} = f(x), \quad x \in \mathbb{R}^n, \qquad (A1) $$

where f is of class C¹. The LaSalle's invariance principle [23] is stated in the following theorem.

**Theorem A1.** Let $\Omega \subset \mathbb{R}^n$ be a bounded closed (compact) set with the property that every solution of (A1) which begins in $\Omega$ remains for all future time in $\Omega$.

Suppose there is a scalar function $V(x)$ which has continuous first partials in $\Omega$ and is such that $\dot{V}(x) \le 0$ in $\Omega$. Let E be the set of all points in $\Omega$ where $\dot{V}(x) = 0$.

Let M be the largest invariant set in E. Then, every solution starting in $\Omega$ approaches M as $t \to +\infty$.
References

1. WHO. Global Tuberculosis Report; World Health Organization: Geneva, Switzerland, 2023. Available online: https://extranet.who.int/tme/generateCSV.asp?ds=notifications (accessed on 1 January 2023).

2. Bernoulli, D. Essai d'une nouvelle analyse de la mortalité causée par la petite vérole, et des avantages de l'inoculation pour la prévenir. *Histoire de l'Acad. Roy. Sci. (Paris) avec Mem 1760*; pp. 1-45. Available online: https://gallica.bnf.fr/ark:/12148/bpt6k3558n/f223.item.r=daniel%20bernoulli (accessed on 1 January 2023).

3. Ross, R. The Prevention of Malaria. *Nature* **1910**, *85*, 263-264. [CrossRef]

4. Hamer, W.H. *The Milroy Lectures on Epidemic Diseases in England: The Evidence of Variability and of Persistency of Type*; The Bedford Press: London, UK, 1906; pp. 733-739.

5. Martini, E. *Berechnungen und Beobachtungen zur Epidemiologie und Bekämpfung der Malaria*; W. Gente: Hamburg, Germany, 1921.

6. Lotka A.J. *Elements of Physical Biology*; Williams and Wilkens: Philadelphia, PA, USA, 1925.

7. Yang, Y.; Tang, S.; Ren, X.; Zhao, H.; Guo, C. Global stability and optimal control for a tuberculosis model with vaccination and treatment. *Discret. Contin. Dyn. Syst. Ser. B* **2016**, *21*, 1009–1022. [CrossRef]

8. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. *Proc. R. Soc. Lond. Ser. A Contain. Pap. Math. Phys. Character* **1927**, *115*, 700–721.

9. Waaler, H.; Geser, A.; Andersen, S. The use of mathematical models in the study of the epidemiology of tuberculosis. *Am. J. Public Health Nation's Health* **1962**, *52*, 1002–1013. [CrossRef] [PubMed]

10. Ullah, S.; Khan, M.A.; Farooq, M.; Gul, T. Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan. *Math. Comput. Simul.* **2019**, *165*, 182–184. [CrossRef]

11. Abdelouahab, M.S.; Arama, A.; Lozi, R. Bifurcation analysis of a model of tuberculosis epidemic with treatment of wider population suggesting a possible role in the seasonality of this disease. *Chaos Interdiscip. J. Nonlinear Sci.* **2021**, *31*, 123125. [CrossRef] [PubMed]

12. Andersen, P.; Doherty, T.M. The success and failure of BCG-implications for novel tuberculosis vaccine. *Nature* **2005**, *3*, 656–662. [CrossRef] [PubMed]

13. Ucakan, Y.; Gulen, S.; Koklu, K. Analysing of Tuberculosis in Turkey through SIR, SEIR and BSEIR Mathematical Models. *Math. Comput. Model. Dyn. Syst.* **2021**, *27*, 179–202. [CrossRef]

14. Egonmwan, A.O.; Okuonghae, D. Mathematical analysis of a tuberculosis model with imperfect vaccine. *Int. J. Biomath.* **2019**, *12*, 1950073. [CrossRef]

15. Revelle, C.S.; Lynn, W.R.; Feldmann, F. Mathematical models for the economic allocation of tuberculosis control activities in developing nations. *Am. Rev. Respir. Dis.* **1967**, *96*, 893–909.

16. Diekmann, O.; Heesterbeek, J.A.P.; Metz, J.A. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. *J. Math. Biol.* **1990**, *28*, 365–382. [CrossRef]

17. Wang, Y.; Cao, J. Global stability of general cholera models with nonlinear incidence and removal rates. *J. Frankl. Inst.* **2015**, *352*, 2464–2485. [CrossRef]

18. Population Growth in Algeria. Available online: https://www.donneesmondiales.com/afrique/algerie/croissance-population.php (accessed on 10 January 2022).

19. Trading Economics. Immunization, BCG (% of One-Year-Old Children) Algeria. Available online: https://www.tradingeconomics.com/algeria/immunization-bcg-percent-of-one-year-old-children-wb-data.html/ (accessed on 20 January 2022).

20. Katelaris, A.L.; Jackson, C.; Southern, J.; Gupta, R.K.; Drobniewski, F.; Lalvani, A.; Lipman, M.; Mangtani, P.; Abubakar, I. Effectiveness of BCG Vaccination Against Mycobacterium tuberculosis Infection in Adults: A Cross-sectional Analysis of a UK-Based Cohort. *J. Infect. Dis.* **2020**, *221*, 146–155. [CrossRef]

21. The World Bank Group. Tuberculosis Treatment Success Rate (% of New Cases)—Algeria. Available online: https://data.worldbank.org/indicator/SH.TBS.CURE.ZS?locations=DZ (accessed on 5 January 2022).

22. Van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Math. Biosci.* **2002**, *180*, 29–48. [CrossRef]

23. LaSalle, J.P. Some extensions of Liapunov's second method. *IRE Trans. Circuit Theory* **1960**, *7*, 520–527. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
References

1. WHO. *Global Tuberculosis Report*; World Health Organization: Geneva, Switzerland, 2023. Available online: https://extranet.who.int/tme/generateCSV.asp?ds=notifications (accessed on 1 January 2023).

2. Bernoulli, D. Essai d'une nouvelle analyse de la mortalité causée par la petite vérole, et des avantages de l'inoculation pour la prévenir. *Histoire de l'Acad. Roy. Sci. (Paris) avec Mem 1760*; pp. 1-45. Available online: https://gallica.bnf.fr/ark:/12148/bpt6k3558n/f223.item.r=daniel%20bernoulli (accessed on 1 January 2023).

3. Ross, R. The Prevention of Malaria. *Nature* **1910**, *85*, 263-264. [CrossRef]

4. Hamer, W.H. *The Milroy Lectures on Epidemic Diseases in England: The Evidence of Variability and of Persistency of Type*; The Bedford Press: London, UK, 1906; pp. 733-739.

5. Martini, E. *Berechnungen und Beobachtungen zur Epidemiologie und Bekämpfung der Malaria*; W. Gente: Hamburg, Germany, 1921.

6. Lotka A.J. *Elements of Physical Biology*; Williams and Wilkens: Philadelphia, PA, USA, 1925.

7. Yang, Y.; Tang, S.; Ren, X.; Zhao, H.; Guo, C. Global stability and optimal control for a tuberculosis model with vaccination and treatment. *Discret. Contin. Dyn. Syst. Ser. B* **2016**, *21*, 1009–1022. [CrossRef]

8. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. *Proc. R. Soc. Lond. Ser. A Contain. Pap. Math. Phys. Character* **1927**, *115*, 700–721.

9. Waaler, H.; Geser, A.; Andersen, S. The use of mathematical models in the study of the epidemiology of tuberculosis. *Am. J. Public Health Nation's Health* **1962**, *52*, 1002–1013. [CrossRef] [PubMed]

10. Ullah, S.; Khan, M.A.; Farooq, M.; Gul, T. Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan. *Math. Comput. Simul.* **2019**, *165*, 182–184. [CrossRef]

11. Abdelouahab, M.S.; Arama, A.; Lozi, R. Bifurcation analysis of a model of tuberculosis epidemic with treatment of wider population suggesting a possible role in the seasonality of this disease. *Chaos Interdiscip. J. Nonlinear Sci.* **2021**, *31*, 123125. [CrossRef] [PubMed]

12. Andersen, P.; Doherty, T.M. The success and failure of BCG-implications for novel tuberculosis vaccine. *Nature* **2005**, *3*, 656–662. [CrossRef] [PubMed]

13. Ucakan, Y.; Gulen, S.; Koklu, K. Analysing of Tuberculosis in Turkey through SIR, SEIR and BSEIR Mathematical Models. *Math. Comput. Model. Dyn. Syst.* **2021**, *27*, 179–202. [CrossRef]

14. Egonmwan, A.O.; Okuonghae, D. Mathematical analysis of a tuberculosis model with imperfect vaccine. *Int. J. Biomath.* **2019**, *12*, 1950073. [CrossRef]

15. Revelle, C.S.; Lynn, W.R.; Feldmann, F. Mathematical models for the economic allocation of tuberculosis control activities in developing nations. *Am. Rev. Respir. Dis.* **1967**, *96*, 893–909.

16. Diekmann, O.; Heesterbeek, J.A.P.; Metz, J.A. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. *J. Math. Biol.* **1990**, *28*, 365–382. [CrossRef]

17. Wang, Y.; Cao, J. Global stability of general cholera models with nonlinear incidence and removal rates. *J. Frankl. Inst.* **2015**, *352*, 2464–2485. [CrossRef]

18. Population Growth in Algeria. Available online: https://www.donneesmondiales.com/afrique/algerie/croissance-population.php (accessed on 10 January 2022).

19. Trading Economics. Immunization, BCG (% of One-Year-Old Children) Algeria. Available online: https://www.tradingeconomics.com/algeria/immunization-bcg-percent-of-one-year-old-children-wb-data.html/ (accessed on 20 January 2022).

20. Katelaris, A.L.; Jackson, C.; Southern, J.; Gupta, R.K.; Drobniewski, F.; Lalvani, A.; Lipman, M.; Mangtani, P.; Abubakar, I. Effectiveness of BCG Vaccination Against Mycobacterium tuberculosis Infection in Adults: A Cross-sectional Analysis of a UK-Based Cohort. *J. Infect. Dis.* **2020**, *221*, 146–155. [CrossRef]

21. The World Bank Group. Tuberculosis Treatment Success Rate (% of New Cases)—Algeria. Available online: https://data.worldbank.org/indicator/SH.TBS.CURE.ZS?locations=DZ (accessed on 5 January 2022).

22. Van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Math. Biosci.* **2002**, *180*, 29–48. [CrossRef]

23. LaSalle, J.P. Some extensions of Liapunov's second method. *IRE Trans. Circuit Theory* **1960**, *7*, 520–527. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.