# Mathematical Modeling and Analysis of TB and COVID-19 Coinfection

Authors: Kassahun Getnet Mekonen; Shiferaw Feyissa Balcha; Legesse Lemecha Obsu; Abdulkadir Hassen
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

Tuberculosis (TB) and coronavirus (COVID-19) are both infectious diseases that globally continue affecting millions of people every year. They have similar symptoms such as cough, fever, and difficulty breathing but differ 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 differential 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. Different simulation cases were performed to supplement the analytical results.


---

## Body

## 1. Introduction

The outbreak of coronavirus (COVID-19) 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 10th 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 byMycobacterium 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–16] and references cited therein. Mathematical models for TB are also extensively studied in [17–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 cholera-schistosomiasis 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 94th day in the absence of public health interventions and on the 138th 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 Section1. 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.


## 2. 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 (LT), active TB-infected individuals (IT), COVID-19 infected individuals with no symptoms but are infectious (EC), COVID-19 infected individuals with clinical symptoms and are infectious (IC), both latent TB and COVID-19 coinfected individuals (LTC), active TB and COVID-19 coinfected individuals (ITC), and recovered population from both TB and COVID-19 (R). With these in mind, the total population at time t, denoted by Nt, is given by
(1)Nt=St+LTt+ITt+ECt+ICt+LTCt+ITCt+Rt.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 (LT) leave to compartment IT 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 (IT) become recovered from the disease at a rate γ while the remaining portion transfers to the coinfection compartment ITC at a rate θ or die due to TB-induced death at a rate δT.The population in compartment (LTC) either progress to coinfected class ITC at a rate ρ or die with COVID-19-induced 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 LTC moves to class IT with a rate of qσ, to compartment IC class rate of pσ, and become recovered at a rate of 1−p+qσ. Furthermore, we assumed that individuals in compartment ITC leaves to compartments IT, IC, or R, respectively, at a rate of nτ, mτ, or 1−m+nτ while coinfection-induced death rate is denoted by δTC. In addition, the COVID-19 exposed individual (EC) has the chance to leave either to compartment LTC, IC, or R, respectively, at a rate of ϵλT, φ, or π. Similarly, the population in compartment (IC) 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.Figure 1
Flow diagram of the model.Based on the flow diagram, it results in systems of the following nonlinear differential equations:(2)dSdt=Λ−λT+λC+μS,dLTdt=λTS−μ+α+ηλC+ωLT,dITdt=αLT+qσLTC+nτITC−μ+θ+δT+γIT,dECdt=λCS−μ+ϵλT+φ+πEC,dICdt=φEC+pσLTC+mτITC−μ+δC+ν+ψIC,dLTCdt=ηλCLT+ϵλTEC−μ+δC+ρ+σLTC,dITCdt=ρLTC+θIT+νIC−μ+δTC+τITC,dRdt=ωLT+πEC+γIT+ψIC+1−p+qσLTC+1−m+nτITC−μR.where(3)λTt=β1NtITt+ITCt,(4)λCt=β2NtECt+ICt+LTCt+ITCt,and with nonnegative initial conditions S0>0, ET0≥0, IT0≥0, EC0≥0, IC0≥0, LTC0≥0, ITC0≥0, and R0≥0. The model parameters are described in Table 1.Table 1
Description of the model parameters.
ParameterBiological description of the parameterΛRate of recruitment to the susceptible individualsβ1Transmission rate of TB diseaseβ2COVID-19 transmission rateαFraction of TB exposed individuals to become infectedλTForce of infection for TB (it is the probability of acquiring infection from a TB infected individual)λCForce of infection for COVID-19 (it is the probability of acquiring infection from a COVID-19 infected individual)νTB infection rate from the COVID-19 infected individualsθCOVID-19 infection rate from the TB infected individualsπRecovery rate of the COVID-19 exposed individualsσRate at which individuals leave theLTH classpRecovery rate ofLTC individuals from TBqRecovery rate ofLTC individuals from COVID-19τRate at which individuals leave the coinfected classITHmTB recovery rate ofITC individualsnCOVID-19 recovery rate ofITC individualsηRate of TB exposed individuals become exposed with COVID-19γRecovery rate of TB infected individualsδTDeath rate due to TB bacteriaδCDeath rate due to the COVID-19 virusδTCDeath 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


### 2.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
(5)Ω=S,LT,IT,EC,IC,LTC,ITC,R∈ℝ+8:0≤Nt≤Λμ.Proof.
All functions on the right-hand side of equation (2) are C1 on ℝ+8 . Thus, by the Picard–Lindelo¨f theorem [33], equations (2) has a unique solution. Let the system of equations in (2) be written in the form x′=fx,t, where x=S,LT,IT,EC,IC,LTC,ITC,R, and f is the right hand sides of the equation (1). Following the results of proposition A.1 in [34], the functions fx,t has the property of
(6)fiS,LT,IT,EC,IC,LTC,ITC,R≥0,whenever x∈0,∞8andxi=0. As there exists a unique solution for the system of equations (2), it follows that xt∈0,∞8 for all t≥0 whenever x0≥0.
The change of total populationNt=S+LT+IT+EC+IC+LTC+ITC+R is governed by
(7)dNdt=Λ−δTIT−δcIc+LTC−δTCITC−μN≤Λ−μN.
The solution for this linear first order ode isNt≤N0e−μt+Λ/μ1−e−μt. Thus, for the initial data 0≤N0, we obtain 0≤Nt≤Λ/μ.
Hence, the solutions of the system of equations in (2) exists, unique, and bounded in a feasible region Ω, which concludes the proof.


## 2.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
(5)Ω=S,LT,IT,EC,IC,LTC,ITC,R∈ℝ+8:0≤Nt≤Λμ.Proof.
All functions on the right-hand side of equation (2) are C1 on ℝ+8 . Thus, by the Picard–Lindelo¨f theorem [33], equations (2) has a unique solution. Let the system of equations in (2) be written in the form x′=fx,t, where x=S,LT,IT,EC,IC,LTC,ITC,R, and f is the right hand sides of the equation (1). Following the results of proposition A.1 in [34], the functions fx,t has the property of
(6)fiS,LT,IT,EC,IC,LTC,ITC,R≥0,whenever x∈0,∞8andxi=0. As there exists a unique solution for the system of equations (2), it follows that xt∈0,∞8 for all t≥0 whenever x0≥0.
The change of total populationNt=S+LT+IT+EC+IC+LTC+ITC+R is governed by
(7)dNdt=Λ−δTIT−δcIc+LTC−δTCITC−μN≤Λ−μN.
The solution for this linear first order ode isNt≤N0e−μt+Λ/μ1−e−μt. Thus, for the initial data 0≤N0, we obtain 0≤Nt≤Λ/μ.
Hence, the solutions of the system of equations in (2) exists, unique, and bounded in a feasible region Ω, which concludes the proof.


## 3. 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.


### 3.1. COVID-19 Submodel

Without considering the infections of people with tuberculosis, the COVID-19 submodel is given by(8)dSdt=Λ−β2NEC+ICS−μS=f1,dECdt=β2NEC+ICS−μ+φ+πEC=f2,dICdt=φEC−μ+δC+ψIC=f3,dRdt=πEC+ψIC−μR=f4.


#### 3.1.1. Local Stability of Equilibrium Points

The equilibrium point of the COVID-19 submodel is a solution to the system of equations:(9)Λ−β2NEC+ICS−μS=0,β2NEC+ICS−μ+φ+πEC=0,φEC−μ+δC+ψIC=0,πEC+ψIC−μR=0.The disease-free equilibrium point is obtained by puttingEC=0 and IC=0 in equation (9). Hence, the disease-free equilibrium point, denoted by E0C, becomes E0C=Λ/μ,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 matrix as given in [35]. Let Fit,x be the input rate of newly infected individuals in the ith compartment and Vi+t,x and Vi−t,x, respectively, be the input rates by other means and rate of transfer of individuals out of a compartment i, where the xi’s are our state variables and Vit,x=Vi−t,x−Vi+t,x.Hence, from equation (8), it follows that
(10)Fi=β2NEC+ICS0,Vi+=0φECandVi−=μ+φ+πECμ+δC+ψIC.The Jacobian matrices ofFx and Vx, respectively, are
(11)F=DFE0C=β2β200andV=DVE0C=μ+φ+π0−φμ+δc+ψ.Hence, the next-generation matrixFV−1 becomes
(12)FV−1=β2μ+φ+π1+φμ+δc+ψβ2μ+δc+ψ00.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(13)R0C=β2μ+φ+πμ+δc+ψ+φμ+δc+ψ.
We use the Jacobean matrix to determine the local stability of the equilibrium points. For the submodel (8), the Jacobean matrix becomes
(14)JS,EC,IC,R=−β2EC+ICEC+IC+RN2−μ−β2S+RSN2−β2S+RSN20β2EC+ICEC+IC+RN2β2S+RSN2−μ+φ+πβ2S+RSN200φ−μ+δC+ψ00πψ−μ.
The Jacobean matrix of the submodel at the disease-free equilibrium pointE0c is given as
(15)JE0C=−μ−β2−β200β2−μ+φ+πβ200φ−μ+δC+ψ00πψ−μ.
Here, the eigenvalues forJE0C are λ1,2=−μ, and the remaining eigenvalues can be easily obtained from the submatrix:
(16)J2=β2−μ+φ+πβ2φ−μ+δC+ψ.
In order to determine the local stability of the disease-free equilibrium point, it is sufficient to show that the trace ofJ2 is less than zero, and its determinant is greater than zero.The trace ofJ2,trJ2=β2−2μ+φ+π+ψ+δC, which is less than zero ifβ2<2μ+φ+π+ψ+δc. Moreover, the determinant of J2,(17)detJ2=−β2−μ+φ+πμ+δC+ψ−φβ2=−β2μ+δC+ψ+φ+μ+φ+πμ+δC+ψ.This value is greater than zero ifβ2μ+δC+ψ+φ/μ+φ+πμ+δC+ψ<1, that is, detJ2>0 if R0C<1 and detJ2<0 if R0C>1. Hence, the disease-free equilibrium point of the COVID-19 submodel is stable if R0C<1 and unstable if R0C>1.The endemic equilibrium point of the COVID-19 submodel is the solution of the system of equation in (18).
(18)Λ−λC+μS=0,λCS−μ+φ+πEC=0,φEC−μ+δC+ψIC=0,πEC+ψIC−μR=0.Now, solving this system of equations in termsλC∗=β2/N∗EC∗+IC∗ one can obtain the following expressions for the endemic equilibrium point ΣC=S∗,EC∗,IC∗,R where
(19)S∗=ΛλC∗+μ,EC∗=λC∗Λc1λC∗+μ,IC∗=λC∗Λφc1c2λC∗+μandR∗=λC∗Λc2π+φψμc1c2λC∗+μ.Here,c1=μ+φ+π, c2=μ+δC+ψ, and N∗=Λ−δCIC∗/μ. Using (19) in the expression for λC∗, it gives λC∗=β2μc2+φ−c1c2μ/c1c2−δCμ. After simplification, we have in compact form
(20)λC∗=μc1c2β2c2+φ/c1c2−1μπc2+φμ+ψ=μc1c2R0C−1μπc2+φμ+ψ.This implies that the force of infectionλC∗ is positive at endemic equilibrium point only if R0C>1. Thus, we have just proved the following theorem.Theorem 3.
IfR0c>1, then the submodel system (8) has a unique endemic equilibrium point.


#### 3.1.2. 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:
(21)LS,EC,IC,R=12S−S∗+EC−EC∗+IC−IC∗+R−R∗2.The functionL 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
(22)dLdt=S−S∗+EC−EC∗+IC−IC∗+R−R∗dSdt+dECdt+dICdt+dRdt=N−Λ−δCIC∗μΛ−δCIC−μN≤N−ΛμΛ−μN≤−Λ−μN2μ≤0.This holds true forR0C>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]
(23)Sp=∂R0C∂ppR0C.In our context, the sensitivity analysis of each parameters for equation (8) becomes
(24)Sβ2=∂R0C∂β2β2R0C=1,Sμ=∂R0C∂μμR0C=−μδCδC+2μ+3ψ+φ+μμ+4ψ+2φ+ψ2ψ+2φ+π+φφ+πδC+μ+ψμ+π+φδC+μ+2ψ+φ,SδC=∂R0C∂δCδCR0C=−δCφδC+μ+ψδC+μ+ψ+φ,Sψ=∂R0C∂ψψR0C=−ψφδC+μ+ψδC+μ+ψ+φ,Sφ=∂R0C∂φφR0C=−φδC+ψ−πμ+π+φδC+μ+ψ+φ,Sπ=∂R0C∂ππR0C=−πμ+π+φ.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.


#### 3.1.4. Bifurcation Analysis for COVID-19 Submodel

Next, we study the solution behavior of equations (8) by taking β2 as the bifurcation parameter. Calculating the value of β2 from R0C=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 JE0c in equation (15) at β2=β2∗ has a left eigenvector (associated with the zero eigenvalue) which is calculated from vTJE0c=0. Here, v=v1,v2,v3,v4, where
(25)v1=0,v2=δC+μ+φ+ψμ+π+φ,v3=1,andv4=0.Similarly, the right eigenvector associated with the zero eigenvalue can be calculated fromJE0cw=0 with w=w1,w2,w3,w4 where
(26)w1=μδC+π+μ+π+φδC+μ+ψπδC+μ+ψ+ψφ,w2=μδC+μ+ψπδC+μ+ψ+ψφ,w3=μφπδC+μ+ψ+ψφ,andw4=1.Now, assume thatfk represents the right-hand side of the kth equation in the COVID-19 submodel (8) and let xk denote the corresponding state variable for k=1,⋯,4. Define
(27)a=∑k,i,j=1nvkwiwj∂2fk∂xi∂xj0,0andb=∑k,i=1nvkwi∂2fk∂xi∂ϕ0,0.The local dynamics of equation (8) near the bifurcation point β2=β∗ are then determined by the signs of two associated constants a and b with ϕ=β2−β∗. Note that, in fk0,0, the first zero corresponds to the DFE, E0C, for the submodel (4). In other words, fk0,ϕ=0, for k=1,⋯,4, if and only if the right-hand sides of equations (8) are equal to zero at E0C. Moreover, from ϕ=β2−β∗, we have ϕ=0 when β2=β∗, which is the second zero component in fk0,0.For the submodel (8), the associated nonzero partial derivatives at the disease-free equilibrium E0C are
(28)∂2f1∂EC2=−2β2∗μΛ,∂2f1∂EC∂IC=−2β2∗μΛ,∂2f1∂EC∂R=β2∗μΛ,∂2f1∂IC2=−2β2∗μΛ,∂2f1∂IC∂EC=−2β2∗μΛ,∂2f1∂IC∂R=β2∗μΛ,∂2f2∂EC2=2β2∗μΛ,∂2f2∂EC∂IC=2β2∗μΛ,∂2f2∂EC∂R=−β2∗μΛ,∂2f2∂IC2=2β2∗μΛ,∂2f2∂IC∂EC=2β2∗μΛ,∂2f2∂IC∂R=−β2∗μΛ.Then, using the above expressions fora, it follows that
(29)a=μ2μ+δC+ψδC+μ+φ+ψΛπδC+μ+ψ+ψφ2μδC+μ+ψ+φπδC+μ+ψ+ψφ−1.For the sign ofb, it can be shown that the associated nonvanishing partial derivatives are
(30)∂2f1∂EC∂β2=−1,∂2f1∂I3∂β2=−1∂2f2∂EC∂β2=1,∂2f2∂I3∂β2=1.Clearly, it also follows from the above expressions that(31)b=μδC+μ+ψ+φ2μ+π+φπδC+μ+ψ+ψφ.From the computations ofa 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.Theorem 4.
If the minimum valueκ=2μδC+μ+ψ+φ/πδC+μ+ψ+ψφ of the reinfection that induces bistability is less than unity, then there exists a forward bifurcation for the COVID-19 submodel (8).


### 3.2. TB Submodel

In the absence of COVID-19 disease, equation (2) automatically reduced to TB-submodel expressed as
(32)dSdt=Λ−β1NITS−μS,dLTdt=β1NITS−μ+α+ωLT,dITdt=αLT−μ+δT+γIT,dRdt=ωLT+γIT−μR.In a similar argument, the disease-free equilibrium point of equation (32) becomes E0T=Λ/μ,0,0,0, and the corresponding basic reproduction number is given by
(33)R0T=β1αμ+α+ωμ+δT+γ.Theorem 5.
The disease-free equilibrium point of the TB submodel (32) is locally asymptotically stable if R0T<1.Proof.
The Jacobean matrix at the disease-free equilibrium point for the TB submodel is given by(34)JE0T=−μ0−β100−μ+α+ωβ100α−μ+δT+γ00ωγ−μ.
Clearly, the two eigenvalues are equal to−μ, while the remaining can be obtained from the reduced matrix
(35)J3=−μ+α+ωβ1α−μ+δT+γ.
Here, the trace of the matrix istrJ3=−2μ+α+ω+δT+γ<0, and its determinant is detJ3=μ+α+ωμ+δT+γ−αβ1. Thus, the determinant is greater than zero if R0T<1. Therefore, the disease-free equilibrium point of the TB submodel is locally asymptotically stable if R0T<1. This completes the proof.
The endemic equilibrium point of the TB-submodel is obtained by solving the system of equations:(36)Λ−λTS−μS=0,λTS−μ+α+ωLT=0,αLT−μ+δT+γIT=0,ωLT+γIT−μR=0,where λT=β1IT/N. Solving this equation, we have
(37)S∗=ΛλT∗+μ,LT∗=λT∗Λk1λT∗+μ,IT∗=αλT∗Λk1k2λT∗+μ,andR∗=λT∗Λωk2+αγk1k2μλT∗+μ.
The simplified expression forλT in terms of S∗, LT∗, IT∗, and R∗ becomes
(38)λT∗=αβ1μ−k1k2μk1k2−αδT=αβ1μ1−1/R0Tk1k2−αδT.
Note that the endemic equilibrium point exists if and only ifR0T>1, since k1k2−αδT>0.Theorem 6.
The endemic equilibrium point of the TB submodel (32) is globally asymptotically stable if R0T>1.Proof.
To proof this theorem, we define the following function:(39)VS,LT,IT,R=12S−S∗+LT−LT∗+IT−IT∗+R−R∗2.
Next, we show that this functionV is a Lyapunov function. Differentiating the function V with respect to time t, we obtain
(40)dVdt=S−S∗+LT−LT∗+IT−IT∗+R−R∗dSdt+dLTdt+dITdt+dRdt=N−Λ−δTIT∗μΛ−δTIT−μN≤N−ΛμΛ−μN≤−Λ−μN2μ≤0.
ForR0T>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.


### 3.3. Analysis of the Coinfection Model

The equilibrium point of the coinfection model (2) is obtained by solving the system of equation (41):
(41)Λ−λT+λC+μS=0,λTS−μ+α+ηλC+ωLT=0,αLT+qσLTC+nτITC−μ+θ+δT+γIT=0,λCS−μ+ϵλT+φ+πEC=0,φEC+pσLTC+mτITC−μ+δC+ν+ψIC=0,ηλCLT+ϵλTEC−μ+δC+ρ+σLTC=0,ρLTC+θIT+νIC−μ+δTC+τITC=0,ωLT+πEC+γIT+ψIC+1−p+qσLTC+1−m+nτITC−μR=0,where λT and λC are as given in equations (3), and (4), respectively. The disease-free equilibrium point (E0) of the system becomes
(42)E0=Λμ,0,0,0,0,0,0,0.In the following, we calculate the basic reproduction number (R0) of the coinfection model using the method of next generation matrix [37]. For this, we use the notation X=LT,IT,EC,IC,LTC,ITC to denote the vector functions:
(43)FX=λTS0λCS000andVX=μ+α+ηλC+ωLTμ+θ+δT+γIT−αLT−qσLTC−nτITCμ+ϵλT+φ+πECμ+ν+δC+ψIC−φEC−pσLTC−mτITCμ+δC+ρ+σLTC−ηλCLT−ϵλTECμ+δTC+τITC−ρLTC−θIT−νIC,representing the appearance of new infections and the transfer of individuals into and out of the infected compartments, respectively. The Jacobian matrices of FX and VX at the disease free equilibrium point E0 are, respectively,
(44)F=DFE0=0β1000β100000000β2β2β2β2000000000000000000,and,V=DVE0=c100000−αc200−qσ−nτ00c300000−φc4−pσ−mτ0000c500−θ0−ν−ρc6,where c1=μ+α+ω, c2=μ+θ+δT+γ, c3=μ+φ+π, c4=μ+ν+δC+ψ, c5=μ+δC+ρ+σ, and c6=μ+τ+δTC. The next-generation matrix FV−1 is then given by
(45)FV−1=αA1c1BA1BφA3c3BA3BA5c5B−β1c4c2+nτB000000αA2c1BA2BA4c3B−β2c2c6+ν−nθτBA6c5B−β2c2c4+mτB000000000000000000,where
(46)A1=−β1c4c6+c4θ−mντ,A2=−β2θc4+mτ,A3=−β1νc2+nτ,A4=−β2c2c4c6+c6φ+νφ−mντ−nθτc4+φ,A5=−β1c4ρc2+ντ+c4qσc6+θ+c2pνσ+νρτnp−mq,A6=−β2c2c4c6+ρ+c2mτρ−ν+c2pσc6+ν+qσθc4+mτ−nθτc4+pσ,B=c2mντ+c4nθτ−c2c6.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 matrixFV−1 are equal to zero. The remaining eigenvalues are obtained from the reduced matrix αA1/c1BφA3/c3BαA2/c1BA4/c3B. Thus, computing the eigenvalues of FV−1 one can easily obtain that
(47)λ1=c1A4+αc3A1−D2c1c3B,λ2=c1A4+αc3A1+D2c1c3B,where D=α2c32A12−2αc1c3A1A2+4φαc1c3A22+c12A42. Hence, the basic reproduction number R0 of the coinfection model (1) is given by
(48)R0=max∣λ1∣,∣λ2∣.


#### 3.3.1. Local Stability Analysis of the Disease-Free Equilibrium PointE0

The local stability analysis of the disease-free equilibrium point for the coinfection model is analyzed using the method of linearization. The Jacobean matrix of the model equation (2) at disease-free equilibrium point E0 represented as
(49)JE0=−μ0−β1−β2−β2−β2−β1−β200−c1β1000β100α−c200qσnτ0000β2−c3β2β2β20000φ−c4pσmτ000000−c50000θ0νρ−c600ωγπψ1−p+qσ1−m+nτ−μ.The local asymptotic stability ofE0 is analyzed based on the sign of the eigenvalues. Expanding the characteristic equation ∣λI8−JE0∣=0 with the first and eighth columns, and the sixth row, we obtain three eigenvalues λ1,2=−μ and λ3=−c5. Further, we obtain the remaining five eigenvalues from the characteristic equation of the reduced matrix:
(50)J5−λI5=−c1−λβ100β1α−c2−λ00nτ00β2−c3−λβ2β200φ−c4−λmτ0θ0ν−c6−λ.Here, the characteristic polynomial is computed using elementary row operations method [38]. Thus, the row reduced form of the matrix becomes
(51)α−c2−λ00nτ0θ0ν−c6−λ00φ−c4−λmτ000β2+β2−c3−λc4+λφβ2+β2−c3−λc4+λφ0000β11+c6+λθ+c2+λαnτ−c2+λc6+λθ.Hence, the characteristic polynomial of our matrix is the product of the diagonal entries times−1k, where k is the number of row swaps. We perform three row swaps, so that the characteristic polynomial is
(52)X5λ=−αθφβ2+β2−c3−λc4+λφβ11+c6+λθ+c2+λαnτ−c2+λc6+λθ.Simplifying it gives(53)X5λ=−β2φ+c4β2−c3+β2−c3−c4λ−λ2,αβ1θ+c6+c1nτθ−c2c6+αβ1+nτθ−c2c6−c1c2+c6λ+c2+c6−c1λ2−λ3.That is, the eigenvalues are the solutions of the characteristic equation:(54)X5λ=λ5+D1λ4+D2λ3+D3λ2+D4λ+D5=0,where
(55)D1=β2+c2+c6−c3+c4+c1,D2=β2φ+c4β2−c3−β2−c3−c4c2+c6−c1+αβ1+nτθ−c2c6−c1c2+c6,D3=αβ1θ+c6+c1nτθ−c2c6−β2φ+c4β2−c3c2+c6−c1−β2−c3−c4αβ1+nτθ−c2c6−c1c2+c6,D4=−β2φ+c4β2−c3αβ1+nτθ−c2c6−c1c2+c6−β2−c3−c4αβ1θ+c6+c1nτθ−c2c6,D5=β2φ+c4β2−c3αβ1θ+c6+c1nτθ−c2c6.Thus, using the Routh-Hurwitz stability criterion, the roots of the equation (54) have negative real parts if the following condition holds.
(56)Di>0,foralli=1,⋯,5,D1D2D3>D32+D12D4,andD1D4−D5D1D2D3−D32−D12D4>D5D1D2−D32+D1D52.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) holds.


#### 3.3.2. Global Stability of Disease-Free EquilibriumE0

Let us rewrite the coinfection model (2) as
(57)dXdt=FX,Z,dZdt=GX,Z,GX,0=0.where X=S,R and Z=LT,IT,EC,IC,LTC,ITC. 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
(58)U0=X0,0,where,X0=Λμ,0.The conditionsH1 and H2 below must be met to guarantee global asymptotically stability:(H1) For dX/dt=FX,0, the equilibrium point U0 is globally stable(H2) GX,Z=AZ−G^X,Z,G^X,Z≥0 for X,Z∈Ω, where A=DZGU0,0 is a Metzler matrix (the nondiagonal entities of A are nonnegative), and Ω is the biological feasible region of the modelFrom the coinfection model (2), we have
(59)dXdt=FX,Z=Λ−λT+λC+μSωLT+πEC+γIT+ψIC+1−p+qσLTC+1−m+nτITC−μR.Hence,(60)FX,0=Λ−μS0,(61)dZdt=GX,Z=λTS−μ+α+ηλC+ωLTαLT+qσLTC+nτITC−μ+θ+δT+γITλCS−μ+ϵλT+φ+πECφEC+pσLTC+mτITC−μ+δC+ν+ψICηλCLT+ϵλTEC−μ+δC+ρ+σLTCρLTC+θIT+νIC−μ+δTC+τITC.Therefore,(62)A=DZGU0,0=−c1β1000β1α−c200qσnτ00β2−c3β2β2β200φ−c4pσmτ0000−c500θ0νρ−c6whichisaMetzlerMatrix.Here,G^X,Z=AZ−GX,Z, and hence
(63)G^X,Z=G^1X,ZG^2X,ZG^3X,ZG^4X,ZG^5X,ZG^6X,Z=ηλCLT+β1IT+ITC−λTS0β2EC+IC+LTC+ITC+ϵλTEC−λCS0−ηλCLT−ϵλTEC0.Thus,G^5X,Z=−ηλCLT−ϵλTEC<0. This shows that condition (H2) is not satisfied. As a result, U0 and then the disease-free equilibrium point E0 cannot be globally asymptotically stable.


### 3.4. 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 numbersR0C in terms of R0T (and vice versa) [40]. Expressing the parameter μ in the equation (33) in terms of R0T, we have
(64)R0T=β1αμ+k1μ+k2,wherek1=α+ω,k2=δT+γ.Solving forμ and simplifying, we obtain
(65)μ=R0Tk1−k22R0T+4αβ1−k1+k2R0T2R0T.Substituting this value ofμ in R0C (13), we obtain R0C in terms of R0T as
(66)R0C=β2R0Th1+h2/2h2/2+φ+ψR0T−k1+k2R0T/2h2/2+δC+ψR0T−k1+k2R0T/2,where h1=R0TδC+φ+ψ−k1+k2/2, and h2=R0Tk1−k22R0T+4αβ1.The partial derivative ofR0C with respect to R0T is then given by
(67)∂R0C∂R0T=β2R0Th1h2+h3+h4h2k1+k2−k1−k22R0T−2αβ12h2h2/2+φ+ψR0T−k1+k2R0T/22h2/2+δC+ψR0T−k1+k2R0T/22,with h3=φ2+π+δC+ψφ+δC+ψ2R0T2, and
(68)h4=R0Tk12+k22R0T+2αβ12−k1+k2δC+φ+ψR0T.Here, we observe that∂R0C/∂R0T is positive. This implies that the expansion of TB infection exacerbates the pandemic of COVID-19.Remark 8.
If∂R0C/∂R0T=0, the expansion of TB in the community has no significant impact on the spread of COVID-19. In the contrary, when ∂R0C/∂R0T<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 expressingR0T in terms of R0c and identifying the sign of the partial derivative of R0T with respect to R0c.


## 3.1. COVID-19 Submodel

Without considering the infections of people with tuberculosis, the COVID-19 submodel is given by(8)dSdt=Λ−β2NEC+ICS−μS=f1,dECdt=β2NEC+ICS−μ+φ+πEC=f2,dICdt=φEC−μ+δC+ψIC=f3,dRdt=πEC+ψIC−μR=f4.


### 3.1.1. Local Stability of Equilibrium Points

The equilibrium point of the COVID-19 submodel is a solution to the system of equations:(9)Λ−β2NEC+ICS−μS=0,β2NEC+ICS−μ+φ+πEC=0,φEC−μ+δC+ψIC=0,πEC+ψIC−μR=0.The disease-free equilibrium point is obtained by puttingEC=0 and IC=0 in equation (9). Hence, the disease-free equilibrium point, denoted by E0C, becomes E0C=Λ/μ,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 matrix as given in [35]. Let Fit,x be the input rate of newly infected individuals in the ith compartment and Vi+t,x and Vi−t,x, respectively, be the input rates by other means and rate of transfer of individuals out of a compartment i, where the xi’s are our state variables and Vit,x=Vi−t,x−Vi+t,x.Hence, from equation (8), it follows that
(10)Fi=β2NEC+ICS0,Vi+=0φECandVi−=μ+φ+πECμ+δC+ψIC.The Jacobian matrices ofFx and Vx, respectively, are
(11)F=DFE0C=β2β200andV=DVE0C=μ+φ+π0−φμ+δc+ψ.Hence, the next-generation matrixFV−1 becomes
(12)FV−1=β2μ+φ+π1+φμ+δc+ψβ2μ+δc+ψ00.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(13)R0C=β2μ+φ+πμ+δc+ψ+φμ+δc+ψ.
We use the Jacobean matrix to determine the local stability of the equilibrium points. For the submodel (8), the Jacobean matrix becomes
(14)JS,EC,IC,R=−β2EC+ICEC+IC+RN2−μ−β2S+RSN2−β2S+RSN20β2EC+ICEC+IC+RN2β2S+RSN2−μ+φ+πβ2S+RSN200φ−μ+δC+ψ00πψ−μ.
The Jacobean matrix of the submodel at the disease-free equilibrium pointE0c is given as
(15)JE0C=−μ−β2−β200β2−μ+φ+πβ200φ−μ+δC+ψ00πψ−μ.
Here, the eigenvalues forJE0C are λ1,2=−μ, and the remaining eigenvalues can be easily obtained from the submatrix:
(16)J2=β2−μ+φ+πβ2φ−μ+δC+ψ.
In order to determine the local stability of the disease-free equilibrium point, it is sufficient to show that the trace ofJ2 is less than zero, and its determinant is greater than zero.The trace ofJ2,trJ2=β2−2μ+φ+π+ψ+δC, which is less than zero ifβ2<2μ+φ+π+ψ+δc. Moreover, the determinant of J2,(17)detJ2=−β2−μ+φ+πμ+δC+ψ−φβ2=−β2μ+δC+ψ+φ+μ+φ+πμ+δC+ψ.This value is greater than zero ifβ2μ+δC+ψ+φ/μ+φ+πμ+δC+ψ<1, that is, detJ2>0 if R0C<1 and detJ2<0 if R0C>1. Hence, the disease-free equilibrium point of the COVID-19 submodel is stable if R0C<1 and unstable if R0C>1.The endemic equilibrium point of the COVID-19 submodel is the solution of the system of equation in (18).
(18)Λ−λC+μS=0,λCS−μ+φ+πEC=0,φEC−μ+δC+ψIC=0,πEC+ψIC−μR=0.Now, solving this system of equations in termsλC∗=β2/N∗EC∗+IC∗ one can obtain the following expressions for the endemic equilibrium point ΣC=S∗,EC∗,IC∗,R where
(19)S∗=ΛλC∗+μ,EC∗=λC∗Λc1λC∗+μ,IC∗=λC∗Λφc1c2λC∗+μandR∗=λC∗Λc2π+φψμc1c2λC∗+μ.Here,c1=μ+φ+π, c2=μ+δC+ψ, and N∗=Λ−δCIC∗/μ. Using (19) in the expression for λC∗, it gives λC∗=β2μc2+φ−c1c2μ/c1c2−δCμ. After simplification, we have in compact form
(20)λC∗=μc1c2β2c2+φ/c1c2−1μπc2+φμ+ψ=μc1c2R0C−1μπc2+φμ+ψ.This implies that the force of infectionλC∗ is positive at endemic equilibrium point only if R0C>1. Thus, we have just proved the following theorem.Theorem 3.
IfR0c>1, then the submodel system (8) has a unique endemic equilibrium point.


### 3.1.2. 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:
(21)LS,EC,IC,R=12S−S∗+EC−EC∗+IC−IC∗+R−R∗2.The functionL 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
(22)dLdt=S−S∗+EC−EC∗+IC−IC∗+R−R∗dSdt+dECdt+dICdt+dRdt=N−Λ−δCIC∗μΛ−δCIC−μN≤N−ΛμΛ−μN≤−Λ−μN2μ≤0.This holds true forR0C>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]
(23)Sp=∂R0C∂ppR0C.In our context, the sensitivity analysis of each parameters for equation (8) becomes
(24)Sβ2=∂R0C∂β2β2R0C=1,Sμ=∂R0C∂μμR0C=−μδCδC+2μ+3ψ+φ+μμ+4ψ+2φ+ψ2ψ+2φ+π+φφ+πδC+μ+ψμ+π+φδC+μ+2ψ+φ,SδC=∂R0C∂δCδCR0C=−δCφδC+μ+ψδC+μ+ψ+φ,Sψ=∂R0C∂ψψR0C=−ψφδC+μ+ψδC+μ+ψ+φ,Sφ=∂R0C∂φφR0C=−φδC+ψ−πμ+π+φδC+μ+ψ+φ,Sπ=∂R0C∂ππR0C=−πμ+π+φ.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.


### 3.1.4. Bifurcation Analysis for COVID-19 Submodel

Next, we study the solution behavior of equations (8) by taking β2 as the bifurcation parameter. Calculating the value of β2 from R0C=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 JE0c in equation (15) at β2=β2∗ has a left eigenvector (associated with the zero eigenvalue) which is calculated from vTJE0c=0. Here, v=v1,v2,v3,v4, where
(25)v1=0,v2=δC+μ+φ+ψμ+π+φ,v3=1,andv4=0.Similarly, the right eigenvector associated with the zero eigenvalue can be calculated fromJE0cw=0 with w=w1,w2,w3,w4 where
(26)w1=μδC+π+μ+π+φδC+μ+ψπδC+μ+ψ+ψφ,w2=μδC+μ+ψπδC+μ+ψ+ψφ,w3=μφπδC+μ+ψ+ψφ,andw4=1.Now, assume thatfk represents the right-hand side of the kth equation in the COVID-19 submodel (8) and let xk denote the corresponding state variable for k=1,⋯,4. Define
(27)a=∑k,i,j=1nvkwiwj∂2fk∂xi∂xj0,0andb=∑k,i=1nvkwi∂2fk∂xi∂ϕ0,0.The local dynamics of equation (8) near the bifurcation point β2=β∗ are then determined by the signs of two associated constants a and b with ϕ=β2−β∗. Note that, in fk0,0, the first zero corresponds to the DFE, E0C, for the submodel (4). In other words, fk0,ϕ=0, for k=1,⋯,4, if and only if the right-hand sides of equations (8) are equal to zero at E0C. Moreover, from ϕ=β2−β∗, we have ϕ=0 when β2=β∗, which is the second zero component in fk0,0.For the submodel (8), the associated nonzero partial derivatives at the disease-free equilibrium E0C are
(28)∂2f1∂EC2=−2β2∗μΛ,∂2f1∂EC∂IC=−2β2∗μΛ,∂2f1∂EC∂R=β2∗μΛ,∂2f1∂IC2=−2β2∗μΛ,∂2f1∂IC∂EC=−2β2∗μΛ,∂2f1∂IC∂R=β2∗μΛ,∂2f2∂EC2=2β2∗μΛ,∂2f2∂EC∂IC=2β2∗μΛ,∂2f2∂EC∂R=−β2∗μΛ,∂2f2∂IC2=2β2∗μΛ,∂2f2∂IC∂EC=2β2∗μΛ,∂2f2∂IC∂R=−β2∗μΛ.Then, using the above expressions fora, it follows that
(29)a=μ2μ+δC+ψδC+μ+φ+ψΛπδC+μ+ψ+ψφ2μδC+μ+ψ+φπδC+μ+ψ+ψφ−1.For the sign ofb, it can be shown that the associated nonvanishing partial derivatives are
(30)∂2f1∂EC∂β2=−1,∂2f1∂I3∂β2=−1∂2f2∂EC∂β2=1,∂2f2∂I3∂β2=1.Clearly, it also follows from the above expressions that(31)b=μδC+μ+ψ+φ2μ+π+φπδC+μ+ψ+ψφ.From the computations ofa 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.Theorem 4.
If the minimum valueκ=2μδC+μ+ψ+φ/πδC+μ+ψ+ψφ of the reinfection that induces bistability is less than unity, then there exists a forward bifurcation for the COVID-19 submodel (8).


## 3.1.1. Local Stability of Equilibrium Points

The equilibrium point of the COVID-19 submodel is a solution to the system of equations:(9)Λ−β2NEC+ICS−μS=0,β2NEC+ICS−μ+φ+πEC=0,φEC−μ+δC+ψIC=0,πEC+ψIC−μR=0.The disease-free equilibrium point is obtained by puttingEC=0 and IC=0 in equation (9). Hence, the disease-free equilibrium point, denoted by E0C, becomes E0C=Λ/μ,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 matrix as given in [35]. Let Fit,x be the input rate of newly infected individuals in the ith compartment and Vi+t,x and Vi−t,x, respectively, be the input rates by other means and rate of transfer of individuals out of a compartment i, where the xi’s are our state variables and Vit,x=Vi−t,x−Vi+t,x.Hence, from equation (8), it follows that
(10)Fi=β2NEC+ICS0,Vi+=0φECandVi−=μ+φ+πECμ+δC+ψIC.The Jacobian matrices ofFx and Vx, respectively, are
(11)F=DFE0C=β2β200andV=DVE0C=μ+φ+π0−φμ+δc+ψ.Hence, the next-generation matrixFV−1 becomes
(12)FV−1=β2μ+φ+π1+φμ+δc+ψβ2μ+δc+ψ00.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(13)R0C=β2μ+φ+πμ+δc+ψ+φμ+δc+ψ.
We use the Jacobean matrix to determine the local stability of the equilibrium points. For the submodel (8), the Jacobean matrix becomes
(14)JS,EC,IC,R=−β2EC+ICEC+IC+RN2−μ−β2S+RSN2−β2S+RSN20β2EC+ICEC+IC+RN2β2S+RSN2−μ+φ+πβ2S+RSN200φ−μ+δC+ψ00πψ−μ.
The Jacobean matrix of the submodel at the disease-free equilibrium pointE0c is given as
(15)JE0C=−μ−β2−β200β2−μ+φ+πβ200φ−μ+δC+ψ00πψ−μ.
Here, the eigenvalues forJE0C are λ1,2=−μ, and the remaining eigenvalues can be easily obtained from the submatrix:
(16)J2=β2−μ+φ+πβ2φ−μ+δC+ψ.
In order to determine the local stability of the disease-free equilibrium point, it is sufficient to show that the trace ofJ2 is less than zero, and its determinant is greater than zero.The trace ofJ2,trJ2=β2−2μ+φ+π+ψ+δC, which is less than zero ifβ2<2μ+φ+π+ψ+δc. Moreover, the determinant of J2,(17)detJ2=−β2−μ+φ+πμ+δC+ψ−φβ2=−β2μ+δC+ψ+φ+μ+φ+πμ+δC+ψ.This value is greater than zero ifβ2μ+δC+ψ+φ/μ+φ+πμ+δC+ψ<1, that is, detJ2>0 if R0C<1 and detJ2<0 if R0C>1. Hence, the disease-free equilibrium point of the COVID-19 submodel is stable if R0C<1 and unstable if R0C>1.The endemic equilibrium point of the COVID-19 submodel is the solution of the system of equation in (18).
(18)Λ−λC+μS=0,λCS−μ+φ+πEC=0,φEC−μ+δC+ψIC=0,πEC+ψIC−μR=0.Now, solving this system of equations in termsλC∗=β2/N∗EC∗+IC∗ one can obtain the following expressions for the endemic equilibrium point ΣC=S∗,EC∗,IC∗,R where
(19)S∗=ΛλC∗+μ,EC∗=λC∗Λc1λC∗+μ,IC∗=λC∗Λφc1c2λC∗+μandR∗=λC∗Λc2π+φψμc1c2λC∗+μ.Here,c1=μ+φ+π, c2=μ+δC+ψ, and N∗=Λ−δCIC∗/μ. Using (19) in the expression for λC∗, it gives λC∗=β2μc2+φ−c1c2μ/c1c2−δCμ. After simplification, we have in compact form
(20)λC∗=μc1c2β2c2+φ/c1c2−1μπc2+φμ+ψ=μc1c2R0C−1μπc2+φμ+ψ.This implies that the force of infectionλC∗ is positive at endemic equilibrium point only if R0C>1. Thus, we have just proved the following theorem.Theorem 3.
IfR0c>1, then the submodel system (8) has a unique endemic equilibrium point.


## 3.1.2. 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:
(21)LS,EC,IC,R=12S−S∗+EC−EC∗+IC−IC∗+R−R∗2.The functionL 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
(22)dLdt=S−S∗+EC−EC∗+IC−IC∗+R−R∗dSdt+dECdt+dICdt+dRdt=N−Λ−δCIC∗μΛ−δCIC−μN≤N−ΛμΛ−μN≤−Λ−μN2μ≤0.This holds true forR0C>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]
(23)Sp=∂R0C∂ppR0C.In our context, the sensitivity analysis of each parameters for equation (8) becomes
(24)Sβ2=∂R0C∂β2β2R0C=1,Sμ=∂R0C∂μμR0C=−μδCδC+2μ+3ψ+φ+μμ+4ψ+2φ+ψ2ψ+2φ+π+φφ+πδC+μ+ψμ+π+φδC+μ+2ψ+φ,SδC=∂R0C∂δCδCR0C=−δCφδC+μ+ψδC+μ+ψ+φ,Sψ=∂R0C∂ψψR0C=−ψφδC+μ+ψδC+μ+ψ+φ,Sφ=∂R0C∂φφR0C=−φδC+ψ−πμ+π+φδC+μ+ψ+φ,Sπ=∂R0C∂ππR0C=−πμ+π+φ.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.


## 3.1.4. Bifurcation Analysis for COVID-19 Submodel

Next, we study the solution behavior of equations (8) by taking β2 as the bifurcation parameter. Calculating the value of β2 from R0C=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 JE0c in equation (15) at β2=β2∗ has a left eigenvector (associated with the zero eigenvalue) which is calculated from vTJE0c=0. Here, v=v1,v2,v3,v4, where
(25)v1=0,v2=δC+μ+φ+ψμ+π+φ,v3=1,andv4=0.Similarly, the right eigenvector associated with the zero eigenvalue can be calculated fromJE0cw=0 with w=w1,w2,w3,w4 where
(26)w1=μδC+π+μ+π+φδC+μ+ψπδC+μ+ψ+ψφ,w2=μδC+μ+ψπδC+μ+ψ+ψφ,w3=μφπδC+μ+ψ+ψφ,andw4=1.Now, assume thatfk represents the right-hand side of the kth equation in the COVID-19 submodel (8) and let xk denote the corresponding state variable for k=1,⋯,4. Define
(27)a=∑k,i,j=1nvkwiwj∂2fk∂xi∂xj0,0andb=∑k,i=1nvkwi∂2fk∂xi∂ϕ0,0.The local dynamics of equation (8) near the bifurcation point β2=β∗ are then determined by the signs of two associated constants a and b with ϕ=β2−β∗. Note that, in fk0,0, the first zero corresponds to the DFE, E0C, for the submodel (4). In other words, fk0,ϕ=0, for k=1,⋯,4, if and only if the right-hand sides of equations (8) are equal to zero at E0C. Moreover, from ϕ=β2−β∗, we have ϕ=0 when β2=β∗, which is the second zero component in fk0,0.For the submodel (8), the associated nonzero partial derivatives at the disease-free equilibrium E0C are
(28)∂2f1∂EC2=−2β2∗μΛ,∂2f1∂EC∂IC=−2β2∗μΛ,∂2f1∂EC∂R=β2∗μΛ,∂2f1∂IC2=−2β2∗μΛ,∂2f1∂IC∂EC=−2β2∗μΛ,∂2f1∂IC∂R=β2∗μΛ,∂2f2∂EC2=2β2∗μΛ,∂2f2∂EC∂IC=2β2∗μΛ,∂2f2∂EC∂R=−β2∗μΛ,∂2f2∂IC2=2β2∗μΛ,∂2f2∂IC∂EC=2β2∗μΛ,∂2f2∂IC∂R=−β2∗μΛ.Then, using the above expressions fora, it follows that
(29)a=μ2μ+δC+ψδC+μ+φ+ψΛπδC+μ+ψ+ψφ2μδC+μ+ψ+φπδC+μ+ψ+ψφ−1.For the sign ofb, it can be shown that the associated nonvanishing partial derivatives are
(30)∂2f1∂EC∂β2=−1,∂2f1∂I3∂β2=−1∂2f2∂EC∂β2=1,∂2f2∂I3∂β2=1.Clearly, it also follows from the above expressions that(31)b=μδC+μ+ψ+φ2μ+π+φπδC+μ+ψ+ψφ.From the computations ofa 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.Theorem 4.
If the minimum valueκ=2μδC+μ+ψ+φ/πδC+μ+ψ+ψφ of the reinfection that induces bistability is less than unity, then there exists a forward bifurcation for the COVID-19 submodel (8).


## 3.2. TB Submodel

In the absence of COVID-19 disease, equation (2) automatically reduced to TB-submodel expressed as
(32)dSdt=Λ−β1NITS−μS,dLTdt=β1NITS−μ+α+ωLT,dITdt=αLT−μ+δT+γIT,dRdt=ωLT+γIT−μR.In a similar argument, the disease-free equilibrium point of equation (32) becomes E0T=Λ/μ,0,0,0, and the corresponding basic reproduction number is given by
(33)R0T=β1αμ+α+ωμ+δT+γ.Theorem 5.
The disease-free equilibrium point of the TB submodel (32) is locally asymptotically stable if R0T<1.Proof.
The Jacobean matrix at the disease-free equilibrium point for the TB submodel is given by(34)JE0T=−μ0−β100−μ+α+ωβ100α−μ+δT+γ00ωγ−μ.
Clearly, the two eigenvalues are equal to−μ, while the remaining can be obtained from the reduced matrix
(35)J3=−μ+α+ωβ1α−μ+δT+γ.
Here, the trace of the matrix istrJ3=−2μ+α+ω+δT+γ<0, and its determinant is detJ3=μ+α+ωμ+δT+γ−αβ1. Thus, the determinant is greater than zero if R0T<1. Therefore, the disease-free equilibrium point of the TB submodel is locally asymptotically stable if R0T<1. This completes the proof.
The endemic equilibrium point of the TB-submodel is obtained by solving the system of equations:(36)Λ−λTS−μS=0,λTS−μ+α+ωLT=0,αLT−μ+δT+γIT=0,ωLT+γIT−μR=0,where λT=β1IT/N. Solving this equation, we have
(37)S∗=ΛλT∗+μ,LT∗=λT∗Λk1λT∗+μ,IT∗=αλT∗Λk1k2λT∗+μ,andR∗=λT∗Λωk2+αγk1k2μλT∗+μ.
The simplified expression forλT in terms of S∗, LT∗, IT∗, and R∗ becomes
(38)λT∗=αβ1μ−k1k2μk1k2−αδT=αβ1μ1−1/R0Tk1k2−αδT.
Note that the endemic equilibrium point exists if and only ifR0T>1, since k1k2−αδT>0.Theorem 6.
The endemic equilibrium point of the TB submodel (32) is globally asymptotically stable if R0T>1.Proof.
To proof this theorem, we define the following function:(39)VS,LT,IT,R=12S−S∗+LT−LT∗+IT−IT∗+R−R∗2.
Next, we show that this functionV is a Lyapunov function. Differentiating the function V with respect to time t, we obtain
(40)dVdt=S−S∗+LT−LT∗+IT−IT∗+R−R∗dSdt+dLTdt+dITdt+dRdt=N−Λ−δTIT∗μΛ−δTIT−μN≤N−ΛμΛ−μN≤−Λ−μN2μ≤0.
ForR0T>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.


## 3.3. Analysis of the Coinfection Model

The equilibrium point of the coinfection model (2) is obtained by solving the system of equation (41):
(41)Λ−λT+λC+μS=0,λTS−μ+α+ηλC+ωLT=0,αLT+qσLTC+nτITC−μ+θ+δT+γIT=0,λCS−μ+ϵλT+φ+πEC=0,φEC+pσLTC+mτITC−μ+δC+ν+ψIC=0,ηλCLT+ϵλTEC−μ+δC+ρ+σLTC=0,ρLTC+θIT+νIC−μ+δTC+τITC=0,ωLT+πEC+γIT+ψIC+1−p+qσLTC+1−m+nτITC−μR=0,where λT and λC are as given in equations (3), and (4), respectively. The disease-free equilibrium point (E0) of the system becomes
(42)E0=Λμ,0,0,0,0,0,0,0.In the following, we calculate the basic reproduction number (R0) of the coinfection model using the method of next generation matrix [37]. For this, we use the notation X=LT,IT,EC,IC,LTC,ITC to denote the vector functions:
(43)FX=λTS0λCS000andVX=μ+α+ηλC+ωLTμ+θ+δT+γIT−αLT−qσLTC−nτITCμ+ϵλT+φ+πECμ+ν+δC+ψIC−φEC−pσLTC−mτITCμ+δC+ρ+σLTC−ηλCLT−ϵλTECμ+δTC+τITC−ρLTC−θIT−νIC,representing the appearance of new infections and the transfer of individuals into and out of the infected compartments, respectively. The Jacobian matrices of FX and VX at the disease free equilibrium point E0 are, respectively,
(44)F=DFE0=0β1000β100000000β2β2β2β2000000000000000000,and,V=DVE0=c100000−αc200−qσ−nτ00c300000−φc4−pσ−mτ0000c500−θ0−ν−ρc6,where c1=μ+α+ω, c2=μ+θ+δT+γ, c3=μ+φ+π, c4=μ+ν+δC+ψ, c5=μ+δC+ρ+σ, and c6=μ+τ+δTC. The next-generation matrix FV−1 is then given by
(45)FV−1=αA1c1BA1BφA3c3BA3BA5c5B−β1c4c2+nτB000000αA2c1BA2BA4c3B−β2c2c6+ν−nθτBA6c5B−β2c2c4+mτB000000000000000000,where
(46)A1=−β1c4c6+c4θ−mντ,A2=−β2θc4+mτ,A3=−β1νc2+nτ,A4=−β2c2c4c6+c6φ+νφ−mντ−nθτc4+φ,A5=−β1c4ρc2+ντ+c4qσc6+θ+c2pνσ+νρτnp−mq,A6=−β2c2c4c6+ρ+c2mτρ−ν+c2pσc6+ν+qσθc4+mτ−nθτc4+pσ,B=c2mντ+c4nθτ−c2c6.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 matrixFV−1 are equal to zero. The remaining eigenvalues are obtained from the reduced matrix αA1/c1BφA3/c3BαA2/c1BA4/c3B. Thus, computing the eigenvalues of FV−1 one can easily obtain that
(47)λ1=c1A4+αc3A1−D2c1c3B,λ2=c1A4+αc3A1+D2c1c3B,where D=α2c32A12−2αc1c3A1A2+4φαc1c3A22+c12A42. Hence, the basic reproduction number R0 of the coinfection model (1) is given by
(48)R0=max∣λ1∣,∣λ2∣.


### 3.3.1. Local Stability Analysis of the Disease-Free Equilibrium PointE0

The local stability analysis of the disease-free equilibrium point for the coinfection model is analyzed using the method of linearization. The Jacobean matrix of the model equation (2) at disease-free equilibrium point E0 represented as
(49)JE0=−μ0−β1−β2−β2−β2−β1−β200−c1β1000β100α−c200qσnτ0000β2−c3β2β2β20000φ−c4pσmτ000000−c50000θ0νρ−c600ωγπψ1−p+qσ1−m+nτ−μ.The local asymptotic stability ofE0 is analyzed based on the sign of the eigenvalues. Expanding the characteristic equation ∣λI8−JE0∣=0 with the first and eighth columns, and the sixth row, we obtain three eigenvalues λ1,2=−μ and λ3=−c5. Further, we obtain the remaining five eigenvalues from the characteristic equation of the reduced matrix:
(50)J5−λI5=−c1−λβ100β1α−c2−λ00nτ00β2−c3−λβ2β200φ−c4−λmτ0θ0ν−c6−λ.Here, the characteristic polynomial is computed using elementary row operations method [38]. Thus, the row reduced form of the matrix becomes
(51)α−c2−λ00nτ0θ0ν−c6−λ00φ−c4−λmτ000β2+β2−c3−λc4+λφβ2+β2−c3−λc4+λφ0000β11+c6+λθ+c2+λαnτ−c2+λc6+λθ.Hence, the characteristic polynomial of our matrix is the product of the diagonal entries times−1k, where k is the number of row swaps. We perform three row swaps, so that the characteristic polynomial is
(52)X5λ=−αθφβ2+β2−c3−λc4+λφβ11+c6+λθ+c2+λαnτ−c2+λc6+λθ.Simplifying it gives(53)X5λ=−β2φ+c4β2−c3+β2−c3−c4λ−λ2,αβ1θ+c6+c1nτθ−c2c6+αβ1+nτθ−c2c6−c1c2+c6λ+c2+c6−c1λ2−λ3.That is, the eigenvalues are the solutions of the characteristic equation:(54)X5λ=λ5+D1λ4+D2λ3+D3λ2+D4λ+D5=0,where
(55)D1=β2+c2+c6−c3+c4+c1,D2=β2φ+c4β2−c3−β2−c3−c4c2+c6−c1+αβ1+nτθ−c2c6−c1c2+c6,D3=αβ1θ+c6+c1nτθ−c2c6−β2φ+c4β2−c3c2+c6−c1−β2−c3−c4αβ1+nτθ−c2c6−c1c2+c6,D4=−β2φ+c4β2−c3αβ1+nτθ−c2c6−c1c2+c6−β2−c3−c4αβ1θ+c6+c1nτθ−c2c6,D5=β2φ+c4β2−c3αβ1θ+c6+c1nτθ−c2c6.Thus, using the Routh-Hurwitz stability criterion, the roots of the equation (54) have negative real parts if the following condition holds.
(56)Di>0,foralli=1,⋯,5,D1D2D3>D32+D12D4,andD1D4−D5D1D2D3−D32−D12D4>D5D1D2−D32+D1D52.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) holds.


### 3.3.2. Global Stability of Disease-Free EquilibriumE0

Let us rewrite the coinfection model (2) as
(57)dXdt=FX,Z,dZdt=GX,Z,GX,0=0.where X=S,R and Z=LT,IT,EC,IC,LTC,ITC. 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
(58)U0=X0,0,where,X0=Λμ,0.The conditionsH1 and H2 below must be met to guarantee global asymptotically stability:(H1) For dX/dt=FX,0, the equilibrium point U0 is globally stable(H2) GX,Z=AZ−G^X,Z,G^X,Z≥0 for X,Z∈Ω, where A=DZGU0,0 is a Metzler matrix (the nondiagonal entities of A are nonnegative), and Ω is the biological feasible region of the modelFrom the coinfection model (2), we have
(59)dXdt=FX,Z=Λ−λT+λC+μSωLT+πEC+γIT+ψIC+1−p+qσLTC+1−m+nτITC−μR.Hence,(60)FX,0=Λ−μS0,(61)dZdt=GX,Z=λTS−μ+α+ηλC+ωLTαLT+qσLTC+nτITC−μ+θ+δT+γITλCS−μ+ϵλT+φ+πECφEC+pσLTC+mτITC−μ+δC+ν+ψICηλCLT+ϵλTEC−μ+δC+ρ+σLTCρLTC+θIT+νIC−μ+δTC+τITC.Therefore,(62)A=DZGU0,0=−c1β1000β1α−c200qσnτ00β2−c3β2β2β200φ−c4pσmτ0000−c500θ0νρ−c6whichisaMetzlerMatrix.Here,G^X,Z=AZ−GX,Z, and hence
(63)G^X,Z=G^1X,ZG^2X,ZG^3X,ZG^4X,ZG^5X,ZG^6X,Z=ηλCLT+β1IT+ITC−λTS0β2EC+IC+LTC+ITC+ϵλTEC−λCS0−ηλCLT−ϵλTEC0.Thus,G^5X,Z=−ηλCLT−ϵλTEC<0. This shows that condition (H2) is not satisfied. As a result, U0 and then the disease-free equilibrium point E0 cannot be globally asymptotically stable.


## 3.3.1. Local Stability Analysis of the Disease-Free Equilibrium PointE0

The local stability analysis of the disease-free equilibrium point for the coinfection model is analyzed using the method of linearization. The Jacobean matrix of the model equation (2) at disease-free equilibrium point E0 represented as
(49)JE0=−μ0−β1−β2−β2−β2−β1−β200−c1β1000β100α−c200qσnτ0000β2−c3β2β2β20000φ−c4pσmτ000000−c50000θ0νρ−c600ωγπψ1−p+qσ1−m+nτ−μ.The local asymptotic stability ofE0 is analyzed based on the sign of the eigenvalues. Expanding the characteristic equation ∣λI8−JE0∣=0 with the first and eighth columns, and the sixth row, we obtain three eigenvalues λ1,2=−μ and λ3=−c5. Further, we obtain the remaining five eigenvalues from the characteristic equation of the reduced matrix:
(50)J5−λI5=−c1−λβ100β1α−c2−λ00nτ00β2−c3−λβ2β200φ−c4−λmτ0θ0ν−c6−λ.Here, the characteristic polynomial is computed using elementary row operations method [38]. Thus, the row reduced form of the matrix becomes
(51)α−c2−λ00nτ0θ0ν−c6−λ00φ−c4−λmτ000β2+β2−c3−λc4+λφβ2+β2−c3−λc4+λφ0000β11+c6+λθ+c2+λαnτ−c2+λc6+λθ.Hence, the characteristic polynomial of our matrix is the product of the diagonal entries times−1k, where k is the number of row swaps. We perform three row swaps, so that the characteristic polynomial is
(52)X5λ=−αθφβ2+β2−c3−λc4+λφβ11+c6+λθ+c2+λαnτ−c2+λc6+λθ.Simplifying it gives(53)X5λ=−β2φ+c4β2−c3+β2−c3−c4λ−λ2,αβ1θ+c6+c1nτθ−c2c6+αβ1+nτθ−c2c6−c1c2+c6λ+c2+c6−c1λ2−λ3.That is, the eigenvalues are the solutions of the characteristic equation:(54)X5λ=λ5+D1λ4+D2λ3+D3λ2+D4λ+D5=0,where
(55)D1=β2+c2+c6−c3+c4+c1,D2=β2φ+c4β2−c3−β2−c3−c4c2+c6−c1+αβ1+nτθ−c2c6−c1c2+c6,D3=αβ1θ+c6+c1nτθ−c2c6−β2φ+c4β2−c3c2+c6−c1−β2−c3−c4αβ1+nτθ−c2c6−c1c2+c6,D4=−β2φ+c4β2−c3αβ1+nτθ−c2c6−c1c2+c6−β2−c3−c4αβ1θ+c6+c1nτθ−c2c6,D5=β2φ+c4β2−c3αβ1θ+c6+c1nτθ−c2c6.Thus, using the Routh-Hurwitz stability criterion, the roots of the equation (54) have negative real parts if the following condition holds.
(56)Di>0,foralli=1,⋯,5,D1D2D3>D32+D12D4,andD1D4−D5D1D2D3−D32−D12D4>D5D1D2−D32+D1D52.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) holds.


## 3.3.2. Global Stability of Disease-Free EquilibriumE0

Let us rewrite the coinfection model (2) as
(57)dXdt=FX,Z,dZdt=GX,Z,GX,0=0.where X=S,R and Z=LT,IT,EC,IC,LTC,ITC. 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
(58)U0=X0,0,where,X0=Λμ,0.The conditionsH1 and H2 below must be met to guarantee global asymptotically stability:(H1) For dX/dt=FX,0, the equilibrium point U0 is globally stable(H2) GX,Z=AZ−G^X,Z,G^X,Z≥0 for X,Z∈Ω, where A=DZGU0,0 is a Metzler matrix (the nondiagonal entities of A are nonnegative), and Ω is the biological feasible region of the modelFrom the coinfection model (2), we have
(59)dXdt=FX,Z=Λ−λT+λC+μSωLT+πEC+γIT+ψIC+1−p+qσLTC+1−m+nτITC−μR.Hence,(60)FX,0=Λ−μS0,(61)dZdt=GX,Z=λTS−μ+α+ηλC+ωLTαLT+qσLTC+nτITC−μ+θ+δT+γITλCS−μ+ϵλT+φ+πECφEC+pσLTC+mτITC−μ+δC+ν+ψICηλCLT+ϵλTEC−μ+δC+ρ+σLTCρLTC+θIT+νIC−μ+δTC+τITC.Therefore,(62)A=DZGU0,0=−c1β1000β1α−c200qσnτ00β2−c3β2β2β200φ−c4pσmτ0000−c500θ0νρ−c6whichisaMetzlerMatrix.Here,G^X,Z=AZ−GX,Z, and hence
(63)G^X,Z=G^1X,ZG^2X,ZG^3X,ZG^4X,ZG^5X,ZG^6X,Z=ηλCLT+β1IT+ITC−λTS0β2EC+IC+LTC+ITC+ϵλTEC−λCS0−ηλCLT−ϵλTEC0.Thus,G^5X,Z=−ηλCLT−ϵλTEC<0. This shows that condition (H2) is not satisfied. As a result, U0 and then the disease-free equilibrium point E0 cannot be globally asymptotically stable.


## 3.4. 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 numbersR0C in terms of R0T (and vice versa) [40]. Expressing the parameter μ in the equation (33) in terms of R0T, we have
(64)R0T=β1αμ+k1μ+k2,wherek1=α+ω,k2=δT+γ.Solving forμ and simplifying, we obtain
(65)μ=R0Tk1−k22R0T+4αβ1−k1+k2R0T2R0T.Substituting this value ofμ in R0C (13), we obtain R0C in terms of R0T as
(66)R0C=β2R0Th1+h2/2h2/2+φ+ψR0T−k1+k2R0T/2h2/2+δC+ψR0T−k1+k2R0T/2,where h1=R0TδC+φ+ψ−k1+k2/2, and h2=R0Tk1−k22R0T+4αβ1.The partial derivative ofR0C with respect to R0T is then given by
(67)∂R0C∂R0T=β2R0Th1h2+h3+h4h2k1+k2−k1−k22R0T−2αβ12h2h2/2+φ+ψR0T−k1+k2R0T/22h2/2+δC+ψR0T−k1+k2R0T/22,with h3=φ2+π+δC+ψφ+δC+ψ2R0T2, and
(68)h4=R0Tk12+k22R0T+2αβ12−k1+k2δC+φ+ψR0T.Here, we observe that∂R0C/∂R0T is positive. This implies that the expansion of TB infection exacerbates the pandemic of COVID-19.Remark 8.
If∂R0C/∂R0T=0, the expansion of TB in the community has no significant impact on the spread of COVID-19. In the contrary, when ∂R0C/∂R0T<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 expressingR0T in terms of R0c and identifying the sign of the partial derivative of R0T with respect to R0c.


## 4. 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.


### 4.1. 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 S0=5000, EC0=1000, IC0=10, R0=1, and the parameter values described in each figure.Figure2 presents the solution of the COVID-19 submodel using the parameter values described in the caption. The basic reproduction number is obtained as R0C=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 2
Solutions of the COVID-19 submodel with parameter valuesΛ=500, β2=0.054, μ=0.0477, φ=0.46, δC=0.0023ψ=0.2065, ω=0.09, and its basic reproduction number R0C=0.2524<1.Figure 3
The solution curves of the susceptible, infected, and recovered populations for different initial values of the COVID-19 only model withR0C=0.2524<1.
(a)
Susceptible(b)
Infected(c)
RecoveredFigure4 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.Figure 4
Solutions of the COVID-19 only model with parameter valuesΛ=500, β2=0.41659, μ=0.0477, φ=0.096, δC=0.0023ψ=0.0650, ω=0.09, and its basic reproduction number R0C=3.2707>1.Figure 5
The stable solution curves of the susceptible, infected, and recovered populations for different initial values of the COVID-19 only model use parameter values of the Figure4 in which its basic reproduction number is greater than unity (R0C=3.2707).
(a)
Susceptible(b)
Infected(c)
RecoveredFigures6 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.Figure 6
The effect of the transmission rate (β2) for the COVID-19 only model.
(a)
Exposed(b)
Infected(c)
RecoveredFigure 7
The effect of the recovery rate of infected individuals (ψ) for the COVID-19 only model.
(a)
Exposed(b)
Infected(c)
Recovered


### 4.2. Numerical Solutions of the TB Submodel

The numerical solutions of the TB submodel with the basic reproduction number ofR0T=2.263 are given in Figure 8. The solutions are executed using the initial conditions of the state variables S0=5000, LT0=100, IT0=10, and R0=1. We observe that all solution curves approaches to their endemic equilibrium point. This supports the analytical results that corresponds to R0T>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).Figure 8
Solutions of the TB-19 model with parameter valuesΛ=500, β1=0.6, μ=0.0477, γ=0.16, δT=0.01, α=0.632, ω=0.09, and its basic reproduction number is R0T=2.263>1.Figure 9
The impact of the contact rate of individuals (β1) for the TB model, with other parameter values described in Figure 8.
(a)
Latent TB(b)
Infected TB(c)
Recovered


### 4.3. 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.Table 2
The initial values of the state variables in the model equation (2).
State variableSLTITECICLTCITCRInitial values100001002103323Table 3
The parameter values of the model equation (2) for our numerical solution.
ParameterValueSourceParameterValueSourceΛ500Assumedμ0.0477[14]β10.6[17]β20.659[14]φ0.02Assumedθ0.001AssumedδC0.023[14]δT0.01[17]ψ0.5Assumedα0.25[17]η0.03Assumedσ0.01Assumedϵ0.03Assumedν0.002Assumedτ0.003Assumedδtc0.2Assumedρ0.02Assumedπ0.05[14]ω0.09[32]γ0.516AssumedThe 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 10
Solutions of the coinfection model (2), with parameter values given in Table 3.Figure 11
The dynamics ofITC with the values of τ given in the legend and other parameter values are given in Table 3.Figure12 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 IT, IC, LTC, and ITC increase as the transmission coefficients increase.Figure 12
Impacts ofτ, ν, and θ on the coinfected classes (LTC+ITC) and on the COVID-19 infected class (IC).
(a)
Coinfected(b)
Infected COVID-19Figure 13
Impacts of contact ratesβ1 and β2 on the transmission dynamics of TB (IT) infected, COVID-19 (IC) infected, and coinfected ones (ITC).
(a)
Coinfected(b)
Infected TB(c)
Infected COVID-19


## 4.1. 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 S0=5000, EC0=1000, IC0=10, R0=1, and the parameter values described in each figure.Figure2 presents the solution of the COVID-19 submodel using the parameter values described in the caption. The basic reproduction number is obtained as R0C=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 2
Solutions of the COVID-19 submodel with parameter valuesΛ=500, β2=0.054, μ=0.0477, φ=0.46, δC=0.0023ψ=0.2065, ω=0.09, and its basic reproduction number R0C=0.2524<1.Figure 3
The solution curves of the susceptible, infected, and recovered populations for different initial values of the COVID-19 only model withR0C=0.2524<1.
(a)
Susceptible(b)
Infected(c)
RecoveredFigure4 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.Figure 4
Solutions of the COVID-19 only model with parameter valuesΛ=500, β2=0.41659, μ=0.0477, φ=0.096, δC=0.0023ψ=0.0650, ω=0.09, and its basic reproduction number R0C=3.2707>1.Figure 5
The stable solution curves of the susceptible, infected, and recovered populations for different initial values of the COVID-19 only model use parameter values of the Figure4 in which its basic reproduction number is greater than unity (R0C=3.2707).
(a)
Susceptible(b)
Infected(c)
RecoveredFigures6 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.Figure 6
The effect of the transmission rate (β2) for the COVID-19 only model.
(a)
Exposed(b)
Infected(c)
RecoveredFigure 7
The effect of the recovery rate of infected individuals (ψ) for the COVID-19 only model.
(a)
Exposed(b)
Infected(c)
Recovered


## 4.2. Numerical Solutions of the TB Submodel

The numerical solutions of the TB submodel with the basic reproduction number ofR0T=2.263 are given in Figure 8. The solutions are executed using the initial conditions of the state variables S0=5000, LT0=100, IT0=10, and R0=1. We observe that all solution curves approaches to their endemic equilibrium point. This supports the analytical results that corresponds to R0T>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).Figure 8
Solutions of the TB-19 model with parameter valuesΛ=500, β1=0.6, μ=0.0477, γ=0.16, δT=0.01, α=0.632, ω=0.09, and its basic reproduction number is R0T=2.263>1.Figure 9
The impact of the contact rate of individuals (β1) for the TB model, with other parameter values described in Figure 8.
(a)
Latent TB(b)
Infected TB(c)
Recovered


## 4.3. 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.Table 2
The initial values of the state variables in the model equation (2).
State variableSLTITECICLTCITCRInitial values100001002103323Table 3
The parameter values of the model equation (2) for our numerical solution.
ParameterValueSourceParameterValueSourceΛ500Assumedμ0.0477[14]β10.6[17]β20.659[14]φ0.02Assumedθ0.001AssumedδC0.023[14]δT0.01[17]ψ0.5Assumedα0.25[17]η0.03Assumedσ0.01Assumedϵ0.03Assumedν0.002Assumedτ0.003Assumedδtc0.2Assumedρ0.02Assumedπ0.05[14]ω0.09[32]γ0.516AssumedThe 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 10
Solutions of the coinfection model (2), with parameter values given in Table 3.Figure 11
The dynamics ofITC with the values of τ given in the legend and other parameter values are given in Table 3.Figure12 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 IT, IC, LTC, and ITC increase as the transmission coefficients increase.Figure 12
Impacts ofτ, ν, and θ on the coinfected classes (LTC+ITC) and on the COVID-19 infected class (IC).
(a)
Coinfected(b)
Infected COVID-19Figure 13
Impacts of contact ratesβ1 and β2 on the transmission dynamics of TB (IT) infected, COVID-19 (IC) infected, and coinfected ones (ITC).
(a)
Coinfected(b)
Infected TB(c)
Infected COVID-19


## 5. 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.



---

*Source: 2449710-2022-03-27.xml*