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

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

Jonner Nainggolan* Moch. Fandi Ansori

Department of Mathematics, Faculty of Mathematics and Natural Sciences, Universitas Cenderawasih, Jayapura 99224, Indonesia

Department of Mathematics, Faculty of Science and Mathematics, Universitas Diponegoro, Semarang 50275, Indonesia

Corresponding Author Email:
jonner2766@gmail.com
Page:
34-46
|
DOI:
https://doi.org/10.18280/mmep.110104
22 July 2023
|
Revised:
20 September 2023
|
Accepted:
30 October 2023
|
Available online:
30 January 2024
| Citation

OPEN ACCESS

Abstract:

This research combines the COVID-19 and BCG vaccination subpopulations to examine the spread of COVID-19 coinfection and tuberculosis (TB) using a compartmental mathematical model. The model analysis yields the non-endemic and endemic equilibrium points in addition to the basic reproduction number. The vaccination variable in the model can reduce the incidence of COVID-19, TB, and coinfection. A sensitivity analysis using elasticity index is conducted and the result is that the natural death rate parameter is the most influential in relation to the accelerated spread of COVID-19 co-infection with tuberculosis. Additionally, we conduct a time-dependent sensitivity analysis to determine how varying parameter values influence each subpopulation. By using this technique, we calculate the sensitivity index after reaching equilibrium of several groups of parameters, and the result is that re-susceptible, immunity rate, symptomatic transition rate of TB, COVID-19 recovery rate, and natural death rate are the most influential for each group of parameters on the dynamics of each subpopulation.

Keywords:

infectious disease, COVID-19, tuberculosis, coinfection, mathematical modeling, sensitivity

1. Introduction

COVID-19 is an infectious disease caused by the SARS-CoV-2 virus that first appeared in Wuhan, China, in December 2019 and then spread to Indonesia on March 2, 2020 [1]. For more than three years, the COVID-19 epidemic has been ongoing. As of June 29, 2023, the disease had claimed the lives of 161,871 people in Indonesia and 6951677 people globally. The Indonesian government declared that the country is free of COVID-19 on June 29, 2023. Co-infection with tuberculosis, heart disease, high blood pressure, HIV/AIDS, asthma, cholera, and diabetes causes the majority of COVID-19 case deaths [1, 2].

COVID-19 has been examined by scientists from all around the world, and researchers in the field of mathematics are helping to control the disease. Mathematicians collaborate with health professionals to combat the spread of the COVID-19 infectious disease [1, 3]. To combat the spread of COVID-19 in the Chinese metropolis of Wuhan, mathematical models are being deployed [3]. The basic reproduction number on the dynamics of COVID-19 distribution is analyzed to determine if the spread is rising or decreasing in Nigeria [4, 5]. If the fundamental reproduction number is smaller than one, the dynamic model of COVID-19 in China with recruitment growth rate logistics spread internationally, and an asymptotically stable non-endemic equilibrium point [6]. The COVID-19 mathematical model study also took into account comorbid and non-comorbid [7], symptomatic and asymptomatic disorders in Indonesia [8, 9].

The omicron form of COVID-19 first appeared in Africa in 2022. This variation spreads more quickly than the alpha, beta, and delta versions [1]. To maximize viral control, the parameters of the omicron variant dispersion model were studied [10]. Strategic interventions are carried out to restrict the spread of COVID-19 by researching preventative and therapy criteria [11].

Vaccination against COVID-19 is effective in preventing its spread since it increases antibodies, lowers the likelihood of transmission, and lessens the virus's significant impact [1]. In Indonesia, mathematical models for the transmission of COVID-19 that take into consideration the vaccine compartment for the Susceptible-Exposed-Infected-Recovered (SEIR) type have been examined [12, 13], by incorporating treatment scenarios [14], or by taking immigration characteristics [15]. Efforts to reduce the number of people who tested positive for COVID-19 and increase the number of people who recovered were given preventive controls such as keeping a safe distance, washing hands before and after activities, independent isolation, care, medication, eating multivitamins, and living a healthy lifestyle [1]. Optimal prevention, treatment, and isolation control can reduce the frequency of infections and enhance the number of recoveries in COVID-19 patients [16, 17].

COVID-19 spreads by droplets and free air in a manner similar to tuberculosis [1]. Researchers have extensively explored the use of mathematical models to combat tuberculosis spread, notably authors in references [18, 19], with a focus on immunological variables, demographics, and the transportation infrastructure in Moscow, Russia. The research [20] reviewed the criteria influencing the growth in dissemination and decrease in the number of infected individuals in the Congo. The research [21] also investigated model studies with respect to children's and adult compartments using fractional models. Several articles [17, 22, 23] investigated the stability of the equilibrium point and the basic reproduction number in models of the spread of tuberculosis types SIR, SEIR, BSEIR, and SVEIRE, with B representing the BCG vaccine compartment.

BCG immunization is helpful in preventing tuberculosis infection [2]. Modeling analysis of tuberculosis spread in relation to vaccination compartment has been investigated in studies [17, 24]. Individuals infected with tuberculosis who are not disciplined in taking anti-tuberculosis drugs may die or develop resistance to anti-tuberculosis drugs [2, 25]. vaccine control and treatment are provided to efforts to enhance immunity, tuberculosis dynamics models [26], and vaccine control and treatment [27].

Individuals infected with concomitant COVID-19 and infectious diseases have a higher risk of dying sooner than those who do not have comorbid diseases [1]. The effectiveness of treating COVID-19 comorbidities with rubella was investigated [28], pneumonia [29], and dengue fever [30]. Papers [31, 32] investigated the stability of the equilibrium point and model parameters for the dynamics of tuberculosis co-infection with HIV/AIDS. The paper [33] discussed efforts to accelerate recovery from co-infection with COVID-19 and cholera by providing optimal treatment control for each disease.

Because the pattern of COVID-19 transmission is nearly identical to that of tuberculosis, persons afflicted with tuberculosis are more easily infected with COVID-19 than individuals who are not infected with tuberculosis [1]. Several researchers, notably, Trajman et al. [34] have analyzed models of the propagation of COVID-19 co-infection with tuberculosis to establish the characteristics that have a major effect on the disease's transmission. Using Pearson's coefficients, the COVID-19 co-infection model for omicron and delta variants with HIV and diabetes, as well as the quarantine compartment, produced prediction values that were close to real data.

Mekonen et al. [35] investigated an eight-compartment model of the spread of COVID-19 co-infection with tuberculosis. Marimuthu et al. [36] constructed a model of COVID-19 co-infection with tuberculosis in India to forecast its transmission pattern in the coming years. Lockdown has hampered tuberculosis patient services, resulting in an increase in tuberculosis dissemination and COVID-19 tuberculosis co-infection [37]. In terms of the influence of the COVID-19 pandemic on TB, the analysis shows that people infected with COVID-19 cause an increase in the number of people infected with tuberculosis and slow down the cure rate [38].

Several researchers have explored the regulation of countermeasures for co-infection with tuberculosis in an effort to minimize the spread of COVID-19 co-infection with tuberculosis. In Indonesia, tuberculosis prevention can prevent roughly 27,878,840 new cases, while tuberculosis treatment can control 5,397,795 new cases [39]. Treatment and care control, as well as prevention in the form of physical separation and isolation in a hospital or independently, can reduce the spread of co-infection with COVID-19 and tuberculosis, both infection with each disease COVID-19 and tuberculosis, or co-infection of the two diseases, and prevention costs are more efficient than treatment costs [40, 41].

On the basis of the research of Mekonen et al. [35] and Inayaturohmat et al. [38] that did not include the vaccination factor, in this paper, we propose a new model of the transmission of co-infection COVID-19 and TB incorporating vaccination compartment, since vaccination is effective in preventing the spread of COVID-19 [1, 12-17] or BCG [2, 17, 24-27]. The goals of this study were to: (1) examine the mathematical model for the spread of COVID-19 co-infection with tuberculosis by focusing on the vaccination compartment, (2) determine the equilibrium point, the basic reproduction number, and analyze the stability of the equilibrium point, and (3) analyze the sensitivity of the model parameters, determining the most sensitive parameter for reducing the spread of COVID-19 co-infection with tuberculosis.

2. Mathematical Model

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 IC. 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 IT. Individuals in the IC 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 ICT. 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

 Parameter Description Value Reference $\Lambda$ Recruitment rate 13700 Estimated $\theta$ Vaccination rates from S to V 0.1 [42] $\varphi$ Re-susceptible rates from V to S 0.000685 [39] $r_1$ Immunity rate due to vaccination 0.001644 [39] $\beta_1$ Symptomatic transition rate of COVID-19 5.5´10-6 Assumed $\beta_2$ Symptomatic transition rate of TB 2 ´ 10-6 Assumed $t_1$ Treatment recovery rate of IC 0.33029 [43] $\mu$ Natural death rate of each subpopulation (65×365)-1 Estimated $\delta$ Quarantine rate 0.13266 [43] $\alpha$ Transfer rate from IC to ICT 0.01 Assumed $\omega$ COVID-19 recovery rate of ICT 0.000389 Assumed $d_1$ COVID-19-death rate 0.002 Assumed $\gamma$ Transfer rate from IT to ICT 0.1 Assumed $\eta$ TB cure rate of ICT 0.002335 Assumed $\rho$ Treatment rate of IT 0.1005 Assumed $d_2$ TB-death rate 0.001088 [39] $\psi$ Recovery rate due to treatment of ICT 0.0055 Assumed $d_3$ COVID-19 & TB coinfection death rate 0.003 Assumed $\tau$ Recovery rate due to treatment of Q 0.11624 [39] $d_4$ COVID-19-death rate of Q 0.001 Assumed $\sigma$ Recovery rate due to treatment of T 0.001 [39] $d_5$ TB-death rate of T 0.000039 [39]

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].

3. Numerical Simulation

To validate the previous section’s analytical results, numerical simulations are performed. We start by simulating the solution of system (1) at different time scales. The second section presents the numerical sensitivity of the basic reproduction number. Third, we analyze the sensitivity of the model's parameters to time and determine which parameter is the most sensitive. To accomplish this, we use the parameter values shown in Table 1 and the initial conditions listed below: $S(0)=267137521190$, $V(0)=58200000$, $I_C(0)= 23601$, $I_T(0)=16000000$, $I_{C T}(0)=9000$, $Q(0)=15000$, $T(0)=6400000$, and $R(0)=280$, with the total population is assumed to be $N=268000000$. All the simulations are performed using MATLAB software. To solve the system of differential equations, we use package function ode45.

3.1 Solution of system

Based on the number of subpopulations, we simulate the solution of system (1) into three groups. Various time scales are used to plot the solution. The susceptible (S), vaccinated (V), and recovered (R) populations make up the first subpopulation group. In a short period of time, the susceptible population is decreasing in subgraph (a) of Figure 2, while the vaccinated population is rapidly increasing and then gradually declining, and the recovered population is slowly increasing. Long-term, as depicted in subgraph (d) of Figure 2, the vaccinated population decreases while the recovered population grows until equilibrium is reached. The second group consists of populations infected with COVID-19 (Ic), co-infected with COVID-19 and tuberculosis (Ict), and quarantined (Q). In subgraphs (b) and (e) of Figure 2, these three subpopulations are declining over both short and long periods of time. The final group consists of TB-infected (It) and TB-treated (T) populations. In a short period of time, the population infected with tuberculosis decreases in subgraph (c) of Figure 2, while the population receiving treatment for tuberculosis increases and then declines. As shown in subgraph (f) of Figure 2, the population receiving TB treatment is decreasing.

(a)

(b)

(c)

## 26.png

(d)

(e)

(f)

Figure 2. (a) and (d) Plots of susceptible S, vaccinated V, and recovered R subpopulations in different time scales; (b) and (e) Plots of $I_C, I_{C T}$, and Q subpopulations in different time scales; (c) and (f) Plots of $I_T$ and $T$ subpopulations in different time scales

Table 2. The value of elasticity index of basic reproduction number

 Parameter $\Lambda$ $\theta$ $\varphi$ $r_1$ $\beta_1$ $\beta_2$ $\Psi_q^{\mathcal{R}_0}$ 1 -4.37e-05 1.26e-05 -1.23e-05 0.37 0.62 Parameter $t_1$ $\mu$ $\delta$ $d_1$ $\rho$ $d_2$ $\Psi_q^{\mathcal{R}_0}$ -0.26 -1.0002 -0.107 -0.0016 -0.62 -0.0066

3.2 Sensitivity of the basic reproduction number

The beginning of disease transmission is associated with the basic reproduction number $\mathcal{R}_0$. To identify the parameter that is most sensitive to the spread of disease, we wish to examine every parameter that appears in the $\mathcal{R}_0$ formula. To accomplish this, we define an elasticity index, also known as the normalized sensitivity index, as follows

$\Psi_q^{\mathcal{R}_0}=\frac{\partial \mathcal{R}_0}{\partial q} \times \frac{q}{\mathcal{R}_0}$

where, $q \in\left\{\Lambda, \theta, \varphi, r_1, \beta_1, \beta_2, t_1, \mu, \delta, d_1, \rho, d_2\right\}$. This tells about the relative change of the beginning disease transmission with respect to the changes of the parameters’ value. Therefore, the highest absolute value of the elasticity index corresponds to the most influential parameters to the diseases. We calculate $\Psi_q^{\mathcal{R}_0}$ using parameters’ value in Table 1. The result is given in Table 2 and also depicted in Figure 3. We can observe that the natural death $\mu$ gives highest proportional impact on $\mathcal{R}_0$, followed by the recruitment rate $\Lambda$.

According to the previous section, the basic reproduction number $\mathcal{R}_0$ can be used to indicate the spread of the COVID19 and TB infections. When $\mathcal{R}_0$ is less than 1, there is no disease transmission. In the meantime, when $\mathcal{R}_0$ exceeds 1, diseases become endemic. Now, we wish to determine the effects of parameter changes on the $\mathcal{R}_0$ value, or, in other words, the sensitivity of $\mathcal{R}_0$ is evaluated. The basic reproduction number $\mathcal{R}_0$ is also can be seen as a function of two parameters. Thus, we can assess the impact of combination of two parameters on $\mathcal{R}_0$ by observing its contour plot. Due to the large number of model parameters, we restrict the analysis to the combination of the COVID-19 and TB symptomatic transition rate parameters $\left(\beta_1\right.$ and $\left.\beta_2\right)$ with the other included parameters in (5). Then, the contour of $\mathcal{R}_0$ is plotted. We examine the contour whose $\mathcal{R}_0$ value is greater than 1.

## 3.png

Figure 3. Elasticity index of basic reproduction number for the case of all parameters

Figure 4 depicts the contour plot of $\mathcal{R}_0$ as a function of various parameter combinations. In comparison to the effects that $\beta_1$ and $\beta_2$ have on $\mathcal{R}_0$ value, all parameters other than $\beta_1$ and $\beta_2$ have a negligible or nonexistent effect on $\mathcal{R}_0$. Simulations demonstrate that as $\beta_1$ and $\beta_2$ values increase, diseases become more endemic. Therefore, the COVID-19 and TB symptomatic transition rate parameters are the most influential disease spread parameters, and their values should be minimized as much as possible through the implementation of control or prevention policies. Comparing parameters $\beta_1$ and $\beta_2$, we find that parameter $\beta_2$ has a greater impact on the value of $\mathcal{R}_0$. Consequently, the rate of transition from asymptomatic to symptomatic TB must be the top priority for regulators.

3.3 Sensitivity analysis of parameters

In this subsection, a sensitivity analysis is conducted to determine the effects of parameter changes on the model compartment over time. Due to the large number of parameters, we group them into five categories based on the similarity of their descriptions. Group 1 is for the susceptible rate parameter and includes $\Lambda$ and $\varphi$. Group 2 is for the coinfected population's treatment rate parameter, which consists of $\theta, r_1$, $t_1, \delta, \rho, \tau$, and $\sigma$. Group 3 is for the symptomatic transition rate parameter, and it is made up of $\beta_1, \beta_2, \alpha$, and $\gamma$. Group 4 is for the death rate parameter, which consists of $d_1, d_2, d_3$, $d_4, d_5$, and $\mu$. Finally, group 5 is for the co-infected population's recovery rate parameter, which consists of $\omega, \eta$, and $\psi$.

## 45.png

Figure 4. Sensitivity analysis of the basic reproduction number $R_0$ affected by symptomatic transition rates $\beta_1$ or $\beta_2$ with combinations of other parameters

The steps for the sensitivity analysis are as follows. Suppose the $k$-th group of parameters is denoted $P_k$, where $k=$ $1,2,3,4,5$. Let the vector of compartments is written as $X=$ $\left(S, V, I_C, I_T, I_{C T}, Q, T, R\right)^{\prime}$ and the vector of right-hand side of (1) is written as $F$. We define a sensitivity function $\Omega_k=\frac{\partial X}{\partial P_k}$. By doing a total differentiation $\Omega_k$ of $t$, we get:

$\frac{d \Omega_k}{d t}=\frac{d}{d t} \frac{\partial X}{\partial P_k}=\frac{\partial}{\partial P_k} \frac{d X}{d t}=\frac{\partial F}{\partial X} \frac{\partial X}{\partial P_k}+\frac{\partial F}{\partial P_k}=J(X) \Omega_k+\frac{\partial F}{\partial P_k}$                            (8)

where, $J(X)$ is the Jacobian matrix of system (1) evaluated at point $X$. Matrix $\Omega_k$ has size $8 \times \mathrm{m}$, where $m=\# P_k$, the matrix $J(X)$ has size $8 \times 8$, and the matrix $\frac{\partial F}{\partial P_k}$ has size $8 \times \mathrm{m}$. Thus, the system of differential equations (8) is solved to get a matrix solution $\Omega_k=\left[s_Y^p\right]_{8 \times m}$, with $s_Y^a=\frac{\partial Y}{\partial a}$ where $Y \in X$ and $a \in$ $P_k$. By solving system (8), we calculate the sensitivity index after reaching equilibrium to observe the most influential parameters on the dynamics of all subpopulations.

3.3.1 Group of parameters 1

There are two parameters for the susceptible rate: $\Lambda$, the recruitment rate, and $\varphi$, the re-susceptibility rate in the vaccinated population. As depicted in Figure 5, the recruitment rate benefits the susceptible, vaccinated, and recovered populations. The re-susceptibility rate has a positive effect on the susceptible population but a negative effect on the vaccinated and recovered populations. The remaining populations are unaffected by the two parameters. After equilibrium has been reached, a comparison of these two parameters using the sensitivity index reveals that the re-susceptible parameter is more sensitive than the recruitment rate.

(a)

(b)

## 53.png

(c)

Figure 5. Sensitivity of parameters (a) $\Lambda$ and (b) $\varphi$ to the populations over time, and sensitivity index of these two parameters

3.3.2 Group of parameters 2

There are seven parameters in group 2 which represent the treatment rate parameter: vaccination rate θ, immunity rate $r_1$, treatment recovery rate $t_1$, quarantine rate δ, tuberculosis treatment rate ρ, COVID-19 recovery rate τ, and tuberculosis recovery rate σ. Figure 6 demonstrates that the vaccination rate has positive effects on the vaccinated and recovered populations, but negative effects on the susceptible population. The rate of immunity has a positive effect on the recovered population but a negative effect on the vaccinated population. The treatment recovery rate has a positive effect on the recovered population but a negative effect on the infected and Q populations. The quarantine rate has positive effects on the Q population, but negative effects on the infected and recovered COVID-19 populations. The treatment rate for tuberculosis has a positive effect on T and TB-infected populations. The COVID-19 recovery rate has a positive effect on the recovered population and a negative effect on the Q population. The tuberculosis recovery rate has a positive effect on the recovered population while having a negative effect on the T population. After reaching equilibrium, a comparison of these parameters using the sensitivity index reveals that the immunity rate is the most sensitive of the group.

(a)

(b)

(c)

(d)

(e)

(f)

## 68.jpg

(g)

(h)

Figure 6. Sensitivity of parameters (a) $\theta$, (b) $r_1$, (c) $t_1$, (d) $\delta$, (e) $\rho$, (f) $\tau$, and (g) $\sigma$ to the populations over time, and (h) their sensitivity inde

3.3.3 Group of parameters 3

There are four parameters in group 3: Symptomatic transition rate of COVID-19 $\beta_1$, symptomatic transition rate of TB $\beta_2$, transfer rate from $I_C$ to $I_{C T} \alpha$, and transfer rate from $I_T$ to $I_{C T}$ $\gamma$. In Figure 7, the symptomatic transition rate of COVID-19 has a positive effect on the populations of $I_C, I_{C T}, Q$, and $R$, but a negative effect on the remaining populations. The symptomatic transition rate of tuberculosis has a positive influence on the $I_T$ and $T$ populations, but a negative influence on the $S$ and $V$ populations. The transfer rate $\alpha$ has a positive effect on the $I_{C T}$ and $T$ populations, but a negative effect on the $I_C$ and $R$ populations. The transfer rate $\gamma$ has favorable effects on the $I_{C T}$ and $R$ populations, but negative effects on the $I_T$ and $T$ populations. According to the sensitivity index at equilibrium, the symptomatic transition rate of tuberculosis is the most sensitive parameter among others.x

(a)

(b)

(c)

(d)

## 75.png

(e)

Figure 7. Sensitivity of parameters (a) $\beta_1$, (b) $\beta_2$, (c) $\alpha$, and (d) $\gamma$ to the populations over time, and (e) their sensitivity index

(a)

(b)

(c)

(d)

(e)

(f)

## 87.png

(g)

Figure 8. Sensitivity of parameters (a) $d_1$, (b) $d_2$, (c) $d_3$, (d) $d_4$, (e) $d_5$, and (f) $\mu$ to the populations over time, and (e) their sensitivity index

3.3.4 Group of parameters 4

There are six death rate parameters, namely $d_1, d_2, d_3, d_4, d_5$, and $\mu$. As shown in Figure 8, it is evident that none of these parameters have a positive effect on populations. The parameter $d_1$ has a negative effect on infected, $Q$, and recovered COVID-19 populations. The parameter $d_2$ has a negative impact on populations that are TB-infected, $T$, and recovered. The parameter $d_3$ has a negative impact on populations of TB infected, coinfected, $T$, and recovered. Parameter $d_4$ has a negative effect on both $Q$ and recovered populations. Parameter $d_5$ has a negative effect on the $T$ and recovered populations. All populations are negatively impacted by the parameter $\mu$. According to the sensitivity index diagram, the natural death rate $\mu$ is the most sensitive parameter, followed by parameter $d_5$.

(a)

(b)

(c)

## 94.png

(d)

Figure 9. Sensitivity of parameters (a) $\omega$, (b) $\eta$, and (c) $\psi$ to the populations over time, and (e) their sensitivity index

3.3.5 Group of parameters 5

There are three parameters in group 5: COVID-19 recovery rate of $I_{C T} \omega$, TB cure rate of $I_{C T} \eta$, and recovery rate of $I_{C T}$ due to treatment $\psi$. Figure 9 illustrates that the parameter $\omega$ has a positive effect on the $S, V, Q$, and $R$ populations, but a negative effect on the $I_T, I_{C T}$, and $T$ populations. Parameter $\eta$ has a positive effect on the $I_T$ and $T$ populations, but a negative effect on the $I_C, I_{C T}$, and $R$ populations. Parameter $\psi$ has a positive effect on the $S, V$, and $R$ populations, but a negative effect on the $I_T, I_{C T}$, and $T$ populations. After their sensitivity index has reached equilibrium, we can observe that the parameter $\omega$ has the greatest impact, followed by parameter $\psi$.

4. Conclusions

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.

Acknowledgement

The authors would like to thank Kemenristek Dikti for giving Higher Education Grants of College Leading Basic Research for Fiscal Year 2023, which were funded by LPPM Universitas Cenderawasih.

References

[1] WHO, Guideline Clinical management of COVID-19 patients: Living guideline, 2022. http://apps.who.int/bookorders.

[2] World Health Organization. (2019). WHO guidelines on tuberculosis infection prevention and control: 2019 update (No. WHO/CDS/TB/2019.1). World Health Organization.

[3] Ndaïrou, F., Area, I., Nieto, J.J., Torres, D.F. (2020). Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan. Chaos, Solitons & Fractals, 135: 109846. https://doi.org/10.1016/j.chaos.2020.109846

[4] Iboi, E.A., Sharomi, O., Ngonghala, C.N., Gumel, A.B. (2020). Mathematical modeling and analysis of COVID-19 pandemic in Nigeria. Mathematical Biosciences and Engineering, 17(6): 7192-7220. https://doi.org/10.3934/MBE.2020369

[5] Kifle, Z.S., Obsu, L.L. (2022). Mathematical modeling for COVID-19 transmission dynamics: A case study in Ethiopia. Results in Physics, 34: 105191. https://doi.org/10.1016/j.rinp.2022.105191

[6] Shen, C.Y. (2020). Logistic growth modelling of COVID-19 proliferation in China and its international implications. International Journal of Infectious Diseases, 96: 582-589. https://doi.org/10.1016/j.ijid.2020.04.085

[7] Nainggolan, J., Ansori, M.F. (2022). Stability and sensitivity analysis of the COVID-19 spread with comorbid diseases. Symmetry, 14(11): 2269. https://doi.org/10.3390/sym14112269

[8] Anggriani, N., Ndii, M.Z., Amelia, R., Suryaningrat, W., Pratama, M.A.A. (2022). A mathematical COVID-19 model considering asymptomatic and symptomatic classes with waning immunity. Alexandria Engineering Journal, 61(1): 113-124. https://doi.org/10.1016/j.aej.2021.04.104

[9] Ahmed, I., Modu, G.U., Yusuf, A., Kumam, P., Yusuf, I. (2021). A mathematical model of Coronavirus Disease (COVID-19) containing asymptomatic and symptomatic classes. Results in Physics, 21: 103776. https://doi.org/10.1016/j.rinp.2020.103776

[10] Khan, M.A., Atangana, A. (2022). Mathematical modeling and analysis of COVID-19: A study of new variant Omicron. Physica A: Statistical Mechanics and its Applications, 599: 127452. https://doi.org/10.1016/j.physa.2022.127452

[11] Teklu, S.W. (2022). Mathematical analysis of the transmission dynamics of COVID-19 infection in the presence of intervention strategies. Journal of Biological Dynamics, 16(1): 640-664. https://doi.org/10.1080/17513758.2022.2111469

[12] Annas, S., Pratama, M.I., Rifandi, M., Sanusi, W., Side, S. (2020). Stability analysis and numerical simulation of SEIR model for pandemic COVID-19 spread in Indonesia. Chaos, Solitons & Fractals, 139: 110072. https://doi.org/10.1016/j.chaos.2020.110072

[13] Nuraini, N., Sukandar, K.K., Hadisoemarto, P., Susanto, H., Hasan, A.I., Sumarti, N. (2021). Mathematical models for assessing vaccination scenarios in several provinces in Indonesia. Infectious Disease Modelling, 6: 1236-1258. https://doi.org/10.1016/j.idm.2021.09.002

[14] Diagne, M.L., Rwezaura, H., Tchoumi, S.Y., Tchuenche, J.M. (2021). A mathematical model of COVID-19 with vaccination and treatment. Computational and Mathematical Methods in Medicine, 2021: 1250129. https://doi.org/10.1155/2021/1250129

[15] Semlali, M., Hattaf, K., El Kettani, M.E. (2022). Modeling and analysis of the dynamics of COVID-19 transmission in presence of immigration and vaccination. Communications in Mathematical Biology and Neuroscience, 2022: 36. https://doi.org/10.28919/cmbn/7270

[16] Nainggolan, J., Harianto, J., Tasman, H. (2023). An optimal control of prevention and treatment of COVID-19 spread in Indonesia. Communications in Mathematical Biology and Neuroscience, 2023: 3. https://doi.org/10.28919/cmbn/7820

[17] Sulayman, F., Abdullah, F.A., Mohd, M.H. (2021). An SVEIRE model of tuberculosis to assess the effect of an imperfect vaccine and other exogenous factors. Mathematics, 9(4): 327. https://doi.org/10.3390/math9040327

[18] Castillo-Chavez, C., Song, B. (2004). Dynamical models of tuberculosis and their applications. Mathematical Biosciences, 1(2): 361-404.

[19] Avilov, K.K., Romanyukha, A.A., Belilovsky, E.M., Borisov, S.E. (2022). Mathematical modelling of the progression of active tuberculosis: Insights from fluorography data. Infectious Disease Modelling, 7(3): 374-386. https://doi.org/10.1016/j.idm.2022.06.007

[20] Kasereka Kabunga, S., Doungmo Goufo, E.F., Ho Tuong, V. (2020). Analysis and simulation of a mathematical model of tuberculosis transmission in Democratic Republic of the Congo. Advances in Difference Equations, 2020(1): 642. https://doi.org/10.1186/s13662-020-03091-0

[21] Fatmawati, M.A.K., Bonyah, E., Hammouch, Z., Shaiful, E.M. (2020). A mathematical model of tuberculosis (TB) transmission with children and adults groups: A fractional model. AIMS Mathematics, 5(4): 2813-2842. http://doi.org/10.3934/math.2020181

[22] Ucakan, Y., Gulen, S., Koklu, K. (2021). Analysing of tuberculosis in Turkey through SIR, SEIR and BSEIR mathematical models. Mathematical and Computer Modelling of Dynamical Systems, 27(1): 179-202. https://doi.org/10.1080/13873954.2021.1881560

[23] Das, K., Murthy, B.S.N., Samad, S.A., Biswas, M.H.A. (2021). Mathematical transmission analysis of SEIR tuberculosis disease model. Sensors International, 2: 100120. https://doi.org/10.1016/j.sintl.2021.100120

[24] Nkamba, L.N., Manga, T.T., Agouanet, F., Mann Manyombe, M.L. (2019). Mathematical model to assess vaccination and effective contact rate impact in the spread of tuberculosis. Journal of Biological Dynamics, 13(1): 26-42. https://doi.org/10.1080/17513758.2018.1563218

[25] Liu, L., Wang, Y. (2014). A mathematical study of a TB model with treatment interruptions and two latent periods. Computational and Mathematical Methods in Medicine, 2014: 932186. https://doi.org/10.1155/2014/932186

[26] Ojo, M.M., Peter, O.J., Goufo, E.F.D., Panigoro, H.S., Oguntolu, F.A. (2023). Mathematical model for control of tuberculosis epidemiology. Journal of Applied Mathematics and Computing, 69(1): 69-87. https://doi.org/10.1007/s12190-022-01734-x

[27] Gao, D.P., Huang, N.J. (2018). Optimal control analysis of a tuberculosis model. Applied Mathematical Modelling, 58: 47-64. https://doi.org/10.1016/j.apm.2017.12.027

[28] Artiono, R., Wintarti, A., Prawoto, B.P., Astuti, Y.P. (2022). Co-infection model for Covid-19 and Rubella with vaccination treatment: stability and threshold. Communications in Mathematical Biology and Neuroscience, 2022: 8. https://doi.org/10.28919/cmbn/6973

[29] Pérez, A.G.C., Oluyori, D.A. (2022). A model for COVID-19 and bacterial pneumonia coinfection with community-and hospital-acquired infections. Mathematical Modelling and Numerical Simulation with Applications, 2(4): 197-210. https://doi.org/10.53391/mmnsa.2022.016

[30] Omame, A., Rwezaura, H., Diagne, M.L., Inyama, S.C., Tchuenche, J.M. (2021). COVID-19 and dengue co-infection in Brazil: optimal control and cost-effectiveness analysis. The European Physical Journal Plus, 136(10): 1090. https://doi.org/10.1140/epjp/s13360-021-02030-6

[31] Sharomi, O., Podder, C.N., Gumel, A.B., Song, B. (2007). Mathematical analysis of the transmission dynamics of HIV/TB coinfection in the presence of treatment. Mathematical Biosciences & Engineering, 5(1): 145-174. https://doi.org/10.3934/mbe.2008.5.145

[32] Shah, N.H., Sheoran, N., Shah, Y. (2020). Dynamics of HIV-TB co-infection model. Axioms, 9(1): 29. https://doi.org/10.3390/axioms9010029

[33] Hezam, I.M., Foul, A., Alrasheedi, A. (2021). A dynamic optimal control model for COVID-19 and cholera co-infection in Yemen. Advances in Difference Equations, 2021(1): 1-30. https://doi.org/10.1186/s13662-021-03271-6

[34] Trajman, A., Felker, I., Alves, L.C., Coutinho, I., Osman, M., Meehan, S.A., Singh, U.B., Schwartz, Y. (2022). The COVID-19 and TB syndemic: The way forward. The International Journal of Tuberculosis and Lung Disease, 26(8): 710-719. https://doi.org/10.5588/ijtld.22.0006

[35] Mekonen, K.G., Balcha, S.F., Obsu, L.L., Hassen, A. (2022). Mathematical modeling and analysis of TB and COVID-19 coinfection. Journal of Applied Mathematics, 2022: 1-20. https://doi.org/10.1155/2022/2449710

[36] Marimuthu, Y., Nagappa, B., Sharma, N., Basu, S., Chopra, K.K. (2020). COVID-19 and tuberculosis: A mathematical model based forecasting in Delhi, India. Indian Journal of Tuberculosis, 67(2): 177-181. https://doi.org/10.1016/j.ijtb.2020.05.006

[37] Cilloni, L., Fu, H., Vesga, J.F., Dowdy, D., Pretorius, C., Ahmedov, S., Nair, S.A., Mosneaga, A., Masini, E., Sahu, S., Arinaminpathy, N. (2020). The potential impact of the COVID-19 pandemic on the tuberculosis epidemic a modelling analysis. EClinicalMedicine, 28: 1-9. https://doi.org/10.1016/j.eclinm.2020.100603

[38] Inayaturohmat, F., Anggriani, N., Supriatna, A.K. (2022). A mathematical model of tuberculosis and COVID-19 coinfection with the effect of isolation and treatment. Frontiers in Applied Mathematics and Statistics, 8: 958081. https://doi.org/10.3389/fams.2022.958081

[39] Rwezaura, H., Diagne, M.L., Omame, A., de Espindola, A.L., Tchuenche, J.M. (2022). Mathematical modeling and optimal control of SARS-CoV-2 and tuberculosis co-infection: A case study of Indonesia. Modeling Earth Systems and Environment, 8(4): 5493-5520. https://doi.org/10.1007/s40808-022-01430-6

[40] Bandekar, S.R., Ghosh, M. (2022). A co-infection model on TB-COVID-19 with optimal control and sensitivity analysis. Mathematics and Computers in Simulation, 200: 1-31. https://doi.org/10.1016/j.matcom.2022.04.001

[41] Goudiaby, M.S., Gning, L.D., Diagne, M.L., Dia, B.M., Rwezaura, H., Tchuenche, J.M. (2022). Optimal control analysis of a COVID-19 and tuberculosis co-dynamics model. Informatics in Medicine Unlocked, 28: 100849. https://doi.org/10.1016/j.imu.2022.100849

[42] Kementerian Kesehatan Republik Indonesia. (2020). Pedoman Kesiapsiagaan Menghadapi Coronavirus Disease (COVID-19). Kementerian Kesehatan RI Direktorat Jenderal Pencegahan dan Pengendalian Penyakit (P2P), pp. 11-12, 46-55.

[43] Liu, S., Bi, Y., Liu, Y. (2020). Modeling and dynamic analysis of tuberculosis in mainland China from 1998 to 2017: The effect of DOTS strategy and further control. Theoretical Biology and Medical Modelling, 17(1): 6. https://doi.org/10.1186/s12976-020-00124-9

[44] Van den Driessche, P., Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1-2): 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6

[45] Hauksdóttir, A.S., Sigurðsson, S.Þ. (2020). Proving Routh’s theorem using the Euclidean algorithm and Cauchy’s theorem. IFAC-PapersOnLine, 53(2): 4460-4467. https://doi.org/10.1016/j.ifacol.2020.12.446