Stability analysis of general delayed HTLV-I dynamics model with mitosis and CTL immunity

: This paper formulates and analyzes a general delayed mathematical model which describe the within-host dynamics of Human T-cell lymphotropic virus class I (HTLV-I) under the e ff ect Cytotoxic T Lymphocyte (CTL) immunity. The models consist of four components: uninfected CD4 + T cells, latently infected cells, actively infected cells and CTLs. The mitotic division of actively infected cells are modeled. We consider general nonlinear functions for the generation, proliferation and clearance rates for all types of cells. The incidence rate of infection is also modeled by a general nonlinear function. These general functions are assumed to be satisfy some suitable conditions. To account for series of events in the infection process and activation of latently infected cells, we introduce two intracellular distributed-time delays into the models: (i) delay in the formation of latently infected cells, (ii) delay in the activation of latently infected cells. We determine a bounded domain for the system’s solutions. We calculate two threshold numbers, the basic reproductive number R 0 and the CTL immunity stimulation number R 1 . We determine the conditions for the existence and global stability of the equilibrium points. the develop two types of treatments, the first type aims to extend the intracellular delay periods, while the second type aims to activate and stimulate the CTL immune response.


Introduction
Human T-cell lymphotropic virus class I (HTLV-I) is one of the most dangerous viruses that infect the human body and causes inflammatory and malignant diseases. It was estimated that about 10 to 25 million humans infected with HTLV-I worldwide [1]. HTLV-I can be transferred via bodily fluids including blood and semen. Further, HTLV-I can be transferred from mother to child via breast-feeding [2]. The high transmission of HTLV-I was reported in Caribbean Islands, Central Africa, Northeast South America and Southwestern Japan [3]. HTLV-I can progress to serious diseases, adult T-cell leukemia (ATL) and HTLV-I-associated myelopathy/tropical spastic paraparesis (HAM/TSP) [4]. HTLV-I is a single-stranded RNA virus which infects the most important component of the immune system, CD4 + T cells [5]. There are two modes of within-host transmission, vertical or mitotic transmission of actively HTLV-I infected cells and horizontal by infectious transmission through direct cell-to-cell touch [6].
Mathematical modeling of HTLV-I was established since 1990 to understand the pathogenesis of the HTLV-I infection inside the host cell that has a critical ramifications for the movement of remedial measures and for the acknowledgment of hazard elements of the advancement of HAM/TSP. A lot of mathematical models were planned by numerous scientists to epitomize the intuitive dynamics withinhost among HTLV-I infection to depict the pathogenesis of HTLV-I related infection [7][8][9][10][11][12].
Stilianakis and Seydel [7] constructed a four-dimensional system of ODEs which describes the within-host HTLV-I dynamics as: (1.1)

2)
A(t) = πL(t) − (δ * + κ)A(t), (1.3) where U, L, A and Q (cells mm −3 ) are the concentrations of healthy or uninfected CD4 + T cells, latently CD4 + T infected cells, actively CD4 + T infected cells and leukemia cells (ATL cells), respectively. The variable t denotes to the time (days). Uninfected cells are created at a constant rate ν. The bilinear incidence term βUA accounts for the rate of infectious via horizontal transmission appears by cell-tocell touch between uninfected and actively infected cells. πL is the activation rate of latently infected cells, and κA is the conversion rate of actively infected cells to ATL cells. Q max is the maximal growth of ATL cells and λ is the maximum proliferation rate of ATL cells. The death rates of the uninfected cells, latently infected cells, actively infected cells and ATL cells are indicated byνU, µL, δ * A and ℓQ, respectively.
As mentioned before HTLV-I has two modes of transmission, horizontal and vertical. Actively infected cells propagate faster than uninfected and latently infected cells. This causes an increasing in the proviral load. As a consequence, mitotic transmission plays an instrumental role in the persistence of HTLV-I infection [3]. Li and Lim [10] developed an HTLV-I infection model by considering both horizontal and vertical transmissions as: (1.5)

6)
where ϵ ∈ (0, 1) is the probability of new HTLV-I infections via horizontal route become latently infected. The vertical transmission is represented in proliferating the actively infected cells at rate ε * A 1 − U+L U max , where U max is the CD4 + T cells carrying capacity. The actively infected cells become latent at rate ωε * A 1 − U+L U max , where ω ∈ (0, 1). It was shown in [13] that U + L < U max . Therefore, in [14], the logistic term ε * A 1 − U+L U max was replaced by an exponential growth term ε * A. Cytotoxic T lymphocyte (CTL) is identified as the tremendous element of human immunity in killing viral infections. It represses viral reproduction and kills the cells which can be tainted by infections. It was mentioned in [13,15] that the CTLs assume a strong part in controlling the HTLV-I related diseases. CTLs can realize and kill the actively infected CD4 + T cells, further, they able to lessen the proviral load. The effect of CTL immunity was included in several models for HTLV-I infection (e.g., [4,[16][17][18][19][20][21][22][23][24][25]).
To incorporate both CTL immunity and mitosis into the HTLV-I dynamics model, Lim and Maini [14] proposed the following system: (1.8)
The authors investigated the global stability for chronic infection equilibrium point via two approaches Lyapunov and geometric.
In the above models, it was assumed that the uninfected CD4 + T cells become latently infected instantaneously. Moreover, the transition from latently infected to actively infected is also instantaneously. However, there exist delays in the infection process and activation of latently infected cells. Therefore, ignoring these time delays is biologically unrealistic. In addition, these time intracellular delays have an important effect on the stability of system [29,30]. HTLV-I infection models have been incorporated with intracellular delay [23,24,[31][32][33] or immune response delay [16,20,21,34,35]. Both intracellular and immune response delays were considered in [25,29,36]. Katri and Ruan [31] modified models (1.1)-(1.4) by including a discrete-time intracellular delay τ 1 which describes the time between the initial infection of uninfected CD4 + T cells to become latently CD4 + T infected cells. Equation (1.2) was replaced bẏ The delayed HTLV-I infection models presented in the literature included either discrete-time delay [16,20,21,24,25,[32][33][34][35][36], or distributed-time delay [23]. All these models neglected the latently infected cells.
In models (1.1)-(1.4), (1.5)-(1.7) and (1.8)-(1.11) and most of the HTLV-I models introduced in the literature, the death rate of the cells is given by linear function of their concentration. However, the death rate of the cells is generally not known.
Our aim of the present paper is to formulate and analyze a general delayed HTLV-I mathematical model under the following generalizations: G1: We assume that the intrinsic growth rate of the uninfected CD4 + T cells is given by a general function Φ(U).
G2: We consider a general incidence rate in the form Λ(U, A). G3: We assume that the death rates of the latently infected cells, actively infected cells and CTLs are given by general nonlinear functions µΨ 1 (L), δ * Ψ 2 (A) and γΨ 3 (C), respectively.
G4: We assume that the transition rate from latently infected phase to actively infected phase is given by a general nonlinear function in the form πΨ 1 (L).
G6: We assume that the replenishment rate of CTLs and the elimination of actively infected cells by CTLs are given by ρΨ 2 (A) Ψ 3 (C) and αΨ 2 (A) Ψ 3 (C), respectively. G7: We include distributed intracellular delay which describes the cell mutates from uninfected to latent.

G8:
We include distributed delay in the activation of latently infected.

HTLV-I dynamics model with delay-distributed
With the above generalizations G1-G8 we propose the following HTLV-I dynamics model with two types of distributed-time delays:

4)
Here, ϖ is a random variable generated from probability distribution functions f i (ϖ), i = 1, 2 over the time interval [0, κ i ], i = 1, 2 where κ i is the upper limit of the delay period. Here f 1 (ϖ)e −ξ 1 ϖ represents the probability of surviving the cells from time t − ϖ to time t. The factor f 2 (ϖ)e −ξ 2 ϖ represents the probability that latently infected cells is survived ϖ time units to become actively infected cells.
Here, Φ, Λ, Ψ i , i = 1, 2, 3 are continuously differentiable functions, in addition they satisfy a set of conditions: for all A > 0. Let κ = max {κ 1 , κ 2 } and assume the initial conditions of systems (2.1)-(2.4) having the following form: We note that the systems (2.1)-(2.4) with initial conditions (2.5) by the fundamental theory of functional differential equations has a unique solution [40].

Remark 1.
If the time delay is not considered, ϖ = 0, all proliferated actively infected by mitosis become latent, ω = 1, the intrinsic growth rate of the uninfected CD4 + T cells is given by linear function, Φ(U) = ν −νU, the infection rate is given by bilinear incidence, Λ(U, A) = βUA, the death rates of the compartments are given by linear functions Ψ i (x) = x, i = 1, 2, 3, and the CTL immune response is given by linear form, ρA, then systems (2.1)-(2.4) will lead to the ODEs systems (1.8)-(1.11). In addition, if the CTL immune response is given by a nonlinear form, ρAC, then systems (2.1)-(2.4) will reduce to the model presented in [1]. Thus, systems (2.1)-(2.4) can be considered a generalization of the systems presented in [14,1]. We mention that functions Φ (U), Λ(U, A) and Ψ i , i = 1, 2, 3, include and generalize many forms presented in the literature. The general distribution functions f i , i = 1, 2 considered in our work allows us to include several special forms of intracellular delays existing in the literature (see [24,[31][32][33]). Finally, models (2.1)-(2.4) can be considered as a generalization of the model presented in [23] which neglected the latently infected cells.

Properties of solutions
Proposition 1. Suppose that Conditions C1-C3 are valid. Then the solutions of (2.6)-(2.9) are nonnegative and ultimately bounded.

Equilibrium points and threshold parameters
Here, we look at the model's equilibrium points and deduce the conditions that allow them to exist. Lemma 1. Suppose that C1-C4 are valid, then there exist two positive threshold numbers R 0 and R 1 with R 0 > R 1 such that: (a) if R 0 ≤ 1, then there exists a single equilibrium point EP 0 = (U 0 , 0, 0, 0). (b) if R 1 ≤ 1 < R 0 , then there exist two equilibrium points EP 0 and EP 1 = (U 1 , L 1 , A 1 , 0). (c) if R 1 > 1, then there exist three equilibrium points EP 0 , EP 1 and EP 2 = (U 2 , L 2 , A 2 , C 2 ).

It is clear from Conditions
∂Λ(U 0 ,0) ∂A > 1, then Γ ′ 1 (U 0 ) < 0 and there exists U 1 ∈ (0, U 0 ) such that Γ 1 (U 1 ) = 0. It follows from Eq (2.14) that, L 1 = Ξ(U 1 ) > 0 and A 1 = Π(U 1 ) > 0. This shows that, an infected equilibrium point without CTL immunity which defines the basic infection reproductive number and decides when the HTLV-I infection will be chronic. The second possibility of Eq (2.13) is Ψ 3 (C) 0 and Ψ 2 (A 2 ) = γ ρ . This yields that We insert the value of A 2 in Eq (2.10) and let Based on Conditions C1 and C2, Γ 2 is a strictly decreasing function of U. Clearly, From Eq (2.12) we get which represents the CTL immunity stimulation number and determines when CTL immunity can be activated. Therefore, C 2 can be written as: Thus, there exists an infected equilibrium point with CTL immunity Conditions C2 and C4 imply that

Global stability analysis
Stability analysis of equilibrium points is at the heart of dynamical system analysis. Stable solutions can only be observed experimentally. Thus, we demonstrate the global asymptotic stability of all equilibrium points in this section by establishing appropriate Lyapunov functions [44][45][46] and applying Lyapunov-LaSalle asymptotic stability theorem (L-LAST) [47][48][49]. We define the function H(x) ≥ 0 for all x > 0 as: We will use the follow equations: Let T j (U, L, A, C) be a Lyapunov function candidate and define a set and Υ ′ j be the largest invariant subset of Υ j . Theorem 1. If R 0 ≤ 1 and Conditions C1-C4 are satisfied, then the equilibrium point EP 0 is globally asymptotically stable (GAS).

Special case of the distributed-time delay
In this section, we take a particular form of the probability distributed functions as: where ∆(.) is the Dirac delta function and τ i , i = 1, 2 are constants. When κ i → ∞, i = 1, 2, we have Hence, models (2.6)-(2.9) will reduce to the following model with discrete-time delays:
(i) If R 0 ≤ 1, then the equilibrium point EP 0 is GAS. (ii) If R 1 ≤ 1 < R 0 , then the infected equilibrium point without CTL immunity EP 1 is GAS. (iii) If R 1 > 1, then the infected equilibrium point with CTL immunity EP 2 is GAS. Proof. (i) We consider the Lyapunov function T 0 as: Following the proof of Theorem 1 we obtain dT 0 dt as: (ii) We define the Lyapunov function T 1 as: Then, we follow the proof of Theorem 2 to get (iii) We define the Lyapunov function T 2 as: Following the proof of Theorem 3 we obtain dT 2 dt as: (3.5) Similar to the proofs of Theorems 1-3 one can complete the proof.

Numerical simulations
To illustrate and validate our theoretical results we present some numerical simulations for systems (3.1)-(3.4) with the following specific shapes of the general functions: • Logistic function for the intrinsic growth rate of uninfected cells • Crowley-Martin incidence rate • Linear death rate of the cells Ψ i (s) = s, i = 1, 2, 3.

2)
3) It was mentioned in [16] that the term ς −ν is not necessary to be positive, therefore we assume that ς <ν. We see that  [10,13,14,26] d dU Also, we can see that Ψ i (s) > 0, for all s > 0 and Ψ i (0) = 0, i = 1, 2, 3. Moreover, Ψ ′ i (s) = 1 > 0, i = 1, 2, 3 for all s ≥ 0. Thus, Condition C3 is satisfied. Finally, we have Thus, Condition C4 is satisfied as well. We can also see that our example satisfied the global stability results (Theorems 1-3) with, R 0 and R 1 given by where U 2 is the positive root of the equation We present some numerical results for models (4.1)-(4.4) to address the stability of equilibrium points. Further, we address the effect of time delays on the dynamical behavior of the system. All the simulations are done using the solver dde23 MATLAB.

Stability of equilibrium points
In this subsection, we illustrate our global stability results given in Theorems 1-3 by showing that from any selected initial conditions, the solution of the HTLV-I infection model will tend to one of the three equilibrium points. Let us select three initial conditions as: We fix the time delays τ 1 and τ 2 to become τ 1 = τ 2 = 0.1. We use the values of some parameters of models (4.1)-(4.4) given in Table 1 and establish the effect of the other parameters, β and ρ on the stability of the three equilibrium points. We investigate three cases as follows: Case 1. (R 0 ≤ 1): Choosing β = 0.0001 and ρ = 0.0001 which gives R 0 = 0.1276 < 1 and R 1 = 0.0332 < 1. Lemma 1 and Theorem 1 prove that the system has a single equilibrium point, EP 0 and it is GAS. Figure 1 shows that for all the three initials IC1-IC3, the concentrations of latently infected cells, actively infected cells and CTLs are tending to zero. In the mean time, the concentration of the uninfected cells tends to its healthy value U 0 = 1288.15. This case means that there is no HTLV-I infection in the body.
Case 2. (R 1 ≤ 1 < R 0 ): We take β = 0.004 and ρ = 0.0001. So, we get R 0 = 5.1054 > 1 and R 1 = 0.0936 < 1. Based on Lemma 1 and Theorem 2, the infected equilibrium point without CTL immunity EP 1 exists. In such case, EP 1 is GAS. Figure 2 shows that the numerical simulations illustrate our theoretical result given in Theorem 2. The system's solutions converge to the equilibrium point EP 1 = (144.65, 343.80, 32.29, 0) for all the three initials IC1-IC3. In this situation, the HTLV-I-patient has non-active CTL immunity.  Figure 3 displays the numerical solutions of the system converge to EP 2 for all the three initials IC1-IC3. The results support the theoretical results presented in Theorem 3. This case corresponds to an HTLV-I chronic infection with active CTL immunity.

Impact of the CTL immunity on the dynamics of HTLV-I
We observe that the values of the parameter α, ρ and γ have no effect on the value of R 0 . It means that, CTL immunity has no effect on the stability of the infection-free equilibrium point EP 0 . Now we address the effect of the proliferation rate constant of CTL ρ on the HTLV-I dynamics. Using the values of the parameters in Table 1 to solve the model (4.1)-(4.4) with initial condition IC2. We choose β = 0.004, τ 1 = τ 2 = 0.1. We can see from Table 2 that R 1 can be increased by increasing the value of ρ. We can see from Figure 4 when ρ is increased, the concentrations of uninfected cells and CTLs are increased, while the concentration of actively infected cells is decreased. Therefore, CTL immunity does not play the role in deleting the virus, but it has an important role in controlling the HTLV-I infection.
4.3. Impact of time delays τ 1 and τ 2 on the dynamics of HTLV-I In this part, we study the effect of vary the time delays τ 1 and τ 2 on the HTLV-1 dynamics. Using the values of the parameters in Table 1 and considering β = 0.004 and ρ = 0.001 to solve the models Consider τ = τ 1 = τ 2 , then R 0 and R 1 can be given as functions of τ as follows: . .
We note that R 0 and R 1 are decreasing functions of τ. Let τ CR 0 and τ CR 1 be such that R 0 τ CR 0 = 1 and R 1 τ CR 1 = 1. Using the values of the parameters in Table 1 and considering β = 0.004 and ρ = 0.001, we obtain τ CR 0 = 0.907187 and τ CR 1 = 0.142356. Therefore, we have the following cases: (a) if τ ≥ 0.907187, then R 0 (τ) ≤ 1 and EP 0 is GAS; (b) if 0.142356 ≤ τ < 0.907187, then R 1 (τ) ≤ 1 < R 0 (τ) and EP 1 is GAS; (c) if τ < 0.142356, then R 1 (τ) > 1 and EP 2 is GAS. It is clear from Figure 5 that as τ increases, the concentration of uninfected CD4 + T cells is increased and the other concentrations are decreased. In Table 3 we calculate the values of R 0 and R 1 as functions of τ. We note that, when τ is increased, the values of both R 0 and R 1 are decreased. We can see from the above discussion that increasing time delays values can have similar effect as antiviral treatments.

Conclusions and discussion
In this work, we proposed and examined a general HTLV-I infection models with CTL immunity and mitotic transmission of actively infected cells. We introduce two intracellular distributed-time delays into the models: (i) delay in the formation of latently infected cells, (ii) delay in the activation of latently infected cells. We constructed the model with more general nonlinear functions for the generation/stimulation and clearance rates of all compartments as well as the incidence rate of infection. We established some conditions on such general functions and determined two threshold numbers R 0 and R 1 which control the existence and global stability of the three equilibrium points. We formulated Lyapunov functions and used L-LAST to establish the global asymptotic stability of the equilibrium points. We proved that • If R 0 ≤ 1, then the infection-free equilibrium point EP 0 is GAS. This result suggests that when R 0 ≤ 1, the HTLV-I infection is predicted to die out regardless of the initial conditions. From a control viewpoint, making R 0 ≤ 1 will be an ideal way but HTLV-I infection is lifelong, and the viruses are rarely cleared.
• If R 1 ≤ 1 < R 0 , then the infected equilibrium point without CTL immunity EP 1 is GAS. This result establishes that when R 1 ≤ 1 < R 0 , the HTLV-I infection with inactive CTL immunity is always established regardless of the initial conditions.
• If R 1 > 1, then the infected equilibrium point with CTL immunity EP 2 is GAS. This result illustrates that when R 1 > 1, the HTLV-I infection with active CTL immunity is always established regardless of the initial conditions. This implies that the infection becomes chronic and the CTL immune response is persistent. In cases R 1 ≤ 1 < R 0 and R 1 > 1, the individuals have a high risk of developing HAM/TSP and can be characterized as the HAM/TSP patients.    We studied the global stability of a general HTLV-I model with discrete-time delays as a special case of the model with distributed-time delays when the probability distributed function is given by Dirac delta function. To support the obtained theoretical results, we conducted numerical simulations for the model with discrete-time delays, logistic form for the intrinsic growth rate of uninfected cells, Crowley-Martin incidence rate and linear death rate of the cells. We discussed the impacts of the proliferation of CTLs on the HTLV-I dynamics. We noted that the CTL immunity does not play the role in deleting the virus from the body, but it has an important role in controlling and suppressing the HTLV-I infection. We showed that the parameter ρ has a significant effect of the dynamical behavior of the HTLV-I. The effect of the time delays on the HTLV-I dynamics was discussed. We established that, when all other parameters are fixed and delays are sufficiently large, R 0 becomes less than one, which makes the infection-free equilibrium point EP 0 globally asymptotically stable. From a biological viewpoint, time delays play positive roles in the HTLV-I infection process in order to eliminate the virus. Sufficiently large time delays makes the HTLV-I development slower, and the HTLV-I is controlled and disappears. This gives us some suggestions on new drugs to prolong the time of formation of latent infected cells, or the time of activation of latent infected cells.  We mention that models (1.8)-(1.11) presented by Lim and Maini [14] considered a linear proliferation rate of the CTLs, ρA. In comparison with our proposed model, models (1.8)-(1.11) admits only two equilibrium points, infection-free equilibrium point EP 0 = ( ν ν , 0, 0, 0) and infected equilibrium point EP 1 = (Ũ 1 ,L 1 ,Ã 1 ,C 1 ). The existence and stability of EP 0 and EP 1 are determined by the basic reproductive numberR 0 = π δ * (π + µ) β ν ν + ε * .  We note that,R 0 does not depends on the proliferation rate constant of the CTLs, ρ. As a result, ρ does not play any role on the existence and stability of both equilibria. In this case, the effect of the CTL immune response on the dynamical behavior of the HTLV-I does not appear explicitly. In our model, by considering a nonlinear proliferation rate of the CTLs, the role of the CTLs in controlling the HTLV-I infection is well presented. In [1], an HTLV-I infection model was presented with a nonlinear proliferation rate of the CTLs, ρAC, however, the intracellular time delay was neglected. In addition, the infection rate was given by bilinear incidence which may not be completely characterize the incidence between uninfected cells and actively infected cells [37]. Overall, our proposed model helps to better understand the dynamical behavior of HTLV-I infection with CTL immunity and time delays.  The model presented in this paper can be extended in many directions such as (i) including the reaction diffusion [55], (ii) considering the stochastic interactions [56], (iii) modeling the coinfection with HTLV-I and other types of viruses.