Accounting for multi-delay effects in an HIV-1 infection model with saturated infection rate, recovery and proliferation of host cells

— Biological models inherently contain delay. Mathematical analysis of a delay-induced model is, however, more difﬁcult compare to its non-delayed counterpart. Difﬁculties multiply if the model contains multiple delays. In this paper, we analyze a realistic HIV-1 infection model in the presence and absence of multiple delays. We consider self-proliferation of CD4 + T cells, nonlinear saturated infection rate and recovery of infected cells due to incomplete reverse transcription in a basic HIV-1 in-host model and incorporate multiple delays to account for successful viral entry and subsequent virus reproduction from the infected cell. Both of delayed and non-delayed system becomes disease-free if the basic reproduction number is less than unity. In the absence of delays, the infected equilibrium is shown to be locally asymptotically stable under some parametric space and unstable otherwise. The system may show unstable oscillatory behaviour in the presence of either delay even when the non-delayed system is stable. The second delay further enhances the instability of the endemic equilibrium which is otherwise stable in the presence of a single delay. Numerical results are shown to be in agreement with the analytical results and reﬂect quite realistic dynamics observed in HIV-1 infected individuals.


I. THE MODEL
Human Immune Deficiency Virus-type 1 (HIV-1 or simply HIV) is a retrovirus which preferably infects CD4 + T lymphocytes and supposed to be the causative agent of AIDS (Acquired Immune Deficiency Syndrome). With the help of T cell receptor, CD4 + T cells recognize the pathogen and interact with other lymphocytes like CD8 + T cells to respond efficiently against the virus so that virus can be either removed or killed. The entry of HIV into a host cell is a complex process [1], [2]. It begins with the adhesion of virus to the host cells and ends with the fusion of its genome into the host cell's cytoplasm. The outer envelope (Env) of HIV contains spikes made up of gp120 and gp41 glycoprotein. The envelope glycoprotein gp120 contains the receptor-binding region and is attached with the receptor of CD4 + T cells. After gp120 binding on CD4 + T molecule, a conformational change occurs in Env followed by dissociation of Env from the viral membrane. The surface of gp120 contains five variable loops which (particularly v3) plays crucial roles in immune evasion and coreceptor binding. Due to dissolution from the viral membrane, amino-terminal hydrophobic domains of gp41 are exposed. This initiates the fusion and HIV enters into CD4 + T cells through six-helix bundle formation [3], [4]. After successful entry of virus in the cell's cytoplasm, information coated in the viral RNA is transcribed into DNA (called provirus) with the help of reverse transcriptase enzyme of virus [5]. Reverse transcriptase, however, is not a good copier [6]. Some observations report that when the virus enters into a CD4 + T cell, viral RNA may not be completely reverse transcribed into DNA [7], [8] and the success depends on the completion time of transcription process [9]. If the cell is activated shortly after entry of the virus into the cell, successful completion of reverse transcription is highly expected [7]. However, if the time is long enough then the partial DNA transcripts may be degraded and the cell may be infection-free [8], [9]. After reverse transcription, provirus moves to the nucleus of CD4 + T cell and integrates with the DNA of the cell with the help of a viral enzyme called integrase and controls over the host cell's machinery [10]. When the infected CD4 + T is activated, the viral genome is transcribed back into RNA, resulting in the production of multiple copies of viral RNA. These RNAs are translated into proteins that require a viral protease to cleave them into active forms. Finally, the mature proteins assemble with the viral RNA to produce new virus particles that bud from the infected cell [11].
Let x(t) be the concentration of susceptible CD4 + T cells in the blood plasma at any time t and y(t) be the concentration of CD4 + T cells in which virus penetration has occurred successfully. We call this later class as infected CD4 + T cells.
Assuming that a proportion of infected CD4 + T cells can be reverted to the uninfected class and v(t) be the concentration of free virus particle in the blood plasma at time t, following mathematical models have been proposed and analyzed [11], [12], [13], [14], [15]: Here f 1 is the demography function of the uninfected CD4 + T cells and f 2 is the incidence function. The parameters d 1 , c and d 2 are the death rate of infected cells, virus reproduction rate and virus clearance rate, respectively. p is the rate at which infected cells revert to the uninfected class due to unsuccessful reverse transcription. In basic HIV models [16], [17], [18], [19], f 1 and f 2 are considered as f 1 = s − µx and f 2 (x, v) = βxv, where s is the constant input rate of CD4 + T cells, µ is its death rate and β is the disease transmission coefficient. It is a well-established fact that CD4 + T cells proliferate upon exposer to HIV antigen through a number of mechanisms, like antigenic stimulation, direct activation by the virus or its products, and homeostatic T cell proliferation as a response to CD4 + T cell depletion [20], [21]. In such a case, density-dependent demography function f 1 (x) = s − µx + rx(1 − x/K) is assumed to be a more realistic one [22], [23]. The parameter r represents the maximum proliferation rate of CD4 + T cells and K is the value where CD4 + T cells proliferation ceases to happen. One of the arguable issues in constructing infection-related model is the representation of the incidence rate. Though most of the epidemic models consider bilinear/mass action incidence, there are many reasons for which this incidence function needs modification. May and Anderson [24] argued that nonlinear incidence should be considered in case of the complex transformation process of parasite infections with indirect life cycles and the case is similar here. Different drawbacks of bilinear/mass action law have also been pointed out, e.g., it does not address the crowding effect when the density of either host cells or virus particle is high [25], [26]. In fact, the incidence becomes unbounded whenever either of x and v becomes large. Many researchers have considered the saturation effect of virus particles only and considered f 2 (x, v) = βxv 1+bv , where β and b are positive constants. Considering virus as predator and CD4 + T cells as its prey, De Boer [27] considered saturation on CD4 + T cells and used f 2 (x, v) = βxv a+x , a > 0 is the half-saturation constant, as the incidence function. Bairagi and Adak [28], [29] considered more general incidence rate of the form f 2 (x, v) = βx n v a n +x n , n ≥ 1. Note that this incidence function does not consider the saturation effect on virus particle. Here we consider a more general incidence function of the form f 2 (x, v) = βxv 1+ax+bv (a, b > 0). It is worth mentioning that this function considers the saturation effect of both CD4 + T cells and virus particle. A similar function has also been used to study the dynamics of other viral infection models [30], [31]. With these assumptions, the model system (1) becomes +py, We consider the initial condition The above model generalizes a large number of HIV-1 infection model. For example, the basic HIV model studied in [16], [17], [18], [19] can be obtained from (2) by setting zero to r, a, b.p.
One can also deduce the basic HIV model with recovery studied in [11], [12] for a = b = r = 0 and the models studied in [14], [15] can be deduced for a = b = 0.
Researchers frequently use delay in the model to take into account various time-consuming events in the biological processes. A single delay has been used in HIV models to represent the incubation delay [25], [26], [32], [33] and virus reproduction delay [34]. An HIV infection model may be more interesting if it considers both delays simultaneously. We here consider a time lag τ 1 between the events a host cell is contacted by a virus particle and the contacted cell becomes infected after various successful biochemical events. We consider another time lag τ 2 which represents the time taken by the virus for completion of its life cycle within the infected cell and subsequent release of the new virus through cell lysis. Local and global stability of the model proposed in (1) has been studied in presence of a single delay τ 1 in [12] with f 1 = s − µx and f 2 = βxv. It has been shown that this delay does not affect the system dynamics and the endemic equilibrium is globally asymptotically stable if the basic reproduction number is greater than 1. Sun and Min [35] analyzed a non-delayed system with f 1 = s − µx and f 2 = βxv x+v and concluded similarly. Zhou et al. [14] has studied the model (1) with f 1 = s − µx + rx(1 − x K ) and f 2 = βxv. It is proved that infection is cleared if the basic reproduction number is less than unity. In the opposite condition, the chronic equilibrium is globally asymptotically stable. There may exist a stable limit cycle around the endemic equilibrium if some additional condition is satisfied. A single delay τ 1 was considered in the model studied by Zhou et al. [14] and instability effect of the delay was observed [15]. Though there exists lots of  [13], [36], [37] which consider multi delays in a basic HIV model because more than one delay increases mathematical complexity significantly. Incorporating two delays τ 1 & τ 2 , our system (2) then reads The model (4) will be analyzed with the initial conditions where τ = max{τ 1 , τ 2 }, τ > 0, and All parameters are assumed to be positive. Parameter descriptions and their values are presented in Table 1. Adak and Bairagi [13] studied the role of immune response in a HIV model when f 2 (x, v) = βx n v 1+ax n , p = 0 in (1). Srivastava et al. [36] considered two delays in a four-dimensional HIV model with f 1 = s−µx and f 2 = βxv, p = 0. They showed that the delay may have both stable and unstable effect on the system dynamics. Pawelek et al. [37] considered the same model and showed the global stability of the interior equilibrium point under some restrictions. They ignored the proliferation character of CD4 + T cells in the presence of antigenic infection, saturation effect of the incidence and recovery of some infected CD4 + T cells due to incomplete reverse transcription. Alshorman et al. [38] considered a two multi-delayed HIV models with f 1 = s − µx and f 2 = βxv, p = 0. One delay is used to represent the time between cell infection and viral production, and the other delay is used to incorporate the time needed for the adaptive immune response to emerge. It is shown that intracellular delay does not change the stability results but immune delay can generate rich dynamics. Multiple delay has also been considered in the fractional-order HIV model [39]. The issue we address here is the role of self proliferation of CD4 + T cells in an in-host HIV-1 ODE model with more general infection rate and recovery of infected cells. The objective is to find how the dynamics change in the presence of multiple delays. In the absence of delays, we show that the disease-free state is locally asymptotically stable if the basic reproduction number is less than unity. In the opposite case, however, the infection always persists. The infected equilibrium is locally asymptotically stable under some parametric space and unstable for some other parametric space. The system may show unstable oscillatory behavior in the presence of either delay even when the nondelayed system is stable. It is further observed that instability is augmented in the presence of multiple delays.
The paper is arranged in the following sequence. In the next section, we analyze the non-delayed system. We show positivity, boundedness and permanence of the system along with the stability results of different equilibrium points. In Section 3, we study the delay-induced system and give local stability results in the presence of single and multiple delays. Numerical results to validate analytical findings are presented in Section 4. The paper ends with a discussion in Section 5.

A. Positivity and boundedness of solutions
Unless stated otherwise, we always assume that self-proliferation rate of healthy CD4 + T cells is greater than its natural death rate, i.e. r > µ. It is also assumed that µK > s, so that CD4 + T cells count decreases if it ever reaches K [40].
PROPOSITION II.1 All solutions of the system (2) are positively invariant and uniformly bounded in , µK > s and r > µ.
We below prove the boundedness of the solutions for all large t. Let V 1 (t) = x(t) + y(t). Differentiating V 1 (t) along the solutions of (2), we

Fig. 1: Sensitivity analysis of system parameters using Partially Ranked Correlation Coefficients (PRCC)
in absence of delays with p-value < 0.001. Lengths of the bar against each parameter measure its sensitivity on the state variables. This diagram shows that c and β are the most sensitive parameters. All parameters have been varied in their reported range mentioned in the Table 1.
, implying that solutions x(t; x(0)) and y(t; y(0)) are uniformly bounded for all large t. Finally, from the third equation of (2), one has lim sup t→+∞ v(t) ≤ cs0 d2d1 . Hence the proposition.

B. Equilibria, stability and permanence
System (2) has two equilibrium points. The disease-free equilibrium E 1 (x 0 , 0, 0) always exists with Then the infected or endemic equilibrium is given by Thus the unique positive root of the quadratic equation is determined as If y * is also positive then E * will exist uniquely. Note that, in view of y * = N (x * ) d1 , we have y * > 0 iff N (x * ) > 0. This last inequality is satisfied if x * < x 0 . Thus, the system (2) possesses a unique interior equilibrium E * whenever x * < x 0 .

C. Basic reproductive number
We determine the basic reproductive number using the next generation matrix method [42]. The Jacobian of (2) at E 1 is given by The infected sub-matrix is One can express J 1 as Then the basic reproduction number R 0 is determined by the spectral radius of the next generation matrix M 1 M −1 2 [42] and evaluated as At E * , we have We, therefore, state the following proposition.

REMARK II.1
The endemic equilibrium Proof: The Jacobian of (2) at E 1 is given by (10). It can be easily verified by a direct calculation that all the eigenvalues of J 0 are negative iff R 0 < 1. Hence the proposition follows.
E. Permanence of the system DEFINITION II.1 System (2) is said to be uniformly persistent if there exists ς > 0 (independent of initial conditions) such that every solution (2) is said to be permanent if there exists a compact region Γ 0 ∈ int Γ L such that every solution of system (2) with initial condition (3) will eventually enter and remain in region Γ 0 .
We further define Metzler matrix and state Perron-Frobenius theorem following Gantmacher [45]. DEFINITION II.3 A square matrix is said to be a Metzler matrix if all its off-diagonal elements are non-negative.
THEOREM II.2 (Perron-Frobenius theorem) Let F be an irreducible Metzler matrix. Then λ F , the eigenvalue of F of largest real part is real and the elements of its associated eigenvector ν F are positive. Moreover, any eigenvector of F with nonnegative elements belongs to span of ν F .
With these results and the concept of persistence in infinite-dimensional system given by Hale and Waltman [44], one can prove the following theorem as in [14]. THEOREM II. 3 Model system (2) is permanent if R 0 > 1, i.e., whenever E * exists.

F. Local stability of the interior equilibrium
After linearizing about E * , one can express system (2) as where and THEOREM II.4 Let R 0 > 1. Then the interior equilibrium E * of (2) is locally asymptotically Proof: The characteristic equation corresponding to M 3 is  Following Routh-Hurwitz conditions, the necessary and sufficient conditions for E * to be locally asymptotically stable are Hence the theorem.

G. Hopf bifurcation analysis of the endemic equilibrium
Suppose E * is locally asymptotically stable. We want to know whether E * will lose its stability when one of the system parameters is smoothly varied. Noting that the virus production rate (c) as one of the most sensitive parameters, we perform this study with respect to c. One can, however, choose another parameter, say β, as the bifurcation parameter.
THEOREM II.5 If c crosses a critical value c * then there exists a Hopf bifurcation around E * when the following conditions hold: Proof: Assume that there exists a critical value c * such that the conditions (i) − (iii) hold. For a Hopf bifurcation to occur at c = c * , the characteristic equation (15) satisfies This equation has three roots where Since l(c * ) = 0 and m(c * ) = B 2 (c * ), we have Solving for l (c * ) from (18) and noting the condition (iii) of Theorem II.5, we have Thus, the transversality condition is satisfied. Therefore, the equilibrium point E * switches its stability (stable to unstable) through Hopf bifurcation at c = c * when the conditions of Theorem II.5 hold. Hence the theorem.

SYSTEM
It is worth mentioning that the system (4) has the same equilibrium points and same expression for the basic reproduction number R 0 as in the case of non-delayed system (2). It is easy to check that for all τ 1 ≥ 0 and τ 2 ≥ 0, E 1 is stable when R 0 < 1 and therefore omitted. We are therefore interested to study the case R 0 > 1, i.e., when the interior equilibrium E * exists, in presence of both delays. Consider the perturbed variables asx = x − x * , y = y − y * ,v = v − v * . Dropping the bars for convenience, system (4) then reads  where the coefficients are as in (14) and with Corresponding linear system is then given by and the corresponding characteristic equation is expressible as (ξ 3 +Āξ 2 +Bξ +C) + (Dξ +Ē)e −ξτ1 One can discuss the following cases in sequel.
Sign of the transversality condition will indicate the direction of Hopf bifurcation. Hence the theorem.

IV. NUMERICAL SIMULATIONS
We considered the default values of the system parameters from their reported range (see Table  1) and simulated the system (2). We first performed a sensitivity analysis using Latin Hypercube Sampling (LHS) method and Partial Ranked Correlation Coefficient (PRCC) technique [58] to determine the most sensitive parameters of the model. As the complexity of the biological models is related to a high degree of uncertainty in estimating the values of various system parameters, uncertainty and sensitivity analysis are essential to interpret the dynamics of biological models. LHS estimates the uncertainty for each input parameter by considering it as a random variable only once in the analysis and by defining probability distribution functions, marginal distributions etc. corresponding to each parameter. PRCC analysis is performed for each input parameter sampled by the LHS scheme and each outcome variable. The bar length against each parameter (see Fig. 1) indicates its sensitivity on the system and the level of significance of the test is given by the p− value, which is obtained to be less than 0.001 implying high accuracy of the LHS-PRCC analysis. It shows that the parameters c and β are the most sensitive parameters out of eleven parameters (excluding delay parameters) of the system. It is noticeable that the values of these two parameters can be changed by antiretroviral drugs used in the treatment of HIV infected patients and therefore it makes sense to be considered as important parameters. Based on this sensitivity analysis, we performed our simulations with respect to these two parameters. All the numerical simulations have been performed using MATLAB R2015a. The time series and bifurcation diagrams have been drawn using solvers ode45 (for the deterministic system) and dde23 (for the delayed system). In Fig. 2, we presented the stability region of different equilibrium points when c and β varied simultaneously. This figure shows that, for any given value of β, the system becomes infection-free when virus reproduction is low and does not cross some critical value of c for which R 0 = 1. Infection, however, persists  Fig. 3.
To demonstrate the effect of single delay τ 2 on the stability of the system, we plotted the stability and instability regions of different equilibrium points in τ 2 − c plane (Fig. 4). It shows the interdependency of this delay with the virus reproduction rate in the dynamic behavior of the system. From this figure, one can determine the critical length of delay τ 2 for some given value of c for which the system will remain in ae state and if it crosses the critical length then the system will be in an unstable state. On the other hand, if one knows the length of delay τ 2 then the value of c can be estimated for obtaining an infection-free system or an infected system with fluctuating or stable populations. To illustrate, we consider c = 800 and compute the critical length of τ 2 as τ * 2 = 0.9877, following Theorem III.2. Therefore, all populations will remain in stable state for τ 2 < 0.9877 and in unstable state if τ 2 > 0.9877 (Fig. 5). In contrast, one can fixed τ 2 (say, τ 2 = 0.6) and then plot a bifurcation diagram to show how the system behavior changes with respect to c (Fig. 6). This figure demonstrates that infection cannot persist if virus production is too low, c < 38. Infection persists in oscillatory state for 38 < c < 663 and in stable state for c > 663. To observe the joint effect of both delays on the system dynamics, following Theorem III.3, we pick up any value of τ 2 , say τ 2 = 0.6, from its stable range 0 < τ 2 < τ * 2 and then compute the critical value of τ 1 as τ * 1 (τ 2 ) = 0.415. Thus, E * will be stable for τ 1 < 0.415 and unstable for τ 1 > 0.415 for the fixed value of τ 2 = 0.6. This behavior of the system is represented by time series solutions (Fig. 7) for τ 1 = 0.37 and τ 1 = 0.45. One can find all such critical values corresponding to the entire stable range of τ 2 ∈ [0, τ * 2 = 0.9877) and obtain the stability region of E * in τ 1 − τ 2 parametric plane (Fig. 8). This figure is important to understand the combined effect of two delays on the system dynamics. It shows that one delay is inversely related to the other. It is worth mentioning that the equilibrium E * loses its stability and population densities fluctuate around it when the value of (τ 1 , τ 2 ) lies above the bifurcation line and it remains stable in the opposite case.

V. DISCUSSION
In this paper, we studied a more general model for the basic HIV-1 in-host infection in presence and absence of multiple delays. We considered self-proliferation of CD4 + T cells, nonlinear saturated infection rate and recovery of infected cells due to incomplete reverse transcription in a basic HIV model. Various cell signalling and biochemical processes are required for a virus to enter into the host cell after attachment to the cell surface and a considerable time is elapsed in accomplishing these activities. We, therefore, incorporated one delay (τ 1 ) in the model system to account for this time. After successful entry, a virus goes through different steps (like uncoating, reverse transcription, circularization, transcription, translation, core particle assembly, final assembly and budding) before producing new virus through cell lysis. To encompass this time, we added another delay (τ 2 ) in the virus growth equation.
The objective was to find how the dynamics of a generalized HIV-1 model change in the presence of multiple delays. Whether infection persists or not it depends on the basic reproduction number of the system which we have determined from the second generation matrix. In particular, we have shown that elimination of infection is possible if the basic reproduction number is less than unity. On the other hand, the infection persists if the basic reproduction number is greater than unity. Local stability of the infected equilibrium of the non-delayed system has been proved under some parametric restrictions.
The novelty of this work lies in the fact that it can reciprocate the experimental observations of many studies [35], [14], [15]. It is reported that virus load in infected individuals may be in steadystate or in oscillatory state [59], [60], [61]. Here we have shown that the endemic equilibrium E * of the non-delayed system may undergo a Hopf bifurcation if the virus production rate, c, crosses a critical value, c * . This implies that the state variables will show stable behaviour if c < c * and oscillatory behaviour if c > c * . Thus, both the steady and fluctuating levels of virus count in an HIV patient may be explained once the Hopf bifurcation point is known. Furthermore, if the non-delayed system is stable, either delay can destabilize the system through Hopf bifurcation if the length of delay crosses some critical value. These observations imply that not only the virus production rate but also the disease transmission delay (τ 1 ) and virus maturation delay (τ 2 ) may be the cause of intra-patient variability of CD4 + T cells and virus load in an HIV patient.
Our simulation results showed four types of dynamic behaviors of the system with respect to the virus reproduction rate. If the virus reproduction is too small then the infection can not establish and the system becomes infection-free. Infection always persists as the virus production crosses the threshold value. The stable coexistence of all population may be possible either for low or for very high virus production rate. coexist in an oscillatory state. Similar complex behavior was observed by Srivastava et al. [36] in a four-dimensional HIV-1 model in the presence of delay. We, however, observed similar behavior in a more generalized three dimensional model in the absence of delay. De Boer [27], Bairagi and Adak [28] and Xu [26] considered two different types of incidence rates with intracellular delay. However, their models could not explain the periodic fluctuations of viral load in HIV-1 patients even in the presence of delay. The model proposed by Zhou et al. [14] shows fluctuations of the endemic equilibrium only in the presence of a delay. Whereas, our model shows oscillations for the variation in a virus production number, even when there is no delay. Moreover, as we have considered a generalized incidence function, our model is able to reproduce the dynamics of various other HIV infection in-host models as a special case of our model. To summarize, we have been able to show that introduction of the single delay makes a system unstable which might be stable in the absence of delay. We have shown an interdependency between virus reproduction delay and virus production rate. Range of virus production for which stable persistence is feasible decreases with the increasing length of virus production delay. A second delay further enhances the instability of the system. The system, which was stable in the presence of a single delay, maybe unstable in the presence of double delay. Thus a basic HIV model that considers nonlinear incidence, self-proliferation of host cells and recovery of infected host cells may show a variety of dynamics even in the absence of delay. Introduction of delay in such a system further increases the complexity in the dynamics. These observations suggest that virus load in infected individuals may be in steady-state or in an oscillatory state. Furthermore, our analysis demonstrates that intra-and inter-patients variabilities of CD4 + T cells and virus particles in case of HIV infec-tion model may be observed through a number of mechanisms and reciprocate the experimental results that there exist considerable fluctuations in the counts of CD4 + T cells and virus particles in the blood of HIV-1 infected individuals.