Global Stability of Delayed SARS-CoV-2 and HTLV-I Coinfection Models within a Host

Abstract

The aim of the present paper is to formulate two new mathematical models to describe the co-dynamics of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and human T-cell lymphotropic virus type-I (HTLV-I) in a host. The models characterizes the interplaying between seven compartments, uninfected ECs, latently SARS-CoV-2-infected ECs, actively SARS-CoV-2-infected ECs, free SARS-CoV-2 particles, uninfected CD4${}^{+}$T cells, latently HTLV-I-infected CD4${}^{+}$T cells and actively HTLV-I-infected CD4${}^{+}$T cells. The models incorporate five intracellular time delays: (i) two delays in the formation of latently SARS-CoV-2-infected ECs and latently HTLV-I-infected CD4${}^{+}$T cells, (ii) two delays in the reactivation of latently SARS-CoV-2-infected ECs and latently HTLV-I-infected CD4${}^{+}$T cells, and (iii) maturation delay of new SARS-CoV-2 virions. We consider discrete-time delays and distributed-time delays in the first and second models, respectively. We first investigate the properties of the model’s solutions, then we calculate all equilibria and study their global stability. The global asymptotic stability is examined by constructing Lyapunov functionals. The analytical findings are supported via numerical simulation. The impact of time delays on the coinfection progression is discussed. We found that, increasing time delays values can have an antiviral treatment-like impact. Our developed coinfection model can contribute to understand the SARS-CoV-2 and HTLV-I co-dynamics and help to select suitable treatment strategies for COVID-19 patients with HTLV-I.

1. Introduction

At the end of 2019, the world witnessed the emergence of a new infectious disease in Wuhan, China, which was called coronavirus disease 2019 (COVID-19). This disease is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Within a few months, SARS-CoV-2 infection spread rapidly in most countries of the world and infected many people. According to the update provided by the World Health Organization (WHO) on 23 October 2022 [1], over 624 million confirmed cases and over 6.5 million deaths have been reported globally. SARS-CoV-2 is transmitted to people when they are exposed to respiratory fluids carrying infectious viral particles. The implementation of preventive measures such as hand washing, using of face masks, physical and social distancing, disinfection of surfaces and getting COVID-19 vaccine can reduce the SARS-CoV-2 transmission. WHO approved eleven vaccines for COVID-19 for emergency use including Novavax, CanSino, Bharat Biotech, Pfizer/BioNTech, Moderna, Serum Institute of India (Novavax formulation), Janssen (Johnson & Johnson), Oxford/AstraZeneca, Serum Institute of India (Oxford/AstraZeneca formulation), Sinopharm, and Sinovac [2].

SARS-CoV-2 is a single-stranded positive-sense RNA virus and infects the ECs of the respiratory tracts. The virus has clinical manifestations including dyspnea, cough, fever, headache, rhinitis, myalgia and sore throat. SARS-CoV-2 can lead to acute respiratory distress syndrome (ARDS), which has high mortality rates, particularly in patients with immunosenescence [3] susceptibility to viral infections [4]. It was reported in [5] that, 94.2% of patients with COVID-19 were also coinfected with several other microorganisms, such as fungi, bacteria and viruses. Important viral copathogens include the respiratory syncytial virus, rhinovirus/enterovirus, influenza A and B viruses (IAV and IBV), metapneumovirus, parainfluenza virus, human immunodeficiency virus (HIV), cytomegalovirus (CMV), dengue virus (DENV), hepatitis B virus (HBV) and Epstein Barr virus (EBV) (see the review paper [6]).

Human T-cell lymphotropic virus type-I (HTLV-I) is a single-stranded RNA virus which infects the essential human system immune cells, CD4${}^{+}$T cells. Therefore, HTLV-I can cause immune dysfunction even in asymptomatic carriers [3]. HTLV-I can leads to two diseases, adult T-cell leukemia (ATL) and HTLV-I-associated myelopathy/tropical spastic paraparesis (HAM/TSP) [7]. Although HTLV-I can cause fatal diseases (ATL and HAM/TSP), most of infected persons remain asymptomatic throughout their lives [3]. An estimation by WHO told that about 5 to 10 million individuals are infected with HTLV-I worldwide [8]. The primary way for HTLV-I transmission is through bodily fluids including: semen, blood and breast milk [9]. The spread of this viral infection that causes significant morbidity and mortality in some countries of the world has increased the interest and awareness of medical and biological scientists.

In [3,10], two cases of COVID-19 patients with HTLV-I infection were reported. These reports highlighted the need for the accumulation of similar cases to illustrate the risk factors for severe illness, the best-in-class antiviral agent, how to manage and prevent secondary infection, and the optimal treatment strategy for patients with SARS-CoV-2 and HTLV-I coinfection. Sajjadi et al. [11] presented a review on the pathogenesis of both SARS-CoV-2 and HTLV-I infections and discussed their similarities in triggering immune responses.

Mathematical Models of Within-Host HTLV-I and SARS-CoV-2 Mono-Infections

During the past decade, mathematical models which describe the within-host dynamics of human viruses have demonstrated their ability to provide useful insight to gain a further understanding the dynamics and mechanisms of the viruses. These models may assist in the development of viral therapies and vaccines as well as the selection of appropriate therapeutic and vaccine strategies. Further, these models are helpful in determine the sufficient number of factors to analyze the experimental results and explain the biological phenomena [12]. Furthermore, mathematical models have been a key in studying viral-viral coinfections and can be particularly useful in examining the interplay between respiratory viruses and chronic viruses.

Hernandez-Vargas and Velasco-Hernandez [13] formulated the basic target cell-limited model for SARS-CoV-2 dynamics within a host. The model contained three compartments, uninfected ECs (X), active SARS-CoV-2-infected ECs (Y) and free SARS-CoV-2 particles (V). The model was formulated as a system of ODEs:

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

(1)

$$\begin{array}{cc}\hfill \dot{Y}\left(t\right)& =\stackrel{\mathrm{SARS}-\mathrm{CoV}-2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\rho V\left(t\right)X\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\xi }_{Y}Y\left(t\right)}},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \dot{V}\left(t\right)& =\stackrel{\mathrm{generation}\mathrm{of}\mathrm{SARS}-\mathrm{CoV}-2}{\overbrace{\eta Y\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\xi }_{V}V\left(t\right)}}.\hfill \end{array}$$

(3)

In the same paper [13], this model was extended to include the eclipse (latent) phase. The eclipse phase is defined as the period of time that elapses between the entry of the virus into the uninfected epithelial cell and the release of virions from that infected cell [14]. There are two methods to include the latent phase in the viral infection models. The first method is to comprise the eclipse phase by introducing it as a time delay in the model [15]. This can be conducted by reformulating Equation (2) as [15]:

$$\dot{Y}\left(t\right)=\rho V(t-\tau )X(t-\tau )-{\xi }_{Y}Y\left(t\right),$$

where $\tau$ is the delay between viral entry into an epithelial cell and the start of production IAV virions. The second method is to define two separate populations of SARS-CoV-2-infected ECs: one population is the latently SARS-CoV-2-infected cells (N) which contains the SARS-CoV-2 virions but not yet producing it; the second population is the actively SARS-CoV-2-infected cells (Y) which produce the SARS-CoV-2 [16]. SARS-CoV-2 infection model with latent infected cells was formulated as [13]:

$$\begin{array}{cc}\hfill \dot{X}\left(t\right)& =-\stackrel{\mathrm{SARS}-\mathrm{CoV}-2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\rho V\left(t\right)X\left(t\right)}},\hfill \\ \hfill \dot{N}\left(t\right)& =\stackrel{\mathrm{SARS}-\mathrm{CoV}-2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\rho V\left(t\right)X\left(t\right)}}-\stackrel{\mathrm{latent}\mathrm{activation}}{\overbrace{\kappa N\left(t\right)}},\hfill \\ \hfill \dot{Y}\left(t\right)& =\stackrel{\mathrm{latent}\mathrm{activation}}{\overbrace{\kappa N\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\xi }_{Y}Y\left(t\right)}},\hfill \\ \hfill \dot{V}\left(t\right)& =\stackrel{\mathrm{generation}\mathrm{of}\mathrm{SARS}-\mathrm{CoV}-2}{\overbrace{\eta Y\left(t\right)}}-\stackrel{\mathrm{death}}{\overbrace{{\xi }_{V}V\left(t\right)}}.\hfill \end{array}$$

Li et al. [17] have considered a SARS-CoV-2 infection model with regeneration and death for the uninfected ECs as:

$$\dot{X}=\delta -{\xi }_{X}X-\rho VX,$$

where $\delta$ and ${\xi }_{X}X$ are the regeneration and death rates of uninfected ECs, respectively. Models presented in [13,17] were extended and modified by including (i) the effect of immune response [18,19,20,21,22,23], (ii) the influence of different drug therapies [24,25], and (iii) the impact of time delay [26].

In very recent works, mathematical models have been formulated to describe the coinfection of COVID-19 with other diseases in epidemiology such as: COVID-19/Dengue [27], COVID-19/Influenza [28], COVID-19/HIV [29], COVID-19/ZIKV [30], COVID-19/Dengue/ HIV [31], COVID-19/Tuberculosis [32], COVID-19/Bacterial [33]. However, modeling of within-host dynamics of SARS-CoV-2 with other pathogens coinfection has been investigated in a few papers: SARS-CoV-2/Bacteria [34], SARS-CoV-2/HIV [35], SARS-CoV-2/malaria [36] and SARS-CoV-2/Influenza A virus [37].

Stability analysis of viral infection models can help researchers to (i) expect the qualitative features of the model for a given set of values of the model’s parameters, (ii) establish the conditions that ensure the persistence or deletion of this infection, and (iii) determine under what conditions the immune system is stimulated against the infection. Stability of time-delay systems has received great interest in different fields (see e.g., [38,39,40,41]). Stability analysis for models describing the within-host dynamics of SARS-CoV-2 infection was studied in [21,22,23,35,37,42]. Hattaf and Yousfi [21] studied a within-host SARS-CoV-2 infection model with cell-to-cell transmission and Cytotoxic T lymphocyte (CTL) immune response. The model included both lytic and nonlytic immune responses. Lyapunov method was used to prove the global stability of the three equilibria of the model. A SARS-CoV-2 infection model with both CTL and antibody immunities was developed and analyzed in [23]. Mathematical analysis of the model presented in [17] was studied in [42]. Both local and global stability analysis of the model’s equilibria were established. Almocera et al. [22] studied the stability of the two-dimensional SARS-CoV-2 dynamics model with immune response presented in [13]. Elaiw et al. [26] studied the global stability of a delayed SARS-CoV-2 dynamics model with logistic growth of the uninfected ECs and antibody immunity. In very recent works, Lyapunov method was used to establish the global stability of coinfection models including, SARS-CoV-2/HIV-1 [35], SARS-CoV-2/Influenza A virus [37] and SARS-CoV-2/malaria [36].

Modeling and analysis of HTLV-I mono-infection have attracted the interest of several mathematician. HTLV-I mono-infection model which describes the interaction of three compartments, healthy (or uninfected) CD4${}^{+}$T cells (U), latently HTLV-I-infected cells (L) and actively HTLV-I-infected cells (A), is given by [43]:

$$\begin{array}{cc}\hfill \dot{U}& =\stackrel{{\mathrm{CD}4}^{+}\mathrm{T}\mathrm{cells}\mathrm{production}}{\overbrace{\gamma }}-\stackrel{\mathrm{death}}{\overbrace{{\xi }_{U}U}}-\stackrel{\mathrm{HTLV}-\mathrm{I}\mathrm{infectious}\mathrm{transmission}}{\overbrace{\pi AU}},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \dot{L}& =\stackrel{\mathrm{HTLV}-\mathrm{I}\mathrm{infectious}\mathrm{transmission}}{\overbrace{\pi AU}}-\stackrel{\mathrm{latent}\mathrm{activation}}{\overbrace{\alpha L}}-\stackrel{\mathrm{death}}{\overbrace{{\xi }_{L}L}},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \dot{A}& =\stackrel{\mathrm{latent}\mathrm{activation}}{\overbrace{\alpha L}}-\stackrel{\mathrm{death}}{\overbrace{{\xi }_A^{*}A}}.\hfill \end{array}$$

(6)

In [44], model (4)–(6) was extended by incorporating a time delay $\tau$ which accounts the time between initial infection of an uninfected target cells to become latently infected cells. Equation (5) was replaced by

$$\dot{L}\left(t\right)=\pi A\left(t-\tau \right)U\left(t-\tau \right)-\left(\alpha +{\xi }_L\right)L\left(t\right).$$

Different biological factors were considered in HTLV-I infection models including (i) Cytotoxic T-Lymphocytes (CTLs) immunity [7,45,46,47,48], (ii) mitotic transmission of actively infected cells [49,50,51,52,53], (iii) intracellular time delay [44,47,48]) or immune response delay [45,46,54], and (iv) reaction and diffusion [55]. HTLV-I has quite similar ways of transmissions as HIV-1. Therefore, we formulated and analyzed some HIV-1/HTLV-I coinfection models [56].

To the best of our knowledge, mathematical modeling of within-host SARS-CoV-2 and HTLV-I coinfection has not been studied before. The objective of this work is to formulate new models for within-host SARS-CoV-2-HTLV-I coinfection with discrete/distributed time delays. We study the properties of the model’s solutions, calculate all equilibrium points, investigate the global stability of equilibria and conduct some numerical simulations. Finally, we discuss the obtained results.

The SARS-CoV-2 and HTLV-I coinfection models presented in this paper can be helpful to describe the co-dynamics of several human viruses. In addition, the models may be used to predict new treatment regimens and strategies for patients who are coinfected with different viruses or multi variants of a virus [57].

2. SARS-CoV-2 and HTLV-I Coinfection Model with Discrete-Time Delays

In this section, we formulate a SARS-CoV-2 and HTLV-I coinfection dynamics model. Let us consider following hypothesis:

Hypothesis 1 (H1).

The model describes the interactions between seven compartments, uninfected ECs (X), latently SARS-CoV-2-infected cells (N), actively SARS-CoV-2-infected cells (Y), free SARS-CoV-2 particles (V), uninfected CD4${}^{+}$T cells (U), latently HTLV-I-infected cells (L) and actively HTLV-I-infected cells (A).

Hypothesis 2 (H2).

The uninfected ECs and CD4${}^{+}$T cells are the targets for SARS-CoV-2 and HTLV-I, respectively, [13,43].

Hypothesis 3 (H3).

The CD4${}^{+}$T cells help the CTLs to kill the actively SARS-CoV-2-infected ECs at rate $\mu YU$ and are proliferated at rate $uYU$ [35].

Hypothesis 4 (H4).

The actively HTLV-I-infected cells are proliferated at rate ${\epsilon }^{*}A$, with a part $\omega {\epsilon }^{*}A$ turn into latent, while the other part $(1-\omega ){\epsilon }^{*}A$ remains active, where $\omega \in (0,1)$ [52,56].

Hypothesis 5 (H5).

There exist two time delays ${\tau }_1$ and ${\tau }_4$ in the formation of latently SARS-CoV-2-infected cells and latently HTLV-I-infected, respectively [26,58,59].

Hypothesis 6 (H6).

There exist two reactivation time delays ${\tau }_2$ and ${\tau }_5$ of the latently SARS-CoV-2-infected cells and latently HTLV-I-infected cells, respectively [26,58,59].

Hypothesis 7 (H7).

There exist a maturation time delay ${\tau }_3$ of SARS-CoV-2 virions [26].

Under hypothesis H1–H7 we propose the model for SARS-CoV-2 and HTLV-I coinfection within a host as:

$$\begin{array}{cc}\hfill \dot{X}\left(t\right)& =\delta -{\xi }_{X}X\left(t\right)-\rho V\left(t\right)X\left(t\right),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \dot{N}\left(t\right)& =\rho e^{-d_1{\tau }_1}V(t-{\tau }_1)X\left(t-{\tau }_1\right)-\left(\kappa +{\xi }_N\right)N\left(t\right),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \dot{Y}\left(t\right)& =\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)-{\xi }_{Y}Y\left(t\right)-\mu Y\left(t\right)U\left(t\right),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \dot{V}\left(t\right)& =\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)-{\xi }_{V}V\left(t\right),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \dot{U}\left(t\right)& =\gamma +uY\left(t\right)U\left(t\right)-{\xi }_{U}U\left(t\right)-\pi A\left(t\right)U\left(t\right),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \dot{L}\left(t\right)& =\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\omega {\epsilon }^{*}A\left(t\right)-\left(\alpha +{\xi }_L\right)L\left(t\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \dot{A}\left(t\right)& =\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)+(1-\omega ){\epsilon }^{*}A\left(t\right)-{\xi }_A^{*}A\left(t\right).\hfill \end{array}$$

(13)

The factors $e^{-d_1{\tau }_1},e^{-d_2{\tau }_2},e^{-d_3{\tau }_3},e^{-d_4{\tau }_4}$ and $e^{-d_5{\tau }_5}$ denote the probability of surviving the compartments during the time delay periods, where, $d_1,d_2,d_3,d_4$ and $d_5$ are constants. Let ${\tau }^{*}=max\left\{{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4,{\tau }_5\right\}$ and assume the initial conditions of system (7)–(13) having the following form:

$$\begin{array}{cc}\hfill X\left(\theta \right)& ={\varphi }_1\left(\theta \right),N\left(\theta \right)={\varphi }_2\left(\theta \right),Y\left(\theta \right)={\varphi }_3\left(\theta \right),V\left(\theta \right)={\varphi }_4\left(\theta \right),\hfill \\ \hfill U\left(\theta \right)& ={\varphi }_5\left(\theta \right),L\left(\theta \right)={\varphi }_6\left(\theta \right),A\left(\theta \right)={\varphi }_7\left(\theta \right),\hfill \\ \hfill {\varphi }_i\left(\theta \right)& \ge 0,\theta \in \left[-{\tau }^{*},0\right],i=1,2,\dots ,7,\hfill \end{array}$$

(14)

where ${\varphi }_i\in$$\mathbb{C}$ and $\mathbb{C}=\mathbb{C}\left(\left[-{\tau }^{*},0\right],{\mathbb{R}}_{\ge 0}\right)$ is the Banach space of continuous mapping the interval $\left[-{\tau }^{*},0\right]$ into ${\mathbb{R}}_{\ge 0}$ with the norm $∥{\varphi }_i∥=\underset{-{\tau }^{*}\le \theta \le 0}{sup}\left|{\varphi }_i\left(\theta \right)\right|$ for ${\varphi }_i\in \mathbb{C}$, $i=1,2,\dots ,7$. By the fundamental theory of functional differential equations, system (7)–(13) with initial conditions (14) has a unique solution [60]. We note that, model (7)–(13) can be seen as a generalization of many models for HTLV-I and SARS-CoV-2 mono-infections presented in the literature.

All parameters of model (7)–(13) are positive. In [52,56], it was proposed that ${\epsilon }^{*}<min\{{\xi }_U,{\xi }_L,{\xi }_A^{*}\}$. Since $0<\omega <1$ and ${\epsilon }^{*}<{\xi }_A^{*}$, thus $(1-\omega ){\epsilon }^{*}<{\xi }_A^{*}$. Denote ${\xi }_A={\xi }_A^{*}-(1-\omega ){\epsilon }^{*}>0$ and $\epsilon =\omega {\epsilon }^{*}.$ We have ${\xi }_A-\epsilon ={\xi }_A^{*}-{\epsilon }^{*}>0$. Therefore, model (7)–(13) becomes:

$$\begin{array}{cc}\hfill \dot{X}\left(t\right)& =\delta -{\xi }_{X}X\left(t\right)-\rho V\left(t\right)X\left(t\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \dot{N}\left(t\right)& =\rho e^{-d_1{\tau }_1}V(t-{\tau }_1)X\left(t-{\tau }_1\right)-\left(\kappa +{\xi }_N\right)N\left(t\right),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \dot{Y}\left(t\right)& =\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)-{\xi }_{Y}Y\left(t\right)-\mu Y\left(t\right)U\left(t\right),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \dot{V}\left(t\right)& =\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)-{\xi }_{V}V\left(t\right),\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \dot{U}\left(t\right)& =\gamma +uY\left(t\right)U\left(t\right)-{\xi }_{U}U\left(t\right)-\pi A\left(t\right)U\left(t\right),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \dot{L}\left(t\right)& =\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\epsilon A\left(t\right)-\left(\alpha +{\xi }_L\right)L\left(t\right),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \dot{A}\left(t\right)& =\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)-{\xi }_{A}A\left(t\right).\hfill \end{array}$$

(21)

Next, we will examine the nonnegativity and boundedness of the system’s solutions.

2.1. Properties of Solutions

Lemma 1.

The solutions of model (15)–(21) with initial conditions (14) are nonnegative and ultimately bounded.

Proof. 

From Equations (15) and (19), we obtain $\dot{X}{\mid }_{X=0}=\delta >0$ and $\dot{U}{\mid }_{U=0}=\gamma >0$. It follows that $X\left(t\right)>0$ and $U\left(t\right)>0$, for all $t\ge 0$. From Equations (16)–(18), (20) and (21) we have

$$\begin{array}{cc}\hfill N\left(t\right)& =e^{-\left(\kappa +{\xi }_N\right)t}{\varphi }_2\left(0\right)+\rho e^{-d_1{\tau }_1}\underset{0}{\overset{t}{\int }}{e}^{-\left(\kappa +{\xi }_N\right)(t-\theta )}V(\theta -{\tau }_1)X\left(\theta -{\tau }_1\right)d\theta \ge 0,\hfill \\ \hfill Y\left(t\right)& =e^{-\underset{0}{\overset{t}{\int }}\left({\xi }_Y+\mu U\left(\varkappa \right)\right)d\varkappa }{\varphi }_3\left(0\right)+\kappa e^{-d_2{\tau }_2}\underset{0}{\overset{t}{\int }}{e}^{-\underset{\theta }{\overset{t}{\int }}\left({\xi }_Y+\mu U\left(\varkappa \right)\right)d\varkappa }N(\theta -{\tau }_2)d\theta \ge 0,\hfill \\ \hfill V\left(t\right)& =e^{-{\xi }_{V}t}{\varphi }_4\left(0\right)+\eta e^{-d_3{\tau }_3}\underset{0}{\overset{t}{\int }}{e}^{-{\xi }_V(t-\theta )}Y\left(\theta -{\tau }_3\right)d\theta \ge 0,\hfill \\ \hfill L\left(t\right)& =e^{-\left(\alpha +{\xi }_L\right)t}{\varphi }_6\left(0\right)+\underset{0}{\overset{t}{\int }}{e}^{-\left(\alpha +{\xi }_L\right)(t-\theta )}\left[\pi e^{-d_4{\tau }_4}A\left(\theta -{\tau }_4\right)U\left(\theta -{\tau }_4\right)+\epsilon A\left(\theta \right)\right]d\theta \ge 0,\hfill \\ \hfill A\left(t\right)& =e^{-{\xi }_{A}t}{\varphi }_7\left(0\right)+\alpha e^{-d_5{\tau }_5}\underset{0}{\overset{t}{\int }}{e}^{-{\xi }_A(t-\theta )}L\left(\theta -{\tau }_5\right)d\theta \ge 0,\hfill \end{array}$$

for all $t\in \left[0,{\tau }^{*}\right].$ Hence, by a recursive argument, we obtain $N\left(t\right),Y\left(t\right),V\left(t\right),L\left(t\right)$,$A\left(t\right)\ge 0$ for all $t\ge 0.$ This guarantees that $\left(X\left(t\right),N\left(t\right),Y\left(t\right),V\left(t\right),U\left(t\right),L\left(t\right),A\left(t\right)\right)\in {\mathbb{R}}_{\ge 0}^7$ for all $t\ge 0$ if $\left(X\left(0\right),N\left(0\right),Y\left(0\right),V\left(0\right),U\left(0\right),L\left(0\right),A\left(0\right)\right)\in {\mathbb{R}}_{\ge 0}^7.$ To investigate the boundedness of the model’s solutions, from Equation (15), we have $\underset{t\to \infty }{lim\; sup}X\left(t\right)\le \frac{\delta }{{\xi }_X}=M_1$. We define

$$\begin{array}{cc}\hfill \Theta \left(t\right)& =e^{-d_1{\tau }_1}X\left(t-{\tau }_1\right)+N\left(t\right)+Y\left(t\right)+\kappa \underset{t-{\tau }_2}{\overset{t}{\int }}{e}^{-d_2\left(t-\theta \right)}N\left(\theta \right)d\theta +\frac{{\xi }_Y}{2\eta }V\left(t\right)+\frac{{\xi }_Y}{2}\underset{t-{\tau }_3}{\overset{t}{\int }}{e}^{-d_3\left(t-\theta \right)}Y\left(\theta \right)d\theta \hfill \\ \hfill & +\frac{\mu }{u}\left[U\left(t\right)+L\left(t\right)+A\left(t\right)+\pi \underset{t-{\tau }_4}{\overset{t}{\int }}{e}^{-d_4\left(t-\theta \right)}A\left(\theta \right)U\left(\theta \right)d\theta +\alpha \underset{t-{\tau }_5}{\overset{t}{\int }}{e}^{-d_5\left(t-\theta \right)}L\left(\theta \right)d\theta \right].\hfill \end{array}$$

Then

$$\begin{array}{cc}\hfill \dot{\Theta }\left(t\right)& =e^{-d_1{\tau }_1}\dot{X}\left(t-{\tau }_1\right)+\dot{N}\left(t\right)+\dot{Y}\left(t\right)-d_2\kappa \underset{t-{\tau }_2}{\overset{t}{\int }}{e}^{-d_2\left(t-\theta \right)}N\left(\theta \right)d\theta +\kappa N\left(t\right)-\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)\hfill \\ \hfill & +\frac{{\xi }_Y}{2\eta }\dot{V}\left(t\right)-d_3\frac{{\xi }_Y}{2}\underset{t-{\tau }_3}{\overset{t}{\int }}{e}^{-d_3\left(t-\theta \right)}Y\left(\theta \right)d\theta +\frac{{\xi }_Y}{2}Y\left(t\right)-\frac{{\xi }_Y}{2}{e}^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)\hfill \\ \hfill & +\frac{\mu }{u}\left[\dot{U}\left(t\right)+\dot{L}\left(t\right)+\dot{A}\left(t\right)-d_4\pi \underset{t-{\tau }_4}{\overset{t}{\int }}{e}^{-d_4\left(t-\theta \right)}A\left(\theta \right)U\left(\theta \right)d\theta +\pi A\left(t\right)U\left(t\right)\right.\hfill \\ \hfill & \left.-\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)-d_5\alpha \underset{t-{\tau }_5}{\overset{t}{\int }}{e}^{-d_5\left(t-\theta \right)}L\left(\theta \right)d\theta +\alpha L\left(t\right)-\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)\right]\hfill \\ \hfill & =e^{-d_1{\tau }_1}\left(\delta -{\xi }_{X}X\left(t-{\tau }_1\right)-\rho V(t-{\tau }_1)X\left(t-{\tau }_1\right)\right)+\rho e^{-d_1{\tau }_1}V(t-{\tau }_1)X\left(t-{\tau }_1\right)\hfill \\ \hfill & -\left(\kappa +{\xi }_N\right)N\left(t\right)+\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)-{\xi }_{Y}Y\left(t\right)-\mu Y\left(t\right)U\left(t\right)-d_2\kappa \underset{t-{\tau }_2}{\overset{t}{\int }}{e}^{-d_2\left(t-\theta \right)}N\left(\theta \right)d\theta \hfill \\ \hfill & +\kappa N\left(t\right)-\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)+\frac{{\xi }_Y}{2\eta }\left(\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)-{\xi }_{V}V\left(t\right)\right)-d_3\frac{{\xi }_Y}{2}\underset{t-{\tau }_3}{\overset{t}{\int }}{e}^{-d_3\left(t-\theta \right)}Y\left(\theta \right)d\theta \hfill \\ \hfill & +\frac{{\xi }_Y}{2}Y\left(t\right)-\frac{{\xi }_Y}{2}{e}^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)+\frac{\mu }{u}\left[\gamma +uY\left(t\right)U\left(t\right)-{\xi }_{U}U\left(t\right)-\pi A\left(t\right)U\left(t\right)\right.\hfill \\ \hfill & +\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\epsilon A\left(t\right)-\left(\alpha +{\xi }_L\right)L\left(t\right)+\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)-{\xi }_{A}A\left(t\right)\hfill \\ \hfill & -d_4\pi \underset{t-{\tau }_4}{\overset{t}{\int }}{e}^{-d_4\left(t-\theta \right)}A\left(\theta \right)U\left(\theta \right)d\theta +\pi A\left(t\right)U\left(t\right)-\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)\hfill \\ \hfill & \left.-d_5\alpha \underset{t-{\tau }_5}{\overset{t}{\int }}{e}^{-d_5\left(t-\theta \right)}L\left(\theta \right)d\theta +\alpha L\left(t\right)-\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)\right]\hfill \\ \hfill & =e^{-d_1{\tau }_1}\delta +\frac{\mu }{u}\gamma -{\xi }_X{e}^{-d_1{\tau }_1}X\left(t-{\tau }_1\right)-{\xi }_{N}N\left(t\right)-\frac{{\xi }_Y}{2}Y\left(t\right)-d_2\kappa \underset{t-{\tau }_2}{\overset{t}{\int }}{e}^{-d_2\left(t-\theta \right)}N\left(\theta \right)d\theta \hfill \\ \hfill & -\frac{{\xi }_Y{\xi }_V}{2\eta }V\left(t\right)-d_3\frac{{\xi }_Y}{2}\underset{t-{\tau }_3}{\overset{t}{\int }}{e}^{-d_3\left(t-\theta \right)}Y\left(\theta \right)d\theta -\frac{\mu }{u}\left[{\xi }_{U}U\left(t\right)+{\xi }_{L}L\left(t\right)+\left({\xi }_A-\epsilon \right)A\left(t\right)\right.\hfill \\ \hfill & \left.+d_4\pi \underset{t-{\tau }_4}{\overset{t}{\int }}{e}^{-d_4\left(t-\theta \right)}A\left(\theta \right)U\left(\theta \right)d\theta +d_5\alpha \underset{t-{\tau }_5}{\overset{t}{\int }}{e}^{-d_5\left(t-\theta \right)}L\left(\theta \right)d\theta \right].\hfill \end{array}$$

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

$$\begin{array}{cc}\hfill \dot{\Theta }\left(t\right)& \le \delta +\frac{\mu }{u}\gamma -\sigma \left[e^{-d_1{\tau }_1}X\left(t-{\tau }_1\right)+N\left(t\right)+Y\left(t\right)+\kappa \underset{t-{\tau }_2}{\overset{t}{\int }}{e}^{-d_2\left(t-\theta \right)}N\left(\theta \right)d\theta +\frac{{\xi }_Y}{2\eta }V\left(t\right)\right.\hfill \\ \hfill & +\frac{{\xi }_Y}{2}\underset{t-{\tau }_3}{\overset{t}{\int }}{e}^{-d_3\left(t-\theta \right)}Y\left(\theta \right)d\theta +\frac{\mu }{u}\left(U\left(t\right)+L\left(t\right)+A\left(t\right)+\pi \underset{t-{\tau }_4}{\overset{t}{\int }}{e}^{-d_4\left(t-\theta \right)}A\left(\theta \right)U\left(\theta \right)d\theta \right.\hfill \\ \hfill & \left.\left.+\alpha \underset{t-{\tau }_5}{\overset{t}{\int }}{e}^{-d_5\left(t-\theta \right)}L\left(\theta \right)d\theta \right)\right]=\delta +\frac{\mu }{u}\gamma -\sigma \Theta \left(t\right),\hfill \end{array}$$

where $\sigma =min\left\{{\xi }_X,{\xi }_N,\frac{{\xi }_Y}{2},d_2,{\xi }_V,d_3,{\xi }_U,{\xi }_L,{\xi }_A^{*}-{\epsilon }^{*},d_4,d_5\right\}$. Hence, $\underset{t⟶\infty }{lim\; sup}\Theta \left(t\right)\le M_2$ for $t\ge 0,$ where $M_2=\frac{\delta +\frac{\mu }{u}\gamma }{\sigma }.$ Since $X,N,Y,V,U,L$ and A are all non-negative, then $\underset{t⟶\infty }{lim\; sup}N\left(t\right)\le M_2,$ $\underset{t⟶\infty }{lim\; sup}Y\left(t\right)\le M_2,$ $\underset{t⟶\infty }{lim\; sup}V\left(t\right)\le M_3,\underset{t⟶\infty }{lim\; sup}U\left(t\right)\le M_4,$$\underset{t⟶\infty }{lim\; sup}L\left(t\right)\le M_4$ and $\underset{t⟶\infty }{lim\; sup}A\left(t\right)\le M_4$, where $M_3=\frac{2\eta }{{\xi }_Y}{M}_2$ and $M_4=\frac{u}{\mu }{M}_2$. Consequently, $X\left(t\right)$, $N\left(t\right)$, $Y\left(t\right),V\left(t\right),U\left(t\right),L\left(t\right)$ and $A\left(t\right)$ are all ultimately bounded. □

According to Lemma 1, it can be shown that the set $\Omega =\{\left(X,N,Y,V,U,L,A\right)\in {\mathbb{C}}_{+}^7:∥X∥\le M_1,$$∥N∥\le M_2,$$∥Y∥\le M_2,$$∥V∥\le M_3,$$∥U∥\le M_4,$$∥L∥\le M_4,$$∥A∥\le M_4\}$ is positively invariant for system (15)–(21).

2.2. Equilibrium Points

To calculate the equilibrium points of the system (15)–(21) we solve the following system:

$$\begin{array}{cc}\hfill 0& =\delta -{\xi }_{X}X-\rho VX,\hfill \\ \hfill 0& =\rho e^{-d_1{\tau }_1}VX-\left(\kappa +{\xi }_N\right)N,\hfill \\ \hfill 0& =\kappa e^{-d_2{\tau }_2}N-{\xi }_{Y}Y-\mu YU,\hfill \\ \hfill 0& =\eta e^{-d_3{\tau }_3}Y-{\xi }_{V}V,\hfill \\ \hfill 0& =\gamma +uYU-{\xi }_{U}U-\pi AU,\hfill \\ \hfill 0& =\pi e^{-d_4{\tau }_4}AU+\epsilon A-\left(\alpha +{\xi }_L\right)L,\hfill \\ \hfill 0& =\alpha e^{-d_5{\tau }_5}L-{\xi }_{A}A.\hfill \end{array}$$

We find that system admits four equilibrium points.

(i) Uninfected equilibrium point, $E{P}_0=\left(X_0,0,0,0,U_0,0,0\right),$ where $X_0=\frac{\delta }{{\xi }_X}$ and $U_0=\frac{\gamma }{{\xi }_U}$.

(ii) HTLV-I mono-infection equilibrium point, $E{P}_1=\left(X_1,0,0,0,U_1,L_1,A_1\right),$ where

$$X_1=\frac{\delta }{{\xi }_X}=X_0,U_1=\frac{U_0}{R_1},L_1=\frac{{\xi }_U{\xi }_A}{\pi \alpha e^{-d_5{\tau }_5}}\left(R_1-1\right),A_1=\frac{{\xi }_U}{\pi }\left(R_1-1\right),$$

where $R_1=\frac{\pi \alpha \gamma e^{-d_4{\tau }_4-d_5{\tau }_5}}{{\xi }_U\left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)\right)}$. Here, $R_1$ is the basic reproduction number of HTLV-I mono-infection. It determines the establishment of HTLV-I infection. Therefore $E{P}_1$ exists when $R_1>1.$

(iii) SARS-CoV-2-mono-infection equilibrium point, $E{P}_2=\left(X_2,N_2,Y_2,V_2,U_2,0,0\right),$ where

$$Y_2=\frac{{\xi }_V}{\eta e^{-d_3{\tau }_3}}{V}_2,U_2=\frac{\eta \gamma e^{-d_3{\tau }_3}}{\eta {\xi }_U{e}^{-d_3{\tau }_3}-u{\xi }_V{V}_2},X_2=\frac{\left(\kappa +{\xi }_N\right)\left({\xi }_Y{Y}_2+\mu U_2{Y}_2\right)}{\kappa \rho e^{-d_1{\tau }_1-d_2{\tau }_2}{V}_2},N_2=\frac{{\xi }_Y{Y}_2+\mu U_2{Y}_2}{\kappa e^{-d_2{\tau }_2}},$$

(22)

and $V_2$ satisfies the following equation:

$$\frac{T_1{V}_2^2+T_2{V}_2+T_3}{\eta \rho \kappa e^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}\left(\eta {\xi }_U{e}^{-d_3{\tau }_3}-u{\xi }_V{V}_2\right)}=0,$$

(23)

where

$$\begin{array}{cc}\hfill T_1& =\rho u{\xi }_Y{\xi }_V^2\left(\kappa +{\xi }_N\right),\hfill \\ \hfill T_2& ={\xi }_V\left(u{\xi }_X{\xi }_Y{\xi }_V\left(\kappa +{\xi }_N\right)-\eta \rho {\xi }_Y{\xi }_U{e}^{-d_3{\tau }_3}\left(\kappa +{\xi }_N\right)-\eta \rho \gamma \mu e^{-d_3{\tau }_3}\left(\kappa +{\xi }_N\right)-\delta \eta \rho \kappa u{e}^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}\right),\hfill \\ \hfill T_3& =e^{-d_3{\tau }_3}\left(\delta {\eta }^2\rho \kappa {\xi }_U{e}^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}-\eta {\xi }_X{\xi }_Y{\xi }_V{\xi }_U\left(\kappa +{\xi }_N\right)-\eta \gamma \mu {\xi }_X{\xi }_V\left(\kappa +{\xi }_N\right)\right).\hfill \end{array}$$

We want to prove that Equation (23) has a positive root. Define a function $F\left(V\right)$ as

$$F\left(V\right)=\frac{T_1{V}^2+T_{2}V+T_3}{\eta \rho \kappa e^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}\left(\eta {\xi }_U{e}^{-d_3{\tau }_3}-u{\xi }_{V}V\right)}.$$

We have

$$\begin{array}{cc}\hfill F\left(0\right)& =\frac{\delta \eta \rho \kappa {\xi }_U{e}^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}-{\xi }_X{\xi }_Y{\xi }_V{\xi }_U\left(\kappa +{\xi }_N\right)-\gamma \mu {\xi }_X{\xi }_V\left(\kappa +{\xi }_N\right)}{\eta \rho \kappa {\xi }_U{e}^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}}\hfill \\ \hfill & =\frac{{\xi }_X{\xi }_V\left(\kappa +{\xi }_N\right)\left({\xi }_Y{\xi }_U+\gamma \mu \right)}{{\xi }_U\eta \rho \kappa e^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}}\left(R_2-1\right),\hfill \end{array}$$

where

$$R_2=\frac{\delta \eta \rho \kappa {\xi }_U{e}^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}}{{\xi }_X{\xi }_V\left(\kappa +{\xi }_N\right)\left({\xi }_Y{\xi }_U+\gamma \mu \right)}.$$

This implies that $F\left(0\right)>0$ when $R_2>1.$ Further,

$$\underset{V⟶{\left(e^{-d_3{\tau }_3}\frac{\eta {\xi }_U}{u{\xi }_V}\right)}^{-}}{lim}F\left(V\right)=-\infty .$$

Furthermore,

$$F^{\prime }\left(V\right)=-\frac{{\xi }_V(\kappa +{\xi }_N)e^{d_1{\tau }_1+d_2{\tau }_2+d_3{\tau }_3}}{\eta \rho \kappa {\left(\eta {\xi }_U-e^{d_3{\tau }_3}u{\xi }_{V}V\right)}^2}\left[\rho \left({\xi }_Y{\left(e^{d_3{\tau }_3}u{\xi }_{V}V-\eta {\xi }_U\right)}^2+\gamma \mu {\eta }^2{\xi }_U\right)+e^{d_3{\tau }_3}\gamma \mu u\eta {\xi }_X{\xi }_V\right].$$

Hence, $F^{\prime }\left(V\right)<0$ for all $V\in \left(0,e^{-d_3{\tau }_3}\frac{\eta {\xi }_U}{u{\xi }_V}\right).$ It follows that there exists a unique $V_2\in \left(0,e^{-d_3{\tau }_3}\frac{\eta {\xi }_U}{u{\xi }_V}\right)$ such that $F\left(V_2\right)=0.$ From Equation (22) we get $Y_2>0,$$U_2>0,$$X_2>0$ and $N_2>0.$ As a result, $E{P}_2$ exists when $R_2>1$. The parameter $R_2$ represents the basic reproduction number of SARS-CoV-2-mono-infection. It determines the establishment of SARS-CoV-2 mono-infection.

(iv) HTLV-I and SARS-CoV-2 coinfection equilibrium point $E{P}_3=(X_3,N_3,Y_3,V_3,U_3,$$L_3,A_3),$ where

$$\begin{array}{cc}{X}_3=\frac{X_0}{R_4},\hfill & N_3=\frac{{\xi }_X{\xi }_V\left({\xi }_Y\alpha \pi e^{-d_4{\tau }_4-d_5{\tau }_5}+\mu \left({\xi }_L{\xi }_A+\alpha e^{-d_5{\tau }_5}\left({\xi }_A-\epsilon \right)\right)\right)}{\eta \rho \kappa \alpha \pi e^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3-d_4{\tau }_4-d_5{\tau }_5}}\left(R_4-1\right),\hfill \\ Y_3=\frac{{\xi }_X{\xi }_V}{\eta \rho e^{-d_1{\tau }_1-d_3{\tau }_3}}\left(R_4-1\right),\hfill & V_3=\frac{{\xi }_X}{\rho e^{-d_1{\tau }_1}}\left(R_4-1\right),\hfill \\ U_3=\frac{{\xi }_L{\xi }_A+\alpha e^{-d_5{\tau }_5}\left({\xi }_A-\epsilon \right)}{\alpha \pi e^{-d_4{\tau }_4-d_5{\tau }_5}},\hfill & L_3=\frac{{\xi }_A\left(u{\xi }_X{\xi }_V+\eta \rho {\xi }_U{e}^{-d_1{\tau }_1-d_3{\tau }_3}\right)}{\eta \rho \alpha \pi e^{-d_1{\tau }_1-d_3{\tau }_3-d_4{\tau }_4-d_5{\tau }_5}}\left(R_3-1\right),\hfill \\ A_3=\frac{\left(u{\xi }_X{\xi }_V+\eta \rho {\xi }_U{e}^{-d_3{\tau }_3}\right)}{\eta \rho \pi e^{-d_1{\tau }_1-d_3{\tau }_3}}\left(R_3-1\right),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill R_3& =\frac{\eta \rho \alpha \pi e^{-d_1{\tau }_1-d_3{\tau }_3-d_4{\tau }_4-d_5{\tau }_5}}{u{\xi }_X{\xi }_V+\eta \rho {\xi }_U{e}^{-d_1{\tau }_1-d_3{\tau }_3}}\times \hfill \\ \hfill & \left(\frac{\gamma }{\left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)\right)}+\frac{\delta \kappa u{e}^{-d_2{\tau }_2}}{\left(\kappa +{\xi }_N\right)\left[{\xi }_Y\alpha \pi e^{-d_4{\tau }_4-d_5{\tau }_5}+\mu \left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)\right)\right]}\right),\hfill \\ \hfill R_4& =\frac{\delta \eta \rho \kappa \alpha \pi e^{-d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3-d_4{\tau }_4-d_5{\tau }_5}}{{\xi }_X{\xi }_V\left(\kappa +{\xi }_N\right)\left[{\xi }_Y\alpha \pi e^{-d_4{\tau }_4-d_5{\tau }_5}+\mu \left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)\right)\right]}.\hfill \end{array}$$

Thus, $E{P}_3$ exists when $R_3>1$ and $R_4>1.$ At this point, $R_3$ and $R_4$ are threshold numbers that determine the occurrence of HTLV-I/SARS-CoV-2 coinfection.

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) if $R_1\le 1$, then the uninfected equilibrium point, $E{P}_0=\left(X_0,0,0,0,U_0,0,0\right)$ is the only equilibrium point, (b) if $R_1>1$, then, in addition to $E{P}_0$, there is an HTLV-I mono-infection equilibrium point, $E{P}_1=\left(X_1,0,0,0,U_1,L_1,A_1\right)$, (c) if $R_2>1$, then, in addition to $E{P}_0$, there is a SARS-CoV-2-mono-infection equilibrium point, $E{P}_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 $E{P}_0$, there is an HTLV-I and SARS-CoV-2 coinfection equilibrium point $E{P}_3=\left(X_3,N_3,Y_3,V_3,U_3,L_3,A_3\right)$.

2.3. Global Stability Analysis

In this section we discuss the global asymptotic stability of the four equilibrium points $E{P}_i$, $i=0,1,2,3$. We follow the work of [61] to build suitable Lyapunov function and apply Lyapunov–LaSalle asymptotic stability theorem (L-LAST) [62,63,64].

Denote $\left(X,N,Y,V,U,L,A\right)=\left(X\left(t\right),N\left(t\right),Y\left(t\right),V\left(t\right),U\left(t\right),L\left(t\right),A\left(t\right)\right)$. Let ${\Delta }_j^{\prime }$ be the largest invariant subset of

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

where, ${\Phi }_j\left(X,N,Y,V,U,L,A\right)$ is a Lyapunov function candidate.

Theorem 1.

If $R_1\le 1$ and $R_2\le 1$, then the uninfected equilibrium point $E{P}_0$ is globally asymptotically stable (GAS).

Proof. 

Define ${\Phi }_0$ as:

$$\begin{array}{cc}\hfill {\Phi }_0& =X_0\mathcal{H}\left(\frac{X}{X_0}\right)+e^{d_1{\tau }_1}N+\frac{\left(\kappa +{\xi }_N\right)}{\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}Y+\frac{\rho X_0}{{\xi }_V}V+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}{U}_0\mathcal{H}\left(\frac{U}{U_0}\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}L+\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}A\hfill \\ \hfill & +{\int }_{t-{\tau }_1}^t\rho V\left(\theta \right)X\left(\theta \right)d\theta +\frac{\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}}{\kappa }{\int }_{t-{\tau }_2}^t\kappa N\left(\theta \right)d\theta +\frac{\rho X_0{e}^{-d_3{\tau }_3}}{{\xi }_V}{\int }_{t-{\tau }_3}^t\eta Y\left(\theta \right)d\theta \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1+d_2{\tau }_2}}{u\kappa }{\int }_{t-{\tau }_4}^t\pi A\left(\theta \right)U\left(\theta \right)d\theta +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}{\int }_{t-{\tau }_5}^t\alpha L\left(\theta \right)d\theta ,\hfill \end{array}$$

where, $\mathcal{H}\left(x\right)=x-1-lnx$. Obviously, ${\Phi }_0\left(X,N,Y,V,U,L,A\right)>0$ for all $X,N,Y,V,U,L,$$A>0$, while ${\Phi }_0(X_0,0,$$0,$$0,$$U_0,0,0)=0$. The derivative of ${\Phi }_0$ w.r.t. t along the solutions of system (15)–(21) is calculated as:

$$\begin{array}{cc}\hfill \frac{d{\Phi }_0}{dt}& =\left(1-\frac{X_0}{X}\right)\dot{X}+e^{d_1{\tau }_1}\dot{N}+\frac{\left(\kappa +{\xi }_N\right)}{\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\dot{Y}+\frac{\rho X_0}{{\xi }_V}\dot{V}+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\left(1-\frac{U_0}{U}\right)\dot{U}\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}\dot{L}+\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}\dot{A}+\rho VX\hfill \\ \hfill & -\rho V\left(t-{\tau }_1\right)X\left(t-{\tau }_1\right)+\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}N-\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}N\left(t-{\tau }_2\right)+\frac{\rho X_0}{{\xi }_V}\eta e^{-d_3{\tau }_3}Y\hfill \\ \hfill & -\frac{\rho X_0}{{\xi }_V}\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\pi e^{d_1{\tau }_1+d_2{\tau }_2}AU-\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\pi e^{d_1{\tau }_1+d_2{\tau }_2}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}L-\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}L\left(t-{\tau }_5\right)\hfill \\ \hfill & =\left(1-\frac{X_0}{X}\right)(\delta -{\xi }_{X}X-\rho VX)+e^{d_1{\tau }_1}\left(\rho e^{-d_1{\tau }_1}V(t-{\tau }_1)X\left(t-{\tau }_1\right)-\left(\kappa +{\xi }_N\right)N\right)\hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\left(\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)-{\xi }_{Y}Y-\mu YU\right)+\frac{\rho X_0}{{\xi }_V}\left(\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)-{\xi }_{V}V\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\left(1-\frac{U_0}{U}\right)\left(\gamma +uYU-{\xi }_{U}U-\pi AU\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}\left(\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\epsilon A-\left(\alpha +{\xi }_L\right)L\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}\left(\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)-{\xi }_{A}A\right)+\rho VX\hfill \\ \hfill & -\rho V\left(t-{\tau }_1\right)X\left(t-{\tau }_1\right)+\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}N-\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}N\left(t-{\tau }_2\right)+\frac{\rho X_0}{{\xi }_V}\eta e^{-d_3{\tau }_3}Y\hfill \\ \hfill & -\frac{\rho X_0}{{\xi }_V}\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\pi e^{d_1{\tau }_1+d_2{\tau }_2}AU\hfill \\ \hfill & -\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\pi e^{d_1{\tau }_1+d_2{\tau }_2}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}L\hfill \\ \hfill & -\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}L\left(t-{\tau }_5\right).\hfill \end{array}$$

Since $\delta ={\xi }_X{X}_0$ and $\gamma ={\xi }_U{U}_0,$ then

$$\begin{array}{cc}\hfill \frac{d{\Phi }_0}{dt}& =-{\xi }_X\frac{{\left(X-X_0\right)}^2}{X}-\frac{\mu {\xi }_U\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1+d_2{\tau }_2}}{u\kappa }\frac{{\left(U-U_0\right)}^2}{U}+\frac{\left(\kappa +{\xi }_N\right)\left({\xi }_Y{\xi }_U+\mu \gamma \right)e^{d_1{\tau }_1+d_2{\tau }_2}}{\kappa {\xi }_U}\left(R_2-1\right)Y\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)\right)}{u\kappa \alpha }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}\left(R_1-1\right)A.\hfill \end{array}$$

Therefore, if $R_1\le 1$ and $R_2\le 1,$ then $\frac{d{\Phi }_0}{dt}\le 0$ for all $X,Y,U,A>0$ and $\frac{d{\Phi }_0}{dt}=0$ when $X=X_0$,$U=U_0$, $\left(R_2-1\right)Y=0$ and $\left(R_1-1\right)A=0$. The solutions of system (15)–(21) converge to ${\Delta }_0^{\prime }$ which comprises elements with $X=X_0$,$U=U_0$, $\left(R_2-1\right)Y=0$ and $\left(R_1-1\right)A=0$. Let consider four cases:

(i) $R_1<1$, $R_2<1$ and $Y=A=0$. Hence, $\dot{Y}=\dot{A}=0$ and Equations (17) and (21) yield

$$\begin{array}{cc}\hfill 0& =\dot{Y}=\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)⟹N\left(t\right)=0,\mathrm{for}\mathrm{all}t,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill 0& =\dot{A}=\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)⟹L\left(t\right)=0,\mathrm{for}\mathrm{all}t.\hfill \end{array}$$

(25)

Further, Equation (15) gives

$$0=\dot{X}=\delta -{\xi }_X{X}_0-\rho V{X}_0⟹V\left(t\right)=0,\mathrm{for}\mathrm{all}t,$$

(26)

(ii) $R_1=R_2=1$. From Equation (26) we have $V\left(t\right)=0$, for all t, and from Equation (18) we get

$$0=\dot{V}=\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)⟹Y\left(t\right)=0,\mathrm{for}\mathrm{all}t.$$

(27)

Equation (19) gives

$$0=\dot{U}=\gamma -{\xi }_U{U}_0-\pi A{U}_0⟹A\left(t\right)=0,\mathrm{for}\mathrm{all}t.$$

(28)

Using Equations (24)–(25) we obtain $N\left(t\right)=L\left(t\right)=0,$ for all $t.$

(iii) $R_1<1$, $R_2=1$ and $A=0$. From Equations (25)–(27) we obtain $L\left(t\right)=V\left(t\right)=Y\left(t\right)=0,$ for all t. Moreover, from Equation (24) we obtain $N\left(t\right)=0,$ for all t.

(iv) $R_1=1$, $R_2<1$ and $Y=0$. From Equations (24), (26) and (28) we obtain $N\left(t\right)=V\left(t\right)=A\left(t\right)=0$, for all t. Finally, Equation (25) implies that $L\left(t\right)=0,$ for all t. In all cases (i)–(iv), we have ${\Delta }_0^{\prime }=\left\{E{P}_0\right\}$. We deduce from L-LAST that, if $R_1\le 1$ and $R_2\le 1$, then $E{P}_0$ is GAS [62,63,64]. □

The result of Theorem 1 suggests that, when $R_1\le 1$ and $R_2\le 1$, both SARS-CoV-2 and HTLV-I infections are predicted to clear regardless of the initial states (any disease stages).

Theorem 2.

If $R_1>1$ and $R_4\le 1$, then the HTLV-I mono-infection equilibrium point $E{P}_1$ is GAS.

Proof. 

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

$$\begin{array}{cc}\hfill {\Phi }_1& =X_1\mathcal{H}\left(\frac{X}{X_1}\right)+e^{d_1{\tau }_1}N+\frac{\left(\kappa +{\xi }_N\right)}{\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}Y+\frac{\rho X_1}{{\xi }_V}V+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}{U}_1\mathcal{H}\left(\frac{U}{U_1}\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}{L}_1\mathcal{H}\left(\frac{L}{L_1}\right)+\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}{A}_1\mathcal{H}\left(\frac{A}{A_1}\right)\hfill \\ \hfill & +{\int }_{t-{\tau }_1}^t\rho V\left(\theta \right)X\left(\theta \right)d\theta +\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}{\int }_{t-{\tau }_2}^{t}N\left(\theta \right)d\theta +\frac{\rho X_1}{{\xi }_V}{e}^{-d_3{\tau }_3}{\int }_{t-{\tau }_3}^t\eta Y\left(\theta \right)d\theta \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\pi A_1{U}_1{\int }_{t-{\tau }_4}^t\mathcal{H}\left(\frac{A\left(\theta \right)U\left(\theta \right)}{A_1{U}_1}\right)d\theta \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}{L}_1{\int }_{t-{\tau }_5}^t\mathcal{H}\left(\frac{L\left(\theta \right)}{L_1}\right)d\theta .\hfill \end{array}$$

Clearly ${\Phi }_1\left(X,N,Y,V,U,L,A\right)>0$ for all $X,N,Y,V,U,L,A>0$ and ${\Phi }_1(X_1,0,$$0,$$0,$$U_1,$$L_1,$$A_1)=0$. Calculate $\frac{d{\Phi }_1}{dt}$ as:

$$\begin{array}{cc}\hfill & \frac{d{\Phi }_1}{dt}=\left(1-\frac{X_1}{X}\right)(\delta -{\xi }_{X}X-\rho VX)+e^{d_1{\tau }_1}\left(\rho e^{-d_1{\tau }_1}V(t-{\tau }_1)X\left(t-{\tau }_1\right)-\left(\kappa +{\xi }_N\right)N\right)\hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\left(\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)-{\xi }_{Y}Y-\mu YU\right)+\frac{\rho X_1}{{\xi }_V}\left(\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)-{\xi }_{V}V\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\left(1-\frac{U_1}{U}\right)\left(\gamma +uYU-{\xi }_{U}U-\pi AU\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}\left(1-\frac{L_1}{L}\right)\times \hfill \\ \hfill & \left(\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\epsilon A-\left(\alpha +{\xi }_L\right)L\right)+\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}}{\alpha u\kappa }\left(1-\frac{A_1}{A}\right)\times \hfill \\ \hfill & \left(\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)-{\xi }_{A}A\right)+\rho VX-\rho V\left(t-{\tau }_1\right)X\left(t-{\tau }_1\right)+\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}N-\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}N\left(t-{\tau }_2\right)\hfill \\ \hfill & +\frac{\rho X_1}{{\xi }_V}{e}^{-d_3{\tau }_3}\eta Y-\frac{\rho X_1}{{\xi }_V}{e}^{-d_3{\tau }_3}\eta Y\left(t-{\tau }_3\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\pi A_1{U}_1\left(\frac{AU}{A_1{U}_1}-\frac{A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)}{A_1{U}_1}\right.\hfill \\ \hfill & \left.+ln\left(\frac{A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)}{AU}\right)\right)+\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}}{u\kappa }{L}_1\left(\frac{L}{L_1}-\frac{L\left(t-{\tau }_5\right)}{L_1}+ln\left(\frac{L\left(t-{\tau }_5\right)}{L}\right)\right).\hfill \end{array}$$

Utilizing the equilibrium point conditions for $E{P}_1$:

$$\begin{array}{cc}\delta ={\xi }_X{X}_1,\hfill & \gamma ={\xi }_U{U}_1+\pi A_1{U}_1,\hfill \\ \pi e^{-d_4{\tau }_4}{A}_1{U}_1=\left(\alpha +{\xi }_L\right)L_1-\epsilon A_1,\hfill & \alpha e^{-d_5{\tau }_5}{L}_1={\xi }_A{A}_1,\hfill \end{array}$$

we obtain

$$\begin{array}{cc}\hfill \frac{d{\Phi }_1}{dt}& =-{\xi }_X\frac{{\left(X-X_1\right)}^2}{X}-\frac{\mu {\xi }_U\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1+d_2{\tau }_2}}{u\kappa }\frac{{\left(U-U_1\right)}^2}{U}\hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)\left[{\xi }_Y\alpha \pi e^{-d_4{\tau }_4-d_5{\tau }_5}+\mu \left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)\right)\right]e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}}{\kappa \alpha \pi }\left(R_4-1\right)Y\hfill \\ \hfill & -\frac{\mu \left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1+d_2{\tau }_2}}{u\kappa }\pi A_1{U}_1\left[\mathcal{H}\left(\frac{U_1}{U}\right)+\mathcal{H}\left(\frac{A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)L_1}{A_1{U}_{1}L}\right)+\mathcal{H}\left(\frac{L\left(t-{\tau }_5\right)A_1}{L_{1}A}\right)\right]\hfill \\ \hfill & -\frac{\mu \left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}}{u\kappa }\epsilon A_1\left[\mathcal{H}\left(\frac{A{L}_1}{A_{1}L}\right)+\mathcal{H}\left(\frac{L\left(t-{\tau }_5\right)A_1}{L_{1}A}\right)\right].\hfill \end{array}$$

Therefore, if $R_4\le 1$, then $\frac{d{\Phi }_1}{dt}\le 0$ for all $X,Y,U,L,A>0$, where $\frac{d{\Phi }_1}{dt}=0$ when $X=X_1,$$Y=0,$ $U=U_1,$$L=L_1$ and $A=A_1$. The solutions of system (15)–(21) are limited ${\Delta }_1^{\prime }$ which comprises elements with $X=X_1$,$U=U_1,$$L=L_1$, $A=A_1$ and $Y=0,$ then $\dot{Y}=0.$ Equation (17) yields

$$0=\dot{Y}=\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)⟹N\left(t\right)=0,\mathrm{for}\mathrm{all}t.$$

Equation (16) gives

$$0=\dot{N}=\rho e^{-d_1{\tau }_1}V(t-{\tau }_1)X_1⟹V\left(t\right)=0,\mathrm{for}\mathrm{all}t.$$

Consequently, ${\Delta }_1^{\prime }=\left\{E{P}_1\right\}$ and then L-LAST implies that $E{P}_1$ is GAS. □

The result of Theorem 2 demonstrates that, when $R_1>1$ and $R_4\le 1$, then the HTLV-I mono-infection is always established regardless of the initial states.

Theorem 3.

If $R_2>1$ and $R_3\le 1$, then the SARS-CoV-2-mono-infection equilibrium point, $E{P}_2$ is GAS.

Proof. 

Define ${\Phi }_2$ as:

$$\begin{array}{cc}\hfill {\Phi }_2& =X_2\mathcal{H}\left(\frac{X}{X_2}\right)+e^{d_1{\tau }_1}{N}_2\mathcal{H}\left(\frac{N}{N_2}\right)+\frac{\left(\kappa +{\xi }_N\right)}{\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}{Y}_2\mathcal{H}\left(\frac{Y}{Y_2}\right)+\frac{\rho X_2}{{\xi }_V}{V}_2\mathcal{H}\left(\frac{V}{V_2}\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}{U}_2\mathcal{H}\left(\frac{U}{U_2}\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}L\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}A+\rho V_2{X}_2{\int }_{t-{\tau }_1}^t\mathcal{H}\left(\frac{V\left(\theta \right)X\left(\theta \right)}{V_2{X}_2}\right)d\theta \hfill \\ \hfill & +\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}{N}_2{\int }_{t-{\tau }_2}^t\mathcal{H}\left(\frac{N\left(\theta \right)}{N_2}\right)d\theta +\frac{\rho X_2}{{\xi }_V}{e}^{-d_3{\tau }_3}\eta Y_2{\int }_{t-{\tau }_3}^t\mathcal{H}\left(\frac{Y\left(\theta \right)}{Y_2}\right)d\theta \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\pi e^{d_1{\tau }_1+d_2{\tau }_2}{\int }_{t-{\tau }_4}^{t}A\left(\theta \right)U\left(\theta \right)d\theta +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}{\int }_{t-{\tau }_5}^{t}L\left(\theta \right)d\theta .\hfill \end{array}$$

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

$$\begin{array}{cc}\hfill \frac{d{\Phi }_2}{dt}& =\left(1-\frac{X_2}{X}\right)(\delta -{\xi }_{X}X-\rho VX)+e^{d_1{\tau }_1}\left(1-\frac{N_2}{N}\right)\left(\rho e^{-d_1{\tau }_1}V(t-{\tau }_1)X\left(t-{\tau }_1\right)-\left(\kappa +{\xi }_N\right)N\right)\hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\left(1-\frac{Y_2}{Y}\right)\left(\kappa e^{-d_2{\tau }_2}N\left(t-{\tau }_2\right)-{\xi }_{Y}Y-\mu YU\right)+\frac{\rho X_2}{{\xi }_V}\left(1-\frac{V_2}{V}\right)\times \hfill \\ \hfill & \left(\eta e^{-d_3{\tau }_3}Y\left(t-{\tau }_3\right)-{\xi }_{V}V\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\left(1-\frac{U_2}{U}\right)\left(\gamma +uYU-{\xi }_{U}U-\pi AU\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}\left(\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\epsilon A-\left(\alpha +{\xi }_L\right)L\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}\left(\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)-{\xi }_{A}A\right)+\rho V_2{X}_2\times \hfill \\ \hfill & \left(\frac{VX}{V_2{X}_2}-\frac{V\left(t-{\tau }_1\right)X\left(t-{\tau }_1\right)}{V_2{X}_2}+ln\left(\frac{V\left(t-{\tau }_1\right)X\left(t-{\tau }_1\right)}{VX}\right)\right)+\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}{N}_2\times \hfill \\ \hfill & \left(\frac{N}{N_2}-\frac{N\left(t-{\tau }_2\right)}{N_2}+ln\left(\frac{N\left(t-{\tau }_2\right)}{N}\right)\right)+e^{-d_3{\tau }_3}\frac{\rho X_2}{{\xi }_V}\eta Y_2\left(\frac{Y}{Y_2}-\frac{Y\left(t-{\tau }_3\right)}{Y_2}+ln\left(\frac{Y\left(t-{\tau }_3\right)}{Y}\right)\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\pi e^{d_1{\tau }_1+d_2{\tau }_2}\left(AU-A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}\left(L-L\left(t-{\tau }_5\right)\right).\hfill \end{array}$$

Utilizing the equilibrium point conditions for $E{P}_2$:

$$\begin{array}{cc}\delta ={\xi }_X{X}_2+\rho V_2{X}_2,\hfill & \rho e^{-d_1{\tau }_1}{V}_2{X}_2=(\kappa +{\xi }_N)N_2,\hfill \\ \kappa e^{-d_2{\tau }_2}{N}_2={\xi }_Y{Y}_2+\mu Y_2{U}_2,\hfill & \eta e^{-d_3{\tau }_3}{Y}_2={\xi }_V{V}_2,\hfill \\ \gamma ={\xi }_U{U}_2-u{Y}_2{U}_2,\hfill \end{array}$$

we obtain

$$\begin{array}{cc}\hfill \frac{d{\Phi }_2}{dt}& =-{\xi }_X\frac{{\left(X-X_2\right)}^2}{X}-e^{d_1{\tau }_1+d_2{\tau }_2}\frac{\mu \gamma \left(\kappa +{\xi }_N\right)}{u\kappa U_2}\frac{{\left(U-U_2\right)}^2}{U}\hfill \\ \hfill & -\rho V_2{X}_2\left[\mathcal{H}\left(\frac{X_2}{X}\right)+\mathcal{H}\left(\frac{V(t-{\tau }_1)X\left(t-{\tau }_1\right)N_2}{V_2{X}_{2}N}\right)+\mathcal{H}\left(\frac{N\left(t-{\tau }_2\right)Y_2}{N_{2}Y}\right)+\mathcal{H}\left(\frac{Y\left(t-{\tau }_3\right)V_2}{Y_{2}V}\right)\right]\hfill \\ \hfill & +e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\left(e^{-d_4{\tau }_4-d_5{\tau }_5}\pi U_2-\frac{{\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)}{\alpha }\right)A.\hfill \end{array}$$

Hence, if $R_3\le 1$, then $E{P}_3$ dose not exist since $A_3\le 0$ and $L_3\le 0.$ This implies that

$$\begin{array}{cc}\hfill \dot{A}\left(t\right)& =\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)-{\xi }_{A}A\le 0,\hfill \\ \hfill \dot{L}\left(t\right)& =\pi e^{-d_4{\tau }_4}A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\epsilon A-\left(\alpha +{\xi }_L\right)L\le 0.\hfill \end{array}$$

It follows that $\left(e^{-d_4{\tau }_4-d_5{\tau }_5}\pi U_2-\frac{{\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)}{\alpha }\right)A\le 0$ for all $A>0.$ Thus, $e^{-d_4{\tau }_4-d_5{\tau }_5}\pi U_2-\frac{{\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5{\tau }_5}\right)}{\alpha }\le 0$. Thus, $\frac{d{\Phi }_2}{dt}\le 0$ for all $X,N,Y,V,U,L,A>0$ and $\frac{d{\Phi }_2}{dt}=0$ when $X=X_2,$ $N=N_2,$$Y=Y_2,$$V=V_2,$$U=U_2$ and $A=0.$ The solutions of the system converge to ${\Delta }_2^{\prime }$ which comprises elements with $A=0.$ It follows that $\dot{A}=0$ and Equation (21) becomes

$$0=\dot{A}=\alpha e^{-d_5{\tau }_5}L\left(t-{\tau }_5\right)⟹L\left(t\right)=0\mathrm{for}\mathrm{all}t.$$

Therefore, ${\Delta }_2^{\prime }=\left\{E{P}_2\right\}.$ L-LAST implies that $E{P}_2$ is GAS. □

The result of Theorem 3 shows that, when $R_2>1$ and $R_3\le 1$, then the SARS-CoV-2-mono-infection is always established regardless of the initial states.

Let us define a parameter $\tilde{R}$ as:

$$\tilde{R}=\frac{\eta \rho \alpha \pi \gamma e^{-d_1{\tau }_1-d_3{\tau }_3-d_4{\tau }_4-d_5{\tau }_5}}{\left(u{\xi }_X{\xi }_V+\eta \rho {\xi }_U{e}^{-d_3{\tau }_3}\right)\left({\xi }_L{\xi }_A+\alpha e^{-d_5{\tau }_5}\left({\xi }_A-\epsilon \right)\right)}.$$

Theorem 4.

If $R_4>1$ and $1<R_3\le 1+\tilde{R}$, then the HTLV-I/SARS-CoV-2 coinfection equilibrium point $E{P}_3$ is GAS.

Proof. 

Define ${\Phi }_3$ as:

$$\begin{array}{cc}\hfill {\Phi }_3& =X_3\mathcal{H}\left(\frac{X}{X_3}\right)+e^{d_1{\tau }_1}{N}_3\mathcal{H}\left(\frac{N}{N_3}\right)+\frac{\left(\kappa +{\xi }_N\right)}{\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}{Y}_3\mathcal{H}\left(\frac{Y}{Y_3}\right)+\frac{\rho X_3}{{\xi }_V}{V}_3\mathcal{H}\left(\frac{V}{V_3}\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}{U}_3\mathcal{H}\left(\frac{U}{U_3}\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}{L}_3\mathcal{H}\left(\frac{L}{L_3}\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}{A}_3\mathcal{H}\left(\frac{A}{A_3}\right)+\rho V_3{X}_3{\int }_{t-{\tau }_1}^t\mathcal{H}\left(\frac{V\left(\theta \right)X\left(\theta \right)}{V_3{X}_3}\right)d\theta \hfill \\ \hfill & +\left(\kappa +{\xi }_N\right)e^{d_1{\tau }_1}{N}_3{\int }_{t-{\tau }_2}^t\mathcal{H}\left(\frac{N\left(\theta \right)}{N_3}\right)d\theta +\frac{\rho X_3}{{\xi }_V}{e}^{-d_3{\tau }_3}\eta Y_3{\int }_{t-{\tau }_3}^t\mathcal{H}\left(\frac{Y\left(\theta \right)}{Y_3}\right)d\theta \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2}\pi A_3{U}_3{\int }_{t-{\tau }_4}^t\mathcal{H}\left(\frac{A\left(\theta \right)U\left(\theta \right)}{A_3{U}_3}\right)d\theta \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{e}^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}{L}_3{\int }_{t-{\tau }_5}^t\mathcal{H}\left(\frac{L\left(\theta \right)}{L_3}\right)d\theta .\hfill \end{array}$$

Clearly, ${\Phi }_3\left(X,N,Y,V,U,L,A\right)>0$ for all $X,N,Y,V,U,L,A>0$ and ${\Phi }_3(X_3,N_3,Y_3,$$V_3,$$U_3,$$L_3,$$A_3)=0$. Calculate $\frac{d{\Phi }_3}{dt}$ as:

$$\begin{array}{cc}\hfill \frac{d{\Phi }_3}{dt}& =\left(1-\frac{X_3}{X}\right)(\delta -{\xi }_{X}X-\rho VX)+e^{d_1{\tau }_1}\left(1-\frac{N_3}{N}\right)\left(e^{-d_1{\tau }_1}\rho V(t-{\tau }_1)X\left(t-{\tau }_1\right)-\left(\kappa +{\xi }_N\right)N\right)\hfill \\ \hfill & +e^{d_1{\tau }_1+d_2{\tau }_2}\frac{\left(\kappa +{\xi }_N\right)}{\kappa }\left(1-\frac{Y_3}{Y}\right)\left(e^{-d_2{\tau }_2}\kappa N\left(t-{\tau }_2\right)-{\xi }_{Y}Y-\mu YU\right)+\frac{\rho X_3}{{\xi }_V}\left(1-\frac{V_3}{V}\right)\times \hfill \\ \hfill & \left(e^{-d_3{\tau }_3}\eta Y\left(t-{\tau }_3\right)-{\xi }_{V}V\right)+e^{d_1{\tau }_1+d_2{\tau }_2}\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\left(1-\frac{U_3}{U}\right)\left(\gamma +uYU-{\xi }_{U}U-\pi AU\right)\hfill \\ \hfill & +e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\left(1-\frac{L_3}{L}\right)\left(e^{-d_4{\tau }_4}\pi A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)+\epsilon A-\left(\alpha +{\xi }_L\right)L\right)\hfill \\ \hfill & +e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4+d_5{\tau }_5}\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa }\left(1-\frac{A_3}{A}\right)\left(e^{-d_5{\tau }_5}\alpha L\left(t-{\tau }_5\right)-{\xi }_{A}A\right)\hfill \\ \hfill & +\rho V_3{X}_3\left(\frac{VX}{V_3{X}_3}-\frac{V\left(t-{\tau }_1\right)X\left(t-{\tau }_1\right)}{V_3{X}_3}+ln\left(\frac{V\left(t-{\tau }_1\right)X\left(t-{\tau }_1\right)}{VX}\right)\right)\hfill \\ \hfill & +e^{d_1{\tau }_1}\left(\kappa +{\xi }_N\right)N_3\left(\frac{N}{N_3}-\frac{N\left(t-{\tau }_2\right)}{N_3}+ln\left(\frac{N\left(t-{\tau }_2\right)}{N}\right)\right)\hfill \\ \hfill & +e^{-d_3{\tau }_3}\frac{\rho X_3}{{\xi }_V}\eta Y_3\left(\frac{Y}{Y_3}-\frac{Y\left(t-{\tau }_3\right)}{Y_3}+ln\left(\frac{Y\left(t-{\tau }_3\right)}{Y}\right)\right)\hfill \\ \hfill & +e^{d_1{\tau }_1+d_2{\tau }_2}\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\pi A_3{U}_3\left(\frac{AU}{A_3{U}_3}-\frac{A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)}{A_3{U}_3}+ln\left(\frac{A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)}{AU}\right)\right)\hfill \\ \hfill & +e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa }{L}_3\left(\frac{L}{L_3}-\frac{L\left(t-{\tau }_5\right)}{L_3}+ln\left(\frac{L\left(t-{\tau }_5\right)}{L}\right)\right)\hfill \end{array}$$

Utilizing the equilibrium point conditions for $E{P}_3$:

$$\begin{array}{cc}\delta ={\xi }_X{X}_3+\rho V_3{X}_3,\hfill & e^{-d_1{\tau }_1}\rho V_3{X}_3=(\kappa +{\xi }_N)N_3,\hfill \\ e^{-d_2{\tau }_2}\kappa N_3={\xi }_Y{Y}_3+\mu Y_3{U}_3,\hfill & e^{-d_3{\tau }_3}\eta Y_3={\xi }_V{V}_3,\hfill \\ \gamma ={\xi }_U{U}_3-u{Y}_3{U}_3+\pi A_3{U}_3,\hfill & e^{-d_4{\tau }_4}\pi A_3{U}_3=\left(\alpha +{\xi }_L\right)L_3-\epsilon A_3,\hfill \\ e^{-d_5{\tau }_5}\alpha L_3={\xi }_A{A}_3,\hfill \end{array}$$

we get

$$\begin{array}{cc}\hfill \frac{d{\Phi }_3}{dt}& =-\frac{{\xi }_X}{X}{\left(X-X_3\right)}^2-\rho V_3{X}_3\left[\mathcal{H}\left(\frac{X_3}{X}\right)+\mathcal{H}\left(\frac{V(t-{\tau }_1)X\left(t-{\tau }_1\right)N_3}{V_3{X}_{3}N}\right)+\mathcal{H}\left(\frac{N\left(t-{\tau }_2\right)Y_3}{N_{3}Y}\right)\right.\hfill \\ \hfill & \left.+\mathcal{H}\left(\frac{Y\left(t-{\tau }_3\right)V_3}{Y_{3}V}\right)\right]-e^{d_1{\tau }_1+d_2{\tau }_2+d_4{\tau }_4}\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\epsilon A_3\left[\mathcal{H}\left(\frac{A{L}_3}{A_{3}L}\right)+\mathcal{H}\left(\frac{L\left(t-{\tau }_5\right)A_3}{L_{3}A}\right)\right]\hfill \\ \hfill & -e^{d_1{\tau }_1+d_2{\tau }_2}\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa }\pi A_3{U}_3\left[\mathcal{H}\left(\frac{U_3}{U}\right)+\mathcal{H}\left(\frac{A\left(t-{\tau }_4\right)U\left(t-{\tau }_4\right)L_3}{A_3{U}_{3}L}\right)+\mathcal{H}\left(\frac{L\left(t-{\tau }_5\right)A_3}{L_{3}A}\right)\right]\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(u{\xi }_X{\xi }_V+e^{-d_3{\tau }_3}\eta \rho {\xi }_U\right)}{e^{-2d_1{\tau }_1-d_2{\tau }_2-d_3{\tau }_3}\eta \rho \kappa u}\frac{{\left(U-U_3\right)}^2}{U}\left(R_3-\tilde{R}-1\right).\hfill \end{array}$$

Moreover, since $1<R_3\le 1+\tilde{R}$ we obtain $\frac{d{\Phi }_3}{dt}\le 0$ for all $X,N,Y,V,U,L,A>0$ and $\frac{d{\Phi }_3}{dt}=0$ when $X=X_3,$ $N=N_3,$$Y=Y_3,$$V=V_3,$$U=U_3,$$L=L_3$ and $A=A_3.$ The solutions of the system converge to ${\Delta }_3^{\prime }.$ Clearly, ${\Delta }_3^{\prime }=\left\{E{P}_3\right\}.$ and thus, L-LAST implies that $E{P}_3$ is GAS. □

The result of Theorem 4 suggests that, when $R_4>1$ and $1<R_3\le 1+\tilde{R}$, then the SARS-CoV-2 and HTLV–I coinfection is always established regardless of the initial states.

3. Model with Distributed-Time Delays

In Section 2, it is assumed that the time delay is constant. It means that, each cell or virus is subject to the same delay during a certain biological process. For example, ${\tau }_3$ is constant means that, the maturation time for all new produced virions is assumed to be same. In reality, the maturation time may vary among immature virions. To accommodate this, one may use a weighted average of all the possible values for the delay time. Therefore, in this section, we consider SARS-CoV-2/HTLV-I coinfection model with five types of distributed-time delays as:

$$\begin{array}{cc}\hfill \dot{X}\left(t\right)& =\delta -{\xi }_{X}X\left(t\right)-\rho V\left(t\right)X\left(t\right),\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \dot{N}\left(t\right)& =\rho \underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }V(t-\varpi )X\left(t-\varpi \right)d\varpi -\left(\kappa +{\xi }_N\right)N\left(t\right),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \dot{Y}\left(t\right)& =\kappa \underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }N\left(t-\varpi \right)d\varpi -{\xi }_{Y}Y\left(t\right)-\mu Y\left(t\right)U\left(t\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \dot{V}\left(t\right)& =\eta \underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }Y\left(t-\varpi \right)d\varpi -{\xi }_{V}V\left(t\right),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \dot{U}\left(t\right)& =\gamma +uY\left(t\right)U\left(t\right)-{\xi }_{U}U\left(t\right)-\pi A\left(t\right)U\left(t\right),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \dot{L}\left(t\right)& =\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }A\left(t-\varpi \right)U\left(t-\varpi \right)d\varpi +\epsilon A\left(t\right)-\left(\alpha +{\xi }_L\right)L\left(t\right),\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \dot{A}\left(t\right)& =\alpha \underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }L\left(t-\varpi \right)d\varpi -{\xi }_{A}A\left(t\right).\hfill \end{array}$$

(35)

Here, $\varpi$ is a random variable generated from probability distribution functions ${\mathsf{f}}_i\left(\varpi \right)$, $i=1,2,\dots ,5$ over the time interval $\left[0,m_i\right]$, $i=1,2,\dots ,5$ where $m_i$ is the upper limit of the delay period. Here ${\mathsf{f}}_i\left(\varpi \right)e^{-d_i\varpi }$ $i=1,2,\dots ,5$ represents the probability of surviving a compartment from time $t-\varpi$ to time t. Functions ${\mathsf{f}}_i\left(\varpi \right)$, $i=1,2,\dots ,5$ satisfy the following conditions:

(i) ${\mathsf{f}}_i\left(\varpi \right)>0$, (ii) $\underset{0}{\overset{m_i}{\int }}{\mathsf{f}}_i\left(\varpi \right)d\varpi =1$, (iii) $\underset{0}{\overset{m_i}{\int }}{e}^{-u\varpi }{\mathsf{f}}_i\left(\varpi \right)d\varpi <\infty ,$$u>0$.

Let us define ${ϝ}_i$ as:

$${ϝ}_i=\underset{0}{\overset{m_i}{\int }}{e}^{-d_i\varpi }{\mathsf{f}}_i\left(\varpi \right)d\varpi ,i=1,2,\dots ,5.$$

Clearly $0<{ϝ}_i\le 1,$$i=1,2,\dots ,5$. We assume that the initial conditions of system (29)–(35) are the same as given by Equation (14) by replacing ${\tau }^{*}$ by $m^{*}=max\left\{m_1,m_2,m_3,m_4,m_5\right\}$.

3.1. Properties of Solutions

Lemma 3.

The solutions of model (29)–(35) satisfying the initial conditions in (14) are nonnegative and ultimately bounded.

Proof. 

From Equations (29) and (33), we obtain $\dot{X}{\mid }_{X=0}=\delta >0$ and $\dot{U}{\mid }_{U=0}=\gamma >0$. It follows that $X\left(t\right)>0$ and $U\left(t\right)>0$, for all $t\ge 0$. From Equations (30)–(32), (34) and (35) we have

$$\begin{array}{cc}\hfill N\left(t\right)& =e^{-\left(\kappa +{\xi }_N\right)t}{\varphi }_2\left(0\right)+\rho \underset{0}{\overset{t}{\int }}{e}^{-\left(\kappa +{\xi }_N\right)(t-\theta )}\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }V(\theta -\varpi )X\left(\theta -\varpi \right)d\varpi d\theta \ge 0,\hfill \\ \hfill Y\left(t\right)& =e^{-\underset{0}{\overset{t}{\int }}\left({\xi }_Y+\mu U\left(\varkappa \right)\right)d\varkappa }{\varphi }_3\left(0\right)+\kappa \underset{0}{\overset{t}{\int }}{e}^{-\underset{\theta }{\overset{t}{\int }}\left({\xi }_Y+\mu U\left(\varkappa \right)\right)d\varkappa }\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }N\left(\theta -\varpi \right)d\varpi d\theta \ge 0,\hfill \\ \hfill V\left(t\right)& =e^{-{\xi }_{V}t}{\varphi }_4\left(0\right)+\eta \underset{0}{\overset{t}{\int }}{e}^{-{\xi }_V(t-\theta )}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }Y\left(\theta -\varpi \right)d\varpi d\theta \ge 0,\hfill \\ \hfill L\left(t\right)& =e^{-\left(\alpha +{\xi }_L\right)t}{\varphi }_6\left(0\right)+\underset{0}{\overset{t}{\int }}{e}^{-\left(\alpha +{\xi }_L\right)(t-\theta )}\left[\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }A\left(\theta -\varpi \right)U\left(\theta -\varpi \right)d\varpi +\epsilon A\left(\theta \right)\right]d\theta \ge 0,\hfill \\ \hfill A\left(t\right)& =e^{-{\xi }_{A}t}{\varphi }_7\left(0\right)+\alpha \underset{0}{\overset{t}{\int }}{e}^{-{\xi }_A(t-\theta )}\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }L\left(\theta -\varpi \right)d\varpi d\theta \ge 0,\hfill \end{array}$$

for all $t\in \left[0,m^{*}\right].$ Hence, by a recursive argument, we obtain $N\left(t\right),Y\left(t\right),V\left(t\right),L\left(t\right)$,$A\left(t\right)\ge 0$ for all $t\ge 0.$ This guarantees that $\left(X\left(t\right),N\left(t\right),Y\left(t\right),V\left(t\right),U\left(t\right),L\left(t\right),A\left(t\right)\right)\in {\mathbb{R}}_{\ge 0}^7.$ From Equation (29), we have $\underset{t\to \infty }{lim\; sup}X\left(t\right)\le \frac{\delta }{{\xi }_X}=M_1$. We define

$$\begin{array}{cc}\hfill \Theta \left(t\right)& =\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }X\left(t-\varpi \right)d\varpi +N\left(t\right)+Y\left(t\right)+\kappa \underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_2\left(t-\theta \right)}N\left(\theta \right)d\theta d\varpi +\frac{{\xi }_Y}{2\eta }V\left(t\right)\hfill \\ \hfill & +\frac{{\xi }_Y}{2}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_3\left(t-\theta \right)}Y\left(\theta \right)d\theta d\varpi +\frac{\mu }{u}\left[U\left(t\right)+L\left(t\right)+A\left(t\right)\right.\hfill \\ \hfill & \left.+\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_4\left(t-\theta \right)}A\left(\theta \right)U\left(\theta \right)d\theta d\varpi +\alpha \underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_5\left(t-\theta \right)}L\left(\theta \right)d\theta d\varpi \right].\hfill \end{array}$$

Then

$$\begin{array}{cc}\hfill \dot{\Theta }\left(t\right)& =\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }\dot{X}\left(t-\varpi \right)d\varpi +\dot{N}\left(t\right)+\dot{Y}\left(t\right)-d_2\kappa \underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_2\left(t-\theta \right)}N\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & +\kappa N\left(t\right)\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)d\varpi -\kappa \underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }N\left(t-\varpi \right)d\varpi +\frac{{\xi }_Y}{2\eta }\dot{V}\left(t\right)\hfill \\ \hfill & -\frac{d_3{\xi }_Y}{2}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_3\left(t-\theta \right)}Y\left(\theta \right)d\theta d\varpi +\frac{{\xi }_Y}{2}Y\left(t\right)\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)d\varpi -\frac{{\xi }_Y}{2}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }Y\left(t-\varpi \right)d\varpi \hfill \\ \hfill & +\frac{\mu }{u}\left[\dot{U}\left(t\right)+\dot{L}\left(t\right)+\dot{A}\left(t\right)-d_4\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_4\left(t-\theta \right)}A\left(\theta \right)U\left(\theta \right)d\theta d\varpi +\pi A\left(t\right)U\left(t\right)\underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)d\varpi \right.\hfill \\ \hfill & -\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }A\left(t-\varpi \right)U\left(t-\varpi \right)d\varpi -d_5\alpha \underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_5\left(t-\theta \right)}L\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & \left.+\alpha L\left(t\right)\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)d\varpi -\alpha \underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }L\left(t-\varpi \right)d\varpi \right]\hfill \\ \hfill & =\delta \underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }d\varpi +\frac{\mu }{u}\gamma -{\xi }_X\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }X\left(t-\varpi \right)d\varpi -{\xi }_{N}N\left(t\right)-\frac{{\xi }_Y}{2}Y\left(t\right)\hfill \\ \hfill & -d_2\kappa \underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_2\left(t-\theta \right)}N\left(\theta \right)d\theta d\varpi -\frac{{\xi }_Y{\xi }_V}{2\eta }V\left(t\right)-d_3\frac{{\xi }_Y}{2}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_3\left(t-\theta \right)}Y\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & -\frac{\mu }{u}\left[{\xi }_{U}U\left(t\right)+{\xi }_{L}L\left(t\right)+\left({\xi }_A-\epsilon \right)A\left(t\right)+d_4\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_4\left(t-\theta \right)}A\left(\theta \right)U\left(\theta \right)d\theta d\varpi \right.\hfill \\ \hfill & \left.+d_5\alpha \underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)\underset{t-\varpi }{\overset{t}{\int }}{e}^{-d_5\left(t-\theta \right)}L\left(\theta \right)d\theta d\varpi \right].\hfill \end{array}$$

We have ${\xi }_A-\epsilon ={\xi }_A^{*}-{\epsilon }^{*}>0$ and $0\le {ϝ}_1\le 1.$ Hence,

$$\dot{\Theta }\left(t\right)\le \delta +\frac{\mu }{u}\gamma -\sigma \Theta \left(t\right),$$

where $\sigma =min\left\{{\xi }_X,{\xi }_N,\frac{{\xi }_Y}{2},d_2,{\xi }_V,d_3,{\xi }_U,{\xi }_L,{\xi }_A^{*}-{\epsilon }^{*},d_4,d_5\right\}$. Hence, $\underset{t⟶\infty }{lim\; sup}\Theta \left(t\right)\le M_2$ for $t\ge 0,$ where $M_2=\frac{\delta +\frac{\mu }{u}\gamma }{\sigma }.$ Since $X,N,Y,V,U,L$ and A are all non-negative, then $\underset{t⟶\infty }{lim\; sup}N\left(t\right)\le M_2,$ $\underset{t⟶\infty }{lim\; sup}Y\left(t\right)\le M_2,$ $\underset{t⟶\infty }{lim\; sup}V\left(t\right)\le M_3,$ $\underset{t⟶\infty }{lim\; sup}U\left(t\right)\le M_4,$$\underset{t⟶\infty }{lim\; sup}L\left(t\right)\le M_4$ and $\underset{t⟶\infty }{lim\; sup}A\left(t\right)\le M_4$, where $M_3=\frac{2\eta }{{\xi }_Y}{M}_2$ and $M_4=\frac{u}{\mu }{M}_2$. Consequently, $X\left(t\right),N\left(t\right),Y\left(t\right),V\left(t\right),U\left(t\right),L\left(t\right)$ and $A\left(t\right)$ are all ultimately bounded. □

According to Lemma 3, it can be shown that the set $\Omega =\{\left(X,N,Y,V,U,L,A\right)\in {\mathbb{C}}_{+}^7:∥X∥\le M_1,$$∥N∥\le M_2,$$∥Y∥\le M_2,$$∥V∥\le M_3,$$∥U∥\le M_4,$$∥L∥\le M_4,$$∥A∥\le M_4\}$ is positively invariant for system (29)–(35).

3.2. Equilibrium Points

To calculate the equilibrium points of the system (29)–(35) we solve the following system:

$$\begin{array}{cc}\hfill 0& =\delta -{\xi }_{X}X-\rho VX,\hfill \\ \hfill 0& =\rho {ϝ}_{1}VX-\left(\kappa +{\xi }_N\right)N,\hfill \\ \hfill 0& =\kappa {ϝ}_{2}N-{\xi }_{Y}Y-\mu YU,\hfill \\ \hfill 0& =\eta {ϝ}_{3}Y-{\xi }_{V}V,\hfill \\ \hfill 0& =\gamma +uYU-{\xi }_{U}U-\pi AU,\hfill \\ \hfill 0& =\pi {ϝ}_{4}AU+\epsilon A-\left(\alpha +{\xi }_L\right)L,\hfill \\ \hfill 0& =\alpha {ϝ}_{5}L-{\xi }_{A}A.\hfill \end{array}$$

We find that system admits four equilibrium points. These equilibrium points can be given as the same as given in previous equilibrium points sub-section by replacing $e^{-d_i{\tau }_i}$ by ${ϝ}_i,i=1,2,\dots ,5.$ The parameters $R_1$, $R_2,$$R_3$ and $R_4$ will be given as:

$$\begin{array}{cc}\hfill R_1& =\frac{\pi \alpha \gamma {ϝ}_4{ϝ}_5}{{\xi }_U\left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon {ϝ}_5\right)\right)},R_2=\frac{\delta \eta \rho \kappa {\xi }_U{ϝ}_1{ϝ}_2{ϝ}_3}{{\xi }_X{\xi }_V\left(\kappa +{\xi }_N\right)\left({\xi }_Y{\xi }_U+\gamma \mu \right)},\hfill \\ \hfill R_3& =\frac{\eta \rho \alpha \pi {ϝ}_1{ϝ}_3{ϝ}_4{ϝ}_5}{u{\xi }_X{\xi }_V+\eta \rho {\xi }_U{ϝ}_1{ϝ}_3}\left(\frac{\gamma }{\left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon {ϝ}_5\right)\right)}+\frac{\delta \kappa u{ϝ}_2}{\left(\kappa +{\xi }_N\right)\left(\alpha \pi {\xi }_Y{ϝ}_4{ϝ}_5+\mu \left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon {ϝ}_5\right)\right)\right)}\right),\hfill \\ \hfill R_4& =\frac{\delta \eta \rho \kappa \alpha \pi {ϝ}_1{ϝ}_2{ϝ}_3{ϝ}_4{ϝ}_5}{{\xi }_X{\xi }_V\left(\kappa +{\xi }_N\right)\left[\alpha \pi {\xi }_Y{ϝ}_4{ϝ}_5+\mu \left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon {ϝ}_5\right)\right)\right]}.\hfill \end{array}$$

The results of Lemma 2 is valid for system (29)–(35).

3.3. Global Stability Analysis

In this subsection, we discuss the global stability of four equilibrium points $E{P}_i$, $i=0,1,2,3$ of system (29)–(35).

Theorem 5.

If $R_1\le 1$ and $R_2\le 1$, then the uninfected equilibrium point $E{P}_0$ is GAS.

Proof. 

Define ${\Phi }_0$ as:

$$\begin{array}{cc}\hfill {\Phi }_0& =X_0\mathcal{H}\left(\frac{X}{X_0}\right)+\frac{1}{{ϝ}_1}N+\frac{\left(\kappa +{\xi }_N\right)}{\kappa {ϝ}_1{ϝ}_2}Y+\frac{\rho X_0}{{\xi }_V}V+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}{U}_0\mathcal{H}\left(\frac{U}{U_0}\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}L\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}A+\frac{1}{{ϝ}_1}\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }{\int }_{t-\varpi }^t\rho V\left(\theta \right)X\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{\kappa {ϝ}_1{ϝ}_2}\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }{\int }_{t-\varpi }^t\kappa N\left(\theta \right)d\theta d\varpi +\frac{\rho X_0}{{\xi }_V}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }{\int }_{t-\varpi }^t\eta Y\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }{\int }_{t-\varpi }^t\pi A\left(\theta \right)U\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }{\int }_{t-\varpi }^t\alpha L\left(\theta \right)d\theta d\varpi .\hfill \end{array}$$

Obviously, ${\Phi }_0\left(X,N,Y,V,U,L,A\right)>0$ for all $X,N,Y,V,U,L,A>0$, while ${\Phi }_0(X_0,0,0,$$0,$$U_0,$$0,0)=0$. We calculate $\frac{d{\Phi }_0}{dt}$ as:

$$\begin{array}{cc}\hfill \frac{d{\Phi }_0}{dt}& =\left(1-\frac{X_0}{X}\right)(\delta -{\xi }_{X}X-\rho VX)+\frac{1}{{ϝ}_1}\left(\rho \underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }V(t-\varpi )X\left(t-\varpi \right)d\varpi -\left(\kappa +{\xi }_N\right)N\right)\hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{\kappa {ϝ}_1{ϝ}_2}\left(\kappa \underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }N\left(t-\varpi \right)d\varpi -{\xi }_{Y}Y-\mu YU\right)\hfill \\ \hfill & +\frac{\rho X_0}{{\xi }_V}\left(\eta \underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }Y\left(t-\varpi \right)d\varpi -{\xi }_{V}V\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}\left(1-\frac{U_0}{U}\right)\left(\gamma +uYU-{\xi }_{U}U-\pi AU\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\left(\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }A\left(t-\varpi \right)U\left(t-\varpi \right)d\varpi +\epsilon A-\left(\alpha +{\xi }_L\right)L\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}\left(\alpha \underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }L\left(t-\varpi \right)d\varpi -{\xi }_{A}A\right)\hfill \\ \hfill & +\frac{1}{{ϝ}_1}\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }\left(\rho VX-\rho V\left(t-\varpi \right)X\left(t-\varpi \right)\right)d\varpi \hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{{ϝ}_1{ϝ}_2}\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }\left(N-N\left(t-\varpi \right)\right)d\varpi +\frac{\rho X_0}{{\xi }_V}\eta \underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }\left(Y-Y\left(t-\varpi \right)\right)d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }\left(AU-A\left(t-\varpi \right)U\left(t-\varpi \right)\right)d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }\left(L-L\left(t-\varpi \right)\right)d\varpi \hfill \\ \hfill & =\left(1-\frac{X_0}{X}\right)(\delta -{\xi }_{X}X)-\frac{\left(\kappa +{\xi }_N\right)}{\kappa {ϝ}_1{ϝ}_2}{\xi }_{Y}Y+\frac{\rho X_0}{{\xi }_V}{ϝ}_3\eta Y+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}\left(1-\frac{U_0}{U}\right)\left(\gamma -{\xi }_{U}U\right)\hfill \\ \hfill & -\frac{\mu \left(\kappa +{\xi }_N\right)}{\kappa {ϝ}_1{ϝ}_2}{U}_{0}Y+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}\pi A{U}_0+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\epsilon A-\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}{\xi }_{A}A.\hfill \end{array}$$

Using $\delta ={\xi }_X{X}_0$ and $\gamma ={\xi }_U{U}_0,$ we get

$$\begin{array}{cc}\hfill \frac{d{\Phi }_0}{dt}& =-{\xi }_X\frac{{\left(X-X_0\right)}^2}{X}-\frac{\mu {\xi }_U\left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}\frac{{\left(U-U_0\right)}^2}{U}+\frac{\left(\kappa +{\xi }_N\right)\left({\xi }_Y{\xi }_U+\mu \gamma \right)}{\kappa {\xi }_U{ϝ}_1{ϝ}_2}\left(R_2-1\right)Y\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon {ϝ}_5\right)\right)}{u\kappa \alpha {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}\left(R_1-1\right)A.\hfill \end{array}$$

Therefore, if $R_1\le 1$ and $R_2\le 1,$ then $\frac{d{\Phi }_0}{dt}\le 0$ for all $X,Y,U,A>0$ and $\frac{d{\Phi }_0}{dt}=0$ when $X=X_0$,$U=U_0$ and $Y=A=0.$ Similar to the proof of Theorem 1 we can show that, ${\Delta }_0^{\prime }=\left\{E{P}_0\right\}$. We deduce from L-LAST that, $E{P}_0$ is GAS. □

Theorem 6.

If $R_1>1$ and $R_4\le 1$, then the HTLV-I mono-infection equilibrium point $E{P}_1$ is GAS.

Proof. 

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

$$\begin{array}{cc}\hfill {\Phi }_1& =X_1\mathcal{H}\left(\frac{X}{X_1}\right)+\frac{1}{{ϝ}_1}N+\frac{\left(\kappa +{\xi }_N\right)}{\kappa {ϝ}_1{ϝ}_2}Y+\frac{\rho X_1}{{\xi }_V}V+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}{U}_1\mathcal{H}\left(\frac{U}{U_1}\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}{L}_1\mathcal{H}\left(\frac{L}{L_1}\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}{A}_1\mathcal{H}\left(\frac{A}{A_1}\right)+\frac{1}{{ϝ}_1}\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }{\int }_{t-\varpi }^t\rho V\left(\theta \right)X\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{{ϝ}_1{ϝ}_2}\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }{\int }_{t-\varpi }^{t}N\left(\theta \right)d\theta d\varpi +\frac{\rho X_1}{{\xi }_V}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }{\int }_{t-\varpi }^t\eta Y\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\pi A_1{U}_1\underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{A\left(\theta \right)U\left(\theta \right)}{A_1{U}_1}\right)d\theta d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}{L}_1\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{L\left(\theta \right)}{L_1}\right)d\theta d\varpi .\hfill \end{array}$$

Calculating $\frac{d{\Phi }_1}{dt}$ and using the equilibrium point conditions for $E{P}_1$:

$$\begin{array}{cc}\hfill \delta & ={\xi }_X{X}_1,\gamma ={\xi }_U{U}_1+\pi A_1{U}_1,\hfill \\ \hfill \pi {ϝ}_4{A}_1{U}_1& =\left(\alpha +{\xi }_L\right)L_1-\epsilon A_1,\alpha {ϝ}_5{L}_1={\xi }_A{A}_1,\hfill \end{array}$$

we obtain

$$\begin{array}{cc}\hfill \frac{d{\Phi }_1}{dt}& =-\frac{{\xi }_X}{X}{\left(X-X_1\right)}^2-\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}\frac{{\xi }_U}{U}{\left(U-U_1\right)}^2\hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)\left[{\xi }_Y\alpha \pi {ϝ}_4{ϝ}_5+\mu \left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon {ϝ}_5\right)\right)\right]}{\kappa \alpha \pi {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}\left(R_4-1\right)Y-\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}\pi A_1{U}_1\left[\mathcal{H}\left(\frac{U_1}{U}\right)\right.\hfill \\ \hfill & \left.+\frac{1}{{ϝ}_4}\underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }\mathcal{H}\left(\frac{A\left(t-\varpi \right)U\left(t-\varpi \right)L_1}{A_1{U}_{1}L}\right)d\varpi +\frac{1}{{ϝ}_5}\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }\mathcal{H}\left(\frac{L\left(t-\varpi \right)A_1}{L_{1}A}\right)d\varpi \right]\hfill \\ \hfill & -\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\epsilon A_1\left[\mathcal{H}\left(\frac{A{L}_1}{A_{1}L}\right)+\frac{1}{{ϝ}_5}\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }\mathcal{H}\left(\frac{L\left(t-\varpi \right)A_1}{L_{1}A}\right)d\varpi \right].\hfill \end{array}$$

Therefore, if $R_4\le 1$, then $\frac{d{\Phi }_1}{dt}\le 0$ for all $X,Y,U,L,A>0$, where $\frac{d{\Phi }_1}{dt}=0$ when $X=X_1,$$Y=0,$ $U=U_1,$$L=L_1$ and $A=A_1$. Similar to the proof of Theorem 2 we can show that ${\Delta }_1^{\prime }=\left\{E{P}_1\right\}$ and then L-LAST implies that $E{P}_1$ is GAS. □

Theorem 7.

If $R_2>1$ and $R_3\le 1$, then the SARS-CoV-2-mono-infection equilibrium point, $E{P}_2$ is GAS.

Proof. 

Define ${\Phi }_2$ as:

$$\begin{array}{cc}\hfill {\Phi }_2& =X_2\mathcal{H}\left(\frac{X}{X_2}\right)+\frac{1}{{ϝ}_1}{N}_2\mathcal{H}\left(\frac{N}{N_2}\right)+\frac{\left(\kappa +{\xi }_N\right)}{\kappa {ϝ}_1{ϝ}_2}{Y}_2\mathcal{H}\left(\frac{Y}{Y_2}\right)+\frac{\rho X_2}{{\xi }_V}{V}_2\mathcal{H}\left(\frac{V}{V_2}\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}{U}_2\mathcal{H}\left(\frac{U}{U_2}\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}L+\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}A\hfill \\ \hfill & +\frac{\rho V_2{X}_2}{{ϝ}_1}\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{V\left(\theta \right)X\left(\theta \right)}{V_2{X}_2}\right)d\theta d\varpi \hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{{ϝ}_1{ϝ}_2}{N}_2\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{N\left(\theta \right)}{N_2}\right)d\theta d\varpi \hfill \\ \hfill & +\frac{\rho X_2}{{\xi }_V}\eta Y_2\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{Y\left(\theta \right)}{Y_2}\right)d\theta d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\pi \underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }{\int }_{t-\varpi }^{t}A\left(\theta \right)U\left(\theta \right)d\theta d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }{\int }_{t-\varpi }^{t}L\left(\theta \right)d\theta d\varpi .\hfill \end{array}$$

Calculating $\frac{d{\Phi }_2}{dt}$ and using the equilibrium point conditions for $E{P}_2$:

$$\begin{array}{cc}\delta ={\xi }_X{X}_2+\rho V_2{X}_2,\hfill & \rho {ϝ}_1{V}_2{X}_2=(\kappa +{\xi }_N)N_2,\hfill \\ \kappa {ϝ}_2{N}_2={\xi }_Y{Y}_2+\mu Y_2{U}_2,\hfill & \eta {ϝ}_3{Y}_2={\xi }_V{V}_2,\hfill \\ \gamma ={\xi }_U{U}_2-u{Y}_2{U}_2,\hfill \end{array}$$

we get

$$\begin{array}{cc}\hfill \frac{d{\Phi }_2}{dt}& =-{\xi }_X\frac{{\left(X-X_2\right)}^2}{X}-\frac{\mu \gamma \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{U}_2}\frac{{\left(U-U_2\right)}^2}{U}\hfill \\ \hfill & -\rho V_2{X}_2\left[\mathcal{H}\left(\frac{X_2}{X}\right)+\frac{1}{{ϝ}_1}\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }\mathcal{H}\left(\frac{V(t-\varpi )X\left(t-\varpi \right)N_2}{V_2{X}_{2}N}\right)d\varpi \right.\hfill \\ \hfill & \left.+\frac{1}{{ϝ}_2}\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }\mathcal{H}\left(\frac{N\left(t-\varpi \right)Y_2}{N_{2}Y}\right)d\varpi +\frac{1}{{ϝ}_3}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }\mathcal{H}\left(\frac{Y\left(t-\varpi \right)V_2}{Y_{2}V}\right)d\varpi \right]\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}\left(\pi {ϝ}_4{ϝ}_5{U}_2-\frac{{\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon {ϝ}_5\right)}{\alpha }\right)A.\hfill \end{array}$$

Similar to the proof of Theorem 3 we can show that ${\Delta }_2^{\prime }=\left\{E{P}_2\right\}.$ L-LAST implies that $E{P}_2$ is GAS. □

We define the parameter $\tilde{R}$ as:

$$\tilde{R}=\frac{\eta \rho \alpha \pi \gamma {ϝ}_1{ϝ}_3{ϝ}_4{ϝ}_5}{\left(u{\xi }_X{\xi }_V+\eta \rho {\xi }_U{ϝ}_3\right)\left({\xi }_L{\xi }_A+\alpha {ϝ}_5\left({\xi }_A-\epsilon \right)\right)}.$$

Theorem 8.

If $R_4>1$ and $1<R_3\le 1+\tilde{R}$, then the HTLV-I/SARS-CoV-2 coinfection equilibrium point $E{P}_3$ is GAS.

Proof. 

Define ${\Phi }_3$ as:

$$\begin{array}{cc}\hfill {\Phi }_3& =X_3\mathcal{H}\left(\frac{X}{X_3}\right)+\frac{1}{{ϝ}_1}{N}_3\mathcal{H}\left(\frac{N}{N_3}\right)+\frac{\left(\kappa +{\xi }_N\right)}{\kappa {ϝ}_1{ϝ}_2}{Y}_3\mathcal{H}\left(\frac{Y}{Y_3}\right)+\frac{\rho X_3}{{\xi }_V}{V}_3\mathcal{H}\left(\frac{V}{V_3}\right)\hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}{U}_3\mathcal{H}\left(\frac{U}{U_3}\right)+\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}{L}_3\mathcal{H}\left(\frac{L}{L_3}\right)+\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{\alpha u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}{A}_3\mathcal{H}\left(\frac{A}{A_3}\right)\hfill \\ \hfill & +\frac{\rho V_3{X}_3}{{ϝ}_1}\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{V\left(\theta \right)X\left(\theta \right)}{V_3{X}_3}\right)d\theta d\varpi \hfill \\ \hfill & +\frac{\left(\kappa +{\xi }_N\right)}{{ϝ}_1{ϝ}_2}{N}_3\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{N\left(\theta \right)}{N_3}\right)d\theta d\varpi \hfill \\ \hfill & +\frac{\rho X_3}{{\xi }_V}\eta Y_3\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{Y\left(\theta \right)}{Y_3}\right)d\theta d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\pi A_3{U}_3\underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{A\left(\theta \right)U\left(\theta \right)}{A_3{U}_3}\right)d\theta d\varpi \hfill \\ \hfill & +\frac{\mu \left(\kappa +{\xi }_N\right)\left(\alpha +{\xi }_L\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4{ϝ}_5}{L}_3\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }{\int }_{t-\varpi }^t\mathcal{H}\left(\frac{L\left(\theta \right)}{L_3}\right)d\theta d\varpi .\hfill \end{array}$$

Calculating $\frac{d{\Phi }_3}{dt}$ and using following equilibrium point conditions for $E{P}_3$:

$$\begin{array}{cc}\delta ={\xi }_X{X}_3+\rho V_3{X}_3,\hfill & {ϝ}_1\rho V_3{X}_3=(\kappa +{\xi }_N)N_3,\hfill \\ {ϝ}_2\kappa N_3={\xi }_Y{Y}_3+\mu Y_3{U}_3,\hfill & {ϝ}_3\eta Y_3={\xi }_V{V}_3,\hfill \\ \gamma ={\xi }_U{U}_3-u{Y}_3{U}_3+\pi A_3{U}_3,\hfill & {ϝ}_4\pi A_3{U}_3=\left(\alpha +{\xi }_L\right)L_3-\epsilon A_3,\hfill \\ {ϝ}_5\alpha L_3={\xi }_A{A}_3,\hfill \end{array}$$

we get

$$\begin{array}{cc}\hfill \frac{d{\Phi }_3}{dt}& =-{\xi }_X\frac{{\left(X-X_3\right)}^2}{X}+\frac{\mu \left(\kappa +{\xi }_N\right)\left(u{\xi }_X{\xi }_V+{ϝ}_3\eta \rho {\xi }_U\right)}{\eta \rho \kappa u{ϝ}_1^2{ϝ}_2{ϝ}_3}\frac{{\left(U-U_3\right)}^2}{U}\left(R_3-1-\tilde{R}\right)\hfill \\ \hfill & -\rho V_3{X}_3\left[\mathcal{H}\left(\frac{X_3}{X}\right)+\frac{1}{{ϝ}_1}\underset{0}{\overset{m_1}{\int }}{\mathsf{f}}_1\left(\varpi \right)e^{-d_1\varpi }\mathcal{H}\left(\frac{V(t-\varpi )X\left(t-\varpi \right)N_3}{V_3{X}_{3}N}\right)d\varpi \right.\hfill \\ \hfill & \left.+\frac{1}{{ϝ}_2}\underset{0}{\overset{m_2}{\int }}{\mathsf{f}}_2\left(\varpi \right)e^{-d_2\varpi }\mathcal{H}\left(\frac{N\left(t-\varpi \right)Y_3}{N_{3}Y}\right)d\varpi +\frac{1}{{ϝ}_3}\underset{0}{\overset{m_3}{\int }}{\mathsf{f}}_3\left(\varpi \right)e^{-d_3\varpi }\mathcal{H}\left(\frac{Y\left(t-\varpi \right)V_3}{Y_{3}V}\right)d\varpi \right]\hfill \\ \hfill & -\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2}\pi A_3{U}_3\left[\mathcal{H}\left(\frac{U_3}{U}\right)+\frac{1}{{ϝ}_4}\underset{0}{\overset{m_4}{\int }}{\mathsf{f}}_4\left(\varpi \right)e^{-d_4\varpi }\mathcal{H}\left(\frac{A\left(t-\varpi \right)U\left(t-\varpi \right)L_3}{A_3{U}_{3}L}\right)d\varpi \right.\hfill \\ \hfill & \left.+\frac{1}{{ϝ}_5}\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }\mathcal{H}\left(\frac{L\left(t-\varpi \right)A_3}{L_{3}A}\right)d\varpi \right]\hfill \\ \hfill & -\frac{\mu \left(\kappa +{\xi }_N\right)}{u\kappa {ϝ}_1{ϝ}_2{ϝ}_4}\epsilon A_3\left[\mathcal{H}\left(\frac{A{L}_3}{A_{3}L}\right)+\frac{1}{{ϝ}_5}\underset{0}{\overset{m_5}{\int }}{\mathsf{f}}_5\left(\varpi \right)e^{-d_5\varpi }\mathcal{H}\left(\frac{L\left(t-\varpi \right)A_3}{L_{3}A}\right)d\varpi \right].\hfill \end{array}$$

Since $1<R_3\le 1+\tilde{R}$, then we obtain $\frac{d{\Phi }_3}{dt}\le 0$ for all $X,N,Y,V,U,L,A>0$ and $\frac{d{\Phi }_3}{dt}=0$ when $X=X_3,$ $N=N_3,Y=Y_3,$$V=V_3,$$U=U_3,$$L=L_3$ and $A=A_3.$ The solutions of the system converge to ${\Delta }_3^{\prime }.$ Clearly, ${\Delta }_3^{\prime }=\left\{E{P}_3\right\},$ and thus L-LAST implies that $E{P}_3$ is GAS. □

Based on the above findings, we summarize the existence and global stability conditions for all equilibrium points in

Table 1. Conditions of existence and global stability of the system’s equilibria.

Equilibrium Point	Existence Conditions	Global Stability Conditions
$E{P}_0=\left(X_0,0,0,0,U_0,0,0\right)$	None	$R_1\le 1$ and $R_2\le 1$
$E{P}_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$
$E{P}_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$
$E{P}_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}$

.

Remark 1.

Let us choose a Dirac delta function $\mathrm{\Pi }\left(.\right)$ as a special form of the function ${\mathsf{f}}_i(.)$ as:

$${\mathsf{f}}_i\left(\varpi \right)=\mathrm{\Pi }\left(\varpi -{\tau }_i\right),{\tau }_i\in \left[0,m_i\right],i=1,2,\dots ,5,$$

and let $m_i$ tends to $\infty ,$ then the properties of $\mathrm{\Pi }\left(.\right)$ implies that:

$${\displaystyle \underset{0}{\overset{\infty }{\int }}}{\mathsf{f}}_i\left(\varpi \right)d\varpi =1,{ϝ}_i=\underset{0}{\overset{\infty }{\int }}{e}^{-d_i\varpi }{\mathsf{f}}_i\left(\varpi \right)d\varpi =e^{-d_i{\tau }_i},i=1,2,\dots ,5.$$

Then, the distributed-time delay model (29)–(35) will lead to the discrete-time delay model (7)–(13).

4. Numerical Simulations

In this section, we present some numerical results for model with discrete-time delays (7)–(13) to illustrate the stability of equilibrium points. Further, we illustrate the effect of time delays on the dynamical behaviors of the HTLV-I/SARS-CoV-2 coinfection. We used the dde23 solver in MATLAB to solve the system of delay differential equations numerically.

4.1. Stability of Equilibrium Points

We solve the system with three initial conditions:

$$\begin{array}{cc}\hfill \mathrm{IC}1& :\left(X\left(\theta \right),N\left(\theta \right),Y\left(\theta \right),V\left(\theta \right),U\left(\theta \right),L\left(\theta \right),A\left(\theta \right)\right)=\left(1,0.0001,0.0002,0.0003,600,150,15\right),\hfill \\ \hfill \mathrm{IC}2& :\left(X\left(\theta \right),N\left(\theta \right),Y\left(\theta \right),V\left(\theta \right),U\left(\theta \right),L\left(\theta \right),A\left(\theta \right)\right)=\left(3,0.0005,0.0006,0.0007,500,200,20\right),\hfill \\ \hfill \mathrm{IC}3& :\left(X\left(\theta \right),N\left(\theta \right),Y\left(\theta \right),V\left(\theta \right),U\left(\theta \right),L\left(\theta \right),A\left(\theta \right)\right)=\left(5,0.001,0.002,0.003,400,250,25\right),\hfill \end{array}$$

where $\theta \in \left[-max\{{\tau }_1,{\tau }_2,\dots ,{\tau }_5\},0\right]$. We fix the time delay parameters ${\tau }_i=0.1$, $i=1,$$2,\dots ,5$.

Table 2. Model parameters.

Parameter	Value	Sources	Parameter	Value	Sources
$\delta$	$0.11$	[26]	u	$0.1$	[35,65]
${\xi }_X$	$0.011$	[26]	${\xi }_U$	$0.012$	[49,51,52]
$\rho$	Varies	Assumed	$\pi$	Varies	Assumed
$\kappa$	$4.08$	[35,66]	$\omega$	$0.9$	[50]
${\xi }_N$	$0.11$	[26]	${\epsilon }^{*}$	$0.011$	[50]
${\xi }_Y$	$0.11$	[35,42]	$\alpha$	$0.003$	[49,51,52,67]
$\mu$	Varies	Assumed	${\xi }_L$	$0.03$	[50,51,52,53]
$\eta$	$0.24$	[35]	${\xi }_A^{*}$	$0.03$	[49,51,52]
${\xi }_V$	Varies	Assumed	$d_i,i=1,\dots ,5$	$1.0$	Assumed
$\gamma$	10	[51,52,68]	${\tau }_i,i=1,\dots ,5$	Varies	Assumed

contains the values of some parameters. We mention that, the values of some parameters of the model are taken from previous studies for SARS-CoV-2 and HTLV-I mono-infections, while other parameters such as the delay parameters ${\tau }_i$, $i=1,2,\dots ,5$ are chosen just to conduct the numerical simulations. To the best of our knowledge, till now there is no available data from SARS-CoV-2 and HTLV-I coinfection patients. Therefore, estimating the parameters of the coinfection model is still open for future investigations.

We vary the parameters, $\rho$, $\mu$, ${\xi }_V$ and $\pi$ to obtain four cases as:

Case 1.

($R_1\le 1$ and $R_2\le 1$): Choosing $\rho =0.8,$ $\mu =1.1,$${\xi }_V=5.5$ and $\pi =0.0001$ which gives $R_1=0.2208<1$ and $R_2=0.0003<1.$ Based on Theorem 1, the equilibrium point, $E{P}_0$ is GAS. This is illustrated in Figure 1, where, the concentrations of the uninfected ECs and uninfected CD4${}^{+}$T cells tends to their healthy values $X_0=10$ and $U_0=833.33$, while, the concentrations of the other compartments tend to zero. This case means the clearance of both SARS-CoV-2 and HTLV-I infections.

Case 2.

($R_1>1$ and $R_4\le 1$): We take $\rho =0.5,$$\mu =1.1,$ ${\xi }_V=5.5$ and $\pi =0.0015.$ So, we obtain $R_1=3.3126>1$ and $R_4=0.0006<1.$ According to Lemma 2 and Theorem 2, the HTLV-I mono-infection equilibrium point, $E{P}_1$ exists and is GAS. Figure 2 shows that the models’s solutions converge to the equilibrium point $E{P}_1=(10,0,0,0,251.56,196.97,18.5)$ for all initials IC1-IC3. In this situation, the patient becomes infected by HTLV-I, while the SARS-CoV-2 infection is cleared.

Case 3.

($R_2>1$ and $R_3\le 1$): We choose $\rho =3,$ $\mu =0.01,$${\xi }_V=0.04$ and $\pi =0.0001$. Then, we calculate $R_2=15.3786>1$ and $R_3=0.2407<1$. Lemma 2 and Theorem 3 state that the SARS-CoV-2-mono-infection equilibrium point, $E{P}_2=\left(0.702,0.022,0.009,0.049,900.46,0,0\right)$ exists and is GAS. Figure 3 displays the numerical solutions of the system converge to $E{P}_2$ for all the three initials IC1-IC3. The results support the theoretical results presented in Theorem 3. This situation represents a patient who infected only with SARS-CoV-2.

Case 4.

($R_4>1$ and $1<R_3\le 1+\tilde{R}$): We consider $\rho =3,$$\mu =0.01,$${\xi }_V=0.6$ and $\pi =0.0015$. So, we obtain $R_4=3.2969>1,$$R_3=3.3113>1$ and $R_3<1+\tilde{R}=4.0299$. Lemma 2 and Theorem 4 state that the HTLV-I/SARS-CoV-2 coinfection $E{P}_3=(3.03,0.017,0.023,0.008,251.56,$$213.49,20.05)$ exists and is GAS. Figure 4 shows that the solutions of the system converge to $E{P}_3$ for all initials IC1-IC3. The results support the theoretical results presented in Theorem 4. In this case, a SARS-CoV-2 and HTLV-I coinfection happens, where a HTLV-I-patient gets contaminated with COVID-19. CD4${}^{+}$T cells are animated to dispense with SARS-CoV-2 disease from the body. In any case, assuming the patient has low CD4${}^{+}$T cell counts, the freedom of SARS-CoV-2 may not be accomplished. This can prompt extreme contamination and passing.

4.2. Impact of Time Delays on the Clearance of Coinfection

In this part, we study the effect of time delay parameters ${\tau }_i,i=1,2,\dots ,5$ on the stability of the uninfected equilibrium point $E{P}_0$. Since $R_1$ and $R_2$ depend on ${\tau }_i,$ $i=1,$$2,\dots ,5,$ then changing the parameters ${\tau }_i$ will change the stability of the uninfected equilibrium point. Using the values of the parameters in Table 2 and considering $\rho =3,$ $\mu =0.01,$ ${\xi }_V=0.06$ and $\pi =0.0015$to solve the model (7)–(13) with initial condition:

$$\mathrm{IC}4:\left(X\right(\theta ),N(\theta ),Y(\theta ),V(\theta ),U(\theta ),L(\theta ),A(\theta \left)\right)=(5,1,0.0012,0.004,400,200,20),$$

where $\theta \in \left[-max{\tau }_i,0\right],i=1,2,\dots ,5.$ Consider $\tau ={\tau }_1={\tau }_2={\tau }_3={\tau }_4={\tau }_5$ and the other parameters are fixed, then $R_1$ and $R_2$ can be given as functions of $\tau$ as follows:

$$R_1\left(\tau \right)=\frac{\pi \alpha \gamma e^{-\left(d_4+d_5\right)\tau }}{{\xi }_U\left({\xi }_L{\xi }_A+\alpha \left({\xi }_A-\epsilon e^{-d_5\tau }\right)\right)},R_2\left(\tau \right)=\frac{\delta \eta \rho \kappa {\xi }_U{e}^{-\left(d_1+d_2+d_3\right)\tau }}{{\xi }_X{\xi }_V\left(\kappa +{\xi }_N\right)\left({\xi }_Y{\xi }_U+\gamma \mu \right)}.$$

We note that, when all other parameters are fixed, then $R_1$ and $R_2$ are decreasing functions of $\tau$. Let ${\tau }_{c{r}_1}$ and ${\tau }_{c{r}_2}$ be such that $R_1\left({\tau }_{c{r}_1}\right)=R_2\left({\tau }_{c{r}_2}\right)=1$. Consequently,

$$\begin{array}{cc}\hfill R_1\left(\tau \right)& \le 1,\mathrm{for}\mathrm{all}\tau \ge {\tau }_{c{r}_1},\hfill \\ \hfill R_2\left(\tau \right)& \le 1,\mathrm{for}\mathrm{all}\tau \ge {\tau }_{c{r}_2}.\hfill \end{array}$$

Therefore, $E{P}_0$ is GAS when $\tau \ge {\tau }_{c{r}_1}$ and $\tau \ge {\tau }_{c{r}_2}$. Using the values of the parameters, we get, ${\tau }_{c{r}_1}=0.692433$ and ${\tau }_{c{r}_2}=0.108309$. It follows that:

(i) If $\tau \ge max\{{\tau }_{c{r}_1},{\tau }_{cr2}\}=0.692433$, then $R_1\left(\tau \right)\le 1$ and $R_2\left(\tau \right)\le 1$ and thus, $E{P}_0$ is GAS.

(ii) If $\tau <0.692433,$ then $R_1\left(\tau \right)>1$ or $R_2\left(\tau \right)>1$ and thus, $E{P}_0$ will lose its stability.

Therefore, sufficiently large time delay can stabilize the system around the equilibrium $E{P}_0$.

It is clear from Figure 5 that, as $\tau$ increases, the concentrations of the uninfected ECs and uninfected CD4${}^{+}$T cells are increased, while the other concentrations of other compartments are decreased. We can see from the above discussion that increasing time delays values can have similar effect as antiviral treatments.

5. Discussion

Coinfection cases with SARS-CoV-2 and HTLV-I were reported in [3,10]. Therefore, it is important to understand the within-host dynamics of this coinfection. In this paper, we developed and examined two SARS-CoV-2 and HTLV-I coinfection models. The models considered the interactions between uninfected ECs, latently SARS-CoV-2-infected ECs, actively SARS-CoV-2-infected ECs, free SARS-CoV-2 particles, uninfected CD4${}^{+}$T cells, latently HTLV-I-infected CD4${}^{+}$T cells and actively HTLV-I-infected CD4${}^{+}$T cells. We introduced five intracellular time delays into the models: (i) two delays in the formation of latently SARS-CoV-2-infected ECs and latently HTLV-I-infected CD4${}^{+}$T cells, (ii) two delays in the reactivation of latently SARS-CoV-2-infected ECs and latently HTLV-I-infected CD4${}^{+}$T cells, and (iii) maturation delay of new SARS-CoV-2 virions. Discrete-time delays and distributed-time delays were incorporated in the first and second models, respectively. We examined the nonnegativity and ultimate boundedness of the solutions. We found that these systems has four equilibria and we proved the following:

(I) The uninfected equilibrium point $E{P}_0$ always exists. It is GAS when $R_1\le 1$ and $R_2\le 1$. In this case, the patient is recovered from both SARS-CoV-2 and HTLV-I. From a control viewpoint, making $R_1\le 1$ and $R_2\le 1$ will be a good strategy. The parameter $R_2$ may be reduced by multiplying the parameters by $\rho$ or $\eta$ by $(1-{ϵ}_1)$ or $(1-{ϵ}_2)$, respectively. Here, ${ϵ}_1\in [0,1]$ and ${ϵ}_2\in [0,1]$ are the effectiveness of the antiviral drugs for blocking the infection and blocking the production of SARS-CoV-2 particles, respectively [20]. Since there is no treatment for HTLV-I infection [11], making $R_1\le 1$ is rarely achieved.

(II) The HTLV-I mono-infection equilibrium point, $E{P}_1$ exists if $R_1>1$. It is GAS when $R_1>1$ and $R_4\le 1$. This case leads to the situation of the patient who infected by HTLV-I, while the SARS-CoV-2 infection is cleared.

(III) The SARS-CoV-2 mono-infection equilibrium point, $E{P}_2$ exists if $R_2>1$. It is GAS when $R_2>1$ and $R_3\le 1$. This case leads to the situation of the patient who has SARS-CoV-2 mono-infection.

(IV) The HTLV-I and SARS-CoV-2 coinfection equilibrium point, $E{P}_3$ exists if $R_3>1$ and $R_4>1$. It is GAS when $R_4>1$ and $1<R_3\le 1+\tilde{R}$. This point represents the situation of SARS-CoV-2 and HTLV-I coinfection patient.

The global stability of equilibria was established using Lyapunov method. We performed numerical simulations and demonstrated that they are in good agreement with the theoretical results. We discussed the impact of time delays on the clearance of coinfection. We concluded that, increasing time delays values can have an antiviral treatment-like impact.

The models developed in this work can be improved by (i) utilizing real data to find a good estimation of the parameters’ values, thus, the theoretical results obtained in this paper need to be tested against empirical findings when real data become available, (ii) considering the mutations of the SARS-CoV-2, (iii) including the stochastic interaction, and (iv) considering the diffusion of cells and viruses. Future studies should evaluate the effect of SARS-CoV-2 on HTLV-1 patients. These research points need further investigations so we leave them to future works.