Co-dynamics of COVID-19 and TB with COVID-19 vaccination and exogenous reinfection for TB: An optimal control application

COVID-19 and Tuberculosis (TB) are among the major global public health problems and diseases with major socioeconomic impacts. The dynamics of these diseases are spread throughout the world with clinical similarities which makes them difficult to be mitigated. In this study, we formulate and analyze a mathematical model containing several epidemiological characteristics of the co-dynamics of COVID-19 and TB. Sufficient conditions are derived for the stability of both COVID-19 and TB sub-models equilibria. Under certain conditions, the TB sub-model could undergo the phenomenon of backward bifurcation whenever its associated reproduction number is less than one. The equilibria of the full TB-COVID-19 model are locally asymptotically stable, but not globally, due to the possible occurrence of backward bifurcation. The incorporation of exogenous reinfection into our model causes effects by allowing the occurrence of backward bifurcation for the basic reproduction number R0 < 1 and the exogenous reinfection rate greater than a threshold (η > η∗). The analytical results show that reducing R0 < 1 may not be sufficient to eliminate the disease from the community. The optimal control strategies were proposed to minimize the disease burden and related costs. The existence of optimal controls and their characterization are established using Pontryagin's Minimum Principle. Moreover, different numerical simulations of the control induced model are carried out to observe the effects of the control strategies. It reveals the usefulness of the optimization strategies in reducing COVID-19 infection and the co-infection of both diseases in the community.


Introduction
Coronavirus disease 2019  and Tuberculosis (TB) are two major infectious diseases that pose a significant public health threat, and their co-infection exacerbates the situation. Both are contagious diseases that attack primarily the lungs. They share plenty of common symptoms such as fever, cough, and difficulty breathing. However, TB have been known with a longer incubation period with a slower onset of the disease (WHO, 2020b). COVID-19 is a communicable disease caused by a family of novel coronavirus called SARS-CoV-2 while TB is a respiratory pathogens caused by Mycobacterium Tuberculosis (MTB), which is transmitted between humans through the respiratory route such as cough, sneeze, speak, kiss or spit and most commonly affects the lungs, but can damage any tissue. TB can spread from active pulmonary TB persons but not latent TB. Worldwide, TB remains a leading cause of death. According to the World Health Organization Global TB Report 2021, approximately 10 million people contracted TB and 1.5 million died from TB in 2020, worldwide (WHO, 2021). COVID-19, on the other hand, about 480 million confirmed cases were reported up to the end of March 2022 and more than 6 million people were died due to COVID-19 worldwide (WHO, 2023). Currently, there is an emerging evidence that patients with latent TB and active TB disease are more likely to get SARS-CoV-2 and become severely infected with COVID-19 (Chen et al., 2020).
In the study of infectious diseases, co-epidemic requires critical understanding of how the diseases are related, prevalence, their prevention and treatment to effectively control these diseases. A co-epidemic occurs when the spread of one infectious disease stimulates that of another. This is just recently occurred with COVID-19 and TB. Different studies have been conducted on the co-interactions of COVID-19 and TB diseases (Chen et al., 2020;Davies, 2020;Motta et al., 2020;Mousquer, Peres, & Fiegenbaum, 2021;Petrone et al., 2021). In particular, Chen et al. (2020), observed that TB can increase the susceptibility to COVID-19 and the severity of its symptoms. In the study by Petrone et al. (2021), it is observed that peoples infected with COVID-19 and TB develop a weak immune response to SARS-COV-2. The work done in Mousquer et al. (2021) shows that COVID-19 and TB co-infected peoples are associated with higher morbidity and deaths. Besides, it is reported that patients having the past history of TB infection have high risk of mortality if they acquired COVID-19 (Davies, 2020). In another work done by Chen et al. (2020), patients having latent and active TB infection would have higher risk of COVID-19 infection with its severity.
Mathematical modeling has become the powerful tool for studying the behaviour of infectious diseases (Egonmwan & Okuonghae, 2019;Omame, Okuonghae, Umana, & Inyama, 2020;Omame, Okuonghae, Nwafor, & Odionyenma, 2021). Different mathematical modeling of COVID-19 with vaccination is studied in (Gumel, Iboi, Ngonghala, & Ngwa, 2020;L. Zhang, Ullah, Al Alwan, Alshehri, & Sumelka, 2021;Diagne, Rwezaura, Tchoumi, & Tchuenche, 2021;Yavuz, Coş ar, Günay, & € Ozdemir, 2021). For instance, Gumel et al. (2020) present a mathematical analysis for the effects of vaccination and nonpharmaceutical interventions on COVID-19 dynamics. Their findings suggests that, adding the impacts of the therapeutic benefits of the vaccination against COVID-19 into the model resulted in a dramatic elimination of the pandemic. Diagne et al. (2021) developed and analyzed a mathematical model of COVID-19 with vaccination and treatment. Their results show the effectiveness of vaccination and treatment in reducing the spread of COVID-19, and more efforts are needed to eliminate the disease. They concluded that the combined or simultaneous use of non-pharmaceutical protective measures such as face masks, hand washing, and social distancing should continue to be encouraged. Furthermore, the effects of vaccination and other constant and time-dependent variable control measures on the dynamics of COVID-19 is investigated in (L. Zhang et al., 2021). Their findings revealed that, compared to constant controls, the time-depended control measures play a vital role in disease eradication. Moreover, mathematical models for the dynamics of COVID-19 and its co-infection with other diseases have been also developed (Hezam, Foul, & Alrasheedi, 2021;Omame, Sene, et al., 2021;Tchoumi, Diagne, Rwezaura, & Tchuenche, 2021).
Optimal controls have many applications in dynamical systems including systems governed by nonlinear ordinary differential equations. It can be seen as a control strategy in control theory (Ross, 2015). Mathematical models with optimal control analysis have become an important tool for understanding disease transmission dynamics and in the process of making decision regarding disease control. For example in the study of Malaria (Cai, Li, Tuncer, Martcheva, & Lashari, 2017;Makinde & Okosun, 2011;Mohammed-Awel & Gumel, 2019), TB (Das, Khajanchi, & Kar, 2020;Kim, Aurelio, & Jung, 2018), COVID-19 (Lemecha Obsu & Feyissa Balcha, 2020;Mondal & Khajanchi, 2022), HIV-TB coinfection (Aggarwal et al., 2020a;Silva & Torres, 2015), Malaria-COVID-19 coinfection , COVID-19 and Diabetes coinfection (Omame, Sene, et al., 2021) and the references there in. Recently, some mathematical modeling of COVID-19 and TB co-infection with optimal control have been proposed and analyzed (Goudiaby et al., 2022;Mekonen, Obsu, & Habtemichael, 2022;Rwezaura, Diagne, Omame, de Espindola, & Tchuenche, 2022). In particular, Goudiaby et al. (2022) developed and studied a deterministic compartmental model for the transmission dynamics of tuberculosis and COVID-19 with optimal control. They suggest that the best strategy of interest to public health policy and decision makers to contain the spread of COVID-19 is to focus on prevention, treatment, and control of co-infection with COVID-19 which also results a lower percentage of total cost in COVID-19 prevention. In Mekonen et al. (2022), the authors have proposed an optimal control model of COVID-19 and TB by introducing four control measures (preventions and treatments for both diseases) to optimally manage the diseases. They showed that the prevalence of the co-infection reduced when concurrently all the control measures were implemented. On the other hand, Rwezaura et al. (2022) proposed the co-infection model for SARS-CoV-2 and TB that incorporates the vaccinated population with optimal control strategies to manage the spread of these two diseases. Their results suggest that either COVID-19 or TB control strategies will significantly mitigate a significant number of new co-infected cases.
However, the spread of COVID-19 and TB are still continuing and claiming more lives and thus, their co-dynamics need further investigation. Based on the descriptions above, the authors in this study are motivated in developing a mathematical model on the co-dynamics of COVID-19 and TB by incorporating COVID-19 vaccination and exogenous reinfection for TB. The importance of considering exogenous reinfection is justified in (Bandera et al., 2001;Caminero et al., 2001;Feng, Castillo-Chavez, & Capurro, 2000;van Rie et al., 1999). We also extend the model to an optimal control model by incorporating time variant controls to mitigate the spread of these two diseases. Compared to the previous COVID-19 and TB co-infection models, our model differ by several ways. For example, compared to the model proposed in (Goudiaby et al., 2022;Mekonen et al., 2022), our model considers the vaccination against COVID-19, TB re-infection after recovery, and exogenous TB re-infection. Moreover, compared to the model given in , our model considers the exogenous TB re-infection.
The rest of this paper has been organized as follows. In Section 2, the details on the formulation of the proposed coinfection model are presented. The invariant region where the model is biological relevant is presented in Section 2.1. Sections 3.1 and 3.2 present and theoretically examine the COVID-19 and TB sub-models, respectively. In Section 3.3, the full COVID-19 and TB co-infection model is analyzed. In Section 4, the extended optimal control model is presented and then analyzed using the well known Pontryagin's Maximum Principle. Moreover, different numerical simulations of control induced model are provided in Section 5. Finally, the paper concludes in Section 6.

Model formulation
In this section, the deterministic compartmental model of the transmission dynamics of both diseases are described. Depending on individuals' disease status, at any time t, the total population N(t) is divided into nine compartments as described in the following Table 1.
The total population is given by NðtÞ ¼ SðtÞ þ VðtÞ þ I c ðtÞ þ R c ðtÞ þ LðtÞ þ I t ðtÞ þ R t ðtÞ þ I cL ðtÞ þ I tc ðtÞ: Individuals in the susceptible class are reduced at the rate, N when infected with COVID-19, and are also reduced at the rate when infected with TB. The forces of infection considered here are of standard incidence type (J. Zhang & Ma, 2003). The parameters b c and b t respectively represent the effective contact rate of COVID-19 and TB infected. The parameter t(t ! 1) accounts for the increased infectivity of co-infected individuals as a result of TB (Aggarwal & Raj, 2021;Rwezaura et al., 2022).
The factor (1 À kr) represents the effect of personal protection against COVID-19 (such as wearing a facial mask, physical distancing, hand washing with soap), where 0 < r < 1 is the efficacy of adopted protection strategy, and 0 < k < 1 is compliance or the fraction of community employing it.
COVID-19 infected individuals may acquire TB infection at the rate ul t and enter into the class I cL . The modification parameter u (u ! 1) is the enhancement factor accounting for the relative infectiousness of individuals acquiring TB following COVID-19 infection. Every individual after acquiring a TB infection, initially moves to the latently infected TB class, that is, in L or I cL , and then enters the active TB class (I t or I tc ) at a progression rate 4 t and h c , respectively. In addition, the latent TB infectives also move to the active TB infectious class, as exogenous reinfection occurs at a rate hl t as a result of recent exposure of latent TB infectives with actively infected individuals. Individuals after recovery from TB again become susceptible to both diseases. They become COVID-19 infected at the rate l c and acquire TB infection at the rate 4l t . Such a transition is used in most TB and its co-dynamics models, where 4(0 4 1) describes TB re-infection rate (Aggarwal et al., 2020a(Aggarwal et al., , 2020bAggarwal & Raj, 2021;Rwezaura et al., 2022). Individuals recovered from COVID-19 gains immunity against re-infection.
However, there is possibility of infection with TB at the rate l t . Latent TB infected individuals may acquire COVID-19 infection at the rate l c and enter into the class I cL . Also, active TB infected individuals may acquire COVID-19 infection at the rate sl c and enter into the class I tc . The modification parameter s(s ! 1) is represents the increasing factors for the relative infectiousness of individuals acquiring COVID-19 following TB infection. The parameters h L and q represents COVID-19 recovery rate for co-infected individuals in I cL and I tc classes, respectively. All other transitions in the model are described by the model flow diagram shown in Fig. 1. The definitions of the model parameters are provided in Table 2 together with their values. From the aforementioned, the formulated model is (1) with initial conditions Sð0Þ ! 0; Vð0Þ ! 0; I c ð0Þ ! 0; R c ð0Þ ! 0; Lð0Þ ! 0; I t ð0Þ ! 0; R t ð0Þ ! 0; I cL ð0Þ ! 0; I tc ð0Þ ! 0: The parameter h in the rate of exogenous reinfection term hl t L represents the reinfection level. When h ¼ 0, model (1) reduces to the coinfection model with out exogenous reinfection. A value of h 2 (0, 1) implies that reinfection is less likely than a new infection.

Positivity and boundedness of solutions
Since each component of the given model (1) represents the human population, we need to show all variables S(t), V(t), I c (t), R c (t), L(t), I t (t), R t (t), I cL (t) and I tc (t) are positive for all time t ! 0. Therefore, we state the following theorem: Theorem 1. All solutions of model (1) with nonnegative initial conditions remain nonnegative for all t > 0.
Proof. For proving the positivity of all the components of the model (1), we assume on the contrary that, there exists a first time t 1 , such that minfXðt 1 Þg ¼ 0 and minfXðtÞg > 0; for all t2½0; t 1 Þ: Here, X(t) ¼ S(t), V(t), I c (t), R c (t), L(t), I t (t), R t (t), I cL (t), I tc (t). Following our assumption, we first let, minfXðt 1 Þg ¼ Sðt 1 Þ: Hence, S(t 1 ) ¼ 0 and S(t) > 0 for all t2½0; t 1 Þ. But, from which we get, which contradicts our assumption S(t 1 ) ¼ 0. Hence, S(t) > 0 for all t ! 0. Similarly, we can prove that all solution components are positive for t ! 0 in all other cases. ,

Invariant region
Theorem 2. The closed region U ¼ ðS; V; I c ; R c ; L; I t ; R t ; I cL ; I tc Þ2R 9 þ : NðtÞ L m : is a positively invariant set for the system (1) and attracts all positive solutions.
Proof. The rate of change of total population for the model (1) is dN dt Thus, it follows that NðtÞ Nð0Þe Àmt þ L m ð1 À e Àmt Þ: More specifically, 0 < NðtÞ L m , if Nð0Þ L m . That is, N(t) is bounded and all solutions beginning in U approach, enter, or remain in U. As a result, the model in (1) can be regarded as being well-posed. ,

COVID-19 sub-model
The COVID-19 sub-model (3) will be studied in the following feasible region Analogous to Theorem 2, U C can be easily proved to be positively invariant.

Basic reproduction numbers
The COVID-19 sub-model (3) has a disease-free equilibrium (DFE) given by We can calculate the basic reproduction number using the next generation matrix operator (Van den Driessche & Watmough, 2002). Thus, following the same approach as in Appendix A, the vaccination-induced reproduction number of the COVID-19 sub-model (3) is given by In the absence of COVID-19 vaccination (n ¼ 0), the basic reproduction number is given by Note that where j ¼ ½mþð1ÀεÞn nþm . Since j < 1, as expected R vC < R 0C . The value of j becomes small with high value of vaccine efficacy (ε2ð0; 1 ). The parameter j which depends on the rate of vaccination and vaccine efficacy represents the effect of COVID-19 vaccine implementation in reducing the vaccine-induced reproduction number.

Stability analysis of DFE
The local stability of DFE ðE 0 c Þ of COVID-19 sub-model (3) is determined using the eigenvalues of the Jacobian matrix at E 0 c , which is given by The corresponding eigenvalues of the Jacobian matrix J E 0 c are l 1 ¼ Àðn þ mÞ; l 2 ¼ Àm; l 3 ¼ ð1 À krÞb c ðm þ ð1 À εÞnÞ n þ m À ðr c þ m þ d 1 Þ; and l 4 ¼ Àm: We can rewrite the eigenvalue l 3 as Since the eigenvalues l 1 , l 2 and l 4 have negative real parts, the local stability of the DFE, ðE 0 c Þ depends on the sign of the eigenvalue l 3 . Thus, the following result have been established.

Lemma 3. (Local stability of DFE:)
The DFE E 0 c of the COVID-19 sub-model (3) is locally asymptotically stable if R vC < 1 and unstable if R vC > 1.

Existence of endemic equilibrium for COVID-19 sub-model
For finding the endemic equilibrium point (EE) point for COVID-19 sub-model (3), we must solve the following system of equations at an arbitrary equilibrium point Rearranging the second equation gives us the following quadratic equation It is important to note that the coefficient A of the quadratic (11) is always positive. If R 0C < 1 (i.e., R vC is also less than unity), then B > 0 and C > 0. In this case, since all coefficients of the quadratic equation (11) are positive, from Routh-Hurwitz criterion there is no positive root(s) of (11). If R vC > 1, then C < 0 and B > 0 or B < 0. In this case, we have a unique positive root for the quadratic (11). Hence, the following result can be established: Theorem 4. The COVID-19 sub-model (3) has only one unique endemic equilibrium if R vC > 1.

Bifurcation analysis
Depending on the threshold parameter R vC , the DFE can interchange its stability with the EE. We want to determine the direction of bifurcation (either forward or backward) using the centre manifold theory as described in Theorem 4.1 from (Castillo-Chavez & Song, 2004). The following variable changes and simplifications are first done in order to apply this theory. Let The Jacobian matrix of the system (13) evaluated at E 0 * c Þ has an obvious simple zero eigenvalue, and all other eigenvalues have a negative sign. Hence, E 0 c is a non-hyperbolic equilibrium, when b c ¼ b * c . Thus, the Castillo-Chavez and Song can be used to analyze the dynamics of the model (13) near b c ¼ b * c . By using the notation in Theorem 4.1 given in (Castillo-Chavez & Song, 2004), we proceed as follows. Right eigenvector: Let w ¼ ðw 1 ; w 2 ; w 3 ; w 4 Þ T be the right eigenvector associated to the zero eigenvalue. i.e., Solving for w 1 , w 2 , w 3 and w 4 , we get the right eigenvector w ¼ ðw 1 ; w 2 ; w 3 ; w 4 Þ T such that The eigenvectors w and v should satisfy the requirement v , w ¼ 1, i.e., v 3 w 3 ¼ 1: For finding the direction of the bifurcation, we need to determine the sign of bifurcation coefficients a and b which are given by Since v 1 ¼ 0, v 2 ¼ 0, and v 4 ¼ 0, we do not need the derivatives of f 1 , f 2 and f 4 . Thus, we need only the second order partial derivatives of f 3 where Therefore, from the second order partial derivatives of f 3 at ðE 0 c ; b * c Þ, the only ones that are nonzero are the following: Now, using the above information, we obtain Since a < 0 and b > 0, by Theorem 4.1 stated in (Castillo-Chavez & Song, 2004), the COVID-19 sub-model (3) exhibit the phenomenon of forward bifurcation at R vC ¼ 1. Hence, the following result follows.
The public health implication for the forward bifurcation is that the COVID-19 disease cannot invade the population for R vC < 1. On the contrary, the disease persist in the community when R vC exceeds unity. (1), we obtain the following TB-only model.

TB sub-model
Analogously to Theorem 2, it can be shown that the region (14). Thus, we consider the dynamics of the TB sub-model (14) in U T .

Disease-free equilibrium and basic reproduction number
System (14) always has the disease-free equilibrium The basic reproduction number for system (14) is calculated in Appendix A and is given by This number is given by the product b t /(r t þ m þ d 2 ), that is, by the product of the average number of susceptibles infected by one infectious individual during his or her effective infectious period and 4 t /(4 t þ m), the fraction of the population which survives the latent period. Therefore, R 0T gives the number of secondary infectious cases produced by an infectious individual during his or her effective infectious period in a population of susceptible individuals.
Remark 1. The basic reproductive number R 0T does not depend on h. However, the number of endemic equilibria depends on the value of reinfection level h. The expression of R 0T here is essentially the same as that of TB sub-model without a exogenous reinfection term (h ¼ 0). Here we show that as the value of h changes system (14) exhibits a backward bifurcation at R 0T ¼ 1.
In the next section, we analyze the local asymptotic stability of the DFE E 0 t of TB sub-model (14).

Stability of the disease-free equilibrium
Lemma 6. (Local stability of DFE:) The DFE E 0 t of the TB sub-model (14) is locally asymptotically stable if R 0T < 1, and unstable otherwise.
Proof. The local stability of DFE ðE 0 t Þ of TB sub-model (14) is obtained using the sign of eigenvalues of the Jacobian matrix at E 0 t . The Jacobian matrix at E 0 t of the TB sub-model (14) is given by It is obvious that l 1 ¼ Àm ¼ l 4 are the two negative eigenvalues of J E 0 t . The remaining eigenvalues of J E 0 t are obtained from the block matrix The characteristic polynomial of J 1E 0 t is given by Because all model parameters are positive, a 1 > 0 and a 2 > 0 if and only if R 0T < 1. Therefore, by Routh-Hurwitz criterion, both roots of the characteristic polynomial (16) have negative real part, if R 0T < 1. Thus, all eigenvalues of the Jacobian matrix J E 0 t have a negative real part if R 0T < 1. If R 0T > 1, then a 3 < 0 and the matrix J E 0 t has at least one eigenvalue with positive real part. Therefore, from this we conclude that the DFE E 0 t of the TB sub-model system (14) is locally asymptotically stable if R 0T < 1 and unstable if R 0T > 1. , The implication of Theorem 6 is that TB can be excluded from the community (for R 0T < 1) if the initial size of the model subpopulation lies in the attractive basin of E 0 t . Since some TB models can exhibit the phenomenon of backward bifurcation (see Sharomi, Podder, Gumel, and Song (2008) and the references therein), it is necessary to determine whether the current TB sub-model undergo this phenomenon.

Existence of endemic equilibrium for TB sub-model
The endemic equilibrium point of the TB sub-model (14) is given by such that while l * t is taken from the positive root of the following third degree polynomial Theorem 7. TB sub-model (14) always has an endemic equilibrium whenever R 0T > 1.
Proof. The proof follows easily from Theorem 2 in (Aldila & Angelina, 2021) and Theorem 3.3 in (Khajanchi, Das, & Kar, 2018). , Theorem 7 implies that there always exists at least one EE if R 0T > 1. However, because Pðl * t Þ is a third degree polynomial, we may have multiple endemic equilibria when R 0T > 1 or R 0T < 1. Thus, in light of the Descartes principles of signs, we examine the maximum number of positive roots of the polynomial (19). Table 3 provides a summary of the results. It can be observed that for R 0T < 1 the polynomial Pðl * t Þ will either have zero or two positive roots. On the other hand, when R 0T > 1, we always have the possibility of having either one or three positive roots.
From Table 3, we can see that the TB sub-model (14) may have an EE even though R 0T < 1. Unfortunately, the results in Table 3 cannot provide a specific condition to guarantee the existence of the EE when R 0T < 1. Therefore, we continue the analysis by determining the possible conditions that guarantee the existence of EE when R 0T < 1, which is indicated by the negative sign of vlt vR0T when R 0T ¼ 1 and l t ¼ 0. Furthermore, if it is satisfied we have two positive endemic equilibria if R 0T < 1.
First, we make A, B, C, and D in Pðl * t Þ to make it a function depending on R 0T by changing b t as , C, and D, and taking the partial derivative of l t with respect to R 0T from Pðl * t Þ (i.e., differentiating Equation (19) implicitly with respect to R 0T ) and evaluating it in R 0T ¼ 1, l t ¼ 0, give us: So, we get the following result.
Theorem 8. TB sub-model (14) has two endemic equilibriums in U T when R 0T < 1 if      Based on Theorems 7 and 8, we propose the possibility of multiple EE points for R 0T < 1; even the DFE is locally asymptotically stable. Therefore, we need to understand these phenomena in more detail. The local stability of this EE is detailed in the next section.

Existence of backward bifurcation
In this section, we employ the same strategy as in Section 3.1 to investigate direction of the bifurcation using the wellknown Castillo-Chavez and Song bifurcation theorem (Castillo-Chavez & Song, 2004). To apply this theory, the following simplification and change of variables are made first of all.
Denoting (14) can be rewritten in the form dx dt ¼ f ðxÞ, where f ðxÞ ¼ ðf 1 ðxÞ; f 2 ðxÞ; f 3 ðxÞ; f 4 ðxÞÞ as follows 8 > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > : . Choosing b * t as the bifurcation parameter and setting R 0T ¼ 1, we get The Jacobian matrix of the system (22) evaluated at E 0 t and b t ¼ b * t is given by Obviously, l 1 ¼ Àm ¼ l 4 are the two negative eigenvalues of J ðE 0 t ;b * t Þ and the remaining eigenvalues are the roots of the characteristic equation * t Þ has a simple zero eigenvalue and other eigenvalues have negative sign. Hence E 0 t is a non-hyperbolic equilibrium, when b t ¼ b * t . Next, we calculate the right and left eigenvectors associated to the zero eigenvalue.Right eigenvector: Let w ¼ ðw 1 ; w 2 ; w 3 ; w 4 Þ T be the right eigenvector associated to the zero eigenvalue. i.e., J ðE 0 t ;b * t Þ ,w ¼ 0: Solving for w 1 , w 2 , w 3 and w 4 , we get the right eigenvector w ¼ ðw 1 ; w 2 ; w 3 ; w 4 Þ T such that The eigenvectors w and v should satisfy the requirement v , w ¼ 1, i.e., v 3 w 3 1 þ r t þ m þ d 2 4 t þ m ¼ 1: Next, to find the direction of the bifurcation, we need to determine the sign of bifurcation coefficients a and b defined as Since v 1 ¼ v 4 ¼ 0, we do not need the derivatives of f 1 and f 4 . Thus, we need only the second-order partial derivatives of f 2 and f 3 where Therefore, the nonzero second-order partial derivatives of f 2 and f 3 at ðE 0 Now, using the above information we obtain the bifurcation coefficients a and b as follows.
Following Theorem 4.1 in Castillo-Chavez and Song that the model (14) experiences the phenomena of backward bifurcation at R 0T ¼ 1 whenever both b > 0 and a ¼ À2 Hence, the following conclusion may be drawn.
Theorem 9. The TB sub-model (14) undergoes backward bifurcation at R Note that the condition h* in equation (23) has the same condition as inequality (21). It is clear from Theorem 9 that the backward bifurcation phenomenon will not occur in the absence of re-infection (i.e., h ¼ 0), because the right-hand side of the inequality in Theorem 9 is non-negative. To further confirm that there is no backward bifurcation in model (14) when no reinfection occurs (i.e., h ¼ 0), one can check the global-asymptotic stability for the DFE.
From an epidemiological point of view, the occurrence of backward/subcritical bifurcation implies that having R 0T < 1, although necessary, is no longer sufficient for disease control.
Having analyzed the dynamics of the two sub-models, the full COVID-19-TB model (1) will now be analysed.

Analysis of the full TB-COVID-19 model
In this section, we analyze the full TB-COVID-19 co-infection model given by (1).

Stability of the disease-free equilibrium
The disease-free equilibrium of the full TB-COVID-19 model (1) is given by ; Ln mðn þ mÞ ; 0; 0; 0; 0; 0; 0; 0 : It is easy to show, using the next generation method (as in the previous sections), that the associated reproduction number for the full TB-COVID-19 model system (1) (denoted by R 0 ) is given by where R vC and R 0T are the basic reproduction numbers corresponding to COVID-19 only and TB-only submodels, given by Eqs. (4) and (15), respectively. This implies that the dynamics of the TB-COVID-19 coinfection will be dominated by the disease with the bigger basic reproduction number. (1) is locally asymptotically stable if the threshold parameter R 0 < 1, and unstable if R 0 > 1.

Lemma 10. The DFE of the TB-COVID-19 model
Proof. The Jacobian matrix for the system (1) evaluated at the DFE E 0 is given by We can easily find the eigenvalues of J E 0 from the following block matrices Obviously from the above block matrices, we obtained the eigenvalues l 1 ¼ Àk 1 , l 2 ¼ Àm, l 3 ¼ Àm, l 4 ¼ k 2 (R vC À 1), l 7 ¼ Àm, l 8 ¼ Àk 5 and l 9 ¼ Àk 6 . The remaining eigenvalues of J E 0 are the roots of the characteristic polynomial l 2 þ a 1 l þ a 2 ¼ 0; The first three and the last three eigenvalues of J E 0 have a negative real part. On the other hand, the fourth eigenvalue and the polynomial (25) will give eigenvalues with a negative real part, if R 0 ¼ maxfR 0C ; R 0T g < 1. Hence, E 0 is locally asymptotically stable if R 0 < 1 and unstable otherwise. , The full TB-COVID-19 model (1) will experience backward bifurcation under the same conditions as the TB sub-model (14), which may experience the phenomena of backward bifurcation (Mtisi, Rwezaura, & Tchuenche, 2009;Mukandavire, Gumel, Garira, & Tchuenche, 2009;Okosun & Makinde, 2013;Tchoumi et al., 2021). Applying the centre manifold theory described in (Castillo-Chavez & Song, 2004), the following result is established.
Theorem 11. If the bifurcation parameters a > 0 and b > 0, the full model (1) will exhibit the phenomenon of backward bifurcation at R 0 ¼ 1.
For the proof, one can follow the same approach as in Section 3.2.4. It should be noted that the occurrence of backward bifurcation prevents the model's equilibria from having global asymptotic stability. That is, both the DFE and endemic equilibria can only be locally asymptotically stable.
In the next section, we formulate and analyze the optimal control problem.

The optimal control model
To examine the possible effects of applying intervention strategies to slow the spread of COVID-19, TB, and their coinfection, we incorporate the time dependent controls in the model (1). We take into account the following model equations for the optimal control problem.
where (i) u 1 (t) represents COVID-19 prevention strategies (such as face mask, hand washing, physical distancing), (ii) u 2 (t) represents time-dependent COVID-19 vaccination rate, (iii) u 3 (t) denotes TB treatment, and l c and l t are as in Section 2 with initial condition given in (2).
The controls u i (,) are bounded between 0 and u i max , with u i max < 1, i ¼ 1, 2, 3. No additional preventative measures are taken to reduce either of diseases when the controls vanish. We assume that u i max is never equal to 1, since it makes the model more realistic from a medical point of view and thus we assumed u i max 1 À d, i ¼ 1, 2, 3, where d ≪ 1 denotes a positive real number. We denote u ¼ (u 1 , u 2 , u 3 ) for our conveniences and later use.
Here, the target is to minimize disease burden and related costs. That is, to reduce the cost of COVID-19 prevention, the cost of COVID-19 vaccination, and the cost of TB treatment for active TB cases. Hence, we formulate the cost functional for the minimization problem as subject to (26) with its initial conditions. The positive coefficients A 1 , A 2 , …, A 4 are the weight constants associated with COVID-19 infection, TB infection, COVID-19 and latent TB co-infection, and COVID-19 and active TB co-infection, respectively. The parameters W 1 , W 2 , and W 3 represent the weight constants that measure the relative cost of the interventions connected to the control u 1 , u 2 , and u 3 , respectively. Larger the values of W 1 , W 2 and W 3 , more expensive is the implementation cost for detection. The fixed constant t f denotes the final intervention time. Moreover, the terms A 1 I c (t), A 2 I t (t), A 3 I cL (t) and A 4 I tc (t) denote, respectively, the cost incurred due to COVID-19 infection, TB infection, COVID-19 and latent TB co-infection, and COVID-19 and active TB co-infection. We take a quadratic form for the cost of control as an approximation to the real nonlinear functional that depends on the assumption that the cost takes a nonlinear form (Kumar, Srivastava, Dong, & Takeuchi, 2020). As a result, we prevent cases of singular or bang-bang optimal control (Fleming & Rishel, 2012;Pontryagin, 1987).
Here, our goal is to find the optimal control u * ¼ ðu * 1 ; u * 2 ; u * 3 Þ of the control functions u(t) ¼ (u 1 (t), u 2 (t), u 3 (t)) such that the associated state trajectories are the solution of the system (26) is an admissible control set. The Pontryagin's Minimum Principle (Fleming & Rishel, 2012) helps to reduce problems (26), (27) and (29) to a problem of minimizing the Hamiltonian H given as: x i ðtÞg i ðt; S; V; I c ; R c ; L; I t ; R t ; I cL ; I tc ; uÞ where g i represents the set of right-hand side functions in (26) and x : ½0; t f /R 9 , x ¼ (x 1 (t), x 2 (t), …x 9 (t)) is a nontrivial absolutely continuous mapping called adjoint vector. In the following, the existence and characterization of the optimal controls are detailed.

Existence of the optimal control
In this subsection, the existence of optimal controls that minimize the cost function in the finite time is discussed. In this case, the boundedness of solutions to control-induced system (26) is very critical in determining the existence and uniqueness of optimal control (Abraha, Al Basir, Obsu, & Torres, 2021).
Theorem 12. Given the objective functional (27) subjected to the control-induced initial value problem (26), there exists optimal control triple u * ¼ ðu * 1 ; u * 2 ; u * 3 Þ in U and corresponding state solutions (S*, V*, I * c , R * c , L*, I * t , R * t , I * cL , I * tc ) such that J ðu * Þ ¼ min U J ðuÞ: Proof. For the existence of the given optimal control u*, we need to verify the following conditions given by Fleming and Rishel (2012).
(i) The set of solutions to system (26) with associated control functions in (29) is nonempty. (ii) The control admissible set U (29) is convex and closed. (iii) Each right-hand side of the state system (26) is continuous, bounded by a sum of the bounded control and the state. (iv) The state system can be written as a linear function of control variables u ¼ (u 1 , u 2 , u 3 ), with the coefficients as functions of time and state variables. (v) The integrand F(t, S, V, I c , R c , L, I t , R t , I cL , I tc , u) in equation (27) is convex with respect to control variables.
(vi) There exist positive constants c 1 , c 2 , c 3 , c 4 and a constant c > 1 such that Fðt; S; V; I c ; R c ; L; I t ; R t ; I cL ; I tc ; uÞ ! c 2 ju 1 j c þ c 3 ju 2 j c þ c 4 ju 3 j c À c 1 The state variables in the proposed model are continuously differentiable. Besides, it is proved in Theorem 2 that the solutions of the state system are continuous and bounded. Furthermore, the right-hand side functions of model (26) satisfy the Lipschitz condition with respect to the state variables, which is follows from the boundedness of the partial derivatives with respect to state variables in the state system. Therefore, the result follows from Picard-Lindel€ of's existence theorem (Coddington & Levinson, 1955). Being a quadratic function of u ¼ (u 1 , u 2 , u 3 ), the integrand in equation (27) is convex on U.
Furthermore, for proving the condition in (vi), it can be seen that, W 3 u 2 3 W 3 as W 3 2 [0, 1], that is, This completes the proof. ,

Characterization of optimal control solutions
The Pontryagin's Maximum Principle (PMP) introduced in (Fleming & Rishel, 2012) provides the necessary optimality condition as stated in the following theorem.
From the boundedness of u * i ðtÞ on (0, 1) and the minimality condition, we have: ð1ÀkrÞb c ðI c þtðI tc þI cL ÞÞ½ðx 3 Àx 1 ÞSþðx 3 Àx 2 Þð1ÀεÞV þðx 3 Àx 7 ÞR t þðx 8 Àx 5 ÞLþðx 9 Àx 6 ÞsI t W 2 N ; if vH vu 1 ¼ 0; Concluding that the optimal control functions u * i ðtÞ can be written in compact form as in (32). , To verify that the optimal control is indeed a minimizer, we need to compute the Hessian matrix of the Hamiltonian H function and check its positive definiteness. Indeed, the Hessian matrix for H w.r.t the control u is given by Clearly, this matrix is a positive definite due to the positive weights W 1 ,W 2 and W 3 . This confirms that the corresponding Hamiltonian function H is convex, and as a result, the optimal control, u ¼ (u 1 , u 2 , u 3 ), is a minimizer (Abraha et al., 2021).
The proposed control induced system (26) is an initial value problem while the adjoint system (31) subject to the transversality conditions x i (t f ) ¼ 0, i ¼ 1, …, 9 is a terminal value problem. Therefore, the boundary value problem given by Theorem 13 need to be solved numerically using forward iteration in the state system, and the adjoint system (31) is integrated through backward iteration (Campos, Silva, & Torres, 2019;Kumar et al., 2020;Lenhart & Workman, 2007;Silva & Torres, 2013).

Simulations of optimal control model
In this section, several numerical simulations of model (26) is performed to support the analytical results. To obtain the optimal solution to the control induced model (26), we apply the forward-backward sweep method presented in (Lenhart & Workman, 2007). This iterative method consists in solving the system of equations in (26) with the guess for the controls over the time interval [0, t f ] using a forward fourth-order Runge-Kutta scheme and the transversality conditions. Then, the adjoint system is solved by applying a backward fourth-order Runge-Kutta scheme to the current iteration solution of the system of equation (26) (For details see references (Campos et al., 2019;Lenhart & Workman, 2007)). The controls are updated by means of a convex combination of the previous controls and the values in the characterization process (32). To repeat the solutions of state variables in (26) and adjoint systems, the updated controls are used. This situation is continued until the values of unknowns at the previous iteration and at the present iteration are close enough.
To assess the cost of each applied strategy, the following weight constants are presumptions.
Also, the cost of TB treatment is assumed to be higher than the cost of implementing COVID-19 prevention control (facemask usage and vaccination).

Initial conditions
In this subsection, we define the initial conditions that correspond to the number of individuals in each class in order to solve model (1). In this regard, the initial total population is assumed to be N(0) ¼ 114, 963, 588 (Kifle & Obsu, 2022). Further, we assume the initial conditions for the state variables used in model simulation as follows: the total number of fully vaccinated individuals against COVID-19 was V(0) ¼ 400, 000, L(0) ¼ 100, 000, I t (0) ¼ 30, 000, I c (0) ¼ 280, 565, I cL (0) ¼ 50, 000, I tc (0) ¼ 5, 000, R c (0) ¼ 263, 587 and R t (0) ¼ 7000. Thus, the initial susceptible population is determined as We start the numerical simulations of optimal control model (26) for the infected populations when no control strategies are applied as shown in Fig. 2. For the numerical simulations, we used the parameter values given in Table 2 which results the basic reproduction number R 0 ¼ maxfR vC ; R 0T g ¼ maxf0:82; 0:93g ¼ 0:93 < 1. As we can observed from Fig. 2, the solutions are not converging to zero. This is because of the existence of backward bifurcation. That is, for R 0 < 1, we can not eliminate the disease from the community. Fig. 3 shows the forward and backward bifurcation diagrams. Fig. 3a shows that the system (3) undergoes the transcritical bifurcation at R vC and changes the stability from the disease free equilibrium (DFE) to the endemic equilibrium (EE). Precisely, the DFE is stable when R vC < 1 and DFE becomes an unstable equilibrium and EE becomes a stable equilibrium for the model (3) when R vC > 1. Further, Fig. 3b shows that the system (14) undergoes the backward bifurcation at R 0T where a stable endemic equilibrium co-exists with a stable disease-free equilibrium for R 0T < 1. Clearly, this phenomenon has significant public health consequences, as it makes the classical requirement of the associated basic reproduction number being less than unity to be necessary, but not sufficient to eradicate the disease.
(a) Forward bifurcation for COVID-19 sub-model (b) Backward bifurcation for TB sub-model The impact of implementing various possible control strategies on the codynamics of COVID-19 and tuberculosis disease is investigated as follows.
5.1.1. Strategy A: prevention against COVID-19 (u 1 s 0) The simulations of the optimal control system (26) when the control strategy COVID-19 prevention (u 1 ) is implemented, are presented in Figs. 4e6. As shown in Fig. 4a, the control u 1 greatly reduce the COVID-19 infected from the community within the intervention time. Further, from Fig. 5a and 6a, it is demonstrated that the COVID-19 prevention strategy have a significant effect on reducing the co-infected individuals. From Fig. 4b and 5b, one can also observe that the COVID-19 prevention method like face-mask usage has an effect in reducing individuals become TB infected. The control profile u 1 representing the effect of prevention against COVID-19 control strategy is depicted in Fig. 6b. When prevention against COVID-19 is the selected strategy, it is observed that the prevention against COVID-19 (u 1 ) should be implemented optimally throughout the intervention period, and then it will start decreasing after about 280 days.

Strategy B: vaccination against COVID-19 (u 2 s 0)
In this case, the COVID-19 vaccination control strategy (u 2 ) is implemented to minimize the disease burden from the population as presented in Fig. 7 and 8. From Fig. 7a, one can observe that when the control u 2 is considered, the optimal strategy provide a reduction of 3 Â 10 5 COVID-19 infected individuals. The corresponding control profile u 2 is also depicted in Fig. 8b. As shown in Fig. 8b, to minimize the cost functional, the optimal control effort u 2 decreases from the maximum to the minimum in the final time. That is, COVID-19 vaccination u 2 should be implemented optimally throughout the intervention period, and then it will start decreasing after about 320 days.

Strategy C: TB treatment (u 3 s 0)
In this case, the TB treatment control strategy u 3 is implemented to observe its effects in minimizing the number of TB infected and co-infected as illustrated in Fig. 9 and 10. From Fig. 9a, one can observe that when the control u 3 is considered, the optimal strategy provide a reduction of 1 Â 10 5 latent TB infected individuals. Similarly, from Fig. 9b, it can be seen that the TB treatment control u 3 have a significant effect on reducing actively infected TB. Further, the effects of TB treatment on reducing the coinfected persons are also depicted in Fig. 10a and b. The corresponding control profile u 3 is also portrayed in Fig. 11, and it can be seen that the TB treatment control u 3 is optimally implemented at the onset of the epidemic, but drastically decreases and only picks up again from day 330 and become effective throughout the simulation period. 5.1.4. Strategy D: COVID-19 prevention and vaccination (u 1 , u 2 s 0) In Fig. 12 and 13, the implementation of control strategy that combines COVID-19 prevention and vaccination (u 1 s 0, u 2 s 0) is implemented. From Fig. 12a as expected, with implementation of this strategy a significant number of COVID-19 infected persons are reduced. Besides, it has a positive impact on reducing the co-infected individuals as shown in Fig. 12b and 13a. The corresponding combined effect of controls u 1 and u 2 are also graphically depicted in Fig. 13b. When the combination of these control strategy are selected, one notes that COVID-19 prevention (control u 1 ) should be implemented optimally and will start decreasing after about one year throughout the intervention period, while COVID-19 vaccination u 2 uptake will start decreasing after about 20 days. 5.1.5. Strategy E: COVID-19 prevention and TB treatment (u 1 , u 3 s 0) In Fig. 14 and 15, the implementation of the control strategy that combines COVID-19 prevention and TB treatment (u 1 s 0, u 3 s 0) is implemented. From Fig. 14a, with implementation of this strategy a significant number of COVID-19 and latent TB co-infected persons are reduced. Besides, it has a positive impact on reducing active TB infected individuals as shown in Fig. 14b. Furthermore, from Fig. 15a one can observe that, when this control strategies are considered, the optimal strategy provides a reduction of 2 Â 10 5 latent TB infected cases. The corresponding combined effect of controls u 1 and u 3 are also graphically depicted in Fig. 15b. When the combination of these control strategy are selected, one notes that COVID-19 prevention (control u 1 ) should be implemented optimally and will start decreasing after about 350 days throughout the intervention period, while TB treatment u 3 uptake will start decreasing after about 60 days.

Conclusions
In this paper, we presented a new mathematical model for the co-dynamics of COVID-19 and TB which incorporates the COVID-19 vaccination and exogenous reinfection for TB. The model's basic properties, including the invariant region, positivity, and boundedness of solutions, are given. Furthermore, theoretical findings such as stability and bifurcation analysis of each sub-models equilibria are discussed. It is shown that the DFE is locally asymptotically stable for both COVID-19 and TB sub-models when the associated basic reproduction numbers R vC < 1 and R 0T < 1, and unstable otherwise. The existence and stability of the endemic equilibrium point was determined for both sub models. Besides, with the help of centre manifold theory, we also proved that the COVID-19 sub-model exhibit forward bifurcation, while the TB sub-model exhibit the backward bifurcation. That is, the incorporation of exogenous reinfection into our model causes effects by allowing the occurrence of backward bifurcation for R 0 < 1 and the exogenous reinfection rate larger than a threshold, (i.e., h > h*). Based on the occurrence of backward bifurcation for TB sub-model, and since the dynamics of COVID-19 and TB co-infection model is driven by that of its sub-model, then the model (1) exhibit the phenomenon of backward bifurcation. Thus, the possible implication of this phenomenon is that reducing R 0 < 1 may not be sufficient to eliminate the disease burden from the community. Further, we observed that when h ¼ 0, system (1) reduces to the coinfection model with out exogenous reinfection. In this case, the backward bifurcation phenomenon will not occur and the disease elimination from the community is possible for R 0 < 1.
The proposed co-infection model is then extended to included three control measures. The conditions for the existence of optimal control and the optimality system for the control induced model are established using Pontryagin's maximum principle. Furthermore, different numerical simulations of the control induced model are also performed to observe the effects of the control strategies by using model parameters given in Table 2. To further understand the dynamical interactions between COVID-19 and TB with and without optimal control, simulations are carried out in two scenarios, one with single control and other the combination of two controls. Results of these two scenarios are compared with the model without control and it reveals that the infected sub-population reduces significantly when combinations of the controls are applied. When COVID-19 vaccination is considered as control strategy, the optimal strategy provide a reduction of significant number of COVID-19 infected cases. Further, the combination of COVID-19 prevention and vaccination control strategy greatly mitigates new co-infection cases.

Data availability
All datasets generated in this study can be recovered from information available within the article.

Declaration of competing interest
The authors declare that there is no conflict of interest. Thus, the next generation matrix FV À1 is given by Taking the dominant eigenvalue of FV À1 , the basic reproductive number for TB sub-model (14) is given by :