*Article* **Kinetics of a Reaction-Diffusion Mtb/SARS-CoV-2 Coinfection Model with Immunity**

**Ali Algarni 1,\*, Afnan D. Al Agha 2, Aisha Fayomi <sup>1</sup> and Hakim Al Garalleh <sup>2</sup>**


**Abstract:** The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and Mycobacterium tuberculosis (Mtb) coinfection has been observed in a number of nations and it is connected with severe illness and death. The paper studies a reaction–diffusion within-host Mtb/SARS-CoV-2 coinfection model with immunity. This model explores the connections between uninfected epithelial cells, latently Mtb-infected epithelial cells, productively Mtb-infected epithelial cells, SARS-CoV-2-infected epithelial cells, free Mtb particles, free SARS-CoV-2 virions, and CTLs. The basic properties of the model's solutions are verified. All equilibrium points with the essential conditions for their existence are calculated. The global stability of these equilibria is established by adopting compatible Lyapunov functionals. The theoretical outcomes are enhanced by implementing numerical simulations. It is found that the equilibrium points mirror the single infection and coinfection states of SARS-CoV-2 with Mtb. The threshold conditions that determine the movement from the monoinfection to the coinfection state need to be tested when developing new treatments for coinfected patients. The impact of the diffusion coefficients should be monitored at the beginning of coinfection as it affects the initial distribution of particles in space.

**Keywords:** tuberculosis; COVID-19; diffusion; coinfection; stability

**MSC:** 35B35; 37N25; 92B05

#### **1. Introduction**

Coronavirus disease 2019 (COVID-19) is a viral disease induced by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and emerged in 2019. Although the number of new cases decreased in the last few months, COVID-19 is continuing its spread around the globe [1]. Following the World Health Organization (WHO) report issued on 1 March 2023, above 758,000,000 affirmed cases and over 6,800,000 deaths have been accounted globally [1]. COVID-19 coinfections with other viral or bacterial diseases are common, which complicates the treatment of COVID-19 [2]. Tuberculosis (TB) is a bacterial infection attributable to Mycobacterium tuberculosis (Mtb). Currently, COVID-19 cooccurring with TB has been declared in a number of nations [3]. As COVID-19 and TB are greatly infectious diseases, understanding Mtb/SARS-CoV-2 coinfection is very crucial for protection and treatment of coinfection.

SARS-CoV-2 is an enveloped RNA virus which is linked with the Betacoronavirus genus [4]. It breaks into the host cell using the angiotensin-converting enzyme 2 (ACE2) receptor [5]. It primarily infects the alveolar epithelial type-II cells of the lungs [6]. Similar to SARS-CoV-2, Mtb infects alveolar epithelial type-II cells through pattern recognition receptors such as toll-like receptors, complement receptors, and CD14 receptors [7]. Thus, the lung is the major infection site for these pathogens. Nevertheless, they can invade cells within different organs [6]. Since SARS-CoV-2 and Mtb infect the same target, this could

**Citation:** Algarni, A.; Al Agha, A.D.; Fayomi, A.; Al Garalleh, H. Kinetics of a Reaction-Diffusion Mtb/SARS-CoV-2 Coinfection Model with Immunity. *Mathematics* **2023**, *11*, 1715. https://doi.org/10.3390/ math11071715

Academic Editor: Patricia J. Y. Wong

Received: 9 March 2023 Revised: 26 March 2023 Accepted: 31 March 2023 Published: 3 April 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

increase the seriousness of disease in coinfected people [4]. Both pathogens are disseminated through respiratory droplets [2]. The most dominant features of Mtb/SARS-CoV-2 co-occurring are fever, cough, and dyspenea [2]. Risk factors in coinfection include age and comorbidities such as diabetes, HIV, and hypertension [2,3]. The immune response in coinfection includes T cells [8]. Specifically, cytotoxic T lymphocytes (CTLs) work on eliminating infected cells from the body. In the mild cases, the immune response can clear both infections. It has been proposed that Mtb/SARS-CoV-2 patients are at higher risk of death and developing severe disease than SARS-CoV-2 patients without Mtb [2,3,8]. Moreover, some studies reported that SARS-CoV-2 infection may cause latent Mtb to become active in coinfected people [3,4]. Other studies also observed that coinfected patients have low lymphocyte counts [2,8]. Thus, understanding the mechanism of coinfection is very important to evolve treatments for coinfected patients.

Mathematical models have been utilized to assist experimental and medical studies of different infections. These models are partitioned into epidemiological and within-host models. Epidemiological systems consider the interactions between individuals at the population level, while within-host systems explore the interplay between pathogens and cells within the host's body. A variety of COVID-19 epidemiological models (see for example, [9–15]) and within-host models (see for example, [16–18]) have been introduced and investigated. Similarly, TB models have been widely studied as epidemiological models (see for example, [19–22]) and within-host models (see for example, [23–27]). Some coinfection models of COVID-19 with other diseases have been developed. For instance, Pinky and Dobrovolny [28] analyzed a model that tests the impact of SARS-CoV-2 coinfection with the influenza virus. Al Agha and Elaiw [29] established a within-host SARS-CoV-2/malaria model with immune response. Elaiw et al. [30] developed a SARS-CoV-2/HIV coinfection model that takes the latent stage of infected epithelial cells (EPCs) into consideration. Elaiw and Al Agha [31] studied a SARS-CoV-2/cancer system with two immune responses.

Many epidemiological models of TB/COVID-19 coinfection have been proposed (see for example, [32–34]). On the other hand, within-host models are not widely considered. In [35], a within-host coinfection model has been formalized using ordinary differential Equations (ODEs). This work develops a reaction–diffusion within-host Mtb/SARS-CoV-2 coinfection model. It depicts the interplay between uninfected EPCs, latently Mtb-infected EPCs, productively Mtb-infected EPCs, SARS-CoV-2-infected EPCs, Mtb particles, SARS-CoV-2 particles, and CTLs. Additionally, this model is formalized using partial differential Equations (PDEs) which count the nonuniform distribution of cells and pathogens with their ability to move. Thus, PDEs are more realistic than ODEs which assume the spatial distribution homogeneity of cells and particles. Using the developed model, we (i) establish the boundedness and nonnegativity of the solutions, (ii) determine all equilibrium points and find the thresholds, (iii) confirm the global stability of each point, and (iv) use numerical simulations to validate the theoretical observations.

The remaining sections are divided as follows. Section 2 represents the model. Section 3 proves the boundedness and nonnegativity of the solutions. Moreover, it recounts all equilibrium points. Section 4 employs Lyapunov functionals to show the global stability of each point. Section 5 implements numerical simulations. The last section provides the conclusion with upcoming works.

#### **2. A Reaction–Diffusion Mtb/SARS-CoV-2 Coinfection Model**

This part describes the model under consideration. In this model, we assume that Mtb and SARS-CoV-2 have the same target, and CTLs kill infected cells at the same rate. The model consists of seven PDEs as follows:

$$\begin{cases} \frac{\partial II(x,t)}{\partial t} = D\_{l\ell} \Delta I + \lambda - \eta\_1 IB - \eta\_2 LV - \varepsilon\_1 I\_{\prime} \\ \frac{\partial L(x,t)}{\partial t} = D\_L \Delta L + \eta\_1 IB - aL\_{\prime} \\ \frac{\partial I^B(x,t)}{\partial t} = D\_{I^B} \Delta I^B + aL - \gamma I^B Z - \varepsilon\_2 I^B, \\ \frac{\partial I^V(x,t)}{\partial t} = D\_{I^V} \Delta I^V + \eta\_2 LV - \gamma I^V Z - \varepsilon\_3 I^V, \\ \frac{\partial B(x,t)}{\partial t} = D\_B \Delta B + \mu\_1 \varepsilon\_2 I^B - \varepsilon\_4 B, \\ \frac{\partial V(x,t)}{\partial t} = D\_V \Delta V + \mu\_2 I^V - \varepsilon\_5 V, \\ \frac{\partial Z(x,t)}{\partial t} = D\_Z \Delta Z + \omega I^B Z + \omega I^V Z - \varepsilon\_6 Z, \end{cases} (1)$$

where the time *<sup>t</sup>* > 0 and the position *<sup>x</sup>* <sup>∈</sup> <sup>Ψ</sup>. The domain <sup>Ψ</sup> is bounded and connected with a smooth boundary *∂*Ψ. The compartments *U*, *L*, *IB*, *IV*, *B*, *V*, and *Z* designate the concentrations of uninfected EPCs, latently Mtb-infected EPCs, productively Mtb-infected EPCs, SARS-CoV-2-infected EPCs, Mtb particles, SARS-CoV-2 particles, and CTLs at (*x*, *t*), respectively. Uninfected EPCs are generated at rate *λ*. Mtb converts healthy EPCs into latently infected cells at rate *η*1*UB*, while SARS-CoV-2 infects the same type of cells at rate *η*2*UV*. Latently Mtb-infected cells become an active producer at rate *aL*. Mtb particles are created at a total production rate *μ*1<sup>2</sup> *IB*. SARS-CoV-2 virions are ejected from SARS-CoV-2-infected cells at rate *μ*<sup>2</sup> *IV*. CTLs remove Mtb and SARS-CoV-2 infected cells at rates *γIBZ* and *γI<sup>V</sup> Z*, respectively. The corresponding stimulation rates are *ωIBZ* and *ωI<sup>V</sup> Z*, respectively. The compartments *U*, *IB*, *IV*, *B*, *V*, and *Z* die at rates 1*U*, <sup>2</sup> *IB*, <sup>3</sup> *IV*, 4*B*, 5*V*, and 6*Z*, respectively. We assume that each compartment *K* diffuses with a diffusion coefficient *DK*. The operator <sup>Δ</sup> <sup>=</sup> *<sup>∂</sup>*<sup>2</sup> *<sup>∂</sup>x*<sup>2</sup> is the Laplacian operator. We presume that all parameters of model (1) are positive. The initial conditions (ICs) of system (1) are

$$\begin{aligned} L(\mathbf{x},0) &= \nu\_1(\mathbf{x}), \quad L(\mathbf{x},0) = \nu\_2(\mathbf{x}), \quad I^B(\mathbf{x},0) = \nu\_3(\mathbf{x}), \quad I^V(\mathbf{x},0) = \nu\_4(\mathbf{x}),\\ B(\mathbf{x},0) &= \nu\_5(\mathbf{x}), \quad V(\mathbf{x},0) = \nu\_6(\mathbf{x}), \quad Z(\mathbf{x},0) = \nu\_7(\mathbf{x}), \quad \mathbf{x} \in \Psi,\end{aligned} \tag{2}$$

where *<sup>ν</sup>i*(*x*) <sup>≥</sup> 0, *<sup>i</sup>* <sup>=</sup> 1, 2, ... , 7, are continuous functions in <sup>Ψ</sup>¯ . The boundary conditions (BCs) of (1) are

$$\frac{\partial I}{\partial \vec{r}} = \frac{\partial L}{\partial \vec{r}} = \frac{\partial I\_B}{\partial \vec{r}} = \frac{\partial I\_V}{\partial \vec{r}} = \frac{\partial B}{\partial \vec{r}} = \frac{\partial V}{\partial \vec{r}} = \frac{\partial Z}{\partial \vec{r}} = 0, \quad t > 0, \quad x \in \partial \Psi,\tag{3}$$

where *<sup>∂</sup> <sup>∂</sup><sup>r</sup>* is the outward normal derivative on *<sup>∂</sup>*Ψ. These Neumann BCs suggest that the boundary is isolated.

In the upcoming parts of the paper and for simplicity, we consider the contraction *K*(*x*, *t*) ≡ *K* for each compartment *K* in model (1).

#### **3. Basic Properties**

This section certifies that the solutions of system (1)–(3) are unique, nonnegative, and bounded. Additionally, it computes all equilibrium points of model (1).

Let S = *BUC* Ψ¯ , R<sup>7</sup> be the set of functions that are bounded and continuous from Ψ¯ to R7. The positive cone S<sup>+</sup> = *BUC* Ψ¯ , R<sup>7</sup> + <sup>⊂</sup> <sup>S</sup> forms a partial order on <sup>S</sup>. Define *f* <sup>S</sup> = sup *<sup>x</sup>*∈Ψ¯ <sup>|</sup> *<sup>f</sup>*(*x*)|. Consequently, (S, ·S) is a Banach lattice [36,37].

**Theorem 1.** *Suppose that DU* = *DL* = *DIB* = *DIV* = *DZ* = *D*1*. Then, model* (1) *with any ICs* (2) *has a unique, nonnegative, and bounded solution on* <sup>Ψ</sup>¯ <sup>×</sup> [0, <sup>+</sup>∞)*.*

**Proof.** For any *<sup>ν</sup>* = (*ν*1, *<sup>ν</sup>*2, *<sup>ν</sup>*3, *<sup>ν</sup>*4, *<sup>ν</sup>*5, *<sup>ν</sup>*6, *<sup>ν</sup>*7)*<sup>T</sup>* <sup>∈</sup> <sup>S</sup>+, we define *<sup>P</sup>* = (*P*1, *<sup>P</sup>*2, *<sup>P</sup>*3, *<sup>P</sup>*4, *<sup>P</sup>*5, *<sup>P</sup>*6, *<sup>P</sup>*7)*<sup>T</sup>* : <sup>S</sup><sup>+</sup> <sup>→</sup> <sup>S</sup> by

$$\begin{cases} P\_1(\nu)(\mathbf{x}) = \lambda - \eta\_1 \nu\_1(\mathbf{x})\nu\_5(\mathbf{x}) - \eta\_2 \nu\_1(\mathbf{x})\nu\_6(\mathbf{x}) - \epsilon\_1 \nu\_1(\mathbf{x}), \\ P\_2(\nu)(\mathbf{x}) = \eta\_1 \nu\_1(\mathbf{x})\nu\_5(\mathbf{x}) - \mu \nu\_2(\mathbf{x}), \\ P\_3(\nu)(\mathbf{x}) = a\nu\_2(\mathbf{x}) - \gamma \nu\_3(\mathbf{x})\nu\_7(\mathbf{x}) - \epsilon\_2 \nu\_3(\mathbf{x}), \\ P\_4(\nu)(\mathbf{x}) = \eta\_2 \nu\_1(\mathbf{x})\nu\_6(\mathbf{x}) - \gamma \nu\_4(\mathbf{x})\nu\_7(\mathbf{x}) - \epsilon\_3 \nu\_4(\mathbf{x}), \\ P\_5(\nu)(\mathbf{x}) = \mu\_1 \epsilon\_2 \nu\_3(\mathbf{x}) - \epsilon\_4 \nu\_5(\mathbf{x}), \\ P\_6(\nu)(\mathbf{x}) = \mu\_2 \nu\_4(\mathbf{x}) - \epsilon\_5 \nu\_6(\mathbf{x}), \\ P\_7(\nu)(\mathbf{x}) = \omega \nu\_3(\mathbf{x})\nu\_7(\mathbf{x}) + \omega \nu\_4(\mathbf{x})\nu\_7(\mathbf{x}) - \epsilon\_6 \nu\_7(\mathbf{x}). \end{cases}$$

We observe that *P* is Lipschitz on S+. Therefore, it is possible to rewrite system (1)–(3) as the abstract DE:

$$\begin{cases} \frac{df}{dt} = Df + P(f)\_{\prime\prime}t > 0\_{\prime\prime} \\ f\_0 = \nu \in \mathbb{S}\_{+\prime} \end{cases}$$

where *J* = (*U*, *L*, *IB*, *IV*, *B*, *V*, *Z*)*<sup>T</sup>* and *D J* = (*DU*Δ*U*, *DL*Δ*L*, *DIB*Δ*IB*, *DIV* Δ*IV*, *DB*Δ*B*, *DV*Δ*V*, *DZ*Δ*Z*)*T*. We can show that

$$\lim\_{h \to 0^+} \frac{1}{h} \operatorname{dist}(\nu + hP(\nu), \mathbb{S}\_+) = 0, \,\,\nu \in \mathbb{S}\_+.$$

Hence, for any *<sup>ν</sup>* <sup>∈</sup> <sup>S</sup>+, model (1)–(3) has a unique nonnegative mild solution for the time interval [0, *Te*).

To verify the boundedness, we consider the function

$$\mathcal{Y}\_1(\mathbf{x}, \mathbf{t}) = \mathcal{U} + L + I^V.$$

*As DU* = *DL* = *DIV* = *D*1, then by utilizing model (1) we obtain

$$\begin{aligned} \frac{\partial Y\_1}{\partial t} - D\_1 \Delta Y\_1 &= \lambda - \epsilon\_1 \mathcal{U} - aL - \gamma I^V Z - \epsilon\_3 I^V \\ &\le \lambda - \epsilon\_1 \mathcal{U} - aL - \epsilon\_3 I^V \\ &\le \lambda - \sigma\_1 \left[ \mathcal{U} + L + I^V \right] \\ &= \lambda - \sigma\_1 \mathcal{Y}\_{1'} \end{aligned}$$

where *σ*<sup>1</sup> = min{*a*, 1, 3}. Thus, Υ<sup>1</sup> satisfies the system

$$\begin{cases} \frac{\partial \mathbf{Y}\_1}{\partial t} - D\_1 \Delta \mathbf{Y}\_1 \le \lambda - \sigma\_1 \mathbf{Y}\_{1\prime}, \\\frac{\partial \mathbf{Y}\_1}{\partial \vec{r}} = \mathbf{0}, \\\mathbf{Y}\_1(x, 0) \ge 0. \end{cases}$$

Assume that Υ01(*t*) satisfies the system

$$\begin{cases} \frac{d\tilde{\mathbf{Y}}\_1(t)}{dt} = \lambda - \sigma\_1 \tilde{\mathbf{Y}}\_1(t),\\ \tilde{\mathbf{Y}}\_1(0) = \max\_{\mathbf{x} \in \mathbb{F}} \mathbf{Y}\_1(\mathbf{x}, 0), \end{cases}$$

which implies that <sup>Υ</sup>01(*t*) <sup>≤</sup> max *<sup>λ</sup> σ*1 , max *<sup>x</sup>*∈Ψ¯ <sup>Υ</sup>1(*x*, 0) 6 . In accord with the comparison principle (CP) [38], we obtain Υ1(*x*, *t*) ≤ Υ01(*t*). As a result, we have

$$\mathcal{Y}\_1(\mathfrak{x}, t) \le \max \left\{ \frac{\lambda}{\mathcal{O}\_1}, \max\_{\mathbf{x} \in \Psi} \mathcal{Y}\_1(\mathfrak{x}, 0) \right\} := Q\_1. \mathfrak{x}$$

This ensures that *U*, *L*, and *I<sup>V</sup>* are bounded. From the third equation of (1), we have

$$\begin{aligned} \frac{\partial I^B}{\partial t} - D\_{I^B} \Delta I^B &= aL - \gamma I^B Z - \epsilon\_2 I^B \\ &\le aL - \epsilon\_2 I^B \\ &\le aQ\_1 - \epsilon\_2 I^B. \end{aligned}$$

We can conclude from the CP [38] that

$$I^B \le \max\left\{ \frac{aQ\_1}{\epsilon\_2}, \max\_{\mathbf{x} \in \mathbb{P}} I^B(\mathbf{x}, 0) \right\} := Q\_2.$$

Thus, *I<sup>B</sup>* is bounded. From the fifth equation of (1), we have

$$\begin{aligned} \frac{\partial B}{\partial t} - D\_B \Delta B &= \mu\_1 \varepsilon\_2 I^B - \varepsilon\_4 B \\ &\leq \mu\_1 \varepsilon\_2 Q\_2 - \varepsilon\_4 B. \end{aligned}$$

Based on the CP [38], we obtain

$$B \le \max\left\{\frac{\mu\_1 \varepsilon\_2 Q\_2}{\epsilon\_4}, \max\_{x \in \Psi} B(x, 0)\right\}.$$

Hence, *B* is bounded. From the sixth equation of model (1), we obtain

$$\begin{aligned} \frac{\partial V}{\partial t} - D\_V \Delta V &= \mu\_2 I^V - \epsilon\_5 V \\ &\leq \mu\_2 Q\_1 - \epsilon\_5 V. \end{aligned}$$

The CP [38] implies that

$$V \le \max\left\{\frac{\mu\_2 Q\_1}{\epsilon\_5}, \max\_{\mathbf{x} \in \Psi} V(\mathbf{x}, 0)\right\} := Q\_3.$$

Thus, *V* is bounded. Finally, we introduce the function

$$\Upsilon\_2(\mathfrak{x}, t) = I^B + I^V + \frac{\gamma}{\omega} Z.$$

Then, we obtain

$$\begin{aligned} \frac{\partial Y\_2}{\partial t} - D\_1 \Delta Y\_2 &= aL + \eta\_2 LU - \varepsilon\_2 I^B - \varepsilon\_3 I^V - \frac{\gamma \varepsilon\_6}{\omega} Z\_2\\ &\leq aQ\_1 + \eta\_2 Q\_1 Q\_3 - \sigma\_2 Y\_{2\prime} \end{aligned}$$

where *σ*<sup>2</sup> = min{2, 3, 6}. By the CP, [38], we obtain

$$Y\_2(x,t) \le \max\left\{\frac{aQ\_1 + \eta\_2 Q\_1 Q\_3}{\sigma\_2}, \max\_{x \in \Psi} Y\_2(x,0)\right\}.$$

This implies that *Z* is bounded. The above results show that all solutions are bounded on <sup>Ψ</sup>¯ <sup>×</sup> [0, *Te*), and so solutions are bounded on <sup>Ψ</sup>¯ <sup>×</sup> [0, <sup>+</sup>∞). This conclusion is derived from the standard theory for semi-linear parabolic Equations [39].

**Proposition 1.** *The conditions* R0*B,* R0*V,* R1*B,* R1*V, and σ exist such that system* (1) *has six equilibrium points:*


$$\frac{\varepsilon\_3}{\varepsilon\_2} + \mathcal{R}\_{1b}, \frac{\omega \lambda \eta\_1 \mu\_1 \varepsilon\_2 \varepsilon\_5}{\varepsilon\_3 \varepsilon\_4 [\omega \varepsilon\_1 \varepsilon\_5 + \eta\_2 \mu\_2 \varepsilon\_6]} + 1 > \frac{\varepsilon\_2}{\varepsilon\_3} + \mathcal{R}\_{1V}, \frac{\mathcal{R}\_{0V}}{\mathcal{R}\_{0B}} > 1, \epsilon\_2 > \epsilon\_3, \text{and } \sigma > 1.$$

**Proof.** The equilibrium points of Equation (1) can be drawn by solving the system:

$$\begin{cases} 0 = \lambda - \eta\_1 LU - \eta\_2 UI - \varepsilon\_1 II, \\ 0 = \eta\_1 UIB - aL, \\ 0 = aL - \gamma I^B Z - \varepsilon\_2 I^B, \\ 0 = \eta\_2 UIV - \gamma I^V Z - \varepsilon\_3 I^V, \\ 0 = \mu\_1 \varepsilon\_2 I^B - \varepsilon\_4 B, \\ 0 = \mu\_2 I^V - \varepsilon\_5 V, \\ 0 = \omega I^B Z + \omega I^V Z - \varepsilon\_6 Z. \end{cases}$$

Then, we obtain the following:


$$\mathcal{U}\_{1} = \frac{\mathfrak{e}\_{4}}{\eta\_{1}\mu\_{1}}, \quad L\_{1} = \frac{\mathfrak{e}\_{1}\mathfrak{e}\_{4}}{a\eta\_{1}\mu\_{1}} \left(\mathcal{R}\_{0\mathcal{B}} - 1\right), \quad I\_{1}^{\mathcal{B}} = \frac{\mathfrak{e}\_{1}\mathfrak{e}\_{4}}{\eta\_{1}\mu\_{1}\mathfrak{e}\_{2}} \left(\mathcal{R}\_{0\mathcal{B}} - 1\right), \quad B\_{1} = \frac{\mathfrak{e}\_{1}}{\eta\_{1}} \left(\mathcal{R}\_{0\mathcal{B}} - 1\right),$$

where <sup>R</sup>0*<sup>B</sup>* <sup>=</sup> *λη*1*μ*<sup>1</sup> 1<sup>4</sup> . We note that *U*<sup>1</sup> is positive, whilst *L*1, *I<sup>B</sup>* <sup>1</sup> , and *B*<sup>1</sup> are positive when <sup>R</sup>0*<sup>B</sup>* <sup>&</sup>gt; 1. Hence, *<sup>E</sup>*<sup>1</sup> exists if <sup>R</sup>0*<sup>B</sup>* <sup>&</sup>gt; 1. The parameter <sup>R</sup>0*<sup>B</sup>* appoints the onset of Mtb infection with inactive CTLs.

(iii) The COVID-19 immune-free equilibrium *E*<sup>2</sup> = *U*2, 0, 0, *I<sup>V</sup>* <sup>2</sup> , 0, *<sup>V</sup>*2, 0 . Its coordinates are written as follows:

$$\mathcal{U}\_2 = \frac{\mathfrak{e}\_3 \mathfrak{e}\_5}{\eta\_2 \mu\_2}, \quad I\_2^V = \frac{\mathfrak{e}\_1 \mathfrak{e}\_5}{\eta\_2 \mu\_2} \left(\mathcal{R}\_{0V} - 1\right), \quad V\_2 = \frac{\mathfrak{e}\_1}{\eta\_2} \left(\mathcal{R}\_{0V} - 1\right),$$

where <sup>R</sup>0*<sup>V</sup>* <sup>=</sup> *λη*2*μ*<sup>2</sup> 13<sup>5</sup> . Thus, *<sup>U</sup>*<sup>2</sup> <sup>&</sup>gt; 0, whilst *<sup>I</sup><sup>V</sup>* <sup>2</sup> <sup>&</sup>gt; 0 and *<sup>V</sup>*<sup>2</sup> <sup>&</sup>gt; 0 when <sup>R</sup>0*<sup>V</sup>* <sup>&</sup>gt; 1. Accordingly, *<sup>E</sup>*<sup>2</sup> is defined when <sup>R</sup>0*<sup>V</sup>* <sup>&</sup>gt; 1. The threshold <sup>R</sup>0*<sup>V</sup>* locates the start of COVID-19 infection, where the CTL immunity is inactive.

(iv) The Mtb equilibrium with immunity *E*<sup>3</sup> = *U*3, *L*3, *I<sup>B</sup>* <sup>3</sup> , 0, *B*3, 0, *Z*<sup>3</sup> , where

$$\begin{array}{ll} \mathcal{U}\_{3} = \frac{\omega \lambda \varepsilon\_{4}}{\omega \varepsilon\_{1} \varepsilon\_{4} + \eta\_{1} \mu\_{1} \varepsilon\_{2} \varepsilon\_{6}}, & \mathcal{L}\_{3} = \frac{\lambda \eta\_{1} \mu\_{1} \varepsilon\_{2} \varepsilon\_{6}}{a \left[\omega \varepsilon\_{1} \varepsilon\_{4} + \eta\_{1} \mu\_{1} \varepsilon\_{2} \varepsilon\_{6}\right]}, & \mathcal{I}\_{3}^{B} = \frac{\varepsilon\_{6}}{\omega}, \\\mathcal{B}\_{3} = \frac{\mu\_{1} \varepsilon\_{2} \varepsilon\_{6}}{\omega \varepsilon\_{4}}, & \mathcal{Z}\_{3} = \frac{\varepsilon\_{2}}{\gamma} \left(\mathcal{R}\_{1B} - 1\right), \end{array}$$

where <sup>R</sup>1*<sup>B</sup>* <sup>=</sup> *ωλη*1*μ*<sup>1</sup> *ω*1<sup>4</sup> + *η*1*μ*12<sup>6</sup> . We note that *U*3, *L*3, *I<sup>B</sup>* <sup>3</sup> and *B*<sup>3</sup> are always positive, while *<sup>Z</sup>*<sup>3</sup> <sup>&</sup>gt; 0 if <sup>R</sup>1*<sup>B</sup>* <sup>&</sup>gt; 1. This implies that *<sup>E</sup>*<sup>3</sup> exists if <sup>R</sup>1*<sup>B</sup>* <sup>&</sup>gt; 1. The threshold <sup>R</sup>1*<sup>B</sup>* sets the activation of CTLs versus Mtb-infected EPCs.

(v) The COVID-19 equilibrium with immunity *E*<sup>4</sup> = *U*4, 0, 0, *I<sup>V</sup>* <sup>4</sup> , 0, *V*4, *Z*<sup>4</sup> . The components are given as

$$\mathcal{U}\_4 = \frac{\omega \lambda \varepsilon\_5}{\omega \varepsilon\_1 \varepsilon\_5 + \eta\_2 \mu\_2 \varepsilon\_6}, \quad \mathcal{I}\_4^V = \frac{\varepsilon\_6}{\omega}, \quad \mathcal{V}\_4 = \frac{\mu\_2 \varepsilon\_6}{\omega \varepsilon\_5}, \quad \mathcal{Z}\_4 = \frac{\varepsilon\_3}{\gamma} (\mathcal{R}\_{1V} - 1),$$

where <sup>R</sup>1*<sup>V</sup>* <sup>=</sup> *ωλη*2*μ*<sup>2</sup> 3[*ω*1<sup>5</sup> + *η*2*μ*26] . We see that *U*4, *I<sup>V</sup>* <sup>4</sup> , *<sup>V</sup>*<sup>4</sup> <sup>&</sup>gt; 0, while *<sup>Z</sup>*<sup>4</sup> <sup>&</sup>gt; 0 if <sup>R</sup>1*<sup>V</sup>* <sup>&</sup>gt; 1. Hence, *<sup>E</sup>*<sup>4</sup> is defined when <sup>R</sup>1*<sup>V</sup>* <sup>&</sup>gt; 1. Here, the threshold <sup>R</sup>1*<sup>V</sup>* defines the stimulation status of CTL immunity versus SARS-CoV-2-infected EPCs.

(vi) The Mtb/SARS-CoV-2 coinfection equilibrium *E*<sup>5</sup> = *U*5, *L*5, *I<sup>B</sup>* <sup>5</sup> , *<sup>I</sup><sup>V</sup>* <sup>5</sup> , *B*5, *V*5, *Z*<sup>5</sup> . The components are defined as

*<sup>U</sup>*<sup>5</sup> <sup>=</sup> (<sup>2</sup> <sup>−</sup> 3)<sup>5</sup> *η*2*μ*2(*σ* − 1) , *<sup>L</sup>*<sup>5</sup> <sup>=</sup> *<sup>η</sup>*1*μ*12345(*ω*1<sup>5</sup> <sup>+</sup> *<sup>η</sup>*2*μ*26) *aω*(*η*1*μ*12<sup>5</sup> − *η*2*μ*24) 2 *ωλη*1*μ*12<sup>5</sup> 34(*ω*1<sup>5</sup> <sup>+</sup> *<sup>η</sup>*2*μ*26) <sup>+</sup> <sup>1</sup> <sup>−</sup> <sup>2</sup> 3 − R1*V* , *I B* <sup>5</sup> <sup>=</sup> 3(*ω*1<sup>5</sup> <sup>+</sup> *<sup>η</sup>*2*μ*26) *ω*(<sup>2</sup> − 3)*η*2*μ*2(*σ* − 1) *ωλη*1*μ*12<sup>5</sup> 34(*ω*1<sup>5</sup> <sup>+</sup> *<sup>η</sup>*2*μ*26) <sup>+</sup> <sup>1</sup> <sup>−</sup> <sup>2</sup> 3 − R1*V* , *I V* <sup>5</sup> <sup>=</sup> 25(*ω*1<sup>4</sup> <sup>+</sup> *<sup>η</sup>*1*μ*126) *ω*(<sup>2</sup> − 3)*η*2*μ*24(*σ* − 1) *ωλη*2*μ*2<sup>4</sup> 25(*ω*1<sup>4</sup> <sup>+</sup> *<sup>η</sup>*1*μ*126) <sup>+</sup> <sup>1</sup> <sup>−</sup> <sup>3</sup> 2 − R1*B* , *<sup>B</sup>*<sup>5</sup> <sup>=</sup> *<sup>μ</sup>*123(*ω*1<sup>5</sup> <sup>+</sup> *<sup>η</sup>*2*μ*26) *ω*(<sup>2</sup> − 3)*η*2*μ*24(*σ* − 1) *ωλη*1*μ*12<sup>5</sup> 34(*ω*1<sup>5</sup> <sup>+</sup> *<sup>η</sup>*2*μ*26) <sup>+</sup> <sup>1</sup> <sup>−</sup> <sup>2</sup> 3 − R1*V* , *<sup>V</sup>*<sup>5</sup> <sup>=</sup> 2(*ω*1<sup>4</sup> <sup>+</sup> *<sup>η</sup>*1*μ*126) *ω*(<sup>2</sup> − 3)*η*24(*σ* − 1) *ωλη*2*μ*2<sup>4</sup> 25(*ω*1<sup>4</sup> <sup>+</sup> *<sup>η</sup>*1*μ*126) <sup>+</sup> <sup>1</sup> <sup>−</sup> <sup>3</sup> 2 − R1*B* , *<sup>Z</sup>*<sup>5</sup> <sup>=</sup> *<sup>η</sup>*1*μ*1235(R0*V*/R0*<sup>B</sup>* <sup>−</sup> <sup>1</sup>) *γη*2*μ*24(*<sup>σ</sup>* <sup>−</sup> <sup>1</sup>) ,

where *<sup>σ</sup>* <sup>=</sup> *<sup>η</sup>*1*μ*12<sup>5</sup> *η*2*μ*2<sup>4</sup> . We see that *L*5, *I<sup>B</sup>* <sup>5</sup> , and *<sup>B</sup>*<sup>5</sup> are positive if *ωλη*1*μ*12<sup>5</sup> 34(*ω*1<sup>5</sup> <sup>+</sup> *<sup>η</sup>*2*μ*26) <sup>+</sup> <sup>1</sup> > <sup>2</sup> 3 <sup>+</sup> <sup>R</sup>1*V*, while *<sup>I</sup><sup>V</sup>* <sup>5</sup> and *<sup>V</sup>*<sup>5</sup> are positive if *ωλη*2*μ*2<sup>4</sup> 25(*ω*1<sup>4</sup> <sup>+</sup> *<sup>η</sup>*1*μ*126) <sup>+</sup> <sup>1</sup> <sup>&</sup>gt; <sup>3</sup> 2 + R1*B*, and *<sup>Z</sup>*<sup>5</sup> <sup>&</sup>gt; 0 if <sup>R</sup>0*<sup>V</sup>* R0*B* <sup>&</sup>gt; 1. In addition, we need the two conditions <sup>2</sup> <sup>&</sup>gt; <sup>3</sup> and *<sup>σ</sup>* <sup>&</sup>gt; 1. Hence, *E*<sup>5</sup> exists when the above conditions are met.

#### **4. Global Properties**

This part is aimed to prove the global stability of all equilibria by adopting correct Lyapunov functionals. The construction of these Lyapunov functionals follows the methods presented in [40–42].

We consider a function Ξ*i*(*U*, *L*, *IB*, *IV*, *B*, *V*, *Z*) and suppose that *χ <sup>i</sup>* is the largest invariant subset of *χ<sup>i</sup>* = (*U*, *<sup>L</sup>*, *<sup>I</sup>B*, *<sup>I</sup>V*, *<sup>B</sup>*, *<sup>V</sup>*, *<sup>Z</sup>*) <sup>|</sup> *<sup>d</sup>*Ξ*<sup>i</sup> dt* <sup>=</sup> <sup>0</sup> 6 , *i* = 0, 1, . . . , 5.

**Theorem 2.** *The equilibrium E*<sup>0</sup> *is globally asymptotically stable (GS) if* R0*<sup>B</sup>* ≤ 1 *and* R0*<sup>V</sup>* ≤ 1*.*

**Proof.** We opt a Lyapunov functional (LF)

$$\Xi\_0(t) = \int\_{\mathbb{T}} \Xi\_0(\mathbf{x}, t) \, d\mathbf{x}, \quad \text{where}$$

$$\tilde{\Xi}\_0 = l l\_0 \left( \frac{lI}{lI\_0} - 1 - \ln \frac{lI}{lI\_0} \right) + L + I^B + I^V + \frac{1}{\mu\_1} B + \frac{\varepsilon\_3}{\mu\_2} V + \frac{\gamma}{\omega} Z.$$

By taking the partial derivative, we obtain

*∂*Ξ0<sup>0</sup> *<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>U</sup>*<sup>0</sup> *U DU*Δ*<sup>U</sup>* <sup>+</sup> *<sup>λ</sup>* <sup>−</sup> *<sup>η</sup>*1*UB* <sup>−</sup> *<sup>η</sup>*2*UV* <sup>−</sup> 1*<sup>U</sup>* + *DL*Δ*L* + *η*1*UB* − *aL* + *DIB*Δ*I B* + *aL* − *γI BZ* <sup>−</sup> <sup>2</sup> *<sup>I</sup> <sup>B</sup>* + *DIV* <sup>Δ</sup>*<sup>I</sup> <sup>V</sup>* <sup>+</sup> *<sup>η</sup>*2*UV* <sup>−</sup> *<sup>γ</sup><sup>I</sup> <sup>V</sup> <sup>Z</sup>* <sup>−</sup> <sup>3</sup> *<sup>I</sup> <sup>V</sup>* + 1 *μ*1 *DB*Δ*B* + *μ*1<sup>2</sup> *I <sup>B</sup>* <sup>−</sup> 4*<sup>B</sup>* <sup>+</sup> <sup>3</sup> *μ*2 *DV*Δ*V* + *μ*<sup>2</sup> *I <sup>V</sup>* <sup>−</sup> 5*<sup>V</sup>* <sup>+</sup> *<sup>γ</sup> ω DZ*Δ*Z* + *ωI BZ* + *ωI <sup>V</sup> <sup>Z</sup>* <sup>−</sup> 6*<sup>Z</sup>* = <sup>1</sup> <sup>−</sup> *<sup>U</sup>*<sup>0</sup> *U* (*λ* − 1*U*) + *<sup>η</sup>*1*U*<sup>0</sup> <sup>−</sup> <sup>4</sup> *μ*1 *B* + *<sup>η</sup>*2*U*<sup>0</sup> <sup>−</sup> 3<sup>5</sup> *μ*2 *<sup>V</sup>* <sup>−</sup> *γ*<sup>6</sup> *<sup>ω</sup> <sup>Z</sup>* <sup>+</sup> <sup>1</sup> <sup>−</sup> *<sup>U</sup>*<sup>0</sup> *U DU*Δ *U* + *DL*Δ*L* + *DIB*Δ*I <sup>B</sup>* + *DIV* <sup>Δ</sup>*<sup>I</sup> <sup>V</sup>* + 1 *μ*1 *DB*Δ*<sup>B</sup>* <sup>+</sup> <sup>3</sup> *μ*2 *DV*Δ*<sup>V</sup>* <sup>+</sup> *<sup>γ</sup> <sup>ω</sup> DZ*Δ*Z*.

$$\begin{aligned} \text{The derivative } \frac{d\Xi\_0}{dt} \text{ is given by} \\ \frac{d\Xi\_0}{dt} = -\varepsilon\_1 \int\_\mathbf{\varUpsilon} \frac{(\varPi - \varPi\_0)^2}{\varPi} \, \mathbf{dx} + \frac{\epsilon\_4}{\mu\_1} (\varPi\_{0\varPi} - 1) \int\_\mathbf{\varUpsilon} \, \mathbf{B} \, \mathbf{dx} + \frac{\xi \mathbf{y} \xi\_5}{\mu\_2} (\varPi\_{0\varPi} - 1) \int\_\mathbf{\varUpsilon} \, \mathbf{J} \, \mathbf{dx} - \frac{\gamma \epsilon\_6}{\omega} \int\_\mathbf{\varUpsilon} \, \mathbf{Z} \, \mathbf{dx} \\ + D\_\mathbf{U} \int\_\mathbf{\varUpsilon} \left( 1 - \frac{\mathbf{U}\_0}{\mathbf{U}} \right) \Delta \mathbf{U} \, \mathbf{dx} + D\_\mathbf{L} \int\_\mathbf{\varUpsilon} \Delta \mathbf{L} \, \mathbf{dx} + D\_{I^3} \int\_\mathbf{\varUpsilon} \Delta \mathbf{P} \, \mathbf{dx} + D\_{I^V} \int\_\mathbf{\varUpsilon} \Delta \mathbf{V}^\mathbf{dx} \\ + \frac{1}{\mu\_1} D\_\mathbf{B} \int\_\mathbf{\varUpsilon} \, \Delta \mathbf{B} \, \mathbf{dx} + \frac{\xi\_3}{\mu\_2} D\_V \int\_\mathbf{\varUpsilon} \, \Delta V \, \mathbf{dx} + \frac{\gamma}{\omega} D\_{\mathbf{Z}} \int\_\mathbf{\varUpsilon} \Delta \mathbf{Z} \, \mathbf{dx}. \end{aligned} \tag{4}$$

Based on the Divergence theorem and Neumann BCs, we obtain

$$\begin{split} 0 &= \int\_{\partial \mathbb{F}} \nabla \Phi \cdot \mathbb{F} \, \mathrm{d}x = \int\_{\mathbb{F}} \mathrm{div} (\nabla \Phi) \, \mathrm{d}x = \int\_{\mathbb{F}} \Delta \Phi \, \mathrm{d}x, \\ 0 &= \int\_{\partial \mathbb{F}} \frac{1}{\Phi} \nabla \Phi \cdot \mathbb{F} \, \mathrm{d}x = \int\_{\mathbb{F}} \mathrm{div} (\frac{1}{\Phi} \nabla \Phi) \, \mathrm{d}x = \int\_{\mathbb{F}} \left[ \frac{\Delta \Phi}{\Phi} - \frac{\|\nabla \Phi\|^2}{\Phi^2} \right] \, \mathrm{d}x, \quad \text{for } \Phi \in \{\mathcal{U}, L, I^B, I^V, B, V, Z\}. \end{split} \tag{5}$$

As a result, the derivative in (4) is altered to

$$\frac{d\Xi\_0}{dt} = -\epsilon\_1 \int\_{\mathbf{Y}} \frac{\left(\mathcal{U} - \mathcal{U}\right)^2}{\mathcal{U}} \, \mathrm{d}x + \frac{\epsilon\_4}{\mu\_1} (\mathcal{R}\_{0B} - 1) \int\_{\mathbf{Y}} \mathcal{B} \, \mathrm{d}x + \frac{\varepsilon\_3 \mathfrak{c} \mathfrak{s}}{\mu\_2} (\mathcal{R}\_{0V} - 1) \int\_{\mathbf{Y}} \mathcal{V} \, \mathrm{d}x - \frac{\gamma \mathfrak{c} \mathfrak{s}}{\omega} \int\_{\mathbf{Y}} \mathcal{Z} \, \mathrm{d}x$$

$$-D\_{\mathcal{U}} \mathcal{U}\_0 \int\_{\mathbf{Y}} \frac{\left\| \mathcal{V} \mathcal{U} \right\|^2}{\mathcal{U}^2} \, \mathrm{d}x.$$

$$\text{We see that } \frac{d\Xi\_0}{dt} \le 0 \text{ if } \mathcal{R}\_{0B} \le 1 \text{ and } \mathcal{R}\_{0V} \le 1. \text{ Furthermore, } \frac{d\Xi\_0}{dt} = 0 \text{ if } \mathcal{U} = \mathcal{U}\_0.$$

and *B* = *V* = *Z* = 0. The solutions approach *χ* <sup>0</sup> that has *<sup>B</sup>* <sup>=</sup> *<sup>V</sup>* <sup>=</sup> 0. Thus, *<sup>∂</sup><sup>B</sup> <sup>∂</sup><sup>t</sup>* <sup>=</sup> 0 and *∂V <sup>∂</sup><sup>t</sup>* <sup>=</sup> 0. According to the fifth and sixth equations of system (1), we acquire *<sup>I</sup><sup>B</sup>* <sup>=</sup> *<sup>I</sup><sup>V</sup>* <sup>=</sup> 0. Therefore, *<sup>∂</sup>I<sup>B</sup> <sup>∂</sup><sup>t</sup>* <sup>=</sup> 0 and the third equation of (1) gives *<sup>L</sup>* <sup>=</sup> 0. Consequently, *<sup>χ</sup>* <sup>0</sup> = {*E*0} and in compliance with LaSalle's invariance principle (LIP) [43], the point *E*<sup>0</sup> is GS if R0*<sup>B</sup>* ≤ 1 and R0*<sup>V</sup>* ≤ 1.

**Theorem 3.** *Let* <sup>R</sup>0*<sup>B</sup>* <sup>&</sup>gt; <sup>1</sup>*. Then, the equilibrium E*<sup>1</sup> *is GS if* <sup>R</sup>0*<sup>V</sup>* R0*B* ≤ 1 *with* R1*<sup>B</sup>* ≤ 1*.*

**Proof.** We adopt an LF

$$\Xi\_1(t) = \int\_{\mathbf{F}} \tilde{\Xi}\_1(\mathbf{x}, t) \, d\mathbf{x}, \quad \text{where}$$

$$\begin{split} \tilde{\Xi}\_1 &= lI\_1\left(\frac{lI}{lI\_1} - 1 - \ln\frac{lI}{lI\_1}\right) + L\_1\left(\frac{L}{L\_1} - 1 - \ln\frac{L}{L\_1}\right) + I\_1^B\left(\frac{I^B}{I\_1^B} - 1 - \ln\frac{I^B}{I\_1^B}\right) + I^V \\ &+ \frac{1}{\mu\_1}B\_1\left(\frac{B}{B\_1} - 1 - \ln\frac{B}{B\_1}\right) + \frac{\varepsilon\_3}{\mu\_2}V + \frac{\gamma}{\omega}Z. \end{split}$$

By calculating the partial derivative, we obtain

$$\begin{split} \frac{\partial \Xi\_{1}}{\partial t} &= \left(1 - \frac{\mathcal{U}\_{1}}{\mathcal{U}}\right) \left(D\_{\mathcal{U}}\Delta\mathcal{U} + \lambda - \eta\_{1}\mathcal{U}B - \eta\_{2}\mathcal{U}\mathcal{V} - \varepsilon\_{1}\mathcal{U}\right) + \left(1 - \frac{\mathcal{L}\_{1}}{\mathcal{L}}\right) \left(D\_{\mathcal{U}}\Delta\mathcal{L} + \eta\_{1}\mathcal{U}B - a\mathcal{L}\right) \\ &+ \left(1 - \frac{\mathcal{I}\_{1}^{B}}{\mathcal{I}^{B}}\right) \left(D\_{\mathcal{V}}\Delta\mathcal{U}^{B} + a\mathcal{L} - \gamma\mathcal{I}^{B}\mathcal{Z} - \varepsilon\_{2}\mathcal{I}^{B}\right) + D\_{\mathcal{V}}\Delta\mathcal{U}^{\mathcal{V}} + \eta\_{2}\mathcal{U}\mathcal{V} - \gamma\mathcal{I}^{\mathcal{V}}\mathcal{Z} - \varepsilon\_{3}\mathcal{I}^{\mathcal{V}} \\ &+ \frac{1}{\mu\_{1}}\left(1 - \frac{B\_{1}}{B}\right) \left(D\_{\mathcal{B}}\Delta\mathcal{B} + \mu\_{1}\varepsilon\_{2}\mathcal{I}^{B} - \epsilon\_{4}\mathcal{B}\right) + \frac{\varepsilon\_{3}}{\mu\_{2}}\left(D\_{\mathcal{V}}\Delta\mathcal{V} + \mu\_{2}\mathcal{I}^{\mathcal{V}} - \epsilon\_{5}\mathcal{V}\right) \\ &+ \frac{\gamma}{\omega}\left(D\_{\mathcal{Z}}\Delta\mathcal{Z} + \omega\mathcal{U}^{B}\mathcal{Z} + \omega\mathcal{U}^{\mathcal{V}}\mathcal{Z} - \epsilon\_{6}\mathcal{Z}\right). \end{split} \tag{6}$$

By employing the equilibrium conditions at *E*1, we obtain

$$\begin{cases} \lambda = \eta\_1 \mathcal{U}\_1 B\_1 + \varepsilon\_1 \mathcal{U}\_1, \\ \eta\_1 \mathcal{U}\_1 B\_1 = a L\_1, \\ a L\_1 = \varepsilon\_2 I\_1^B, \\ \epsilon\_2 I\_1^B = \frac{\varepsilon\_4}{\mu\_1} B\_1. \end{cases}$$

Hence, the derivative in (6) can be simplified to

$$\begin{split} \frac{\partial \Xi\_{1}}{\partial t} &= \left(1 - \frac{lL\_{1}}{\mathcal{U}}\right) (\varepsilon\_{1}lI\_{1} - \varepsilon\_{1}lI) + \eta\_{1}lI\_{1}B\_{1}\left(4 - \frac{lI\_{1}}{\mathcal{U}} - \frac{lIL\_{1}B}{\mathcal{U}\_{1}L\_{1}B} - \frac{LI\_{1}^{B}}{L\_{1}I^{B}} - \frac{l^{B}B\_{1}}{I\_{1}^{B}B}\right) + \left(\eta\_{2}\Delta I\_{1} - \frac{\varepsilon\_{2}\eta\_{2}}{\mu\_{2}}\right)V \\ &+ \left(\gamma I\_{1}^{B} - \frac{\gamma\varepsilon\_{6}}{\omega}\right)Z + \left(1 - \frac{lI\_{1}}{\mathcal{U}}\right)D\_{\mathcal{U}}\Delta I + \left(1 - \frac{L\_{1}}{L}\right)D\_{\mathcal{U}}\Delta L + \left(1 - \frac{l\_{1}^{B}}{I\_{B}^{B}}\right)D\_{\mathcal{U}}\Delta I^{B} \\ &+ D\_{\mathcal{U}}\Delta I^{V} + \frac{1}{\mu\_{1}}\left(1 - \frac{B\_{1}}{B}\right)D\_{B}\Delta B + \frac{\varepsilon\_{3}}{\mu\_{2}}D\_{V}\Delta V + \frac{\gamma}{\omega}D\_{Z}\Delta Z. \end{split}$$

By using (5), *<sup>d</sup>*Ξ<sup>1</sup> *dt* is given by

*d*Ξ<sup>1</sup> *dt* <sup>=</sup> <sup>−</sup> <sup>1</sup> Ψ (*U* − *U*1) 2 *<sup>U</sup>* <sup>d</sup>*<sup>x</sup>* <sup>+</sup> *<sup>η</sup>*1*U*1*B*<sup>1</sup> Ψ <sup>4</sup> <sup>−</sup> *<sup>U</sup>*<sup>1</sup> *<sup>U</sup>* <sup>−</sup> *UL*1*<sup>B</sup> U*1*LB*<sup>1</sup> <sup>−</sup> *LI<sup>B</sup>* 1 *<sup>L</sup>*<sup>1</sup> *<sup>I</sup><sup>B</sup>* <sup>−</sup> *<sup>I</sup>BB*<sup>1</sup> *IB* 1 *B* <sup>d</sup>*<sup>x</sup>* <sup>+</sup> 3<sup>5</sup> *μ*2 R0*<sup>V</sup>* R0*B* − 1 Ψ *V* d*x* <sup>+</sup> *<sup>γ</sup>*(*ω*1<sup>4</sup> <sup>+</sup> *<sup>η</sup>*1*μ*126) *ωη*1*μ*1<sup>2</sup> (R1*<sup>B</sup>* − 1) Ψ *Z* d*x* − *DUU*<sup>1</sup> Ψ *U*<sup>2</sup> *<sup>U</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DLL*<sup>1</sup> Ψ *L*<sup>2</sup> *<sup>L</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* − *DIB I B* 1 Ψ *IB*<sup>2</sup> *<sup>I</sup>B*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DBB*<sup>1</sup> *μ*1 Ψ *B*<sup>2</sup> *<sup>B</sup>*<sup>2</sup> <sup>d</sup>*x*.

In this situation, *<sup>d</sup>*Ξ<sup>1</sup> *dt* <sup>≤</sup> 0 if <sup>R</sup>0*<sup>V</sup>* R0*B* <sup>≤</sup> 1 and <sup>R</sup>1*<sup>B</sup>* <sup>≤</sup> 1. In addition, *<sup>d</sup>*Ξ<sup>1</sup> *dt* <sup>=</sup> 0 when *U* = *U*1, *L* = *L*1, *I<sup>B</sup>* = *I<sup>B</sup>* <sup>1</sup> , *B* = *B*1, while *V* = *Z* = 0. The solutions tend to *χ* <sup>1</sup> with *V* = 0 and therefore *<sup>∂</sup><sup>V</sup> <sup>∂</sup><sup>t</sup>* <sup>=</sup> 0. The sixth equation of (1) yields *<sup>I</sup><sup>V</sup>* <sup>=</sup> 0. Hence, *<sup>χ</sup>* <sup>1</sup> = {*E*1} and *E*<sup>1</sup> is GS when <sup>R</sup>0*<sup>B</sup>* <sup>&</sup>gt; 1, <sup>R</sup>0*<sup>V</sup>* R0*B* ≤ 1 and R1*<sup>B</sup>* ≤ 1 according to LIP [43].

**Theorem 4.** *Let* <sup>R</sup>0*<sup>V</sup>* <sup>&</sup>gt; <sup>1</sup>*. The equilibrium E*<sup>2</sup> *is GS if* <sup>R</sup>0*<sup>B</sup>* R0*V* ≤ 1 *and* R1*<sup>V</sup>* ≤ 1*.*

**Proof.** We pick an LF

$$
\Xi\_2(t) = \int\_{\mathbb{F}} \vec{\Xi}\_2(\mathfrak{x}, t) \, \mathrm{d}\mathfrak{x}, \quad \text{where}
$$

Ξ0<sup>2</sup> =*U*<sup>2</sup> *U U*<sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>U</sup> U*<sup>2</sup> + *L* + *I <sup>B</sup>* + *I V* 2 *IV IV* 2 <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>I</sup><sup>V</sup> IV* 2 + 1 *μ*1 *<sup>B</sup>* <sup>+</sup> <sup>3</sup> *μ*2 *V*2 *V V*2 <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>V</sup> V*2 <sup>+</sup> *<sup>γ</sup> <sup>ω</sup> <sup>Z</sup>*. Then, *<sup>∂</sup>*Ξ0<sup>2</sup> *<sup>∂</sup><sup>t</sup>* is computed as *∂*Ξ0<sup>2</sup> *<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>U</sup>*<sup>2</sup> *U DU*Δ*<sup>U</sup>* <sup>+</sup> *<sup>λ</sup>* <sup>−</sup> *<sup>η</sup>*1*UB* <sup>−</sup> *<sup>η</sup>*2*UV* <sup>−</sup> 1*<sup>U</sup>* + *DL*Δ*L* + *η*1*UB* − *aL* + *DIB*Δ*I <sup>B</sup>* + *aL* − *γI BZ* <sup>−</sup> <sup>2</sup> *<sup>I</sup> <sup>B</sup>* + <sup>1</sup> <sup>−</sup> *<sup>I</sup><sup>V</sup>* 2 *IV DIV* Δ*I <sup>V</sup>* <sup>+</sup> *<sup>η</sup>*2*UV* <sup>−</sup> *<sup>γ</sup><sup>I</sup> <sup>V</sup> <sup>Z</sup>* <sup>−</sup> <sup>3</sup> *<sup>I</sup> V* + 1 *μ*1 *DB*Δ*B* + *μ*1<sup>2</sup> *I <sup>B</sup>* <sup>−</sup> 4*<sup>B</sup>* <sup>+</sup> <sup>3</sup> *μ*2 <sup>1</sup> <sup>−</sup> *<sup>V</sup>*<sup>2</sup> *V DV*Δ*<sup>V</sup>* <sup>+</sup> *<sup>μ</sup>*<sup>2</sup> *<sup>I</sup> <sup>V</sup>* <sup>−</sup> 5*<sup>V</sup>* <sup>+</sup> *<sup>γ</sup> ω DZ*Δ*Z* + *ωI BZ* + *ωI <sup>V</sup> <sup>Z</sup>* <sup>−</sup> 6*<sup>Z</sup>* . (7)

By considering the equilibrium conditions at *E*<sup>2</sup>

$$\begin{cases} \lambda = \eta\_2 \mathcal{U}\_2 V\_2 + \epsilon\_1 \mathcal{U}\_{2\prime} \\ \eta\_2 \mathcal{U}\_2 V\_2 = \epsilon\_3 I\_2^V \\ \epsilon\_3 I\_2^V = \frac{\epsilon\_3 \epsilon\_5}{\mu\_2} V\_{2\prime} \end{cases}$$

the derivative in (7) becomes

$$\begin{split} \frac{\partial \Xi\_2}{\partial t} &= \left(1 - \frac{lL\_2}{lI}\right) (\varepsilon\_1 l L\_2 - \varepsilon\_1 l I) + \eta\_2 l L\_2 V\_2 \left(3 - \frac{lL\_2}{lI} - \frac{lII\_2^V V}{lL\_2 I^V V\_2} - \frac{l^V V\_2}{l\_2^V V}\right) + \left(\eta\_1 l L\_2 - \frac{\varepsilon\_4}{\mu\_1}\right) B \\ &+ \left(\gamma I\_2^V - \frac{\gamma \varepsilon\_6}{\omega}\right) Z + \left(1 - \frac{lL\_2}{lI}\right) D\_{lI} \Delta l I + D\_L \Delta L + D\_{I^B} \Delta I^B + \left(1 - \frac{l^V}{I^V}\right) D\_{I^V} \Delta I^V \\ &+ \frac{1}{\mu\_1} D\_B \Delta B + \frac{\varepsilon\_3}{\mu\_2} \left(1 - \frac{V\_2}{V}\right) D\_{V} \Delta V + \frac{\gamma}{\omega} D\_Z \Delta Z. \end{split}$$

By using (5), the derivative of Ξ2(*t*) is expressed as

*d*Ξ<sup>2</sup> *dt* <sup>=</sup> <sup>−</sup> <sup>1</sup> Ψ (*U* − *U*2) 2 *<sup>U</sup>* <sup>d</sup>*<sup>x</sup>* <sup>+</sup> *<sup>η</sup>*2*U*2*V*<sup>2</sup> Ψ <sup>3</sup> <sup>−</sup> *<sup>U</sup>*<sup>2</sup> *<sup>U</sup>* <sup>−</sup> *U I<sup>V</sup>* <sup>2</sup> *V U*<sup>2</sup> *IVV*<sup>2</sup> <sup>−</sup> *<sup>I</sup>VV*<sup>2</sup> *IV* <sup>2</sup> *V* <sup>d</sup>*<sup>x</sup>* <sup>+</sup> <sup>4</sup> *μ*1 <sup>R</sup>0*<sup>B</sup>* R0*V* − 1 Ψ *B* d*x* <sup>+</sup> *<sup>γ</sup>*(*ω*1<sup>5</sup> <sup>+</sup> *<sup>η</sup>*2*μ*26) *ωη*2*μ*<sup>2</sup> (R1*<sup>V</sup>* − 1) Ψ *Z* d*x* − *DUU*<sup>2</sup> Ψ *U*<sup>2</sup> *<sup>U</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* − *DIV I V* 2 Ψ *IV*<sup>2</sup> *<sup>I</sup>V*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DVV*2<sup>3</sup> *μ*2 Ψ *V*<sup>2</sup> *<sup>V</sup>*<sup>2</sup> <sup>d</sup>*x*. We note that *<sup>d</sup>*Ξ<sup>2</sup> *dt* <sup>≤</sup> 0 if <sup>R</sup>0*<sup>B</sup>* R0*V* <sup>≤</sup> 1, and <sup>R</sup>1*<sup>V</sup>* <sup>≤</sup> 1. Moreover, *<sup>d</sup>*Ξ<sup>2</sup> *dt* <sup>=</sup> 0 when *<sup>U</sup>* <sup>=</sup> *<sup>U</sup>*2, 

*I<sup>V</sup>* = *I<sup>V</sup>* <sup>2</sup> , *V* = *V*<sup>2</sup> and *B* = *Z* = 0. The solutions approach *χ* <sup>2</sup>, which has an element with *<sup>B</sup>* <sup>=</sup> 0 and hence *<sup>∂</sup><sup>B</sup> <sup>∂</sup><sup>t</sup>* <sup>=</sup> 0. From the fifth equation of (1), we have *<sup>I</sup><sup>B</sup>* <sup>=</sup> 0. Consequently, *∂I<sup>B</sup> <sup>∂</sup><sup>t</sup>* <sup>=</sup> 0 and thus *<sup>L</sup>* <sup>=</sup> 0 according to the third equation of (1). Thereupon, *<sup>χ</sup>* <sup>2</sup> = {*E*2} and *<sup>E</sup>*<sup>2</sup> is GS when <sup>R</sup>0*<sup>V</sup>* <sup>&</sup>gt; 1, <sup>R</sup>0*<sup>B</sup>* R0*V* ≤ 1 and R1*<sup>V</sup>* ≤ 1 as attributed to LIP [43].

**Theorem 5.** *Assume that* <sup>R</sup>1*<sup>B</sup>* <sup>&</sup>gt; <sup>1</sup>*. Then, the equilibrium <sup>E</sup>*<sup>3</sup> *is GS if λωη*2*μ*2<sup>4</sup> 25[*ω*1<sup>4</sup> + *η*1*μ*126] + <sup>1</sup> <sup>≤</sup> <sup>3</sup> 2 + R1*B.*

**Proof.** We pick an LF

$$\Xi\_3(t) = \int\_{\mathbf{Y}} \tilde{\Xi}\_3(\mathbf{x}, t) \, d\mathbf{x}, \quad \text{where}$$

$$\begin{aligned} \tilde{\Xi}\_3 &= lI\_3\left(\frac{\underline{U}}{\underline{U}\_3} - 1 - \ln\frac{\underline{U}}{\underline{U}\_3}\right) + lJ\_3\left(\frac{\underline{L}}{\underline{L}\_3} - 1 - \ln\frac{\underline{L}}{\underline{L}\_3}\right) + l^B\_3 \left(\frac{\underline{I}^B}{\underline{I}\_3^B} - 1 - \ln\frac{\underline{I}^B}{\underline{I}\_3^B}\right) + \underline{I}^V\\ &+ \left(\frac{1}{\mu\_1} + \frac{\gamma Z\_3}{\mu\_1 \varepsilon\_2}\right)B\_3\left(\frac{B}{\underline{B}\_3} - 1 - \ln\frac{B}{\underline{B}\_3}\right) + \left(\frac{\epsilon\_3}{\mu\_2} + \frac{\gamma Z\_3}{\mu\_2}\right)V + \frac{\gamma}{\omega}Z\_3\left(\frac{Z}{\underline{Z}\_3} - 1 - \ln\frac{Z}{\underline{Z}\_3}\right), \end{aligned}$$
Then,  $\frac{\partial \Xi\_3}{\partial t}$  is written as

*∂*Ξ0<sup>3</sup> *<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>U</sup>*<sup>3</sup> *U DU*Δ*<sup>U</sup>* <sup>+</sup> *<sup>λ</sup>* <sup>−</sup> *<sup>η</sup>*1*UB* <sup>−</sup> *<sup>η</sup>*2*UV* <sup>−</sup> 1*<sup>U</sup>* + <sup>1</sup> <sup>−</sup> *<sup>L</sup>*<sup>3</sup> *L DL*Δ*<sup>L</sup>* <sup>+</sup> *<sup>η</sup>*1*UB* <sup>−</sup> *aL* + <sup>1</sup> <sup>−</sup> *<sup>I</sup><sup>B</sup>* 3 *IB DIB*Δ*I <sup>B</sup>* <sup>+</sup> *aL* <sup>−</sup> *<sup>γ</sup><sup>I</sup> BZ* <sup>−</sup> <sup>2</sup> *<sup>I</sup> B* + *DIV* Δ*I <sup>V</sup>* <sup>+</sup> *<sup>η</sup>*2*UV* <sup>−</sup> *<sup>γ</sup><sup>I</sup> <sup>V</sup> <sup>Z</sup>* <sup>−</sup> <sup>3</sup> *<sup>I</sup> V* + 1 *μ*1 <sup>+</sup> *<sup>γ</sup>Z*<sup>3</sup> *μ*1<sup>2</sup> <sup>1</sup> <sup>−</sup> *<sup>B</sup>*<sup>3</sup> *B DB*Δ*<sup>B</sup>* <sup>+</sup> *<sup>μ</sup>*1<sup>2</sup> *<sup>I</sup> <sup>B</sup>* <sup>−</sup> 4*<sup>B</sup>* + <sup>3</sup> *μ*2 <sup>+</sup> *<sup>γ</sup>Z*<sup>3</sup> *μ*2 *DV*Δ*<sup>V</sup>* <sup>+</sup> *<sup>μ</sup>*<sup>2</sup> *<sup>I</sup> <sup>V</sup>* <sup>−</sup> 5*<sup>V</sup>* <sup>+</sup> *<sup>γ</sup> ω* <sup>1</sup> <sup>−</sup> *<sup>Z</sup>*<sup>3</sup> *Z DZ*Δ*<sup>Z</sup>* <sup>+</sup> *<sup>ω</sup><sup>I</sup> BZ* + *ωI <sup>V</sup> <sup>Z</sup>* <sup>−</sup> 6*<sup>Z</sup>* . (8)

By utilizing the equilibrium requirements at *E*<sup>3</sup> to add the terms of Equation (8)

$$\begin{cases} \lambda = \eta\_1 \mathcal{U}\_3 B\_3 + \epsilon\_1 \mathcal{U}\_{3\prime} \\ \eta\_1 \mathcal{U}\_3 B\_3 = a \mathcal{L}\_{3\prime} \\ a \mathcal{L}\_3 = \gamma I\_3^B Z\_3 + \epsilon\_2 I\_3^B \mathcal{L}\_{3\prime} \\ \epsilon\_2 I\_3^B = \frac{\varepsilon\_4}{\mu\_1} B\_{3\prime} \\ \gamma I\_3^B Z\_3 = \frac{\gamma \varepsilon\_6}{\omega} Z\_{3\prime} \end{cases}$$

we obtain

$$\begin{split} \frac{\partial \Xi\_{3}}{\partial t} &= \left(1 - \frac{lL\_{3}}{lI}\right) (\varepsilon\_{1}l\_{3} - \varepsilon\_{1}lI) + \eta\_{1}l\_{3}B\_{3} \left(4 - \frac{lL\_{3}}{lI} - \frac{lL\_{3}B}{lL\_{3}L\_{3}} - \frac{Ll\_{3}^{B}}{L\_{3}I^{B}} - \frac{I^{B}B\_{3}}{l\_{3}^{B}B}\right) \\ &+ \left(\eta\_{2}lL\_{3} - \frac{\varepsilon\_{3}\varepsilon\_{5}}{\mu\_{2}} - \frac{\gamma\varepsilon\_{5}Z\_{3}}{\mu\_{2}}\right)V + \left(1 - \frac{lL\_{3}}{lI}\right)D\_{l}\Delta lI + \left(1 - \frac{L\_{3}}{L}\right)D\_{L}\Delta L + \left(1 - \frac{l^{B}\_{3}}{I^{B}}\right)D\_{I}\Delta I^{B} \\ &+ D\_{I^{B}}\Delta I^{V} + \left(\frac{1}{\mu\_{1}} + \frac{\gamma Z\_{3}}{\mu\_{1}\varepsilon\_{2}}\right)\left(1 - \frac{B\_{3}}{B}\right)D\_{B}\Delta B + \left(\frac{\varepsilon\_{3}}{\mu\_{2}} + \frac{\gamma Z\_{3}}{\mu\_{2}}\right)D\_{V}\Delta I \\ &+ \frac{\gamma}{\omega}\left(1 - \frac{Z\_{3}}{Z}\right)D\_{Z}\Delta Z. \end{split}$$

By using (5), the derivative of Ξ3(*t*) is presented as

*d*Ξ<sup>3</sup> *dt* <sup>=</sup> <sup>−</sup> <sup>1</sup> Ψ (*U* − *U*3) 2 *<sup>U</sup>* <sup>d</sup>*<sup>x</sup>* <sup>+</sup> *<sup>η</sup>*1*U*3*B*<sup>3</sup> Ψ <sup>4</sup> <sup>−</sup> *<sup>U</sup>*<sup>3</sup> *<sup>U</sup>* <sup>−</sup> *UL*3*<sup>B</sup> U*3*LB*<sup>3</sup> <sup>−</sup> *LI<sup>B</sup>* 3 *<sup>L</sup>*<sup>3</sup> *<sup>I</sup><sup>B</sup>* <sup>−</sup> *<sup>I</sup>BB*<sup>3</sup> *IB* 3 *B* d*x* <sup>+</sup> 2<sup>5</sup> *μ*2 *λωη*2*μ*2<sup>4</sup> 25[*ω*1<sup>4</sup> + *η*1*μ*126] <sup>+</sup> <sup>1</sup> <sup>−</sup> <sup>3</sup> 2 − R1*B* Ψ *V* d*x* − *DUU*<sup>3</sup> Ψ *U*<sup>2</sup> *<sup>U</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* − *DLL*<sup>3</sup> Ψ *L*<sup>2</sup> *<sup>L</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DIB <sup>I</sup> B* 3 Ψ *IB*<sup>2</sup> *<sup>I</sup>B*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DBB*<sup>3</sup> 1 *μ*1 <sup>+</sup> *<sup>γ</sup>Z*<sup>3</sup> *μ*1<sup>2</sup> Ψ *B*<sup>2</sup> *<sup>B</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *<sup>γ</sup>DZZ*<sup>3</sup> *ω* Ψ *Z*<sup>2</sup> *<sup>Z</sup>*<sup>2</sup> <sup>d</sup>*x*.

We see that *<sup>d</sup>*Ξ<sup>3</sup> *dt* <sup>≤</sup> 0 if *λωη*2*μ*2<sup>4</sup> 25[*ω*1<sup>4</sup> + *η*1*μ*126] <sup>+</sup> <sup>1</sup> <sup>≤</sup> <sup>3</sup> 2 + R1*B*. In addition, it is possible to show that *<sup>d</sup>*Ξ<sup>3</sup> *dt* <sup>=</sup> 0 when (*U*, *<sup>L</sup>*, *<sup>I</sup>B*, *<sup>I</sup>V*, *<sup>B</sup>*, *<sup>V</sup>*, *<sup>Z</sup>*)=(*U*3, *<sup>L</sup>*3, *<sup>I</sup><sup>B</sup>* <sup>3</sup> , 0, *B*3, 0, *Z*3). Then, *χ* <sup>3</sup> <sup>=</sup> {*E*3} and in reference to LIP [43], *<sup>E</sup>*<sup>3</sup> is GS when <sup>R</sup>1*<sup>B</sup>* <sup>&</sup>gt; 1 and *λωη*2*μ*2<sup>4</sup> 25[*ω*1<sup>4</sup> + *η*1*μ*126] + <sup>1</sup> <sup>≤</sup> <sup>3</sup> 2 + R1*B*.

**Theorem 6.** *Let* <sup>R</sup>1*<sup>V</sup>* <sup>&</sup>gt; <sup>1</sup>*. Thereupon, the equilibrium <sup>E</sup>*<sup>4</sup> *is GS if λωη*1*μ*12<sup>5</sup> 34[*ω*1<sup>5</sup> + *η*2*μ*26] + 1 ≤ 2 3 + R1*V.*

**Proof.** We nominate an LF

$$\Xi\_4(t) = \int\_{\mathbb{T}} \tilde{\Xi}\_4(\mathbf{x}, t) \, d\mathbf{x}, \quad \text{where}$$

$$\begin{split} \tilde{\Xi}\_4 &= \mathcal{U}\_4 \left( \frac{\mathcal{U}}{\mathcal{U}\_4} - 1 - \ln \frac{\mathcal{U}}{\mathcal{U}\_4} \right) + L + I^B + I\_4^V \left( \frac{I^V}{I\_4^V} - 1 - \ln \frac{I^V}{I\_4^V} \right) + \left( \frac{1}{\mu\_1} + \frac{\gamma Z\_4}{\mu\_1 \varepsilon\_2} \right) B^\perp \\ &+ \left( \frac{\varepsilon\_3}{\mu\_2} + \frac{\gamma Z\_4}{\mu\_2} \right) V\_4 \left( \frac{V}{V\_4} - 1 - \ln \frac{V}{V\_4} \right) + \frac{\gamma}{\omega} Z\_4 \left( \frac{Z}{Z\_4} - 1 - \ln \frac{Z}{Z\_4} \right). \end{split}$$

By computing the partial derivative, we obtain

*∂*Ξ<sup>4</sup> *<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>U</sup>*<sup>4</sup> *U DU*Δ*<sup>U</sup>* <sup>+</sup> *<sup>λ</sup>* <sup>−</sup> *<sup>η</sup>*1*UB* <sup>−</sup> *<sup>η</sup>*2*UV* <sup>−</sup> 1*<sup>U</sup>* + *DL*Δ*L* + *η*1*UB* − *aL* + *DIB*Δ*I <sup>B</sup>* + *aL* − *γI BZ* <sup>−</sup> <sup>2</sup> *<sup>I</sup> <sup>B</sup>* + <sup>1</sup> <sup>−</sup> *<sup>I</sup><sup>V</sup>* 4 *IV DIV* Δ*I <sup>V</sup>* <sup>+</sup> *<sup>η</sup>*2*UV* <sup>−</sup> *<sup>γ</sup><sup>I</sup> <sup>V</sup> <sup>Z</sup>* <sup>−</sup> <sup>3</sup> *<sup>I</sup> V* + 1 *μ*1 <sup>+</sup> *<sup>γ</sup>Z*<sup>4</sup> *μ*1<sup>2</sup> *DB*Δ*<sup>B</sup>* <sup>+</sup> *<sup>μ</sup>*1<sup>2</sup> *<sup>I</sup> <sup>B</sup>* <sup>−</sup> 4*<sup>B</sup>* + <sup>3</sup> *μ*2 <sup>+</sup> *<sup>γ</sup>Z*<sup>4</sup> *μ*2 <sup>1</sup> <sup>−</sup> *<sup>V</sup>*<sup>4</sup> *V DV*Δ*<sup>V</sup>* <sup>+</sup> *<sup>μ</sup>*<sup>2</sup> *<sup>I</sup> <sup>V</sup>* <sup>−</sup> 5*<sup>V</sup>* <sup>+</sup> *<sup>γ</sup> ω* <sup>1</sup> <sup>−</sup> *<sup>Z</sup>*<sup>4</sup> *Z DZ*Δ*<sup>Z</sup>* <sup>+</sup> *<sup>ω</sup><sup>I</sup> BZ* + *ωI <sup>V</sup> <sup>Z</sup>* <sup>−</sup> 6*<sup>Z</sup>* . (9)

By considering the equilibrium conditions at *E*<sup>4</sup>

$$\begin{cases} \lambda = \eta\_2 \mathcal{U}\_4 V\_4 + \epsilon\_1 \mathcal{U}\_4 \\ \eta\_2 \mathcal{U}\_4 V\_4 = \gamma I\_4^V Z\_4 + \epsilon\_3 I\_4^V \\ \epsilon\_3 I\_4^V = \frac{\epsilon\_3 \epsilon\_5}{\mu\_2} V\_{4\prime} \\ \gamma I\_4^V Z\_4 = \frac{\gamma \epsilon\_6}{\omega} Z\_{4\prime} \end{cases}$$

the derivative in (9) is transformed to

$$\begin{split} \frac{\partial \Xi\_{4}}{\partial t} &= \left(1 - \frac{\mathcal{U}\_{4}}{\mathcal{U}}\right) (\varepsilon\_{1} \mathcal{U}\_{4} - \varepsilon\_{1} \mathcal{U}) + \eta\_{2} \mathcal{U}\_{4} V\_{4} \left(3 - \frac{\mathcal{U}\_{4}}{\mathcal{U}} - \frac{\mathcal{U} \mathcal{U}\_{4} \mathcal{V}}{\mathcal{U}\_{4} \mathcal{V} \mathcal{V}\_{4}} - \frac{\mathcal{V} \mathcal{V}\_{4}}{\mathcal{U}\_{4} \mathcal{V}}\right) + \left(\eta\_{1} \mathcal{U}\_{4} - \frac{\varepsilon\_{4}}{\mu\_{1}} - \frac{\gamma \varepsilon\_{4} \mathcal{Z}\_{4}}{\mu\_{1} \varepsilon\_{2}}\right) B \mathcal{V}\_{4} \\ &+ \left(1 - \frac{\mathcal{U}\_{4}}{\mathcal{U}}\right) D\_{\mathcal{U}} \Delta \mathcal{U} + D\_{\mathcal{L}} \Delta L + D\_{\mathcal{I}^{B}} \Delta I^{B} + \left(1 - \frac{\mathcal{V}\_{4}}{\mathcal{V}}\right) D\_{\mathcal{V}^{B}} \Delta \mathcal{V}^{V} + \left(\frac{1}{\mu\_{1}} + \frac{\gamma \mathcal{Z}\_{4}}{\mu\_{1} \varepsilon\_{2}}\right) D\_{\mathcal{B}} \Delta B \\ &+ \left(\frac{\varepsilon\_{3}}{\mu\_{2}} + \frac{\gamma \mathcal{Z}\_{4}}{\mu\_{2}}\right) \left(1 - \frac{V\_{4}}{V}\right) D\_{\mathcal{V}} \Delta V + \frac{\gamma}{\omega} \left(1 - \frac{\mathcal{Z}\_{4}}{\mathcal{Z}}\right) D\_{\mathcal{Z}} \Delta \mathcal{Z}. \end{split}$$

By using (5), the derivative of Ξ4(*t*) has the form

*d*Ξ<sup>4</sup> *dt* <sup>=</sup> <sup>−</sup> <sup>1</sup> Ψ (*U* − *U*4) 2 *<sup>U</sup>* <sup>d</sup>*<sup>x</sup>* <sup>+</sup> *<sup>η</sup>*2*U*4*V*<sup>4</sup> Ψ <sup>3</sup> <sup>−</sup> *<sup>U</sup>*<sup>4</sup> *<sup>U</sup>* <sup>−</sup> *U I<sup>V</sup>* <sup>4</sup> *V U*<sup>4</sup> *IVV*<sup>4</sup> <sup>−</sup> *<sup>I</sup>VV*<sup>4</sup> *IV* <sup>4</sup> *V* d*x* <sup>+</sup> 3<sup>4</sup> *μ*1<sup>2</sup> *λωη*1*μ*12<sup>5</sup> 34[*ω*1<sup>5</sup> + *η*2*μ*26] <sup>+</sup> <sup>1</sup> <sup>−</sup> <sup>2</sup> 3 − R1*V* Ψ *B* d*x* − *DUU*<sup>4</sup> Ψ *U*<sup>2</sup> *<sup>U</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* − *DIV I V* 4 Ψ *IV*<sup>2</sup> *<sup>I</sup>V*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DVV*<sup>4</sup> <sup>3</sup> *μ*2 <sup>+</sup> *<sup>γ</sup>Z*<sup>4</sup> *μ*2 Ψ *V*<sup>2</sup> *<sup>V</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *<sup>γ</sup>DZZ*<sup>4</sup> *ω* Ψ *Z*<sup>2</sup> *<sup>Z</sup>*<sup>2</sup> <sup>d</sup>*x*. It follows that *<sup>d</sup>*Ξ<sup>4</sup> *dt* <sup>≤</sup> 0 if *λωη*1*μ*12<sup>5</sup> 34[*ω*1<sup>5</sup> + *η*2*μ*26] <sup>+</sup> <sup>1</sup> <sup>≤</sup> <sup>2</sup> 3 <sup>+</sup> <sup>R</sup>1*V*. In addition, *<sup>d</sup>*Ξ<sup>4</sup> *dt* <sup>=</sup> <sup>0</sup> when (*U*, *L*, *IB*, *IV*, *B*, *V*, *Z*)=(*U*4, 0, 0, *I<sup>V</sup>* <sup>4</sup> , 0, *V*4, *Z*4). Hence, *χ* <sup>4</sup> = {*E*4} and *E*<sup>4</sup> is GS if

$$\mathcal{R}\_{1V} > 1 \text{ and } \frac{\lambda \omega \eta\_1 \mu\_1 \varepsilon\_2 \varepsilon\_5}{\varepsilon\_3 \varepsilon\_4 [\omega \varepsilon\_1 \varepsilon\_5 + \eta\_2 \mu\_2 \varepsilon\_6]} + 1 \le \frac{\varepsilon\_2}{\varepsilon\_3} + \mathcal{R}\_{1V} \text{ based on LIP [43].} \quad \Box$$

**Theorem 7.** *Suppose that λωη*2*μ*2<sup>4</sup> 25[*ω*1<sup>4</sup> + *η*1*μ*126] <sup>+</sup> <sup>1</sup> > <sup>3</sup> 2 <sup>+</sup> <sup>R</sup>1*B, λωη*1*μ*12<sup>5</sup> 34[*ω*1<sup>5</sup> + *η*2*μ*26] + <sup>1</sup> > <sup>2</sup> 3 + R1*V,* R0*V* R0*B* <sup>&</sup>gt; <sup>1</sup>*,* <sup>2</sup> <sup>&</sup>gt; 3*, and <sup>σ</sup>* <sup>&</sup>gt; <sup>1</sup>*. Then, the equilibrium E*<sup>5</sup> *is GS.*

**Proof.** We start with an LF

$$
\Xi\_{\mathbb{S}}(t) = \int\_{\Psi} \underline{\Xi}\_{\mathbb{S}}(\mathfrak{x}, t) \, \mathrm{d}x,\quad \text{where}
$$

Ξ0<sup>5</sup> =*U*<sup>5</sup> *U U*<sup>5</sup> <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>U</sup> U*<sup>5</sup> + *L*<sup>5</sup> *L L*5 <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>L</sup> L*5 + *I B* 5 *IB IB* 5 <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>I</sup><sup>B</sup> IB* 5 + *I V* 5 *IV IV* 5 <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>I</sup><sup>V</sup> IV* 5 + 1 *μ*1 <sup>+</sup> *<sup>γ</sup>Z*<sup>5</sup> *μ*1<sup>2</sup> *B*5 *B B*5 <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>B</sup> B*5 + <sup>3</sup> *μ*2 <sup>+</sup> *<sup>γ</sup>Z*<sup>5</sup> *μ*2 *V*5 *V V*5 <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>V</sup> V*5 <sup>+</sup> *<sup>γ</sup> <sup>ω</sup> <sup>Z</sup>*<sup>5</sup> *Z Z*5 <sup>−</sup> <sup>1</sup> <sup>−</sup> ln *<sup>Z</sup> Z*5 .

By computing the partial derivative, we obtain

*∂*Ξ<sup>5</sup> *<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>U</sup>*<sup>5</sup> *U DU*Δ*<sup>U</sup>* <sup>+</sup> *<sup>λ</sup>* <sup>−</sup> *<sup>η</sup>*1*UB* <sup>−</sup> *<sup>η</sup>*2*UV* <sup>−</sup> 1*<sup>U</sup>* + <sup>1</sup> <sup>−</sup> *<sup>L</sup>*<sup>5</sup> *L DL*Δ*<sup>L</sup>* <sup>+</sup> *<sup>η</sup>*1*UB* <sup>−</sup> *aL* + <sup>1</sup> <sup>−</sup> *<sup>I</sup><sup>B</sup>* 5 *IB DIB*Δ*I <sup>B</sup>* <sup>+</sup> *aL* <sup>−</sup> *<sup>γ</sup><sup>I</sup> BZ* <sup>−</sup> <sup>2</sup> *<sup>I</sup> B* + <sup>1</sup> <sup>−</sup> *<sup>I</sup><sup>V</sup>* 5 *IV DIV* Δ*I <sup>V</sup>* <sup>+</sup> *<sup>η</sup>*2*UV* <sup>−</sup> *<sup>γ</sup><sup>I</sup> <sup>V</sup> <sup>Z</sup>* <sup>−</sup> <sup>3</sup> *<sup>I</sup> V* + 1 *μ*1 <sup>+</sup> *<sup>γ</sup>Z*<sup>5</sup> *μ*1<sup>2</sup> <sup>1</sup> <sup>−</sup> *<sup>B</sup>*<sup>5</sup> *B DB*Δ*<sup>B</sup>* <sup>+</sup> *<sup>μ</sup>*1<sup>2</sup> *<sup>I</sup> <sup>B</sup>* <sup>−</sup> 4*<sup>B</sup>* + <sup>3</sup> *μ*2 <sup>+</sup> *<sup>γ</sup>Z*<sup>5</sup> *μ*2 <sup>1</sup> <sup>−</sup> *<sup>V</sup>*<sup>5</sup> *V DV*Δ*<sup>V</sup>* <sup>+</sup> *<sup>μ</sup>*<sup>2</sup> *<sup>I</sup> <sup>V</sup>* <sup>−</sup> 5*<sup>V</sup>* <sup>+</sup> *<sup>γ</sup> ω* <sup>1</sup> <sup>−</sup> *<sup>Z</sup>*<sup>5</sup> *Z DZ*Δ*<sup>Z</sup>* <sup>+</sup> *<sup>ω</sup><sup>I</sup> BZ* + *ωI <sup>V</sup> <sup>Z</sup>* <sup>−</sup> 6*<sup>Z</sup>* .

At equilibrium, the following conditions are satisfied:

$$\begin{cases} \lambda = \eta\_1 \mathcal{U}\_5 B\_5 + \eta\_2 \mathcal{U}\_5 V\_5 + \varepsilon\_1 \mathcal{U}\_5 \tau\_2 \\ \eta\_1 \mathcal{U}\_5 B\_5 = a \mathcal{L}\_5 \\ a \mathcal{L}\_5 = \gamma I\_5^B Z\_5 + \varepsilon\_2 I\_5^B \\ \eta\_2 \mathcal{U}\_5 V\_5 = \gamma I\_5^V Z\_5 + \varepsilon\_3 I\_5^V \\ \varepsilon\_2 I\_5^B = \frac{\varepsilon\_4}{\mu\_1} B \mathfrak{s}\_{\prime} \\ \varepsilon\_3 I\_5^V = \frac{\varepsilon\_3 \varepsilon\_5}{\mu\_2} V\_5 \\ \gamma I\_5^B Z \mathfrak{s} + \gamma I\_5^V Z \mathfrak{s} = \frac{\gamma \varepsilon\_6}{\omega} Z \mathfrak{s}. \end{cases}$$

By using the above conditions with (5), the time derivative of Ξ5(*t*) is written as

*d*Ξ<sup>5</sup> *dt* <sup>=</sup> <sup>−</sup> <sup>1</sup> Ψ (*U* − *U*5) 2 *<sup>U</sup>* <sup>d</sup>*<sup>x</sup>* <sup>+</sup> *<sup>η</sup>*1*U*5*B*<sup>5</sup> Ψ <sup>4</sup> <sup>−</sup> *<sup>U</sup>*<sup>5</sup> *<sup>U</sup>* <sup>−</sup> *UL*5*<sup>B</sup> U*5*LB*<sup>5</sup> <sup>−</sup> *LI<sup>B</sup>* 5 *<sup>L</sup>*<sup>5</sup> *<sup>I</sup><sup>B</sup>* <sup>−</sup> *<sup>I</sup>BB*<sup>5</sup> *IB* 5 *B* d*x* + *η*2*U*5*V*<sup>5</sup> Ψ <sup>3</sup> <sup>−</sup> *<sup>U</sup>*<sup>5</sup> *<sup>U</sup>* <sup>−</sup> *U I<sup>V</sup>* <sup>5</sup> *V U*<sup>5</sup> *IVV*<sup>5</sup> <sup>−</sup> *<sup>I</sup>VV*<sup>5</sup> *IV* <sup>5</sup> *V* d*x* − *DUU*<sup>5</sup> Ψ *U*<sup>2</sup> *<sup>U</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DLL*<sup>5</sup> Ψ *L*<sup>2</sup> *<sup>L</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* − *DIB I B* 5 Ψ *IB*<sup>2</sup> *<sup>I</sup>B*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DIV <sup>I</sup> V* 5 Ψ *IV*<sup>2</sup> *<sup>I</sup>V*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *DBB*<sup>5</sup> 1 *μ*1 <sup>+</sup> *<sup>γ</sup>Z*<sup>5</sup> *μ*1<sup>2</sup> Ψ *B*<sup>2</sup> *<sup>B</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* − *DVV*<sup>5</sup> <sup>3</sup> *μ*2 <sup>+</sup> *<sup>γ</sup>Z*<sup>5</sup> *μ*2 Ψ *V*<sup>2</sup> *<sup>V</sup>*<sup>2</sup> <sup>d</sup>*<sup>x</sup>* <sup>−</sup> *<sup>γ</sup>DZZ*<sup>5</sup> *ω* Ψ *Z*<sup>2</sup> *<sup>Z</sup>*<sup>2</sup> <sup>d</sup>*x*.

Thus, *<sup>d</sup>*Ξ<sup>5</sup> *dt* <sup>≤</sup> 0 and *<sup>d</sup>*Ξ<sup>5</sup> *dt* <sup>=</sup> 0 when (*U*, *<sup>L</sup>*, *<sup>I</sup>B*, *<sup>I</sup>V*, *<sup>B</sup>*, *<sup>V</sup>*, *<sup>Z</sup>*)=(*U*5, *<sup>L</sup>*5, *<sup>I</sup><sup>B</sup>* <sup>5</sup> , *<sup>I</sup><sup>V</sup>* <sup>5</sup> , *B*5, *V*5, *Z*5). This implies that *χ* <sup>5</sup> = {*E*5} and *E*<sup>5</sup> is GS when it exists in regard to LIP [43].

#### **5. Numerical Simulations**

In this part, we implement numerical simulations using MATLB PDE solver (pdepe) to validate the theoretical observations attained in the previous parts. This solver solves initial boundary value problems for systems of PDEs in one spatial variable *x* and time *t*. The domain of *x* is provided as Ψ = [0, 2] with step sizes Δ*x* = 0.02 and Δ*t* = 0.1. The ICs of system (1) are determined as the following:

*U*(*x*, 0) = 105(1 + 0.2 cos2(*πx*)), *L*(*x*, 0) = 104(1 + 0.2 cos2(*πx*)), *I <sup>B</sup>*(*x*, 0) = 103(1 + 0.2 cos2(*πx*)), *I <sup>V</sup>*(*x*, 0) = 103(1 + 0.2 cos2(*πx*)), *B*(*x*, 0) = 500(1 + 0.2 cos2(*πx*)), *V*(*x*, 0) = 500(1 + 0.2 cos2(*πx*)), *Z*(*x*, 0) = 0.1(1 + 0.2 cos2(*πx*)).

> To present the global stability of the equilibria of system (1), the results are divided into six cases. In each case, we change the values of *η*1, *η*2, and *ω* while keeping all other values as shown in Table 1. These cases are stated as follows:


#### *5.1. The Movement from the Monoinfection to the Coinfection State*

From the results above, we see that increasing the infection rate of EPCs by SARS-CoV-2, *η*2, forces the system to move from Mtb monoinfection state to SARS-CoV-2/Mtb coinfection state. In other words, *E*<sup>3</sup> loses its stability and *E*<sup>5</sup> becomes GS. Similarly, increasing the infection rate by Mtb, *η*1, pushes the system from SARS-CoV-2 monoinfection state to the coinfection state. In this case, *E*<sup>4</sup> loses its stability and *E*<sup>5</sup> becomes GS. Therefore, the values of these parameters need to be controlled as they have a powerful effect in converting the system from the monoinfection zone to the coinfection zone.

#### *5.2. The Impact of the Diffusion Coefficients On Coinfection*

To test the impact of the diffusion coefficients in model (1) on the behavior of the solutions, we change the values of the coefficients considered in case (vi) to *DU* = *DL* = *DIB* = *DIV* <sup>=</sup> *DB* <sup>=</sup> *DV* <sup>=</sup> *DZ* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−5. We observe from Figure 7 that the effect of this change appears at the initial times, while the final solutions are not affected. Thus, the diffusion coefficients do not affect the robustness of the global stability of the solutions. Therefore, the impact of these coefficients should be monitored at the beginning of coinfection as it affects the distribution of particles in space.


**Table 1.** Parameters' values of system (1).

(**g**) CTLs **Figure 1.** The numerical results of system (1) for *<sup>η</sup>*<sup>1</sup> <sup>=</sup> 2.5 <sup>×</sup> <sup>10</sup><sup>−</sup>9, *<sup>η</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>11, and *<sup>ω</sup>* <sup>=</sup> <sup>8</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>3. The uninfected equilibrium *E*<sup>0</sup> = <sup>4</sup> <sup>×</sup> <sup>10</sup>5, 0, 0, 0, 0, 0, 0 is GS.

(**g**) CTLs **Figure 3.** The numerical results of system (1) for *<sup>η</sup>*<sup>1</sup> <sup>=</sup> 2.5 <sup>×</sup> <sup>10</sup>−9, *<sup>η</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−9, and *<sup>ω</sup>* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−8. The equilibrium *E*<sup>2</sup> = 8571.43, 0, 0, 391,429, 0, 4.56667 <sup>×</sup> <sup>10</sup>8, 0 is GS.

**Figure 4.** The numerical results of system (1) for *<sup>η</sup>*<sup>1</sup> <sup>=</sup> 2.5 <sup>×</sup> <sup>10</sup><sup>−</sup>7, *<sup>η</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>11, and *<sup>ω</sup>* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>4. The equilibrium *E*<sup>3</sup> = (117,514, 7062.15, 1000, 0, 96,153.8, 0, 4.64972) is GS.

**Figure 6.** The numerical results of system (1) for *<sup>η</sup>*<sup>1</sup> <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>−7, *<sup>η</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−9, and *<sup>ω</sup>* <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>−6. The equilibrium *E*<sup>5</sup> = 27,125.6, 5712.6, 4380.44, 45,619.6, 421,196, 5.322 <sup>×</sup> 107, 0.043 is GS.

**Figure 7.** The impact of changing the diffusion coefficients in case (vi) to 1 <sup>×</sup> <sup>10</sup>−5. The initial distributions of the solutions are affected, while the global stability is not affected.

#### **6. Conclusions and Future Works**

There is an emerging evidence that the COVID-19 patients who have Mtb are more likely to develop acute disease and die [2,3,8]. Therefore, understanding Mtb/SARS-CoV-2 coinfection is critical to treat this group of patients. Here, we introduced a reaction–diffusion within-host Mtb/SARS-CoV-2 model. It counts the connections between uninfected EPCs, latently Mtb-infected EPCs, productively Mtb-infected EPCs, SARS-CoV-2-infected EPCs, Mtb particles, SARS-CoV-2 virions, and CTLs. It owns six equilibrium points as the following:


$$\begin{array}{ll}\text{(vi)} & \text{The } & \text{Mtb/SAR-CovV-2} \text{ coinfection equilibrium } E\_5 \text{ exists and it is CS if} \\ & \frac{\omega \lambda \eta\_2 \mu\_2 \varepsilon\_4}{\varepsilon\_2 \varepsilon\_5 [\omega \epsilon\_1 \epsilon\_4 + \eta\_1 \mu\_1 \varepsilon\_2 \varepsilon\_6]} + 1 > \frac{\varepsilon\_3}{\varepsilon\_2} + \mathcal{R}\_{1B} \frac{\omega \lambda \eta\_1 \mu\_1 \varepsilon\_2 \varepsilon\_5}{\varepsilon\_3 \varepsilon\_4 [\omega \epsilon\_1 \epsilon\_5 + \eta\_2 \mu\_2 \varepsilon\_6]} + 1 > \frac{\varepsilon\_2}{\varepsilon\_3} + \mathcal{R}\_{1V}, \\ & \frac{\mathcal{R}\_{0V}}{\mathcal{R}\_{0B}} > 1, \; \varepsilon\_2 > \varepsilon\_3, \text{ and } \sigma > 1. \text{ Here, the patient with a single infection becomes } \\ & \text{infected with both SAR-CovV-2 and Mtb.} \end{array}$$

We found that the numerical computations are quite congruous with the theoretical contributions. The equilibrium points reflect three states: the healthy state, the monoinfection state, and the coinfection state. The threshold parameters defined in Proposition 1 determine the locomotion between these states. Thus, the values of parameters in model (1) should be selected with caution. In addition, the global stability of the solutions of model (1) is robust against the values of the diffusion coefficients. However, the initial distributions of particles are affected by the selection of these values. Thus, it should be monitored as it may affect the initial status of the coinfected patients. In fact, Mtb/SARS-CoV-2 coinfection is a disease that needs to be further investigated and requires more awareness in high-TB burden regions such as India, Indonesia, and China [2]. Understanding the dynamics of coinfection will help develop new treatments, find better ways to treat coinfected patients, or recommend preventive measures for coinfected patients. The main limitation of this work is that we did not acquire real data to estimate the values of parameters in system (1). We gathered the values from SARS-CoV-2 monoinfection models or Mtb monoinfection models. Furthermore, we proved the boundedness only for the case when *DU* = *DL* = *DIB* = *DIV* = *DZ*. In addition, we assumed that CTLs kill infected cells at the same rate constant. Therefore, this work could be polished by (i) utilizing real data to obtain an estimation of the values of parameters in system (1) when the data on coinfection become available, (ii) proving the boundedness for different diffusion coefficients, (iii) analyzing the model with different killing rates of CTLs, (iv) counting the time delays inherent in the latent stage or other responses, (v) adding the role of antibodies in eliminating SARS-CoV-2 or Mtb particles, (vi) using fractional derivatives to study model (1) [45,46], (vii) performing a sensitivity analysis for the threshold parameters to identify the most sensitive parameters in the model [47], (viii) considering mutations that can generate more aggressive variants of

SARS-CoV-2 and their effect on coinfection dynamics [48], and (ix) developing a multiscale model to connect within-host dynamics with between-hosts dynamics and gain a better comprehension of the coinfection mechanism.

**Author Contributions:** Conceptualization, A.A., A.F. and H.A.G.; Methodology, A.A., A.D.A.A., A.F. and H.A.G.; Formal analysis, A.D.A.A.; Investigation, A.A. and A.F.; Writing—original draft, A.D.A.A.; Writing—review & editing, A.D.A.A. and H.A.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project was funded by the Deanship Scientific Research (DSR), King Abdulaziz University, Jeddah, under the Grant No. (IFPIP:699-130-1443).

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** This project was funded by the Deanship Scientific Research (DSR), King Abdulaziz University, Jeddah, under the Grant No. (IFPIP:699-130-1443). The authors acknowledge with thanks DSR for technical and financial support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
