Stability of an HTLV-HIV coinfection model with multiple delays and CTL-mediated immunity.

In the literature, several mathematical models have been formulated and developed to describe the within-host dynamics of either human immunodeficiency virus (HIV) or human T-lymphotropic virus type I (HTLV-I) monoinfections. In this paper, we formulate and analyze a novel within-host dynamics model of HTLV-HIV coinfection taking into consideration the response of cytotoxic T lymphocytes (CTLs). The uninfected CD 4 + T cells can be infected via HIV by two mechanisms, free-to-cell and infected-to-cell. On the other hand, the HTLV-I has two modes for transmission, (i) horizontal, via direct infected-to-cell touch, and (ii) vertical, by mitotic division of active HTLV-infected cells. It is well known that the intracellular time delays play an important role in within-host virus dynamics. In this work, we consider six types of distributed-time delays. We investigate the fundamental properties of solutions. Then, we calculate the steady states of the model in terms of threshold parameters. Moreover, we study the global stability of the steady states by using the Lyapunov method. We conduct numerical simulations to illustrate and support our theoretical results. In addition, we discuss the effect of multiple time delays on stability of the steady states of the system.

Introduction

During the past several decades many human viruses and their associated diseases, such as human immunodeficiency virus (HIV), hepatitis C virus (HCV), hepatitis B virus (HBV), dengue virus, human T-lymphotropic virus type I (HTLV-I), and recently coronavirus, have been recognized. Human body can be infected by more that one virus at the same time such as HTLV-HIV, coronavirus/influenza, HCV-HIV, HBV-HIV, HCV-HBV, and malaria-HIV. HTLV and HIV are universal public health matters. HTLV and HIV are two viruses which infect most effective immune cells, cells. Adult T-cell leukemia (ATL) and HTLV-I-associated myelopathy/tropical spastic paraparesis (HAM/TSP) are the last stage of HTLV-I infection. Chronic HIV infection leads to acquired immunodeficiency syndrome (AIDS). Both HTLV and HIV have the same ways of transmission such as sharing contaminated needles and unprotected sexual contact with infected partners. Over the last 10 years HTLV-HIV coinfection has been widely documented (see e.g. [1-3], and [4]).

Mathematical models

Mathematical models of HIV and HTLV-I dynamics have become efficient tools to biological and medical scientists. These models can provide a deeper understanding of within-host virus dynamics and assist in predicting the impact of antiviral drug efficacy on viral infection progression (see e.g. [5-16]).

HIV monoinfection model: The standard HIV dynamics model under the effect of cytotoxic T lymphocytes (CTLs) has been formulated by Nowak and Bangham [17] as follows: where , , , and are the concentrations of uninfected cells, active HIV-infected cells, free HIV particles, and HIV-specific CTLs, respectively, and t is the time. η refers to the generation rate of the uninfected cells. The uninfected cells are infected via HIV particles (free-to-cell infection) at rate . The HIV-infected cells produce HIV particles at rate bI. The stimulation rate of effective HIV-specific CTLs due to the presence of HIV-infected cells is defined by . The term accounts for the killing rate of HIV-infected cells due to its specific CTLs. The four compartments S, I, V, and have normal death rates ϱS, aI, εV, and , respectively. Several extensions on model (1) have been accomplished (see e.g. [18-20]).

HTLV-I monoinfection model: The within-host dynamics of HTLV-I has been mathematically modeled in several papers [21-24]. CTL immunity has been included into the HTLV-I dynamics models in many works (see e.g. [25-33]). Lim and Maini [28] have formulated a model for HTLV-I dynamics under the consideration of CTL immunity and mitotic division of active HTLV-infected cells as follows: where , , and are the concentrations of latent HTLV-infected cells, active HTLV-infected cells, and HTLV-specific CTLs at time t, respectively. The term denotes the infected-to-cell contact rate between HTLV-infected cells and uninfected cells (horizontal transmission). The active HTLV-infected cells transmit vertically to latent compartment at rate (mitotic transmission), where . The HTLV-specific CTLs kill the active HTLV-infected cells at rate and are stimulated at rate . The term ψE denotes the activation rate of latent HTLV-infected cells. The death rates of E, Y, and are given by ωE, , and , respectively.

HTLV-HIV coinfection model: Elaiw and AlShamrani [34] have recently formulated an HTLV-HIV coinfection model as follows: where is the concentration of latent HIV-infected cells. The term describes the infection rate of uninfected cells by HIV-infected cells. λL and γL are the activation and death rates of latent HIV-infected cells. The parameter represents the part of newly HIV-infected cells that becomes active, and the other part enters a latent stage. The parameter refers to the part of newly HTLV-infected cells that become latent.

Intracellular delay plays a crucial role in within-host virus dynamics and is defined as the time lapse between viral entry a cell and its production. In case of HIV, it has been estimated that the time between the HIV enters a target cell until producing new HIV particles is about 0.9 days [35]. Time delay has also an important effect in HTLV-I infection. Several works have been devoted to developing mathematical models with time delays to describe the dynamics of HIV (see e.g. [36-43]) and HTLV (see e.g. [44-53]).

Our aim is to take model (3) to further destination by incorporating multiple intracellular time delays and mitotic transmission. We study the fundamental and global properties of the system, then we present numerical simulation. The outcomes of this paper will help clinicians to estimate the suitable time to start the treatment. Our model may be helpful to study different coinfections such as influenza-coronavirus, HCV-HIV, HBV-HIV, and malaria-HIV. It is interesting to note that fractional-order differential equations (FODEs) have been widely studied in several works (see e.g. [54-57]). Modeling and analysis of HIV dynamics with FODEs have been investigated in many papers (see e.g. [58-60]). Clinicians can use the information (in terms of behavior predictions) of fractional-order systems to fit patients’ data with the most appropriate noninteger-order index. As a future work, our coinfection model can be formulated as a system of FODEs.

The multiple delays model

In this section, we extend system (3) by taking under consideration multiple types of distributed-time delays and mitosis of active HTLV-infected cells. We achieve this goal by considering the following system of delay differential equations (DDEs): The factor represents the probability that uninfected cells contacted by HIV particles or active HIV-infected cells at time survived time units and become latent infected at time t. The term is the probability that uninfected cells contacted by HIV particles or active HIV-infected cells at time survived time units and become actively infected at time t. The term is the probability that latent HIV-infected cells survived time units before transmitted to be active at time t. Moreover, the factor demonstrates the probability that the initial infection of uninfected cells and the HTLV-infected cells at time completing all the intracellular processes that are required for it to become latent HTLV-infected cells at time t. Further, the probability that latent HTLV-infected cells survived time units before transmitted to active HTLV-infected cells at time t is given by the factor . Furthermore, the term refers to the probability that new immature HIV particles at time lost time units and become mature at time t. Here , , are positive constants. The delay parameter is randomly taken from a probability distribution function over the time interval , , where is the limit superior of this delay period. The function , satisfies and where . Let us denote where . Thus , .

According to [28], we assume that . This yields . Let and . Then system (4) becomes The initial conditions of system (5) are given by: where , , and is the Banach space of continuous functions mapping the interval into with norm for . Therefore, system (5) with initial conditions (6) has a unique solution by using the standard theory of functional differential equations [61, 62].

Well-posedness of solutions

Proposition 1

All solutions of system (5) with initial conditions (6) are nonnegative and ultimately bounded.

Proof

From the first equation of system (5), we have , then for all . Moreover, the rest of equations of system (5) give us the following: Therefore, for all . Thus, by a recursive argument, we get for all . Hence, the solutions of system (5) with initial conditions (6) satisfy for all . Next, we establish the boundedness of the model’s solutions. The nonnegativity of the model’s solution implies that . To show the ultimate boundedness of , we let Then where . It follows that , where . Since and are nonnegative, then . Further, we let Then we obtain where . It follows that , where . Since , and are nonnegative, then and , where . Furthermore, we let Then Since , then where . It follows that . Since , , , and , then , , and , where and . Finally, from the sixth equation of system (5), we have This implies that , where . □

According to Proposition 1, we can show that the region is positively invariant with respect to system (5).

Steady states analysis

In this section, we calculate all possible steady states of the model and derive the threshold parameters which guarantee the existence of the steady states. Let us define which will be used throughout the paper. Let be any steady state of system (5) satisfying the following equations: We find that system (5) has eight possible steady states.

(i) Infection-free steady state, , where . In this case, the body is free from HTLV and HIV.

(ii) Persistent HIV monoinfection steady state with ineffective immune response, , where and is the basic HIV monoinfection reproduction number for system (5) and is defined as follows: where The parameter determines whether or not a persistent HIV infection can be established. In fact, measures the average number of secondary HIV infected generation caused by an existing free HIV particle due to free-to-cell transmission, while measures the average numbers of secondary HIV infected generation caused by living active HIV-infected cells due to infected-to-cell transmission. The steady state describes the case of persistent HIV monoinfection without immune response.

(iii) Persistent HTLV monoinfection steady state with ineffective immune response, , where and is the basic HTLV monoinfection reproduction number for system (5) and is defined as follows: The parameter decides whether or not a persistent HTLV infection can be established. The steady state describes a persistent HTLV monoinfection without immune response.

We mention that and state the threshold dynamics of infection-free equilibrium and can be calculated by different methods such as (a) the next-generation matrix method of van den Driessche and Watmough [63], (b) local stability of the infection-free equilibrium , and (c) the existence of the chronic HIV and HTLV monoinfection equilibria with inactive immune response. In the present paper, we derive and by method (c).

(iv) Persistent HIV monoinfection steady state with only effective HIV-specific CTL, , where and is the HIV-specific CTL immunity reproduction number in case of HIV monoinfection. The parameter determines whether or not the HIV-specific CTL immune response is effective in the absence of HTLV.

(v) Persistent HTLV monoinfection steady state with only effective HTLV-specific CTL, , where and is the HTLV-specific CTL immunity reproduction number in case of HTLV monoinfection and is stated as follows: The parameter determines whether or not the HTLV-specific CTL immune response is effective in the absence of HIV.

(vi) Persistent HTLV-HIV coinfection steady state with only effective HIV-specific CTL, , where where Here, the parameter is the HTLV infection reproduction number in the presence of HIV infection and determines whether or not HIV-infected patients could be dually infected with HTLV.

(vii) Persistent HTLV-HIV coinfection steady state with only effective HTLV-specific CTL, , where and is the HIV infection reproduction number in the presence of HTLV infection. It is clear that determines whether or not HTLV-infected patients could be dually infected with HIV.

(viii) Persistent HTLV-HIV coinfection steady state with effective HIV-specific CTL and HTLV-specific CTL, , where and The parameter is the competed HIV-specific CTL immunity reproduction number in case of HTLV-HIV coinfection. The parameter is the competed HTLV-specific CTL immunity reproduction number in case of HTLV-HIV coinfection. Clearly, exists when and .

Global stability analysis

In this section, we use the Lyapunov method to show the global asymptotic stability of the model’s steady states. For formation of Lyapunov functionals, we follow the works [64, 65]. Denote , where .

Let a function and be the largest invariant subset of We define a function .

Theorem 1

If and , then is globally asymptotically stable (GAS).

We define a Lyapunov functional as follows: Clearly, for all , and . We calculate along the solutions of model (5) as follows: Summing the terms of Eq. (9), we obtain Using , we obtain Since and , then . Therefore, for all ; moreover, when . The solutions of system (5) converge to . The set includes elements with . Then and the first and fifth equations of system (5) become which give for all t. In addition, we have , and from the third equation of system (5) we have which yields for all t and hence . Applying Lyapunov–LaSalle asymptotic stability (LLAS) theorem [66-68], we get that is GAS. □

The following equalities are needed in the next theorems: Further,

Theorem 2

Let , , and , then is GAS.

Define a functional as follows: Calculate as follows: Summing the terms of Eq. (12), we get The steady state conditions for are given by Then we get Further, we obtain Using the equalities given by (10) in case of , we get Therefore, Eq. (13) becomes Since and , then for all . Moreover, when and . The solutions of system (5) converge to which includes elements that satisfy and i.e. for all . If , then from Eq. (14) we get , , and for all t. Further, for each element of , we have and then . The fifth equation of system (5) becomes which provides for all t, and hence . Therefore, using LLAS theorem we get that is GAS. □

Theorem 3

If , , and , then is GAS.

Define We calculate as follows: By collecting the terms of Eq. (15), we get Using the steady state conditions for we obtain Using the equalities given by (11) in case of , we get Therefore, Eq. (16) becomes Thus, if and , then for all . Moreover, when . The solutions of system (5) converge to which includes elements with . Then we have , and the first equation of system (5) becomes which yields for all t. Moreover, we have and from the third equation of system (5) we get which implies that for all t. Therefore, . Applying LLAS theorem, we get is GAS. □

Theorem 4

Let and , then is GAS.

Define a functional as follows: We calculate as follows: Collecting the terms of Eq. (18), we derive Using the steady state conditions for we get Further, we obtain Using the equalities given by (10) in case of , we get Therefore, Eq. (19) becomes Hence, if , then does not exist since and . In this case Now we want to find the value S̄ such that, for all , we get and . Let us consider This happens when . Clearly, for all , where occurs at and . The solutions of system (5) converge to which includes elements satisfying and i.e. for all . If , then from Eq. (20) we get , , and for all t. Thus, contains elements with , , , and then , . The third and fifth equations of system (5) become which yield and for all t. Therefore, . Applying LLAS theorem, we get is GAS. □

Theorem 5

If and , then is GAS.

Let Calculate as follows: Summing the terms of Eq. (21), we get Using the steady state conditions for we obtain Using the equalities given by (11) in case of , we get Therefore, Eq. (22) becomes Hence, if , then for all , where occurs at , , , and . The trajectories of system (5) converge to which includes elements with , , , and then . The first and fifth equations of system (5) become which imply that and for all t. Moreover, we have , then the third equation of system (5) becomes which yields for all t, and then . Applying LLAS theorem, we get is GAS. □

Theorem 6

If , , and , then is GAS.

Define Calculate as follows: Summing the terms of Eq. (23), we get Using the steady state conditions for we obtain Moreover, we get Using the equalities given by (10) and (11) in case of , we get Therefore, Eq. (24) becomes Hence, if , then for all . One can show that when . The solutions of model (5) tend to which includes elements with , and then . The third equation of system (5) becomes which yields for all t, and hence . Applying LLAS theorem, we get is GAS. □

Theorem 7

If , , and , then is GAS.

Define Calculate as follows: Collecting the terms of Eq. (25), we obtain Using the steady state conditions for we get Moreover, we get Using the equalities given by (10) and (11) in case of , we get Therefore, Eq. (26) becomes Hence, if , then for all . Similar to the previous theorems, one can show that at . The solutions of system (5) reach which contains elements with , , and then . The fifth equation of system (5) becomes which yields for all t, and hence . Applying LLAS theorem, we get is GAS. □

Theorem 8

If and , then is GAS.

Consider Calculate as follows: Summing the terms of Eq. (27), we get Using the steady state conditions for we get Moreover, we get Using the equalities given by (10) and (11) in case of , we get Therefore, Eq. (28) becomes Hence, for all . Similar to the previous theorems, one can show that when . The solutions of system (5) converge to which includes elements with . Then . The third and fifth equations of system (5) become which yield and for all t, and hence . Applying LLAS theorem, we get is GAS. □

Numerical simulations

In this section, we perform numerical simulations to illustrate the results of Theorems 1–8. Moreover, we study the influence of time delays on the dynamical behavior of the system. Let us choose a Dirac delta function as a special form of the kernel as follows: Let , then we get Thus, model (4) reduces to For model (29), the threshold parameters are given by where To solve system (29) numerically, we fix the values of some parameters (see Table 1) and the others will be varied.

Table 1

The data of model (29)

Parameter	Value	Parameter	Value	Parameter	Value
η	10	b	5	ω	0.03
ϱ	0.01	π_1	0.1	ψ	0.003
ϑ_i, i = 1,2,3	Varied	π_2	0.1	ℏ_1	0.2
β	0.7	μ_1	0.2	ℏ_2	0.3
a	0.5	μ_2	0.2	ℏ_3	0.4
φ	0.2	ε	2	ℏ_4	0.5
K	0.9	γ	0.1	ℏ_5	0.6
r^∗	0.008	σ_i, i = 1,2	Varied	ℏ_6	0.9
δ^∗	0.05	λ	0.2	ξ_i, i = 1,2,…,6	Varied

Stability of the steady states

In this subsection, we select the delay parameters as , , , , , . Besides, we choose the following three different initial conditions for system (29):

Initial-1: ,

Initial-2: ,

Initial-3: , where .

Choosing different values of , , , , , and under the above initial conditions leads to the following sets:

Set 1 (Stability of ): , , , , and . For this set of parameters, we have and . Figure 1 illustrates that the trajectories starting different initials converge to the steady state . This supports the global stability result of Theorem 1. Here, a healthy state will be reached where both viruses are absent.

Figure 1

Solutions of system (29) when and

Set 2 (Stability of ): , , , , and . With such a choice we get , , and hence . The steady state exists with . The stability of given in Theorem 2 is shown in Fig. 2. This leads to the case where HIV monoinfection is chronic but with an ineffective CTL immunity.

Figure 2

Solutions of system (29) when , , and

Set 3 (Stability of ): , , , , and . Then we calculate , , and then and . We can see from Fig. 3 that the system’s solutions tend to , which is compatible with Theorem 3. This case means that an HTLV monoinfection is chronic with an ineffective CTL immunity.

Figure 3

Solutions of system (29) when , , and

Set 4 (Stability of ): , , , , and . Then we calculate and . Figure 4 shows that the trajectories starting with different states tend to . Therefore, is GAS, and this supports Theorem 4. Hence, an HIV monoinfection is chronic with effective HIV-specific CTL immunity.

Figure 4

Solutions of system (29) when and

Set 5 (Stability of ): , , , , and . Then we calculate and , and exists with . We observe from Fig. 5 that the system’s trajectories tend to and it is GAS. Here, an HTLV monoinfection is chronic with effective HTLV-specific CTL immunity.

Figure 5

Solutions of system (29) when and

Set 6 (Stability of ): , , , , and . Then we calculate , , and . The numerical solutions of the system drawn in Fig. 6 confirm that exists and is GAS. This case leads to a chronic coinfection with HTLV and HIV where the HIV-specific CTL immunity is effective while the HTLV-specific CTL immunity is ineffective.

Figure 6

Solutions of system (29) when , , and

Set 7 (Stability of ): , , , , and . We compute , , and . According to these values, we obtain that exists. The numerical solutions of our system plotted in Fig. 7 show that is GAS (Theorem 7). This case leads to a chronic coinfection with HTLV and HIV where the HTLV-specific CTL immunity is effective and the HIV-specific CTL immunity is not working.

Figure 7

Solutions of system (29) when , , and

Set 8 (Stability of ): , , , , and . These data give and . According to these values, the steady state exists. Figure 8 illustrates that the solutions of the system initiating with three different states tend to . In this case, a chronic coinfection with HTLV and HIV is reached where both immune responses are well working.

Figure 8

Solutions of system (29) when and

Effect of time delays on the HTLV-HIV dynamics

In this part we vary the delay parameters , , and fix the parameters , , , , and . Since and given by Eqs. (30) and (31) depend on , , then changing the parameters will change the stability of steady states. Let us consider the following situations:

,

, , , , , and ,

, , , , , and ,

, , , , , and .

With these values we solve system (29) under the following initial condition:

Initial-4: , where , .

From Fig. 9 we observe that the presence of time delays can increase the number of uninfected CD4+ T cells and decrease the number of other compartments. Table 2 presents the values and for selected values of , . It is clear that and are decreased when are increased, and thus the stability of can be changed. Let us calculate the critical value of the time delay that changes the stability of . Without loss of generality, we let the parameters and fix , , and write and as functions of and , respectively, as follows: To force the threshold parameters and to satisfy and , we choose , where is the solution of and Therefore, if and , then is GAS. Let us choose the value and and compute , as , . It follows that:

Figure 9

Impact of delay parameters , , on the behavior of solution trajectories of system (29)

Table 2

The variation of and with respect to the delay parameters

Delay parameters	ℜ_1	ℜ_2
ℓ_1=ℓ_2=ℓ_3=ℓ_4=ℓ_5=ℓ_6=0	2.790	3.690
ℓ_1=0.3, ℓ_2=0.4, ℓ_3=0.5, ℓ_4=0.6, ℓ_5=0.7, and ℓ_6=0.8	1.408	1.787
ℓ_1=0.4, ℓ_2=0.5, ℓ_3=0.6, ℓ_4=0.7, ℓ_5=0.8, and ℓ_6=0.9	1.280	1.600
ℓ_1=0.6, ℓ_2=0.7, ℓ_3=0.8, ℓ_4=0.9, ℓ_5=1, and ℓ_6=1.2	1.009	1.283
ℓ_1=1, ℓ_2=1.5, ℓ_3=2, ℓ_4=2.5, ℓ_5=3, and ℓ_6=3.5	0.368	0.173
ℓ_1=2, ℓ_2=3, ℓ_3=4, ℓ_4=5, ℓ_5=6, and ℓ_6=7	0.188	0.008
ℓ_1=3, ℓ_2=4, ℓ_3=5, ℓ_4=6, ℓ_5=7, and ℓ_6=8	0.136	0.003
ℓ_1=4, ℓ_2=6, ℓ_3=8, ℓ_4=9, ℓ_5=10, and ℓ_6=11	0.072	0.1 × 10−3
ℓ_1=6, ℓ_2=7, ℓ_3=9, ℓ_4=10, ℓ_5=11, and ℓ_6=12	0.052	0.3 × 10−4
ℓ_1=10, ℓ_2=11, ℓ_3=12, ℓ_4=13, ℓ_5=14, and ℓ_6=15	0.016	1.2 × 10−6

(i) If and , then , and is GAS.

(ii) If or , then or and will lose its stability.