Mathematical Modeling and Sensitivity Analysis of COVID-19 and Tuberculosis Coinfection with Vaccination

COVID-19 and TB are both airborne infectious diseases. The model investigated in this paper is the co-infection model of COVID-19 and tuberculosis. The population is divided into eight subpopulations. Individuals in the susceptible subpopulation (S) are those who are still healthy. Vaccination subpopulation (V), consisting of healthy individuals immunized against the COVID-19 and BCG vaccines. Individuals who have been vaccinated and are immune to COVID-19 and TB join the R subpopulation. Individuals who are immune to COVID-19 but not to TB, or who are immune to TB but not to COVID-19, are classified as subpopulation V. Individuals who have been immunized against COVID-19 and BCG are vulnerable to COVID-19 and TB again and may revert to subpopulation S. Individuals who are still healthy come into contact with COVID-19 infected individuals and become infected with COVID-19, where they then enter the COVID-19 infected subpopulation, as denoted by I_C. Individuals who are still healthy come into contact with TB-infected individuals and become infected with TB, where they then enter the TB-infected subpopulation, as denoted by I_T. Individuals in the I_C subpopulation come into contact with TB-infected individuals and become infected with TB, where they then enter the co-infected subpopulation with COVID-19 and TB, as denoted by I_CT. Individuals infected with COVID-19 who have not recovered after being given vitamins are classified as Q. Individuals who have TB but are not infected with COVID-19 will enter the tuberculosis treatment subpopulation, denoted by T. Individuals infected with COVID or TB, or both, will be treated and have a chance to recover before joining the recovered subpopulation, denoted by R.

Some additional assumptions for the model in this study are as follows: (1) Constant recruitment into the S subpopulation. (2) Susceptible individuals who are vaccinated and immune will join the R subpopulation. (3) Individuals infected with COVID-19 or TB may die as a result of the disease. (4) Each subpopulation may die naturally. (5) Infected individuals may recover as a result of treatment. (6) Each subpopulation is homogeneous. Figure 1 depicts the spread of COVID-19 and TB coinfection.

1111111.png

Figure 1. Transmission diagram of COVID-19 and tuberculosis coinfection model

Based on Figure 1, the mathematical model of the COVID-19 and TB coinfection is as follows:

$\begin{aligned} & \frac{d S}{d t}=\Lambda+\varphi V-\left(\theta+\xi_1+\xi_2+\mu\right) S \\ & \frac{d V}{d t}=\theta S-\left(\varphi+r_1+\mu\right) V \\ & \frac{d I_C}{d t}=\xi_1 S+\omega I_{C T}-\left(\alpha \xi_2+t_1+\delta+d_1+\mu\right) I_C \\ & \frac{d I_T}{d t}=\xi_2 S+\eta I_{C T}-\left(\gamma \xi_1+\rho+d_2+\mu\right) I_T \\ & \frac{d I_{C T}}{d t}=\alpha \xi_2 I_C+\gamma \xi_1 I_T-\left(\omega+\eta+\psi+d_3+\mu\right) I_{C T} \\ & \frac{d Q}{d t}=\delta I_C-\left(\tau+d_4+\mu\right) Q \\ & \frac{d T}{d t}=\rho I_T-\left(\sigma+d_5+\mu\right) T \\ & \frac{d R}{d t}=r_1 V+t_1 I_C+\tau Q+\sigma T+\psi I_{C T}-\mu R\end{aligned}$                           (1)

where, $\xi_1=\frac{\beta_1\left(I_C+I_{C T}\right)}{N}$, and $\xi_2=\frac{\beta_2\left(I_T+T+I_{C T}\right)}{N}$, and with initial conditions:

$\begin{aligned} & S(0)>0, V(0)>0, I_C(0) \geq 0, I_T(0) \geq 0 \\ & I_{C T}(0) \geq 0, Q(0) \geq 0, T(0) \geq 0, R(0) \geq 0\end{aligned}$          (2)

The changes in the total population per unit time is given by:

$\frac{d N}{d t}=\Lambda-\mu N-\left(d_1+d_2+d_3+d_4+d_5\right)$                 (3)

Notation, description, and value of parameters are presented in Table 1.

All variables in model (1) have non-negative solutions with non-negative initial conditions (2) in the following bounded region,

$\Upsilon=\left\{\left(S, V, I_C, I_T, I_{C T}, Q, T, R\right) \in R_{+}^8, N \leq \frac{\Lambda}{\mu}\right\}$                (4)

It will be shown that the solutions of the model with non-negative initial conditions will always be non-negative for each t > 0, and it is stated in Theorem 1 below. The theorem can be interpreted in biologically meanings that population size is never be negative.

Theorem 1

Based on the initial conditions in Eq. (3), at time t, all subpopulations, $S(t), V(t), I_C(t), I_T(t), I_{C T}(t), Q(t), T(t), R(t)$ are nonnegative.

Proof.

Define $\zeta=\sup \left\{S(t), V(t)>0, I_C(t)>0, I_T(0)>0, I_{C T}(t)>0, Q(t)>0, T(t)>0, R(t)>0\right\}$. Because $S(t), V(t), I_C(t)$, $I_T(t), I_{C T}(t), Q(t), T(t), R(t)$ are continuous, then $\zeta>0$.

If $\zeta=+\infty$, then the positivity prevails. If $0<\zeta<+\infty$, then $S(\zeta)=0, V(\zeta)=0, I_T(\zeta)=0, I_{C T}(\zeta)=0, Q(\zeta)=0, T(\zeta)=0$, and $R(\zeta)=0$.

From the first differential equation of system (1), $\frac{d S}{d t}=\Lambda+\varphi V-\left(\theta+\xi_1+\xi_2+\mu\right) S$, by integrating both sides, we get $S(\zeta)=$ $Q_1 S(0)+Q_1 \int_0^\zeta e^{-\int\left(\theta+\xi_1+\xi_2+\mu\right) d t}(\Lambda+\varphi V) d t$.

Because $e^{-\int\left(\theta+\xi_1+\xi_2+\mu\right) d t}>0$ and $S(0)>0$, and also from definition of $\zeta$, we have $V(t)>0, I_C(t)>0, I_T(t)>0, I_{C T}(t)>$ $0, Q(t)>0, T(t)>0, R(t)>0$, then $S(\zeta)>0$, and hence $S(\zeta) \neq 0$.

In the similar way, for the second until the eighth differential equation of system (1), we will obtain $(t)>0, I_C(t)>0, I_T(t)>$ $0, I_{C T}(t)>0, Q(t)>0, T(t)>0$, and $R(t)>0$. Since $\zeta$ is infinite, this means that $\zeta=+\infty$, then all the solutions of the coinfection model of COVID-19 and TB model (1) are non-negative.

Table 1. Parameter values for COVID-19 and TB coinfection transmission

2.1 Nonendemic equilibrium point

The non-endemic point of the model (1) is obtained if the rate of spread of each subpopulation is constant and there are no COVID-19 and TB diseases. Non-endemic point is as follows:

$E_0=\left(S^*, V^*, 0,0,0,0,0,0, R^*\right)$

where,

$\begin{aligned} S^* & =\frac{\left(\varphi+r_1+\mu\right) \Lambda}{\mu\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)} \\ V^* & =\frac{\theta \Lambda}{\mu\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)} \\ R^* & =\frac{\theta r_1 \Lambda}{\mu\left(\mu\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)\right)} .\end{aligned}$

2.2 The basic reproduction number

The basic reproduction number of an infectious disease spread model is an important parameter for determining whether disease spread is increasing or decreasing. The basic reproduction number is commonly represented by $\mathcal{R}_0$. If $\mathcal{R}_0>1$ then the spread of the epidemic increases, and vice versa if $\mathcal{R}_0>1$ then the spread of the epidemic then slows or stops. The Jacobian matrix of individuals who are in contact between active and susceptible infected individuals is obtained using model (1) and the non-endemic equilibrium point, namely:

$F=\left[\begin{array}{cccc}q_1 & 0 & q_1 & 0 \\ 0 & q_2 & q_2 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0\end{array}\right]$

where, $q_1=\frac{\beta_1 \Lambda\left(\varphi+r_1+\mu\right)}{\mu N\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)}$ and $q_2=\frac{\beta_2 \Lambda\left(\varphi+r_1+\mu\right)}{\mu N\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)}$. Whereas, the Jacobian matrix of non-contact individuals between subpopulations and out is:

$V=\left[\begin{array}{cccc}t_1+\delta+d_1+\mu & 0 & -\omega & 0 \\ 0 & \rho+d_2+\mu & -\eta & 0 \\ 0 & 0 & \omega+\eta+\psi+d_3+\mu & 0 \\ 0 & -\rho & 0 & \sigma+d_5+\mu\end{array}\right]$

One method for determining $\mathcal{R}_0$ of an epidemic spread is to use the next generation matrix approach [44]. From this we get:

$\mathcal{R}_0=\mathcal{R}_{0 C}+\mathcal{R}_{0 T}$                  (5)

where, $\mathcal{R}_{0 C}=\frac{\beta_1 \Lambda\left(\varphi+r_1+\mu\right)}{\mu N\left(t_1+\delta+d_1+\mu\right)\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)}$ is the basic reproduction number of COVID-19 spread and $\mathcal{R}_{0 T}=\frac{\beta_2 \Lambda\left(\varphi+r_1+\mu\right)}{\mu N\left(\rho+d_2+\mu\right)\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)}$ is the basic reproduction number of TB spread.

From $\mathcal{R}_{0 C}$, it can be shown that the COVID-19 vaccine can reduce the number of individuals infected with COVID-19 at a rate of $p r_1$, resulting in a decrease in the basic reproduction number of COVID-19 spread.

From $\mathcal{R}_{0 C}$, it can be shown that the COVID-19 vaccine can reduce the number of individuals infected with COVID-19 at a rate of $p r_1$, resulting in a decrease in the basic reproduction number of COVID-19 spread.

$\mathcal{R}_{0 C}=\frac{\beta_1 \Lambda\left(\varphi+p r_1+\mu\right)}{\mu N\left(t_1+\delta+d_1+\mu\right)\left(\mu+p r_1+\varphi\right)+\theta\left(p r_1+\mu\right)}$

It will be demonstrated that the COVID vaccine can reduce the incidence of COVID-19 infection. Consider the basic reproduction number of the spread of COVID-19 in the absence of COVID-19 vaccination:

$\mathcal{R}_{0 C 1}=\frac{\beta_1 \Lambda(\varphi+\mu)}{\mu N\left(t_1+\delta+d_1+\mu\right)(\mu+\varphi)+\theta \mu}$

By substituting the values of the parameters in Table 1, we obtain

$\frac{\beta_1 \Lambda\left(\varphi+p r_1+\mu\right)}{\mu N\left(t_1+\delta+d_1+\mu\right)\left(\mu+p r_1+\varphi\right)+\theta\left(p r_1+\mu\right)}<\frac{\beta_1 \Lambda(\varphi+\mu)}{\mu N\left(t_1+\delta+d_1+\mu\right)(\mu+\varphi)+\theta \mu}$

or

$\mathcal{R}_{0 C}<\mathcal{R}_{0 C 1}$                  (6)

In the similar way, consider the basic reproduction number for the spread of tuberculosis in the absence of BCG vaccination:

$\mathcal{R}_{0 T}=\frac{\beta_2 \Lambda(\varphi+\mu)}{\mu N\left(\rho+d_2+\mu\right)(\mu+\varphi)+\theta \mu}$

It will be demonstrated that BCG vaccination can reduce the spread of tuberculosis at a rate of $(1-p) r_1$, so the basic reproduction number for the spread of tuberculosis becomes

$\mathcal{R}_{0 T 1}=\frac{\beta_2 \Lambda\left(\varphi+(1-p) r_1+\mu\right)}{\mu N\left(\rho+d_2+\mu\right)\left(\mu+(1-p) r_1+\varphi\right)+\theta\left((1-p) r_1+\mu\right)}$

The basic reproduction number of tuberculosis spread without BCG vaccination is calculated by entering the parameter values in Table 1,

$\frac{\beta_2 \Lambda\left(\varphi+(1-p) r_1+\mu\right)}{\mu N\left(\rho+d_2+\mu\right)\left(\mu+(1-p) r_1+\varphi\right)+\theta\left((1-p) r_1+\mu\right)}<\frac{\beta_2 \Lambda(\varphi+\mu)}{\mu N\left(\rho+d_2+\mu\right)(\mu+\varphi)+\theta \mu}$

or

$\mathcal{R}_{0 T}<\mathcal{R}_{0 T 1}$                   (7)

Based on Eqs. (6) and (7), it can be concluded that vaccination can reduce the number of COVID-19, TB, and COVID-19 and TB co-infections.

The following is a theorem regarding the local stability of nonendemic equilibrium point. The interpretation of the theorem is that the diseases’ infection will be vanished at certain conditions.

Theorem 2

The nonendemic equilibrium point $E_0$S of the COVID-19 and TB co-infection model is asymptotically stable locally if $\mathcal{R}_0<1$ and vice versa.

Proof

Consider the Jacobian matrix model of Eq. (1) at the nonendemic equilibrium point $E_0$.

$J\left(E_0\right)=\left[\begin{array}{cccccccc}-\theta-\mu & \varphi & -S_1 & -S_2 & -S_1-S_2 & 0 & -S_2 & 0 \\ \theta & -h_1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & S_1-h_2 & 0 & S_1+\omega & 0 & 0 & 0 \\ 0 & 0 & 0 & S_2-h_3 & S_2+\eta & 0 & S_2 & 0 \\ 0 & 0 & 0 & 0 & -h_4 & 0 & 0 & 0 \\ 0 & 0 & \delta & 0 & 0 & -h_5 & 0 & 0 \\ 0 & 0 & 0 & \rho & 0 & 0 & -h_6 & 0 \\ 0 & v_1 & t_1 & 0 & \psi & \tau & \sigma & -\mu\end{array}\right]$

where,

$\begin{aligned} & S_1=\frac{\left(\varphi+r_1+\mu\right) \beta_1 \Lambda}{\mu N\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)}, \\ & S_2=\frac{\left(\varphi+r_1+\mu\right) \beta_2 \Lambda}{\mu N\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)}, \\ & h_1=\varphi+r_1+\mu, \\ & h_2=t_1+\delta+d_1+\mu, \\ & h_3=\rho+d_2+\mu, \\ & h_4=\omega+\eta+\psi+d_3+\mu, \\ & h_5=\tau+d_4+\mu, \\ & h_6=\sigma+d_5+\mu .\end{aligned}$

The eigenvalues are obtained from the Jacobian matrix,

$\begin{aligned} & \lambda_1=-\mu, \lambda_2=-\theta-\mu, \\ & \lambda_3=-\left(\tau+d_4+\mu\right), \\ & \lambda_4=-\left(\omega+\eta+\psi+d_3+\mu\right), \\ & \lambda_5=-\frac{\left(\varphi+r_1+\mu\right) \beta_1 \Lambda}{\mu N\left(\mu+r_1+\varphi\right)+\theta\left(r_1+\mu\right)}-t_1+\delta+d_1+\mu, \\ & \lambda_6=-\left(\varphi+r_1+\mu\right), \\ & \lambda_7=\frac{1}{2}\left(-h_3-h_6+S_2\right)+ \frac{1}{2} \sqrt{S_2^2+\left(-2 h_3+2 h_6+4 \rho\right) S_2+\left(h_3-h_6\right)^2}, \\ & \lambda_8=\frac{1}{2}\left(-h_3-h_6+S_2\right)- \frac{1}{2} \sqrt{S_2^2+\left(-2 h_3+2 h_6+4 \rho\right) S_2+\left(h_3-h_6\right)^2} .\end{aligned}$

The value of $\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_5, \lambda_6$ and $\lambda_8$ are clearly negative, meanwhile for $\lambda_7$, By entering the parameter values listed in Table 1, $\lambda_7$ becomes negative as well. Because all the eigen values of the Jacobian matrix at the point of nonendemic equilibrium are negative, corresponding to $\mathcal{R}_0<1$, and vice versa. Therefore, the nonendemic equilibrium point $E_0$  in the COVID-19 and TB co-infection model (1) is asymptotically stable locally, and vice versa.

2.3 Endemic equilibrium point

When the number of each subpopulation remains constant and the disease spreads, endemic equilibrium is reached as follows

$E_1=\left(S^{* *}, V^{* *}, E_C^{* *}, E_T^{* *}, E_{C T}^{* *}, I_C^{* *}, I_T^{* *}, I_{C T}^{* *}, R^*\right)$

where,

$S^{* *}=\frac{k_2 \Lambda}{k_1 k_2-\theta \varphi}$

$V^{* *}=\frac{\theta \Lambda}{k_1 k_2-\theta \varphi}$

$I_C^{* *}=\frac{k_2 \Lambda \xi_1\left(\eta \gamma \xi_1-\alpha \omega \xi_2-k_4 k_5\right)}{\left(\alpha k_4 \omega \xi_2+\eta \gamma k_3 \xi_1-k_3 k_4 k_5\right)\left(k_1 k_2-\theta \varphi\right)}$,

$I_T^{* *}=\frac{\xi_2 k_2 \Lambda\left(\alpha \eta \xi_1-\alpha \omega \xi_2+k_3 k_5\right)}{\left(\alpha k_4 \omega \xi_2+\eta \gamma k_3 \xi_1-k_3 k_4 k_5\right)\left(k_1 k_2-\theta \varphi\right)}$,

$I_{C T}^{* *}=\frac{\xi_2 k_2 \Lambda \xi_1\left(\alpha k_4+\gamma k_3\right)}{\left(\alpha k_4 \omega \xi_2+\eta \gamma k_3 \xi_1-k_3 k_4 k_5\right)\left(k_1 k_2-\theta \varphi\right)}$

$Q^{* *}=\frac{\delta k_2 \Lambda \xi_1\left(\eta \gamma \xi_1-\gamma \omega \xi_2-k_4 k_5\right)}{k_6\left(\left[\alpha k_4 \omega \xi_2+\eta \gamma k_3 \xi_1-k_3 k_4 k_5\right]\left(k_1 k_2-\theta \varphi\right)\right)}$

$T^*=\frac{\xi_2 \rho k_2 \Lambda\left(\alpha \eta \xi_1-\alpha \omega \xi_2 \mp k_3 k_5\right)}{k_7\left(\left[\alpha k_4 \omega \xi_2+\eta \gamma k_3 \xi_1-k_3 k_4 k_5\right]\left(k_1 k_2-\theta \varphi\right)\right)}$

$\begin{aligned} & R^*=\frac{\left(\left[\left(\left[-\gamma t_1 \eta \xi_1\left(\left[k_3 \psi+\omega t_1\right] \gamma+\psi \alpha k_4\right) \xi_2+k_4 k_5 t_1\right] \xi_1 k_7+\rho \xi_2 \sigma\left(\alpha \eta \xi_1-\alpha \omega \xi_2+k_4 k_5\right)\right) k_6-\right.\right.}{\mu k_6 k_7\left(\left[\alpha k_4 \omega \xi_2+\eta \gamma k_3 \xi_1-k_3 k_4 k_5\right]\left(\theta \varphi-k_1 k_2\right)\right)} \\ & \frac{\left.\left.-\tau \xi_1 \delta k_7\left(\eta \gamma \xi_1-\gamma \omega \xi_2-k_4 k_5\right]\right) k_2-k_7 r_1\left[\xi_1 \gamma k_3 \eta+k_4\left(\alpha \omega \xi_2-k_3 k_5\right)\right] \theta k_6\right) \Lambda}{\mu k_6 k_7\left(\left[\alpha k_4 \omega \xi_2+\eta \gamma k_3 \xi_1-k_3 k_4 k_5\right]\left(\theta \varphi-k_1 k_2\right)\right)}\end{aligned}$

$\begin{aligned} & \text {with} k_1=\theta+\xi_1+\xi_2+\mu, \\ & k_2=\varphi+r_1+\mu \\ & k_3=\alpha \xi_2+t_1+\delta+d_1+\mu, \\ & k_4=\gamma \xi_1+\rho+d_2+\mu, \\ & k_5=\omega+\eta+\psi+d_3+\mu, \\ & k_6=\tau+d_4+\mu, \text { and } \\ & k_7=\sigma+d_5+\mu .\end{aligned}$

The following theorem states about the stability of the endemic equilibrium, and in biologically meanings, it means that in some conditions the diseases will continue to spread forever.

Theorem 3

If $\mathcal{R}_0>1$, then the endemic equilibrium $E_2$ is locally asymptotic stable.

Proof.

Based on the system of Eq. (1) the Jacobian Matrix from the system of Eq. (1) above at point $E_2$

$J\left(E_2^*\right)=\left[\begin{array}{cccccccc}-a_{11} & \varphi & -S_1^{* *} & -S_2^{* *} & -S_1^{* *}-S_2^{* *} & 0 & -S_2^{* *} & 0 \\ \theta & -h_1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \xi_1^{* *} & 0 & a_{33} & -\frac{\alpha \beta_2 I_C^{* *}}{N} & a_{35} & 0 & -\frac{\alpha \beta_2 I_C^{* *}}{N} & 0 \\ \xi_2^{* *} & 0 & -\frac{\gamma \beta_1 I_T^{* *}}{N} & a_{44} & a_{45} & 0 & S_2^{* *} & 0 \\ 0 & 0 & a_{53} & a_{54} & a_{55} & 0 & \frac{\alpha \beta_2 I_C^{* *}}{N} & 0 \\ 0 & 0 & \delta & 0 & 0 & -h_5 & 0 & 0 \\ 0 & 0 & 0 & \rho & 0 & 0 & -h_6 & 0 \\ 0 & r_1 & t_1 & 0 & \psi & \tau & \sigma & -\mu\end{array}\right]$

where,

$a_{11}=\theta+\xi_1^{* *}+\xi_2^{* *}+\mu, S_1^{* *}=\frac{\beta_1 k_2 \Lambda}{\left(k_1 k_2-\theta \varphi\right) N}$,

$S_2^{* *}=\frac{\beta_2 k_2 \Lambda}{\left(k_1 k_2-\theta \varphi\right) N}, \xi_1^{* *}=\frac{\beta_1\left(I_C^{* *}+I_{C T}^{* *}\right)}{N}$

$\xi_2^{* *}=\frac{\beta_2\left(I_T^{* *}+T^{* *}+I_{C T}^{* *}\right)}{N}, a_{33}=S_1^{* *}-\alpha \xi_2^{* *}-h_2$

$a_{35}=S_1^{* *}-\frac{\alpha \beta_2 I_C^{* *}}{N}+\omega, a_{44}=S_2^{* *}-\gamma \xi_1^{* *}-h_3$

$a_{45}=S_2^{* *}-\frac{\gamma \beta_1 I_T^{* *}}{N}+\eta, a_{53}=\alpha \xi_2^{* *}+\frac{\gamma \beta_1 I_T^{* *}}{N}$

$a_{54}=\gamma \xi_1^{* *}+\frac{\alpha \beta_2 I_C^{* *}}{N}, a_{55}=\frac{\alpha \beta_2 I_C^{* *}}{N}+\frac{\gamma \beta_1 I_T^{* *}}{N}-h_4$

The characteristic equation of $J\left(E_1^*\right)$ is $h(\lambda)=g(\lambda)(\lambda+$ $\left.\tau+d_4+\mu\right)(\lambda+\mu)$, where $g(\lambda)=\left(c_0 \lambda^6+c_1 \lambda^5+c_2 \lambda^4+\right.$ $\left.c_3 \lambda^3+c_4 \lambda^2+c_5 \lambda+c_6\right)$. We get $\lambda_7=-\tau-d_4-\mu<0$, $\lambda_8=-\mu<0$ from polynomial $h(\lambda)$, and $\lambda_i, i=$ $1,2,3,4,5,6$ will be negative if $c_j>0$, where $j=$ $0,1,2,3,4,5,6, \mathcal{R}_0>1, c_1 c_2>c_0 c_3, \quad c_1\left(c_2 c_3+c_0 c_5\right)>$ $c_1^2 c_4+c_0 c_3^2$, and $c_1 c_2 c_4>c_0\left(c_1 c_6+c_2 c_5\right)$. Because the form of the coefficients of the polynomial $g(\lambda)$ is very long, checking the Routh-Hurwitz criterion is accomplished numerically. Using the value of parameters provided in Table 1, we compute the characteristic equation of the Jacobian matrix $J\left(E_2^*\right)$ by employing the Maple software. We get $c_1=$ $0.3091 ; c_2=0.1167 ; c_3=0.0195 ; c_4=0.0013 ; c_5=$ $0.000019 ; c_6=0.00000000083 ; c_1 c_2=0.0955 ;$ $c_0 c_3=0.0195 ;$ $c_1\left(c_2 c_3+c_0 c_5\right)=0.00071 ;$ $c_1^2 c_4+c_0 c_3^2=0.0005$; $c_1 c_2 c_4=0.000048 ; c_0\left(c_1 c_6+c_2 c_5\right)=0.0000023$. It satisfies the condition that the local asymptotically stable endemic equilibrium point exists based on the Routh-Hurwitz criteria [45].

The new model investigated in this study is a development of the work of Mekonen et al. [35] and Inayaturohmat et al. [38] in the form of a dynamic model of co-infection with COVID-19 and tuberculosis expressed in a system of differential equations that takes the vaccination compartment into account. The reproduction number study results indicate that the natural date rate is the proportion of parameters that significantly influence the spread of COVID-19 and TB infection. The re-susceptible rate is the most sensitive of a group of susceptible rate factors in the sensitivity analysis. The immunity rate is the most sensitive of a group of coinfected populations, according to treatment rate criteria. A collection of symptomatic transition rate metrics, the tuberculosis symptomatic transition rate is the most sensitive. The natural death rate is the most vulnerable of the death rate metrics. The COVID-19 recovery rate is the most sensitive, according to the recovery rate criteria for a group of co-infected populations.

Based on the findings of this study, future research could create a mathematical model that takes into consideration the characteristics of recruitment into the infected compartment or reinfection from the recovered compartment in COVID-19 co-infection with tuberculosis. Models can also be created by include exposed compartments or compartments resistant to anti-tuberculosis medications in COVID-19 co-infection with tuberculosis, as well as comorbid dynamic models of COVID-9 disease with tuberculosis and other infectious diseases. Model development can also be accomplished by providing optimal control or by employing a fractional, stochastic differential equation system technique, and the model must be validated using real-world data.