www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Accounting for multi-delay effects in an HIV-1 infection model with saturated infection rate, recovery and proliferation of host cells Debadatta Adak1, Nandadulal Bairagi2, Robert Hakl3 1Department of Applied Mathematics Maharaja Bir Bikram University, Agartala, India dev.adak.math95325@gmail.com 2Centre for Mathematical Biology and Ecology Department of Mathematics, Jadavpur University, Kolkata, India nbairagi@yahoo.co.in 3Institute of Mathematics Czech Academy of Sciences, Branch in Brno, Brno, Czech Republic hakl@drs.ipm.cz Received: 11 November 2019, accepted: 29 December 2020, published: 31 December 2020 Abstract— Biological models inherently contain delay. Mathematical analysis of a delay-induced model is, however, more difficult compare to its non-delayed counterpart. Difficulties 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 con- sider 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 reflect quite realistic dynamics observed in HIV-1 infected individuals. Keywords-HIV model, saturated incidence, self- proliferation, recovery, multiple delays, stability, bifurcation Copyright: © 2020 Adak et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Debadatta Adak, Nandadulal Bairagi, Robert Hakl, Accounting for multi-delay effects in an HIV-1 infection model with saturated infection rate, recovery and proliferation of host cells, Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 1 of 20 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... 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 im- mune evasion and coreceptor binding. Due to dis- solution 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 cy- toplasm, 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]: ẋ = f1(x) −f2(x,v) + py, ẏ = f2(x,v) −d1y −py, (1) v̇ = cy −d2v, Here f1 is the demography function of the un- infected CD4+T cells and f2 is the incidence function. The parameters d1,c and d2 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], f1 and f2 are considered as f1 = s − µx and f2(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 pro- liferation as a response to CD4+T cell depletion [20], [21]. In such a case, density-dependent de- mography function f1(x) = s − µx + rx(1 − x/K) is assumed to be a more realistic one [22], Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 2 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... [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 f2(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 f2(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 f2(x,v) = βxnv an+xn ,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 f2(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 dx dt = s−µx+rx ( 1− x K ) − βxv 1+ax+bv +py, dy dt = βxv 1 + ax + bv −d1y −py, (2) dv dt = cy −d2v. We consider the initial condition x(0) > 0, y(0) ≥ 0, v(0) ≥ 0. (3) 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 de- duced 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 in- cubation 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 be- tween the events a host cell is contacted by a virus particle and the contacted cell becomes infected af- ter various successful biochemical events. We con- sider 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 f1 = s−µx and f2 = β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 f1 = s − µx and f2 = βxvx+v and concluded similarly. Zhou et al. [14] has studied the model (1) with f1 = s−µx + rx(1 − xK ) and f2 = β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 Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 3 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... TABLE I: Parameter descriptions with their values and sources. Parameter Description Reported range References Default value s Constant inpute 0-10 cells mm−3 [50], [51] 10 rate of CD4+T cells µ Death rate of 0.0014-0.03 day−1 [54], [53] 0.01 susceptible CD4+T cells r Proliferation rate 0.03-3 day−1 [50], [52] 3 of CD4+T cells K CD4+T cell density 600-1600 cells mm−3 [52], [23] 1500 where proliferation stops β Disease transmission 0.00025-0.5 [50], [52] 0.0006 coefficient mm−3 virion day−1 a Positive constant ... ... 0.007 b Positive constant ... ... 0.00001 p Recovery rate of 0.01-1.3 day−1 [55], [56] 0.01 infected CD4+T cells d1 Removal rate of 0.16-1 day−1 [57], [12] 1 infected CD4+T cells c Virus production rate constant 26-4000 virion day−1 [56], [57] 800 d2 Removal rate of virus 2-5 day −1 [52], [57], [53] 2.5 τ1 Virus transmission delay 1 day [51] variable τ2 Virus reproduction delay 2 days variable literature with a single delay, there are very few works [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 dx(t) dt =s−µx(t) + rx(t)(1 − x(t) K ) − βx(t)v(t) 1 + ax(t) + bv(t) + py(t), dy(t) dt = βx(t− τ1)v(t− τ1) 1 + ax(t− τ1) + bv(t− τ1) −d1y(t) −py(t), dv(t) dt =cy(t− τ2) −d2v(t). (4) The model (4) will be analyzed with the initial conditions x(θ) = %1(θ),y(θ) = %2(θ),v(θ) = %3(θ), θ ∈ [−τ, 0], (5) where τ = max{τ1,τ2}, τ > 0, and % = (%1,%2,%3) ∈ R3+ with %1(θ) ≥ 0, %2(θ) ≥ 0, %3(θ) ≥ 0. 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 f2(x,v) = βxnv 1+axn , p = 0 in (1). Srivastava et al. [36] considered two delays in a four-dimensional HIV model with f1 = s−µx and f2 = β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, satura- tion 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 f1 = s−µx and f2 = βxv,p = 0. One delay is used to repre- sent the time between cell infection and viral pro- duction, and the other delay is used to incorporate the time needed for the adaptive immune response Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 4 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... 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 non- delayed 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 per- manence 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. II. ANALYSIS OF THE NON-DELAYED SYSTEM 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 ΓL when t is large, where ΓL = { (x,y,v) ∈ R3+ | 0 s and r > µ. Proof: From the second and third equations of the system (2), it follows that if y(0) = 0 and v(0) = 0 then y(t) = 0 and v(t) = 0 for t ∈ R+. We therefore assume that y(0) + v(0) > 0. (7) Then from the second equation of (2), we obtain y(t) = e−(d1+p)ty(0) + e−(d1+p)t ∫ t 0 e(d1+p)s βx(s)v(s) 1 + ax(s) + bv(s) ds for t > 0, (8) which together with (7) implies that there exists T > 0 such that y(t) > 0 for t ∈ (0,T). Conse- quently, the third equation of (2) yields v(t) > 0 for t ∈ (0,T). We will show that x(t) > 0 for t ∈ (0,T). (9) Indeed, assume that x has a zero in (0,T) and let t0 ∈ (0,T) be such that x(t0) = 0, x(t) > 0 for t ∈ [0, t0). Then, obviously, ẋ(t0) ≤ 0. On the other hand, the first equation in (2) yields ẋ(t0) = s + py(t0) > 0, a contradiction. Thus (9) holds. Now let T0 > 0 be such that y(T0) = 0, y(t) > 0 for t ∈ (0,T0). Then, from the above-proven, we have x(t) > 0, v(t) > 0 for t ∈ (0,T0) and so, from (8) we get y(T0) > 0, a contradiction. Consequently, y(t) > 0 for t > 0 and therefore also x(t) > 0 and v(t) > 0 for t > 0. We below prove the boundedness of the solu- tions for all large t. Let V1(t) = x(t) + y(t). Differentiating V1(t) along the solutions of (2), we Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 5 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... Parameters s r K a b p c P R C C -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 x y v µ β d 1 d 2 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. have V̇1(t) = ẋ(t) + ẏ(t) = s + [(r −µ) + d1]x(t) − rx2(t) K −d1[x(t) + y(t)] ≤ s0 −d1V1(t), where s0 = s + K[(r −µ) + d1]2 4r . Therefore, lim sup t→+∞ V1(t) ≤ s0d1 , implying that so- lutions 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 E1(x0, 0, 0) always exists with x0 = K2r [(r − µ) + √ (r −µ)2 + 4sr K ]. For simplicity, we write N(x) = s−µx+rx ( 1− x K ) . Then the infected or endemic equilibrium is given by E∗ = (x∗,y∗,v∗), where v∗ = cy∗ d2 , y∗ = N(x∗) d1 and x∗ is the positive root of A1x 2 + A2x + A3 = 0, where A1 = rbc d1d2K (> 0), A2 = βc d2(d1 + p) − (r −µ)bc d1d2 −a, A3 = − ( sbc d1d2 + 1 ) (< 0). Thus the unique positive root of the quadratic equation is determined as x∗ = −A2 + √ A22 − 4A1A3 2A1 > 0. Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 6 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... c0 200 400 600 β 0.00025 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 Unstable E * Stable E * Stable E * Stable E 1 Fig. 2: Stability regions of the disease-free equilibrium E1 and infected equilibrium E∗ in c−β plane. Here s = 10, µ = 0.01, r = 3, K = 1500, a = 0.007, b = 0.00001, p = 0.01, d1 = 1, d2 = 2.5, τ1 = τ2 = 0. 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∗ < x0. Thus, the system (2) possesses a unique interior equilibrium E∗ whenever x∗ < x0. C. Basic reproductive number We determine the basic reproductive number using the next generation matrix method [42]. The Jacobian of (2) at E1 is given by J0 =   r −µ− 2rx0 K p − βx0 1+ax0 0 −(d1 + p) βx01+ax0 0 c −d2   . (10) The infected sub-matrix is J1 = ( −(d1 + p) βx01+ax0 c −d2 ) . One can express J1 as J1 = ( 0 βx0 1+ax0 0 0 ) − ( d1 + p 0 −c d2 ) = M1 −M2. Then the basic reproduction number R0 is deter- mined by the spectral radius of the next generation matrix M1M −1 2 [42] and evaluated as R0 = βcx0 d2(1 + ax0)(d1 + p) . (11) At E∗, we have βcx∗v∗ 1 + ax∗ + bv∗ = (d1 + p)y ∗ or, equivalently, βcx∗ d2(d1+p) [ 1+ax∗+ bc d1d2 ( s−µx∗+rx∗(1−x∗ K ) )] = 1. Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 7 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... Define M(x) = βcx d2(d1+p) [ 1+ax+ bc d1d2 ( s−µx+rx(1− x K ) )], P(x) = 1 M(x) , ∀x > 0, (12) so that M(0) = 0, M(x∗) = 1, M(x0) = R0 and P ′(x) = −d2(d1+p) βc ( bcs d1d2 1 x2 + 1 x2 + rbc d1d2K ) < 0 for all x > 0. Thus, P(x) is monotonically decreasing for all x > 0, which in turn implies that M(x) is monotonically increasing for all x > 0. Therefore, x0 > x∗ ⇒ M(x0) > M(x∗) ⇒ R0 > 1. We, therefore, state the following proposition. REMARK II.1 The endemic equilibrium E∗(x∗,y∗,v∗) of system (2) exists uniquely if x∗ < x0, or equivalently R0 > 1. D. Local stability of disease-free equilibrium THEOREM II.1 The disease-free equilibrium E1 is locally asymptotically stable if R0 < 1. Proof: The Jacobian of (2) at E1 is given by (10). It can be easily verified by a direct calcula- tion that all the eigenvalues of J0 are negative iff R0 < 1. Hence the proposition follows. E. Permanence of the system DEFINITION II.1 System (2) is said to be uni- formly persistent if there exists ς > 0 (indepen- dent of initial conditions) such that every solution (x(t),y(t),v(t)) with initial condition (3) satisfies lim inf t→+∞ x(t)≥ς, lim inf t→+∞ y(t)≥ς, lim inf t→+∞ v(t)≥ς. DEFINITION II.2 System (2) is said to be perma- nent 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 non- negative 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 theo- rem as in [14]. THEOREM II.3 Model system (2) is permanent if R0 > 1, i.e., whenever E∗ exists. F. Local stability of the interior equilibrium After linearizing about E∗, one can express system (2) as dX dt = M3X(t), (13) where X(t) = [x(t), y(t), v(t)]T , M3 =   m11 m12 m13n21 m22 n23 0 l32 m33   , and  m11 = r−µ− 2rx∗ K − βv∗(1+bv∗) (1+ax∗+bv∗)2 , m12 = p, m13 = − βx∗(1 + ax∗) (1 + ax∗ + bv∗)2 , n21 = βv∗(1 + bv∗) (1 + ax∗ + bv∗)2 , m22 = −(d1 + p) n23 = βx∗(1 + ax∗) (1 + ax∗ + bv∗)2 , l32 = c, m33 = −d2. (14) THEOREM II.4 Let R0 > 1. Then the interior equilibrium E∗ of (2) is locally asymptotically stable if B1 > 0, B3 > 0 and B1B2 − B3 > 0, where Bi (i = 1, 2, 3) are defined bellow. Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 8 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... 0 500 1000 1,000 1,200 1,400 (ii) 0 500 1000 0 200 400 600 800 1000 1200 0 500 1000 0 5000 10000 15000 t 0 100 200 0 100 200 300 400 (iii) 0 100 200 0 200 400 600 800 1000 0 100 200 0 1 2 3 4 x 10 4 t 0 200 400 0 10 20 30 40 (iv) 0 200 400 0 20 40 60 0 200 400 0 5000 10000 15000 t 0 500 1000 1200 1400 1600 1500 x (i) 0 1000 2000 0 0.01 0.02 0.03 0.04 0.05 y 0 1000 2000 0 0.5 1 1.5 2 t v c = 34 c = 250 c = 800c = 30 Fig. 3: Time series solutions of system (2) for some particular values of c. The first column shows that infection is eradicated if virus reproduction is too low (c = 30). Subsequent columns demonstrate stable coexistence (c = 34), oscillatory coexistence (c = 250) and again stable coexistence of infection with increasing virus reproduction (c = 800). Here β = 0.0006 and other parameters are as in Fig. 2. Proof: The characteristic equation corre- sponding to M3 is ξ3 + B1ξ 2 + B2ξ + B3 = 0, (15) where B1 = −(m11 + m22 + m33), B2 = m11m22 + m22m33 + m33m11 −n23l32 −m12n21, B3 = (m11n23 −m13n21)l32 +(m12n21 −m11m22)m33. (16) We have B1B2−B3 =−m211m22−m 2 11m33−m 2 22m11 −m222m33 −m 2 33m11 −m 2 33m22 +m11m12n21 + m22n23l32 + m22m12n21 +m33n23l32 + m13n21l32 − 2m11m22m33. Following Routh-Hurwitz conditions, the neces- sary and sufficient conditions for E∗ to be locally asymptotically stable are B1 > 0, B3 > 0, B1B2 −B3 > 0. Hence the theorem. G. Hopf bifurcation analysis of the endemic equi- librium 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. Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 9 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... THEOREM II.5 If c crosses a critical value c∗ then there exists a Hopf bifurcation around E∗ when the following conditions hold: (i) B1(c∗) > 0, B3(c∗) > 0, (ii) B1(c∗)B2(c∗) −B3(c∗) = 0. (iii) B ′ 3(c∗)−B ′ 1(c∗)B2(c∗)−B1(c∗)B ′ 2(c∗) 6= 0. 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 [ξ2(c∗) + B2(c∗)][ξ(c∗) + B1(c∗)] = 0. (17) This equation has three roots ξ1(c∗) = i √ B2(c∗), ξ2(c∗) = −i √ B2(c∗) and ξ3(c∗) = −B1(c∗) < 0. For the existence of Hopf bifurcation at c = c∗, one needs to verify the transversality condition[ Re(ξ(c)) dc ] c=c∗ 6= 0. Note that, for all c, the roots of (15) are in general of the form ξ1(c) = l(c) + i m(c), ξ2(c) = l(c) − i m(c) and ξ3 = −B1(c). Substituting ξj(c) = l(c) ± i m(c), (j = 1, 2) in (15) and calculating the derivative, one obtains K̂(c∗) l ′ (c) − L̂(c) m ′ (c) + M̂(c) = 0, L̂(c∗) l ′ (c) + K̂(c) m ′ (c) + N̂(c) = 0, (18) where K̂(c) = 3l2(c) + 2B1(c)l(c) + B2(c) − 3m2(c), L̂(c) = 6l(c)m(c) + 2B1(c)m(c), M̂(c) = l2(c)B ′ 1(c) + B ′ 2(c)l(c) + B ′ 3(c) −B ′ 1(c)m 2(c), N̂(c) = 2l(c)m(c)B ′ 1(c) + m(c)B ′ 2(c). Since l(c∗) = 0 and m(c∗) = √ B2(c∗), we have K̂(c∗) = −2B2(c∗), L̂(c∗) = 2B1(c∗) √ B2(c∗), M̂(c∗) = B ′ 3(c∗) −B ′ 1(c∗)B2(c∗), N̂(c∗) = B ′ 2(c∗) √ B2(c∗). Solving for l ′ (c∗) from (18) and noting the condi- tion (iii) of Theorem II.5, we have [ Re(ξ1,2(c)) dc ] c=c∗ = l ′ (c∗) = − L̂(c∗)N̂(c∗) + K̂(c∗)M̂(c∗) K̂2(c∗) + L̂2(c∗) = B ′ 3(c∗) −B ′ 1(c∗)B2(c∗) −B1(c∗)B ′ 2(c∗) 2[B22(c∗) + B 2 1(c∗)] 6= 0. Thus, the transversality condition is satisfied. Therefore, the equilibrium point E∗ switches its stability (stable to unstable) through Hopf bifur- cation at c = c∗ when the conditions of Theorem II.5 hold. Hence the theorem. III. ANALYSIS OF THE DELAY-INDUCED SYSTEM It is worth mentioning that the system (4) has the same equilibrium points and same expression for the basic reproduction number R0 as in the case of non-delayed system (2). It is easy to check that for all τ1 ≥ 0 and τ2 ≥ 0, E1 is stable when R0 < 1 and therefore omitted. We are therefore interested to study the case R0 > 1, i.e., when the interior equilibrium E∗ exists, in presence of both delays. Consider the perturbed variables as x̄ = x − x∗, ȳ = y − y∗, v̄ = v − v∗. Dropping the bars for convenience, system (4) then reads dx dt =m11x(t) + m12y(t) + m13v(t) + F1, dy dt =n21x(t− τ1) + m22y(t) + n23v(t− τ1) + F2, dv dt =l32y(t− τ2) + m33v(t) + F3, (19) Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 10 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... τ 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 c 0 200 400 600 800 1,000 1,200 Stable E * Unstable E * Stable E 1 Fig. 4: Stability regions of the disease-free equilibrium E1 and infected equilibrium E∗ in τ2−c plane. Parameters are s = 10, µ = 0.01, r = 3, K = 1500, a = 0.007, b = 0.00001, p = 0.01, d1 = 1, d2 = 2.5, τ1 = 0 with β = 0.0006. where the coefficients are as in (14) and F1 =α21x 2(t) + α22v 2(t) + α23x(t)v(t) + α24x 3(t) + α25v 3(t) + α26x 2(t)v(t) + ..... , F2 =β21x 2(t− τ1) + β22v2(t− τ1) + β23x(t− τ1)v(t− τ1) + β24x 3(t− τ1) + β25v3(t− τ1) + β26x 2(t− τ1)v(t− τ1) + ..... , F3 =0, (20) with α21 = − r K + aβv∗(1 + bv∗) (1 + ax∗ + bv∗)3 , β21 = − aβv∗(1 + bv∗) (1 + ax∗ + bv∗)3 , α22 = bβx∗(1 + ax∗) (1 + ax∗ + bv∗)3 = −β22, α23 = − (1 + ax∗ + bv∗ + 2abx∗v∗)β (1 + ax∗ + bv∗)3 =−β23, α24 = − a2βv∗(1 + bv∗) (1 + ax∗ + bv∗)4 = −β24, α25 = − b2βx∗(1 + ax∗) (1 + ax∗ + bv∗)4 = −β25, α26 = aβ(1 + ax∗ + 2abx∗v∗ − b2v∗2) (1 + ax∗ + bv∗)4 =−β26. Corresponding linear system is then given by dx dt = m11x(t) + m12y(t) + m13v(t), dy dt = n21x(t− τ1) + m22y(t) + n23v(t− τ1), dv dt = l32y(t− τ2) + m33v(t), (21) and the corresponding characteristic equation is expressible as (ξ3 + Āξ2 + B̄ξ + C̄) + (D̄ξ + Ē)e−ξτ1 + (F̄ξ + Ḡ)e−(τ1+τ2) = 0, (22) Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 11 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... where Ā = −(m11 + m22 + m33), B̄ = m11m22 + m22m33 + m33m11, C̄ = −m11m22m33, D̄ = −m12n21, Ē = m12n21m33, F̄ = −n23l32, Ḡ = (m11n23 −m13n21)l32. (23) One can discuss the following cases in sequel. CASE 1. τ1 = τ2 = 0 THEOREM III.1 The interior equilibrium E∗ is locally asymptotically stable in absence of delays if R0 > 1, B1 > 0, B3 > 0 and B1B2 −B3 > 0, where B1, B2, B3 are given by (16). Proof: In this case, (22) takes the form ξ3 +Āξ2 +(B̄+D̄+F̄)ξ+(C̄+Ē+Ḡ) = 0. (24) One can notice that (23) and (16) are identical as Ā = −(m11 + m22 + m33) = B1, B̄ + D̄ + F̄ = m11m22 +m22m33 +m33m11−n23l32−m12n21 = B2 and C̄ + Ē + Ḡ = (m11n23 − m13n21)l32 + (m12n21 − m11m22)m33 = B3. Hence the result follows. CASE 2. τ2 > 0, τ1 = 0 THEOREM III.2 Assume that the conditions of Theorem III.1 are satisfied and the interior equi- librium E∗ is stable when τ1 = 0 = τ2. Then the following results are true. (i) E∗ is stable for all τ2 ≥ 0 if (28) has no positive root. (ii) If (28) has one positive root then there exists a critical value τ∗2 of τ2 such that E ∗ is stable for 0 < τ2 < τ∗2 and unstable for τ2 > τ ∗ 2 . System undergoes a Hopf bifurcation around E∗ at τ2 = τ ∗ 2 . (iii) If (28) has k positive roots, k = 2, 3, then there exists k critical values of τ2 given by τ∗2,j, 1 ≤ j ≤ k ; j ∈ N, where E ∗ will change its stability and a Hopf bifurcation will occur at each of these critical points. Proof: If τ2 > 0, τ1 = 0 then (22) has the form ξ3 + Āξ2 + (B̄ + D̄)ξ + (C̄ + Ē) + (F̄ξ + Ḡ)e−ξτ2 = 0. (25) We investigate if (25) has a pair of purely imag- inary roots of the form ξ = ±iω2, ω2 > 0, for some parametric conditions. Putting ξ = iω2, we obtain β1 cos(ω2τ2) + β2 sin(ω2τ2) −β3 + i (β2 cos(ω2τ2) −β1 sin(ω2τ2) −β4) = 0, where β1 = Ḡ, β2 = F̄ω2, β3 = Āω 2 2 − (C̄ + Ē), β4 = ω 3 2 −ω2(B̄ + D̄). (26) Separating real and imaginary parts, we get β1 cos(ω2τ2) + β2 sin(ω2τ2) = β3, β2 cos(ω2τ2) −β1 sin(ω2τ2) = β4. Eliminating τ2, these equations yield F2(ω2) = ω62 + K1 ω 4 2 + K2 ω 2 2 + K3 = 0, (27) where K1 = Ā2 − 2(B̄ + D̄),K2 = (B̄ + D̄)2 − 2Ā(C̄ + Ē) − F̄2 and K3 = (C̄ + Ē)2 − Ḡ2. We reduce the degree of Eq. (27) by setting h2 = ω22 and obtain F2(h2) = h32 + K1 h 2 2 + K2 h2 + K3 = 0. (28) If F2(h2) = 0 has no positive root then no change in stability will occur for τ2 > 0. Since E∗ was stable at τ2 = 0, it will remain stable for all τ2 > 0. If F2(h2) has one positive root h2 = h+2 (say) and the corresponding positive ω2 is ω∗2 = √ h+2 , then Eq. (25) will have a pair of purely imaginary roots of the form ξ = ±iω∗2 and the equilibrium E∗ will undergo a Hopf bifurcation at τ2 = τ∗2 , i.e., a family of periodic solutions will bifurcate from E∗ as τ2 crosses the value τ∗2 , where τ∗2 = 1 ω∗2 arccos ( β1β3 + β2β4 β21 + β 2 2 ) . (29) Here α1, α2, α3, α4 will be calculated from (26) with ω2 = ω∗2. Moreover, the transversality condition is given by[ Re ( dξ dτ2 )−1] ξ=iω∗2, τ2=τ ∗ 2 = F ′ 2(ω ∗2 2 ) ω∗22 (Ḡ 2 + F̄2ω∗22 ) 6= 0. (30) Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 12 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... 0 100 200 300 400 500 0 10 20 30 40 50 (ii) 0 100 200 300 400 500 0 50 100 0 100 200 300 400 500 0 1 2 3 4 x 10 4 t 0 200 400 600 800 1000 0 5 10 15 20 x (i) 0 200 400 600 800 1000 0 20 40 60 y 0 200 400 600 800 1000 0 0.5 1 1.5 2 x 10 4 t v τ 2 = 1.1τ 2 = 0.95 Fig. 5: Time evolutions of system (4) for different values of τ2. First panel shows stable behavior of the system populations when τ2 = 0.95(< τ∗2 ) and the second panel depicts the unstable behavior for τ2 = 1.1(> τ ∗ 2 ). Parameters are as in Fig. 4 with β = 0.0006, c = 800 and τ ∗ 2 = 0.9877 . If F2(h2) has k positive roots, k = 2, 3, given by h2 = h + 2,j, 1 ≤ j ≤ k ; j ∈ N, we will get k positive values of ω2 given by ω∗2,j = √ h+2,j , 1 ≤ j ≤ k ; j ∈ N. For each of these critical val- ues of ω∗2, Eq. (25) will have a pair of purely imaginary roots of the form ξ = ±iω∗2,j and the equilibrium E∗ will undergo a Hopf bifurcation at τ2 = τ ∗ 2,j , 1 ≤ j ≤ k ; j ∈ N, where τ∗2,j = 1 ω∗2,j arccos ( β1β3 + β2β4 β21 + β 2 2 ) , 1 ≤ j ≤ k ; j ∈ N, k = 2, 3. (31) As before, β1, β2, β3, β4 are calculated from (26) with ω2 = ω∗2,j , 1 ≤ j ≤ k ; j ∈ N. The transversality condition is given by[ Re ( dξ dτ2 )−1] ξ=iω∗2,j, τ2=τ ∗ 2,j = G ′ 1(ω ∗2 2,j) ω∗22,j(Ē + Ḡ) 2 + ω∗41,j(D̄ + F̄) 2 6= 0, 1 ≤ j ≤ k; j ∈ N,k = 2, 3. This completes the proof. CASE 3. τ1 > 0, τ2 > 0 THEOREM III.3 Assume that conditions of Theo- rem III.1 and either of Theorem III.2 (ii) or III.2 (iii) are satisfied. Corresponding to each positive root ω = ω̂1 of Eq. (35), there exists a critical value τ∗1 (τ2) of τ1 such that E ∗ will change its Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 13 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... 0 100 200 300 400 500 600 700 800 900 1000 0 200 400 600 800 1000 1200 1400 1600 c y τ 2 = 0.6 Fig. 6: Bifurcation diagram of the system (4) with c as the bifurcation parameter when τ2 = 0.6. Other parameters are as in Fig. 4. It demonstrates that infection can not persist if c < 38 and persists in oscillatory state for 38 < c < 663 and in stable state for c > 663. stability at each τ1 = τ∗1 (τ2), where τ2 is fixed and taken from the interval where E∗ was stable in the case τ2 > 0, τ1 = 0. A Hopf bifurcation will occur at τ1 = τ∗1 (τ2). Proof: Putting ξ = iω in (22) and equating real and imaginary parts, we have Ē cos(ωτ1) + D̄ω sin(ωτ1) + Ḡ cos(ω(τ1 + τ2)) + F̄ω sin(ω(τ1 + τ2)) + ( C̄ − Āω2 ) = 0, D̄ω cos(ωτ1) − Ē sin(ωτ1) + F̄ω cos(ω(τ1 + τ2)) − Ḡ sin(ω(τ1 + τ2)) + ( B̄ω −ω3 ) = 0. (32) Here τ2 is fixed in any of its stable range stated in Theorem III.2 (ii) or III.2 (iii), and τ1 is considered as a free parameter. From (32), we then get κ1 cos(ωτ1) + κ2 sin(ωτ1) = κ3, κ2 cos(ωτ1) −κ1 sin(ωτ1) = κ4, (33) where  κ1 = Ē + Ḡ cos(ωτ2) + F̄ω sin(ωτ2), κ2 = D̄ω − Ḡ sin(ωτ2) + F̄ω cos(ωτ2), κ3 = (Āω 2 − C̄), κ4 = ω 3 − B̄ω. (34) Eq. (33) yields J(ω) =ω6 + G1ω4 + G2ω2 + G3 + 2G4 sin(ωτ2) + 2G5 cos(ωτ2) = 0, (35) Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 14 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... where   G1 = Ā 2 − 2B̄, G2 = B̄ 2 − D̄2 − F̄2 − 2ĀC̄, G3 = C̄ 2 − (Ē2 + Ḡ2), G4 = (D̄Ḡ− ĒF̄)ω, G5 = −ĒḠ− D̄F̄ω2. One can note that if C̄2 < (Ē + Ḡ)2 then Eq. (35) will have at least one positive root ω̂1. For this ω̂1, (22) will have a pair of purely imaginary roots of the form ξ = ±iω̂1 and E∗ will change its stability. A Hopf bifurcation will occur at τ1 = τ∗1 (τ2), where τ∗2 (τ1) = 1 ω̂1 arccos ( κ1κ3 + κ2κ4 κ21 + κ 2 2 ) , (36) κ1,κ2,κ3,κ4 have to be calculated from (34) with ω = ω̂1. The transversality condition is also verified as Re [( dξ dτ2 )−1] ω=ω̂1, τ1=τ ∗ 1 (τ2) = �7(�1 + �3 + �5) + �8(�2 + �4 + �6) �27 + � 2 8 6= 0, where  �1 = B̄ − 3ω2, �2 = 2Āω, �3 = (D̄ − Ēτ1) cos(ωτ1) −D̄ωτ1 sin(ωτ1), �4 = −(D̄ − Ēτ1) sin(ωτ1) −D̄ωτ1 cos(ωτ1), �5 = (F̄ − Ḡ(τ1 + τ2)) cos(ω(τ1 + τ2)) −F̄ω(τ1 + τ2) sin(ω(τ1 + τ2)), �6 = −(F̄ − Ḡ(τ1 + τ2)) sin(ω(τ1 + τ2)) −F̄ω(τ1 + τ2) cos(ω(τ1 + τ2)), �7 = Ēω sin(ωτ1) − D̄ω2 cos(ωτ1) +Ḡ cos(ω(τ1 +τ2))+F̄ω sin(ω(τ1 +τ2)), �8 = Ēω sin(ωτ1) + D̄ω 2 sin(ωτ1) −Ḡ sin(ω(τ1 +τ2))+F̄ω cos(ω(τ1 +τ2)). Sign of the transversality condition will indicate the direction of Hopf bifurcation. Hence the theo- rem. 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 per- formed a sensitivity analysis using Latin Hyper- cube 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 esti- mating 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 parame- ter 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 de- lay 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 bi- furcation 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 R0 = 1. Infection, however, persists Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 15 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... 0 200 400 600 800 1000 0 5 10 15 20 x (i) 0 200 400 600 800 1000 0 20 40 60 y 0 200 400 600 800 1000 0 0.5 1 1.5 2 x 10 4 t v 0 100 200 300 400 500 0 5 10 15 20 (ii) 0 100 200 300 400 500 0 20 40 60 0 100 200 300 400 500 0 0.5 1 1.5 2 x 10 4 t τ 1 = 0.45 τ 2 = 0.6 τ 1 = 0.37 τ 2 = 0.6 Fig. 7: Time series representation of populations of the system (4) for τ1 = 0.37(< τ∗1 ), depicting stability of E∗ (first panel) and the same for τ1 = 0.45(> τ∗1 ), depicting instability of E ∗. Here β = 0.0006, c = 800, τ2 = 0.6 and other parameters are as in Fig. 4. for all higher values once c crosses the critical value. Persistence of host cells and virus particles in a stable state or oscillatory state depends on the value of virus reproduction rate. Populations coexist in a stable state for a small range of c after crossing the critical value and coexist in the oscillatory state for a longer intermediate range of c. All populations, however, coexist in a stable state when c is significantly large. Such behavior of the solutions for some fixed values of c are presented in 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 interdepen- dency 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 popula- tions 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 Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 16 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... 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 sat- urated 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 steady- state 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 repro- duction 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 coexis- tence of all population may be possible either for low or for very high virus production rate. For intermediate virus production rate populations Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 17 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... τ 1 0 0.2 0.4 0.6 0.8 1.0414 1.2 τ 2 0 0.2 0.4 0.6 0.8 0.9877 1.1 Unstable E * Stable E * Fig. 8: Stability and instability regions of the infected equilibrium E∗ in τ1−τ2 plane. The bifurcation line separates the plane into two stability zones. Parameters are as in Fig. 4 with c = 800 and β = 0.0006. 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 dif- ferent 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 per- sistence 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 consid- ers 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 oscilla- tory state. Furthermore, our analysis demonstrates that intra-and inter-patients variabilities of CD4+T cells and virus particles in case of HIV infec- Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 18 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... 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. ACKNOWLEDGEMENTS The research of Nandadulal Bairagi was sup- ported by the CSIR, Ref. NO.: 25(0294)/18/EMR- II. The research of Robert Hakl was supported by RVO 67985840. REFERENCES [1] C. B. Wilen, J. C. Tilton, and R. W. Doms, HIV: Cell binding and entry, Cold Spring Harb Perspect Med doi: 10.1101/cshperspect.a006866, 2012. [2] J. A. Levy, Pathogenesis of HIV infection, Microbiol Rev, vol. 57, pp. 183-289, 1993. [3] B. S. Stein et al., pH-independent HIV entry into CD4- positive T cells via virus envelope fusion to the plasma membrane, Cell, vol. 49, no. 5, pp. 659-668, 1987. [4] E. A. Berger, P. M. Murphy, J. M. Farber, Chemokine receptors as HIV-1 coreceptors: roles in viral entry, tropism, and disease, Annu Rev Immunol., vol. 17, pp. 657-700, 1991. [5] Y. H. Zheng, N. Lovsin and B. M. Peterlin, Newly iden- tified host factors modulate hiv replication, Immunol. Lett., vol. 97, pp. 225-234, 2005. [6] J. A. Levy, HIV and the pathogenesis of AIDS (second edition) American Society for Microbiology, Washing- ton, D. C. Press, 1998. [7] J. A. Zack et al., HIV-1 entry into quiescent primary lymphocytes: Molecular analysis reveals a labile latent viral structure, Cell, vol. 61, pp. 213-222, 1990. [8] J. A. Zack et al., Incompletely reverse-transcribed human immunodeficiency virus type 1 genomes in quiescent cells can function as intermediates in the retroviral cycle, J. Virol., vol. 66, pp. 1717-1725, 1992. [9] P. Essunger and A. S. Perelson, Modelling HIV infection of CD4+T-cell subpopulations, J. Theo. Biol., vol. 170, pp. 367-391, 1994. [10] J. A. Hiscox, RNA viruses: hijacking the dynamic nucle- olus, Nature Reviews Microbiology, vol. 5, pp. 119-127, 2007. [11] L. Rong, M. A. Gilchrist, Z. Feng, A. S. Perelson, Modeling within host HIV-1 dynamics and the evolution of drug resistance: Trade offs between viral enzyme function and drug susceptibility, J. Theoret. Biol., vol. 247, pp. 804-818, 2007. [12] P. K. Srivastava, and P. Chandra, Modeling the dynam- ics of HIV and CD4+T cells during primary infection, Nonlinear Ana. RWA., vol. 11, pp. 612-618, 2010. [13] D. Adak and N. Bairagi, Bifurcation analysis of a multi- delayed HIV model in presence of immune response and understanding of in-host viral dynamics, Math. Anal. Applied Sciences, vol. 42, no. 12, pp. 4256-4272, 2019. [14] X. Zhou, X. Song and X. Shi, A differential equation model of HIV infection of CD4+T with cure rate, J. Math. Anal. Appl., vol. 342, pp. 1342-1355, 2008. [15] J. Yang, X. Wang and F. Zhang, A Differential Equa- tion Model of HIV Infection of CD4+T-Cells with de- lay, Discrete Dynamics in Nature and Society, doi:10.1155/2008/903678, 2008. [16] A. V. M. Herz et al., Viral dynamics in vivo; limitations on estimates of intracellular delay and virus decay, Proc. Natl. Acad. Sci. USA, vol. 93, pp. 7247-7251, 1996. [17] A. Perelson et al., HIV-1 dynamics in vivo: Virion clear- ance rate, infected cell life-span, and viral generation time, Science, vol. 271, pp. 1582-1586, 1996. [18] S. Bonhoeffer et al., Virus dynamics and drug therapy, Proc. Natl. Acad. Sci. USA, vol. 94, pp. 6971-6976, 1997. [19] M. A. Nowak and R. M. May, Virus Dynamics: Mathe- matical Principles of Immunology and Virology, Oxford University, Oxford, 2000. [20] S. Fisson et al., Continuous activation of autoreactive CD4+CD25+ regulatory T cells in the steady state, J Exp Med, vol. 198, pp.737-746, 2003. [21] M. L. Chan et al., Limited CD4+ T cell proliferation leads to preservation of CD4+T cell counts in SIV- infected sooty mangabeys, Proc. R. Soc. B., vol. 277, pp. 3773-3781, 2010. [22] H. Shu, L. Wang, J. Watmough, Global stability of a nonlinear viral infection model with infinitely distributed intracellular delays and CTL immune responses, SIAM J. Appl. Math., vol. 73, pp. 1280-1302, 2013. [23] L. Wang and M. Y. Li, Mathematical analysis of the global dynamics of a model for HIV infection of CD4+T cells, Math. Biosci., vol. 200, pp. 44-57, 2006. [24] R. M. May and R. M. Anderson, Population biology of infectious disease: II, Nature, vol. 280, pp. 455-461, 1979. [25] D. Li and W. Ma, Asymptotic properties of a HIV-1 infection model with time delay, J. Math. Anal. Appl., vol. 335, pp.683-691, 2007. [26] R. Xu, Global stability of an HIV-1 infection model with saturation infection and intracellular delay, J. Math. Anal. Appl., vol. 375, pp. 75-81, 2011. [27] R. J. De Boer, Understanding the failure of CD8+T-cell vaccination against Simian/Human Immunodeficiency Virus, J. Viro., vol. 81, pp. 2838-2848, 2007. [28] N. Bairagi and D. Adak, Global analysis of HIV-1 dynamics with Hill type infection rate and intracellular delay, Applied Math. Modelling, vol. 38, pp. 5047- 5066, 2014. [29] D. Adak and N. Bairagi, Analysis and computation of multi-pathways and multi-delays HIV-1 infection model, Applied Math. Modelling, vol. 54, pp. 517-537, 2018. [30] Y. Zhang and Z. Xu, Dynamics of a diffusive HBV model Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 19 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 D. Adak, N. Bairagi, R. Hakl, Accounting for multi-delay effects in an HIV-1 infection model with ... with delayed Beddington-DeAngelis response, Nonlinear Anal. RWA., vol. 15, pp. 118-139, 2014. [31] H. Ansari and M. Hesaaraki, A with-in host dengue infection model with immune response and Beddington- DeAngelis incidence rate, Appl. Math., vol. 3, pp. 177- 184, 2012. [32] M. Y. Li and H. Shu, Global dynamics of an in-host viral model with intracellular delay, Bull. Math. Biol., vol. 72, pp. 1492-1505, 2010. [33] Y. Wang, Y. Zhou, J. Wub, J. Heffernan, Oscillatory viral dynamics in a delayed HIV pathogenesis model, Math. Biosci., vol. 219, pp. 104-112, 2009. [34] H. Zhu and X. Zou, Impact of delays in cell infection and virus production on HIV-1 dynamics, Math. Medi. and Bio., vol. 25, pp. 99-112, 2008. [35] Q. Sun and L. Min, Dynamics analysis and simulation of a modified HIV infection model with a saturated infection rate, Computa. Mat. Methods in Medi., Article ID 145162, vol. 2014. [36] P. K. Srivastava, M. Banerjee and P. Chandra, A primary infection model for HIV and immune response with two discrete time delays, Diff. Equ. Dyn. Syst., vol. 18, pp. 385-399, 2010. [37] K. A. Pawelek, S. Liu, F. Pahlevani and L. Rong, A model of HIV-1 infection with two time delays: Math- ematical analysis and comparison with patient data, Math. Bios., vol. 235, pp. 98-109, 2012. [38] A. Alshorman, X. Wang, M. M. Joseph andL. Rong, Analysis of HIV models with two time delays, J. Biol. Dynamics, vol. 11(sup1), pp. 40-64, 2017. [39] N. Sweilam, A. M. Seham, S. Shatta and D. Baleanu, Numerical study for a novel variable-order multiple time delay awareness programs mathematical model, Appl. Numer. Math., vol. 158, pp. 212-235, 2020. [40] A. S. Perelson, D. E. Kirschner, D. E. Boer, Dynamics of HIV infection of CD4+T cells, Math. Bios., vol. 114, pp. 81-125, 1993. [41] N. Nagumo, Uber die lage der integralkurven gewn- licherdifferantialgleichungen, Proc. Phys. Math. Soc. Jpn., vol. 24, pp. 551, 1942. [42] P. Van Den Driessche, and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for com- partmental models of disease transmission, Math. Bios., vol. 180, pp. 2948, 2002. [43] J. P. LaSalle, The stability of dynamical systems, in: Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, PA, 1976. [44] J. K. Hale and P. Waltman, Persistence in infinite- dimensional systems, SIAM J. Math. Anal., vol. 20, pp. 388-396, 1989. [45] F. R. Gantmacher, The Theory of Matrices, Chelsea Publ. Co., New York, 1959. [46] M. Y. Li and L. Wang, Global stability in some SEIR epidemic models, Mathematical Approaches for Emerging and Reemerging Infectious Diseases: Models, Methods, and Theory, Springer, New York, pp. 295-311, 2002. [47] G. H. Butler and P. Waltman, Persistence in dynamical systems, Proc. Amer. Math. Soc., vol. 96, pp. 425-430, 1986. [48] M. W. Hirch, Systems of differential equations which are competitive or cooperative, IV: structrual stabilitites in three dimensional systems, SIAM J Math. Anal., vol. 21, pp. 1225-1234, 1990. [49] X. Wang and X. Song, Global stability and periodic solution of a model for HIV infection of CD4+T cells, Appl. Math. Compu., vol. 189, no. 2, pp. 1331-1340, 2007. [50] P. De Leenheer and H. L. Smith, Virus dynamics: a global analysis, SIAM J. Appl. Math., vol. 63, pp. 1313-1327, 2003. [51] N. M. Dixit, A. S. Perelson, Complex patterns of viral load decay under antiretroviral therapy: influence of pharmacokinetics and intracellular delay, J. Theor. Biol., vol. 226, pp. 95, 2004. [52] R. V. Culshaw and S. Ruan, A delay-differential equa- tion model of HIV infection of CD4+T-cells, Math. Biosci., vol. 165, pp. 27-39, 2000. [53] P. Nelson, J. Murray, A. Perelson, A model of HIV-1 pathogenesis that includes an intracellular delay, Math. Biosci., vol. 163, pp. 201, 2000. [54] R. Culshaw, S. Ruan, R. Spiteri, Optimal HIV treatment by maximising immune response, J. Math. Biol., vol. 48, pp. 545-562, 2004. [55] L. Rong, M. A. Gilchrist, Z. Feng and A. S. Perelson, Modeling within-host HIV-1 dynamics and the evolution of drug resistance: trade-offs between viral enzyme func- tion and drug susceptibility, J. Theo. Biol., vol. 247, no. 4, pp. 804-818, 2007. [56] B. Buonomo and C. Vargas-De-Leon, Global stability for an HIV-1 infection model including an eclipse stage of infected cells, J. Math. Anal. Appl., vol. 385, pp. 709-720, 2012. [57] X. Zhou, X. Song and X. Shi, A differential equation model of HIV infection of CD4+T cells with cure rate, J. Math. Anal. Appl., vol. 342, pp. 1342-1355, 2008. [58] S. M. Blower and H. Dowlatabadi, Sensitivity and uncertainty analysis of complex models of disease trans- mission: an HIV model, as an example, Int. Stat. Rev., vol. 62, pp. 229-243, 1994. [59] J. M. Raboud et al., Quantification of the variation due to laboratory and physiologic sources in CD4 lymphocyte counts of clinically stable HIV-infected individuals, J. Acqu. Immu. Deffi. Syndr. Hum. Retro., vol. 10, pp. 67- 73, 1995. [60] B. G. Williams et al., HIV infection, antiretroviral therapy, and CD4 cell count distributions in African populations, J. Infect. Diss., vol. 194, pp. 1450-1458, 2006. [61] A. C. Crampin et al., Normal range of CD4 cell counts and temporal changes in two HIV negative Malawian populations, Open AIDS J., vol. 5, pp. 74-79, 2011. Biomath 9 (2020), 2012297, http://dx.doi.org/10.11145/j.biomath.2020.12.297 Page 20 of 20 http://dx.doi.org/10.11145/j.biomath.2020.12.297 The model Analysis of the non-delayed system Positivity and boundedness of solutions Equilibria, stability and permanence Basic reproductive number Local stability of disease-free equilibrium Permanence of the system Local stability of the interior equilibrium Hopf bifurcation analysis of the endemic equilibrium Analysis of the delay-induced system Numerical simulations Discussion References