Stability analysis for a HIV model with cell-to-cell transmission, two immune responses and induced apoptosis

1.   Introduction

Human T-cell leukemia virus (HTLV) is a family of retroviruses that predominantly target human T cells, a vital component of the immune system responsible for immune response regulation. HTLV is significant due to its link to various diseases and its capacity to persist in the human body over extended periods. Currently, four strains of this virus HTLV-1, HTLV-2, HTLV-3, and HTLV-4 are known to infect humans. Of these, only HTLV-1 and HTLV-2 have been associated with particular health issues ^[1]. These two strains share similar biological features and modes of transmission. HTLV-1 primarily targets CD4$^{+}$ T cells (also known as helper T cells) and is associated with two major diseases: adult T-cell leukemia (ATL) and a neurological disorder known as HTLV-1-associated myelopathy or tropical spastic paraparesis (HAM/TSP) ^[2]. On the other hand, HTLV-2 primarily infects CD8$^{+}$ T cells, commonly referred to as cytotoxic T lymphocytes (CTLs), which are essential for eliminating cells infected by viruses ^[2]. This strain has been connected to tropical spastic paraparesis as well as peripheral neuropathy ^[2]. People co-infected with HTLV-1 and HTLV-2 frequently face physical symptoms like chronic pain and weakness, along with social and psychological difficulties, such as depression and the stigma associated with carrying these viruses ^[3,4]. The successful transmission of HTLV-1 and HTLV-2 to target cells necessitates direct interaction between cells (cell-to-cell) ^[5]. Both viruses rely on the envelope glycoproteins to mediate cell attachment and entry ^[5]. HTLV is spread through three primary pathways: from parent to child (e.g., during birth or breastfeeding), through parenteral exposure (such as infected blood transfusions, organ transplants, or sharing needles), and via sexual contact ^[5]. As of 2012, it was estimated that between five and ten million people globally were infected with HTLV-1. Regions with high prevalence included South America, the Caribbean, Southwest Japan, sub-Saharan Africa, the Middle East, and Australo Melanesia ^[6]. In contrast, HTLV-2 infections were significantly less common by 2015, with estimates ranging from 670,000 to 890,000 cases ^[7]. Most HTLV-2 cases were reported in the United States, particularly among Native American communities and intravenous drug users. Brazil, the country with the second-highest HTLV-2 prevalence, showed a similar distribution pattern. Infections with these viruses also aggravate other illnesses, such as tuberculosis and strongyloidiasis, which are widespread in endemic regions, adding to the challenges and expenses of their treatment ^[4,8]. Although HTLV-1 and HTLV-2 co-infection is relatively uncommon, it has been identified in certain high-risk populations, especially among individuals who use injection drugs. The health consequences of these co-infections are not yet fully understood, highlighting the need for additional research to clarify their effects on disease progression and overall well-being. At present, no effective treatment exists for infections caused by HTLV-1 and HTLV-2, highlighting the importance of prevention and control measures to limit the dissemination of these viruses and reduce their impact on public health, especially in endemic regions ^[4].

The study of mathematical models offers a powerful and insightful approach to understanding the dynamics of viral infections within a host. This method enhances our comprehension of the mechanisms underlying diseases caused by different viruses. Recently, there has been an increased focus on mathematical models of HTLV-1 dynamics within a host, as they help to uncover the complex interactions between the virus, host cells, and the immune response. The model describing within-host HTLV-1 dynamics, incorporating the CTL-mediated immune response, is expressed as follows ^[9]:

Here,

represent, respectively, the concentrations of healthy CD4$^{+}$ T cells, HTLV-1-infected CD4$^{+}$ T cells, and healthy CD8$^{+}$ T cells at time $t$. This mathematical model has been further expanded and refined by numerous researchers. Lim and Maini ^[10] developed a model to study HTLV-1 dynamics, factoring in CTL response and cell division. Pan et al. ^[11] created a model to describe HTLV dynamics, incorporating CTL response and time delays. Wang et al. ^[12] constructed a detailed HTLV-1 infection model that integrates nonlinear CTL responses (both lytic and nonlytic), a nonlinear incidence rate, distributed delays, and CTL dysfunction. Bera et al. ^[13] examined an HTLV-1 infection model that includes delayed CTL response. In another study, Wang et al. ^[14] investigated HTLV-1 dynamics using a model with two distinct delays: one for intracellular processes and another for CTL immune response. Papers ^[15,16,17] explored HTLV-1 dynamics with CTL response and cell division. Chen et al. ^[18] conducted a global dynamical analysis of an HTLV-1 infection model, incorporating logistic growth of CD4$^{+}$ T cells and nonlinear CTL response. Then the model was extended by taking into account delay CTL response ^[19] and environment noise ^[20]. The model proposed by ^[20] includes the impacts of reverse-transcriptase inhibitors and IL-2 immunotherapy, leading to the determination of an optimal therapeutic strategy. Wang and Ma ^[21] integrated CTL immunity and cell division into a diffusive model of HTLV infection. In contrast, research on HTLV-2 through mathematical modeling has been limited, with much less attention dedicated to its study.

Recent studies have focused on developing mathematical models to examine the co-infection dynamics of HTLV-1 with other viruses, such as HIV-1 (see, for instance, ^{[22,23,24,25,26]}) and SARS-CoV-2 ^[27]. A recent study by ^[28] developed and analyzed a co-infection model for HIV-1 and HTLV-2. As far as we can ascertain from a comprehensive review of the literature, no previous work has performed a dynamical analysis of HTLV-1 and HTLV-2 co-infection. Although various studies have explored HTLV-1 individually, we have not encountered any research that develops and examines a mathematical model specifically addressing the HTLV-1 and HTLV-2 co-infection dynamics.

This research focuses on formulating and examining a mathematical model that captures the dynamics of HTLV-1 and HTLV-2 co-infection within a host. The novelty of this study lies in the following key aspects:

B1. A novel co-infection model has been developed to describe the within-host interactions of HTLV-1 and HTLV-2, capturing their simultaneous presence in the host system.

B2. The model is structured based on the distinct cellular targets of each virus: HTLV-1 primarily infects CD4$^{+}$ T cells, whereas HTLV-2 predominantly targets CD8$^{+}$ T cells.

B3. A rigorous analysis of the model's solutions confirms their non-negativity and boundedness, ensuring both mathematical consistency and biological relevance.

B4. Four threshold parameters are established, which fully determine the conditions for the existence and stability of the model's equilibrium points.

B5. The global stability of each equilibrium is examined using the Lyapunov function approach.

B6. Theoretical results are substantiated through numerical simulations.

B7. A sensitivity analysis is conducted on the basic reproduction numbers for HTLV-1 ($R_{1}$) and HTLV-2 ($R_{2}$), assessing their dependence on key model parameters.

This methodology provides a comprehensive framework for examining the co-dynamics of HTLV variants and their impact on the host. HTLV-1 and HTLV-2 were chosen for this study due to their significant public health implications, their association with severe diseases, and their epidemiological relevance in HIV co-infection. Understanding their dynamics is crucial for enhancing disease control measures and developing effective intervention strategies. Moreover, the proposed model can be adapted to investigate the competitive transmission dynamics of different COVID-19 strains, such as Omicron and Delta.

The structure of this paper is as follows: Sections 2 and 3 introduce the formulation of the co-infection model, analyze the non-negativity and boundedness of its solutions, and derive the equilibrium points along with the threshold parameters. Section 4 explores the global stability of the equilibria. In Section 5, numerical simulations are conducted to validate the theoretical results. Finally, Section 6 provides a summary of the study's findings, discusses their implications, and suggests directions for future research.

2.   Model formulation

This section provides a detailed description of the proposed model. The formulation of the model is based on the following assumptions:

${\rm{A1}}.$ The model represents four populations: healthy CD4$^{+}$ T cells $(U)$, HTLV-1-infected CD4$^{+}$ T cells $(H)$, healthy CD8$^{+}$ T cells $(C)$, and HTLV-2-infected CD8$^{+}$ T cells $(M)$. The compartments $U$, $H$, $C$, and $M$ have respective mortality rates of $\rho_{U}{U}$, $\rho_{H}{H}$, $\rho_{C}{C}$, and $\rho_{M}{M}$.

${\rm{A2}}.$ Healthy CD4$^{+}$ T cells, which serve as the main targets for HTLV-1, are produced at a constant rate $\alpha$. These cells can be infected by HTLV-1 through a cell-to-cell transmission mechanism at a rate of $\eta_{UH}{U}{H}$ ^[9] (see Eq (2.1)).

${\rm{A3}}.$ HTLV-1-infected CD4$^{+}$ T cells are generated at a rate of $\eta_{UH}{U}{H}$. These infected cells are eliminated through immune-mediated destruction by CD8$^{+}$ T cells, at a rate represented by $\chi HC$ ^[9] (see Eq (2.2)).

${\rm{A4}}.$ Healthy CD8$^{+}$ T cells, the primary targets of HTLV-2, are produced at a constant rate $\gamma$ (self-regulating immune response) and proliferate in response to the presence of HTLV-1-infected CD4$^{+}$ T cells at a rate of $\sigma{H}{C}$ (predator prey-like immune response) ^[9]. These cells can become infected through direct contact with HTLV-2-infected CD8$^{+}$ T cells via a cell-to-cell mechanism, at a rate described by $\eta_{CM}{C}{M}$ ^[5,28] (see Eq (2.3)).

${\rm{A5}}.$ HTLV-2-infected CD8$^{+}$ T cells are generated at a rate of $\eta_{CM}{C}{M}$ ^[28] (see Eq (2.4)).

Figure 1 illustrates the diagram that outlines the dynamics of HTLV-1 and HTLV-2 co-infection.

According to assumptions A1–A5, the structure of the proposed models (2.1)–(2.4) is outlined as follows:

The descriptions and values of the variables and parameters are shown in Table 1.

Remark 1. We mention that we have used bilinear incidence rate for the infection rate. This form is commonly used in mathematical virology due to its simplicity and analytical tractability while capturing the fundamental interaction between target cells and viruses or infected cells (see, e.g., ^[9,10,11,17]). Specifically, the terms $\eta_{UH}UH$ and $\eta_{CM}CM$ represent the infection rates proportional to the product of interacting populations, aligning with the standard mass-action principle. Some modifed versions for the infection rate have been included in the HTLV-1 models, such as:

● Saturated incidence: $\frac{\eta_{{U}{H}}{U}{H}}{1+{mH}}$, where $m\geq0$. This incidence form can be suitable when the concentration of infected T cells becomes significantly high ^[35].

● Holling type-Ⅰ: $\frac{\eta_{{U}{H}}{U}{H}}{1+n{U}}$ ^[16].

● Beddington-DeAngelis incidence: $\frac{\eta_{{U}{H}}{U}{H} }{1+n{U+mH}}$, where $n, m\geq0$ ^[36].

● Crowley Martin incidence: $\frac{\eta_{{U}{H}}{U}{H} }{(1+n{U)(1+mH)}}$ ^[12].

● Nonlinear incidence: $\eta_{{U}{H}}U^{p}H^{q}$, where $p$ and $q$ are positive constants ^[21]. Another form, presented in ^[37,38] as $\Psi(H)U$, where $\Psi$ is a nonlinear function. In ^[12], the bilinear incidence $\eta_{UH} UH$ was adjusted to $\frac{\eta_{UH}UH}{1+q_{0}C}$ to incorporate the impact of the non-lytic CTL (NL-CTL) response, which suppresses viral replication through soluble mediators. Here, $q_{0}$ represents the effectiveness of the NL-CTL response.

● Standard incidence: $\frac{\eta_{{U}{H}}{U}{H}}{{U+H}}$, ^[39]. This form has been generalized in ^[40] as $\frac{\eta_{{U}{H}}{U}{H}}{({U+H)} ^{\varepsilon}}$, where $0\leq\varepsilon\leq$1.

● General incidence: $\Psi(U, H)H$, where $\Psi (U, H)$ is a general function that is the average number of healthy CD4$^{+}$ T cells that are infected by unit HTLV-1-infected CD4 $^{+}$ T cells per unit time ^[21]. In ^[41], a general incidence in the form $\Psi(U, H)$ is considered.

Due to the limited experimental data specifically quantifying HTLV-1 and HTLV-2 co-infection transmission rates, we adopted this bilinear incidence as a reasonable first approximation.

3.   Preliminaries

This section analyzes the essential qualitative characteristics of the systems (2.1)–(2.4), including the solutions' nonnegativity and boundedness of solutions. Moreover, each equilibrium is identified together with its corresponding threshold number.

3.1.   Non-negativity and boundedness of solutions

In this subsection, we establish the well-posedness of the models (2.1)–(2.4) by demonstrating that the solutions remain nonnegative and bounded over time.

Lemma 1. Solutions of the systemS (2.1)–(2.4) are nonnegative and bounded.

Proof. We have

Hence, in conformity with [42, Proposition B.7]

for any $t\geq0$, when

To illustrate the boundedness of the solutions. Let's establish the definition of $\psi(t)$ as:

Next, we obtain

where

Thus,

if

Consequently

if

where

This completes the proof. □

3.2.   Equilibria and thresholds

Define

and

and introduce the following four indices, which will serve as threshold parameters, denoted by

and given by

It is crucial to highlight that $R_{1}$ indicates the basic reproduction number for HTLV-1 single infection and represents the number of new CD4$^{+}$ T cells infected with HTLV-1 that are generated from a single HTLV-1-infected CD4$^{+}$ T cell during the early phase of HTLV-1 single infection. The parameter $R_{2}$ stands for the basic reproduction number for HTLV-2 single infection and represents the quantity of newly infected CD8$^{+}$ T cells carrying HTLV-2 that originate from one HTLV-2-infected CD8$^{+}$ T cell in the initial phases of HTLV-2 single infection. The threshold values $R_{3}$ and $R_{4}$ serve as indicators of the likelihood of co-infection with HTLV-1 and HTLV-2.

Lemma 2. For the models (2.1)–(2.4), four equilibrium points (EPs) exist, such that

${\bf{(I)}}$ Disease-free equilibrium, $EP_{0}$, is always presented, where

${\bf{(II)}}$ If $R_{1} > 1$, an equilibrium for HTLV-1 single infection,

will emerge in addition to $EP_{0}$.

${\bf{(III)}}$ If $R_{2} > 1$, an equilibrium for HTLV-2 single infection,

will emerge in addition to $EP_{0}$.

${\bf{(IV)}}$ If $R_{3} > 1$ and $R_{4} > 1$, an equilibrium HTLV-1/HTLV-2 co-infection,

will emerge in addition to $EP_{0}$.

Proof. The EPs of models (2.1)–(2.4) satisfy the following:

Solving the algebraic system (3.2), we obtain four equilibrium points as follows:

(1) Disease-free equilibrium,

(2) HTLV-1 single infection equilibrium,

where

and $H_{1}$ fulfills the following:

where

Let us define a function $\mathcal{D}(H)$ as follows:

Note that, $\mathcal{D}$ is continuous on $\left[0, \frac{\rho_{C}}{\sigma }\right)$. We have

The presence of an HTLV-1 single infection is determined by evaluating the parameter $R_{1}$. Since

if $R_{1} > 1$ in addition to

there exists $H_{1}$ such that

and satisfies

Consequently, we obtain ${U}_{1} > 0, \text{ and }{C}_{1} > 0$.

(3) HTLV-2 single infection equilibrium,

where

The persistence of an HTLV-2 single infection can be ascertained by assessing the parameter $R_{2}$.

(4) HTLV-1/HTLV-2 co-infection equilibrium,

where

This completes the proof. □

4.   Global stability

This section aims to examine the global asymptotic stability of all EPs in the models (2.1)–(2.4) using the Lyapunov approach, as proposed in the work conducted by ^[43]. We define a function

Moreover, the arithmetic mean-geometric mean inequality is utilized to establish the proofs of Theorems 1–4 as:

Let $\mathcal{L}_{k}$ be the candidate Lyapunov function and define $\mathcal{H}_{k}^{\prime}$ as the largest invariant set of

Theorem 1. The disease-free equilibrium $EP_{0}$ is globally asymptotically stable (GAS) if $R_{1}\leq1$ and $R_{2}\leq 1$; otherwise, it is unstable.

Proof. To verify (ⅰ), we define $\mathcal{L}_{0}({U}, {H}, {C}, {M})$ as:

Obviously,

for any

and

The derivative of $\mathcal{L}_{0}$ along the solutions of systems (2.1)–(2.4) can be computed as:

By substituting the equations stated in models (2.1)–(2.4), we derive

By gathering the terms and substituting $\alpha$ with $\rho_{U} {U}_{0}$ and $\gamma$ with $\rho_{C}{C}_{0}$, we derive the following expression:

Then,

Ultimately, we obtain

Thus,

satisfies if $R_{1} \leq1\text{ and }R_{2}\leq1$. Moreover,

when

The solutions of the system approach $\mathcal{H} _{0}^{\prime}$ ^[44]. Every element in $\mathcal{H}_{0}^{\prime}$ satisfies $U = U_{0}$, $C = C_{0}$,

As a result, four cases arise:

${\bf{(I)}}$ $R_{1} = 1$ and $R_{2} = 1$. Then from Eq (2.1) we obtain

From Eq (2.3) we have

Hence,

${\bf{(II)}}$ $R_{1} < 1$ and $R_{2} < 1$. Then, Eq (4.2) implies that $H = M = 0$ and, hence,

${\bf{(III)}}$ $R_{1} = 1$ and $R_{2} < 1$. Then, Eq (4.2) suggests $M = 0$ and Eq (4.3) gives $H = 0$. Thus,

${\bf{(IV)}}$ $R_{1} < 1$ and $R_{2} = 1$. Eq (4.2) gives $H = 0$ while Eq (4.4) implies $M = 0$. Consequently,

By LaSalle's invariance principle (LP) ^[45], $EP_{0}$ is GAS.

Now, let us establish the instability of $EP_{0}$ if $R_{1} > 1$ and/or $R_{2} > 1$. First of all, we construct the Jacobian matrix

of models (2.1)–(2.4) as:

Next, we calculate the characteristic equation at the equilibrium point $EP_{0}$ as:

where $\Lambda$ represents the eigenvalue, $I$ represents the identity matrix, and

Then matrix $J$ has the following eigenvalues:

Clearly, $\Lambda_{3} > 0$ and $\Lambda_{4} > 0$ , when $R_{1} > 1$ and $R_{2} > 1$ . It follows that if either $R_{1} > 1$ , $R_{2} > 1$ , or both, then $EP_{0}$ is unstable.

This completes the proof. □

Theorem 2. HTLV-1 single infection equilibrium $EP_{1}$ is GAS if $R_{1} > 1$ and $R_{4}\leq1$.

Proof. Construct $\mathcal{L}_{1}({U}, {H}, {C}, {M})$ as:

Clearly,

for any

and

Calculating $\frac{d\mathcal{L}_{1}}{dt}$ as:

Collecting terms results to

By using the subsequent equilibrium conditions

we obtain

It follows that

We will now demonstrate that if

as:

Thus, by applying inequality (4.1), it follows that

for all

In addition,

if

The solutions of models (2.1)–(2.4) approach $\mathcal{H}_{1}^{\prime}$, where

and

As a result, two cases arise:

${\bf{(I)}}$ $C_{1} = C_{3}$, then from Eq (2.1)

From Eq (2.3), we have

Then,

${\bf{(II)}}$ $C_{1} < C_{3}$, then Eq (4.7) implies that $M = 0$ and Eq (4.8) gives $H = H_{1}$. Hence,

Thus, by LP, $EP_{1}$ is GAS. □

Theorem 3. HTLV-2 single infection equilibrium, $EP_{2}$ is GAS if $R_{2} > 1$ and $R_{3}\leq1$, and if $R_{3} > 1$, then $EP_{2}$ is unstable.

Proof. Define $\mathcal{L}_{2}({U}, {H}, {C}, {M})$ as:

Evidently,

for any

and

Calculating $\frac{d\mathcal{L}_{2}}{dt}$ as:

Collecting the above terms leads to

Utilizing the equilibrium conditions

We obtain

Then, if $R_{3}\leq1$ and by using inequality (4.1), we obtain that

for any

Moreover,

if

The model's solutions converge to $\mathcal{H}_{2}^{\prime}$ where

and

As a result, two cases arise:

${\bf{(I)}}$ $R_{3} = 1$. From Eq (2.1), we have

Equation (2.3) implies that

Hence,

${\bf{(II)}}$ $R_{3} < 1$. Then, from Eq (4.10), we obtain $H = 0$. In addition, Eq (4.12) gives $M = M_{2}$ and, hence,

Consequently, by LP, $EP_{2}$ is GAS.

We demonstrate that if $R_{3} > 1$, then $EP_{2}$ is unstable. By applying the Jacobian matrix in Eq (4.5), we can derive the characteristic equation at the equilibrium $EP_{2}$ as follows:

$\Gamma$ is continuous on $[0, \infty)$,

and

Since $R_{2} > 1$ and $R_{3} > 1$, then $\Gamma(0) < 0$. Hence, $\Gamma(\Lambda)$ has a positive root, and thus $EP_{2}$ is unstable.

This completes the proof. □

Theorem 4. HTLV-1/HTLV-2 co-infection equilibrium, $EP_{3}$ is GAS if $R_{3} > 1$ and $R_{4} > 1$.

Proof. Define $\mathcal{L}_{3}({U}, {H}, {C}, {M})$ as:

Calculating $\frac{d\mathcal{L}_{3}}{dt}$ as:

Then we obtain

Using the equilibrium conditions

We obtain

It follows that

Therefore, if $R_{3} > 1$ and $R_{4} > 1$, and by applying inequality (4.1), we deduce that

for all

In addition,

if

Solutions of models (2.1)–(2.4) converge to $\mathcal{H^{\prime}}_{4}$, where

It follows that

and then Eqs (2.1) and (2.3) give

Thus, by LP

and $EP_{3}$ is GAS.

This completes the proof. □

Table 2 summarizes the conditions for the existence and global stability of each equilibrium in the models (2.1)–(2.4).

Remark 2. Investigating the memory effects and hereditary properties on our model's behavior using fractional differential equations (FDEs) offers a compelling research direction. FDEs are particularly effective in representing memory effects and non-local interactions, which are crucial in biological ^[46] and epidemiological systems ^[47,48]. While our study focuses on the integer-order version, extending the model to a fractional-order system could offer deeper insights into the long-term behavior of the system, particularly in capturing the persistence of infection or addiction-related processes. In recent years, the Lyapunov method has gained significant attention for studying the global stability of FDE models ^[49]. The Lyapunov functions constructed in this section are essential for examining a fractional model describing HTLV-1 and HTLV-2 co-dynamics.

5.   Numerical simulations

The suggested model systems (2.1)–(2.4) are subjected to numerical simulation in order to acquire insights into dynamical features, hence enhancing our comprehension of how control approaches impact the dynamics of infectious disease transmission. First, the stability characteristics of the infectious model are examined.

5.1.   Numerical simulations for systems (2.1)–(2.4)

This subsection provides numerical simulations based on the parameter values listed in Table 1 to visually represent the analytical results established in Theorems 1–4. We employed the ode45 solver in MATLAB to numerically solve the ordinary differential equation model. To establish the global stability of the system's equilibria, we show that its solutions converge to a specific equilibrium point, irrespective of the initial conditions. Given the scarcity of real-world data on initial infections with HTLV-1 and HTLV-2, we consider three distinct initial values as follows:

We perform the numerical simulation by choosing the parameters $\eta_{UH}$ and $\eta_{CM}$ in the following scenarios:

${\bf{Scenario-1:}}$ $\eta_{UH}$ and $\eta_{CM}$ are set at $0.005$ and $0.0005$, respectively. Our findings indicate that with these values

The paths shown in Figure 2 originating from the three points all converge to the equilibrium point

This demonstrates that $EP_{0}$ is GAS as provided in Theorem 1. Consequently, the HTLV-1 and HTLV-2 will be eradicated by this process.

${\bf{Scenario-2:}}$ $\eta_{UH}$ and $\eta_{CM}$ are set at $0.01$ and $0.0005$, respectively. This means that

and

The results depicted in Figure 3 illustrate how the solutions reach the equilibrium point

In this scenario, Theorem 2 aligns with the findings. This situation exemplifies the effects of HTLV-1 infection without HTLV-2 on an individual health. It is evident that HTLV-1 has caused a decrease in CD4$^{+}$ T cell counts and an increase in CD8$^{+}$ T cell levels.

${\bf{Scenario-3:}}$ The values of $\eta_{UH}$ and $\eta_{CM}$ are 0.001 and 0.002, respectively. Next, we calculate

and

It is clear that the criteria outlined in Table 2 are satisfactorily met. In Figure 4 we can observe how the solutions converge towards the equilibrium point

thereby confirming Theorem 3. This scenario illustrates the progression when an individual is solely infected with HTLV-2. Even though CD4$^{+}$ T cell levels remain normal, it is evident that HTLV-2 has caused a decrease in CD8$^{+}$ T cell counts.

${\bf{Scenario-4:}}$ $\eta_{UH}$ is $0.003$, while $\eta_{CM}$ is $0.004$. After that, we calculate

and

The conditions specified in Table 2 are clearly met. As depicted in Figure 5, the solutions approach the equilibrium

and confirm Theorem 4. This scenario illustrates the consequences of co-infection with both HTLV-1 and HTLV-2 viruses in an individual. The quantities of healthy CD4$^{+}$ T cells and CD8$^{+}$ T cells decline, resulting in a weakening of the patient's immune system.

To provide additional verification, an analysis of the local stability is conducted for all equilibria. The Jacobian matrix $\mathcal{J}$, which depends on the variables $U$, $H$, $C$, and $M$, is calculated in Eq (4.5). To determine the stability of each equilibrium, we compute the eigenvalues $\lambda_{k}$, where $k = 1, \ldots, 4$ of $\mathcal{J}$. An equilibrium point is considered stable if all eigenvalues have a real part less than zero, denoted as

for all $k = 1, 2, 3, 4$. By calculating the nonnegative EPs and utilizing the parameter values specified in Scenarios 1–4, we can determine the eigenvalues associated with all equilibria. Table 3 provided a comprehensive overview of the positive equilibria and the real part of the eigenvalues.

5.2.   Sensitivity analysis

Sensitivity analysis is important for understanding systems by showing how model outcomes are affected by different parameters ^[50]. It helps enhance our grasp of model research. Sensitivity analysis evaluates how biological reactions change when parameters are adjusted, helping pinpoint the factors influencing model results ^[51]. Several methods of sensitivity analysis were introduced for biological models ^[51]. We apply derivative-based sensitivity analysis to our system. By computing the derivatives with respect to model parameters, we can analytically determine the indices. In our study, we examine the effect of $R_{1}$ and $R_{2}$ on the stability of disease-free infection equilibrium using sensitivity analysis. The normalized forward sensitivity index for $R_{i}$ (where $k = 1, 2$) is expressed as follows:

where $\beta$ represents a given parameter.

5.2.1.   Sensitivity analysis for $R_{1}$

Applying form (5.1) and using the parameter values provided in Table 1, we calculated the sensitivity indices of $R_{1}$ with respect to each parameter, as presented in Table 4.

The sensitivity index values for $R_{1}$ with respect to the parameters are given in Table 1. Based on the presented sensitivity indices in Table 4, we observe that:

● The parameters $\alpha$, $\eta_{UH}$, and $\rho_{C}$ show positive sensitivity indices, suggesting that changes in these parameters will lead to corresponding increases or decreases in the basic reproduction number $R_{1}$ for HTLV-1 mono-infection. Among them, $\alpha$ and $\eta_{UH}$ exhibit the highest positive sensitivity indices.

● On the other hand, $\rho_{U}$, $\chi$, $\gamma$, and $\rho_{H}$ have a negative influence on $R_{1}$. As a result, an increase in these values will lead to a reduction in the value of $R_{1}$. Among these parameters, $\rho_{U}$, $\chi$, and $\gamma$ show greater significance compared to $\rho_{H}$.

5.2.2.   Sensitivity analysis for $R_{2}$

Using Eq (5.1), we can calculate the sensitivity indices of $R_{2}$ with respect to each parameter as shown in Table 5. It is clear that the sensitivity analysis of $R_{2}$ is independent of the parameter values, as the sensitivity to any parameter is consistently either $1$ or $-1$. The signs presented in Table 5 help us understand how each parameter contributes to the sensitivity analysis as:

● The values of parameters $\gamma$ and $\eta_{CM}$ positively influence the growth of HTLV-2 in the body, suggesting their role in its proliferation.

● In contrast, parameters $\rho_{C}$ and $\rho_{M}$ are associated with reducing the transmission rate of HTLV-2 within humans.

5.3.   HTLV-1 and HTLV-2 co-dynamics under the effect of stimulation rate of CD8$^{+}$ T cells

This subsection examines the impact of the stimulated rate constant of CD8$^{+}$ T cells, denoted as $\sigma$, on the system dynamics described by (2.1)–(2.4). In order to investigate the impact of CD8$^{+}$ T cells on the model's solutions, we hold the values of

and

while varying the parameter $\sigma$. By choosing the following initial point:

From Figure 6, we observe that as $\sigma$ increases, the quantities of healthy CD4$^{+}$ T and CD8$^{+}$ T cells remain increased. Moreover, the number of HTLV-1-infected CD4$^{+}$ T cells decreases.

It is also important to note that an increase in CD8$^{+}$ T cells will increase the number of HTLV-2-infected CD8$^{+}$ T cells. Clearly, enhanced CD8$^{+}$ T cell stimulation helps control HTLV-1 infection while also promoting the spread of HTLV-2 infection. Due to the fact that $R_{1}$ and $R_{2}$ are independent of $\sigma$, increasing $\sigma$ does not result in the attainment of $EP_{0}$. To boost CD8$^{+}$ T cell stimulation, the parameter $\sigma$ can be substituted with $(1+\epsilon_{\text{{IT}} })\sigma$, where $\epsilon_{\text{{IT}}}\in\lbrack 0, 1]$ denotes the drug efficacy of immunotherapy (IT), which enhances the proliferation of CD8$^{+}$ T cells in infected individuals through the subcutaneous administration of IL-2, promoting their activation and differentiation ^[20].

5.4.   Impact of infection rates $\eta_{UH}$ and $\eta_{CM}$ on the HTLV1 and HTLV-2 co-dynamics

In this part, we explore the impact of infection rates $\eta_{UH}$ and $\eta_{CM}$ on the HTLV1 and HTLV-2 co-dynamics. To demonstrate the impact of $\eta_{UH}$ on the HTLV1 and HTLV-2 co-dynamics, we set

We consider different values of $\eta_{UH}$ as

and numerically solve systems (2.1)–(2.4) together with the following initial point:

Figure 7 reveals a correlation in which an increase in $\eta_{UH}$ leads to a reduction in the number of healthy CD4$^{+}$ T cells, while the behavior of healthy CD8$^{+}$ T cells does not change too much. At the same time, populations of HTLV-1-infected CD4$^{+}$ T cells and HTLV-2-infected CD8$^{+}$ T cells increase. In this context, raising $\eta_{UH}$ may increase the risk of HTLV-2 infection. To reduce the HTLV-1 infection, the parameter $\eta_{UH}$ can be replaced with $(1-\epsilon_{\text{RTI} })\eta_{UH}$, where $\epsilon_{\text{RTI}}\in\lbrack0, 1]$ represents the drug efficacy of reverse-transcriptase inhibitors (RTI) like zidovudine, which effectively block the virus from spreading throughout the body ^[20,52].

To demonstrate the effect of $\eta_{CM}$ on the HTLV1 and HTLV-2 co-dynamics, we set

We examine different values of

and numerically solve systems (2.1)–(2.4) with the following initial point:

Analysis of Figure 8 shows a correlation where an increase in $\eta_{CM}$ results in a decrease in both healthy CD4$^{+}$ T cells and healthy CD8$^{+}$ T cells. At the same time, the populations of HTLV-1-infected CD4$^{+}$ T cells and HTLV-2-infected CD8$^{+}$ T cells increase. In this context, raising $\eta_{CM}$ may increase the HTLV-1 progression.

6.   Discussion

Since HTLV-1 and HTLV-2 are closely linked to various health concerns and share transmission routes, they are especially common in high-risk populations, particularly among individuals who use injection drugs. Therefore, understanding the within-host co-dynamics of HTLV-1 and HTLV-2 is crucial. In this paper, we develop a model to describe the dynamics of HTLV-1 and HTLV-2 co-infection. A summary of the equilibria in our model is as follows:

$({\rm{I}})$ Disease-free equilibrium. $EP_{0}$: this equilibrium always exists and is GAS when

In this scenario, both HTLV-1 and HTLV-2 will be eradicated from the system. While there is no definitive cure for HTLV infections, treatment mainly focuses on managing associated conditions such as ATL and HAM/TSP using chemotherapy, corticosteroids, interferon, and supportive care. If antiviral drugs for HTLV-1 and HTLV-2 become available, it may be possible to achieve

by adjusting the parameters $\eta_{UH}$ and $\eta_{CM}$. Specifically, these parameters can be replaced with $\left(1-\epsilon_{\text{RTI-HTLV-1}}\right) \eta_{UH}$ and $\left(1-\epsilon_{\text{RTI-HTLV-2}}\right) \eta_{UH}$ where $\epsilon_{\text{RTI-HTLV-1}}\in\lbrack0, 1]$, and $\epsilon_{\text{RTI-HTLV-2}}\in\lbrack0, 1]$ represent the efficacy of RTI drugs for HTLV-1 and HTLV-2, respectively. Then $R_{1}\ $ and $R_{2}\ $become

Therefore, $R_{1} \leq1$ when

$R_{2} \leq1$ when

where

In this context, $\epsilon_{\text{RTI-HTLV-1}}^{\min}$ and $\epsilon_{\text{RTI-HTLV-2}}^{\min}$ represent the minimum drug efficacies necessary to stabilize $EP_{0}$ and eliminate the co-infection from the host. By treating drug efficacies as control variables, optimal control theory can be applied to develop treatment strategies that minimize both the cost of the drugs and their associated side effects ^[24].

$({\rm{II}})$ HTLV-1 single infection equilibrium point. $EP_{1}$: this equilibrium exists if $R_{1} > 1$ and is GAS if $R_{4}\leq1$. This scenario represents an individual infected only with HTLV-1. This likely occurs when $\epsilon_{\text{RTI-HTLV-1}}$ is insufficient to eliminate the HTLV-1 infection, whereas $\epsilon_{\text{RTI-HTLV-2}}$ is effective in clearing HTLV-2.

$({\rm{III}})$ HTLV-2 single infection equilibrium point. $EP_{2}$: this equilibrium exists if $R_{2} > 1$ and is GAS if $R_{3}\leq1$. It describes a situation where an individual is infected only with HTLV-2. This situation may arise when $\epsilon_{\text{RTI-HTLV-2}}$ fails to eradicate HTLV-2, whereas $\epsilon_{\text{RTI-HTLV-1}}$ is sufficient to eliminate HTLV-1.

$({\rm{IV}})$ HTLV-1 and HTLV-2 co-infection equilibrium point. $EP_{3}$: this equilibrium point exists and is GAS if

It represents an individual co-infected with both HTLV-1 and HTLV-2. Our findings indicate that co-infection with both HTLV-1 and HTLV-2 leads to a decrease in the number of healthy CD4$^{\mathbf{+}}$ and CD8$^{\mathbf{+}}$ T cells, while the number of HTLV-1-infected CD4$^{\mathbf{+}}$ T cells and HTLV-2-infected CD8$^{\mathbf{+}}$ T cells increases, ultimately compromising the patient's immune response. This may elevate the risk associated with both infections.

A key limitation of our study is the inability to estimate the model's parameter values using real-world data. This challenge arises due to several factors: (ⅰ) the scarcity of available data on HTLV-1 and HTLV-2 co-infection; (ⅱ) the limited number of studies on this topic, making comparisons less reliable; and (ⅲ) the difficulty in obtaining clinical data from patients infected with both viruses.

7.   Conclusions and future perspectives

In this study, we examined a mathematical model that captures the population dynamics of two HTLV strains, HTLV-1 and HTLV-2, which infect distinct target cells specifically CD4$^{+}$ T cells and CD8$^{+}$ T cells, both crucial to the human immune system. Our model captures the interactions among four components: healthy CD4$^{+}$ T cells, healthy CD8$^{+}$ T cells, HTLV-1-infected CD4$^{+}$ T cells, and HTLV-2-infected CD8$^{+}$ T cells. We have provided preliminary results regarding the boundedness and non-negativity of the solutions to the model. We then determined that the models possess four EPs. Through the application of LaSalle's invariance principle and the construction of suitable Lyapunov functions, we identified four threshold parameters ($R_{1}, R_{2}, R_{3}, R_{4}$) that determine the global stability of the equilibrium points. To validate these theoretical results, numerical simulations were carried out. We found, a close match between the numerical and theoretical outcomes. A sensitivity analysis was conducted to evaluate the impact of the model's parameters on the basic reproduction numbers $R_{1}$ and $R_{2}$. The key findings emerged from the analysis:

● The model examined encompasses several clinical scenarios, including a patient who has recovered from both HTLV-1 and HTLV-2 infections, an individual with chronic HTLV-1 mono-infection, another with chronic HTLV-2 mono-infection, and a patient with persistent co-infection of both viruses.

● Enhanced stimulation of CD8$^{+}$ T cells helps control HTLV-1 infection, while simultaneously promoting the spread of HTLV-2 infection. Numerous studies provide valuable insights into diverse approaches for boosting CD8$^{\mathbf{+}}$ T cell activation and their potential use in therapy (see, e.g., ^[53,54]). CD8$^{+}$ T cell stimulation can be enhanced through immunotherapy, which involves promoting the proliferation of CD8$^{+}$ T cells in infected individuals via subcutaneous injection of IL-2, thereby activating and differentiating the T cells ^[20].

● The coinfection with HTLV-1 and HTLV-2 may increase the risk of both viruses. Immune system dysregulation caused by both viruses could contribute to a higher inflammatory state, possibly worsening neurological symptoms.

Modeling the interactions between HTLV-1 and HTLV-2 has provided crucial insights into their pathogenesis and has also played a role in shaping guidelines for more effective treatment strategies for co-infection.

Our proposed model can be expanded in several directions:

(ⅰ) Incorporating a generalized incidence function $\Psi(U, H)$ that encompasses various forms, such as saturated incidence, Beddington-DeAngelis incidence, and Crowley-Martin incidence ^[41];

(ⅱ) Extending the model using partial differential equations to account for cellular mobility;

(ⅲ) Employing fractional differential equations to capture the influence of immunological memory ^[48]. Future research could explore incorporating the impact of various drug therapies into the model.

Moreover, we aim to compare the model's outcomes with data from infected patients to validate its predictions.

Author contributions

E. A. Almohaimeed: conceptualization, formal analysis; A. M. Elaiw: investigation, formal analysis; A. D. Hobiny: methodology, writing original draft. All authors have read and agreed to the published version of the manuscript.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The researchers would like to thank the Deanship of Graduate Studies and Scientific Research at Qassim University for financial support (QU-APC-2025).

Conflict of interest

The authors declare that they have no conflicts of interest.