Mathematical modeling of HIV/HTLV co-infection with CTL-mediated immunity

1.   Introduction

During the last decades different dangerous viruses have been recognized which attack the human body and causes many fatal diseases. As an example of these viruses, the human immunodeficiency virus (HIV) which is the causative agent for acquired immunodeficiency syndrome (AIDS). According to global health observatory (GHO, 2018) data of HIV/AIDS published by WHO ^[1] that says, globally, about 37.9 million HIV-infected people in 2018, 1.7 million newly HIV-infected and 770,000 HIV-related death in the same year. HIV is a retrovirus that infects the susceptible CD$4^{+}$T cells which play a central role in immune system defence. During the last decades, mathematical modeling of within-host HIV infection has witnessed a significant development ^[2]. Nowak and Bangham ^[3] have introduced the basic HIV infection model which describes the interaction between three compartments, susceptible CD$4^{+}$T cells ($S$), actively HIV-infected cells ($I$) and free HIV particles ($V$). Latent viral reservoirs remain one of the major hurdles for eradicating the HIV by current antiviral therapy. Latently HIV-infected cells include HIV virions but do not produce them until they become activated. Mathematical modeling of HIV dynamics with latency can help in predicting the effect of antiviral drug efficacy on HIV progression ^[4]. Rong and Perelson ^[5] have incorporated the latently infected cells in the basic HIV model presented in ^[3] as:

where $S = S(t),$ $L = L(t)$, $I = I(t)\ $and $V = V(t)$ are the concentrations of susceptible CD$4^{+}$T cells, latently HIV-infected cells, actively HIV-infected cells and free HIV particles at time $t$, respectively. The susceptible CD$4^{+}$T cells are produced at specific constant rate $\rho$. The HIV virions can replicate using virus-to-cell (VTC) transmission. The term $\eta_{1}SV$ refers to the rate at which new infectious appears by VTC contact between free HIV particles and susceptible CD$4^{+}$T cells. Latently HIV-infected cells are transmitted to be active at rate $\lambda L$. The free HIV particles are generated at rate $bI$. The natural death rates of the susceptible CD$4^{+}$T cells, latently HIV-infected cells, actively HIV-infected cells and free HIV particles are given by $\alpha S$, $\gamma L,$ $aI$ and $\varepsilon V$, respectively. A fraction $\beta \in(0, 1)$ of new HIV-infected cells will be active, and the remaining part $1-\beta$ will be latent. During the last decades, mathematical modeling and analysis of HIV mono-infection with both latently and actively HIV-infected cells have witnessed a significant development ^{[6,7,8,9,10,11,12]}.

Model (1.1) assumed that the HIV can only spread by VTC transmission. However, several works have reported that there is another mode of transmission called cell-to-cell (CTC) where the HIV can be transmitted directly from an infected cell to a healthy CD$4^{+}$T cell through the formation of virological synapses ^[13]. Sourisseau et al. ^[14] have shown that CTC transmission plays an efficient role in the HIV replication. Sigal et al. ^[15] have demonstrated the importance of CTC transmission in the HIV infection process during the antiviral treatment. Iwami et al. ^[13] have shown that about 60% of HIV infections are due to CTC transmission. In addition, CTC transmission can increase the HIV fitness by 3.9 times and decrease the production time of HIV particles by 0.9 times ^[16]. HIV dynamics model with latency and both VTC and CTC transmissions is given by ^[17,18]:

where, the term $\eta_{2}SI$ refers to the rate at which new infectious appears by CTC contact between HIV-infected cells and susceptible CD$4^{+}$T cells.

Another example of the dangerous human viruses is called Human T-lymphotropic virus type Ⅰ (HTLV-I) which can lead to two diseases, adult T-cell leukemia (ATL) and HTLV-I-associated myelopathy/tropical spastic paraparesis (HAM/TSP). The discovery of the first human retrovirus HTLV-I is back to 1980, and after 3 years the HIV was determined ^[19]. HTLV-I is global epidemic that infects about 10-25 million persons ^[20]. The infection is endemic in the Caribbean, southern Japan, the Middle East, South America, parts of Africa, Melanesia and Papua New Guinea ^[21]. HTLV-I is a provirus that targets the susceptible CD$4^{+}$T cells. HTLV-I can spread to susceptible CD$4^{+}$T cells from CTC through the virological synapse. HTLV-infected cells can be divided into two kinds based on the presence of Tax inside the cell or not: (i) Tax$^{\bf{-}}$, or latently HTLV-infected cells are resting CD4$^{+}$T cells that contain a provirus and do not express Tax, and (ii) Tax$^{\bf{+}}$, or actively HTLV-infected cells are activated provirus-carrying CD4$^{\bf{+}}$T cells that do express Tax ^[22]. During the primary infection stage of HTLV-I, the proviral load can reach high level, approximately 30-50% ^[23]. Unlike in the case of HIV infection, however, only a small percentage of infected individuals develop the disease and 2-5% percent of HTLV-I carriers develop symptoms of ATL and another 0.25-3% develop HAM/TSP ^[24]. Stilianakis and Seydel ^[25] have formulated an HTLV-I model to describe the interaction of susceptible CD$4^{+}$T cells, latently HTLV-infected cells, Tax-expressing HTLV-infected cells (actively HTLV-infected cells) and leukemia cells (ATL cells) as:

where $S = S(t),$ $E = E(t),$ $Y = Y(t)$ and $Z = Z(t)$ are the concentrations of susceptible CD$4^{+}$T cells, latently HTLV-infected cells, Tax-expressing HTLV-infected cells and ATL cells, at time $t$, respectively. In contrast of HIV, the transmission of HTLV-I can be only from CTC that is the HTLV virions can only survive inside the host CD$4^{+}$T cells and cannot be detectable in the plasma. The rate at which new infectious appears by CTC contact between Tax-expressing HTLV-infected cells and susceptible CD$4^{+}$T cells is assumed to be $\eta_{3}SY$. The natural death rate of the latently HTLV-infected cells, Tax-expressing HTLV-infected cells and ATL cells are represented by $\omega E$, $\delta Y$ and $\theta Z$, respectively. The term $\psi E$ accounts for the rate of latently HTLV-infected cells that become Tax-expressing HTLV-infected cells. $\vartheta Y$ is the transmission rate at which Tax-expressing HTLV-infected cells convert to ATL cells. The logistic term $\ell Z\left(1-\frac{Z}{Z_{\max}}\right)$ denotes the proliferation rate of the ATL cells, where $Z_{\max}$ is the maximal concentration that ATL cells can grow. The parameter $\ell$ is the maximum proliferation rate constant of ATL cells. Many researchers have been concerned to study mathematical modeling and analysis of HTLV-I mono-infection in several works ^[26,27,28].

Cytotoxic T lymphocytes (CTLs) are recognized as the significant component of the human immune response against viral infections. CTLs inhibit viral replication and kill the cells which are infected by viruses. In fact, CTLs are necessary and universal to control HIV infection ^[29]. During the recent years, great efforts have been made to formulate and analyze the within-host HIV mono-infection models under the influence of CTL immune response (see e.g. ^[2,3]). In ^[30], latently HIV-infected cells have been included in the HIV dynamics models with CTL immune response. In case of HTLV-I infection, it has been reported in ^[31] that the CTLs play an effective role in controlling such infection. CTLs can recognize and kill the Tax-expressing HTLV-infected cells, moreover, they can reduce the proviral load. In the literature, several mathematical models have been proposed to describe the dynamics of HTLV-I under the effect of CTL immune response (see e.g. ^{[21,32,33,34,35,36]}). In ^[20,37,38], HTLV-I dynamics models have been presented by incorporating latently HTLV-infected cells and CTL immune response.

Simultaneous infection by HIV and HTLV-I and the etiology of their pathogenic and disease outcomes have become a global health matter over the past 10 years ^[39]. It is commonly that HIV/HTLV-I co-infection can be endemic in areas where individuals experience high risk attitudes; such as unprotected sexual contact and unsafe injection practices; that cause transmission of contaminated body fluids between individuals. This shed a light on the importance of studying HIV/HTLV-I co-infection ^[40]. Although CD$4^{+}$T cells are the major targets of both HIV and HTLV-I, however, these viruses present a different biological behavior that causes diverse impacts on host immunity and ultimately lead to numerous clinical diseases ^[41]. It has been reported that the HTLV-I co-infection rate among HIV infected patients as increase as 100 to 500 times in comparison with the general population ^[42]. In seroepidemiologic studies, it has been recorded that HIV-infected patients are more exposure to be co-infected with HTLV-I, and vice versa compared to the general population ^[43]. HIV/HTLV-I co-infection is usually found in individuals of specific ethnic or who belonged to geographic origins where these viruses are simultaneously endemic ^[44]. As an example, the co-infection rates in individuals living in Bahia have reached 16% of HIV-infected patients ^[45]. The prevalence of dual infection with HIV and HTLV-I has become more widely in several geographical regions throughout the world such as South America, Europe, the Caribbean, Bahia (Brazil), Mozambique (Africa), and Japan ^{[39,43,45,46,47]}. HIV and HTLV-I dual infection appears to have an overlap on the course of associated clinical outcomes with both viruses ^[43]. Several reports have concluded that HIV/HTLV-I co-infected patients were found to have an increase of CD$4^{+}$T cells count in comparison with HIV mono-infected patients, although there is no evident to result in a better immune response ^[41,48]. Indeed, simultaneous infected patients by both viruses with CD$4^{+}$T counts greater than 200 cells/mm$^{3}$ are more exposure to have other opportunistic infections as compared with HIV mono-infected patients who have similar CD$4^{+}$T counts ^[48]. Studies have reported that higher mortality and shortened survival rates were accompany with co-infected individuals more than mono-infected individuals ^[46]. Considering the natural history of HIV, many researchers have noted that co-infection with HIV and HTLV-I can accelerate the clinical progression to AIDS. On the other hand, HIV can adjust HTLV-I expression in co-infected individuals which leads them to a higher risk of developing HTLV-I related diseases such as ATL and TSP/HAM ^[42,43,46].

Great efforts have been made to develop and analyze mathematical models of HIV and HTLV-I mono-infections, however, modeling of HIV/HTLV-I co-infection has not been studied. In fact, such co-infection modeling and its analysis will be needed to help clinicians on estimating the appropriate time to initiate treatment in co-infected patients. Therefore, the aim of the present paper is to formulate a new HIV/HTLV-I co-infection model. We show that the model is well-posed by establishing that the solutions of the model are nonnegative and bounded. We derive a set of threshold parameters which govern the existence and stability of the equilibria of the model. Global stability of all equilibria is proven by constructing suitable Lyapunov functions and utilizing Lyapunov-LaSalle asymptotic stability theorem. We conduct some numerical simulations to illustrate the theoretical results.

The results of this work, such as co-infection model and its analysis will help clinicians estimate the appropriate time for patients with co-infection to begin treatment. On the other hand, this study, from a certain point of view, illustrate the complexity of this co-infection model and the model is helpful to clinic treatment. It is worth mentioning, if we look at research perspectives, that appropriate developments of the model presented in this paper can be focused on the within host modeling of the competition between COVID19 virus and the immune system by a complex dynamics described in ^[49]. This dynamics which occurs, in human lungs, once the virus, after contagion, has gone over the biological barriers which protect each individuals, see ^[47].

2.   Model formulation

We set up an ordinary differential equation model that describes the change of concentrations of eight compartments with respect to time $t$; susceptible CD$4^{+}$T cells $S(t)$, latently HIV-infected cells $L(t)$, actively HIV-infected cells $I(t)$, latently HTLV-infected cells $E(t)$, Tax-expressing HTLV-infected cells $Y(t)$, free HIV particles $V(t)$, HIV-specific CTLs $C^{I}(t)$ and HTLV-specific CTLs $C^{Y}(t)$. The dynamics of HIV/HTLV-I co-infection is schematically shown in the transfer diagram given in Figure 1. Our proposed model is given by the following form:

where $(S, L, I, E, Y, V, C^{I}, C^{Y}) = (S(t), L(t), I(t), E(t), Y(t), V(t), C^{I} (t), C^{Y}(t))$. The term $\mu_{1}C^{I}I$ is the killing rate of actively HIV-infected cells due to their specific immunity. The term $\mu_{2}C^{Y}Y$ is the killing rate of Tax-expressing HTLV-infected cells due to their specific immunity. The proliferation and death rates for both effective HIV-specific CTLs and HTLV-specific CTLs are given by $\sigma_{1}C^{I}I$, $\sigma_{2} C^{Y}Y$, $\pi_{1}C^{I}$ and $\pi_{2}C^{Y}$, respectively. All remaining parameters have the same biological meaning as explained in the previous section. All parameters and their definitions are summarized in Table 1.

3.   Preliminaries

Let $\Omega_{j} > 0$, $j = 1, ..., 5$ and define

Proposition 1. The compact set $\Theta$ is positively invariant for system (2.1).

Proof. We have

This ensures that, $(S(t), L(t), I(t), E(t), Y(t), V(t), C^{I}(t), C^{Y} (t))\in \mathbb{R}_{\geq0}^{8}$ for all $t\geq0$ when $(S(0), L(0), I(0), E(0), Y(0), V(0), C^{I}(0), C^{Y}(0))\in \mathbb{R}_{\geq0}^{8}$. To show the boundedness of all state variables, we let

Then

where $\phi = \min \{ \alpha, \gamma, \frac{a}{2}, \omega, \delta, \varepsilon, \pi _{1}, \pi_{2}\}$. Hence, $0\leq \Psi(t)\leq \Omega_{1}$ if $\Psi(0)\leq \Omega _{1}$ for $t\geq0,$ where $\Omega_{1} = \frac{\rho}{\phi}.$ Since $S$, $L,$ $I,$ $E,$ $Y,$ $V,$ $C^{I}$, and $C^{Y}$ are all nonnegative then $0\leq S(t), L(t), I(t)\leq \Omega_{1}$, $0\leq E(t), Y(t)\leq \Omega_{2}$, $0\leq V(t)\leq \Omega_{3}$, $0\leq C^{I}(t)\leq \Omega_{4}$, $0\leq C^{Y}(t)\leq \Omega_{5}$ if $S(0)+L(0)+I(0)+\frac{1}{\varphi}\left(E(0)+Y(0)\right) +\frac{a}{2b}V(0)+\frac{\mu_{1}}{\sigma_{1}}C^{I}(0)+\frac{\mu_{2}} {\varphi \sigma_{2}}C^{Y}(0)\leq \Omega_{1}$, where $\Omega_{2} = \varphi \Omega_{1}$, $\Omega_{3} = \dfrac{2b\Omega_{1}}{a}$, $\Omega_{4} = \dfrac {\sigma_{1}\Omega_{1}}{\mu_{1}}$ and $\Omega_{5} = \dfrac{\varphi \sigma _{2}\Omega_{1}}{\mu_{2}}$.

4.   Equilibria

In this section, we derive eight threshold parameters which guarantee the existence of the equilibria of the model. Let $(S, L, I, E, Y, V, C^{I}, C^{Y})$ be any equilibrium of system (2.1) satisfying the following equations:

The straightforward calculation finds that system (2.1) admits eight equilibria.

(ⅰ) Infection-free equilibrium, $Đ _{0} = (S_{0}, 0, 0, 0, 0, 0, 0, 0)$, where $S_{0} = \rho/\alpha$. This case describes the situation of healthy state where both HIV and HTLV are absent.

(ⅱ) Chronic HIV mono-infection equilibrium with inactive immune response, $Đ _{1} = (S_{1}, L_{1}, I_{1}, 0, 0, V_{1}, 0, 0),$ where

Therefore, $Đ _{1}$ exists when

At the equilibrium $Đ _{1}$ the chronic HIV mono-infection persists while the immune response is unstimulated. The basic HIV mono-infection reproductive ratio for system (2.1) is defined as:

where

The parameter $\Re_{1}$ determines whether or not a chronic HIV infection can be established. In fact, $\Re_{11}$ measures the average number of secondary HIV infected generation caused by an existing free HIV particles, while $\Re_{12}$ measures the average number of secondary HIV infected generation caused by an HIV-infected cell. Therefore, $\Re_{11}$ and $\Re_{12}$ are the basic HIV mono-infection reproductive ratio corresponding to VTC and CTC infections, respectively. In terms of $\Re_{1}$, we can write

(ⅲ) Chronic HTLV mono-infection equilibrium with inactive immune response, $Đ _{2} = (S_{2}, 0, 0, E_{2}, Y_{2}, 0, 0, 0)$, where

Therefore, $Đ _{2}$ exists when

At the equilibrium $Đ _{2}$ the chronic HTLV mono-infection persists while the immune response is unstimulated. The basic HTLV mono-infection reproductive ratio for system (2.1) is defined as:

The parameter $\Re_{2}$ decides whether or not a chronic HTLV infection can be established. In terms of $\Re_{2}$, we can write

Remark 1. We note that both $\Re_{1}$ and $\Re_{2}$ does not depend of parameters $\sigma_{i}$, $\pi_{i}$ and $\mu_{i}$, $i = 1, 2$. Therefore, without treatment CTLs will not able to clear HIV or HTLV-I from the body.

(ⅳ) Chronic HIV mono-infection equilibrium with only active HIV-specific CTL, $Đ _{3} = (S_{3}, L_{3}, I_{3}, 0, 0, V_{3}, C_{3}^{I}, 0)$, where

We note that $Đ _{3}$ exists when $\dfrac{\sigma_{1}\rho \left(\beta \gamma+\lambda \right) \left(\eta_{1}b+\eta_{2}\varepsilon \right) }{a\left(\gamma+\lambda \right) \left[ \pi_{1}\left(\eta_{1}b+\eta _{2}\varepsilon \right) +\alpha \varepsilon \sigma_{1}\right] } > 1$. The HIV-specific CTL-mediated immunity reproductive ratio in case of HIV mono-infection is stated as:

Thus, $C_{3}^{I} = \dfrac{a}{\mu_{1}}(\Re_{3}-1)$. The parameter $\Re_{3}$ determines whether or not the HIV-specific CTL-mediated immune response is stimulated in the absent of HTLV infection.

(ⅴ) Chronic HTLV mono-infection equilibrium with only active HTLV-specific CTL, $Đ _{4} = (S_{4}, 0, 0, E_{4}, Y_{4}, 0, 0, C_{4}^{Y})$, where

We note that $Đ _{4}$ exists when $\dfrac{\sigma_{2}\rho \varphi \eta_{3}\psi }{\delta \left(\psi+\omega \right) \left(\pi_{2}\eta_{3}+\alpha \sigma _{2}\right) } > 1$. The HTLV-specific CTL-mediated immunity reproductive ratio in case of HTLV mono-infection is stated as:

Thus, $C_{4}^{Y} = \dfrac{\delta}{\mu_{2}}(\Re_{4}-1)$. The parameter $\Re_{4}$ determines whether or not the HTLV-specific CTL-mediated immune response is stimulated in the absent of HIV infection.

(ⅵ) Chronic HIV/HTLV co-infection equilibrium with only active HIV-specific CTL, $Đ _{5} = (S_{5}, L_{5}, I_{5}, E_{5}, Y_{5}, V_{5}, C_{5}^{I}, 0)$, where

We note that $Đ _{5}$ exists when $\Re_{1}/\Re_{2} > 1$ and $\frac{\rho \varphi \varepsilon \eta_{3}\sigma_{1}\psi}{\delta \left(\psi+\omega \right) \left[ \pi_{1}\left(\eta_{1}b+\eta_{2}\varepsilon \right) +\alpha \varepsilon \sigma_{1}\right] } > 1$. The HTLV infection reproductive ratio in the presence of HIV infection is stated as:

Thus, $E_{5} = \dfrac{\delta \left[ \pi_{1}\left(\eta_{1}b+\eta_{2} \varepsilon \right) +\alpha \varepsilon \sigma_{1}\right] }{\varepsilon \eta _{3}\sigma_{1}\psi}\left(\Re_{5}-1\right),$ $Y_{5} = \dfrac{\pi_{1}\left(\eta_{1}b+\eta_{2}\varepsilon \right) +\alpha \varepsilon \sigma_{1} }{\varepsilon \eta_{3}\sigma_{1}}\left(\Re_{5}-1\right)$. The parameter $\Re_{5}$ determines whether or not HIV-infected patients could be co-infected with HTLV.

(ⅶ) Chronic HIV/HTLV co-infection equilibrium with only active HTLV-specific CTL, $Đ _{6} = (S_{6}, L_{6}, I_{6}, E_{6}, Y_{6}, V_{6}, 0, C_{6}^{Y})$, where

We note that $Đ _{6}$ exists when $\Re_{2}/\Re_{1} > 1$ and $\frac{\rho \sigma_{2}\left(\beta \gamma+\lambda \right) \left(\eta_{1}b+\eta _{2}\varepsilon \right) }{a\varepsilon \left(\gamma+\lambda \right) \left(\pi_{2}\eta_{3}+\alpha \sigma_{2}\right) } > 1$. The HIV infection reproductive ratio in the presence of HTLV infection is stated as:

Thus, $L_{6} = \dfrac{a\varepsilon \left(1-\beta \right) \left(\pi_{2} \eta_{3}+\alpha \sigma_{2}\right) }{\sigma_{2}\left(\beta \gamma +\lambda \right) \left(\eta_{1}b+\eta_{2}\varepsilon \right) }\left(\Re_{6}-1\right),$ $I_{6} = \dfrac{\varepsilon \left(\pi_{2}\eta_{3} +\alpha \sigma_{2}\right) }{\sigma_{2}\left(\eta_{1}b+\eta_{2} \varepsilon \right) }\left(\Re_{6}-1\right)$, $V_{6} = \dfrac{b\left(\pi_{2}\eta_{3}+\alpha \sigma_{2}\right) }{\sigma_{2}\left(\eta_{1} b+\eta_{2}\varepsilon \right) }\left(\Re_{6}-1\right)$. The parameter $\Re_{6}$ determines whether or not HTLV-infected patients could be co-infected with HIV.

(ⅷ) Chronic HIV/HTLV co-infection equilibrium with active HIV-specific CTL and HTLV-specific CTL, $Đ _{7} = (S_{7}, L_{7}, I_{7}, E_{7}, Y_{7}, V_{7}, C_{7}^{I}, C_{7}^{Y})$, where

It is obvious that $Đ _{7}$ exists when $\dfrac{\sigma_{1}\sigma_{2} \rho \left(\beta \gamma+\lambda \right) \left(\eta_{1}b+\eta_{2} \varepsilon \right) }{a\left(\gamma+\lambda \right) \left[ \pi_{1} \sigma_{2}\left(\eta_{1}b+\eta_{2}\varepsilon \right) +\pi_{2}\eta _{3}\varepsilon \sigma_{1}+\alpha \varepsilon \sigma_{1}\sigma_{2}\right] } > 1$ and $\dfrac{\varphi \eta_{3}\varepsilon \sigma_{1}\sigma_{2}\rho \psi} {\delta \left(\psi+\omega \right) \left[ \pi_{1}\sigma_{2}\left(\eta _{1}b+\eta_{2}\varepsilon \right) +\pi_{2}\eta_{3}\varepsilon \sigma_{1} +\alpha \varepsilon \sigma_{1}\sigma_{2}\right] } > 1$. Now we define

Clearly, $Đ _{7}$ exists when $\Re_{7} > 1$ and $\Re_{8} > 1$ and we can write $C_{7}^{I} = \dfrac{a}{\mu_{1}}\left(\Re_{7}-1\right)$ and $C_{7}^{Y} = \dfrac{\delta}{\mu_{2}}\left(\Re_{8}-1\right).$ The parameter $\Re_{7}$ refers to the competed HIV-specific CTL-mediated immunity reproductive ratio in case of HIV/HTLV co-infection. On the other hand, the parameter $\Re_{8}$ refers to the competed HTLV-specific CTL-mediated immunity reproductive ratio case of HIV/HTLV co-infection.

The eight threshold parameters are given as follows:

According to the above discussion, we sum up the existence conditions for all equilibria in Table 2.

5.   Global stability analysis

In this section we prove the global asymptotic stability of all equilibria by constructing Lyapunov functional following the method presented in ^[50]. We will use the arithmetic-geometric mean inequality

which yields

Theorem 1. If $\Re_{1}\leq1$ and $\Re_{2}\leq1$, then $Đ _{0}$ is globally asymptotically stable (G.A.S).

Proof. We construct a Lyapunov function $\Phi_{0}(S, L, I, E, Y, V, C^{I}, C^{Y})$ as:

where,

Clearly, $\Phi_{0}(S, L, I, E, Y, V, C^{I}, C^{Y}) > 0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$, and $\Phi_{0}(S_{0}, 0, 0, 0, 0, 0, 0, 0) = 0$. Calculating $\frac{d\Phi _{0}}{dt}$ along the solutions of system (2.1) as:

Using $S_{0} = \rho/\alpha$, we obtain

Therefore, $\frac{d\Phi_{0}}{dt}\leq0$ for all $S, I, Y, C^{I}, C^{Y} > 0$ and $\frac{d\Phi_{0}}{dt} = 0$ when $S = S_{0}$ and $I = Y = C^{I} = C^{Y} = 0.$ Define $\Upsilon_{0} = \left \{ (S, L, I, E, Y, V, C^{I}, C^{Y}):\frac {d\Phi_{0}}{dt} = 0\right \}$ and let $\Upsilon_{0}^{\prime}$ be the largest invariant subset of $\Upsilon_{0}$. The solutions of system (2.1) converge to $\Upsilon_{0}^{^{\prime}}$. The set $\Upsilon _{0}^{^{\prime}}$ includes elements with $S = S_{0}$ and $I = Y = C^{I} = C^{Y} = 0$, and hence $\dot{S} = \dot{Y} = 0$. The first and fifth equations of system (2.1) yield

Thus, $V(t) = E(t) = 0$ for all $t$. In addition, we have $\dot{I} = 0$ and from the third equation of system (2.1) we obtain

which yields $L(t) = 0$ for all $t$. Therefore, $\Upsilon_{0}^{\prime} = \left \{ Đ _{0}\right \}$ and by applying Lyapunov-LaSalle asymptotic stability theorem ^[51,52,53] we get that $Đ _{0}$ is G.A.S.

Theorem 2. Let $\Re_{1} > 1$, $\Re_{2}/\Re_{1}\leq1$ and $\Re_{3}\leq 1$, then $Đ _{1}$ is G.A.S.

Proof. Define a function $\Phi_{1}(S, L, I, E, Y, V, C^{I}, C^{Y})$ as:

Calculating $\frac{d\Phi_{1}}{dt}$ as:

Using the equilibrium conditions for $Đ _{1}$, we get

Then, we obtain

Therefore, Eq (5.5) becomes

Since $\Re_{2}/\Re_{1}\leq1$ and $\Re_{3}\leq1$, then using inequalities (5.1)–(5.3) we get $\frac{d\Phi_{1}}{dt}\leq0$ for all $S, L, I, Y, V, C^{I}, C^{Y} > 0$. Moreover, $\frac{d\Phi_{1}}{dt} = 0$ when $S = S_{1}$, $L = L_{1}$, $I = I_{1}$, $V = V_{1}$ and $Y = C^{I} = C^{Y} = 0.$ The solutions of system (2.1) converge to $\Upsilon_{1}^{\prime}$ the largest invariant subset of $\Upsilon_{1} = \left \{ (S, L, I, E, Y, V, C^{I}, C^{Y}):\frac{d\Phi_{1}}{dt} = 0\right \}$. The set $\Upsilon_{1}^{\prime}$ includes $Y = 0$, and then $\dot{Y} = 0$. The fifth equation of system (2.1) implies

which yields $E(t) = 0$ for all $t.$ Hence, $\Upsilon_{1}^{\prime} = \{Đ_1\}$ and $Đ _1$ is G.A.S using Lyapunov-LaSalle asymptotic stability theorem.

Theorem 3. If $\Re_{2} > 1$, $\Re_{1}/\Re_{2}\leq1$ and $\Re_{4}\leq1$, then $Đ _{2}$ is G.A.S.

Proof. We define $\Phi_{2}(S, L, I, E, Y, V, C^{I}, C^{Y})$ as:

We calculate $\frac{d\Phi_{2}}{dt}$ as:

Using the equilibrium conditions for $Đ _{2}$:

we obtain

Since $\Re_{1}/\Re_{2}\leq1$ and $\Re_{4}\leq1$, then using inequality (5.4) we get $\frac{d\Phi_{2}}{dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$. In addition, $\frac{d\Phi_{2}}{dt} = 0$ when $S = S_{2},$ $E = E_{2},$ $Y = Y_{2}$ and $I = C^{I} = C^{Y} = 0.$ Define $\Upsilon_{2} = \left \{ (S, L, I, E, Y, V, C^{I}, C^{Y}):\frac{d\Phi_{2}}{dt} = 0\right \}$ and let $\Upsilon_{2}^{\prime}$ be the largest invariant subset of $\Upsilon_{2}$. The solutions of system (2.1) converge to $\Upsilon_{2}^{\prime}$ which includes elements with $S = S_{2},$ $Y = Y_{2},$ $I = 0$, then $\dot{S} = 0$. The first equation of system (2.1) gives

From conditions (5.7) we get $V(t) = 0$ for all $t$. Moreover, we have $\dot{I} = 0$ and from the third equation of system (2.1) we obtain

which yields $L(t) = 0$ for all $t$. Therefore, $\Upsilon_{2}^{\prime} = \{ Đ _2 \}$. By applying Lyapunov-LaSalle asymptotic stability theorem we get that $Đ _{2}$ is G.A.S.

Theorem 4. Let $\Re_{3} > 1$ and $\Re_{5}\leq1$, then $Đ _{3}$ is G.A.S.

Proof. Define a function $\Phi_{3}(S, L, I, E, Y, V, C^{I}, C^{Y})$ as:

We calculate $\frac{d\Phi_{3}}{dt}$ as:

Using the equilibrium conditions for $Đ _{3}$:

we obtain

Hence, if $\Re_{5}\leq1$, then using inequalities (5.1)–(5.3) we get $\frac{d\Phi_{3}}{dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$. Moreover, $\frac{d\Phi_{3}}{dt} = 0$ at $S = S_{3}$, $L = L_{3},$ $I = I_{3},$ $V = V_{3}$ and $Y = C^{Y} = 0.$ The solutions of system (2.1) converge to $\Upsilon_{3}^{\prime}$ the largest invariant subset of $\Upsilon_{3} = \left \{ (S, L, I, E, Y, V, C^{I}, C^{Y}):\frac{d\Phi_{3}} {dt} = 0\right \}$. The set $\Upsilon_{3}^{\prime}$ contains elements with $S = S_{3}$, $L = L_{3},$ $I = I_{3},$ $V = V_{3},$ $Y = 0$, and then $\dot{I} = \dot{Y} = 0$. The third and fifth equations of system (2.1) give

which yield $C^{I}(t) = C_{3}^{I}$ and $E(t) = 0$ for all $t.$ Therefore, $\Upsilon_{3}^{\prime} = \{ Đ _{3} \}$. By applying Lyapunov-LaSalle asymptotic stability theorem we get that $Đ _{3}$ is G.A.S.

Theorem 5. If $\Re_{4} > 1$ and $\Re_{6}\leq1$, then $Đ _{4}$ is G.A.S.

Proof. Define $\Phi_{4}(S, L, I, E, Y, V, C^{I}, C^{Y})$ as:

Calculating $\frac{d\Phi_{4}}{dt}$ as:

Using the equilibrium conditions for $Đ _{4}$:

We obtain

If $\Re_{6}\leq1$, then using inequality (5.4) we get $\frac{d\Phi _{4}}{dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$, where $\frac{d\Phi_{4} }{dt} = 0$ at $S = S_{4}$, $E = E_{4},$ $Y = Y_{4}$ and $I = C^{I} = 0.$ The solutions of system (2.1) converge to $\Upsilon_{4}^{\prime}$ the largest invariant subset of $\Upsilon_{4} = \left \{ (S, L, I, E, Y, V, C^{I}, C^{Y}):\frac{d\Phi_{4}}{dt} = 0\right \}$. The set $\Upsilon _{4}^{\prime}$ contains elements with $S = S_{4}$, $E = E_{4},$ $Y = Y_{4}$, $I = 0$, and then $\dot{S} = \dot{Y} = 0$. The first and fifth equations of system (2.1) imply

which yield $V(t) = 0$ and $C^{Y}(t) = C_{4}^{Y}$ for all $t$. Since $\dot{I} = 0$, then from the third equation of system (2.1) we obtain

which yields $L(t) = 0$ for all $t.$ Therefore, $\Upsilon_{4}^{\prime} = \left \{ Đ _{4}\right \}$. By applying Lyapunov-LaSalle asymptotic stability theorem we obtain that $Đ _{4}$ is G.A.S.

Theorem 6. If $\Re_{5} > 1$, $\Re_{8}\leq1$ and $\Re_{1}/\Re_{2} > 1,$ then $Đ _{5}$ is G.A.S.

Proof. Define $\Phi_{5}(S, L, I, E, Y, V, C^{I}, C^{Y})$ as:

Calculating $\frac{d\Phi_{5}}{dt}$ as:

Using the equilibrium conditions for $Đ _{5}$:

We obtain

Then, Eq (5.8) will be reduced to the form

If $\Re_{8}\leq1$, then using inequalities (5.1)–(5.4) we get $\frac{d\Phi_{5}}{dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$, where $\frac{d\Phi_{5}}{dt} = 0$ at $S = S_{5}$, $L = L_{5},$ $I = I_{5},$ $E = E_{5},$ $Y = Y_{5},$ $V = V_{5}$ and $C^{Y} = 0.$ The solutions of system (2.1) converge to $\Upsilon_{5}^{\prime}$ the largest invariant subset of $\Upsilon_{5} = \left \{ (S, L, I, E, Y, V, C^{I}, C^{Y}):\frac{d\Phi_{5}} {dt} = 0\right \}$. The set $\Upsilon_{5}^{\prime}$ contains elements with $S = S_{5}$, $L = L_{5}$, $I = I_{5},$ $V = V_{5}$, and then $\dot {I} = 0$. Third equation of system (2.1) implies

which gives $C^{I}(t) = C_{5}^{I}$ for all $t.$ Therefore, $\Upsilon_{5} ^{\prime} = \{ Đ _{5} \}$. Applying Lyapunov-LaSalle asymptotic stability theorem we get $Đ _{5}$ is G.A.S.

Theorem 7. If $\Re_{6} > 1$, $\Re_{7}\leq1$ and $\Re_{2}/\Re_{1} > 1$, then $Đ _{6}$ is G.A.S.

Proof. Define $\Phi_{6}(S, L, I, E, Y, V, C^{I}, C^{Y})$ as:

Calculating $\frac{d\Phi_{6}}{dt}$ as:

Using the equilibrium conditions for $Đ _{6}$:

We obtain

Then, Eq (5.9) will be reduced to the form

Therefore, if $\Re_{7}\leq1$, then using inequalities (5.1)–(5.4) we get $\frac{d\Phi_{6}}{dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0,$ where $\frac{d\Phi_{6}}{dt} = 0$ at $S = S_{6}$, $L = L_{6},$ $I = I_{6},$ $E = E_{6},$ $Y = Y_{6},$ $V = V_{6}$ and $C^{I} = 0.$ Define $\Upsilon _{6} = \left \{ (S, L, I, E, Y, V, C^{I}, C^{Y}):\frac{d\Phi_{6}}{dt} = 0\right \}$ and let $\Upsilon_{6}^{\prime}$ be the largest invariant subset of $\Upsilon_{6}$. The solutions of system (2.1) tend to $\Upsilon_{6}^{\prime}$ which includes elements with $E = E_{6},$ $Y = Y_{6}$, and then $\dot{Y} = 0$. The fifth equation of system (2.1) implies

which ensures that $C^{Y}(t) = C_{6}^{Y}$ for all $t.$ Therefore, $\Upsilon _{6}^{\prime} = \left \{ Đ _{6}\right \}$. Applying Lyapunov-LaSalle asymptotic stability theorem we get $Đ _{6}$ is G.A.S.

Theorem 8. If $\Re_{7} > 1$ and $\Re_{8} > 1$, then $Đ _{7}$ is G.A.S.

Proof. Define $\Phi_{7}(S, L, I, E, Y, V, C^{I}, C^{Y})$ as:

Calculating $\frac{d\Phi_{7}}{dt}$ and after collecting terms we get

Using the equilibrium conditions for $Đ _{7}$:

We obtain

Therefore, using inequalities (5.1)–(5.4) we get $\frac {d\Phi_{7}}{dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$. In addition we have $\frac{d\Phi_{7}}{dt} = 0$ at $S = S_{7}$, $L = L_{7},$ $I = I_{7},$ $E = E_{7},$ $Y = Y_{7}$ and $V = V_{7}.$ The solutions of system (2.1) converge to $\Upsilon_{7}^{\prime}$ the largest invariant subset of $\Upsilon_{7} = \left \{ (S, L, I, E, Y, V, C^{I}, C^{Y}):\frac{d\Phi_{7}} {dt} = 0\right \}$. The set $\Upsilon_{7}^{\prime}$ contains elements with $S = S_{7}$, $L = L_{7},$ $I = I_{7},$ $E = E_{7},$ $Y = Y_{7}$ and $V = V_{7}$. Then $\dot{I} = \dot{Y} = 0$ and from the third and fifth equations of system (2.1) we get

which ensure that $C^{I}(t) = C_{7}^{I}$ and $C^{Y}(t) = C_{7}^{Y}$ for all $t.$ Therefore, $\Upsilon_{7}^{\prime} = \left \{ Đ _{7}\right \}$. Applying Lyapunov-LaSalle asymptotic stability theorem we get $Đ _{7}$ is G.A.S.

In Table 3, we summarize the global stability results given in Theorems 1–8.

6.   Numerical results and discussions

In this section, we illustrate the results of Theorems 1–8 by performing numerical simulations. Moreover, we study the effect of HTLV-I infection on the HIV mono-infected individuals by making a comparison between the dynamics of HIV mono-infection and HIV/HTLV-I co-infection. Otherwise, we investigate the influence of HIV infection on the HTLV-I mono-infected individuals by conducting a comparison between the dynamics of HTLV-I mono-infection and HIV/HTLV-I co-infection.

To solve system (2.1) numerically we fix the values of some parameters (see Table 4) and the others will be varied.

6.1.   Stability of the equilibria

In this subsection, we choose the following three different initial conditions for system (2.1):

Initial-1 :$(S(0), L(0), I(0), E(0), Y(0), V(0), C^{I}(0), C^{Y} (0)) = (600, 1.5, 1.5, 30, 0.3, 5, 1, 3)$,

Initial-2:$(S(0), L(0), I(0), E(0), Y(0), V(0), C^{I}(0), C^{Y} (0)) = (500, 1, 1, 20, 0.2, 2, 2, 2),$

Initial-3:$(S(0), L(0), I(0), E(0), Y(0), V(0), C^{I}(0), C^{Y} (0)) = (300, 0.5, 0.5, 10, 0.1, 1.5, 3, 1)$.

Choosing selected values of $\eta_{1}$, $\eta_{2}$, $\eta_{3}$, $\sigma_{1}$ and $\sigma_{2}$ under the above initial conditions leads to the following scenarios:

Scenario 1 (Stability of ${\mathit{Đ}} _{0}$): $\eta _{1} = \eta_{2} = 0.0001,$ $\eta_{3} = 0.001$ and $\sigma_{1} = \sigma _{2} = 0.2$. For this set of parameters, we have $\Re_{1} = 0.68 < 1$ and $\Re _{2} = 0.23 < 1$. Figure 2 displays that the trajectories initiating with Initial-1, Initial-2 and Initial-3 reach the equilibrium $Đ _{0} = (1000, 0, 0, 0, 0, 0, 0, 0)$. This shows that $Đ _{0}$ is G.A.S according to Theorem 1. In this situation both HIV and HTLV will be died out.

Scenario 2 (Stability of ${\mathit{Đ}} _{1}$): $\eta _{1} = 0.0005,$ $\eta_{2} = 0.0003,$ $\eta_{3} = 0.0005,$ $\sigma _{1} = 0.003$ and $\sigma_{2} = 0.2$. With such choice we get $\Re_{2} = 0.12 < 1 < 3.02 = \Re_{1}$, $\Re_{3} = 0.49 < 1$ and hence $\Re_{2}/\Re_{1} = 0.04 < 1$. Therefore, the conditions in Table 2 is verified. In fact, the equilibrium point $Đ _{1}$ exists with $Đ _{1} = \left(331.63, 9.11, 13, 0, 0, 32.51, 0, 0\right)$. Figure 3 displays that the trajectories initiating with Initial-1, Initial-2 and Initial-3 tend to $Đ _{1}$. Therefore, the numerical results support Theorem 2. This case corresponds to a chronic HIV mono-infection but with unstimulated CTL-mediated immune response.

Scenario 3 (Stability of ${\mathit{Đ}} _{2}$): $\eta _{1} = 0.0001,$ $\eta_{2} = 0.0002,$ $\eta_{3} = 0.01$, $\sigma_{1} = 0.001$ and $\sigma_{2} = 0.05$. Then, we calculate $\Re_{1} = 0.88 < 1 < 2.31 = \Re_{2}$, $\Re_{4} = 0.77 < 1$ and then $\Re_{1}/\Re_{2} = 0.38 < 1$. Hence, the conditions in Table 2 is satisfied. The numerical results show that $Đ _{2} = \left(433.33, 0, 0, 87.18, 1.31, 0, 0, 0\right)$ exists. Figure 4 illustrates that the trajectories initiating with Initial-1, Initial-2 and Initial-3 tend to $Đ _{2}$. Thus, the numerical results consistent with Theorem 3. This situation leads to a persistent HTLV mono-infection with unstimulated CTL-mediated immune response.

Scenario 4 (Stability of ${\mathit{Đ}} _{3}$): $\eta _{1} = 0.001,$ $\eta_{2} = 0.0001,$ $\eta_{3} = 0.005$ and $\sigma _{1} = \sigma_{2} = 0.01$. Then, we calculate $\Re_{3} = 1.41 > 1$ and $\Re _{5} = 0.32 < 1$. Table 2 and Figure 5 show that the trajectories initiating with Initial-1, Initial-2 and Initial-3 tend to $Đ _{3} = \left(277.78, 9.85, 10, 0, 0, 25, 1.01, 0\right)$. Therefore, $Đ _{3}$ is G.A.S and this agrees with Theorem 4. Hence, a chronic HIV mono-infection with HIV-specific CTL-mediated immune response is attained.

Scenario 5 (Stability of ${\mathit{Đ}} _{4}$): $\eta _{1} = 0.0007,$ $\eta_{2} = 0.0001,$ $\eta_{3} = 0.1$, $\sigma_{1} = 0.05$ and $\sigma_{2} = 0.3$. Then, we calculate $\Re_{4} = 5.33 > 1$ and $\Re_{6} = 0.83 < 1$. According to Table 2, $Đ _{4}$ exists with $Đ _{4} = \left(230.77, 0, 0,118.34, 0.33, 0, 0, 4.33\right)$. In Figure 6, we show that the trajectories initiating with Initial-1, Initial-2 and Initial-3 tend to $Đ _{4}$ and then it is G.A.S which agrees with Theorem 5. Hence, a chronic HTLV mono-infection with HTLV-specific CTL-mediated immune response is attained.

Scenario 6 (Stability of ${\mathit{Đ}} _{5}$): $\eta _{1} = 0.001,$ $\eta_{2} = 0.0001,$ $\eta_{3} = 0.01$, $\sigma_{1} = 0.05$ and $\sigma_{2} = 0.08$. Then, we calculate $\Re_{5} = 1.52 > 1$, $\Re_{8} = 0.83 < 1$ and $\Re_{1}/\Re_{2} = 2.19 > 1$. Table 2 and the numerical results demonstrated in Figure 7 show that $Đ _{5} = \left(433.33, 3.07, 2, 52.51, 0.79, 5, 2.98, 0\right)$ exists and it is G.A.S and this agrees with Theorem 6. As a result, a chronic co-infection with HIV and HTLV is attained where the HIV-specific CTL-mediated immune response is active and the HTLV-specific CTL-mediated immune response is unstimulated.

Scenario 7 (Stability of ${\mathit{Đ}} _{6}$): $\eta _{1} = 0.0006,$ $\eta_{2} = 0.0001,$ $\eta_{3} = 0.04$, $\sigma_{1} = 0.01$ and $\sigma_{2} = 0.5$. We compute $\Re_{6} = 1.73 > 1$, $\Re_{7} = 0.92 < 1$ and $\Re_{2}/\Re_{1} = 2.97 > 1$. Based on the conditions in Table 2, the equilibrium $Đ _{6} = \left(321.26, 5.75, 8.2, 39.54, 0.2, 20.51, 0, 1.97\right)$ exists. Moreover, the numerical results plotted in Figure 8 show that $Đ _{6}$ is G.A.S and this illustrates Theorem 7. As a result, a chronic co-infection with HIV and HTLV is attained where the HTLV-specific CTL-mediated immune response is active and the HIV-specific CTL-mediated immune response is unstimulated.

Scenario 8 (Stability of ${\mathit{Đ}} _{7}$): $\eta _{1} = 0.0006,$ $\eta_{2} = 0.0002,$ $\eta_{3} = 0.04$, $\sigma_{1} = 0.05$ and $\sigma_{2} = 0.5$. These data give $\Re_{7} = 1.55 > 1$ and $\Re_{8} = 4.31 > 1$. According to Table 2, the equilibrium $Đ _{7}$ exists. Figure 9 illustrates that the trajectories initiating with Initial-1, Initial-2 and Initial-3 tend to $Đ _{7} = \left(467.29, 2.17, 2, 57.51, 0.2, 5, 1.36, 3.31\right)$. The numerical results displayed in Figure 9 show that $Đ _{7}$ is G.A.S based on Theorem 8. In this case, a chronic co-infection with HIV and HTLV is attained where both HIV-specific CTL-mediated and HTLV-specific CTL-mediated immune responses are working.

To further confirmation, we calculate the Jacobian matrix $J = J(S, L, I, E, Y, V, C^{I}, C^{Y})$ of system (2.1) as in the following form:

Then, we calculate the eigenvalues $\lambda_{i},$ $i = 1, 2, ..., 8$ of the matrix $J$ at each equilibrium. The examined steady will be locally stable if all its eigenvalues satisfy the following condition:

We use the parameters $\eta_{1}$, $\eta_{2}$, $\eta_{3}$, $\sigma_{1}$ and $\sigma_{2}$ the same as given above to compute all positive equilibria and the corresponding eigenvalues. From the scenarios 1–8, we present in Table 5 the positive equilibria, the real parts of the eigenvalues and whether the equilibrium is locally stable or unstable.

6.2.   Comparison results

In this subsection, we study the influence of HTLV-I infection on HIV mono-infection dynamics, and how affect the HIV infection on the dynamics of HTLV-I mono-infection as well.

6.2.1.   Impact of HTLV-I infection on HIV mono-infection dynamics

To investigate the effect of HTLV-I infection on HIV mono-infection dynamics, we make a comparison between model (2.1) and the following HIV mono-infection model:

We fix parameters $\eta_{1} = 0.0006$, $\eta_{2} = 0.0001$, $\sigma_{1} = 0.05$, and $\sigma_{2} = 0.5$ and consider the following initial condition:

Initial-4: $(S(0), L(0), I(0), E(0), Y(0), V(0), C^{I}(0), C^{Y}(0)) = \left(600, 2.4, 1.8, 60, 0.2, 4.5, 1.8, 3.5\right)$.

We choose two values of the parameter $\eta_{3}$ as $\eta_{3} = 0.04$ (HIV/HTLV-I co-infection), and $\eta_{3} = 0.0$ (HIV mono-infection). It can be seen from Figure 10 that when the HIV mono-infected individual is co-infected with HTLV-I then the concentrations of susceptible CD$4^{+}$T cells, latently HIV-infected cells and HIV-specific CTLs are decreased. Although, the concentration of free HIV particles tend to the same value in both HIV mono-infection and HIV/HTLV-I co-infection. Indeed, such observation are compatible with the study that has been performed by Vandormael et al. in 2017 ^[54]. The researchers have not found any worthy differences in the concentration of HIV virus particles in comparison between HIV mono-infected and HIV/HTLV-I co-infected patients.

6.2.2.   Impact of HIV infection on HTLV-I mono-infection dynamics

To investigate the effect of HIV infection on HTLV-I mono-infection dynamics, we make a comparison between model (2.1) and the following HTLV-I mono-infection model:

We fix parameters $\eta_{3} = 0.01$; $\sigma_{1} = 0.05$, and $\sigma_{2} = 0.5$ and consider the following initial condition:

Initial-5: $(S(0), L(0), I(0), E(0), Y(0), V(0), C^{I}(0), C^{Y} (0)) = (700, 4, 2, 21, 0.198, 5, 4.5, 0.6)$.

We choose two values of the parameters $\eta_{1},$ $\eta_{2}$ as $\eta_{1} =$ $0.001,$ $\eta_{2} = 0.0002$ (HIV/HTLV-I co-infection), and $\eta_{1} = \eta _{2} = 0.0$ (HTLV-I mono-infection). It can be seen from Figure 11 that when the HTLV-I mono-infected individual is co-infected with HIV then the concentrations of susceptible CD$4^{+}$T cells, latently HTLV-infected cells and HTLV-specific CTLs are decreased. Although, the concentration of Tax-expressing HTLV-infected cells tend to the same value in both HTLV-I mono-infection and HIV/HTLV-I co-infection.

6.2.3.   Effect of CTL immnue response

As we discussed in Section 1 that CTLs have significant important in controlling HIV and HTLV-I mono-infections by killing infected cells. Model (2.1) in the absence of CTL immune response leads to a model with competition between HIV and HTLV-I on CD$4^{+}$T cells:

The system has only three equilibria, infection-free equilibrium, $\overline{Đ }_{0} = (S_{0}, 0, 0, 0, 0, 0)$, chronic HIV mono-infection equilibrium, $\overline{Đ}_{1} = (S_{1}, L_{1}, I_{1}, 0, 0, V_{1})$ and chronic HTLV mono-infection equilibrium, $\overline{Đ}_{2} = (S_{2}, 0, 0, E_{2}, Y_{2}, 0)$, where $S_{0}$, $S_{1}$, $L_{1}$, $I_{1}$, $V_{1}$, $S_{2}$, $E_{2}$ and $Y_{2}$ are given in Section 4. The existence of the these three equilibria is determined by two threshold parameters $\Re_{1}$ and $\Re_{2}$ which are defined in Section 4.

Corollary 1. For system (6.3), the following statements hold true.

(ⅰ) If $\Re_{1}\leq1$ and $\Re_{2}\leq1$, then $\overline{Đ}_{0}$ is G.A.S.

(ⅱ) If $\Re_{1} > 1$ and $\Re_{2}/\Re_{1}\leq1$, then $\overline{Đ}_{1}$ is G.A.S.

(ⅲ) If $\Re_{2} > 1$ and $\Re_{1}/\Re_{2}\leq1$, then $\overline{Đ}_{2}$ is G.A.S.

Therefore, the system will tend to one of the three equilibria $\overline{Đ}_{0}$, $\overline{Đ}_{1}$ and $\overline{Đ}_{2}$. The above result says that in the absence of immune response, the competition between HIV and HTLV-I consuming common resources, only one type of viruses with maximum basic reproductive ratio can survive. However, in our proposed model (2.1) involving HIV- and HTLV-specific CTLs, HIV and HTLV-I coexist at equilibrium. We can consider this situation as follows. Since CTL immune responses suppress viral progression, the competition between HIV and HTLV-I is also suppressed and the coexistence of HIV and HTLV-I is occurred ^[55].

7.   Conclusion

This research work formulates a mathematical model which describes the within host dynamics of HIV/HTLV-I co-infection. The model incorporated the effect of HIV-specific CTLs and HTLV-specific CTLs. HIV has two predominant infection modes: the classical VTC infection and CTC spread. The HTLV-I has two ways of transmission, (ⅰ) horizontal transmission via direct CTC contact, and (ⅱ) vertical transmission through mitotic division of Tax-expressing HTLV-infected cells. We first proved that the model is well-posed by showing that the solutions are nonnegative and bounded. We derived eight threshold parameters that governed the existence and stability of the eight equilibria of the model. We constructed appropriate Lyapunov functions and applied Lyapunov-LaSalle asymptotic stability theorem to prove the global asymptotic stability of all equilibria. We conducted numerical simulations to support and clarify our theoretical results. We studied the effect of HIV infection on HTLV-I mono-infection dynamics and vice versa. The model analysis suggested that co-infected individuals with both viruses will have smaller number of healthy CD$4^{+}$T cells in comparison with HIV or HTLV-I mono-infected individuals. We discussed the influence of CTL immune response on the co-infection dynamics.

Our model can be extended in many directions:

● In our model (2.1), we assumed that susceptible CD$4^{+}$T cells are produced at a constant rate $\rho$ and have a linear death rate $\alpha S$. It would be more reasonable to consider the density dependent production rate. One possibility is to assume a logistic growth for the susceptible CD$4^{+}$T cells in the absence of infection. Moreover, the model assumed bilinear incidence rate of infections, $\eta_{1}SV$, $\eta_{2}SI$ and $\eta_{3}SY$. However, such bilinear form may not describe the virus dynamics during the full course of infection. Therefore, it is reasonable to consider other forms of the incidence rate such as: saturated incidence, Beddington-DeAngelis incidence and general incidence ^[56,57,58].

● Model (2.1) assumed that once susceptible CD$4^{+}$T cell is contacted by an HIV or HIV-infected or HTLV-infected cell it becomes latently or actively infected instantaneously. However, such process needs time. Intracellular time delay has a significant effect on the virus dynamics. Delayed viral infection models have been constructed and analyzed in several works (see, e.g., ^{[59,60,61,62,63,64]}).

● Model (2.1) assumes that cells and viruses are equally distributed in the domain with no spatial variations. Taking into account spatial variations in the case of HIV/HTLV-I co-infection will be significant ^[65,66].

We leave these extensions as a future project.

Acknowledgment

This project was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, Saudi Arabia under grant no. (KEP-PhD-20-130-41). The authors, therefore, acknowledge with thanks DSR technical and financial support.

Conflict of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.