Global Dynamics of an HTLV-I and SARS-CoV-2 Co-Infection Model with Diffusion

Abstract

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a novel respiratory virus that causes coronavirus disease 2019 (COVID-19). Symptoms of COVID-19 range from mild to severe illness. It was observed that disease progression in COVID-19 patients depends on their immune response, especially in elderly patients whose immune system suppression may put them at increased risk of infection. Human T-cell lymphotropic virus type-I (HTLV-I) attacks the CD4⁺ T cells (T cells) of the immune system and leads to immune dysfunction. Co-infection with HTLV-I and SARS-CoV-2 has been reported in recent studies. Modeling HTLV-I and SARS-CoV-2 co-infection can be a helpful tool to understand the in-host co-dynamics of these viruses. The aim of this study was to construct a model that characterizes the in-host dynamics of HTLV-I and SARS-CoV-2 co-infection. By considering the mobility of the viruses and cells, the model is represented by a system of partial differential equations (PDEs). The system contains two independent variables, time t and position x, and seven dependent variables for representing the densities of healthy epithelial cells (ECs), latent SARS-CoV-2-infected ECs, active SARS-CoV-2-infected ECs, SARS-CoV-2, healthy T cells, latent HTLV-I-infected T cells and active HTLV-I-infected T cells. We first studied the fundamental properties of the solutions of the system, then deduced all steady states and proved their global properties. We examined the global stability of the steady states by constructing appropriate Lyapunov functions. The analytical results were illustrated by performing numerical simulations. We discussed the effect of HTLV-I infection on COVID-19 progression. The results suggest that patients with HTLV-I have a weakened immune response; consequently, their risk of COVID-19 infection may be increased.

1. Introduction

In recent years, several types of viruses that target the human body have spread, some of which cause serious diseases and may even lead to death. The spread of these diseases has a significant effect on the global health burden. Examples of such viruses are human immunodeficiency virus class-1 (HIV-1), human T-cell lymphotropic virus class-I (HTLV-I), hepatitis B (and C) viruses, influenza virus, Ebola virus, chikungunya virus, dengue virus, Zika virus, and Middle East Respiratory Syndrome coronavirus. On December 2019, a novel respiratory virus appeared, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). SARS-CoV-2 is the causative agent of coronavirus disease 2019 (COVID-19). This pandemic swept around the world and caused the death of millions of people [1]. Since the beginning of this pandemic, the efforts of researchers and scientists from several disciplines have been directed to find different ways to confront this pandemic, such as the synthesis of vaccines and antiviral drugs [2]. On the other hand, governments implemented several measures to reduce the spread of the pandemic.

SARS-CoV-2 is an RNA virus from the Coronaviridae family. This virus is transmitted to individuals when exposed to the fluids of the respiratory system that contain viruses. After entering the human body, the virus targets epithelial cells (ECs) in the respiratory tract [3]. the symptoms of COVID-19 include rhinitis, headache, myalgia, cough, dyspnea, fever, and sore throat. The body’s immunity is an important element in confronting emerging viruses. The progression of the infection in COVID-19 patients depends on the response of the immune system to the virus, especially in elderly patients, whose immunosuppression may predispose them to an increased risk of infection [4].

HTLV-I is a fetal virus that attacks humans. It is a blood-borne virus and sexually transmitted infection [5]. HTLV-I retrovirus attacks CD4${}^{+}$T cells (T cells,) which are one the most effective components of the adaptive immune response. This leads to immune dysfunction in HTLV-I patients [6]. Infection by HTLV-I can cause a variety of symptoms that lead to fatal diseases; the most well-known diseases associated with HTLV-I infection are HTLV-I-associated myelopathy/tropical spastic paraparesis (HAM/TSP) and adult T-cell leukemia (ATL). Asquith and Bangham [7] presented an experimental study explaining how HTLV-I infection occurs and persists. T cells infected with HTLV-I are, in fact, divided into two distinct compartments, latent and active. Latent HTLV-I-infected cells involve the viral DNA but are incapable of generating it, so until they are activated, they cannot infect susceptible T cells. Active HTLV-I-infected T cells may then convert to ATL cells.

B cells, cytotoxic T lymphocytes (CTLs), and T cells are central elements of adaptive immunity. B cells and CTLs are both detected in the blood of COVID-19 patients about one week after the onset of symptoms [8]. B cells produce antibodies to neutralize SARS-CoV-2 [9]. CTL is an essential component controlling HTLV-I related diseases by killing viral-infected cells [7,10]. T cells are essential for the processing of both CTLs and B cells [8]. Sajjadi et al. [11] provided an overview of the pathogenesis of both HTLV-I and SARS-CoV-2 infections and on their similarities in stimulating the immune response.

Co-infections of COVID-19 and other diseases have become the interest of many scientists and researchers. Zhu et al. [12] reported that 94.2% of individuals infected with COVID-19 were also co-infected with many other types of micro-organisms, such as viruses, bacteria, and fungi. Co-infections with HTLV-I and SARS-CoV-2 were reported [6,13,14,15].

1.1. Mathematical Models of In-Host HTLV-I and SARS-CoV-2 Single Infections

Mathematical modeling of in-host viral dynamics has attracted the attention of researchers. Mathematical models and their analyses are useful tools for better understanding in-host virus dynamics, enabling the prediction of viral progression and the development of anti-viral treatment strategies.

1.1.1. SARS-CoV-2 Single-Infection Model

Hernandez-Vargas and Velasco-Hernandez [16] presented a SARS-CoV-2 single-infection model with limited target cells as:

$$\begin{array}{cc}\hfill \frac{dX\left(t\right)}{dt}& =-\stackrel{\mathrm{SARS}-\mathrm{CoV}-2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\sigma V\left(t\right)X\left(t\right)}},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{dN\left(t\right)}{dt}& =\stackrel{\mathrm{SARS}-\mathrm{CoV}-2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\sigma V\left(t\right)X\left(t\right)}}-\stackrel{\mathrm{activation}\mathrm{of}\mathrm{latent}\mathrm{infected}\mathrm{ECs}}{\overbrace{\varrho N\left(t\right)}},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{dY\left(t\right)}{dt}& =\stackrel{\mathrm{activation}\mathrm{of}\mathrm{latent}\mathrm{infected}\mathrm{ECs}}{\overbrace{\varrho N\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{Y}Y\left(t\right)}},\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{dV\left(t\right)}{dt}& =\stackrel{\mathrm{production}\mathrm{of}\mathrm{SARS}-\mathrm{CoV}-2}{\overbrace{\rho Y\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{V}V\left(t\right)}},\hfill \end{array}$$

(4)

where t is the time; $X\left(t\right)$, $N\left(t\right)$, and $Y\left(t\right)$ are the densities of healthy, latent, and active ECs, respectively; and $V\left(t\right)$ is the concentration of SARS-CoV-2 particles. In [17], the death and reproduction of healthy ECs were considered. Formulated models [16,17] were modified and extended by taking into account the influence of the immune response [18,19,20,21,22,23,24], drug therapies [25,26,27], and time delay [28]. The stability of in-host COVID-19 single-infection models was investigated [22,23,24,28,29].

1.1.2. HTLV-I Single-Infection Models

The latent HTLV-I single-infection model (ignoring the ATL) is given as [30,31]:

$$\begin{array}{cc}\hfill \frac{dU\left(t\right)}{dt}& =\stackrel{\mathrm{generation}\mathrm{of}\mathrm{healthy}\mathrm{T}\mathrm{cells}}{\overbrace{\varkappa }}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{U}U\left(t\right)}}-\stackrel{\mathrm{HTLV}-\mathrm{I}\mathrm{infectious}\mathrm{transmission}}{\overbrace{\tau A\left(t\right)U\left(t\right)}},\hfill \\ \hfill \frac{dL\left(t\right)}{dt}& =\stackrel{\mathrm{HTLV}-\mathrm{I}\mathrm{infectious}\mathrm{transmission}}{\overbrace{\tau A\left(t\right)U\left(t\right)}}-\stackrel{\mathrm{activation}\mathrm{of}\mathrm{latent}\mathrm{infected}\mathrm{T}\mathrm{cells}}{\overbrace{\lambda L\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{L}L\left(t\right)}},\hfill \\ \hfill \frac{dA\left(t\right)}{dt}& =\stackrel{\mathrm{activation}\mathrm{of}\mathrm{latent}\mathrm{infected}\mathrm{T}\mathrm{cells}}{\overbrace{\lambda L\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_A^{*}A\left(t\right)}},\hfill \end{array}$$

where $U\left(t\right)$, $L\left(t\right)$, and $A\left(t\right)$ are the densities of healthy, latent, and active T cells, respectively. Different biological processes were included in the modeling of the infection of HTLV-I by incorporating (i) CTL immunity [32,33,34,35,36], (ii) active T cells mitosis [37,38,39,40,41], and (iii) time delay [33,42,43,44,45].

1.1.3. HTLV-I and SARS-CoV-2 Co-Infection Model

From 2020 to 2022, some epidemiological models were developed to describe the between-host transmission of COVID-19 co-infection with other diseases such as bacteria/COVID-19 [46], tuberculosis/COVID-19 [47], Dengue/COVID-19 [48], Dengue/HIV/COVID-19 [49], ZIKV/COVID-19 [50], HIV/COVID-19 [51], and influenza/COVID-19 [52]. On the other hand, virological models that describe the in-host dynamics of SARS-CoV-2 with co-infection with other micro-organisms have been investigated in recent studies: influenza A virus/SARS-CoV-2 [53,54], malaria/SARS-CoV-2 [55], HIV/SARS-CoV-2 [56], and bacteria/ SARS-CoV-2 [57]. To formulate an HTLV-I and SARS-CoV-2 co-infection model, Elaiw et al. [58] considered the following factors:

Factor 1

The model contains seven interplaying populations: healthy ECs (X), latent SARS-CoV-2-infected ECs (N), active SARS-CoV-2-infected ECs (Y), SARS-CoV-2 particles (V), healthy T cells (U), latent HTLV-I-infected T cells (L), and active HTLV-I-infected T cells (A).

Factor 2

The healthy T cells and ECs are the targets of HTLV-I and SARS-CoV-2, respectively [16,30].

Factor 3

The healthy T cells help the CTLs to destroy the active SARS-CoV-2-infected ECs at rate $\xi YU$, which expand at rate $\vartheta YU$ [56].

Factor 4

The active HTLV-I-infected cells are replicated with rate ${\gamma }^{*}A$, some of which $\phi {\gamma }^{*}A$ are converted into latent infected, and the remainder $\left(1-\phi \right){\gamma }^{*}A$ stay active, where, $\phi \in (0,1)$ [40,59].

Based on these factors, the HTLV-I and SARS-CoV-2 co-infection model was described by a system of ODEs as [58]:

$$\begin{array}{cc}\hfill \frac{dX\left(t\right)}{dt}& =\stackrel{\mathrm{generation}\mathrm{of}\mathrm{healthy}\mathrm{ECs}}{\overbrace{\upsilon }}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{X}X\left(t\right)}}-\stackrel{\mathrm{SARS}-\mathrm{CoV}-2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\sigma V\left(t\right)X\left(t\right)}},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \frac{dN\left(t\right)}{dt}& =\stackrel{\mathrm{SARS}-\mathrm{CoV}-2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\sigma V\left(t\right)X\left(t\right)}}-\stackrel{\mathrm{activation}\mathrm{of}\mathrm{latent}\mathrm{infected}\mathrm{ECs}}{\overbrace{\varrho N\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{N}N\left(t\right)}},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{dY\left(t\right)}{dt}& =\stackrel{\mathrm{activation}\mathrm{of}\mathrm{latent}\mathrm{infected}\mathrm{ECs}}{\overbrace{\varrho N\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{Y}Y\left(t\right)}}-\stackrel{\mathrm{killing}\mathrm{of}\mathrm{active}\mathrm{infected}\mathrm{ECs}}{\overbrace{\xi Y\left(t\right)U\left(t\right)}},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{dV\left(t\right)}{dt}& =\stackrel{\mathrm{production}\mathrm{of}\mathrm{SARS}-\mathrm{CoV}-2}{\overbrace{\rho Y\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{V}V\left(t\right)}},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \frac{dU\left(t\right)}{dt}& =\stackrel{\mathrm{generation}\mathrm{of}\mathrm{healthy}\mathrm{T}\mathrm{cells}}{\overbrace{\varkappa }}+\stackrel{\mathrm{proliferation}\mathrm{of}\mathrm{healthy}\mathrm{T}\mathrm{cells}}{\overbrace{\vartheta Y\left(t\right)U\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{U}U\left(t\right)}}-\stackrel{\mathrm{HTLV}-\mathrm{I}\mathrm{infectious}\mathrm{transmission}}{\overbrace{\tau U\left(t\right)A\left(t\right)}},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{dL\left(t\right)}{dt}& =\stackrel{\mathrm{HTLV}-\mathrm{I}\mathrm{infectious}\mathrm{transmission}}{\overbrace{\tau U\left(t\right)A\left(t\right)}}+\stackrel{\mathrm{mitotic}\mathrm{transmission}}{\overbrace{\phi {\gamma }^{*}A\left(t\right)}}-\stackrel{\mathrm{activation}\mathrm{of}\mathrm{latent}\mathrm{infected}\mathrm{T}\mathrm{cells}}{\overbrace{\lambda L\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_{L}L\left(t\right)}},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \frac{dA\left(t\right)}{dt}& =\stackrel{\mathrm{activation}\mathrm{of}\mathrm{latent}\mathrm{infected}\mathrm{T}\mathrm{cells}}{\overbrace{\lambda L\left(t\right)}}+\stackrel{\mathrm{proliferation}\mathrm{of}\mathrm{active}\mathrm{infected}\mathrm{T}\mathrm{cells}}{\overbrace{\left(1-\phi \right){\gamma }^{*}A\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\beta }_A^{*}A\left(t\right)}}.\hfill \end{array}$$

(11)

The above models are based on the assumption that the T cells and viruses are homogeneously distributed in the body. It was reported [60,61,62] that T cells move from high-concentration to low-concentration places. The movement of epithelial cells was reported [63]. Therefore, spatial structure is important when describing the dynamics of SARS-CoV-2 and HTLV-I co-infection in a host. Wang and Ma [64] formulated a diffusive HTLV-I single-infection model with time delay, mitosis, and CTL response. SARS-CoV-2 infection dynamics models with reaction–diffusion were studied [65,66].

The aim of this study was to construct and analyze a mathematical model for the in-host dynamics of HTLV-I and SARS-CoV-2 co-infection with diffusion. The co-infection model is described by a set of PDEs. We established the fundamental properties of the system’s solutions, computed all steady states, and proved their global asymptotic stability. Furthermore, we performed numerical simulations to support the theoretical findings. We discussed the effect of HTLV-I infection on the progression of COVID-19 infection.

2. Model Formulation

We developed an HTLV-I and SARS-CoV-2 co-infection model by considering Factors 1–4; in addition, we assumed that cells and viruses are mobile in the environment. The model is then given by a system of PDEs as follows:

$$\begin{array}{cc}\hfill \frac{\partial X\left(x,t\right)}{\partial t}& ={\delta }_X\Delta X\left(x,t\right)+\upsilon -{\beta }_{X}X\left(x,t\right)-\sigma X\left(x,t\right)V\left(x,t\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \frac{\partial N\left(x,t\right)}{\partial t}& ={\delta }_N\Delta N\left(x,t\right)+\sigma X\left(x,t\right)V\left(x,t\right)-(\varrho +{\beta }_N)N\left(x,t\right),\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \frac{\partial Y\left(x,t\right)}{\partial t}& ={\delta }_Y\Delta Y\left(x,t\right)+\varrho N\left(x,t\right)-{\beta }_{Y}Y\left(x,t\right)-\xi Y\left(x,t\right)U\left(x,t\right),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{\partial V\left(x,t\right)}{\partial t}& ={\delta }_V\Delta V\left(x,t\right)+\rho Y\left(x,t\right)-{\beta }_{V}V\left(x,t\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \frac{\partial U\left(x,t\right)}{\partial t}& ={\delta }_U\Delta U\left(x,t\right)+\varkappa +\vartheta Y\left(x,t\right)U\left(x,t\right)-{\beta }_{U}U\left(x,t\right)-\tau U\left(x,t\right)A\left(x,t\right),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \frac{\partial L\left(x,t\right)}{\partial t}& ={\delta }_L\Delta L\left(x,t\right)+\tau U\left(x,t\right)A\left(x,t\right)+\phi {\gamma }^{*}A\left(x,t\right)-\left(\lambda +{\beta }_L\right)L\left(x,t\right),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \frac{\partial A\left(x,t\right)}{\partial t}& ={\delta }_A\Delta A\left(x,t\right)+\lambda L\left(x,t\right)+\left(1-\phi \right){\gamma }^{*}A\left(x,t\right)-{\beta }_A^{*}A\left(x,t\right),\hfill \end{array}$$

(18)

where $x\in \mathrm{\Omega },$$t>0.$ The spatial domain $\mathrm{\Omega }\subset {\mathbb{R}}^m,m\ge 1$ is connected and bounded by a smooth boundary $\partial \mathrm{\Omega }$ and $x=\left(x_1,x_2,\dots ,x_m\right)$. $\Delta =\frac{{\partial }^2}{\partial x^2}$ is the Laplacian operator, while ${\delta }_{\mathrm{\Lambda }}$ is the diffusion coefficient corresponding to compartment $\mathrm{\Lambda }$ of the system. All the parameters given in model (12)–(18) are positive.

The initial conditions are given by

$$\begin{array}{cc}\hfill X\left(x,0\right)& ={\mathrm{\Gamma }}_1\left(x\right),N\left(x,0\right)={\mathrm{\Gamma }}_2\left(x\right),Y\left(x,0\right)={\mathrm{\Gamma }}_3\left(x\right),V\left(x,0\right)={\mathrm{\Gamma }}_4\left(x\right)\hfill \\ \hfill U\left(x,0\right)& ={\mathrm{\Gamma }}_5\left(x\right),L\left(x,0\right)={\mathrm{\Gamma }}_6\left(x\right),A\left(x,0\right)={\mathrm{\Gamma }}_7\left(x\right),x\in \mathrm{\Omega }.\hfill \end{array}$$

(19)

Functions ${\mathrm{\Gamma }}_i\left(x\right)\ge 0,$$i=1,2,\dots ,7$ are continuous. Furthermore, the following homogeneous Neumann boundary conditions are considered:

$$\frac{\partial X}{\partial \overrightarrow{\eta }}=\frac{\partial N}{\partial \overrightarrow{\eta }}=\frac{\partial Y}{\partial \overrightarrow{\eta }}=\frac{\partial V}{\partial \overrightarrow{\eta }}=\frac{\partial U}{\partial \overrightarrow{\eta }}=\frac{\partial L}{\partial \overrightarrow{\eta }}=\frac{\partial A}{\partial \overrightarrow{\eta }}=0,t>0,x\in \partial \mathrm{\Omega },$$

(20)

where $\frac{\partial }{\partial \overrightarrow{\eta }}$ is the outward normal derivative on the boundary $\partial \mathrm{\Omega }$.

In [40,59], it was assumed that ${\gamma }^{*}<min\{{\beta }_U,{\beta }_L,{\beta }_A^{*}\}$. Because ${\gamma }^{*}<{\beta }_A^{*}$ and $0<\phi <1$, $(1-\phi ){\gamma }^{*}<{\beta }_A^{*}$. Let us assume that ${\beta }_A={\beta }_A^{*}-(1-\phi ){\gamma }^{*}>0$ and $\gamma =\phi {\gamma }^{*}$. We have ${\beta }_A-\gamma ={\beta }_A^{*}-{\gamma }^{*}>0$. Then, system (12)–(18) becomes:

$$\begin{array}{cc}\hfill \frac{\partial X\left(x,t\right)}{\partial t}& ={\delta }_X\Delta X\left(x,t\right)+\upsilon -{\beta }_{X}X\left(x,t\right)-\sigma V\left(x,t\right)X\left(x,t\right),\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \frac{\partial N\left(x,t\right)}{\partial t}& ={\delta }_N\Delta N\left(x,t\right)+\sigma V\left(x,t\right)X\left(x,t\right)-(\varrho +{\beta }_N)N\left(x,t\right),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \frac{\partial Y\left(x,t\right)}{\partial t}& ={\delta }_Y\Delta Y\left(x,t\right)+\varrho N\left(x,t\right)-{\beta }_{Y}Y\left(x,t\right)-\xi Y\left(x,t\right)U\left(x,t\right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \frac{\partial V\left(x,t\right)}{\partial t}& ={\delta }_V\Delta V\left(x,t\right)+\rho Y\left(x,t\right)-{\beta }_{V}V\left(x,t\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \frac{\partial U\left(x,t\right)}{\partial t}& ={\delta }_U\Delta U\left(x,t\right)+\varkappa +\vartheta Y\left(x,t\right)U\left(x,t\right)-{\beta }_{U}U\left(x,t\right)-\tau A\left(x,t\right)U\left(x,t\right),\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \frac{\partial L\left(x,t\right)}{\partial t}& ={\delta }_L\Delta L\left(x,t\right)+\tau A\left(x,t\right)U\left(x,t\right)+\gamma A\left(x,t\right)-\left(\lambda +{\beta }_L\right)L\left(x,t\right),\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \frac{\partial A\left(x,t\right)}{\partial t}& ={\delta }_A\Delta A\left(x,t\right)+\lambda L\left(x,t\right)-{\beta }_{A}A\left(x,t\right).\hfill \end{array}$$

(27)

3. Properties of Solutions

Lemma 1.

Assume that ${\delta }_X={\delta }_N={\delta }_Y={\delta }_V={\delta }_U={\delta }_L={\delta }_A=\tilde{\delta }$. Then, system (21)–(27) has a unique, non-negative, and bounded solution defined on $\overline{\mathrm{\Omega }}\times [0,+\infty )$ for any given initial datum satisfying (19).

Proof. 

Let $\mathcal{F}=BUC$$\left(\overline{\mathrm{\Omega }},{\mathbb{R}}^7\right)$ be the set of all bounded and uniformly continuous functions from $\overline{\mathrm{\Omega }}$ to ${\mathbb{R}}^7,$ and let ${\mathcal{F}}_{+}=BUC$$\left(\overline{\mathrm{\Omega }},{\mathbb{R}}_{+}^7\right)\subset \mathcal{F}.$ Then, the positive cone ${\mathcal{F}}_{+}$ induces a partial order on $\mathcal{F}$. Suppose that the norm is defined by ${∥\varphi ∥}_{\mathcal{F}}=\underset{x\in \overline{\mathrm{\Omega }}}{sup}\left|\varphi \left(x\right)\right|,$ where $\left|.\right|$ is the Euclidean norm on ${\mathbb{R}}^7$. This implies that the space $(\mathcal{F},{∥.∥}_{\mathcal{F}})$ is a Banach lattice [67,68].

For any initial data $\mathrm{\Gamma }={({\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_2,{\mathrm{\Gamma }}_3,{\mathrm{\Gamma }}_4,{\mathrm{\Gamma }}_5,{\mathrm{\Gamma }}_6,{\mathrm{\Gamma }}_7)}^{\prime }\in$${\mathcal{F}}_{+}$, we define $\mathrm{\Theta }=\left({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,{\mathrm{\Theta }}_3,{\mathrm{\Theta }}_4,\right.$${\left.{\mathrm{\Theta }}_5,{\mathrm{\Theta }}_6,{\mathrm{\Theta }}_7\right)}^{\prime }:{\mathcal{F}}_{+}\to \mathcal{F}$ by

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_1\left(\mathrm{\Gamma }\right)\left(x\right)& =\upsilon -{\beta }_X{\mathrm{\Gamma }}_1\left(x\right)-\sigma {\mathrm{\Gamma }}_4\left(x\right){\mathrm{\Gamma }}_1\left(x\right),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_2\left(\mathrm{\Gamma }\right)\left(x\right)& =\sigma {\mathrm{\Gamma }}_4\left(x\right){\mathrm{\Gamma }}_1\left(x\right)-(\varrho +{\beta }_N){\mathrm{\Gamma }}_2\left(x\right),\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_3\left(\mathrm{\Gamma }\right)\left(x\right)& =\varrho {\mathrm{\Gamma }}_2\left(x\right)-{\beta }_Y{\mathrm{\Gamma }}_3\left(x\right)-\xi {\mathrm{\Gamma }}_3\left(x\right){\mathrm{\Gamma }}_5\left(x\right),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_4\left(\mathrm{\Gamma }\right)\left(x\right)& =\rho {\mathrm{\Gamma }}_3\left(x\right)-{\beta }_V{\mathrm{\Gamma }}_4\left(x\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_5\left(\mathrm{\Gamma }\right)\left(x\right)& =\varkappa +\vartheta {\mathrm{\Gamma }}_3\left(x\right){\mathrm{\Gamma }}_5\left(x\right)-{\beta }_U{\mathrm{\Gamma }}_5\left(x\right)-\tau {\mathrm{\Gamma }}_7\left(x\right){\mathrm{\Gamma }}_5\left(x\right),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_6\left(\mathrm{\Gamma }\right)\left(x\right)& =\tau {\mathrm{\Gamma }}_7\left(x\right){\mathrm{\Gamma }}_5\left(x\right)+\gamma {\mathrm{\Gamma }}_7\left(x\right)-\left(\lambda +{\beta }_L\right){\mathrm{\Gamma }}_6\left(x\right),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_7\left(\mathrm{\Gamma }\right)\left(x\right)& =\lambda {\mathrm{\Gamma }}_6\left(x\right)-{\beta }_A{\mathrm{\Gamma }}_7\left(x\right).\hfill \end{array}$$

(34)

Clearly, $\mathrm{\Theta }$ is locally Lipschitz on ${\mathcal{F}}_{+}$. We can rewrite system (28)–(34) with initial and boundary conditions (19) and (20), respectively, as the following abstract functional DE:

$$\left\{\begin{array}{c}\frac{d\mathrm{\Lambda }}{dt}=\mathrm{Ð}\mathrm{\Lambda }+\mathrm{\Theta }\left(\mathrm{\Lambda }\right),t>0,\\ \mathrm{\Lambda }\left(0\right)=\mathrm{\Gamma }\in {\mathcal{F}}_{+},\end{array}\right.$$

where $\mathrm{\Lambda }={\left(X,N,Y,V,U,L,A\right)}^{\prime }$ and Ð$\mathrm{\Lambda }=\left({\delta }_X\Delta X,{\delta }_N\Delta N,{\delta }_Y\Delta Y,{\delta }_V\Delta V,{\delta }_U\Delta U,{\delta }_L\Delta L,\right.$${\left.{\delta }_A\Delta A\right)}^{\prime }.$ It is possible to show that $\underset{h\to 0^{+}}{lim}\frac{1}{h}$dist$\left(\mathrm{\Gamma }\left(0\right)+h\mathrm{\Theta }\left(\mathrm{\Gamma }\right),{\mathcal{F}}_{+}\right)=0,$ for all $\mathrm{\Gamma }\in {\mathcal{F}}_{+}.$ It follows from [67,68,69] that, for any $\mathrm{\Gamma }\in {\mathcal{F}}_{+}$, system (21)–(27) with initial and boundary conditions (19)–(20), has a unique, non-negative, and mild solution $\left(X\right(x,t),N(x,t),Y(x,t),V(x,t),$$U(x,t),L(x,t),A(x,t))$ defined on $\overline{\mathrm{\Omega }}\times [0,{\chi }_m)$, where $[0,{\chi }_m)$ is the maximal existence time interval.

To prove the boundedness of all state variables, we define

$$\mathrm{\Psi }(x,t)=X(x,t)+N(x,t)+Y(x,t)+\frac{{\beta }_Y}{2\rho }V(x,t)+\frac{\xi }{\vartheta }\left(U(x,t)+L(x,t)+A(x,t)\right).$$

Because ${\delta }_X={\delta }_N={\delta }_Y={\delta }_V={\delta }_U={\delta }_L={\delta }_A=\tilde{\delta }$, using Equations (21)–(27), we obtain

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Psi }(x,t)}{\partial t}-\tilde{\delta }\Delta \mathrm{\Psi }(x,t)& =\upsilon +\frac{\xi }{\vartheta }\varkappa -{\beta }_{X}X(x,t)-{\beta }_{N}N(x,t)-\frac{{\beta }_Y}{2}Y(x,t)-\frac{{\beta }_Y{\beta }_V}{2\rho }V(x,t)\hfill \\ \hfill & -\frac{\xi {\beta }_U}{\vartheta }U(x,t)-\frac{\xi {\beta }_L}{\vartheta }L(x,t)-\frac{\xi \left({\beta }_A-\gamma \right)}{\vartheta }A(x,t).\hfill \end{array}$$

We have ${\beta }_A-\gamma ={\beta }_A^{*}-{\gamma }^{*}>0$. Hence,

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Psi }(x,t)}{\partial t}-\tilde{\delta }\Delta \mathrm{\Psi }(x,t)& =\upsilon +\frac{\xi }{\vartheta }\varkappa -{\beta }_{X}X(x,t)-{\beta }_{N}N(x,t)-\frac{{\beta }_Y}{2}Y(x,t)-\frac{{\beta }_Y{\beta }_V}{2\rho }V(x,t)\hfill \\ \hfill & -\frac{\xi {\beta }_U}{\vartheta }U(x,t)-\frac{\xi {\beta }_L}{\vartheta }L(x,t)-\frac{\xi \left({\beta }_A^{*}-{\gamma }^{*}\right)}{\vartheta }A(x,t)\hfill \\ \hfill & \le \upsilon +\frac{\xi }{\vartheta }\varkappa -\varsigma \left[X(x,t)+N(x,t)+Y(x,t)+\frac{{\beta }_Y}{2\rho }V(x,t)\right.\hfill \\ \hfill & \left.+\frac{\xi }{\vartheta }\left(U(x,t)+L(x,t)+A(x,t)\right)\right]=\upsilon +\frac{\xi }{\vartheta }\varkappa -\varsigma \mathrm{\Psi }(x,t),\hfill \end{array}$$

where $\varsigma =min\left\{{\beta }_X,{\beta }_N,\frac{{\beta }_Y}{2},{\beta }_V,{\beta }_U,{\beta }_L,{\beta }_A^{*}-{\gamma }^{*}\right\}$. Hence, $\mathrm{\Psi }(x,t)$ satisfies

$$\left\{\begin{array}{c}\frac{\partial \mathrm{\Psi }(x,t)}{\partial t}-\tilde{\delta }\Delta \mathrm{\Psi }(x,t)\le \upsilon +\frac{\xi }{\vartheta }\varkappa -\varsigma \mathrm{\Psi }(x,t),\hfill \\ \frac{\partial \mathrm{\Psi }}{\partial \overrightarrow{\eta }}=0,\hfill \\ \mathrm{\Psi }(x,0)\ge 0.\hfill \end{array}\right.$$

Let $\tilde{\mathrm{\Psi }}\left(t\right)$ be a solution of the ODE:

$$\left\{\begin{array}{c}\frac{d\tilde{\mathrm{\Psi }}\left(t\right)}{dt}=\upsilon +\frac{\xi }{\vartheta }\varkappa -\varsigma \tilde{\mathrm{\Psi }}\left(t\right),\hfill \\ \tilde{\mathrm{\Psi }}\left(0\right)=\underset{x\in \overline{\mathrm{\Omega }}}{max}\mathrm{\Psi }(x,0).\hfill \end{array}\right.$$

This yields $\tilde{\mathrm{\Psi }}\left(t\right)\le max\left\{\underset{x\in \overline{\mathrm{\Omega }}}{max}\mathrm{\Psi }(x,0),\frac{\upsilon +\frac{\xi }{\vartheta }\varkappa }{\varsigma }\right\}.$ Applying the comparison principle [70], we have $\mathrm{\Psi }(x,t)\le \tilde{\mathrm{\Psi }}\left(t\right).$ Hence, we obtain

$$\mathrm{\Psi }(x,t)\le max\left\{\underset{x\in \overline{\mathrm{\Omega }}}{max}\mathrm{\Psi }(x,0),\frac{\upsilon +\frac{\xi }{\vartheta }\varkappa }{\varsigma }\right\},$$

which ensures the boundedness of $X(x,t),N(x,t),Y(x,t)$,$V(x,t),U(x,t),L(x,t)$ and $A(x,t)$ on $\overline{\mathrm{\Omega }}\times [0,{\chi }_m)$. From the standard theory for semi-linear parabolic PDEs, we obtain ${\chi }_m=+\infty$ [71]. □

4. Steady States and Thresholds

To find the steady states, we use

$$\begin{array}{cc}\hfill 0& =\upsilon -{\beta }_{X}X-\sigma VX,\hfill \\ \hfill 0& =\sigma VX-(\varrho +{\beta }_N)N,\hfill \\ \hfill 0& =\varrho N-{\beta }_{Y}Y-\xi YU,\hfill \\ \hfill 0& =\rho Y-{\beta }_{V}V,\hfill \\ \hfill 0& =\varkappa +\vartheta YU-{\beta }_{U}U-\tau AU,\hfill \\ \hfill 0& =\tau AU+\gamma A-\left(\lambda +{\beta }_L\right)L,\hfill \\ \hfill 0& =\lambda L-{\beta }_{A}A,\hfill \end{array}$$

and solve this algebraic system. We found that the system admits four steady states:

(i) Healthy steady state: ${\mathrm{\Xi }}_0=\left(X_0,0,0,0,U_0,0,0\right),$ where $X_0=\frac{\upsilon }{{\beta }_X}$ and $U_0=\frac{\varkappa }{{\beta }_U}$.

(ii) HTLV-I single-infection steady state, ${\mathrm{\Xi }}_1=\left(X_1,0,0,0,U_1,L_1,A_1\right),$ where

$$\begin{array}{cc}\hfill X_1& =X_0=\frac{\upsilon }{{\beta }_X},U_1=\frac{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}{\lambda \tau }=\frac{U_0}{R_1},\hfill \\ \hfill L_1& =\frac{{\beta }_U{\beta }_A}{\tau \lambda }\left[\frac{\lambda \varkappa \tau }{{\beta }_U\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}-1\right]=\frac{{\beta }_U{\beta }_A}{\tau \lambda }\left(R_1-1\right),\hfill \\ \hfill A_1& =\frac{{\beta }_U}{\tau }\left[\frac{\lambda \varkappa \tau }{{\beta }_U\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}-1\right]=\frac{{\beta }_U}{\tau }\left(R_1-1\right),\hfill \end{array}$$

where

$$R_1=\frac{\lambda \varkappa \tau }{{\beta }_U\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}.$$

(35)

Here, $R_1$ is the basic HTLV-I single-infection reproduction number. It decides the establishment of HTLV-I single infection. Clearly, $X_1$ is always positive. Additionally, because ${\beta }_A-\gamma >0$, $U_1$ and $R_1$ are always positive, while $L_1$ and $A_1$ are positive if $R_1>1.$ Therefore, ${\mathrm{\Xi }}_1$ exists when $R_1>1.$

(iii) SARS-CoV-2 single-infection steady state, ${\mathrm{\Xi }}_2=\left(X_2,N_2,Y_2,V_2,U_2,0,0\right),$ where

$$Y_2=\frac{{\beta }_V}{\rho }{V}_2,U_2=\frac{\rho \varkappa }{\rho {\beta }_U-\vartheta {\beta }_V{V}_2},X_2=\frac{\left(\varrho +{\beta }_N\right)\left(Y_2{\beta }_Y+U_2{Y}_2\xi \right)}{\varrho \sigma V_2},N_2=\frac{Y_2{\beta }_Y+U_2{Y}_2\xi }{\varrho },$$

(36)

and $V_2$ satisfies the following equation:

$$\frac{T_1{V}_2^2+T_2{V}_2+T_3}{\rho \sigma \varrho \left(\rho {\beta }_U-\vartheta {\beta }_V{V}_2\right)}=0,$$

(37)

where

$$\begin{array}{cc}\hfill T_1& =\sigma {\beta }_Y{\beta }_V^2\vartheta \left(\varrho +{\beta }_N\right),\hfill \\ \hfill T_2& ={\beta }_X{\beta }_Y{\beta }_V^2\vartheta \left(\varrho +{\beta }_N\right)-\rho \sigma {\beta }_Y{\beta }_V{\beta }_U\left(\varrho +{\beta }_N\right)-\rho \sigma \varkappa \xi {\beta }_V\left(\varrho +{\beta }_N\right)-\upsilon \rho \sigma \varrho \vartheta {\beta }_V,\hfill \\ \hfill T_3& =\upsilon {\rho }^2\sigma \varrho {\beta }_U-\rho {\beta }_X{\beta }_Y{\beta }_V{\beta }_U\left(\varrho +{\beta }_N\right)-\rho \varkappa \xi {\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right).\hfill \end{array}$$

Now, we show that Equation (37) has a positive solution. Let us define

$$F\left(V\right)=\frac{T_1{V}^2+T_{2}V+T_3}{\rho \sigma \varrho \left(\rho {\beta }_U-\vartheta {\beta }_{V}V\right)}.$$

We have

$$F\left(0\right)=\frac{\upsilon {\rho }^2\sigma \varrho {\beta }_U-\rho {\beta }_X{\beta }_Y{\beta }_V{\beta }_U\left(\varrho +{\beta }_N\right)-\rho \varkappa \xi {\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)}{{\rho }^2\sigma \varrho {\beta }_U}=\frac{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y{\beta }_U+\varkappa \xi \right)}{\rho \sigma \varrho {\beta }_U}\left(R_2-1\right),$$

where

$$R_2=\frac{\upsilon \rho \sigma \varrho {\beta }_U}{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y{\beta }_U+\varkappa \xi \right)}.$$

(38)

Hence, $F\left(0\right)>0$ when $R_2>1.$ Moreover,

$$\underset{V⟶{\epsilon }^{-}}{lim}F\left(V\right)=-\infty ,$$

where $\epsilon =\frac{\rho {\beta }_U}{\vartheta {\beta }_V}$. Furthermore,

$$F^{\prime }\left(V\right)=-\frac{{\beta }_V(\varrho +{\beta }_N)}{\varrho \rho \sigma {\left(\rho {\beta }_U-V\vartheta {\beta }_V\right)}^2}\left[\sigma {\beta }_Y{\left(V\vartheta {\beta }_V-\rho {\beta }_U\right)}^2+\varkappa \xi {\beta }_X\vartheta \rho {\beta }_V+\varkappa \xi \sigma {\rho }^2{\beta }_U\right].$$

This shows that $F^{\prime }\left(V\right)<0$ for all $V\in (0,\epsilon )$. Then, there exists a unique $V_2\in (0,\epsilon )$, such that $F\left(V_2\right)=0$. From Equation (36), we have $Y_2>0,$$U_2>0,$$X_2>0$ and $N_2>0$. As a result, ${\mathrm{\Xi }}_2$ exists when $R_2>1.$ Here, $R_2$ denotes the basic SARS-CoV-2 single-infection reproduction number, and it decides the occurrence of SARS-CoV-2 single infection.

(iv) HTLV-I and SARS-CoV-2 co-infection steady state ${\mathrm{\Xi }}_3=\left(X_3,N_3,Y_3,V_3,U_3,L_3,A_3\right)$, where

$$\begin{array}{cc}\hfill X_3& =\frac{{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}{\rho \sigma \varrho \lambda \tau },\hfill \\ \hfill N_3& =\frac{{\beta }_X{\beta }_V\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}{\rho \sigma \varrho \lambda \tau }\left[\frac{\upsilon \rho \sigma \varrho \lambda \tau }{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}-1\right],\hfill \\ \hfill Y_3& =\frac{{\beta }_X{\beta }_V}{\rho \sigma }\left[\frac{\upsilon \rho \sigma \varrho \lambda \tau }{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}-1\right],\hfill \\ \hfill V_3& =\frac{{\beta }_X}{\sigma }\left[\frac{\upsilon \rho \sigma \varrho \lambda \tau }{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}-1\right],\hfill \\ \hfill U_3& =\frac{1}{\lambda \tau }\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right),\hfill \\ \hfill L_3& =\frac{{\beta }_A\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)}{\rho \sigma \tau \lambda }\times \hfill \\ \hfill & \left[\frac{\rho \sigma \lambda \tau }{\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U}\left(\frac{\varkappa }{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}+\frac{\upsilon \varrho \vartheta }{\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}\right)-1\right],\hfill \\ \hfill A_3& =\frac{\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)}{\rho \sigma \tau }\times \hfill \\ \hfill & \left[\frac{\rho \sigma \lambda \tau }{\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U}\left(\frac{\varkappa }{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}+\frac{\upsilon \varrho \vartheta }{\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}\right)-1\right].\hfill \end{array}$$

It follows that because ${\beta }_A-\gamma >0$, then $X_3$ and $U_3$ are always positive, while $N_3>0,$$Y_3>0$ and $V_3>0$ if $\frac{\upsilon \rho \sigma \varrho \lambda \tau }{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}>1$. On the other hand, $L_3>0$ and $A_3>0$ when $\frac{\rho \sigma \lambda \tau }{\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U}\left(\frac{\varkappa }{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}+\frac{\upsilon \varrho \vartheta }{\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}\right)>1$.

Therefore, we can rewrite the components of ${\mathrm{\Xi }}_3$ as

$$\begin{array}{cc}\hfill X_3& =\frac{X_0}{R_4},N_3=\frac{{\beta }_X{\beta }_V\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}{\rho \sigma \varrho \lambda \tau }\left(R_4-1\right),\hfill \\ \hfill Y_3& =\frac{{\beta }_X{\beta }_V}{\rho \sigma }\left(R_4-1\right),V_3=\frac{{\beta }_X}{\sigma }\left(R_4-1\right),\hfill \\ \hfill U_3& =\frac{1}{\lambda \tau }\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right),L_3=\frac{{\beta }_A\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)}{\rho \sigma \tau \lambda }\left(R_3-1\right),\hfill \\ \hfill A_3& =\frac{\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)}{\rho \sigma \tau }\left(R_3-1\right),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill R_3& =\frac{\rho \sigma \lambda \tau }{\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U}\left(\frac{\varkappa }{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}+\frac{\upsilon \varrho \vartheta }{\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}\right),\hfill \\ \hfill R_4& =\frac{\upsilon \rho \sigma \varrho \lambda \tau }{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right)}.\hfill \end{array}$$

Thus, ${\mathrm{\Xi }}_3$ exists when $R_3>1$ and $R_4>1.$ At this point, $R_3$ and $R_4$ are threshold numbers that determine the HTLV-I and SARS-CoV-2 coexistence.

Now we summarize the above results in the following Lemma.

Lemma 2.

There exist four threshold numbers $R_i$, $i=1,2,3,4$ such that:

(a)

The healthy steady state, ${\mathrm{\Xi }}_0=\left(X_0,0,0,0,U_0,0,0\right)$, always exists;

(b)

if $R_1>1$, then, in addition to ${\mathrm{\Xi }}_0$, there is an HTLV-I single-infection steady state, ${\mathrm{\Xi }}_1=\left(X_1,0,0,0,\right.$$\left.U_1,L_1,A_1\right)$;

(c)

If $R_2>1$, then, in addition to ${\mathrm{\Xi }}_0$, there is a SARS-CoV-2 single-infection steady state, ${\mathrm{\Xi }}_2=\left(X_2,N_2,Y_2,V_2,U_2,0,0\right)$;

(d)

If $R_3>1$ and $R_4>1$ then, in addition to ${\mathrm{\Xi }}_0$, there is anHTLV-I and SARS-CoV-2 co-infection steady state ${\mathrm{\Xi }}_3=\left(X_3,N_3,Y_3,V_3,U_3,L_3,A_3\right)$.

5. Global Stability

In this section, we examine the global stability of the four steady states ${\mathrm{\Xi }}_i$, $i=0,1,2,3$. We followed [72] to design Lyapunov function and used LaSalle’s invariance principle (LIP) [73]. We used the following geometric mean–arithmetic mean inequality:

$$\sqrt[n]{\left(s_1\right)\left(s_2\right)\dots \left(s_n\right)}\le \frac{s_1+s_2+\dots +s_n}{n},s_i\ge 0,i=1,2,\dots ,n.$$

(39)

Next, we denote $\left(X,N,Y,V,U,L,A\right)=\left(X\left(x,t\right),N\left(x,t\right),Y\left(x,t\right),V\left(x,t\right),U\left(x,t\right),L\left(x,t\right),\right.$$\left.A\left(x,t\right)\right)$. Let ${\mathrm{\Phi }}_j\left(x,t\right)=$${\mathrm{\Phi }}_j\left(X,N,Y,V,U,L,A\right)$ and define

$${\widehat{\mathrm{\Phi }}}_j\left(t\right)={\int }_{\mathrm{\Omega }}{\mathrm{\Phi }}_j\left(x,t\right)dx,j=0,1,2,3.$$

Let ${\mathrm{\Pi }}_j^{\prime }$ be the largest invariant subset of

$${\mathrm{\Pi }}_j=\left\{\left(X,N,Y,V,U,L,A\right):\frac{d{\widehat{\mathrm{\Phi }}}_j}{dt}=0\right\},j=0,1,2,3.$$

Theorem 1.

If $R_1\le 1$ and $R_2\le 1$, then, the healthy steady state ${\mathrm{\Xi }}_0$ is globally asymptotically stable (GAS).

Proof. 

We define ${\mathrm{\Phi }}_0\left(x,t\right)$ as:

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_0\left(x,t\right)& =X_0\left(\frac{X}{X_0}-1-ln\left(\frac{X}{X_0}\right)\right)+N+\frac{\varrho +{\beta }_N}{\varrho }Y+\frac{\sigma X_0}{{\beta }_V}V\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{U}_0\left(\frac{U}{U_0}-1-ln\left(\frac{U}{U_0}\right)\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }L+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }A.\hfill \end{array}$$

Clearly, for all $X,N,Y,V,U,L,A>0$, we have ${\mathrm{\Phi }}_0>0$; moreover, ${\mathrm{\Phi }}_0(X_0,0,0,0,U_0,0,$$0)=0$. We calculate $\frac{\partial {\mathrm{\Phi }}_0}{\partial t}$ as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_0}{\partial t}& =\left(1-\frac{X_0}{X}\right)\frac{\partial X}{\partial t}+\frac{\partial N}{\partial t}+\frac{\varrho +{\beta }_N}{\varrho }\frac{\partial Y}{\partial t}+\frac{\sigma X_0}{{\beta }_V}\frac{\partial V}{\partial t}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_0}{U}\right)\frac{\partial U}{\partial t}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{\partial L}{\partial t}+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\frac{\partial A}{\partial t}.\hfill \end{array}$$

From Equations (21)–(27), we write

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_0}{\partial t}& =\left(1-\frac{X_0}{X}\right)({\delta }_X\Delta X+\upsilon -{\beta }_{X}X-\sigma VX)+{\delta }_N\Delta N+\sigma VX-(\varrho +{\beta }_N)N\hfill \\ \hfill & +\frac{\varrho +{\beta }_N}{\varrho }\left({\delta }_Y\Delta Y+\varrho N-{\beta }_{Y}Y-\xi YU\right)+\frac{\sigma X_0}{{\beta }_V}\left({\delta }_V\Delta V+\rho Y-{\beta }_{V}V\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_0}{U}\right)\left({\delta }_U\Delta U+\varkappa +\vartheta YU-{\beta }_{U}U-\tau AU\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left({\delta }_L\Delta L+\tau AU+\gamma A-\left(\lambda +{\beta }_L\right)L\right)+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left({\delta }_A\Delta A+\lambda L-{\beta }_{A}A\right).\hfill \end{array}$$

(40)

Collecting the terms of (40), we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_0}{\partial t}& =\left(1-\frac{X_0}{X}\right)(\upsilon -{\beta }_{X}X)-\frac{\varrho +{\beta }_N}{\varrho }{\beta }_{Y}Y+\frac{\sigma X_0}{{\beta }_V}\rho Y+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_0}{U}\right)\left(\varkappa -{\beta }_{U}U\right)\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }Y{U}_0+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A{U}_0+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A-\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_{A}A\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_0}{X}\right)\Delta X+{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y+{\delta }_V\frac{\sigma X_0}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_0}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A.\hfill \end{array}$$

Because $\upsilon ={\beta }_X{X}_0$ and $\varkappa ={\beta }_U{U}_0,$ then

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_0}{\partial t}& =-{\beta }_X\frac{{\left(X-X_0\right)}^2}{X}+\left(\frac{\sigma X_0}{{\beta }_V}\rho -\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y-\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }{U}_0\right)Y-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_0\right)}^2}{U}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau U_0+\gamma -\frac{\left(\lambda +{\beta }_L\right)}{\lambda }{\beta }_A\right)A+{\delta }_X\left(1-\frac{X_0}{X}\right)\Delta X+{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y\hfill \\ \hfill & +{\delta }_V\frac{\sigma X_0}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_0}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A\hfill \\ \hfill & =-{\beta }_X\frac{{\left(X-X_0\right)}^2}{X}+\left(\frac{\sigma \upsilon }{{\beta }_X{\beta }_V}\rho -\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y-\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }\frac{\varkappa }{{\beta }_U}\right)Y-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_0\right)}^2}{U}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau \frac{\varkappa }{{\beta }_U}-\frac{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}{\lambda }\right)A+{\delta }_X\left(1-\frac{X_0}{X}\right)\Delta X+{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y\hfill \\ \hfill & +{\delta }_V\frac{\sigma X_0}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_0}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A\hfill \\ \hfill & =-{\beta }_X\frac{{\left(X-X_0\right)}^2}{X}+\frac{\left(\varrho +{\beta }_N\right)\left({\beta }_Y{\beta }_U+\xi \varkappa \right)}{\varrho {\beta }_U}\left(\frac{\upsilon \rho \sigma \varrho {\beta }_U}{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left({\beta }_Y{\beta }_U+\xi \varkappa \right)}-1\right)Y\hfill \\ \hfill & -\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_0\right)}^2}{U}+\frac{\xi \left(\varrho +{\beta }_N\right)\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}{\lambda \vartheta \varrho }\left(\frac{\tau \lambda \varkappa }{{\beta }_U\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}-1\right)A\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_0}{X}\right)\Delta X+{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y+{\delta }_V\frac{\sigma X_0}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_0}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A.\hfill \end{array}$$

From Equations (35) and (38), we write

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_0}{\partial t}& =-{\beta }_X\frac{{\left(X-X_0\right)}^2}{X}+\frac{\left(\varrho +{\beta }_N\right)\left({\beta }_Y{\beta }_U+\xi \varkappa \right)}{\varrho {\beta }_U}\left(R_2-1\right)Y-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_0\right)}^2}{U}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}{\lambda \vartheta \varrho }\left(R_1-1\right)A+{\delta }_X\left(1-\frac{X_0}{X}\right)\Delta X+{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y\hfill \\ \hfill & +{\delta }_V\frac{\sigma X_0}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_0}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A.\hfill \end{array}$$

Consequently, we calculate the derivative of ${\widehat{\mathrm{\Phi }}}_0\left(x,t\right)$ with respect to t as:

$$\begin{array}{cc}\hfill \frac{d{\widehat{\mathrm{\Phi }}}_0}{dt}& ={\int }_{\mathrm{\Omega }}\frac{\partial {\mathrm{\Phi }}_0}{\partial t}dx=-{\beta }_X{\int }_{\mathrm{\Omega }}\frac{{\left(X-X_0\right)}^2}{X}dx\hfill \\ \hfill & +\frac{\left(\varrho +{\beta }_N\right)\left({\beta }_Y{\beta }_U+\xi \varkappa \right)}{\varrho {\beta }_U}\left(R_2-1\right){\int }_{\mathrm{\Omega }}Ydx-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{\int }_{\mathrm{\Omega }}\frac{{\left(U-U_0\right)}^2}{U}dx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}{\lambda \vartheta \varrho }\left(R_1-1\right){\int }_{\mathrm{\Omega }}Adx+{\delta }_X{\int }_{\mathrm{\Omega }}\left(1-\frac{X_0}{X}\right)\Delta Xdx+{\delta }_N{\int }_{\mathrm{\Omega }}\Delta Ndx\hfill \\ \hfill & +{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }{\int }_{\mathrm{\Omega }}\Delta Ydx+{\delta }_V\frac{\sigma X_0}{{\beta }_V}{\int }_{\mathrm{\Omega }}\Delta Vdx+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{\int }_{\mathrm{\Omega }}\left(1-\frac{U_0}{U}\right)\Delta Udx\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{\int }_{\mathrm{\Omega }}\Delta Ldx+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\int }_{\mathrm{\Omega }}\Delta Adx.\hfill \end{array}$$

(41)

Divergence theorem and Neumann boundary conditions imply that

$$\begin{array}{cc}\hfill 0& ={\int }_{\partial \mathrm{\Omega }}\nabla \mathrm{\Lambda }.\overrightarrow{\eta }dx={\int }_{\mathrm{\Omega }}div\left(\nabla \mathrm{\Lambda }\right)dx={\int }_{\mathrm{\Omega }}\Delta \mathrm{\Lambda }dx,\mathrm{for}\mathrm{\Lambda }\in \{X,N,Y,V,U,L,A\}.\hfill \\ \hfill 0& ={\int }_{\partial \mathrm{\Omega }}\frac{1}{\mathrm{\Lambda }}\nabla \mathrm{\Lambda }.\overrightarrow{\eta }dx={\int }_{\mathrm{\Omega }}div\left(\frac{1}{\mathrm{\Lambda }}\nabla \mathrm{\Lambda }\right)dx={\int }_{\mathrm{\Omega }}\left[\frac{1}{\mathrm{\Lambda }}\Delta \mathrm{\Lambda }-\frac{{∥\nabla \mathrm{\Lambda }∥}^2}{{\mathrm{\Lambda }}^2}\right]dx,\mathrm{for}\mathrm{\Lambda }\in \{X,N,Y,V,U,L,A\}.\hfill \end{array}$$

(42)

Hence,

$${\int }_{\mathrm{\Omega }}\frac{1}{\mathrm{\Lambda }}\Delta \mathrm{\Lambda }dx={\int }_{\mathrm{\Omega }}\frac{{∥\nabla \mathrm{\Lambda }∥}^2}{{\mathrm{\Lambda }}^2}dx\ge 0,\mathrm{where}\mathrm{\Lambda }=X,U>0.$$

(43)

Therefore, using Equations (42) and (43) in Equation (41) gives

$$\begin{array}{cc}\hfill \frac{d{\widehat{\mathrm{\Phi }}}_0}{dt}& =-{\beta }_X{\int }_{\mathrm{\Omega }}\frac{{\left(X-X_0\right)}^2}{X}dx+\frac{\left(\varrho +{\beta }_N\right)\left({\beta }_Y{\beta }_U+\xi \varkappa \right)}{\varrho {\beta }_U}\left(R_2-1\right){\int }_{\mathrm{\Omega }}Ydx-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{\int }_{\mathrm{\Omega }}\frac{{\left(U-U_0\right)}^2}{U}dx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}{\lambda \vartheta \varrho }\left(R_1-1\right){\int }_{\mathrm{\Omega }}Adx-{\delta }_X{X}_0{\int }_{\mathrm{\Omega }}\frac{{∥\nabla X∥}^2}{X^2}dx\hfill \\ \hfill & -\frac{\xi {\delta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{U}_0{\int }_{\mathrm{\Omega }}\frac{{∥\nabla U∥}^2}{U^2}dx.\hfill \end{array}$$

Therefore, if $R_1\le 1$ and $R_2\le 1,$ then $\frac{d{\widehat{\mathrm{\Phi }}}_0}{dt}\le 0$ for all $X,Y,U,A>0$ and $\frac{d{\widehat{\mathrm{\Phi }}}_0}{dt}=0$ when $X(x,t)=X_0$,$U(x,t)=U_0$ and $Y(x,t)=A(x,t)=0.$ The solutions of system (21)–(27) converge to ${\mathrm{\Pi }}_0^{\prime }$, which includes elements with $X(x,t)=X_0$,$U(x,t)=U_0$ and $Y(x,t)=A(x,t)=0;$ then, $\frac{\partial Y}{\partial t}=\frac{\partial A}{\partial t}=\Delta Y=\Delta A=0.$, Equations (23) and (27) yield

$$\begin{array}{cc}\hfill 0& =\frac{\partial Y}{\partial t}=\varrho N(x,t)⟹N\left(x,t\right)=0,\mathrm{for}\mathrm{all}\left(x,t\right),\hfill \\ \hfill 0& =\frac{\partial A}{\partial t}=\lambda L(x,t)⟹L\left(x,t\right)=0,\mathrm{for}\mathrm{all}\left(x,t\right).\hfill \end{array}$$

Furthermore, Equation (22) yields

$$0=\frac{\partial N}{\partial t}=\sigma X_{0}V\left(x,t\right)⟹V\left(x,t\right)=0,\mathrm{for}\mathrm{all}\left(x,t\right).$$

Therefore, ${\mathrm{\Pi }}_0^{\prime }=\left\{{\mathrm{\Xi }}_0\right\}.$ We deduce from LIP that ${\mathrm{\Xi }}_0$ is GAS [73]. □

The result of Theorem 1 shows that if there exist control parameters (e.g., drug therapies), which make $R_1\le 1$ and $R_2\le 1$, then both HTLV-I and SARS-CoV-2 are removed from the body regardless of the initial states.

Theorem 2.

If $R_1>1$ and $R_4\le 1$, then the HTLV-I single-infection steady state ${\mathrm{\Xi }}_1$ is GAS.

Proof. 

Let ${\mathrm{\Phi }}_1$ be defined as:

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_1& =X_1\left(\frac{X}{X_1}-1-ln\left(\frac{X}{X_1}\right)\right)+N+\frac{\varrho +{\beta }_N}{\varrho }Y+\frac{\sigma X_1}{{\beta }_V}V\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{U}_1\left(\frac{U}{U_1}-1-ln\left(\frac{U}{U_1}\right)\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{L}_1\left(\frac{L}{L_1}-1-ln\left(\frac{L}{L_1}\right)\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{A}_1\left(\frac{A}{A_1}-1-ln\left(\frac{A}{A_1}\right)\right).\hfill \end{array}$$

Clearly, for all $X,N,Y,V,U,L,A>0$, we have ${\mathrm{\Phi }}_1>0$ and ${\mathrm{\Phi }}_1\left(X_1,0,0,0,U_1,L_1,A_1\right)=0$. We calculate $\frac{\partial {\mathrm{\Phi }}_1}{\partial t}$ as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_1}{\partial t}& =\left(1-\frac{X_1}{X}\right)\frac{\partial X}{\partial t}+\frac{\partial N}{\partial t}+\frac{\varrho +{\beta }_N}{\varrho }\frac{\partial Y}{\partial t}+\frac{\sigma X_1}{{\beta }_V}\frac{\partial V}{\partial t}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\frac{\partial U}{\partial t}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_1}{L}\right)\frac{\partial L}{\partial t}+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_1}{A}\right)\frac{\partial A}{\partial t}.\hfill \end{array}$$

From Equations (21)–(27), we write

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_1}{\partial t}& =\left(1-\frac{X_1}{X}\right)({\delta }_X\Delta X+\upsilon -{\beta }_{X}X-\sigma VX)+{\delta }_N\Delta N+\sigma VX-(\varrho +{\beta }_N)N\hfill \\ \hfill & +\frac{\varrho +{\beta }_N}{\varrho }\left({\delta }_Y\Delta Y+\varrho N-{\beta }_{Y}Y-\xi YU\right)+\frac{\sigma X_1}{{\beta }_V}\left({\delta }_V\Delta V+\rho Y-{\beta }_{V}V\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\left({\delta }_U\Delta U+\varkappa +\vartheta YU-{\beta }_{U}U-\tau AU\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_1}{L}\right)\left({\delta }_L\Delta L+\tau AU+\gamma A-\left(\lambda +{\beta }_L\right)L\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_1}{A}\right)\left({\delta }_A\Delta A+\lambda L-{\beta }_{A}A\right).\hfill \end{array}$$

(44)

Equation (44) is simplified as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_1}{\partial t}& =\left(1-\frac{X_1}{X}\right)\left(\upsilon -{\beta }_{X}X\right)-\frac{\varrho +{\beta }_N}{\varrho }{\beta }_{Y}Y+\frac{\sigma X_1}{{\beta }_V}\rho Y+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\left(\varkappa -{\beta }_{U}U\right)\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }Y{U}_1+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A{U}_1-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau AU\frac{L_1}{L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A\frac{L_1}{L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\lambda +{\beta }_L\right)L_1-\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\vartheta \varrho }L\frac{A_1}{A}\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_{A}A+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_A{A}_1+{\delta }_X\left(1-\frac{X_1}{X}\right)\Delta X+{\delta }_N\Delta N\hfill \\ \hfill & +{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y+{\delta }_V\frac{\sigma X_1}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_1}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_1}{A}\right)\Delta A\hfill \\ \hfill & =\left(1-\frac{X_1}{X}\right)\left(\upsilon -{\beta }_{X}X\right)+\left(\frac{\sigma X_1}{{\beta }_V}\rho -\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y-\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }{U}_1\right)Y\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\left(\varkappa -{\beta }_{U}U\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_1{U}_1\frac{A}{A_1}-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_1{U}_1\frac{AU{L}_1}{A_1{U}_{1}L}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_1\frac{A}{A_1}-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_1\frac{A{L}_1}{A_{1}L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\lambda +{\beta }_L\right)L_1\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\vartheta \varrho }{L}_1\frac{L{A}_1}{L_{1}A}-\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_A{A}_1\frac{A}{A_1}+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_A{A}_1\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_1}{X}\right)\Delta X+{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y+{\delta }_V\frac{\sigma X_1}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_1}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_1}{A}\right)\Delta A.\hfill \end{array}$$

Using the steady-state conditions for ${\mathrm{\Xi }}_1$:

$$\begin{array}{cc}\hfill \upsilon & ={\beta }_X{X}_1,\hfill \\ \hfill \varkappa & ={\beta }_U{U}_1+\tau A_1{U}_1,\hfill \\ \hfill \tau A_1{U}_1& =\left(\lambda +{\beta }_L\right)L_1-\gamma A_1,\hfill \\ \hfill \lambda L_1& ={\beta }_A{A}_1,\hfill \end{array}$$

we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_1}{\partial t}& =-{\beta }_X\frac{{\left(X-X_1\right)}^2}{X}+\left(\frac{\sigma X_1}{{\beta }_V}\rho -\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y-\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }{U}_1\right)Y-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_1\right)}^2}{U}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\tau A_1{U}_1+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_1{U}_1\frac{A}{A_1}-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_1{U}_1\frac{AU{L}_1}{A_1{U}_{1}L}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_1\frac{A}{A_1}-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_1\frac{A{L}_1}{A_{1}L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau A_1{U}_1+\gamma A_1\right)\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau A_1{U}_1+\gamma A_1\right)\frac{L{A}_1}{L_{1}A}-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau A_1{U}_1+\gamma A_1\right)\frac{A}{A_1}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau A_1{U}_1+\gamma A_1\right)\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_1}{X}\right)\Delta X+{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y+{\delta }_V\frac{\sigma X_1}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_1}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_1}{A}\right)\Delta A.\hfill \end{array}$$

Collect the terms as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_1}{\partial t}& =-{\beta }_X\frac{{\left(X-X_1\right)}^2}{X}-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_1\right)}^2}{U}\hfill \\ \hfill & +\frac{\left(\varrho +{\beta }_N\right)\left[{\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right]}{\varrho \lambda \tau }\left(\frac{\upsilon \rho \sigma \varrho \lambda \tau }{{\beta }_X{\beta }_V\left(\varrho +{\beta }_N\right)\left[{\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right]}-1\right)Y\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_1{U}_1\left(3-\frac{U_1}{U}-\frac{AU{L}_1}{A_1{U}_{1}L}-\frac{L{A}_1}{L_{1}A}\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_1\left(2-\frac{A{L}_1}{A_{1}L}-\frac{L{A}_1}{L_{1}A}\right)\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_1}{X}\right)\Delta X+{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y+{\delta }_V\frac{\sigma X_1}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_1}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_1}{A}\right)\Delta A.\hfill \end{array}$$

(45)

Equation (45) can be written as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_1}{\partial t}& =-{\beta }_X\frac{{\left(X-X_1\right)}^2}{X}-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_1\right)}^2}{U}\hfill \\ \hfill & +\frac{\left(\varrho +{\beta }_N\right)\left[{\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right]}{\varrho \lambda \tau }\left(R_4-1\right)Y\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_1{U}_1\left(3-\frac{U_1}{U}-\frac{AU{L}_1}{A_1{U}_{1}L}-\frac{L{A}_1}{L_{1}A}\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_1\left(2-\frac{A{L}_1}{A_{1}L}-\frac{L{A}_1}{L_{1}A}\right)+{\delta }_X\left(1-\frac{X_1}{X}\right)\Delta X\hfill \\ \hfill & +{\delta }_N\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\Delta Y+{\delta }_V\frac{\sigma X_1}{{\beta }_V}\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_1}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_1}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_1}{A}\right)\Delta A.\hfill \end{array}$$

Then, $\frac{d{\widehat{\mathrm{\Phi }}}_1}{dt}$ is calculated as:

$$\begin{array}{cc}\hfill \frac{d{\widehat{\mathrm{\Phi }}}_1}{dt}& ={\int }_{\mathrm{\Omega }}\frac{\partial {\mathrm{\Phi }}_1}{\partial t}dx=-{\beta }_X{\int }_{\mathrm{\Omega }}\frac{{\left(X-X_1\right)}^2}{X}dx-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{\int }_{\mathrm{\Omega }}\frac{{\left(U-U_1\right)}^2}{U}dx\hfill \\ \hfill & +\frac{\left(\varrho +{\beta }_N\right)\left[{\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right]}{\varrho \lambda \tau }\left(R_4-1\right){\int }_{\mathrm{\Omega }}Ydx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_1{U}_1{\int }_{\mathrm{\Omega }}\left(3-\frac{U_1}{U}-\frac{AU{L}_1}{A_1{U}_{1}L}-\frac{L{A}_1}{L_{1}A}\right)dx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_1{\int }_{\mathrm{\Omega }}\left(2-\frac{A{L}_1}{A_{1}L}-\frac{L{A}_1}{L_{1}A}\right)dx+{\delta }_X{\int }_{\mathrm{\Omega }}\left(1-\frac{X_1}{X}\right)\Delta Xdx+{\delta }_N{\int }_{\mathrm{\Omega }}\Delta Ndx\hfill \\ \hfill & +{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }{\int }_{\mathrm{\Omega }}\Delta Ydx+{\delta }_V\frac{\sigma X_1}{{\beta }_V}{\int }_{\mathrm{\Omega }}\Delta Vdx+\frac{\xi {\delta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{\int }_{\mathrm{\Omega }}\left(1-\frac{U_1}{U}\right)\Delta Udx\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{\int }_{\mathrm{\Omega }}\left(1-\frac{L_1}{L}\right)\Delta Ldx+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\int }_{\mathrm{\Omega }}\left(1-\frac{A_1}{A}\right)\Delta Adx.\hfill \end{array}$$

Using Equations (42) and (43), we obtain

$$\begin{array}{cc}\hfill \frac{d{\widehat{\mathrm{\Phi }}}_1}{dt}& =-{\beta }_X{\int }_{\mathrm{\Omega }}\frac{{\left(X-X_1\right)}^2}{X}dx-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{\int }_{\mathrm{\Omega }}\frac{{\left(U-U_1\right)}^2}{U}dx\hfill \\ \hfill & +\frac{\left(\varrho +{\beta }_N\right)\left[{\beta }_Y\lambda \tau +\xi \left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)\right]}{\varrho \lambda \tau }\left(R_4-1\right){\int }_{\mathrm{\Omega }}Ydx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_1{U}_1{\int }_{\mathrm{\Omega }}\left(3-\frac{U_1}{U}-\frac{AU{L}_1}{A_1{U}_{1}L}-\frac{L{A}_1}{L_{1}A}\right)dx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_1{\int }_{\mathrm{\Omega }}\left(2-\frac{A{L}_1}{A_{1}L}-\frac{L{A}_1}{L_{1}A}\right)dx-{\delta }_X{X}_1{\int }_{\mathrm{\Omega }}\frac{{∥\nabla X∥}^2}{X^2}dx\hfill \\ \hfill & -{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{U}_1{\int }_{\mathrm{\Omega }}\frac{{∥\nabla U∥}^2}{U^2}dx-{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{L}_1{\int }_{\mathrm{\Omega }}\frac{{∥\nabla L∥}^2}{L^2}dx\hfill \\ \hfill & -{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{A}_1{\int }_{\mathrm{\Omega }}\frac{{∥\nabla A∥}^2}{A^2}dx.\hfill \end{array}$$

Using inequality (39), we obtain

$$3\le \frac{U_1}{U}+\frac{AU{L}_1}{A_1{U}_{1}L}+\frac{L{A}_1}{L_{1}A},\forall U,L,A>0,$$

$$2\le \frac{A{L}_1}{A_{1}L}+\frac{L{A}_1}{L_{1}A},\forall L,A>0.$$

Therefore, if $R_4\le 1$, then $\frac{d{\widehat{\mathrm{\Phi }}}_1}{dt}\le 0$ for all $X,Y,U,L,A>0$, where $\frac{d{\widehat{\mathrm{\Phi }}}_1}{dt}=0$ at $X(x,t)=X_1,$ $Y(x,t)=0,$$U(x,t)=U_1,$$L(x,t)=L_1$, and $A(x,t)=A_1$. The solutions to system (21)–(27) tend to ${\mathrm{\Pi }}_1^{\prime }$, which has elements such that $X(x,t)=X_1,$$U(x,t)=U_1,$$L(x,t)=L_1$$A(x,t)=A_1$ and $Y(x,t)=0,$ then $\frac{\partial Y}{\partial t}=\Delta Y=0.$ Equation (23) yields

$$0=\frac{\partial Y}{\partial t}=\varrho N(x,t)⟹N\left(x,t\right)=0,\mathrm{for}\mathrm{all}\left(x,t\right).$$

Equation (22) gives

$$0=\frac{\partial N}{\partial t}=\sigma X_{1}V\left(x,t\right)⟹V\left(x,t\right)=0,\mathrm{for}\mathrm{all}\left(x,t\right).$$

Then, ${\mathrm{\Pi }}_1^{\prime }=\left\{{\mathrm{\Xi }}_1\right\}$, and LIP confirms the global stability of ${\mathrm{\Xi }}_1$. □

Theorem 2 suggests that if the model’s parameters are adjusted such that $R_1>1$ and $R_4\le 1$, then the SARS-CoV-2 infection will be extinct, and the patient will have chronic HTLV-I single infection.

Theorem 3.

If $R_2>1$ and $R_3\le 1$, then the SARS-CoV-2 single-infection steady state, ${\mathrm{\Xi }}_2$, is GAS.

Proof. 

We define ${\mathrm{\Phi }}_2$ as:

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_2& =X_2\left(\frac{X}{X_2}-1-ln\left(\frac{X}{X_2}\right)\right)+N_2\left(\frac{N}{N_2}-1-ln\left(\frac{N}{N_2}\right)\right)\hfill \\ \hfill & +\frac{\varrho +{\beta }_N}{\varrho }{Y}_2\left(\frac{Y}{Y_2}-1-ln\left(\frac{Y}{Y_2}\right)\right)+\frac{\sigma X_2}{{\beta }_V}{V}_2\left(\frac{V}{V_2}-1-ln\left(\frac{V}{V_2}\right)\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{U}_2\left(\frac{U}{U_2}-1-ln\left(\frac{U}{U_2}\right)\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }L+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }A.\hfill \end{array}$$

Clearly, for all $X,N,Y,V,U,L,A>0$, we have ${\mathrm{\Phi }}_2>0$ and ${\mathrm{\Phi }}_2\left(X_2,N_2,Y_2,V_2,U_2,0,0\right)=0$. Calculate $\frac{\partial {\mathrm{\Phi }}_2}{\partial t}$ as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_2}{\partial t}& =\left(1-\frac{X_2}{X}\right)\frac{\partial X}{\partial t}+\left(1-\frac{N_2}{N}\right)\frac{\partial N}{\partial t}+\frac{\varrho +{\beta }_N}{\varrho }\left(1-\frac{Y_2}{Y}\right)\frac{\partial Y}{\partial t}+\frac{\sigma X_2}{{\beta }_V}\left(1-\frac{V_2}{V}\right)\frac{\partial V}{\partial t}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\frac{\partial U}{\partial t}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{\partial L}{\partial t}+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\frac{\partial A}{\partial t}.\hfill \end{array}$$

From Equations (21)–(27), we write

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_2}{\partial t}& =\left(1-\frac{X_2}{X}\right)({\delta }_X\Delta X+\upsilon -{\beta }_{X}X-\sigma VX)+\left(1-\frac{N_2}{N}\right)\left({\delta }_N\Delta N+\sigma VX-(\varrho +{\beta }_N)N\right)\hfill \\ \hfill & +\frac{\varrho +{\beta }_N}{\varrho }\left(1-\frac{Y_2}{Y}\right)\left({\delta }_Y\Delta Y+\varrho N-{\beta }_{Y}Y-\xi YU\right)+\frac{\sigma X_2}{{\beta }_V}\left(1-\frac{V_2}{V}\right)\left({\delta }_V\Delta V+\rho Y-{\beta }_{V}V\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\left({\delta }_U\Delta U+\varkappa +\vartheta YU-{\beta }_{U}U-\tau AU\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left({\delta }_L\Delta L+\tau AU+\gamma A-\left(\lambda +{\beta }_L\right)L\right)+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left({\delta }_A\Delta A+\lambda L-{\beta }_{A}A\right).\hfill \end{array}$$

Collect the terms as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_2}{\partial t}& =\left(1-\frac{X_2}{X}\right)\left(\upsilon -{\beta }_{X}X\right)-\sigma VX\frac{N_2}{N}+(\varrho +{\beta }_N)N_2-\left(\varrho +{\beta }_N\right)N\frac{Y_2}{Y}-\frac{\varrho +{\beta }_N}{\varrho }{\beta }_{Y}Y+\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y{Y}_2\hfill \\ \hfill & +\frac{\varrho +{\beta }_N}{\varrho }\xi Y_{2}U+\frac{\sigma X_2}{{\beta }_V}\rho Y-\frac{\sigma X_2}{{\beta }_V}\rho Y\frac{V_2}{V}+\sigma X_2{V}_2+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\left(\varkappa -{\beta }_{U}U\right)\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }Y{U}_2+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A{U}_2+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A-\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_{A}A\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_2}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_2}{N}\right)\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_2}{Y}\right)\Delta Y+{\delta }_V\frac{\sigma X_2}{{\beta }_V}\left(1-\frac{V_2}{V}\right)\Delta V\hfill \\ \hfill & +{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A\hfill \\ \hfill & =\left(1-\frac{X_2}{X}\right)\left(\upsilon -{\beta }_{X}X\right)-\sigma VX\frac{N_2}{N}+(\varrho +{\beta }_N)N_2-\left(\varrho +{\beta }_N\right)N\frac{Y_2}{Y}+\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y{Y}_2\hfill \\ \hfill & +\left(\frac{\sigma X_2}{{\beta }_V}\rho -\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y-\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }{U}_2\right)Y+\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\xi Y_{2}U-\frac{\sigma X_2}{{\beta }_V}\rho Y\frac{V_2}{V}+\sigma X_2{V}_2\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\left(\varkappa -{\beta }_{U}U\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A{U}_2+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A-\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_{A}A\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_2}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_2}{N}\right)\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_2}{Y}\right)\Delta Y+{\delta }_V\frac{\sigma X_2}{{\beta }_V}\left(1-\frac{V_2}{V}\right)\Delta V\hfill \\ \hfill & +{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A.\hfill \end{array}$$

Using the steady-state conditions for ${\mathrm{\Xi }}_2$:

$$\begin{array}{cc}\hfill \upsilon & ={\beta }_X{X}_2+\sigma V_2{X}_2,\hfill \\ \hfill \sigma V_2{X}_2& =(\varrho +{\beta }_N)N_2,\hfill \\ \hfill \varrho N_2& ={\beta }_Y{Y}_2+\xi Y_2{U}_2,\hfill \\ \hfill \rho Y_2& ={\beta }_V{V}_2,\hfill \\ \hfill \varkappa & ={\beta }_U{U}_2-\vartheta Y_2{U}_2,\hfill \end{array}$$

we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_2}{\partial t}& =-{\beta }_X\frac{{\left(X-X_2\right)}^2}{X}+\sigma V_2{X}_2-\sigma V_2{X}_2\frac{X_2}{X}-\sigma V_2{X}_2\frac{VX{N}_2}{V_2{X}_{2}N}+\sigma V_2{X}_2-\sigma V_2{X}_2\frac{N{Y}_2}{N_{2}Y}\hfill \\ \hfill & +\sigma V_2{X}_2-\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\xi Y_2{U}_2+\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\xi Y_2{U}_2\frac{U}{U_2}-\sigma V_2{X}_2\frac{Y{V}_2}{Y_{2}V}+\sigma V_2{X}_2\hfill \\ \hfill & -\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_2\right)}^2}{U}-\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }{Y}_2{U}_2\left(1-\frac{U_2}{U}\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau U_2-\frac{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}{\lambda }\right)A\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_2}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_2}{N}\right)\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_2}{Y}\right)\Delta Y+{\delta }_V\frac{\sigma X_2}{{\beta }_V}\left(1-\frac{V_2}{V}\right)\Delta V\hfill \\ \hfill & +{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A.\hfill \end{array}$$

Collect the terms as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_2}{\partial t}& =-{\beta }_X\frac{{\left(X-X_2\right)}^2}{X}+\sigma V_2{X}_2\left(4-\frac{X_2}{X}-\frac{VX{N}_2}{V_2{X}_{2}N}-\frac{N{Y}_2}{N_{2}Y}-\frac{Y{V}_2}{Y_{2}V}\right)-\frac{\varrho +{\beta }_N}{\varrho }\xi Y_2{U}_2\left(2-\frac{U}{U_2}-\frac{U_2}{U}\right)\hfill \\ \hfill & -\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_2\right)}^2}{U}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau U_2-\frac{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}{\lambda }\right)A+{\delta }_X\left(1-\frac{X_2}{X}\right)\Delta X\hfill \\ \hfill & +{\delta }_N\left(1-\frac{N_2}{N}\right)\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_2}{Y}\right)\Delta Y+{\delta }_V\frac{\sigma X_2}{{\beta }_V}\left(1-\frac{V_2}{V}\right)\Delta V\hfill \\ \hfill & +{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A.\hfill \end{array}$$

We have

$$\begin{array}{cc}\hfill & -\frac{\varrho +{\beta }_N}{\varrho }\xi Y_2{U}_2\left(2-\frac{U}{U_2}-\frac{U_2}{U}\right)-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_2\right)}^2}{U}\hfill \\ \hfill & =\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\frac{\xi Y_2}{U}{\left(U-U_2\right)}^2-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_2\right)}^2}{U}\hfill \\ \hfill & =\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_2\right)}^2}{U}\left(\vartheta Y_2-{\beta }_U\right)\hfill \\ \hfill & =-\frac{\xi \varkappa \left(\varrho +{\beta }_N\right)}{\vartheta \varrho U_2}\frac{{\left(U-U_2\right)}^2}{U}.\hfill \end{array}$$

Then, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_2}{\partial t}& =-{\beta }_X\frac{{\left(X-X_2\right)}^2}{X}+\sigma V_2{X}_2\left(4-\frac{X_2}{X}-\frac{VX{N}_2}{V_2{X}_{2}N}-\frac{N{Y}_2}{N_{2}Y}-\frac{Y{V}_2}{Y_{2}V}\right)\hfill \\ \hfill & -\frac{\xi \varkappa \left(\varrho +{\beta }_N\right)}{\vartheta \varrho U_2}\frac{{\left(U-U_2\right)}^2}{U}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau U_2-\frac{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}{\lambda }\right)A\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_2}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_2}{N}\right)\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_2}{Y}\right)\Delta Y\hfill \\ \hfill & +{\delta }_V\frac{\sigma X_2}{{\beta }_V}\left(1-\frac{V_2}{V}\right)\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_2}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\Delta L\hfill \\ \hfill & +{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\Delta A.\hfill \end{array}$$

Computing $\frac{d{\widehat{\mathrm{\Phi }}}_2}{dt}$ and using Equations (42) and (43), we obtain

$$\begin{array}{cc}\hfill \frac{d{\widehat{\mathrm{\Phi }}}_2}{dt}& ={\int }_{\mathrm{\Omega }}\frac{\partial {\mathrm{\Phi }}_2}{\partial t}dx=-{\beta }_X{\int }_{\mathrm{\Omega }}\frac{{\left(X-X_2\right)}^2}{X}dx\hfill \\ \hfill & +\sigma V_2{X}_2{\int }_{\mathrm{\Omega }}\left(4-\frac{X_2}{X}-\frac{VX{N}_2}{V_2{X}_{2}N}-\frac{N{Y}_2}{N_{2}Y}-\frac{Y{V}_2}{Y_{2}V}\right)dx\hfill \\ \hfill & -\frac{\xi \varkappa \left(\varrho +{\beta }_N\right)}{\vartheta \varrho U_2}{\int }_{\mathrm{\Omega }}\frac{{\left(U-U_2\right)}^2}{U}dx+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau U_2-\frac{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}{\lambda }\right){\int }_{\mathrm{\Omega }}Adx\hfill \\ \hfill & -{\delta }_X{X}_2{\int }_{\mathrm{\Omega }}\frac{{∥\nabla X∥}^2}{X^2}dx-{\delta }_N{N}_2{\int }_{\mathrm{\Omega }}\frac{{∥\nabla N∥}^2}{N^2}dx-{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }{Y}_2{\int }_{\mathrm{\Omega }}\frac{{∥\nabla Y∥}^2}{Y^2}dx\hfill \\ \hfill & -{\delta }_V\frac{\sigma }{{\beta }_V}{X}_2{V}_2{\int }_{\mathrm{\Omega }}\frac{{∥\nabla V∥}^2}{V^2}dx+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{U}_2{\int }_{\mathrm{\Omega }}\frac{{∥\nabla U∥}^2}{U^2}dx.\hfill \end{array}$$

Hence, if $R_3\le 1$, then ${\mathrm{\Xi }}_3$ does not exist because $A_3\le 0$ and $L_3\le 0.$ Because both the ODE system (5)–(11) and the PDE system (21)–(27) have the same steady states, then ${\mathrm{\Xi }}_3$ does not exist for system (5)–(11). This implies that

$$\begin{array}{cc}\hfill \frac{dA\left(t\right)}{dt}& =\lambda L\left(t\right)-{\beta }_{A}A\left(t\right)\le 0,\hfill \\ \hfill \frac{dL\left(t\right)}{dt}& =\tau A\left(t\right)U\left(t\right)+\gamma A\left(t\right)-\left(\lambda +{\beta }_L\right)L\left(t\right)\le 0.\hfill \end{array}$$

It follows that $\left(\tau U\left(t\right)-\frac{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}{\lambda }\right)A\left(t\right)\le 0$ for all $A\left(t\right)>0$; thus, $\tau U_2-\frac{{\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)}{\lambda }\le 0$. By using inequality (39), we obtain

$$\frac{X_2}{X}+\frac{VX{N}_2}{V_2{X}_{2}N}+\frac{N{Y}_2}{N_{2}Y}+\frac{Y{V}_2}{Y_{2}V}\ge 4,\forall X,N,Y,V>0.$$

Thus, $\frac{d{\widehat{\mathrm{\Phi }}}_2}{dt}\le 0$ for all $X,N,Y,V,U,A>0$ and $\frac{d{\widehat{\mathrm{\Phi }}}_2}{dt}=0$ at $X\left(x,t\right)=X_2,$$N\left(x,t\right)=N_2,$$Y\left(x,t\right)=Y_2,$$V\left(x,t\right)=V_2,$$U\left(x,t\right)=U_2$ and $A\left(x,t\right)=0.$ The solutions of the model tend to ${\mathrm{\Pi }}_2^{\prime }$ for which $A\left(x,t\right)=0.$ It follows that $\frac{\partial A}{\partial t}=\Delta A=0$, and Equation (27) becomes

$$0=\frac{\partial A}{\partial t}=\lambda L\left(x,t\right)⟹L\left(x,t\right)=0\mathrm{for}\mathrm{all}\left(x,t\right).$$

Therefore, ${\mathrm{\Pi }}_2^{\prime }=\left\{{\mathrm{\Xi }}_2\right\}$. LIP implies the global stability of ${\mathrm{\Xi }}_2$. □

Theorem 3 suggests that if the model’s parameters are controlled such that $R_2>1$ and $R_3\le 1$, then the HTLV-I infection will be extinct, and the patient will have SARS-CoV-2 single infection.

We define

$$\tilde{R}=\frac{\rho \sigma \lambda \tau \varkappa }{\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}.$$

Theorem 4.

If $R_4>1$ and $1<R_3\le 1+\tilde{R}$, then the HTLV-I and SARS-CoV-2 co-infection steady state ${\mathrm{\Xi }}_3$ is GAS.

Proof. 

We define ${\mathrm{\Phi }}_3$ as:

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_3& =X_3\left(\frac{X}{X_3}-1-ln\left(\frac{X}{X_3}\right)\right)+N_3\left(\frac{N}{N_3}-1-ln\left(\frac{N}{N_3}\right)\right)\hfill \\ \hfill & +\frac{\varrho +{\beta }_N}{\varrho }{Y}_3\left(\frac{Y}{Y_3}-1-ln\left(\frac{Y}{Y_3}\right)\right)+\frac{\sigma X_3}{{\beta }_V}{V}_3\left(\frac{V}{V_3}-1-ln\left(\frac{V}{V_3}\right)\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{U}_3\left(\frac{U}{U_3}-1-ln\left(\frac{U}{U_3}\right)\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{L}_3\left(\frac{L}{L_3}-1-ln\left(\frac{L}{L_3}\right)\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{A}_3\left(\frac{A}{A_3}-1-ln\left(\frac{A}{A_3}\right)\right).\hfill \end{array}$$

Clearly, for all $X,N,Y,V,U,L,A>0$, we have ${\mathrm{\Phi }}_3>0$ and ${\mathrm{\Phi }}_3(X_3,N_3,Y_3,V_3,U_3,L_3,$$A_3)=0$. Calculate $\frac{\partial {\mathrm{\Phi }}_3}{\partial t}$ as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_3}{\partial t}& =\left(1-\frac{X_3}{X}\right)\frac{\partial X}{\partial t}+\left(1-\frac{N_3}{N}\right)\frac{\partial N}{\partial t}+\frac{\varrho +{\beta }_N}{\varrho }\left(1-\frac{Y_3}{Y}\right)\frac{\partial Y}{\partial t}+\frac{\sigma X_3}{{\beta }_V}\left(1-\frac{V_3}{V}\right)\frac{\partial V}{\partial t}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\frac{\partial U}{\partial t}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_3}{L}\right)\frac{\partial L}{\partial t}+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_3}{A}\right)\frac{\partial A}{\partial t}.\hfill \end{array}$$

Substituting from Equations (21)–(27), we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_3}{\partial t}& =\left(1-\frac{X_3}{X}\right)({\delta }_X\Delta X+\upsilon -{\beta }_{X}X-\sigma VX)+\left(1-\frac{N_3}{N}\right)\left({\delta }_N\Delta N+\sigma VX-(\varrho +{\beta }_N)N\right)\hfill \\ \hfill & +\frac{\varrho +{\beta }_N}{\varrho }\left(1-\frac{Y_3}{Y}\right)\left({\delta }_Y\Delta Y+\varrho N-{\beta }_{Y}Y-\xi YU\right)+\frac{\sigma X_3}{{\beta }_V}\left(1-\frac{V_3}{V}\right)\left({\delta }_V\Delta V+\rho Y-{\beta }_{V}V\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\left({\delta }_U\Delta U+\varkappa +\vartheta YU-{\beta }_{U}U-\tau AU\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_3}{L}\right)\left({\delta }_L\Delta L+\tau AU+\gamma A-\left(\lambda +{\beta }_L\right)L\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_3}{A}\right)\left({\delta }_A\Delta A+\lambda L-{\beta }_{A}A\right).\hfill \end{array}$$

Collecting the terms, we have

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_3}{\partial t}& =\left(1-\frac{X_3}{X}\right)\left(\upsilon -{\beta }_{X}X\right)-\sigma VX\frac{N_3}{N}+(\varrho +{\beta }_N)N_3-\left(\varrho +{\beta }_N\right)N\frac{Y_3}{Y}-\frac{\varrho +{\beta }_N}{\varrho }{\beta }_{Y}Y\hfill \\ \hfill & +\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y{Y}_3+\frac{\varrho +{\beta }_N}{\varrho }\xi Y_{3}U+\frac{\sigma X_3}{{\beta }_V}\rho Y-\frac{\sigma X_3}{{\beta }_V}\rho Y\frac{V_3}{V}+\sigma X_3{V}_3\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\left(\varkappa -{\beta }_{U}U\right)-\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }Y{U}_3+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A{U}_3\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau AU\frac{L_3}{L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A\frac{L_3}{L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\lambda +{\beta }_L\right)L_3\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\vartheta \varrho }L\frac{A_3}{A}-\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_{A}A+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_A{A}_3\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_3}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_3}{N}\right)\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_3}{Y}\right)\Delta Y\hfill \\ \hfill & +{\delta }_V\frac{\sigma }{{\beta }_V}{X}_3\left(1-\frac{V_3}{V}\right)\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_3}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_3}{A}\right)\Delta A.\hfill \end{array}$$

(46)

Equation (46) can be written as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_3}{\partial t}& =\left(1-\frac{X_3}{X}\right)\left(\upsilon -{\beta }_{X}X\right)-\sigma VX\frac{N_3}{N}+(\varrho +{\beta }_N)N_3-\left(\varrho +{\beta }_N\right)N\frac{Y_3}{Y}+\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y{Y}_3\hfill \\ \hfill & +\left(\frac{\sigma X_3}{{\beta }_V}\rho -\frac{\varrho +{\beta }_N}{\varrho }{\beta }_Y-\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }{U}_3\right)Y+\frac{\varrho +{\beta }_N}{\varrho }\xi Y_{3}U-\frac{\sigma X_3}{{\beta }_V}\rho Y\frac{V_3}{V}+\sigma X_3{V}_3\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\left(\varkappa -{\beta }_{U}U\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A{U}_3-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau AU\frac{L_3}{L}\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A\frac{L_3}{L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\lambda +{\beta }_L\right)L_3-\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\vartheta \varrho }L\frac{A_3}{A}\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_{A}A+\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{\beta }_A{A}_3+{\delta }_X\left(1-\frac{X_3}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_3}{N}\right)\Delta N\hfill \\ \hfill & +{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_3}{Y}\right)\Delta Y+{\delta }_V\frac{\sigma }{{\beta }_V}{X}_3\left(1-\frac{V_3}{V}\right)\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_3}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_3}{A}\right)\Delta A.\hfill \end{array}$$

Utilizing the steady-state conditions for ${\mathrm{\Xi }}_3$:

$$\begin{array}{cc}\hfill \upsilon & ={\beta }_X{X}_3+\sigma V_3{X}_3,\hfill \\ \hfill \sigma V_3{X}_3& =(\varrho +{\beta }_N)N_3,\hfill \\ \hfill \varrho N_3& ={\beta }_Y{Y}_3+\xi Y_3{U}_3,\hfill \\ \hfill \rho Y_3& ={\beta }_V{V}_3,\hfill \\ \hfill \varkappa & ={\beta }_U{U}_3-\vartheta Y_3{U}_3+\tau A_3{U}_3,\hfill \\ \hfill \tau A_3{U}_3& =\left(\lambda +{\beta }_L\right)L_3-\gamma A_3,\hfill \\ \hfill \lambda L_3& ={\beta }_A{A}_3,\hfill \end{array}$$

we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_3}{\partial t}& =-{\beta }_X\frac{{\left(X-X_3\right)}^2}{X}+\sigma V_3{X}_3-\sigma V_3{X}_3\frac{X_3}{X}-\sigma V_3{X}_3\frac{VX{N}_3}{V_3{X}_{3}N}+\sigma V_3{X}_3-\sigma V_3{X}_3\frac{N{Y}_3}{N_{3}Y}\hfill \\ \hfill & +\sigma V_3{X}_3-\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\xi Y_3{U}_3+\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\xi Y_3{U}_3\frac{U}{U_3}-\sigma V_3{X}_3\frac{Y{V}_3}{Y_{3}V}+\sigma X_3{V}_3-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_3\right)}^2}{U}\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }{Y}_3{U}_3\left(1-\frac{U_3}{U}\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_3{U}_3\left(1-\frac{U_3}{U}\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_3{U}_3\frac{A}{A_3}\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_3{U}_3\frac{AU{L}_3}{A_3{U}_{3}L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_3\frac{A}{A_3}-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_3\frac{A{L}_3}{A_{3}L}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau A_3{U}_3+\gamma A_3\right)\hfill \\ \hfill & -\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau A_3{U}_3+\gamma A_3\right)\frac{L{A}_3}{L_{3}A}-\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau A_3{U}_3+\gamma A_3\right)\frac{A}{A_3}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(\tau A_3{U}_3+\gamma A_3\right)\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_3}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_3}{N}\right)\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_3}{Y}\right)\Delta Y+{\delta }_V\frac{\sigma }{{\beta }_V}{X}_3\left(1-\frac{V_3}{V}\right)\Delta V\hfill \\ \hfill & +{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\Delta U+{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_3}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_3}{A}\right)\Delta A.\hfill \end{array}$$

Collect the terms as:

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_3}{\partial t}& =-{\beta }_X\frac{{\left(X-X_3\right)}^2}{X}+\sigma V_3{X}_3\left(4-\frac{X_3}{X}-\frac{VX{N}_3}{V_3{X}_{3}N}-\frac{N{Y}_3}{N_{3}Y}-\frac{Y{V}_3}{Y_{3}V}\right)-\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\xi Y_3{U}_3\left(2-\frac{U}{U_3}-\frac{U_3}{U}\right)\hfill \\ \hfill & -\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_3\right)}^2}{U}+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_3{U}_3\left(3-\frac{U_3}{U}-\frac{AU{L}_3}{A_3{U}_{3}L}-\frac{L{A}_3}{L_{3}A}\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_3\left(2-\frac{A{L}_3}{A_{3}L}-\frac{L{A}_3}{L_{3}A}\right)+{\delta }_X\left(1-\frac{X_3}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_3}{N}\right)\Delta N\hfill \\ \hfill & +{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_3}{Y}\right)\Delta Y+{\delta }_V\frac{\sigma }{{\beta }_V}{X}_3\left(1-\frac{V_3}{V}\right)\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_3}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_3}{A}\right)\Delta A.\hfill \end{array}$$

We have

$$\begin{array}{cc}\hfill & -\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\xi Y_3{U}_3\left(2-\frac{U}{U_3}-\frac{U_3}{U}\right)-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_3\right)}^2}{U}\hfill \\ \hfill & =\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\frac{\xi Y_3}{U}{\left(U-U_3\right)}^2-\frac{\xi {\beta }_U\left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\frac{{\left(U-U_3\right)}^2}{U}\hfill \\ \hfill & =\frac{\xi \left(\varrho +{\beta }_N\right)}{\varrho }\frac{{\left(U-U_3\right)}^2}{U}\left(Y_3-\frac{{\beta }_U}{\vartheta }\right)\hfill \\ \hfill & =\frac{\xi \left(\varrho +{\beta }_N\right)\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)}{\rho \sigma \varrho \vartheta }\frac{{\left(U-U_3\right)}^2}{U}\left(R_3-\frac{\rho \sigma \lambda \tau \varkappa }{\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)\left({\beta }_L{\beta }_A+\lambda \left({\beta }_A-\gamma \right)\right)}-1\right)\hfill \\ \hfill & =\frac{\xi \left(\varrho +{\beta }_N\right)\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)}{\rho \sigma \varrho \vartheta }\frac{{\left(U-U_3\right)}^2}{U}\left(R_3-\tilde{R}-1\right).\hfill \end{array}$$

Then, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathrm{\Phi }}_3}{\partial t}& =-{\beta }_X\frac{{\left(X-X_3\right)}^2}{X}+\sigma V_3{X}_3\left(4-\frac{X_3}{X}-\frac{VX{N}_3}{V_3{X}_{3}N}-\frac{N{Y}_3}{N_{3}Y}-\frac{Y{V}_3}{Y_{3}V}\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)}{\rho \sigma \varrho \vartheta }\frac{{\left(U-U_3\right)}^2}{U}\left(R_3-\tilde{R}-1\right)\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_3{U}_3\left(3-\frac{U_3}{U}-\frac{AU{L}_3}{A_3{U}_{3}L}-\frac{L{A}_3}{L_{3}A}\right)+\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_3\left(2-\frac{A{L}_3}{A_{3}L}-\frac{L{A}_3}{L_{3}A}\right)\hfill \\ \hfill & +{\delta }_X\left(1-\frac{X_3}{X}\right)\Delta X+{\delta }_N\left(1-\frac{N_3}{N}\right)\Delta N+{\delta }_Y\frac{\left(\varrho +{\beta }_N\right)}{\varrho }\left(1-\frac{Y_3}{Y}\right)\Delta Y\hfill \\ \hfill & +{\delta }_V\frac{\sigma }{{\beta }_V}{X}_3\left(1-\frac{V_3}{V}\right)\Delta V+{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{U_3}{U}\right)\Delta U\hfill \\ \hfill & +{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\left(1-\frac{L_3}{L}\right)\Delta L+{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }\left(1-\frac{A_3}{A}\right)\Delta A.\hfill \end{array}$$

Computing $\frac{d{\widehat{\mathrm{\Phi }}}_3}{dt}$ and using Equations (42) and (43), we obtain

$$\begin{array}{cc}\hfill \frac{d{\widehat{\mathrm{\Phi }}}_3}{dt}& ={\int }_{\mathrm{\Omega }}\frac{\partial {\mathrm{\Phi }}_3}{\partial t}dx=-{\beta }_X{\int }_{\mathrm{\Omega }}\frac{{\left(X-X_3\right)}^2}{X}dx+\sigma V_3{X}_3{\int }_{\mathrm{\Omega }}\left(4-\frac{X_3}{X}-\frac{VX{N}_3}{V_3{X}_{3}N}-\frac{N{Y}_3}{N_{3}Y}-\frac{Y{V}_3}{Y_{3}V}\right)dx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)\left(\vartheta {\beta }_X{\beta }_V+\rho \sigma {\beta }_U\right)}{\rho \sigma \varrho \vartheta }\left(R_3-\tilde{R}-1\right){\int }_{\mathrm{\Omega }}\frac{{\left(U-U_3\right)}^2}{U}dx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\tau A_3{U}_3{\int }_{\mathrm{\Omega }}\left(3-\frac{U_3}{U}-\frac{AU{L}_3}{A_3{U}_{3}L}-\frac{L{A}_3}{L_{3}A}\right)dx\hfill \\ \hfill & +\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }\gamma A_3{\int }_{\mathrm{\Omega }}\left(2-\frac{A{L}_3}{A_{3}L}-\frac{L{A}_3}{L_{3}A}\right)dx-{\delta }_X{X}_3{\int }_{\mathrm{\Omega }}\frac{{∥\nabla X∥}^2}{X^2}dx-{\delta }_N{N}_3{\int }_{\mathrm{\Omega }}\frac{{∥\nabla N∥}^2}{N^2}dx\hfill \\ \hfill & -\frac{{\delta }_Y\left(\varrho +{\beta }_N\right)}{\varrho }{Y}_3{\int }_{\mathrm{\Omega }}\frac{{∥\nabla Y∥}^2}{Y^2}dx-{\delta }_V\frac{\sigma }{{\beta }_V}{X}_3{V}_3{\int }_{\mathrm{\Omega }}\frac{{∥\nabla V∥}^2}{V^2}dx-{\delta }_U\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{U}_3{\int }_{\mathrm{\Omega }}\frac{{∥\nabla U∥}^2}{U^2}dx\hfill \\ \hfill & -{\delta }_L\frac{\xi \left(\varrho +{\beta }_N\right)}{\vartheta \varrho }{L}_3{\int }_{\mathrm{\Omega }}\frac{{∥\nabla L∥}^2}{L^2}dx-{\delta }_A\frac{\xi \left(\varrho +{\beta }_N\right)\left(\lambda +{\beta }_L\right)}{\lambda \vartheta \varrho }{A}_3{\int }_{\mathrm{\Omega }}\frac{{∥\nabla A∥}^2}{A^2}dx.\hfill \end{array}$$

Using inequality (39), we obtain

$$\begin{array}{cc}\hfill 4& \le \frac{X_3}{X}+\frac{VX{N}_3}{V_3{X}_{3}N}+\frac{N{Y}_3}{N_{3}Y}+\frac{Y{V}_3}{Y_{3}V},\forall X,N,Y,V>0,\hfill \\ \hfill 3& \le \frac{U_3}{U}+\frac{AU{L}_3}{A_3{U}_{3}L}+\frac{L{A}_3}{L_{3}A},\forall U,L,A>0,\hfill \\ \hfill 2& \le \frac{A{L}_3}{A_{3}L}+\frac{L{A}_3}{L_{3}A},\forall L,A>0.\hfill \end{array}$$

Because $1<R_3\le 1+\tilde{R},$ we have $\frac{d{\widehat{\mathrm{\Phi }}}_3}{dt}\le 0$; moreover, $\frac{d{\widehat{\mathrm{\Phi }}}_3}{dt}=0$ when $X\left(x,t\right)=X_3,$ $N\left(x,t\right)=N_3, Y\left(x,t\right)=Y_3,$$V\left(x,t\right)=V_3,$$U\left(x,t\right)=U_3,$$L\left(x,t\right)=L_3$ and $A\left(x,t\right)=A_3$. The solutions of the system tend to ${\mathrm{\Pi }}_3^{\prime }.$ Clearly, ${\mathrm{\Pi }}_3^{\prime }=\left\{{\mathrm{\Xi }}_3\right\};$ hence, LIP confirms the global stability of ${\mathrm{\Xi }}_3$. □

Theorem 4 suggests that if $R_4>1$ and $1<R_3\le 1+\tilde{R}$, then the HTLV-I and SARS-CoV-2 co-infection will be established regardless of the initial states.

The results of Section 3 and Section 4 are summarized in

Table 1. Existence and stability conditions.

Steady State	Existence Conditions	Global Stability Conditions
${\mathrm{\Xi }}_0=\left(X_0,0,0,0,U_0,0,0\right)$	None	$R_1\le 1$ and $R_2\le 1$
${\mathrm{\Xi }}_1=\left(X_1,0,0,0,U_1,L_1,A_1\right)$	$R_1>1$	$R_1>1$ and $R_4\le 1$
${\mathrm{\Xi }}_2=\left(X_2,N_2,Y_2,V_2,U_2,0,0\right)$	$R_2>1$	$R_2>1$ and $R_3\le 1$
${\mathrm{\Xi }}_3=\left(X_3,N_3,Y_3,V_3,U_3,L_3,A_3\right)$	$R_3>1$ and $R_4>1$	$R_4>1$ and $1<R_3\le 1+\tilde{R}$

.

6. Numerical Simulations

This section presents some numerical results for models (12)–(18) to illustrate the stability of the steady states. In addition, we present a comparison of the results on the impact of HTLV-I infection on the dynamics of SARS-CoV-2 single infection. To numerically solve the system of PDEs, we use the solver PDEPE in MATLAB (see the code given in the link given in [66]: https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/math10224390/s1 (accessed on 1 December 2022)).

6.1. Stability of Steady States

We solve the system with the following initial conditions:

$$\begin{array}{cc}\hfill X\left(x,0\right)& =5\left[1+0.2{cos}^2\left(2\pi x\right)\right],\hfill \\ \hfill N\left(x,0\right)& =0.0001\left[1+0.2{cos}^2\left(2\pi x\right)\right],\hfill \\ \hfill Y\left(x,0\right)& =0.0002\left[1+0.2{cos}^2\left(2\pi x\right)\right],\hfill \\ \hfill V\left(x,0\right)& =0.0003\left[1+0.2{cos}^2\left(2\pi x\right)\right],\hfill \\ \hfill U\left(x,0\right)& =500\left[1+0.2{cos}^2\left(2\pi x\right)\right],\hfill \\ \hfill L\left(x,0\right)& =150\left[1+0.2{cos}^2\left(2\pi x\right)\right],\hfill \\ \hfill A\left(x,0\right)& =50\left[1+0.2{cos}^2\left(2\pi x\right)\right],x\in \left[0,1\right],\hfill \end{array}$$

(47)

and homogeneous Neumann boundary conditions

$$\frac{\partial X}{\partial \overrightarrow{\eta }}=\frac{\partial N}{\partial \overrightarrow{\eta }}=\frac{\partial Y}{\partial \overrightarrow{\eta }}=\frac{\partial V}{\partial \overrightarrow{\eta }}=\frac{\partial U}{\partial \overrightarrow{\eta }}=\frac{\partial L}{\partial \overrightarrow{\eta }}=\frac{\partial A}{\partial \overrightarrow{\eta }}=0,t>0,x\in \left[0,1\right].$$

(48)

Note that we did not have real data from HTLV-I and SARS-CoV-2 co-infected patients; therefore, the initial conditions were arbitrarily selected. Our global stability results guarantee that the solutions of the model converge to one of the four steady states regardless of the selected initial conditions.

Table 2. Model parameters.

Parameter	Value	Source	Parameter	Value	Source
$\upsilon$	$0.11$	[28]	$\phi$	$0.9$	[38]
${\beta }_X$	$0.011$	[28]	${\gamma }^{*}$	$0.011$	[38]
$\sigma$	Varies	Assumed	$\lambda$	$0.003$	[7,37,39,40]
$\varrho$	$4.08$	[53,56]	${\beta }_L$	$0.03$	[38,39,40,41]
${\beta }_N$	$0.11$	[28]	${\beta }_A^{*}$	$0.03$	[37,39,40]
${\beta }_Y$	$0.11$	[29,56]	${\delta }_X$	$0.01$	Assumed
$\xi$	Varies	Assumed	${\delta }_N$	$0.01$	Assumed
$\rho$	$0.24$	[29,56]	${\delta }_Y$	$0.02$	Assumed
${\beta }_V$	Varies	Assumed	${\delta }_V$	$0.02$	Assumed
$\varkappa$	10	[39,40,74]	${\delta }_U$	$0.1$	Assumed
$\vartheta$	$0.1$	[56,75]	${\delta }_L$	$0.01$	Assumed
${\beta }_U$	$0.012$	[37,39,40]	${\delta }_A$	$0.2$	Assumed
$\tau$	Varies	Assumed

lists the values of parameters. We selected four sets of the parameters $(\sigma ,\xi ,{\beta }_V,\tau )$ to obtain the following plans.

Plan (I) (stability of${\mathrm{\Xi }}_0$): Choosing $(\sigma ,\xi ,{\beta }_V,\tau )=(1,1.5,6,0.0002)$ gives $R_1=0.5411<1$ and $R_2=0.0003<1.$ Based on Theorem 1, the steady state, ${\mathrm{\Xi }}_0$, is GAS. As shown in Figure 1, the numbers of healthy ECs and T cells tend to their normal values $X_0=10$ and $U_0=833.333$, respectively. Moreover, the other populations will be extinct. This plan leads to the clearance of both HTLV-I and COVID-19. In fact, condition $R_2\le 1$ can be achieved by considering two classes of anti-SARS-CoV-2 drug therapies that reduce the parameters $\rho$ and $\sigma$ by replacing them with $(1-{\omega }_1)\rho$ and $(1-{\omega }_2)\sigma$, respectively. Here, $0\le {\omega }_1\le 1$ and $0\le {\omega }_2\le 1$ are the treatment efficacies of anti-SARS-CoV-2 drug therapies for blocking the infection and production of SARS-CoV-2 particles, respectively [20]. However, $R_1\le 1$ is difficult to achieve because no treatment is currently recommended for individuals with HTLV-1 infection.

Plan (II) (stability of${\mathrm{\Xi }}_1$): We take $(\sigma ,\xi ,{\beta }_V,\tau )=(0.5,1.5,6,0.0025)$. So, we obtain $R_1=6.7641>1$ and $R_4=0.0011<1.$ Lemma 2 and Theorem 2 show the existence and global stability of ${\mathrm{\Xi }}_1$. Figure 2 demonstrates that the numerical solutions reach the steady state ${\mathrm{\Xi }}_1=(10,0,0,0,123.2,266.53,27.67)$. This plan leads the co-infected patient to be infected with HTLV- I only, while the COVID-19 infection disappears. Because $R_4$ depends on the parameters $\rho$ and $\sigma$, $R_4\le 1$ can be achieved by applying treatments with efficacies ${\omega }_1$ and ${\omega }_2$.

Plan (III) (stability of${\mathrm{\Xi }}_2$): We choose $(\sigma ,\xi ,{\beta }_V,\tau )=(5,0.015,0.06,0.0002)$. Then, we calculate $R_2=15.444>1$ and $R_3=0.5769<1$. Based on Lemma 2 and Theorem 3, ${\mathrm{\Xi }}_2$ exists, and it is GAS. The numerical solution plotted in Figure 3 converges to ${\mathrm{\Xi }}_2=\left(0.69,0.024,0.007,0.03,888.28,0,0\right)$. This plan leads the co-infected patient to being infected with SARS-CoV-2 only. This case is practically difficult to be achieved because of the unavailability of HTLV-I treatment.

Plan (IV) (stability of${\mathrm{\Xi }}_3$): We consider $(\sigma ,\xi ,{\beta }_V,\tau )=(5,0.015,0.6,0.0025)$. So, we obtain $R_4=9.9463>1,$$R_3=6.9035>1$ and $R_3<1+\tilde{R}=7.4676$. From Lemma 2 and Theorem 4, we have that the HTLV-I and SARS-CoV-2 co-infection steady state ${\mathrm{\Xi }}_3=\left(1.005,0.024,0.049,0.02,123.2,285.49,29.64\right)$ exists, and it is GAS. The numerical solution of the system converges to ${\mathrm{\Xi }}_3$ (see Figure 4). This plan suggests that the patient will still be co-infected with HTLV-I and SARS-CoV-2. Co-infection with these two diseases for a long time may expose patients to a significant deterioration in their health condition, which may lead to death.

6.2. Comparison of Results

Next, we investigated the impact of HTLV-I infection on the progression of COVID-19 infection by presenting a comparison between the dynamics of SARS-CoV-2 single infection, and HTLV-I and SARS-CoV-2 co-infection. To establish the effect of HTLV-I infection on COVID-19 infection progression, we compared the solutions of system (12)–(18) with those of the following system for SARS-CoV-2 single infection:

$$\begin{array}{cc}\hfill \frac{\partial X\left(x,t\right)}{\partial t}& ={\delta }_X\Delta X\left(x,t\right)+\upsilon -{\beta }_{X}X\left(x,t\right)-\sigma V\left(x,t\right)X\left(x,t\right),\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill \frac{\partial N\left(x,t\right)}{\partial t}& ={\delta }_N\Delta N\left(x,t\right)+\sigma V\left(x,t\right)X\left(x,t\right)-(\varrho +{\beta }_N)N\left(x,t\right),\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill \frac{\partial Y\left(x,t\right)}{\partial t}& ={\delta }_Y\Delta Y\left(x,t\right)+\varrho N\left(x,t\right)-{\beta }_{Y}Y\left(x,t\right),\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill \frac{\partial V\left(x,t\right)}{\partial t}& ={\delta }_V\Delta V\left(x,t\right)+\rho Y\left(x,t\right)-{\beta }_{V}V\left(x,t\right).\hfill \end{array}$$

(52)

We fixed the parameters $(\sigma ,\xi ,{\beta }_V,\tau )=(5,0.015,0.6,0.0025)$ and solved systems (12)–(18) and (49)–(52) with initial conditions (47) and boundary conditions (48). We can see from Figure 5 that the presence of HTLV-I reduces the numbers of healthy ECs and T cells, while it increases the numbers of SARS-CoV-2 particles and SARS-CoV-2 infected ECs. This means that HTLV-I infection weakens the immune response and increases COVID-19 disease progression.

7. Conclusions

Mathematical models are considered important tools that can be used for better understanding the in-host dynamics of human viral co-infections under the impact of the immune response. Cases of HTLV-I and SARS-CoV-2 co-infection were recorded [6,13,14,15]. A mathematical model for HTLV-I and SARS-CoV-2 co-infection was developed and analyzed [58]. The vast majority of this model is based on an insufficient approach where the populations (cells and viruses) are homogeneously distributed in the body. In this study, we proposed and examined an HTLV-I and SARS-CoV-2 co-infection model with diffusion. We determined four threshold numbers, $R_1,R_2,R_3$, and $R_4$, which determine the existence and global stability of the four steady states. We formulated Lyapunov functions and applied LIP to establish the global stability of the steady states. We proved the following: (a) ii $R_1\le 1$ and $R_2\le 1$, then the infection-free steady state is GAS; (b) if $R_1>1$ and $R_4\le 1$, then the HTLV-I single infection steady state is GAS; (c) if $R_2>1$ and $R_3\le 1$, then the SARS-CoV-2 single-infection steady state is GAS; (d) if $R_4>1$ and $1<R_3\le 1+\tilde{R}$, then the HTLV-I and SARS-CoV-2 co-infection steady state is GAS. To validate the theoretical results, we performed numerical simulations for the PDE model. We discussed the impact of HTLV-I infection on COVID-19 progression. We found that the presence of HTLV-I inhibits the immune response and increases the progression of COVID-19 infection. This observation agrees with that of [6], who reported that HTLV-I causes immune dysfunction even in asymptomatic carriers. Therefore, HTLV-I may increase the risk of COVID-19 infection.

In the future, the co-infection model presented in this article can be extended by (i) incorporating the intracellular and immune response time delays, (ii) considering viral mutations, and (iii) incorporating the memory effect by formulating the co-infection model by fractional differential equations.