Mathematical Modeling and Analysis of TB and COVID-19 Coinfection

Tuberculosis (TB) and coronavirus (COVID-19) are both infectious diseases that globally continue a ﬀ ecting millions of people every year. They have similar symptoms such as cough, fever, and di ﬃ culty breathing but di ﬀ er in incubation periods. This paper introduces a mathematical model for the transmission dynamics of TB and COVID-19 coinfection using a system of nonlinear ordinary di ﬀ erential equations. The well-posedness of the proposed coinfection model is then analytically studied by showing properties such as the existence, boundedness, and positivity of the solutions. The stability analysis of the equilibrium points of submodels is also discussed separately after computing the basic reproduction numbers. In each case, the disease-free equilibrium points of the submodels are proved to be both locally and globally stable if the reproduction numbers are less than unity. Besides, the coinfection disease-free equilibrium point is proved to be conditionally stable. The sensitivity and bifurcation analysis are also studied. Di ﬀ erent simulation cases were performed to supplement the analytical results.


Introduction
The outbreak of coronavirus  was first reported by the World Health Organization (WHO) as pneumonia of unknown cause on December 31, 2019, observed in Wuhan, China. However, on the 10 th of January 2020, the virus was determined to be the new family of novel coronavirus with the same category as the severe acute respiratory syndrome virus (SARS-CoV) and the Middle Eastern Respiratory Syndrome virus (MERS-CoV) that occurred, respectively, in Asia and Saudi Arabia [1]. COVID-19 is a respiratory virus that transmits mainly through droplets of saliva generated from an infected person through coughing or sneezing. The rapid expansion of the virus forced the WHO to declare it as a pandemic on March 11, 2020 [2]. The WHO report of COVID-19 indicated that on 10 August 2021, more than 202.14 million people were infected with the virus, and over 4.28 million have died. Different variants of COVID-19 occur in different countries, namely, α, β, γ, and δ. The δ variant is the highest contagion, which is 40 − 60% more transmissible than the α variant and with a higher risk of hospitalization [3].
Tuberculosis (TB) is a communicable disease which is caused by Mycobacterium tuberculosis. It is a major cause of death from a single infectious agent. The bacteria that cause TB are usually spread when an infected person coughs, sneezes, or any other forceful expiratory maneuver that shears respiratory secretions. It mostly affects the lungs but can also affect other organs such as bones, brain, kidneys, and glands [4]. About a quarter of the worldwide population is infected with Mycobacterium tuberculosis and thus at risk of developing TB disease. Globally, approximately 10 million people fell ill with TB, and 1.45 million people died from TB in 2018 [5].
Both TB and COVID-19 are infectious diseases that are transmitted mainly via close contact. The growing clinical evidence suggests that TB diseases are associated with COVID-19 outcomes, including an approximately two to three-fold increase in mortality and a 25% relative decrease in the possibility of recovery for COVID-19 coinfection with current TB disease [6]. Studies on COVID-19 and TB coinfection are limited. It was forecast that people infected with both TB and COVID-19 may have poorer treatment outcomes [7]. TB patients should take precautions as suggested by health professionals to take care of COVID-19 and continue their TB treatment as directed. People infected with COVID-19 and TB show similar symptoms like cough, fever, and difficulty breathing [7,8]. The coinfection of the two diseases exists when individuals are infected with both the COVID-19 and TB at the same time.
Modeling the transmission dynamics of diseases is a way to formulate what is known about the transmissions and explore all possible features of the disease with mathematical techniques [9]. The dynamical model study is used to investigate the vital parameters, predict future trends, and evaluate control measures to provide decisive information for decision making [10]. Recently, many mathematical models have been developed to study the transmission dynamics of COVID-19, see for example [11][12][13][14][15][16] and references cited therein. Mathematical models for TB are also extensively studied in [17][18][19][20][21] and references cited therein. The studies in [22] described the effect of speeding the detection time of COVID-19 infected individuals and reducing the transmission rate using both deterministic and stochastic models. In [23], the authors developed a mathematical to describe the dynamics of COVID-19 by taking into account all modes of transmission and prosocial awareness. Moreover, an optimal control strategy was proposed in [24] that helps to reduce the contact between TB infectious and exposed individuals.
The coinfection of various diseases is widely studied in the mathematical literature, for instance, TB and HIV/AIDS coinfection [4,25], Hepatitis E and HIV [26], and choleraschistosomiasis coinfection in [27]. On the contrary, the study on COVID-19 and other coinfection are at its infant stage, and related mathematical literature is scant. For instance, COVID-19 and malaria coinfection are discussed in [28,29], while TB and COVID-19 coinfection presented in [30,31]. The studies in [30] attempted to estimate the number of TB patients co-infected with COVID-19 with and without public health interventions in New Delhi, India. Their study showed that the peak of COVID-19 and TB coinfected patients would occur on the 94 th day in the absence of public health interventions and on the 138 th day with interventions. Besides, the coinfection of the MERS-CoV and TB model is analyzed in [31]. In their model, the authors studied the bifurcation analysis of the model and performed a sensitivity analysis of the basic reproductive number for the model parameters. To the best of our knowledge, mathematical models considering the coinfection of TB and COVID-19 are rare and they prompt us for further investigation.
In this work, we propose a new mathematical model for TB and COVID-19 coinfection by classifying the total population into eight compartments. The paper is organized as follows. A detailed introduction to both diseases has been discussed in Section 1. The model formulation of COVID-19 and TB coinfection is investigated in Section 2. The analysis of the full model is given in Section 3, and the mathematical results of the submodels are also investigated. In Section 4, numerical results are obtained and discussed with detailed for each submodel as well as the coinfection model, while our work is concluded in Section 5.

Model Formulation
In this section, we formulate the governing mathematical model for the codynamics of TB and COVID-19 by dividing the total population into eight mutually exclusive compartments denoted as susceptible individuals (S), latent TB infected individuals who have no symptoms of TB disease and are not infectious (L T ), active TB-infected individuals (I T ), COVID-19 infected individuals with no symptoms but are infectious (E C ), COVID-19 infected individuals with clinical symptoms and are infectious (I C ), both latent TB and COVID-19 coinfected individuals (L TC ), active TB and COVID-19 coinfected individuals (I TC ), and recovered population from both TB and COVID-19 (R). With these in mind, the total population at time t, denoted by NðtÞ, is given by We assumed that the susceptible population is increased by the recruitment at a rate Λ. All population in each compartment suffer from natural death rate μ. Susceptible individuals acquire TB through contact with active TB patients by a force of infection λ T as expressed in equation (3). In this expression, β 1 denotes the effective contact rate for TB infection. The latent TB infected individuals are assumed to be asymptomatic that does not transmit the disease [32]. Similarly, susceptible individuals acquire COVID-19 infection following effective contact with people infected with COVID-19 at a force of infection λ C defined as in equation (4). Here, β 2 represents the effective contact rate for COVID-19 disease. Moreover, we assumed that individuals in latent TB compartment (L T ) leave to compartment ðI T Þ at a rate α, and to latent TB and COVID-19 coinfected class at a force of infection ηλ C and some portion recovered at a treatment rate ω. Besides, Individuals in TB-infected class (I T ) become recovered from the disease at a rate γ while the remaining portion transfers to the coinfection compartment ðI TC Þ at a rate θ or die due to TB-induced death at a rate δ T .
The population in compartment (L TC ) either progress to coinfected class ðI TC Þ at a rate ρ or die with COVID-19induced death rate δ C . The remaining individuals are assumed to be transferred to either compartment at a constant multiple of σ as illustrated in Figure 1. That is, the population in class L TC moves to class I T with a rate of qσ, to compartment I C class rate of pσ, and become recovered at a rate of ð1 − ðp + qÞÞσ. Furthermore, we assumed that individuals in compartment I TC leaves to compartments I T , I C , or R, respectively, at a rate of nτ, mτ, or ð1 − ðm + nÞÞτ while 2 Journal of Applied Mathematics coinfection-induced death rate is denoted by δ TC . In addition, the COVID-19 exposed individual (E C ) has the chance to leave either to compartment L TC , I C , or R, respectively, at a rate of ϵλ T , φ, or π. Similarly, the population in compartment (I C ) is recovered at a constant rate of ψ or transferred to co-infection compartment at a rate of ν. The disease induced death rate in this compartment is denoted by δ c . The flow diagram of the proposed model is illustrated in Figure 1. Based on the flow diagram, it results in systems of the following nonlinear differential equations: where and with nonnegative initial conditions Sð0Þ > 0, E T ð0Þ ≥ 0, I T ð0Þ ≥ 0, E C ð0Þ ≥ 0, I C ð0Þ ≥ 0, L TC ð0Þ ≥ 0, I TC ð0Þ ≥ 0, and Rð0Þ ≥ 0. The model parameters are described in Table 1.

Positivity and Boundedness of Solutions.
As we are dealing with human populations, all solutions must be positive and bonded in a feasible region. To assure these, we have the following theorem.
Theorem 1. The solutions of the system of equation (2) are positive, unique, and bounded in the region Proof. All functions on the right-hand side of equation (2) are C 1 on ℝ 8 + . Thus, by the Picard-Lindel€ of theorem [33], equations (2) has a unique solution. Let the system of equations in (2) be written in the form x′ = f ðx, tÞ, where x = ðS , L T , I T , E C , I C , L TC , I TC , RÞ, and f is the right hand sides of the equation (1). Following the results of proposition A.1 in [34], the functions f ðx, tÞ has the property of whenever x ∈ ½0,∞Þ 8 andx i = 0. As there exists a unique solution for the system of equations (2), it follows that xðtÞ ∈ ½0,∞Þ 8 for all t ≥ 0 whenever xð0Þ ≥ 0. The change of total population NðtÞ = S + L T + I T + E C + I C + L TC + I TC + R is governed by The solution for this linear first order ode is NðtÞ ≤ Nð 0Þe −μt + ðΛ/μÞð1 − e −μt Þ. Thus, for the initial data 0 ≤ Nð0Þ, we obtain 0 ≤ NðtÞ ≤ Λ/μ: Hence, the solutions of the system of equations in (2) exists, unique, and bounded in a feasible region Ω, which concludes the proof.

Model Analysis
To understand the dynamics of the proposed model, we find the equilibrium points of the system and investigate the dynamics of the equilibrium points. The analysis will be done by investigating the behavior of the submodels for TB and COVID-19 as well as the coinfection model.

COVID-19
Submodel. Without considering the infections of people with tuberculosis, the COVID-19 submodel is given by 3.1.1. Local Stability of Equilibrium Points. The equilibrium point of the COVID-19 submodel is a solution to the system of equations: The disease-free equilibrium point is obtained by putting E C = 0 and I C = 0 in equation (9). Hence, the disease-free equilibrium point, denoted by E C 0 , becomes E C 0 = ððΛ/μÞ, 0, 0, 0Þ.
The basic reproduction number of the submodel defined as the average number of secondary infections produced by a single COVID-19 infected individual in a completely susceptible population. It is computed using the next generation Rate of TB exposed individuals become exposed with COVID-19 γ Recovery rate of TB infected individuals δ T Death rate due to TB bacteria δ C Death rate due to the COVID-19 virus δ TC Death rate due to the coinfection of both diseases μ Natural death rate of the individuals ρ Transfer rate of TB and COVID-19 exposed individuals to the co-infected class φ Transfer rate of COVID-19 exposed individuals to become infected ψ Recovery rate of COVID-19 infected individual ϵ Rate of COVID-19 exposed individuals to become exposed with TB ω Recovery rate of TB latent infected individuals 4 Journal of Applied Mathematics matrix as given in [35]. Let F i ðt, xÞ be the input rate of newly infected individuals in the i th compartment and V + i ð t, xÞ and V − i ðt, xÞ, respectively, be the input rates by other means and rate of transfer of individuals out of a compartment i, where the x i 's are our state variables and V i ðt, xÞ = V − i ðt, xÞ − V + i ðt, xÞ: Hence, from equation (8), it follows that The Jacobian matrices of FðxÞ and V ðxÞ, respectively, are Hence, the next-generation matrix FV −1 becomes The eigenvalues of the next generation matrix are found to be λ 1 = ðβ 2 /μ + φ + πÞð1 + ðφ/μ + δ c + ψÞÞ, and λ 2 = 0. The largest spectral radius of the next generation matrix denotes the basic reproduction number for the submodel. Thus, we have the following result.

Theorem 2. The basic reproduction number of the COVID-19 sub-submodel is
We use the Jacobean matrix to determine the local stability of the equilibrium points. For the submodel (8), the Jacobean matrix becomes The Jacobean matrix of the submodel at the disease-free equilibrium point E c 0 is given as Here, the eigenvalues for JðE C 0 Þ are λ 1,2 = −μ, and the remaining eigenvalues can be easily obtained from the submatrix: In order to determine the local stability of the disease-free equilibrium point, it is sufficient to show that the trace of J 2 is less than zero, and its determinant is greater than zero.

Journal of Applied Mathematics
The trace of J 2 ,trJ 2 This value is greater than zero if β 2 ðμ + δ C + ψ + φÞ/ðμ The endemic equilibrium point of the COVID-19 submodel is the solution of the system of equation in (18).
Now, solving this system of equations in terms λ * C = ð   Journal of Applied Mathematics β 2 /N * ÞðE * C + I * C Þ one can obtain the following expressions for the endemic equilibrium point Here, This implies that the force of infection λ * C is positive at endemic equilibrium point only if R C 0 > 1. Thus, we have just proved the following theorem.
Theorem 3. If R c 0 > 1, then the submodel system (8) has a unique endemic equilibrium point.

Global Stability of the Endemic Equilibrium Point.
The global stability analysis of the endemic equilibrium point Σ C is analyzed using Lyapunov function. For this, we define the following: The function L is greater than zero, and it is equal to zero at the endemic equilibrium point Σ C . Differentiating the function with respect to time, we obtain This holds true for R C 0 > 1 which corresponds the existence of Σ C . Thus, dL/dt < 0 implies that the function is a strictly Lyapunov function which indicates that the endemic equilibrium point Σ C is globally asymptotically stable. Epidemically, this result implies that COVID-19 will continue to exist in the human population for a long duration.
3.1.3. Sensitivity Analysis of the Submodel. In this section, we perform the sensitivity analysis of the model parameters for the COVID-19 submodel given in (8). The sensitivity of a parameter p reflects how the model behavior responds to a small change in a parameter value, and it is defined as [35] In our context, the sensitivity analysis of each parameters for equation (8) becomes Here, we can observe that the contact rate β 2 positively influences the spread of COVID-19. Further, if δ C + ψ − π < 0, the transfer rate φ from exposed compartment to infected class has a positive impact on the spread of the virus. That is, increasing the value of these parameters will increase the endemic. The remaining parameters, μ, ψ, δ C , and π, have negative impact, which means increasing the value of these parameters will decrease the number of people infected with COVID-19. However, in the study of sensitivity analysis, it is unethically increasing human mortality rate to control disease epidemic and hence do not considered.   Submodel. Next, we study the solution behavior of equations (8) by taking β 2 as the bifurcation parameter. Calculating the value of β 2 from R C 0 = 1, i.e., ðβ 2 /μ + φ + πÞðμ + δ C + ψ + φ/μ + δ C + ψÞ = 1, we have β * 2 = ðμ + φ + πÞðμ + δ C + ψÞ/μ + δ C + ψ + φ. Following the result given in Theorem 4 [36], we can calculate the eigenvalues of the Jacobin matrix at the disease-free equilibrium point by substituting β * 2 . Thus, substituting β 2 = β * 2 in equations (15), it gives zero eigenvalue. This means that the Jacobean matrix JðE c 0 Þ in equation (15) at β 2 = β * 2 has a left eigenvector (associated with the zero eigenvalue) Now, assume that f k represents the right-hand side of the k th equation in the COVID-19 submodel (8) and let x k denote the corresponding state variable for k = 1, ⋯, 4.

Define
The local dynamics of equation (8)   For the submodel (8), the associated nonzero partial derivatives at the disease-free equilibrium E C 0 are Then, using the above expressions for a, it follows that For the sign of b, it can be shown that the associated nonvanishing partial derivatives are Clearly, it also follows from the above expressions that   Journal of Applied Mathematics From the computations of a and b, we observe that b is always positive, and the sign of the bifurcation coefficient a depends on the minimum value of reinfection that induces bistability (κ). Thus, using Theorem 4 of [36], the following result is proved.

TB Submodel.
In the absence of COVID-19 disease, equation (2) automatically reduced to TB-submodel expressed as In a similar argument, the disease-free equilibrium point of equation (32) becomes E T 0 = ððΛ/μÞ, 0, 0, 0Þ, and the corresponding basic reproduction number is given by Theorem 5. The disease-free equilibrium point of the TB submodel (32) is locally asymptotically stable if R T 0 < 1.
Proof. The Jacobean matrix at the disease-free equilibrium point for the TB submodel is given by Clearly, the two eigenvalues are equal to −μ, while the remaining can be obtained from the reduced matrix Here, the trace of the matrix is trðJ 3 Þ = −ð2μ + α + ω + δ T + γÞ < 0, and its determinant is det ðJ 3 Þ = ðμ + α + ωÞðμ + δ T + γÞ − αβ 1 . Thus, the determinant is greater than zero if R T 0 < 1. Therefore, the disease-free equilibrium point of the TB submodel is locally asymptotically stable if R T 0 < 1. This completes the proof.
The endemic equilibrium point of the TB-submodel is obtained by solving the system of equations: where λ T = β 1 I T /N. Solving this equation, we have The simplified expression for λ T in terms of S * , L * T , I * T , and R * becomes Note that the endemic equilibrium point exists if and only if R T 0 > 1, since k 1 k 2 − αδ T > 0.

Journal of Applied Mathematics
Proof. To proof this theorem, we define the following function: Next, we show that this function V is a Lyapunov function. Differentiating the function V with respect to time t, we For R T 0 > 1, the endemic equilibrium point exists. Thus, dV/dt < 0 implies that the function V is strictly Lyapunov function which deduces that the endemic equilibrium point is globally asymptotically stable.  12 Journal of Applied Mathematics system of equation (41): where λ T and λ C are as given in equations (3), and (4), respectively. The disease-free equilibrium point (E 0 ) of the system becomes In the following, we calculate the basic reproduction number (R 0 ) of the coinfection model using the method of next generation matrix [37]. For this, we use the notation X = ðL T , I T , E C , I C , L TC , I TC Þ to denote the vector functions: representing the appearance of new infections and the transfer of individuals into and out of the infected compartments,   Table 3.  Table 3. where The largest spectral radius of the next generation matrix denotes the basic reproduction number of the coinfection model. Clearly, four of the eigenvalues of the matrix FV −1 are equal to zero. The remaining eigenvalues are obtained from the reduced matrix : Thus, computing the eigenvalues of FV −1 one can easily obtain that where p . Hence, the basic reproduction number R 0 of the coinfection model (1) is given by 14 Journal of Applied Mathematics equation (2) at disease-free equilibrium point E 0 represented as The local asymptotic stability of E 0 is analyzed based on the sign of the eigenvalues. Expanding the characteristic equation |λI 8 − JðE 0 Þ | = 0 with the first and eighth columns, and the sixth row, we obtain three eigenvalues λ 1,2 = −μ and λ 3 = −c 5 . Further, we obtain the remaining five eigenvalues from the characteristic equation of the reduced matrix: Here, the characteristic polynomial is computed using 15 Journal of Applied Mathematics elementary row operations method [38]. Thus, the row reduced form of the matrix becomes Hence, the characteristic polynomial of our matrix is the product of the diagonal entries times ð−1Þ k , where k is the number of row swaps. We perform three row swaps, so that the characteristic polynomial is Simplifying it gives That is, the eigenvalues are the solutions of the characteristic equation: where Thus, using the Routh-Hurwitz stability criterion, the roots of the equation (54) have negative real parts if the fol-lowing condition holds.
Hence, we have proved the following result.
Theorem 7. The disease-free equilibrium point of the coinfection model (2) is locally asymptotically stable if the condition given in equation (56) where X = ðS, RÞ and Z = ðL T , I T , E C , I C , L TC , I TC Þ. Here, the components of X ∈ ℝ 2 + denotes the number of uninfected individuals while Z ∈ ℝ 6 + represents the number of infected ones [39]. The disease-free equilibrium is denoted now as The conditions ðH 1 Þ and ðH 2 Þ below must be met to guarantee global asymptotically stability: (H 1 ) For dX/dt = FðX, 0Þ, the equilibrium point U 0 is globally stable (H 2 ) GðX, ZÞ = AZ −ĜðX, ZÞ,ĜðX, ZÞ ≥ 0 for ðX, ZÞ ∈ Ω , where A = D Z GðU 0 , 0Þ is a Metzler matrix (the nondiagonal entities of A are nonnegative), and Ω is the biological feasible region of the model From the coinfection model (2), we have Thus,Ĝ 5 ðX, ZÞ = −ηλ C L T − ϵλ T E C < 0. This shows that condition (H 2 ) is not satisfied. As a result, U 0 and then the disease-free equilibrium point E 0 cannot be globally asymptotically stable.

Impacts of TB on COVID-19.
To analyze the impact of TB disease on COVID-19 (and vice versa), we began by expressing the basic reproduction numbers R C 0 in terms of R T 0 (and vice versa) [40]. Expressing the parameter μ in the equation (33) in terms of R T 0 , we have Solving for μ and simplifying, we obtain Substituting this value of μ in R C 0 (13), we obtain R C 0 in terms of R T 0 as where h 1 = R T 0 ððδ C + φ + ψÞ − k 1 + k 2 /2Þ, and h 2 = R T 0 ð ðk 1 − k 2 Þ 2 R T 0 + 4αβ 1 Þ. The partial derivative of R C 0 with respect to R T 0 is then given by with h 3 = ðφ 2 + ðπ + δ C + ψÞφ + ðδ C + ψÞ 2 ÞðR T 0 Þ 2 , and Here, we observe that ∂R C 0 /∂R T 0 is positive. This implies that the expansion of TB infection exacerbates the pandemic of COVID-19.
Remark 8. If ∂R C 0 /∂R T 0 = 0, the expansion of TB in the community has no significant impact on the spread of COVID-19. In the contrary, when ∂R C 0 /∂R T 0 < 0, the increase in TB epidemic will negatively influence the spread of COVID-19.
Similarly, the impact of COVID-19 for TB cases can be addressed by expressing R T 0 in terms of R c 0 and identifying the sign of the partial derivative of R T 0 with respect to R c 0 .

Numerical Simulation
In the previous sections, we have discussed the analytical behaviors of the model equation (2) and the submodels. In the following, we perform numerical solutions of submodels as well as the coinfection model. The solutions of the system equations are integrated using the ode45 solver in MATLAB.
The parameter values are also described in respective subsections.

Numerical Solutions of the COVID-19
Submodel. The analytical findings of the submodel (8) are illustrated in section 3.1. We now present the numerical solutions using initial conditions for state variables Sð0Þ = 5000, E C ð0Þ = 1000, I C ð0Þ = 10, Rð0Þ = 1, and the parameter values described in each figure. Figure 2 presents the solution of the COVID-19 submodel using the parameter values described in the caption. The basic reproduction number is obtained as R C 0 = 0:2524, which implies the disease free equilibrium point is stable, and no endemic equilibrium point exists. Using different initial conditions for each compartment, all solution curves converge to the disease-free equilibrium point as seen in Figure 3. In Figure 2, the susceptible populations approach to the value of Λ/μ, while the other compartments go to zero, which supports our analytical findings in the previous sections. Figure 4 illustrates the time serious plot of the COVID-19 submodel, for which its basic reproduction number is greater than unity. In the analytical findings, we have shown that the disease-free equilibrium point is unstable, and the endemic equilibrium point is stable. Using different initial conditions for each compartment, all solution curves converge to the endemic equilibrium point as seen for the susceptible, infected, and recovered plots of Figure 5. As time goes, in Figures 5(a)-5(c), the solutions close to the endemic equilibrium point given in equation (19) which agrees with the analytical properties. Figures 6 and 7 show the effects of the COVID-19 transmission rate and the recovery rate of the COVID-19 infected individuals. The other parameter values used in the simulation of the figures are Λ = 500, μ = 0:0477, φ = 0:096, δ C = 0:0023, ω = 0:09, and in Figure 6, β 2 = 0:41659 while in Figure 7, ψ = 0:065. When the transmission rates are increasing, the populations in each compartment (exposed, infected, and recovered) become increasing accordingly. If there are low contact rates, then the exposed and infected individuals are low and as a result, the number of recovered individuals is also low.

Numerical
Solutions of the TB Submodel. The numerical solutions of the TB submodel with the basic reproduction number of R T 0 = 2:263 are given in Figure 8. The solutions are executed using the initial conditions of the state variables Sð0Þ = 5000, L T ð0Þ = 100, I T ð0Þ = 10, and Rð0Þ = 1. We observe that all solution curves approaches to their endemic equilibrium point. This supports the analytical results that corresponds to R T 0 > 1, where the endemic equilibrium point of the TB submodel is globally asymptotically stable. Figure 9 shows the impact of the contact rate (β 1 ) of TB endemic. As it increases, the number of individuals in the compartments, latent TB, active TB, and recovered ones is also increasing. If there is a low contact rate, then the transmission of TB disease will decrease as shown in Figures 9(a)-9(c).

Numerical Solutions of the Coinfection Model.
In this subsection, we present the numerical solutions of the coinfection model (2). In doing so, we use initial values of the state variables for the coinfection model given in Table 2, and all parameter values are illustrated in Table 3.
The time series graphs for the numerical solutions of the coinfection model (2) are plotted in Figure 10 with parameter values given in Table 3. The impact of treatment rate for the coinfected individuals with both diseases is illustrated in Figure 11. It was shown that an increase in the treatment rate for the disease will decrease the number of individuals in that class. Figure 12 displays the impact of transferring rates to coinfected with TB and COVID-19 from each actively infected individual of the two diseases. In Figure 13, we have shown that the impact of contact rates of the disease with infected compartments of system (2). Thus, we conclude that the states I T , I C , L TC , and I TC increase as the transmission coefficients increase.

Conclusion
In this study, a nonlinear deterministic mathematical modeling for the coinfection of TB and COVID-19 is proposed to examine the possible spread of these diseases. The biological meaningfulness of the model is proved by showing the existence, uniqueness, positivity, and boundedness of solutions in a given region. Then, the equilibria are computed for the coinfection model and each submodels independently. The stability analysis of these equilibria is also presented with the help of the calculated respective basic reproduction numbers. The analytical results in both submodels reveal that the disease-free equilibrium points are stable if the basic reproduction numbers are less than one otherwise unstable. Furthermore, there exists a stable endemic equilibrium point for the basic reproduction numbers greater than one. Besides, the stability analysis of coinfection disease-free equilibrium point is investigated via corresponding basic reproduction number. Different simulation cases were performed to supplement the analytical results, and it is observed to be in good agreement. Furthermore, the simulation result reveals that an increase in contact rate exacerbated the coinfection of both diseases, while a decrease in the effective contact rates and increase in treatments would minimize the spread of co-infection diseases.

Data Availability
All the necessary data was included in the main text.

Conflicts of Interest
The authors would declare that they have no known competing interests that could have appeared to influence this work.