Analysis of an HTLV/HIV dual infection model with diffusion

1.   Introduction

Since past decades humanity has been under attack by many viruses such as hepatitis C virus (HCV), human immunodeficiency virus (HIV), hepatitis B virus (HBV), human T-lymphotropic virus type Ⅰ (HTLV-I), dengue virus and lastly coronavirus. Both HTLV-I and HIV have similar ways of transmission from infected individual to uninfected one. HTLV-I and HIV primarily target the same host, CD$4^{+}$T cells. HTLV-I can lead to adult T-cell leukemia (ATL) and HTLV-I-associated myelopathy/tropical spastic paraparesis (HAM/TSP), while HIV causes acquired immunodeficiency syndrome (AIDS). Viral infection models have become an indispensable tool to biological researchers, where they can improve the understanding of a within-host virus dynamics and help in predicting the effect of antiviral drug efficacy on disease's progression ^[1]. In 1996, Nowak and Bangham ^[2] have presented an important HIV dynamics model. After that, this model has been extended in many works (see e.g., ^{[3,4,5,6,7,8,9,10,11,13,14,15,16,17]}). All of the above mentioned models are given by ordinary or delay differential equations under the assumption that the cells and viruses are well mixed. Wang and Wang ^[18] have extended the model presented in ^[2] by incorporating spatial dependence as:

where $S(x, t)$, $I(x, t)$ and $V(x, t)$ are the concentrations of uninfected cells, active infected cells and free virus particles at position $x = \left(x_{1}, x_{2}, ..., x_{m}\right)$ and time $t$. The parameter $\rho$ represents the creation rate of the uninfected cells. The free virus particles infect the uninfected cells via free-to-cell (FTC) transmission at rate $\varkappa_{1} SV$. The infected cells produce viruses at rate $bI$. The uninfected cells, active infected cells and free virus particles are die with rates $\alpha S$, $aI$ and $\varepsilon V$, respectively. Here, $d_{V}$ is the diffusion coefficient and $\Delta$ is the Laplacian operator. In ^[19], Kang et al. have studied a four-dimensional diffusive viral infection model with Crowley-Martin infection rate. Model (1.1) has been extended by including different factors; (i) time delay ^[19,20], (ii) different forms of the incidence rate ^[19,20], (iii) Cytotoxic T lymphocytes (CTLs) immune response ^[19], and (iv) both CTL and humoral immune responses ^[21].

In model (1.1), it has been assumed that the virus can only infect the target cell via FTC contact. In case of HIV, the infected cell can infect the target cell via direct infected-to-cell (ITC) contact ^[22]. Wang et al. ^[23] have extended model (1.1) by incorporating ITC transmission as:

where the term $\varkappa_{2}SI$ represents the ITC incidence rate. This model has been generalized in ^[24] by including time delay and general FTC incidence rate function in the form $f(S, V)$. Model (1.2) assumes that the virus particles can move based on Fickian diffusion, while the cells can not. In ^{[25,26,27,28,29,30,31]}, it has been assumed that the viruses, uninfected cells, infected cells and immune cells can diffuse.

Modeling and analysis of HTLV-I single infection have been addressed in several works ^{[32,33,34,35,36]}. The effect of CTLs on HTLV-I dynamics has been addressed in many works (see e.g. ^{[37,38,39,40,41,42,43]}). Lim and Maini ^[37] have proposed an HTLV-I dynamics model with mitotic division of active HTLV-infected cells and CTL immunity. Both mitotic division of active HTLV-infected cells and CTL immune response have been included in the HTLV-I dynamics in many papers (see e.g., ^[37,44,45]). All of these HTLV dynamics models did not include the diffusion of the viruses and cells. Wang and Ma ^[46] have introduced a diffusive HTLV-I infection model with mitotic division of active infected cells and CTL immune response:

where $E(x, t)$, $Y(x, t)$ and $C^{Y}(x, t)$ are the concentrations of latent HTLV-infected cells, active HTLV-infected cells and HTLV-specific CTLs at position $x$ and time $t$, respectively. The uninfected CD$4^{+}$T cells become HTLV-infected cells due to ITC contact at rate $\varkappa_{3}SY$ (horizontal transmission). The fraction $\varphi \in(0, 1)$ represents the probability of new HTLV infections via horizontal transmission could be enter a latent period. The term $r^{\ast}Y$ (vertical transmission) represents the rate at which active HTLV-infected cells become latent. The terms $\omega E$ and $\delta^{\ast}Y$ denote the death rates of the latent and active HTLV-infected cells, respectively. The latent HTLV-infected cells are activated with rate $\psi E$. The active HTLV-infected cells are killed by their specific CTLs at rate $\mu_{2}C^{Y}Y$. The linear term $\sigma_{2}Y$ represents the expansion rate of HTLV-specific CTLs. The HTLV-specific CTLs decay at rate $\pi_{2}C^{Y}$.

During the last decades HTLV-I and HIV dual infection has been extensively reported. It has been discovered that the simultaneous infection by the two viruses affects the pathogenic development and influences the outcomes for associated chronic diseases ^[47]. In fact, concurrent infections with HTLV-I and HIV have occurred frequently in areas where peoples living at high risk activities such as needle injection sharing and unprotected sexual relationships. In addition, HTLV/HIV dual infection have documented in specific geographic regions where both retroviruses become endemic ^[48], and among those who belonged to a specific ethnic as well. For instance, the dual infection rates in peoples living in some parts of Brazil have reached 16% of HIV-infected patients ^[49]. In a recent work, it has been estimated that the HIV single infected patients are more exposure to be dually infected with HTLV-I at a higher rate initiating from 100 to 500 times in comparison to the uninfected peoples ^[50]. Moreover, some seroepidemiologic studies have reported that HTLV-infected patients are at risk to have a concurrent infection with HIV, and vice versa compared to those who are infection-free from the general population ^[48]. HTLV-I and HIV are mainly attack the CD$4^{+}$T cells and lead to immune dysfunctional as well, however, they also conflict no doubt with respect to the etiology of their pathogenic and clinical outcomes ^[51]. HTLV-I and HIV dual infection appears to have an overlap on the course of associated clinical outcomes with both viruses ^[48]. Many researchers have reported that HIV infected individuals who are possibly infected with HTLV-I simultaneously can potentially associated with clinical progression to AIDS. In contrast, HIV can modify HTLV-I expression in dual infected patients which leads them to a higher risk of developing HTLV-I related diseases such as TSP/HAM and ATL ^[48,50].

While many efforts have been made to investigate mathematical modeling and analysis of both HTLV-I and HIV single infection, almost none have focused on the modeling of HTLV/HIV dual infection dynamics. The only exceptions are the very recent works presented in ^[52,53], however, in these papers the diffusion of the viruses and cells has been neglected. Therefore, the contributions of the present paper can be stated as follows: (i) formulate an HTLV/HIV dual infection model taking into account the mobility of cells and viruses, (ii) study the basic properties of the proposed model, (iii) calculate all possible steady states and derive their existence conditions, (iv) study the global stability of all steady states using Lyapunov-LaSalle asymptotic stability theorem, (v) perform some numerical simulations to illustrate the theoretical results.

Since an individual can be infected with two or more viruses in the same time, our model may be helpful to study different dual infections such as Coronavirus/Influenza, HIV/HCV, HIV/HBV and HIV/Malaria.

2.   Model formulation

We set up a partial differential equation model that describes the change of concentrations of eight compartments with respect to position $x$ and time $t$; uninfected CD$4^{+}$T cells $S(x, t)$, latent HIV-infected cells $L(x, t)$, active HIV-infected cells $I(x, t)$, latent HTLV-infected cells $E(x, t)$, active HTLV-infected cells $Y(x, t)$, free HIV particles $V(x, t)$, HIV-specific CTLs $C^{I}(x, t),$ and HTLV-specific CTLs $C^{Y}(x, t)$. We consider the following factors:

(F1) The uninfected CD$4^{+}$T cells are the main target of each of HTLV-I and HIV;

(F2) There exist latent HIV-infected and HTLV-infected cells;

(F3) Bilinear specific CTL immune response for both HTLV-I and HIV;

(F4) The HIV can spread when an uninfected CD$4^{+}$T cell is contacted with free HIV particle (FTC infection) or active HIV-infected cell (ITC infection);

(F5) HTLV-I can be transmitted via two routes, (i) horizontal transmission via direct ITC touch via virological synapse, and (ii) vertical transmission by mitotic division of active HTLV-infected cells.

(F6) Spatial diffusion for all compartments.

Taking into account factors (F1)-(F6) we propose the following PDEs model:

where $x\in \Gamma,$ $t > 0$. A fraction $\beta \in(0, 1)$ of new HIV-infected cells will be active, and the remaining part $1-\beta$ will be latent. Latent HIV-infected cells are transmitted to be active at rate $\lambda L$ and die at rate $\gamma L$. The term $\mu_{1}C^{I}I$ is the killing rate of active HIV-infected cells due to their specific immunity. The term $\kappa r^{\ast }Y,$ $\kappa \in \left(0, 1\right)$ represents the rate active HTLV-infected cells that transmit to latent HTLV-infected cells and get-away from the immune system ^[45]. The expansion rate for HIV-specific CTLs and HTLV-specific CTLs are represented by $\sigma_{1}C^{I}I$ and $\sigma_{2} C^{Y}Y$, respectively. The HIV-specific CTLs decay at rate $\pi_{1}C^{I}$. All remaining parameters have the same biological meaning as explained in Section 1. The spatial domain $\Gamma \subset \mathbb{R}^{m},$ $m\geq1$ is connected and bounded with a smooth boundary $\partial \Gamma,$ while $d_{\mathcal{U}}$ is the diffusion coefficient where $\mathcal{U}\in \{S, L, I, E, Y, V, C^{I}, C^{Y}\}$. The initial conditions are given by

where $\mathcal{G}_{i}(x)$, $i = 1, ..., 8,$ are non-negative and continuous functions. In addition, we take the following homogeneous Neumann boundary conditions:

where $\dfrac{\partial}{\partial \mathcal{\vec{V}}}$ is the outward normal derivative on the boundary $\partial \Gamma$. These boundary conditions indicate that cells and viruses cannot cross the isolated boundary ^[54].

We assume that $r^{\ast} < \min \left \{ \alpha, \omega, \delta^{\ast}\right \}$ ^[37]. It follows that $\left(1-\kappa \right) r^{\ast} < \delta^{\ast}$ and then

Let $\delta = \delta^{\ast}-\left(1-\kappa \right) r^{\ast}$ and $r = \kappa r^{\ast}$. Therefore, model (2.1) can be written as:

3.   Well-posedness of solutions

Proposition 1. Assume that $d_{S} = d_{L} = d_{I} = d_{E} = d_{Y} = d_{V} = d_{C^{I}} = d_{C^{Y}} = \tilde{d}$. Then, model (2.4) with any initial satisfying (2.2) has a unique, non-negative and bounded solution defined on $\bar{\Gamma}\times \lbrack0, +\infty)$.

Proof. We denote $\mathcal{X} = BUC\left(\bar{\Gamma}, \mathbb{R} ^{8}\right)$ the set of all bounded and uniformly continuous functions from $\bar{\Gamma}$ to $\mathbb{R}^{8}$, with norm $\Vert \theta \Vert_{\mathcal{X} } = \sup \limits_{x\in \bar{\Gamma}}|\theta(x)|$. Define the positive cone $\mathcal{X}_{+} = BUC\left(\bar{\Gamma}, \mathbb{R}_{+}^{8}\right) \mathbb{\subset}\mathcal{X}$ which induces a partial order on $\mathcal{X}$. This shows that the space $\left(\mathcal{X}, \Vert \cdot \Vert_{\mathcal{X} }\right)$ is a Banach lattice ^[55,56].

For any initial data $\mathcal{G} = (\mathcal{G}_{1}, \mathcal{G}_{2}, \mathcal{G}_{3}, \mathcal{G}_{4}, \mathcal{G}_{5}, \mathcal{G}_{6}, \mathcal{G}_{7}, \mathcal{G}_{8})^{T}\in \mathcal{X}_{+}$, we define $H = (H_{1}, H_{2}, H_{3}, H_{4}, H_{5}, H_{6}, H_{7}, H_{8})^{T}:\mathcal{X} _{+}\rightarrow \mathcal{X}$ by

We note that $H$ is locally Lipschitz on $\mathcal{X}_{+}$. System (2.4) with initial conditions (2.2) and boundary conditions (2.3) can be rewritten as the following abstract functional differential equation

where $\mathcal{\bar{U}} = (S, L, I, E, Y, V, C^{I}, C^{Y})^{T}$ and $\Theta \mathcal{\bar{U}} = \left(d_{S}\Delta S, d_{L}\Delta L, d_{I}\Delta I, d_{E}\Delta E, d_{Y}\Delta Y, d_{V}\Delta V, d_{C^{I}}\Delta C^{I}, \right. \left. d_{C^{Y}}\Delta C^{Y}\right) ^{T}.$ It is possible to show that

It follows from ^[55,56,57] that for any $\mathcal{G} \in \mathcal{X}_{+}$, system (2.4) with (2.2)-(2.3) has a unique non-negative mild solution $(S\left(x, t\right), L\left(x, t\right), I\left(x, t\right), E\left(x, t\right), Y\left(x, t\right), V\left(x, t\right), C^{I}\left(x, t\right), C^{Y}\left(x, t\right))$ defined on $\bar{\Gamma}\times \lbrack0, \mathcal{T}_{m})$, where $[0, \mathcal{T}_{m})$ is the maximal existence time interval on which the solution exists. In addition, this solution also is a classical solution for the given problem.

To prove the boundedness of solutions, we define

Since $d_{S} = d_{L} = d_{I} = d_{E} = d_{Y} = d_{V} = d_{C^{I}} = d_{C^{Y}} = \tilde{d}$, then using system (2.4) we obtain

We have $\delta-r = \delta^{\ast}-r^{\ast} > 0$. Hence,

where $\phi = \min \{ \alpha, \gamma, \frac{a}{2}, \omega, \delta^{\ast}-r^{\ast }, \varepsilon, \pi_{1}, \pi_{2}\}$. Thus, $\Psi(x, t)$ satisfies the following system

Let $\widetilde{\Psi}(t)$ be a solution of the following ODE

This gives that $\widetilde{\Psi}(t)\leq \max \left \{ \frac{\rho}{\phi }, \underset{x\in \bar{\Gamma}}{\max}\Psi(x, 0)\right \}$. On the basis of comparison principle (see ^[58]), we obtain $\Psi(x, t)\leq \widetilde{\Psi}(t)$. Then, we get

which implies that $S\left(x, t\right),$ $L\left(x, t\right),$ $I\left(x, t\right),$ $E\left(x, t\right),$ $Y\left(x, t\right),$ $V\left(x, t\right),$ $C^{I}\left(x, t\right),$ and $C^{Y}\left(x, t\right)$ are bounded on $\bar{\Gamma}\times \lbrack0, \mathcal{T}_{m})$. We deduce from the standard theory for semi-linear parabolic systems that $\mathcal{T} _{m} = +\infty$ ^[59]. This shows that solution $\left(S\left(x, t\right), L\left(x, t\right), I\left(x, t\right), E\left(x, t\right), Y\left(x, t\right), V\left(x, t\right), C^{I}\left(x, t\right), C^{Y}\left(x, t\right) \right)$ is defined for all $x\in \Gamma,$ $t > 0$ and also is unique and non-negative.

4.   Steady state analysis

In this section, we calculate the steady states and derive the threshold parameters which guarantee their existence. The steady states of system (2.4) satisfying the following equations:

We find that system (2.4) has eight possible steady states.

(i) Infection-free steady state, ${\text{Ð}}_{0} = (S_{0}, 0, 0, 0, 0, 0, 0, 0)$, where $S_{0} = \rho/\alpha$. In this case the body is free from HIV and HTLV.

(ii) Persistent HIV single infection steady state with an ineffective immune response, ${\text{Ð}}_{1} = (S_{1}, L_{1}, I_{1}, 0, 0, V_{1}, 0, 0),$ where

The parameter $\Re_{1}$ represents the basic HIV single infection reproduction number for system (2.4) and is defined as:

The parameter $\Re_{1}$ decides whether or not a persistent HIV single infection can be established. It is clear that at the steady state ${\text{Ð}}_{1}$ the HIV single infection persists with ineffective immune response.

(iii) Persistent HTLV single infection steady state with an ineffective immune response, ${\text{Ð}}_{2} = (S_{2}, 0, 0, E_{2}, Y_{2}, 0, 0, 0)$, where

The parameter $\Re_{2}$ denotes the basic HTLV single infection reproduction number for system (2.4) and is defined as:

The parameter $\Re_{2}$ decides whether or not a persistent HTLV single infection can be established. At the steady state ${\text{Ð}}_{2}$ the HTLV single infection persists with an ineffective immune response.

(iv) Persistent HIV single infection steady state with only effective HIV-specific CTL, ${\text{Ð}}_{3} = (S_{3}, L_{3}, I_{3}, 0, 0, V_{3}, C_{3}^{I}, 0)$, where

and

is the HIV-specific CTL immunity reproduction number in case of HIV single infection. The parameter $\Re_{3}$ determines whether or not the HIV-specific CTL immune response is effective in the absence of HTLV. At the steady state ${\text{Ð}}_{3}$ the HIV single infection persists with an effective immune response.

(v) Persistent HTLV single infection steady state with only effective HTLV-specific CTL, ${\text{Ð}}_{4} = (S_{4}, 0, 0, E_{4}, Y_{4}, 0, 0, C_{4}^{Y})$, where

The HTLV-specific CTL immunity reproduction number in case of HTLV single infection is stated as:

The parameter $\Re_{4}$ determines whether or not the HTLV-specific CTL immune response is effective in the absence of HIV. At the steady state ${\text{Ð}}_{4}$ the HTLV single infection persists with an effective immune response.

(vi) Persistent HTLV/HIV dual infection steady state with only effective HIV-specific CTL, ${\text{Ð}}_{5} = (S_{5}, L_{5}, I_{5}, E_{5}, Y_{5}, V_{5}, C_{5}^{I}, 0)$, where

The HTLV infection reproduction number in the presence of HIV infection is stated as:

It is obvious that the parameter $\Re_{5}$ determines whether or not HIV-infected patients could be dually infected with HTLV. At the steady state ${\text{Ð}}_{5}$ the HTLV/HIV dual infection persists with only HIV-specific CTL immune response.

(vii) Persistent HTLV/HIV dual infection steady state with only effective HTLV-specific CTL, ${\text{Ð}}_{6} = (S_{6}, L_{6}, I_{6}, E_{6}, Y_{6}, V_{6}, 0, C_{6} ^{Y})$, where

The HIV infection reproduction number in the presence of HTLV infection is stated as:

It is clear that the parameter $\Re_{6}$ determines whether or not HTLV-infected patients could be dually infected with HIV. At the steady state ${\text{Ð}}_{6}$ the HTLV/HIV dual infection persists with only HTLV-specific CTL immune response.

(viii) Persistent HTLV/HIV dual infection steady state with effective HIV-specific CTL and HTLV-specific CTL, ${\text{Ð}}_{7} = (S_{7}, L_{7}, I_{7}, E_{7}, Y_{7}, V_{7}, C_{7}^{I}, C_{7}^{Y})$, where

and

The parameter $\Re_{7}$ is the competed HIV-specific CTL immunity reproduction number in case of HTLV/HIV dual infection. The parameter $\Re _{8}$ is the competed HTLV-specific CTL immunity reproduction number in case of HTLV/HIV dual infection. Clearly, ${\text{Ð}}_{7}$ exists when $\Re_{7} > 1$ and $\Re_{8} > 1$.

5.   Global stability analysis

In this section, we analyze the global asymptotic stability of all steady states by Lyapunov method. For constructing Lyapunov functions we follow the works in ^[60,61]. To prove Theorems 1–8 we need the arithmetic-geometric mean inequality

which yields

Consider a function $\Phi_{j}(S, L, I, E, Y, V, C^{I}, C^{Y})$ and define

Let $\Upsilon_{j}^{^{\prime}}$ be the largest invariant subset of

We define a function

Numann boundary conditions (2.3) and Divergence Theorem imply that

for $\mathcal{U}\in \{S, L, I, E, Y, V, C^{I}, C^{Y}\}.$ Thus, we obtain

For convenience, we drop the input notation i.e., $(S, L, I, E, Y, V, C^{I}, C^{Y}) = (S(x, t), L(x, t), I(x, t), E(x, t), Y(x, t), V(x, t), C^{I}(x, t), C^{Y}(x, t))$.

Theorem 1. If $\Re_{1}\leq1$ and $\Re_{2}\leq1$, then ${\text{Ð}}_{0}$ is globally asymptotically stable (GAS).

Proof. Define $\Phi_{0}(x, t)$ as:

Clearly, $\hat{\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 $\hat{\Phi}_{0}(S_{0}, 0, 0, 0, 0, 0, 0, 0) = 0$. We calculate $\frac{\partial \Phi_{0}}{\partial t}$ along the solutions of model (2.4) as:

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

Consequently, we calculate $\frac{d\hat{\Phi}_{0}}{dt}$ as follows:

Using equality (5.5), Eq (5.6) is reduced to the following form

Therefore, $\frac{d\hat{\Phi}_{0}}{dt}\leq0$ for all $S, I, Y, C^{I}, C^{Y} > 0$ and $\frac{d\hat{\Phi}_{0}}{dt} = 0$ with equality holding when $(S, I, Y, C^{I}, C^{Y}) = (S_{0}, 0, 0, 0, 0)$. The solutions of system (2.4) converge to $\Upsilon_{0}^{^{\prime}}$. The elements of $\Upsilon_{0}^{^{\prime}}$ satisfy $(S, I, Y, C^{I}, C^{Y}) = (S_{0}, 0, 0, 0, 0)$ and then $\frac{\partial S}{\partial t} = \frac{\partial Y}{\partial t} = \Delta S = \Delta Y = 0$. The first and fifth equations of system (2.4) reduce to

This yields $V = E = 0$. Further, we have $\frac{\partial I}{\partial t} = \Delta I = 0$, then the third equation of system (2.4) becomes

which provides that $L = 0$. Hence, $\Upsilon_{0}^{^{\prime}} = \left \{ {\text{Ð}}_{0}\right \}$ and by applying Lyapunov-LaSalle asymptotic stability theorem ^[62,63,64] we get that ${\text{Ð}}_{0}$ is GAS.

Theorem 2. Let $\Re_{1} > 1$, $\Re_{2}/\Re_{1}\leq1$ and $\Re_{3}\leq 1$, then ${\text{Ð}}_{1}$ is GAS.

Proof. Define a function $\Phi_{1}(x, t)$ as:

Calculating $\frac{\partial \Phi_{1}}{\partial t}$ as:

The steady state conditions of ${\text{Ð}}_{1}$ imply that

Then, we obtain

Therefore, Eq (5.7) becomes

Calculating the time derivative of $\hat{\Phi}_{1}(t)$ and using equality (5.5) to get

Since $\Re_{2}/\Re_{1}\leq1$ and $\Re_{3}\leq1$, then utilizing inequalities (5.1)–(5.3) we obtain $\frac{d\hat{\Phi}_{1}}{dt}\leq0$ for all $S, L, I, Y, V, C^{I}, C^{Y} > 0$. Further, $\frac{d\hat{\Phi}_{1}}{dt} = 0$ at $(S, L, I, V, Y, C^{I}, C^{Y}) = (S_{1}, L_{1}, I_{1}, V_{1}, 0, 0, 0).$ The solutions of system (2.4) tend to $\Upsilon_{1}^{\prime}$ which contains elements with $Y = 0$, and hence $\frac{\partial Y}{\partial t} = \Delta Y = 0$. The fifth equation of system (2.4) reduces to

which yields $E = 0$. Hence, $\Upsilon_{1}^{\prime} = \left \{ {\text{Ð}} _{1} \right \}$ and ${\text{Ð}}_{1}$ is GAS by using Lyapunov-LaSalle asymptotic stability theorem.

Theorem 3. If $\Re_{2} > 1$, $\Re_{1}/\Re_{2}\leq1$ and $\Re_{4}\leq1$, then ${\text{Ð}}_{2}$ is GAS.

Proof. Let $\Phi_{2}(x, t)$ be defined as:

We calculate $\frac{\partial \Phi_{2}}{\partial t}$ as:

Using the steady state conditions for ${\text{Ð}}_{2}$:

we obtain

Therefore, we take the time derivative of $\hat{\Phi}_{2}(t)$ along the positive solutions of (2.4) and use equality (5.5) to find

Thus, if $\Re_{1}/\Re_{2}\leq1$ and $\Re_{4}\leq1$, then from inequality (5.4) we obtain $\frac{d\hat{\Phi}_{2}}{dt}\leq0$ for all $S, I, E, Y, C^{I}, C^{Y} > 0$. In addition, $\frac{d\hat{\Phi}_{2}}{dt} = 0$ at $(S, E, Y, I, C^{I}, C^{Y}) = (S_{2}, E_{2}, Y_{2}, 0, 0, 0)$. The solutions of model (2.4) tend to $\Upsilon_{2}^{^{\prime}}$ which satisfy $(S, Y, I) = (S_{2}, Y_{2}, 0)$ and hence $\frac{\partial S}{\partial t} = \Delta S = 0$. The first equation of system (2.4) becomes

From conditions (5.9) we get $V = 0$. Moreover, we have $\frac{\partial I}{\partial t} = \Delta I = 0$, then the third equation of system (2.4) becomes

which provides $L = 0$. Thus, $\Upsilon_{2}^{^{\prime}} = \left \{ {\text{Ð}} _{2} \right \}$ and by applying Lyapunov-LaSalle asymptotic stability theorem we get that ${\text{Ð}}_{2}$ is GAS.

Theorem 4. If $\Re_{3} > 1$ and $\Re_{5}\leq1$, then ${\text{Ð}}_{3}$ is GAS.

Proof. Define a function $\Phi_{3}(x, t)$ as:

We calculate $\frac{\partial \Phi_{3}}{\partial t}$ as:

Using the steady state conditions for ${\text{Ð}}_{3}$:

we obtain

After taking the derivative of $\hat{\Phi}_{3}(t)$ with respect to time $t$ and using equality (5.5), Eq (5.11) will take the form

Since $\Re_{5}\leq1$, then from inequalities (5.1)–(5.3) we obtain $\frac{d\hat{\Phi}_{3}}{dt}\leq0$ for all $S, L, I, Y, V, C^{I}, C^{Y} > 0$. In addition $\frac{d\hat{\Phi}_{3}}{dt} = 0$ when $(S, L, I, V, C^{I}, Y, C^{Y}) = (S_{3}, L_{3}, I_{3}, V_{3}, C_{3}^{I}, 0, 0)$. The trajectories of system (2.4) tend to $\Upsilon_{3}^{^{\prime}}$ which has elements satisfying $Y = 0$. Hence, $\frac{\partial Y}{\partial t} = \Delta Y = 0$ and the fifth equation of system (2.4) becomes

which yields $E = 0\ $and hence, $\Upsilon_{3}^{^{\prime}} = \left \{ {\text{Ð}} _{3} \right \}$. Applying Lyapunov-LaSalle asymptotic stability theorem we get ${\text{Ð}}_{3}$ is GAS.

Theorem 5. If $\Re_{4} > 1$ and $\Re_{6}\leq1$, then ${\text{Ð}}_{4}$ is GAS.

Proof. Define $\Phi_{4}(x, t)$ as:

Calculating $\frac{\partial \Phi_{4}}{\partial t}$ as:

Using the steady state conditions for ${\text{Ð}}_{4}$:

We obtain

Calculating $\frac{d\hat{\Phi}_{4}}{dt}$ and using equality (5.5) we to obtain

Hence, if $\Re_{6}\leq1$, then from inequality (5.4) we obtain $\frac{d\hat{\Phi}_{4}}{dt}\leq0$ for all $S, I, E, Y, C^{I}, C^{Y} > 0$. Moreover, $\frac{d\hat{\Phi}_{4}}{dt} = 0$ at $(S, E, Y, C^{Y}, I, C^{I}) = (S_{4}, E_{4}, Y_{4}, C_{4}^{Y}, 0, 0)$. The solutions of model (2.4) tend to $\Upsilon_{4}^{^{\prime}}$ which has elements satisfying $(S, Y, I) = (S_{4}, Y_{4}, 0)$, and then $\frac{\partial S}{\partial t} = \Delta S = 0$. Further, the first equation of system (2.4) becomes

which yields $V = 0$. Furthermore, we have $\frac{\partial I}{\partial t} = \Delta I = 0$ and the third equation of system (2.4) reduces to

which gives $L = 0\ $and hence, $\Upsilon_{4}^{^{\prime}} = \left \{ {\text{Ð}}_{4}\right \}$. Applying Lyapunov-LaSalle asymptotic stability theorem we get ${\text{Ð}}_{4}$ is GAS.

Theorem 6. If $\Re_{5} > 1$, $\Re_{8}\leq1$ and $\Re_{1}/\Re_{2} > 1,$ then ${\text{Ð}}_{5}$ is GAS.

Proof. Define $\Phi_{5}(x, t)$ as:

Calculating $\frac{\partial \Phi_{5}}{\partial t}$ as:

Using the steady state conditions for ${\text{Ð}}_{5}$:

We obtain

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

Calculating $\frac{d\hat{\Phi}_{5}}{dt}$ along the solution trajectories of system (2.4) and using equality (5.5) to get

Hence, if $\Re_{8}\leq1$, then from inequalities (5.1)–(5.4) we obtain $\frac{d\hat{\Phi}_{5}}{dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$. We have also $\frac{d\hat{\Phi}_{5}}{dt} = 0$ at $(S, L, I, E, Y, V, C^{I}, C^{Y}) = (S_{5}, L_{5}, I_{5}, E_{5}, Y_{5}, V_{5}, C_{5}^{I}, 0)$. The trajectories of system (2.4) converge to $\Upsilon_{5}^{^{\prime}}\ $and hence, $\Upsilon_{5}^{^{\prime}} = \left \{ {\text{Ð}} _{5} \right \}$. Applying Lyapunov-LaSalle asymptotic stability theorem we get ${\text{Ð}}_{5}$ is GAS.

Theorem 7. If $\Re_{6} > 1$, $\Re_{7}\leq1$ and $\Re_{2}/\Re_{1} > 1$, then ${\text{Ð}}_{6}$ is GAS.

Proof. Define $\Phi_{6}(x, t)$ as:

Calculating $\frac{\partial \Phi_{6}}{\partial t}$ as:

Using the steady state conditions for ${\text{Ð}}_{6}$:

We obtain

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

We calculate $\frac{d\hat{\Phi}_{6}}{dt}$ along the solution trajectories of system (2.4) and then we use equality (5.5) to obtain

Hence, if $\Re_{7}\leq1$, then using inequalities (5.1)–(5.4) we get $\frac{d\hat{\Phi}_{6}}{dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$, where $\frac{d\hat{\Phi}_{6}}{dt} = 0$ when $(S, L, I, E, Y, V, C^{Y}, C^{I}) = (S_{6}, L_{6}, I_{6}, E_{6}, Y_{6}, V_{6}, C_{6}^{Y}, 0)$. The solutions of model (2.4) tend to $\Upsilon_{6}^{^{\prime}}\ $and hence $\Upsilon _{6}^{^{\prime}} = \left \{ {\text{Ð}}_{6}\right \}$. Applying Lyapunov-LaSalle asymptotic stability theorem we get ${\text{Ð}}_{6}$ is GAS.

Theorem 8. If $\Re_{7} > 1$ and $\Re_{8} > 1$, then ${\text{Ð}}_{7}$ is GAS.

Proof. Define $\Phi_{7}(x, t)$ as:

Calculating $\frac{\partial \Phi_{7}}{\partial t}$ as:

Using the steady state conditions for ${\text{Ð}}_{7}$:

We obtain

Calculating $\frac{d\hat{\Phi}_{7}}{dt}$ and using equality (5.5) we obtain

Inequalities (5.1)–(5.4) imply that $\frac{d\hat{\Phi}_{7}} {dt}\leq0$ for all $S, L, I, E, Y, V, C^{I}, C^{Y} > 0$. Moreover, $\frac{d\hat{\Phi }_{7}}{dt} = 0$ when $(S, L, I, E, Y, V, C^{I}, C^{Y}) = (S_{7}, L_{7}, I_{7}, E_{7}, Y_{7}, V_{7}, C_{7}^{I}, C_{7}^{Y}).$ The solutions of model (2.4) converge to $\Upsilon_{7}^{^{\prime}} = \left \{ {\text{Ð}}_{7}\right \}$. Applying Lyapunov-LaSalle asymptotic stability theorem we get ${\text{Ð}}_{7}$ is GAS.

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

6.   Numerical simulations

In this section, we numerically show the global stability of steady states using the values of the parameters given in Table 2. Moreover, we present comparison between single and dual infections. We choose the spatial domain as $\Gamma = [0, 2]$ with a step size $0.02$. The step size for time is given by $0.1$. Further, we choose the following initial conditions for system (2.1):

In addition, we consider the homogeneous Neumann boundary conditions:

6.1.   Stability of the steady states

In this subsection, we select different values of $\varkappa_{1}$, $\varkappa_{2}$, $\varkappa_{3}$, $\sigma_{1},$ and $\sigma_{2}$ under the above initial and boundary conditions which leads to the following strategies:

Strategy 1 (Stability of ${\bf{ Ð}}_{0}$): $\varkappa _{1} = \varkappa_{2} = 0.0001,$ $\varkappa_{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 1 shows that the solution of system (2.1) converges the steady state ${\text{Ð}}_{0} = (1000, 0, 0, 0, 0, 0, 0, 0)$. This shows that ${\text{Ð}}_{0}$ is GAS according to Theorem 1. In this case both HTLV-I and HIV will be cleared.

Strategy 2 (Stability of ${\bf{ Ð}}_{1}$): $\varkappa_{1} = 0.0005,$ $\varkappa_{2} = 0.0003,$ $\varkappa_{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$. Theorem 2 implies that ${\text{Ð}}_{1} = \left(331.63, 9.11, 13, 0, 0, 32.51, 0, 0\right)$ is GAS. This will lead to the situation of persistent HIV single infection but with an ineffective CTL immune response.

Strategy 3 (Stability of ${\bf{ Ð}} _{2}$): $\varkappa_{1} = 0.0001,$ $\varkappa_{2} = 0.0002,$ $\varkappa_{3} = 0.01$, $\sigma_{1} = 0.001,$ and $\sigma_{2} = 0.05$. Then, we calculate $\Re _{1} = 0.88 < 1 < 2.33 = \Re_{2}$, $\Re_{4} = 0.78 < 1$ and then $\Re_{1}/\Re_{2} = 0.38 < 1$. The numerical results show that ${\text{Ð}}_{2} = \left(428, 0, 0, 88.74, 1.34, 0, 0, 0\right)$ exists and is GAS according to Theorem 3. It means that, a persistent HTLV single infection with an ineffective CTL immune response will be reached.

Strategy 4 (Stability of ${\bf{ Ð}} _{3}$): $\varkappa_{1} = 0.001,$ $\varkappa_{2} = 0.0001,$ $\varkappa_{3} = 0.005,$ and $\sigma_{1} = \sigma_{2} = 0.01$. Then, we calculate $\Re_{3} = 1.41 > 1$ and $\Re_{5} = 0.32 < 1$. The numerical results show that ${\text{Ð}}_{3} = \left(277.78, 9.85, 10, 0, 0, 25, 1.01, 0\right)$ is GAS based on Theorem 4. Hence, a persistent HIV single infection with effective HIV-specific CTL immune response is attained.

Strategy 5 (Stability of ${\bf{ Ð}} _{4}$): $\varkappa_{1} = 0.0007,$ $\varkappa_{2} = 0.0001,$ $\varkappa_{3} = 0.1$, $\sigma_{1} = 0.05,$ and $\sigma_{2} = 0.3$. Then, we calculate $\Re_{4} = 5.38 > 1$ and $\Re_{6} = 0.83 < 1$. According to these data ${\text{Ð}}_{4}$ exists with ${\text{Ð}}_{4} = \left(230.77, 0, 0,118.53, 0.33, 0, 0, 4.34\right)$ and it is GAS based on Theorem 5. In this case, a persistent HTLV single infection with effective HTLV-specific CTL immunity is reached.

Strategy 6 (Stability of ${\bf{ Ð}} _{5}$): $\varkappa_{1} = 0.001,$ $\varkappa_{2} = 0.0001,$ $\varkappa_{3} = 0.01$, $\sigma_{1} = 0.05,$ and $\sigma_{2} = 0.08$. Then, we calculate $\Re_{5} = 1.53 > 1$, $\Re_{8} = 0.84 < 1$ and $\Re_{1}/\Re_{2} = 2.17 > 1$. The numerical results show that ${\text{Ð}}_{5} = \left(428, 3.03, 2, 54.21, 0.82, 5, 2.91, 0\right)$ exists and it is GAS and this supports Theorem 6. This case leads to a persistent dual infection with HTLV and HIV where the HIV-specific CTL immunity is effective and the HTLV-specific CTL immunity is ineffective.

Strategy7 (Stability of ${\bf{ Ð}} _{6}$): $\varkappa_{1} = 0.0006,$ $\varkappa_{2} = 0.0001,$ $\varkappa_{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.99 > 1$. The numerical outcomes show that ${\text{Ð}}_{6} = \left(321.26, 5.75, 8.2, 39.65, 0.2, 20.51, 0, 1.98\right)$ is GAS which support Theorem 7. This situation leads to a persistent dual infection with HTLV and HIV where the HTLV-specific CTL immunity is effective and the HIV-specific CTL immunity is not working.

Strategy 8 (Stability of ${\bf{ Ð}} _{7}$): $\varkappa_{1} = 0.0006,$ $\varkappa_{2} = 0.0002,$ $\varkappa_{3} = 0.04$, $\sigma_{1} = 0.05,$ and $\sigma_{2} = 0.5$. These data give $\Re_{7} = 1.55 > 1$ and $\Re_{8} = 4.35 > 1$. Figure 2 illustrates that ${\text{Ð}}_{7} = \left(467.29, 2.17, 2, 57.62, 0.2, 5, 1.36, 3.33\right)$ is GAS which confirms Theorem 8. In this case, a persistent dual infection with HTLV and HIV is reached where both the immune response is well working.

6.2.   Comparison study

In this subsection, we compare between single and dual infections dynamics

Influence of HTLV infection on the dynamics of HIV single infection

To study the effect of HTLV infection on the dynamics of HIV single infection, we make a comparison between model (2.1) and the following HIV single infection model:

We fix the parameters $\varkappa_{1} = 0.0006$, $\varkappa_{2} = 0.0001$, $\sigma_{1} = 0.05$, and $\sigma_{2} = 0.5$ and consider initial conditions (6.1) and boundary conditions (6.2). We choose $\varkappa _{3} = 0.04$ (HTLV/HIV dual infection). Figure 3 shows that when an individual who has only HIV infection is dually infected with HTLV then the numbers of uninfected (and latent) CD$4^{+}$T cells and HIV-specific CTLs are declined. In contrast, the numbers of free HIV particles in both HIV single infection and HTLV/HIV dual infection limits to a same value. In fact, this observation is consistent with the recent study ^[73], where it has found that there is no worthy differences in the concentration of HIV particles in comparison between HIV single infected and HTLV/HIV dual infected patients.

Influence of HIV infection on the dynamics of HTLV single infection

To see the effect of HIV infection on the dynamics of HTLV single infection, we perform a comparison between model (2.1) and the following HTLV single infection model:

We fix parameters $\varkappa_{3} = 0.01$; $\sigma_{1} = 0.05$, and $\sigma _{2} = 0.5$ and consider initial conditions (6.1) and boundary conditions (6.2). We choose $\varkappa_{1} =$ $0.001$ and $\varkappa_{2} = 0.0002$ (HTLV/HIV dual infection). Figure 5 displays the solutions of two systems (2.1) and (6.4). We observe that the concentrations of uninfected CD$4^{+}$T cells, latent HTLV-infected cells and HTLV-specific CTLs are smaller in case of dual infection than that of HTLV single infection. In contrast, the concentration of active HTLV-infected cells reaches the same value in both HTLV single and HTLV/HIV dual infections.

7.   Conclusions

This work proposes and investigates a within host HTLV/HIV dual infection model taking into account the mobility of viruses and cells. The model was given by 8-dimentional nonlinear PDEs which describe the evolution of eight compartments with respect to position and time; uninfected CD4$^{+}$T cells, latent HIV-infected cells, active HIV-infected cells, latent HTLV-infected cells, active HTLV-infected cells, free HIV particles, HIV-specific CTLs, and HTLV-specific CTLs. We considered two ways of HIV transmission, free-to-cell and infected-to-cell. We also included two directions of HTLV transmission, horizontal via infected-to-cell contact, and vertical transmission through mitosis of active HTLV-infected cells. We first showed the existence of global solutions and the boundedness of the model's solutions. We showed that the model has eight steady states and their existence and stability are governed by eight threshold parameters. The global asymptotic stability of all steady states was investigated by formulating suitable Lyapunov functions and utilizing Lyapunov-LaSalle asymptotic stability theorem. We conducted some numerical simulations to clearify the theoretical results. We made a comparison between the dynamical behavior of dual HTLV/HIV infection and single HTLV (or HIV) infection. We found that HTLV/HIV dual infected patients have less uninfected CD$4^{+}$T cells counts in comparison with HTLV or HIV single infected patients.

Conflict of interest

The authors declare that they have no conflict interests.