Elaiw et al. Advances in Difference Equations (2018) 2018:85 https://doi.org/10.1186/s13662-018-1523-0

RESEARCH

Open Access

Effect of cellular reservoirs and delays on the global dynamics of HIV A. M. Elaiw1* , E. K. Elnahary2 and A. A. Raezah3 *

Correspondence: [email protected] 1 Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia Full list of author information is available at the end of the article

Abstract We investigate general HIV infection models with three types of infected cells: latently infected cells, long-lived productively infected cells, and short-lived productively infected cells. We incorporate three discrete or distributed time delays into the models. Moreover, we consider the effect of humoral immunity on the dynamical behavior of the HIV. The HIV-target incidence rate, production/proliferation, and removal rates of the cells and HIV are represented by general nonlinear functions. We show that the solutions of the proposed model are nonnegative and ultimately bounded. We derive two threshold parameters which fully determine the existence and stability of the three steady states of the model. Using Lyapunov functionals, we establish the global stability of the steady states of the model. The theoretical results are confirmed by numerical simulations. MSC: 34D20; 34D23; 37N25; 92B05 Keywords: HIV infection; HAART; Global stability; Humoral immune response; Latency; Viral reservoirs

1 Introduction Mathematical modeling and analysis of within-host human immunodeficiency virus (HIV) dynamics have become one of the hot topics during the last decades [1–31]. These works can help researchers better understand the HIV dynamical behavior and provide new suggestions for clinical treatment. Most of the mathematical models presented in the literature have focused on modeling the interaction between three main compartments: uninfected CD4+ T cells (s), infected cells (y), and free HIV particles (p). Other models have differentiated between latent and active infected cells by introducing a new variable (w) for the latently infected cells [32–37]. In [38], an HIV mathematical model has been presented by considering three types of infected cells: latently infected cells (w), short-lived productively infected cells (y), and long-lived productively infected cells (u) as follows:   s(t) – (1 – ε1 )(λ1 + λ2 + λ3 )s(t)p(t), s˙(t) = ρ – β1 s(t) + ωs(t) 1 – smax

(1)

˙ = (1 – ε1 )λ1 s(t)p(t) – (a1 + β2 )w(t), w(t)

(2)

y˙ (t) = (1 – ε1 )λ2 s(t)p(t) + a1 w(t) – β3 y(t),

(3)

© The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 2 of 36

˙ = (1 – ε1 )λ3 s(t)p(t) – β4 u(t), u(t)

(4)

p˙ (t) = (1 – ε2 )Nβ3 y(t) + (1 – ε2 )Mβ4 u(t) – β5 p(t),

(5)

where ρ is the creation rate of the uninfected CD4+ T cells, βi , i = 1, . . . , 5, are the death rate constants of the five compartments s, w, y, u, and p, respectively. The model incorporates reverse transcriptase inhibitor (RTI) with efficacy ε1 and protease inhibitor (PI) with efficacy ε2 , where ε1 , ε2 ∈ [0, 1]. The parameters ω and smax are the maximum proliferation rate and the maximum level of concentration of the uninfected CD4+ T cells, respectively. The latently infected cells are activated at rate a1 w. The parameters N and M are the average numbers of HIV particles generated in the lifetime of the short-lived and long-lived infected cells, respectively. The uninfected CD4+ T cells become infected with infectivity λ1 + λ2 + λ3 . Elaiw et al. [27] have generalized the above model by incorporating the humoral immune response and considering general nonlinear functions for the generation and removal rates of all compartments:     s˙(t) = π s(t) – (1 – ε1 )(λ1 + λ2 + λ3 )χ s(t), p(t) ,     ˙ = (1 – ε1 )λ1 χ s(t), p(t) – (a1 + β2 )g1 w(t) , w(t)       y˙ (t) = (1 – ε1 )λ2 χ s(t), p(t) + a1 g1 w(t) – β3 g2 y(t) ,     ˙ = (1 – ε1 )λ3 χ s(t), p(t) – β4 g3 u(t) , u(t)     p˙ (t) = (1 – ε2 )Nβ3 g2 y(t) + (1 – ε2 )Mβ4 g3 u(t)       – β5 g4 p(t) – qg4 p(t) g5 x(t) ,       x˙ (t) = rg4 p(t) g5 x(t) – β6 g5 x(t) ,

(6) (7) (8) (9)

(10) (11)

where x represents the concentration of the B cells. π , χ , gi , i = 1, . . . , 5, are general nonlinear functions. Model (6)–(11) assumes that, once the HIV contacts a CD4+ T cell, it becomes infected in the same time. Neglecting the time delays is an unrealistic assumption (see, e.g., [29, 30]). The aim of this paper is to propose HIV infection models which improve model (6)–(11) by taking into account three time delays, discrete or distributed. We derive two threshold parameters and present some mild sufficient conditions for the existence and global stability of the steady states of the models.

2 HIV dynamics model with discrete delays We formulate a nonlinear HIV dynamics model with latent reservoirs, humoral immunity, and discrete time delays:     s˙ = π s(t) – (1 – ε1 )(λ1 + λ2 + λ3 )χ s(t), p(t) ,     ˙ = (1 – ε1 )λ1 e–μ1 τ1 χ s(t – τ1 ), p(t – τ1 ) – (a1 + β2 )g1 w(t) , w       y˙ = (1 – ε1 )λ2 e–μ2 τ2 χ s(t – τ2 ), p(t – τ2 ) + a1 g1 w(t) – β3 g2 y(t) ,     u˙ = (1 – ε1 )λ3 e–μ3 τ3 χ s(t – τ3 ), p(t – τ3 ) – β4 g3 u(t) ,

(12) (13) (14) (15)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 3 of 36

          p˙ = (1 – ε2 )Nβ3 g2 y(t) + (1 – ε2 )Mβ4 g3 u(t) – β5 g4 p(t) – qg4 p(t) g5 x(t) ,       x˙ = rg4 p(t) g5 x(t) – β6 g5 x(t) .

(16) (17)

All the parameters are positive. Here, τ1 is the time between viral entry and latent infection (i.e., the integration of viral DNA into cell’s DNA has finished), while τ2 and τ3 are the times between viral entry and viral production form short-lived productively infected and long-lived productively infected cells, respectively. The factor e–μi τi accounts for the loss of target cells during the delay period of length τi , where μi is a constant. Let us define λi = (1 – ε1 )λi , i = 1, 2, 3, N = (1 – ε2 )N, M = (1 – ε2 )M, and λ = λ1 + λ2 + λ3 . Functions χ , π , gi , i = 1, . . . , 5, are continuously differentiable and satisfy the following hypotheses: (H1). (i) There exists s0 such that π(s0 ) = 0, π(s) > 0 for s ∈ [0, s0 ); (ii) π  (s) < 0 for s ∈ (0, ∞); (iii) There are b > 0 and b > 0 such that π(s) ≤ b – bs for s ∈ [0, ∞). (H2). (i) χ(s, p) > 0 and χ(0, p) = χ(s, 0) = 0 for s, p ∈ (0, ∞); (ii) ∂χ(s,p) > 0, ∂χ(s,p) > 0, and ∂χ(s,0) > 0 for all s, p ∈ (0, ∞); ∂s ∂p ∂p (iii) ( ∂χ(s,0) ) > 0 for s ∈ (0, ∞). ∂p (H3). (i) gj (η) > 0 for η ∈ (0, ∞), gj (0) = 0, j = 1, . . . , 5; (ii) gj (η) > 0 for η ∈ (0, ∞), j = 1, 2, 3, 5, g4 (η) > 0, for η ∈ [0, ∞); (iii) there are αj > 0, j = 1, . . . , 5, such that gj (η) ≥ αj η for η ∈ [0, ∞). ∂ χ(s,p) (H4). (i) ∂p ( g4 (p) ) ≤ 0 for p ∈ (0, ∞). We consider system (12)–(17) with the initial conditions: s(θ ) = ϕ1 (θ ),

w(θ ) = ϕ2 (θ ),

y(θ ) = ϕ3 (θ ),

u(θ ) = ϕ4 (θ ),

p(θ ) = ϕ5 (θ ),

x(θ ) = ϕ6 (θ ),

ϕj (θ ) ≥ 0, θ ∈ [–ς, 0],   ϕj (θ ) ∈ C [–ς, 0], R6≥0 ,

(18)

j = 1, . . . , 6,

where ς = max{τ1 , τ2 , τ3 } and C is the Banach space of continuous functions mapping the interval [–ς, 0] into R6≥0 . Then the uniqueness of the solution for t > 0 is guaranteed [39].

2.1 Preliminaries Lemma 1 Let hypotheses (H1)–(H3) be valid, then the solutions of system (12)–(17) are nonnegative and ultimately bounded. ˙ = L(k(t)), where k = (s, w, y, Proof Let us write system (12)–(17) in the matrix form k(t) u, p, x)T , L = (L1 , L2 , L3 , L4 , L5 , L6 )T , and ⎛

⎞ L1 (k(t)) ⎜ ⎟ ⎜L2 (k(t))⎟ ⎜ ⎟   ⎜L3 (k(t))⎟ ⎜ ⎟, L k(t) = ⎜ ⎟ ⎜L4 (k(t))⎟ ⎜ ⎟ ⎝L5 (k(t))⎠ L6 (k(t))

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 4 of 36

⎞ π(s(t)) – λχ(s(t), p(t)) ⎟ ⎜ λ1 e–μ1 τ1 χ(s(t – τ1 ), p(t – τ1 )) – (a1 + β2 )g1 (w(t)) ⎟ ⎜ ⎟ ⎜ ⎜ λ2 e–μ2 τ2 χ(s(t – τ2 ), p(t – τ2 )) + a1 g1 (w(t)) – β3 g2 (y(t)) ⎟ ⎟. L=⎜ ⎟ ⎜ λ3 e–μ3 τ3 χ(s(t – τ3 ), p(t – τ3 )) – β4 g3 (u(t)) ⎟ ⎜ ⎟ ⎜ ⎝Nβ3 g2 (y(t)) + Mβ4 g3 (u(t)) – β5 g4 (p(t)) – qg4 (p(t))g5 (x(t))⎠ rg4 (p(t))g5 (x(t)) – β6 g5 (x(t)) ⎛

We have   Li k(t) k (t)∈R6 ≥ 0, i

>0

i = 1, . . . , 6.

Using Lemma 2 in [40], the solutions of system (12)–(17) with the initial states (18) satisfy X(t) ∈ R6≥0 for all t ≥ 0. The nonnegativity of the model’s solution implies that lim supt→∞ s(t) ≤ M1 , where M1 = bb¯ . Let T(t) = Ne–μ1 τ1 s(t – τ1 ) + Ne–μ2 τ2 s(t – τ2 ) + Me–μ3 τ3 s(t – τ3 ) 1 q + Nw(t) + Ny(t) + Mu(t) + p(t) + x(t), 2 2r then

    ˙ = Ne–μ1 τ1 π s(t – τ1 ) – λχ s(t – τ1 ), p(t – τ1 ) T(t)

    + Ne–μ2 τ2 π s(t – τ2 ) – λχ s(t – τ2 ), p(t – τ2 )

    + Me–μ3 τ3 π s(t – τ3 ) – λχ s(t – τ3 ), p(t – τ3 )    

+ N λ1 e–μ1 τ1 χ s(t – τ1 ), p(t – τ1 ) – (a1 + β2 )g1 w(t)  

    + N λ2 e–μ2 τ2 χ s(t – τ2 ), p(t – τ2 ) + a1 g1 w(t) – β3 g2 y(t)    

+ M λ3 e–μ3 τ3 χ s(t – τ3 ), p(t – τ3 ) – β4 g3 u(t)           1 Nβ3 g2 y(t) + Mβ4 g3 u(t) – β5 g4 p(t) – qg4 p(t) g5 x(t) 2      q  rg4 p(t) g5 x(t) – β6 g5 x(t) + 2r   

≤ Ne–μ1 τ1 b – bs(t – τ1 ) + Ne–μ2 τ2 b – bs(t – τ2 ) + Me–μ3 τ3 b – bs(t – τ3 ) +

1 1 1 q – Nβ2 α1 w(t) – Nβ3 α2 y(t) – Mβ4 α3 u(t) – β5 α4 p(t) – β6 α5 p(t) 2 2 2 2r   ≤ b Ne–μ1 τ1 + Ne–μ2 τ2 + Me–μ3 τ3  – σ Ne–μ1 τ1 s(t – τ1 ) + Ne–μ2 τ2 s(t – τ2 ) + Me–μ3 τ3 s(t – τ3 ) 1 q + Nw(t) + Ny(t) + Mu(t) + p(t) + x(t) 2 2r



≤ b(2N + M) – σ T(t), where σ = min{b, β2 α1 , 12 β3 α2 , 12 β4 α3 , β5 α4 , β6 α5 }. Then lim supt→∞ T(t) ≤ nonnegativity of the system’s variables implies that lim sup w(t) ≤ t→∞

b(2N + M) = M2 , Nσ

b(2N+M) . σ

The

Elaiw et al. Advances in Difference Equations (2018) 2018:85

lim sup y(t) ≤

b(2N + M) = M2 , Nσ

lim sup u(t) ≤

b(2N + M) = M3 , Mσ

lim sup p(t) ≤

2b(2N + M) = M4 , σ

lim sup x(t) ≤

2rb(2N + M) = M5 . qσ

t→∞

t→∞

t→∞

t→∞

Therefore, s(t), w(t), y(t), u(t), p(t), and x(t) are ultimately bounded.

Page 5 of 36



According to Lemma 1, we can show that the region   = (s, w, y, u, p, x) ∈ C 6 : s ≤ M1 , w ≤ M2 , y ≤ M2 ,  u ≤ M3 , p ≤ M4 , x ≤ M5

(19)

is positively invariant with respect to system (12)–(17), where φ = supt→∞ φ(t). Lemma 2 Suppose that hypotheses (H1)–(H4) are valid (12)–(17), then there exist two bifurcation parameters R0 and R1 with R0 > R1 > 0 such that (i) if R0 ≤ 1, then there exists only one steady state 0 ; (ii) if R1 ≤ 1 < R0 , then there exist two steady states 0 and 1 ; (iii) if R1 > 1, then there exist three steady states 0 , 1 , and 2 . Proof Let (s, w, y, u, p, x) be any steady state of (12)–(17) satisfying the following equations: 0 = π(s) – λχ(s, p),

(20)

0 = λ1 e–μ1 τ1 χ(s, p) – (a1 + β2 )g1 (w),

(21)

0 = λ2 e–μ2 τ2 χ(s, p) + a1 g1 (w) – β3 g2 (y),

(22)

0 = λ3 e–μ3 τ3 χ(s, p) – β4 g3 (u),

(23)

0 = Nβ3 g2 (y) + Mβ4 g3 (u) – β5 g4 (p) – qg4 (p)g5 (x),

(24)

0 = rg4 (p)g5 (x) – β6 g5 (x).

(25)

From Eq. (25) we have two possible solutions: g5 (x) = 0 and g4 (p) = β6 /r. The first possibility g5 (x) = 0 implies that x = 0. Hypothesis (H3) implies that gi–1 , i = 1, . . . , 5, exist, strictly increasing and gi–1 (0) = 0. Let us define  ψ1 (s) = g1–1 

 λ1 e–μ1 τ1 π(s) , λ(a1 + β2 )

 a1 λ1 e–μ1 τ1 + (a1 + β2 )λ2 e–μ2 τ2 π(s) , β3 λ(a1 + β2 )    –μ3 τ3  λ3 e γ π(s) , ψ4 (s) = g4–1 π(s) , ψ3 (s) = g3–1 β4 λ λ

ψ2 (s) = g2–1

(26)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 6 of 36

where γ=

N(a1 λ1 e–μ1 τ1 + (a1 + β2 )λ2 e–μ2 τ2 ) + Mλ3 (a1 + β2 )e–μ3 τ3 . β5 (a1 + β2 )

It follows from Eqs. (20)–(24) that w = ψ1 (s),

y = ψ2 (s),

u = ψ3 (s),

p = ψ4 (s).

(27)

Obviously, ψi (s) > 0 for s ∈ [0, s0 ) and ψi (s0 ) = 0, i = 1, . . . , 4. From Eqs. (20), (26), and (27) we obtain     γ χ s, ψ4 (s) – g4 ψ4 (s) = 0.

(28)

Equation (28) admits a solution s = s0 which gives the infection-free steady state 0 (s0 , 0, 0, 0, 0, 0). Let     1 (s) = γ χ s, ψ4 (s) – g4 ψ4 (s) = 0. It is clear from hypotheses (H1) and (H2) that 1 (0) = –g4 (ψ4 (0)) < 0 and 1 (s0 ) = 0. Moreover, 1 (s0 ) = γ



 ∂χ(s0 , 0) ∂χ(s0 , 0)  + ψ4 (s0 ) – g4 (0)ψ4 (s0 ). ∂s ∂p

We note from hypothesis (H2) that 1 (s0 ) = ψ4 (s0 )g4 (0)



∂χ(s0 ,0) ∂s

= 0. Then

 γ ∂χ(s0 , 0) – 1 . g4 (0) ∂p

From Eq. (26), we get 1 (s0 ) =

  γ  γ ∂χ(s0 , 0) π (s0 )  –1 . λ g4 (0) ∂p

Hence, from hypothesis (H1), we have π  (s0 ) < 0. Therefore, if g γ(0) ∂χ(s∂p0 ,0) > 1, then 1 (s0 ) < 4 0, and there exists s1 ∈ (0, s0 ) such that 1 (s1 ) = 0. Hypotheses (H1)–(H3) imply that w1 = ψ1 (s1 ) > 0,

y1 = ψ2 (s1 ) > 0,

u1 = ψ3 (s1 ) > 0,

p1 = ψ4 (s1 ) > 0.

(29)

It means that a humoral-inactivated infection steady state 1 (s1 , w1 , y1 , u1 , p1 , 0) exists when g γ(0) ∂χ(s∂p0 ,0) > 1. Let us define the basic infection reproduction number as follows: 4

R0 =

γ g4 (0)

∂χ(s0 , 0) . ∂p

The other solution of Eq. (25) is g4 (p2 ) = βr6 , which yields p2 = g4–1 ( βr6 ) > 0. Substitute p = p2 in Eq. (20) and let 2 (s) = π(s) – λχ(s, p2 ) = 0. According to hypotheses (H1) and (H2),

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 7 of 36

2 is strictly decreasing, 2 (0) = π(0) > 0 and 2 (s0 ) = –λχ(s0 , p2 ) < 0. Thus, there exists unique s2 ∈ (0, s0 ) such that 2 (s2 ) = 0. It follows from Eqs. (24) and (27) that  w2 = ψ1 (s2 ) > 0,  x2 = g5–1

y2 = ψ2 (s2 ) > 0,

u2 = ψ3 (s2 ) > 0,

p2 = g4–1

  χ(s2 , p2 ) β5 γ –1 . q g4 (p2 )

β6 r

 > 0,

2 ,p2 ) > 1. Now we define the humoral immune response activation Thus, x2 > 0 when γ χ(s g4 (p2 ) number as follows:

R1 = γ

χ(s2 , p2 ) . g4 (p2 )

If R1 > 1, then x2 = g5–1 ( βq5 (R1 –1)) > 0, and there exists a humoral-activated infection steady state 2 (s2 , w2 , y2 , u2 , p2 , x2 ). Clearly, from hypotheses (H2) and (H4), we have R1 = γ

χ(s2 , p) γ ∂χ(s2 , 0) γ ∂χ(s0 , 0) χ(s2 , p2 ) ≤ γ lim+ =  <  = R0 . p→0 g4 (p2 ) g4 (p) g4 (0) ∂p g4 (0) ∂p



We will use the following equalities throughout the paper:   χ(s(t – τ1 ), p(t – τ1 )) ln χ(s, p)     ˆ g1 (w)χ(s(t – τ1 ), p(t – τ1 )) χ(ˆs, pˆ ) = ln + ln g1 (w)χ(ˆs, pˆ ) χ(s, pˆ )       g4 (ˆp)g2 (y) g2 (ˆy)g1 (w) g4 (p)χ(s, pˆ ) + ln + ln , + ln ˆ g4 (ˆp)χ(s, p) g4 (p)g2 (ˆy) g2 (y)g1 (w)   χ(s(t – τ2 ), p(t – τ2 )) ln χ(s, p)     g2 (ˆy)χ(s(t – τ2 ), p(t – τ2 )) χ(ˆs, pˆ ) = ln + ln g2 (y)χ(ˆs, pˆ ) χ(s, pˆ )     g4 (ˆp)g2 (y) g4 (p)χ(s, pˆ ) + ln + ln , g4 (p)g2 (ˆy) g4 (ˆp)χ(s, p)   χ(s(t – τ3 ), p(t – τ3 )) ln χ(s, p)     ˆ – τ3 ), p(t – τ3 )) χ(ˆs, pˆ ) g3 (u)χ(s(t + ln = ln g3 (u)χ(ˆs, pˆ ) χ(s, pˆ )     g3 (u)g4 (ˆp) g4 (p)χ(s, pˆ ) + ln + ln . ˆ 4 (p) g3 (u)g g4 (ˆp)χ(s, p)

(30)

2.2 Global properties The following theorems investigate the global stability of the steady states of system (12)– (17). Let us define the function F : (0, ∞) → [0, ∞) as F(z) = z – 1 – ln z. Denote (s, w, y, u, p, x) = (s(t), w(t), y(t), u(t), p(t), x(t)).

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 8 of 36

Theorem 1 If R0 ≤ 1 and hypotheses (H1)–(H4) are valid, then 0 is globally asymptotically stable (GAS). Proof Define a Lyapunov functional V0 as follows:  V0 = s – s0 – 

s

lim+

s0 p→0

χ(s0 , p) dη + 1 w + 2 y + 3 u χ(η, p)

τ1

  λ1 e–μ1 τ1 χ s(t – θ ), p(t – θ ) dθ

τ2

  λ2 e–μ2 τ2 χ s(t – θ ), p(t – θ ) dθ

τ3

  λ3 e–μ3 τ3 χ s(t – θ ), p(t – θ ) dθ + 4 p + 5 x,

+ 1 0

 + 2

0

 + 3

0

where λ = λ1 1 e–μ1 τ1 + λ2 2 e–μ2 τ2 + λ3 3 e–μ3 τ3 , 2 = N4 ,

3 = M4 ,

(a1 + β2 )1 = a1 2 ,

(31)

q4 = r5 .

The solution of Eqs. (31) is given by 1 =

a1 Nλ , γβ5 (a1 + β2 )

We calculate

dV0 dt

2 =

Nλ , γβ5

3 =

Mλ , γβ5

4 =

λ , γβ5

5 =

qλ . rγβ5

(32)

along the trajectories of (12)–(17) as follows:

   χ(s0 , p)  dV0 = 1 – lim+ π(s) – λχ(s, p) p→0 χ(s, p) dt  –μ τ    + 1 λ1 e 1 1 χ s(t – τ1 ), p(t – τ1 ) – (a1 + β2 )g1 (w)     + 2 λ2 e–μ2 τ2 χ s(t – τ2 ), p(t – τ2 ) + a1 g1 (w) – β3 g2 (y)     + 3 λ3 e–μ3 τ3 χ s(t – τ3 ), p(t – τ3 ) – β4 g3 (u)    + 1 λ1 e–μ1 τ1 χ(s, p) – χ s(t – τ1 ), p(t – τ1 )    + 2 λ2 e–μ2 τ2 χ(s, p) – χ s(t – τ2 ), p(t – τ2 )    + 3 λ3 e–μ3 τ3 χ(s, p) – χ s(t – τ3 ), p(t – τ3 )   + 4 Nβ3 g2 (y) + Mβ4 g3 (u) – β5 g4 (p) – qg4 (p)g5 (x)   + 5 rg4 (p)g5 (x) – β6 g5 (x) . Collecting terms of Eq. (33) and using π(s0 ) = 0, we obtain    dV0  χ(s0 , p) = π(s) – π(s0 ) 1 – lim+ p→0 χ(s, p) dt   χ(s0 , p) λχ(s, p) lim – 4 β5 g4 (p) – 5 β6 g5 (x) + g4 (p) p→0+ χ(s, p)     χ(s0 , p) ≤ π(s) – π(s0 ) 1 – lim+ p→0 χ(s, p)

(33)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 9 of 36



 λχ(s, p) χ(s0 , p) + lim+ lim – 4 β5 g4 (p) – 5 β6 g5 (x) p→0 g4 (p) p→0+ χ(s, p)     ∂χ(s0 , 0)/∂p = π(s) – π(s0 ) 1 – ∂χ(s, 0)/∂p   ∂χ(s0 , 0) λ – 1 g4 (p) – 5 β6 g5 (x) +  4 β5 4 β5 g4 (0) ∂p     ∂χ(s0 , 0)/∂p = π(s) – π(s0 ) 1 – + 4 β5 (R0 – 1)g4 (p) – 5 β6 g5 (x). ∂χ(s, 0)/∂p By hypotheses (H1) and (H2), we obtain 

  ∂χ(s0 , 0)/∂p ≤ 0. π(s) – π(s0 ) 1 – ∂χ(s, 0)/∂p 

0 Therefore, if R0 ≤ 1, then dV ≤ 0 for s, p, x > 0. Clearly, dt invariance principle (LIP), we get that 0 is GAS.

dV0 dt

= 0 at 0 . Applying LaSalle’s 

Lemma 3 If R0 > 1 and hypotheses (H1)–(H4) are valid, then sgn(R1 – 1) = sgn(p1 – p2 ) = sgn(s2 – s1 ). Proof Using hypotheses (H1) and (H2), for s1 , s2 , p1 , p2 > 0, we get   (s1 – s2 ) π(s2 ) – π(s1 ) > 0,   (s2 – s1 ) χ(s2 , p2 ) – χ(s1 , p2 ) > 0,   (p2 – p1 ) χ(s1 , p2 ) – χ(s1 , p1 ) > 0,

(34) (35) (36)

and from hypothesis (H4), we obtain 

χ(s1 , p2 ) χ(s1 , p1 ) (p1 – p2 ) – g2 (p2 ) g2 (p1 )

 > 0.

(37)

First, we show that sgn(p1 – p2 ) = sgn(s2 – s1 ). Suppose that sgn(p2 – p1 ) = sgn(s2 – s1 ). Using the steady state conditions of 1 and 2 , we obtain

 π(s2 ) – π(s1 ) = λ χ(s2 , p2 ) – χ(s1 , p1 )   

 = λ χ(s2 , p2 ) – χ(s1 , p2 ) + χ(s1 , p2 ) – χ(s1 , p1 ) . Therefore, from inequalities (34)–(36) we obtain sgn(s2 – s1 ) = sgn(s1 – s2 ), which is a contradiction; hence, sgn(p1 – p2 ) = sgn(s2 – s1 ). Using Eq. (29) and the definition of R1 , we get  R1 – 1 = γ

χ(s2 , p2 ) χ(s1 , p1 ) – g4 (p2 ) g4 (p1 )



  χ(s1 , p2 ) χ(s1 , p1 ) 1  – χ(s2 , p2 ) – χ(s1 , p2 ) + . =γ g4 (p2 ) g4 (p2 ) g4 (p1 ) 

Thus, from Eqs. (35) and (37) we obtain sgn(R1 – 1) = sgn(p1 – p2 ).



Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 10 of 36

Theorem 2 If R1 ≤ 1 < R0 and hypotheses (H1)–(H4) are valid, then 1 is GAS. Proof Let    w χ(s1 , p1 ) g1 (w1 ) dη + 1 w – w1 – dη s1 χ(η, p1 ) w1 g1 (η)      y  u g2 (y1 ) g3 (u1 ) dη + 3 u – u1 – dη + 2 y – y1 – y1 g2 (η) u1 g3 (η)   τ1  χ(s(t – θ ), p(t – θ )) –μ1 τ1 dθ F + 1 λ1 χ(s1 , p1 )e χ(s1 , p1 ) 0   τ2  χ(s(t – θ ), p(t – θ )) F dθ + 2 λ2 χ(s1 , p1 )e–μ2 τ2 χ(s1 , p1 ) 0   τ3  χ(s(t – θ ), p(t – θ )) –μ3 τ3 dθ F + 3 λ3 χ(s1 , p1 )e χ(s1 , p1 ) 0    p g4 (p1 ) + 4 p – p1 – dη + 5 x. p1 g4 (η) 

s

V1 = s – s1 –

Calculating

dV1 dt

along the solutions of (12)–(17), we obtain

   χ(s1 , p1 )  dV1 = 1– π(s) – λχ(s, p) dt χ(s, p1 )     g1 (w1 )  –μ1 τ1  λ1 e χ s(t – τ1 ), p(t – τ1 ) – (a1 + β2 )g1 (w) + 1 1 – g1 (w)     g2 (y1 )  –μ2 τ2  λ2 e χ s(t – τ2 ), p(t – τ2 ) + a1 g1 (w) – β3 g2 (y) + 2 1 – g2 (y)     g3 (u1 )  –μ3 τ3  + 3 1 – χ s(t – τ3 ), p(t – τ3 ) – β4 g3 (u) λ3 e g3 (u)    + 1 λ1 e–μ1 τ1 χ(s, p) – χ s(t – τ1 ), p(t – τ1 )   χ(s(t – τ1 ), p(t – τ1 )) –μ1 τ1 ln + 1 λ1 χ(s1 , p1 )e χ(s, p)    –μ2 τ2 χ(s, p) – χ s(t – τ2 ), p(t – τ2 ) + 2 λ 2 e   χ(s(t – τ2 ), p(t – τ2 )) –μ2 τ2 + 2 λ2 χ(s1 , p1 )e ln χ(s, p)    –μ3 τ3 χ(s, p) – χ s(t – τ3 ), p(t – τ3 ) + 3 λ 3 e   χ(s(t – τ3 ), p(t – τ3 )) –μ3 τ3 ln + 3 λ3 χ(s1 , p1 )e χ(s, p)    g4 (p1 )  + 4 1 – Nβ3 g2 (y) + Mβ4 g3 (u) – β5 g4 (p) – qg4 (p)g5 (x) g4 (p)   + 5 rg4 (p)g5 (x) – β6 g5 (x) . Collecting terms of Eq. (38) and applying the conditions of the steady state 1 π(s1 ) = λχ(s1 , p1 ),

(38)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 11 of 36

(a1 + β2 )g1 (w1 ) = λ1 e–μ1 τ1 χ(s1 , p1 ),   2 β3 g2 (y1 ) = 1 λ1 e–μ1 τ1 + 2 λ2 e–μ2 τ2 χ(s1 , p1 ), β4 g3 (u1 ) = λ3 e–μ3 τ3 χ(s1 , p1 ),

4 β5 g4 (p1 ) = λχ(s1 , p1 ),

we get      χ(s1 , p1 ) dV1  χ(s1 , p1 ) = π(s) – π(s1 ) 1 – + λχ(s1 , p1 ) 1 – dt χ(s, p1 ) χ(s, p1 )   g4 (p) χ(s, p) – + λχ(s1 , p1 ) χ(s, p1 ) g4 (p1 ) – 1 λ1 e–μ1 τ1 χ(s1 , p1 )

g1 (w1 )χ(s(t – τ1 ), p(t – τ1 )) + 1 λ1 e–μ1 τ1 χ(s1 , p1 ) g1 (w)χ(s1 , p1 )

– 2 λ2 e–μ2 τ2 χ(s1 , p1 )

g2 (y1 )χ(s(t – τ2 ), p(t – τ2 )) g2 (y)χ(s1 , p1 )

g2 (y1 )g1 (w) g2 (y)g1 (w1 )   + 1 λ1 e–μ1 τ1 + 2 λ2 e–μ2 τ2 χ(s1 , p1 )

– 1 λ1 e–μ1 τ1 χ(s1 , p1 )

g3 (u1 )χ(s(t – τ3 ), p(t – τ3 )) χ(s1 , p1 )g3 (u)   χ(s(t – τ1 ), p(t – τ1 )) –μ3 τ3 –μ1 τ1 χ(s1 , p1 ) + 1 λ1 e χ(s1 , p1 ) ln + 3 λ 3 e χ(s, p)   χ(s(t – τ2 ), p(t – τ2 )) + 2 λ2 e–μ2 τ2 χ(s1 , p1 ) ln χ(s, p)   χ(s(t – τ3 ), p(t – τ3 )) –μ3 τ3 χ(s1 , p1 ) ln + 3 λ 3 e χ(s, p) – 3 λ3 e–μ3 τ3 χ(s1 , p1 )

  g2 (y)g4 (p1 ) – 1 λ1 e–μ1 τ1 + 2 λ2 e–μ2 τ2 χ(s1 , p1 ) g2 (y1 )g4 (p) g3 (u)g4 (p1 ) g3 (u1 )g4 (p)   β6 g5 (x). + λχ(s1 , p1 ) + r5 g4 (p1 ) – r

– 3 λ3 e–μ3 τ3 χ(s1 , p1 )

(39)

ˆ = w1 , yˆ = y1 , and pˆ = p1 , we can obtain Using inequalities (30) with sˆ = s1 , w       dV1  g4 (p) χ(s, p1 ) χ(s1 , p1 ) χ(s, p) = π(s) – π(s1 ) 1 – + λχ(s1 , p1 ) – 1– dt χ(s, p1 ) χ(s, p1 ) g4 (p1 ) χ(s, p)      g4 (p)χ(s, p1 ) χ(s1 , p1 ) – λχ(s1 , p1 ) F +F χ(s, p1 ) g4 (p1 )χ(s, p)      g2 (y1 )g1 (w) g1 (w1 )χ(s(t – τ1 ), p(t – τ1 )) +F – 1 λ1 e–μ1 τ1 χ(s1 , p1 ) F g1 (w)χ(s1 , p1 ) g2 (y)g1 (w1 )   g2 (y1 )χ(s(t – τ2 ), p(t – τ2 )) – 2 λ2 e–μ2 τ2 χ(s1 , p1 )F g2 (y)χ(s1 , p1 )      g3 (u)g4 (p1 ) g3 (u1 )χ(s(t – τ3 ), p(t – τ3 )) +F – 3 λ3 e–μ3 τ3 χ(s1 , p1 ) F g3 (u)χ(s1 , p1 ) g3 (u1 )g4 (p)

Elaiw et al. Advances in Difference Equations (2018) 2018:85



–μ1 τ1

– 1 λ 1 e

–μ2 τ2

+ 2 λ 2 e

Page 12 of 36



  + r5 g4 (p1 ) – g4 (p2 ) g5 (x).



g2 (y)g4 (p1 ) χ(s1 , p1 )F g2 (y1 )g4 (p)



Hypotheses (H1), (H2), (H4), Lemma 3, and the condition R1 ≤ 1 imply that   χ(s1 , p1 ) π(s) – π(s1 ) 1 – ≤ 0, χ(s, p1 )    g4 (p) χ(s, p1 ) χ(s, p) – 1– ≤ 0, χ(s, p1 ) g4 (p1 ) χ(s, p)





g4 (p1 ) – g4 (p2 ) ≤ 0. It follows that, for all s, y, p, x > 0, we have

dV1 dt

≤ 0 and

dV1 dt

= 0 at 1 . By LIP 1 is GAS. 

Theorem 3 If R1 > 1 and hypotheses (H1)–(H4) are valid, then 2 is GAS. Proof Define    w χ(s2 , p2 ) g1 (w2 ) dη + 1 w – w2 – dη V2 = s – s2 – s2 χ(η, p2 ) w2 g1 (η)      y  u g2 (y2 ) g3 (u2 ) dη + 3 u – u2 – dη + 2 y – y2 – y2 g2 (η) u2 g3 (η)   τ1  χ(s(t – θ ), p(t – θ )) –μ1 τ1 dθ F + 1 λ1 χ(s2 , p2 )e χ(s2 , p2 ) 0   τ2  χ(s(t – θ ), p(t – θ )) F dθ + 2 λ2 χ(s2 , p2 )e–μ2 τ2 χ(s2 , p2 ) 0   τ3  χ(s(t – θ ), p(t – θ )) –μ3 τ3 dθ F + 3 λ3 χ(s2 , p2 )e χ(s2 , p2 ) 0      p  x g4 (p2 ) g5 (x2 ) dη + 5 x – x2 – dη . + 4 p – p2 – p2 g4 (η) x2 g5 (η) 

Calculating

dV2 dt

s

along the solutions of model (12)-(17), we get

   χ(s2 , p2 )  dV2 = 1– π(s) – λχ(s, p) dt χ(s, p2 )     g1 (w2 )  –μ1 τ1  λ1 e χ s(t – τ1 ), p(t – τ1 ) – (a1 + β2 )g1 (w) + 1 1 – g1 (w)     g2 (y2 )  –μ2 τ2  λ2 e χ s(t – τ2 ), p(t – τ2 ) + a1 g1 (w) – β3 g2 (y) + 2 1 – g2 (y)     g3 (u2 )  –μ3 τ3  + 3 1 – χ s(t – τ3 ), p(t – τ3 ) – β4 g3 (u) λ3 e g3 (u)    + 1 λ1 e–μ1 τ1 χ(s, p) – χ s(t – τ1 ), p(t – τ1 )   χ(s(t – τ1 ), p(t – τ1 )) + 1 λ1 e–μ1 τ1 χ(s2 , p2 ) ln χ(s, p)    + 2 λ2 e–μ2 τ2 χ(s, p) – χ s(t – τ2 ), p(t – τ2 )

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 13 of 36

  χ(s(t – τ2 ), p(t – τ2 )) + 2 λ 2 e χ(s2 , p2 ) ln χ(s, p)    + 3 λ3 e–μ3 τ3 χ(s, p) – χ s(t – τ3 ), p(t – τ3 )   χ(s(t – τ3 ), p(t – τ3 )) –μ3 τ3 + 3 λ 3 e χ(s2 , p2 ) ln χ(s, p)    g4 (p2 )  Nβ3 g2 (y) + Mβ4 g3 (u) – β5 g4 (p) – qg4 (p)g5 (x) + 4 1 – g4 (p)    g5 (x2 )  rg4 (p)g5 (x) – β6 g5 (x) . + 5 1 – g5 (x) –μ2 τ2

(40)

Collecting terms of Eq. (40) and using the steady state conditions for 2 : π(s2 ) = λχ(s2 , p2 ), (a1 + β2 )g1 (w2 ) = λ1 e–μ1 τ1 χ(s2 , p2 ),   2 β3 g2 (y2 ) = 1 λ1 e–μ1 τ1 + 2 λ2 e–μ2 τ2 χ(s2 , p2 ), β4 g3 (u2 ) = λ3 e–μ3 τ3 χ(s2 , p2 ), 4 β5 g4 (p) = λχ(s2 , p2 )

4 β5 g4 (p2 ) = λχ(s2 , p2 ) – q4 g4 (p2 )g5 (x2 ),

g4 (p) – q4 g4 (p)g5 (x2 ), g4 (p2 )

we obtain      χ(s2 , p2 ) χ(s2 , p2 ) dV2  = π(s) – π(s2 ) 1 – + λχ(s2 , p2 ) 1 – dt χ(s, p2 ) χ(s, p2 )   g4 (p) χ(s, p) – + λχ(s2 , p2 ) χ(s, p2 ) g4 (p2 ) – 1 λ1 e–μ1 τ1 χ(s2 , p2 )

g1 (w2 )χ(s(t – τ1 ), p(t – τ1 )) + 1 λ1 e–μ1 τ1 χ(s2 , p2 ) g1 (w)χ(s2 , p2 )

– 2 λ2 e–μ2 τ2 χ(s2 , p2 )

g2 (y2 )χ(s(t – τ2 ), p(t – τ2 )) g2 (y)χ(s2 , p2 )

g2 (y2 )g1 (w) g2 (y)g1 (w2 )   + 1 λ1 e–μ1 τ1 + 2 λ2 e–μ2 τ2 χ(s2 , p2 )

– 1 λ1 e–μ1 τ1 χ(s2 , p2 )

g3 (u2 )χ(s(t – τ3 ), p(t – τ3 )) g3 (u)χ(s2 , p2 )   χ(s(t – τ1 ), p(t – τ1 )) + 3 λ3 e–μ3 τ3 χ(s2 , p2 ) + 1 λ1 e–μ1 τ1 χ(s2 , p2 ) ln χ(s, p)   χ(s(t – τ2 ), p(t – τ2 )) + 2 λ2 e–μ2 τ2 χ(s2 , p2 ) ln χ(s, p)   χ(s(t – τ3 ), p(t – τ3 )) + 3 λ3 e–μ3 τ3 χ(s2 , p2 ) ln χ(s, p)

– 3 λ3 e–μ3 τ3 χ(s2 , p2 )

  g2 (y)g4 (p2 ) – 1 λ1 e–μ1 τ1 + 2 λ2 e–μ2 τ2 χ(s2 , p2 ) g2 (y2 )g4 (p) – 3 λ3 e–μ3 τ3 χ(s2 , p2 )

g3 (u)g4 (p2 ) + λχ(s2 , p2 ). g3 (u2 )g4 (p)

(41)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 14 of 36

ˆ = w2 , yˆ = y2 , and pˆ = p2 , we can get By inequalities (30) with sˆ = s2 , w       dV2  g4 (p) χ(s, p2 ) χ(s2 , p2 ) χ(s, p) = π(s) – π(s2 ) 1 – + λχ(s2 , p2 ) – 1– dt χ(s, p2 ) χ(s, p2 ) g4 (p2 ) χ(s, p)      χ(s2 , p2 ) g4 (p)χ(s, p2 ) +F – λχ(s2 , p2 ) F χ(s, p2 ) g4 (p2 )χ(s, p)      g2 (y2 )g1 (w) g1 (w2 )χ(s(t – τ1 ), p(t – τ1 )) –μ1 τ1 +F χ(s2 , p2 ) F – 1 λ 1 e g1 (w)χ(s2 , p2 ) g2 (y)g1 (w2 )   g2 (y2 )χ(s(t – τ2 ), p(t – τ2 )) – 2 λ2 e–μ2 τ2 χ(s2 , p2 )F g2 (y)χ(s2 , p2 )      g3 (u2 )χ(s(t – τ3 ), p(t – τ3 )) g3 (u)g4 (p2 ) –μ3 τ3 – 3 λ 3 e χ(s2 , p2 ) F +F g3 (u)χ(s2 , p2 ) g3 (u2 )g4 (p)     g2 (y)g4 (p2 ) . – 1 λ1 e–μ1 τ1 + 2 λ2 e–μ2 τ2 χ(s2 , p2 )F g2 (y2 )g4 (p) According to hypotheses (H1), (H2), and (H4), we get implies that 2 is GAS.

dV2 dt

≤ 0 and

dV2 dt

= 0 at 2 . LIP 

3 Model with delay-distributed In the next model, we consider a general delay-distributed HIV infection model with humoral immunity as follows:     s˙ = π s(t) – λχ s(t), p(t) ,  h1     ˙ = λ1 w f1 (τ )e–μ1 τ χ s(t – τ ), p(t – τ ) dτ – (a1 + β2 )g1 w(t) ,

(42) (43)

0



h2

y˙ = λ2

      f2 (τ )e–μ2 τ χ s(t – τ ), p(t – τ ) dτ + a1 g1 w(t) – β3 g2 y(t) ,

(44)

    f3 (τ )e–μ3 τ χ s(t – τ ), p(t – τ ) dτ – β4 g3 u(t) ,

(45)

0

 u˙ = λ3

h3

0

          p˙ = Nβ3 g2 y(t) + Mβ4 g3 u(t) – β5 g4 p(t) – qg4 p(t) g5 x(t) ,       x˙ = rg4 p(t) g5 x(t) – β6 g5 x(t) ,

(46) (47)

where fi (τ )e–μi τ over the time interval [0, hi ], i = 1, 2, 3, represent the probabilities that uninfected cells contacted by HIV at time t – τ survived τ time units and became infected at time t. The probability distribution function fi (τ ) is assumed to satisfy fi (τ ) > 0 and 



hi

hi

fi (τ ) dτ = 1, 0

fi (η)eυη dη < ∞,

i = 1, 2, 3,

0

where υ is a positive constant. Let us denote i (τ ) = fi (τ )e–μi τ and Fi = 0 < Fi ≤ 1, i = 1, 2, 3.

 hi 0

i (τ ) dτ ; thus

Lemma 4 Let hypotheses (H1)–(H3) be valid, then the solutions of system (42)–(47) are nonnegative and ultimately bounded.

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 15 of 36

Proof The nonnegativity of solutions of system (42)–(47) can easily be shown as given in Lemma 1. From Eq. (42) we have that lim supt→∞ s(t) ≤ bb¯ = M1 . Let 



h1 0



h2

1 (τ )s(t – τ ) dτ + N

G(t) = N

h3

2 (τ )s(t – τ ) dτ + M 0

3 (τ )s(t – τ ) dτ 0

q 1 + Nw(t) + Ny(t) + Mu(t) + p(t) + x(t), 2 2r then ˙ =N G(t)



h1 0



    1 (τ ) π s(t – τ ) – λχ s(t – τ ), p(t – τ ) dτ h2

    2 (τ ) π s(t – τ ) – λχ s(t – τ ), p(t – τ ) dτ

h3

    3 (τ ) π s(t – τ ) – λχ s(t – τ ), p(t – τ ) dτ

+N 0

 +M

0

  + N λ1   + N λ2

    1 (τ )χ s(t – τ ), p(t – τ ) dτ – (a1 + β2 )g1 w(t)

h2

      2 (τ )χ s(t – τ ), p(t – τ ) dτ + a1 g1 w(t) – β3 g2 y(t)

0

0





h1



h3

+ M λ3



     3 (τ )χ s(t – τ ), p(t – τ ) dτ – β4 g3 u(t)

0

          1 + Nβ3 g2 y(t) + Mβ4 g3 u(t) – β5 g4 p(t) – qg4 p(t) g5 x(t) 2      q  rg4 p(t) g5 x(t) – β6 g5 x(t) + 2r  h1  h2



 ≤N 1 (τ ) b – bs(t – τ ) dτ + N 2 (τ ) b – bs(t – τ ) dτ 0



0 h3

+M

 3 (τ ) b – bs(t – τ ) dτ

0

1 1 1 q – Nβ2 α1 w(t) – Nβ3 α2 y(t) – Mβ4 α3 u(t) – β5 α4 p(t) – β6 α5 x(t) 2 2 2 2r   h1  h2 ≤ b(NF1 + NF2 + MF3 ) – σ N 1 (τ )s(t – τ ) dτ + N 2 (τ )s(t – τ ) dτ  +M 0

0 h3

0

q 1 3 (τ )s(t – τ ) dτ + Nw(t) + Ny(t) + Mu(t) + p(t) + x(t) 2 2r

≤ b(2N + M) – σ G(t). Hence, lim supt→∞ G(t) ≤

b(2N+M) σ

and

lim sup w(t) ≤ M2 ,

lim sup y(t) ≤ M2 ,

lim sup p(t) ≤ M4 ,

lim sup x(t) ≤ M5 ,

t→∞

t→∞

t→∞

t→∞

lim sup u(t) ≤ M3 , t→∞



Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 16 of 36

where σ , M1 , . . . , M5 are given in the previous section. Therefore, s(t), w(t), y(t), u(t), p(t), and x(t) are ultimately bounded. Moreover, the set  defined by (19) is also positively invariant with respect to system (42)–(47). 

3.1 Steady states Lemma 5 Suppose that hypotheses (H1)–(H4) are valid (42)–(47), then there exist two bifurcation parameters R0 and R1 with R0 > R1 > 0 such that (i) if R0 ≤ 1, then there exists only one steady state 0 ; (ii) if R1 ≤ 1 < R0 , then there exist two steady states 0 and 1 ; (iii) if R1 > 1, then there exist three steady states 0 , 1 , and 2 . Proof Let (s, w, y, u, p, x) be any steady state of (42)–(47) satisfying the following equations: 0 = π(s) – λχ(s, p),

(48)

0 = λ1 F1 χ(s, p) – (a1 + β2 )g1 (w),

(49)

0 = λ2 F2 χ(s, p) + a1 g1 (w) – β3 g2 (y),

(50)

0 = λ3 F3 χ(s, p) – β4 g3 (u),

(51)

0 = Nβ3 g2 (y) + Mβ4 g3 (u) – β5 g4 (p) – qg4 (p)g5 (x),

(52)

0 = rg4 (p)g5 (x) – β6 g5 (x).

(53)

From Eq. (53) we obtain two possible solutions: g5 (x) = 0 and g4 (p) = β6 /r. First, we consider the case g5 (x) = 0, then from hypothesis (H3) we have x = 0. Hypothesis (H3) implies that gi–1 , i = 1, . . . , 5, exist, strictly increasing and gi–1 (0) = 0. Let us define 

   λ1 F 1 a1 λ1 F1 + (a1 + β2 )λ2 F2 π(s) , α2 (s) = g2–1 π(s) , λ(a1 + β2 ) β3 λ(a1 + β2 )     λ3 F 3 ζ π(s) , α4 (s) = g4–1 π(s) , α3 (s) = g3–1 β4 λ λ α1 (s) = g1–1

(54)

where ζ=

N(a1 λ1 F1 + (a1 + β2 )λ2 F2 ) + Mλ3 (a1 + β2 )F3 . β5 (a1 + β2 )

From Eqs. (48)–(53) we can get w = α1 (s),

y = α2 (s),

u = α3 (s),

p = α4 (s).

(55)

Obviously, αi (s) > 0 for s ∈ [0, s0 ) and αi (s0 ) = 0, i = 1, . . . , 4. From Eqs. (48), (54), and (55) we obtain     ζ χ s, α4 (s) – g4 α4 (s) = 0.

(56)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 17 of 36

Equation (56) admits a solution s = s0 which gives the infection-free steady state 0 (s0 , 0, 0, 0, 0, 0). Let     1 (s) = ζ χ s, α4 (s) – g4 α4 (s) = 0. By hypotheses (H1) and (H2) we have 1 (0) = –g4 (α4 (0)) < 0 and 1 (s0 ) = 0. Moreover, 1 (s0 ) = ζ



 ∂χ(s0 , 0) ∂χ(s0 , 0) + α4 (s0 ) – g4 (0)α4 (s0 ). ∂s ∂p

From hypothesis (H2) we note that 1 (s0 ) = α4 (s0 )g4 (0)



∂χ(s0 ,0) ∂s

= 0. Then

 ζ ∂χ(s0 , 0) –1 . g4 (0) ∂p

From Eq. (54), we get 1 (s0 ) =

  ζ  ζ ∂χ(s0 , 0) π (s0 )  –1 . λ g4 (0) ∂p

Therefore, using hypothesis (H1), we get π  (s0 ) < 0. So that, if g  ζ(0) ∂χ(s∂p0 ,0) > 1, then 1 (s0 ) < 4 0 and there exists s1 ∈ (0, s0 ) such that 1 (s1 ) = 0. Hypotheses (H1)–(H3) imply that w1 = α1 (s1 ) > 0,

y1 = α2 (s1 ) > 0,

u1 = α3 (s1 ) > 0,

p1 = α4 (s1 ) > 0.

(57)

It means that a humoral-inactivated infection steady state 1 (s1 , w1 , y1 , u1 , p1 , 0) exists when g  ζ(0) ∂χ(s∂p0 ,0) > 1. Now we can define 4

R0 =

ζ ∂χ(s0 , 0) . g4 (0) ∂p

The other solution of Eq. (53) is g4 (p2 ) = βr6 , which yields p2 = g4–1 ( βr6 ) > 0. Substitute p = p2 in Eq. (48) and let 2 (s) = π(s) – λζ (s, p2 ) = 0. According to hypotheses (H1) and (H2), 2 is strictly decreasing, 2 (0) = π(0) > 0, and 2 (s0 ) = –λζ (s0 , p2 ) < 0. Thus, there exists a unique s2 ∈ (0, s0 ) such that 2 (s2 ) = 0. It follows from Eqs. (52) and (55) that  w2 = α1 (s2 ) > 0,  x2 = g5–1

y2 = α2 (s2 ) > 0,

u2 = α3 (s2 ) > 0,

p2 = g4–1

  χ(s2 , p2 ) β5 ζ –1 . q g4 (p2 )

β6 r

 > 0,

2 ,p2 ) Thus, x2 > 0 when ζ χ(s > 1. Let us define the parameter R1 as follows: g4 (p2 )

R1 = ζ

χ(s2 , p2 ) . g4 (p2 )

If R1 > 1, then x2 = g5–1 ( βq5 (R1 –1)) > 0, and there exists a humoral-activated infection steady state 2 (s2 , w2 , y2 , u2 , p2 , x2 ).

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 18 of 36

Clearly, from hypotheses (H2) and (H4), we have R1 = ζ

ζ ∂χ(s2 , 0) ζ ∂χ(s0 , 0) χ(s2 , p2 ) χ(s2 , p) ≤ ζ lim+ =  <  = R0 . p→0 g4 (p2 ) g4 (p) g4 (0) ∂p g4 (0) ∂p



3.2 Global properties Theorem 4 If R0 ≤ 1 and hypotheses (H1)–(H4) hold true for system (42)–(47), then 0 is GAS. Proof Define 

s

χ(s0 , p) dη + κ1 w + κ2 y + κ3 u χ(η, p)  h1  τ   + κ1 λ 1 1 (τ ) χ s(t – θ ), p(t – θ ) dθ dτ

W0 = s – s 0 –

lim

+ s0 p→0

0





h2

0 τ

  χ s(t – θ ), p(t – θ ) dθ dτ

τ

  χ s(t – θ ), p(t – θ ) dθ dτ

2 (τ )

+ κ2 λ 2 0





h3

+ κ3 λ 3

0

3 (τ ) 0

0

+ κ4 p + κ5 x,

(58)

where λ = λ 1 κ1 F 1 + λ 2 κ2 F 2 + λ 3 κ3 F 3 , κ2 = Nκ4 ,

κ3 = Mκ4 ,

(a1 + β2 )κ1 = a1 κ2 , qκ4 = rκ5 .

(59)

The solution of Eqs. (59) is given by κ1 =

a1 Nλ , ζβ5 (a1 + β2 )

We evaluate

dW0 dt

Nλ , ζβ5

κ2 =

κ3 =

Mλ , ζβ5

κ4 =

λ , ζβ5

κ5 =

qλ . rζβ5

(60)

along the solutions of (42)-(47) as follows:

   h1    dW0 χ(s0 , p)  1 (τ )χ s(t – τ ), p(t – τ ) dτ = 1 – lim+ π(s) – λχ(s, p) + κ1 λ1 p→0 χ(s, p) dt 0  h2   – κ1 (a1 + β2 )g1 (w) + κ2 λ2 2 (τ )χ s(t – τ ), p(t – τ ) dτ + a1 κ2 g1 (w)  – κ2 β3 g2 (y) + κ3 λ3 

0 h3

  3 (τ )χ s(t – τ ), p(t – τ ) dτ – κ3 β4 g3 (u)

0 h1

   1 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) dτ

h2

   2 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) dτ

h3

   3 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) dτ

+ κ1 λ 1 0

 + κ2 λ 2

0

 + κ3 λ 3

0

Elaiw et al. Advances in Difference Equations (2018) 2018:85

  + κ4 Nβ3 g2 (y) + Mβ4 g3 (u) – β5 g4 (p) – qg4 (p)g5 (x)   + κ5 rg4 (p)g5 (x) – β6 g5 (x) .

Page 19 of 36

(61)

Collecting terms of Eq. (61) and using π(s0 ) = 0, we obtain    χ(s0 , p) dW0  = π(s) – π(s0 ) 1 – lim+ p→0 χ(s, p) dt   χ(s0 , p) λχ(s, p) lim+ – κ4 β5 g4 (p) – κ5 β6 g5 (x) + g4 (p) p→0 χ(s, p)     χ(s0 , p) ≤ π(s) – π(s0 ) 1 – lim+ p→0 χ(s, p)   λχ(s, p) χ(s0 , p) lim+ – κ4 β5 g4 (p) – κ5 β6 g5 (x) + lim+ p→0 g4 (p) p→0 χ(s, p)     ∂χ(s0 , 0)/∂p = π(s) – π(s0 ) 1 – ∂χ(s, 0)/∂p   λ ∂χ(s0 , 0) + κ 4 β5 – 1 g4 (p) – κ5 β6 g5 (x) κ4 β5 g4 (0) ∂p     ∂χ(s0 , 0)/∂p + κ4 β5 (R0 – 1)g4 (p) – κ5 β6 g5 (x). = π(s) – π(s0 ) 1 – ∂χ(s, 0)/∂p

(62)

From hypotheses (H1) and (H2), we have 

  ∂χ(s0 , 0)/∂p ≤ 0. π(s) – π(s0 ) 1 – ∂χ(s, 0)/∂p 

Therefore, if R0 ≤ 1, then

dW0 dt

≤ 0. Thus, 0 is GAS.



Theorem 5 If R1 ≤ 1 < R0 and hypotheses (H1)–(H4) are valid for (42)–(47), then 1 is GAS. Proof We introduce    w χ(s1 , p1 ) g1 (w1 ) dη + κ1 w – w1 – dη W1 = s – s 1 – s1 χ(η, p1 ) w1 g1 (η)      y  u g2 (y1 ) g3 (u1 ) dη + κ3 u – u1 – dη + κ2 y – y1 – y1 g2 (η) u1 g3 (η)   h1  τ  χ(s(t – θ ), p(t – θ )) dθ dτ 1 (τ ) F + κ1 λ1 χ(s1 , p1 ) χ(s1 , p1 ) 0 0   h2  τ  χ(s(t – θ ), p(t – θ )) dθ dτ + κ2 λ2 χ(s1 , p1 ) 2 (τ ) F χ(s1 , p1 ) 0 0   h3  τ  χ(s(t – θ ), p(t – θ )) dθ dτ 3 (τ ) F + κ3 λ3 χ(s1 , p1 ) χ(s1 , p1 ) 0 0    p g4 (p1 ) dη + κ5 x. + κ4 p – p1 – p1 g4 (η) 

s

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Evaluating

dW1 dt

Page 20 of 36

along the trajectories of (42)–(47), we have

   χ(s1 , p1 )  dW1 = 1– π(s) – λχ(s, p) dt χ(s, p1 )    h1    g1 (w1 ) λ1 1 (τ )χ s(t – τ ), p(t – τ ) dτ – (a1 + β2 )g1 (w) + κ1 1 – g1 (w) 0    h2    g2 (y1 ) λ2 2 (τ )χ s(t – τ ), p(t – τ ) dτ + a1 g1 (w) – β3 g2 (y) + κ2 1 – g2 (y) 0     h3   g3 (u1 ) 3 (τ )χ s(t – τ ), p(t – τ ) dτ – β4 g3 (u) λ3 + κ3 1 – g3 (u) 0   h1   1 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) + χ(s1 , p1 ) + κ1 λ 1 0

  χ(s(t – τ ), p(t – τ )) dτ × ln χ(s, p)   h2   2 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) + χ(s1 , p1 ) + κ2 λ 2 0

  χ(s(t – τ ), p(t – τ )) dτ × ln χ(s, p)   h3   3 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) + χ(s1 , p1 ) + κ3 λ 3 0

  χ(s(t – τ ), p(t – τ )) dτ × ln χ(s, p)    g4 (p1 )  Nβ3 g2 (y) + Mβ4 g3 (u) – β5 g4 (p) – qg4 (p)g5 (x) + κ4 1 – g4 (p)   + κ5 rg4 (p)g5 (x) – β6 g5 (x) . Collecting terms of Eq. (63) and applying the conditions of the steady state 1 : π(s1 ) = λχ(s1 , p1 ), (a1 + β2 )g1 (w1 ) = λ1 F1 χ(s1 , p1 ), β4 g3 (u1 ) = λ3 F3 χ(s1 , p1 ),

κ2 β3 g2 (y1 ) = (κ1 λ1 F1 + κ2 λ2 F2 )χ(s1 , p1 ),

κ4 β5 g4 (p1 ) = λχ(s1 , p1 ),

we get    χ(s1 , p1 ) dW1  = π(s) – π(s1 ) 1 – dt χ(s, p1 )     χ(s1 , p1 ) χ(s, p) g4 (p) + λχ(s1 , p1 ) 1 – + λχ(s1 , p1 ) – χ(s, p1 ) χ(s, p1 ) g4 (p1 )  h1 χ(s(t – τ ), p(t – τ ))g1 (w1 ) dτ + κ1 λ1 F1 χ(s1 , p1 ) 1 (τ ) – κ1 λ1 χ(s1 , p1 ) χ(s1 , p1 )g1 (w) 0  h2 χ(s(t – τ ), p(t – τ ))g2 (y1 ) dτ 2 (τ ) – κ2 λ2 χ(s1 , p1 ) χ(s1 , p1 )g2 (y) 0

(63)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

– κ1 λ1 F1 χ(s1 , p1 ) 

Page 21 of 36

g2 (y1 )g1 (w) + (κ1 λ1 F1 + κ2 λ2 F2 )χ(s1 , p1 ) g2 (y)g1 (w1 ) h3

χ(s(t – τ ), p(t – τ )g3 (u1 )) dτ + κ3 λ3 F3 χ(s1 , p1 ) χ(s1 , p1 )g3 (u) 0    h1 χ(s(t – τ ), p(t – τ )) + κ1 λ1 χ(s1 , p1 ) 1 (τ ) ln dτ χ(s, p) 0    h2 χ(s(t – τ ), p(t – τ )) dτ 2 (τ ) ln + κ2 λ2 χ(s1 , p1 ) χ(s, p) 0    h3 χ(s(t – τ ), p(t – τ )) dτ + κ3 λ3 χ(s1 , p1 ) 3 (τ ) ln χ(s, p) 0 – κ3 λ3 χ(s1 , p1 )

3 (τ )

g2 (y)g4 (p1 ) g3 (u)g4 (p1 ) – κ3 λ3 F3 χ(s1 , p1 ) g2 (y1 )g4 (p) g3 (u1 )g4 (p)   β6 g5 (x). + λχ(s1 , p1 ) + rκ5 g4 (p1 ) – r – (κ1 λ1 F1 + κ2 λ2 F2 )χ(s1 , p1 )

ˆ = w1 , yˆ = y1 , pˆ = p1 , and τ1 = τ2 = τ3 = τ , we can obtain Using inequalities (30) with sˆ = s1 , w    χ(s1 , p1 ) dW1  = π(s) – π(s1 ) 1 – dt χ(s, p1 )    g4 (p) χ(s, p1 ) χ(s, p) – 1– + λχ(s1 , p1 ) χ(s, p1 ) g4 (p1 ) χ(s, p)      g4 (p)χ(s, p1 ) χ(s1 , p1 ) – λχ(s1 , p1 ) F +F χ(s, p1 ) g4 (p1 )χ(s, p)    h1 g1 (w1 )χ(s(t – τ ), p(t – τ )) 1 (τ )F dτ – κ1 λ1 χ(s1 , p1 ) g1 (w)χ(s1 , p1 ) 0    h2 g2 (y1 )χ(s(t – τ ), p(t – τ )) dτ 2 (τ )F – κ2 λ2 χ(s1 , p1 ) g2 (y)χ(s1 , p1 ) 0    h3 g3 (u1 )χ(s(t – τ ), p(t – τ )) – κ3 λ3 χ(s1 , p1 ) 3 (τ )F g3 (u)χ(s1 , p1 ) 0   g2 (y1 )g1 (w) – κ1 λ1 F1 χ(s1 , p1 )F g2 (y)g1 (w1 )   g2 (y)g4 (p1 ) – (κ1 λ1 F1 + κ2 λ2 F2 )χ(s1 , p1 )F g2 (y1 )g4 (p)     g3 (u)g4 (p1 ) + rκ5 g4 (p1 ) – g4 (p2 ) g5 (x). – κ3 λ3 F3 χ(s1 , p1 )F g3 (u1 )g4 (p) Using hypotheses (H1), (H2), (H4), and Lemma 3, we get occurs at 1 . By LIP, 1 is GAS.

dW1 dt

≤ 0, where the equality 

Theorem 6 If R1 > 1 and hypotheses (H1)–(H4) are valid for (42)–(47), then 2 is GAS. Proof Define 

s

W2 = s – s 2 – s2

   w χ(s2 , p2 ) g1 (w2 ) dη + κ1 w – w2 – dη χ(η, p2 ) w2 g1 (η)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 22 of 36

    u g2 (y2 ) g3 (u2 ) dη + κ3 u – u2 – dη y2 g2 (η) u2 g3 (η)   h1  τ  χ(s(t – θ ), p(t – θ )) 1 (τ ) F + κ1 λ1 χ(s2 , p2 ) dθ dτ χ(s2 , p2 ) 0 0   h2  τ  χ(s(t – θ ), p(t – θ )) dθ dτ 2 (τ ) F + κ2 λ2 χ(s2 , p2 ) χ(s2 , p2 ) 0 0   h3  τ  χ(s(t – θ ), p(t – θ )) dθ dτ 3 (τ ) F + κ3 λ3 χ(s2 , p2 ) χ(s2 , p2 ) 0 0      p  x g4 (p2 ) g5 (x2 ) dη + κ5 x – x2 – dη . + κ4 p – p2 – p2 g4 (η) x2 g5 (η)   + κ2 y – y2 –

Calculating

dW2 dt

y

along the solutions of model (42)–(47), we get

   χ(s2 , p2 )  dW2 = 1– π(s) – λχ(s, p) dt χ(s, p2 )   h1     g1 (w2 ) λ1 1 (τ )χ s(t – τ ), p(t – τ ) dτ – (a1 + β2 )g1 (w) + κ1 1 – g1 (w) 0   h2     g2 (y2 ) λ2 2 (τ )χ s(t – τ ), p(t – τ ) dτ + a1 g1 (w) – β3 g2 (y) + κ2 1 – g2 (y) 0     h3   g3 (u2 ) 3 (τ )χ s(t – τ ), p(t – τ ) dτ – β4 g3 (u) λ3 + κ3 1 – g3 (u) 0   h1   1 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) + χ(s2 , p2 ) + κ1 λ 1 0

  χ(s(t – τ ), p(t – τ )) dτ × ln χ(s, p)   h2   2 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) + χ(s2 , p2 ) + κ2 λ 2 0

  χ(s(t – τ ), p(t – τ )) dτ × ln χ(s, p)   h3   3 (τ ) χ(s, p) – χ s(t – τ ), p(t – τ ) + χ(s2 , p2 ) + κ3 λ 3 0

  χ(s(t – τ ), p(t – τ )) dτ × ln χ(s, p)    g4 (p2 )  + κ4 1 – Nβ3 g2 (y) + Mβ4 g3 (u) – β5 g4 (p) – qg4 (p)g5 (x) g4 (p)    g5 (x2 )  + κ5 1 – rg4 (p)g5 (x) – β6 g5 (x) . g5 (x) Collecting terms of Eq. (64) and applying the steady state conditions for 2 : π(s2 ) = λχ(s2 , p2 ), (a1 + β2 )g1 (w2 ) = λ1 F1 χ(s2 , p2 ), β4 g3 (u2 ) = λ3 F3 χ(s2 , p2 ),

κ2 β3 g2 (y2 ) = (κ1 λ1 F1 + κ2 λ2 F2 )χ(s2 , p2 ),

κ4 β5 g4 (p2 ) = λχ(s2 , p2 ) – qκ4 g4 (p2 )g5 (x2 ),

(64)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

κ4 β5 g4 (p) = λχ(s2 , p2 )

Page 23 of 36

g4 (p) – qκ4 g4 (p)g5 (x2 ) g4 (p2 )

we obtain      χ(s2 , p2 ) χ(s2 , p2 ) dW2  = π(s) – π(s2 ) 1 – + λχ(s2 , p2 ) 1 – dt χ(s, p2 ) χ(s, p2 )   χ(s, p) g4 (p) + λχ(s2 , p2 ) – χ(s, p2 ) g4 (p2 )  h1 g1 (w2 )χ(s(t – τ ), p(t – τ )) dτ 1 (τ ) – κ1 λ1 χ(s2 , p2 ) g1 (w)χ(s2 , p2 ) 0  h2 g2 (y2 )χ(s(t – τ ), p(t – τ )) dτ 2 (τ ) + κ1 λ1 F1 χ(s2 , p2 ) – κ2 λ2 χ(s2 , p2 ) g2 (y)χ(s2 , p2 ) 0 – κ1 λ1 F1 χ(s2 , p2 ) 

g2 (y2 )g1 (w) + (κ1 λ1 F1 + κ2 λ2 F2 )χ(s2 , p2 ) g2 (y)g1 (w2 ) h3

g3 (u2 )χ(s(t – τ ), p(t – τ )) dτ + κ3 λ3 F3 χ(s2 , p2 ) g3 (u)χ(s2 , p2 ) 0    h1 χ(s(t – τ ), p(t – τ )) dτ 1 (τ ) ln + κ1 λ1 χ(s2 , p2 ) χ(s, p) 0    h2 χ(s(t – τ ), p(t – τ )) + κ2 λ2 χ(s2 , p2 ) 2 (τ ) ln dτ χ(s, p) 0    h3 χ(s(t – τ ), p(t – τ )) dτ 3 (τ ) ln + κ3 λ3 χ(s2 , p2 ) χ(s, p) 0 – κ3 λ3 χ(s2 , p2 )

3 (τ )

– (κ1 λ1 F1 + κ2 λ2 F2 )χ(s2 , p2 ) – κ3 λ3 F3 χ(s2 , p2 )

g2 (y)g4 (p2 ) g2 (y2 )g4 (p)

g3 (u)g4 (p2 ) + λχ(s2 , p2 ). g3 (u2 )g4 (p)

(65)

ˆ = w2 , yˆ = y2 , pˆ = p2 , and τ1 = τ2 = τ3 = τ , we can obtain Using inequalities (30) with sˆ = s2 , w    χ(s2 , p2 ) dW2  = π(s) – π(s2 ) 1 – dt χ(s, p2 )    g4 (p) χ(s, p2 ) χ(s, p) – 1– + λχ(s2 , p2 ) χ(s, p2 ) g4 (p2 ) χ(s, p)      χ(s2 , p2 ) g4 (p)χ(s, p2 ) +F – λχ(s2 , p2 ) F χ(s, p2 ) g4 (p2 )χ(s, p)    h1 g1 (w2 )χ(s(t – τ ), p(t – τ )) dτ 1 (τ )F – κ1 λ1 χ(s2 , p2 ) g1 (w)χ(s2 , p2 ) 0    h2 g2 (y2 )χ(s(t – τ ), p(t – τ )) dτ 2 (τ )F – κ2 λ2 χ(s2 , p2 ) g2 (y)χ(s2 , p2 ) 0    h3 g3 (u2 )χ(s(t – τ ), p(t – τ )) dτ – κ3 λ3 χ(s2 , p2 ) 3 (τ )F g3 (u)χ(s2 , p2 ) 0   g2 (y2 )g1 (w) – κ1 λ1 F1 χ(s2 , p2 )F g2 (y)g1 (w2 )

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 24 of 36



g2 (y)g4 (p2 ) – (κ1 λ1 F1 + κ2 λ2 F2 )χ(s2 , p2 )F g2 (y2 )g4 (p)   g3 (u)g4 (p2 ) – κ3 λ3 F3 χ(s2 , p2 )F . g3 (u2 )g4 (p) According to hypotheses (H1), (H2), and (H4), we get that 2 is GAS.



dW2 dt

≤ 0. Applying LIP, one can show 

4 Numerical simulations We now perform some computer simulations on the following application:   (1 – ε1 )λs(t)p(t) s(t) s˙(t) = ρ – β1 s(t) + ωs(t) 1 – – , smax 1 + ηp(t)

(66)

˙ = w(t)

(1 – ε1 )λ1 e–μ1 τ1 s(t – τ1 )p(t – τ1 ) – (a1 + β2 )w(t), 1 + ηp(t – τ1 )

(67)

y˙ (t) =

(1 – ε1 )λ2 e–μ2 τ2 s(t – τ2 )p(t – τ2 ) + a1 w(t) – β3 y(t), 1 + ηp(t – τ2 )

(68)

˙ = u(t)

(1 – ε1 )λ3 e–μ3 τ3 s(t – τ3 )p(t – τ3 ) – β4 u(t), 1 + ηp(t – τ3 )

(69)

p˙ (t) = (1 – ε2 )Nβ3 y(t) + (1 – ε2 )Mβ4 u(t) – β5 p(t) – qp(t)x(t),

(70)

x˙ (t) = rp(t)x(t) – β6 x(t).

(71)

We assume that ω < β1 . In this application, we consider the following specific forms of the general functions:   s(t) π s(t) = ρ – β1 s(t) + ωs(t) 1 – , smax 



gi (θ ) = θ ,

  s(t)p(t) , χ s(t), p(t) = 1 + ηp(t)

i = 1, . . . , 5.

First we verify hypotheses (H1)–(H4) for the chosen forms, then we solve the system using MATLAB. Clearly, π(0) = ρ > 0 and π(s0 ) = 0, where    smax 4ρω . s0 = ω – β1 + (ω – β1 )2 + 2ω smax We have π  (s) = –β1 + ω –

2ωs < 0. smax

(72)

Clearly, π(s) > 0 for s ∈ [0, s0 ) and

π(s) = ρ – (β1 – ω)s – ω

s2 smax

≤ ρ – (β1 – ω)s.

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 25 of 36

Then hypothesis (H1) is satisfied. We also have χ(s, p) > 0, χ(0, p) = χ(s, 0) = 0 for s, p ∈ (0, ∞), and p ∂χ(s, p) = , ∂s 1 + ηp

∂χ(s, p) s = , ∂p (1 + ηp)2

Then ∂χ(s,p) > 0, ∂χ(s,p) > 0, and ∂s ∂p satisfied. In addition,

∂χ(s,0) ∂p

∂χ(s, 0) = s. ∂p

> 0 for s, p ∈ (0, ∞). Therefore, hypothesis (H1) is

∂χ(s, 0) sp ≤ sp = p , 1 + ηp ∂p   ∂χ(s, 0)  = 1 > 0 for all s > 0. ∂p

χ(s, p) =

It follows that (H2) is satisfied. Clearly, hypothesis (H3) holds true. Moreover,   –ηs ∂ χ(s, p) = < 0. ∂p g4 (p) (1 + ηp)2 Therefore, hypothesis (H4) holds true and Theorems 1–3 are applicable. The parameters R0 and R1 for this application are given by

R0 =

(1 – ε1 )(1 – ε2 ){N(a1 λ1 e–μ1 τ1 + (a1 + β2 )λ2 e–μ2 τ2 ) + Mλ3 e–μ3 τ3 (a1 + β2 )} s0 , β5 (a1 + β2 )

R1 =

(1 – ε1 )(1 – ε2 ){N(a1 λ1 e–μ1 τ1 + (a1 + β2 )λ2 e–μ2 τ2 ) + Mλ3 e–μ3 τ3 (a1 + β2 )} s2 . β5 (a1 + β2 ) 1 + ηp2

Remark There are several forms of the general function χ(s, p) where (H1)–(H4) can be satisfied such as: sp , (i) Holling-type incidence χ(s, p) = 1+η 1s (ii) Beddington–DeAngelis incidence χ(s, p) = 1+η1sps+η2 p , sp (iii) Crowley–Martin incidence χ(s, p) = (1+η1 s)(1+η , 2 p) (iv) Hill-type incidence χ(s, p) =

sm p . ηm +sm

Now we are ready to perform some numerical simulations for system (66)–(71). The data of system (66)–(71) are provided in Table 1. We let τ1 = τ2 = τ3 = τ . •Effect of the parameters λi and r on the stability of the steady states To discuss our global results, we choose three different initial conditions: IC1: (s(0), w(0), y(0), u(0), p(0), x(0)) = (900, 5, 5, 15, 3, 3). Table 1 The data of example (66)–(71) Parameter

Value

Parameter

Value

Parameter

Value

Parameter

Value

ρ β1 ω

10 0.01 0.0001 1200 0.5

β2 β3 β4 β5 β6

0.02 0.36 0.031 3.0 0.1

μ1 μ2 μ3 η

1 1 1 0.01 0.2

N M

6 3 Varied Varied Varied

smax q

a1

ε1 , ε2 , τ λi , i = 1, 2, 3 r

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 26 of 36

Figure 1 The concentration of uninfected CD4+ T cells for system (66)–(71)

Figure 2 The concentration of latently infected cells for system (66)–(71)

IC2: (s(0), w(0), y(0), u(0), p(0), x(0)) = (700, 7, 8, 30, 5, 5). IC3: (s(0), w(0), y(0), u(0), p(0), x(0)) = (500, 15, 18, 60, 12, 7). Let us address three cases for the parameters λi , i = 1, 2, 3, and r. We assume that ε1 = ε2 = 0 (there is no treatment) and τi = 0 (there is no time delay). Case (I): Choose λi = 0.0000625 and r = 0.005, which gives R0 = 0.3016 < 1 and R1 = 0.1917 < 1. Therefore, based on Lemma 2 and Theorem 1, the system has unique steady state, that is, 0 , and it is GAS. As we can see from Figs. 1–6, the concentration of the uninfected CD4+ T cells is increased and approaches its normal value before infection,

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 27 of 36

Figure 3 The concentration of short-lived productively infected cells for system (66)–(71)

Figure 4 The concentration of long-lived productively infected cells for system (66)–(71)

that is, s0 = 1001.98; while concentrations of the other compartments converge to zero for all the three initial conditions. As a result, the HIV-1 is removed from the plasma. Case (II): By taking λi = 0.000625 and r = 0.005. For these values, R1 = 0.6095 < 1 < R0 = 3.0163. Consequently, based on Lemma 2 and Theorem 2, the humoral-inactivated infection steady state 1 is positive and is GAS. Figures 1–6 confirm that the numerical results support the theoretical results presented in Theorem 2. It can be observed that the variables of the model eventually converge to 1 = (366.317, 9.72758, 11.3488, 69.0344,

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 28 of 36

Figure 5 The concentration of free virus particles for system (66)–(71)

Figure 6 The concentration of B cells for system (66)–(71)

10.3112, 0.0) for all the three initial conditions. This case corresponds to a chronic HIV-1 infection in the absence of immune response. Case (III): λi = 0.000625 and r = 0.05. Then we calculate R0 = 3.0163 > 1 and R1 = 2.1648 > 1. According to Lemma 2 and Theorem 3, the humoral-activated infection steady state 2 is positive and is GAS. We can see from Figs. 1–6 that there is a consistency between the numerical results and theoretical results of Theorem 3. The states of the system converge to 2 = (734.586, 4.09194, 4.77393, 29.0396, 2.0, 7.01238) for all the three initial

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 29 of 36

Figure 7 The concentration of uninfected CD4+ T cells for system (66)–(71)

Figure 8 The concentration of latently infected cells for system (66)–(71)

conditions. In this case the humoral immune response is activated and can control the disease. •Effect of the HAART on the HIV dynamics We take τi = 0, λi = 0.000625, and r = 0.05. We choose the initial conditions (s(0), w(0), y(0), u(0), p(0), x(0)) = (850, 5, 6, 17, 1.5, 5). In Figs. 7–12 we show the effect of the drug efficacy parameters ε1 and ε2 on the HIV dynamics. Also, we can observe that, as the drug efficacy parameters ε1 and ε2 are increased, the concentration of uninfected cells is increased, while the concentrations of free virus particles and the three types of infected cells are decreased. Table 2 shows that the values of R0 and R1 are decreased as ε1 and ε2 are increased. Let us define the overall HAART effect as εe = ε1 + ε2 – ε1 ε2 [13]. If εe = 0, then the HAART has no effect, if εe = 1, the HIV growth is completely halted. Consequently, the

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 30 of 36

Figure 9 The concentration of short-lived productively infected cells for system (66)–(71)

Figure 10 The concentration of long-lived productively infected cells for system (66)–(71)

parameter R0 is given by R0 (εe ) =

(1 – εe ){N(a1 λ1 e–μ1 τ1 + (a1 + β2 )λ2 e–μ2 τ2 ) + Mλ3 e–μ3 τ3 (a1 + β2 )} s0 . β5 (a1 + β2 )

Since the goal is to clear the viruses from the body, we have to determine the drug efficacy that makes R0 (εe ) ≤ 1 for system (66)–(71). In this case, we get the critical drug efficacy (i.e., the efficacy needed in order to stabilize the system around the uninfected steady state). For the model (66)–(71), 0 is GAS when R0 (εe ) ≤ 1, i.e., εecrit ≤ εe < 1,

  R0 (0) – 1 . εecrit = max 0, R0 (0)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 31 of 36

Figure 11 The concentration of free virus particles for system (66)–(71)

Figure 12 The concentration of B cells for system (66)–(71)

Table 2 The values of steady states R0 and R1 for model (66)–(71) with different values of ε1 and ε2 Drug

Steady states

R0

R1

ε1 = 0, ε2 = 0 ε1 = 0.1, ε2 = 0.2 ε1 = 0.3, ε2 = 0.4 ε1 = 0.5, ε2 = 0.6

2 = (734.586, 4.09194, 4.77393, 29.0396, 2.0, 7.01238) 2 = (734.586, 3.68275, 4.29654, 26.1356, 2.0, 3.36892) 1 = (801.457, 2.14802, 2.50603, 15.244, 1.36614, 0) 0 = (1001.98, 0, 0, 0, 0, 0)

3.01635 2.17177 1.26687 0.60327

2.16484 1.55868 0.909233 0.432968

Using the data in Table 1, we have εecrit = 0.668473. •Effect of the time delay on the stability of the system Choosing ε1 = ε2 = 0, λi = 0.000625, and r = 0.05. The initial conditions are considered to be (s(0), w(0), y(0), u(0), p(0), x(0)) = (850, 2, 3, 17, 1, 5). Figures 13–18 and Table 3 show

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 32 of 36

Figure 13 The concentration of uninfected CD4+ T cells for system (66)–(71)

Figure 14 The concentration of latently infected cells for system (66)–(71)

the effect of the time delay parameter τ on the stability of 0 , 1 , and 2 . Clearly, the parameter τ has similar effect as the drug efficacy parameters ε1 and ε2 .

5 Conclusion In this paper, we have proposed and analyzed two general nonlinear HIV infection models with humoral immune response. We have considered three types of infected cells: latently infected cells, long-lived productively infected cells, and short-lived productively infected cells. We have incorporated three discrete or distributed time delays into the models. We have considered more general nonlinear functions for the HIV-target incidence rate, pro-

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 33 of 36

Figure 15 The concentration of short-lived productively infected cells for system (66)–(71)

Figure 16 The concentration of long-lived productively infected cells for system (66)–(71)

duction/proliferation, and removal rates of the cells and HIV. We have derived a set of conditions on these general functions and have determined two threshold parameters: the basic reproduction number R0 and the humoral immune response activation number R1 . We have proved the nonnegativity and ultimate boundedness of the model’s solutions and the existence and stability of the model’s steady states. Using Lyapunov functionals, we have established the global stability of the three steady states of the models. We have presented an example and performed some numerical simulations to support our theoretical results.

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 34 of 36

Figure 17 The concentration of free virus particles for system (66)–(71)

Figure 18 The concentration of B cells for system (66)–(71)

Table 3 The values of steady states R0 and R1 for model (66)–(71) with different values of τ Drug

Steady states

R0

R1

τ τ τ τ τ

2 = (734.586, 4.09194, 4.77393, 29.0396, 2.0, 7.01238) 2 = (734.586, 2.48189, 2.89554, 17.6134, 2.0, 1.89241) 1 = (826.24, 1.09341, 1.27564, 7.75968, 1.15901, 0) 0 = (1001.98, 0, 0, 0, 0, 0) 0 = (1001.98, 0, 0, 0, 0, 0)

3.01635 1.82951 1.22636 0.673038 0.247597

2.16484 1.31304 0.880158 0.483041 0.177701

=0 = 0.5 = 0.9 = 1.5 = 2.5

Acknowledgements This article was funded by the Deanship of Scientific Research (DSR) at King Abdulaziz University, Jeddah. The authors, therefore, acknowledge with thanks DSR for technical and financial support.

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 35 of 36

Competing interests The authors declare that they have no competing interests. Authors’ contributions All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript. Author details 1 Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia. 2 Department of Mathematics, Faculty of Science, Sohag University, Sohag, Egypt. 3 Department of Mathematics, Faculty of Science, King Khalid University, Abha, Saudi Arabia.

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Received: 13 December 2017 Accepted: 13 February 2018 References 1. Wodarz, D., Nowak, M.A.: Mathematical models of HIV pathogenesis and treatment. BioEssays 24, 1178–1187 (2002) 2. Rong, L., Perelson, A.S.: Modeling HIV persistence, the latent reservoir, and viral blips. J. Theor. Biol. 260, 308–331 (2009) 3. Nowak, M.A., Bangham, C.R.M.: Population dynamics of immune responses to persistent viruses. Science 272, 74–79 (1996) 4. Nowak, M.A., May, R.M.: Virus Dynamics: Mathematical Principles of Immunology and Virology. Oxford University Press, Oxford (2000) 5. Hattaf, K., Yousfi, N.: A generalized virus dynamics model with cell-to-cell transmission and cure rate. Adv. Differ. Equ. 2016, 174 (2016) 6. Wang, J., Teng, Z., Miao, H.: Global dynamics for discrete-time analog of viral infection model with nonlinear incidence and CTL immune response. Adv. Differ. Equ. 2016, 143 (2016) 7. Kang, C., Miao, H., Chen, X., Xu, J., Huang, D.: Global stability of a diffusive and delayed virus dynamics model with Crowley–Martin incidence function and CTL immune response. Adv. Differ. Equ. 2017, 324 (2017) 8. Elaiw, A.M., Raezaha, A.A., Shehata, A.M.: Stability of general virus dynamics models with both cellular and viral infections. J. Nonlinear Sci. Appl. 10, 1538–1560 (2017) 9. Li, X., Fu, S.: Global stability of a virus dynamics model with intracellular delay and CTL immune response. Math. Methods Appl. Sci. 38, 420–430 (2015) 10. Perelson, A.S., Nelson, P.W.: Mathematical analysis of HIV-1 dynamics in vivo. SIAM Rev. 41, 3–44 (1999) 11. Roy, P.K., Chatterjee, A.N., Greenhalgh, D., Khan, Q.J.A.: Long term dynamics in a mathematical model of HIV-1 infection with delay in different variants of the basic drug therapy model. Nonlinear Anal., Real World Appl. 14, 1621–1633 (2013) 12. Yang, Y., Xu, Y.: Global stability of a diffusive and delayed virus dynamics model with Beddington–DeAngelis incidence function and CTL immune response. Comput. Math. Appl. 71(4), 922–930 (2016) 13. Huang, Y., Rosenkranz, S.L., Wu, H.: Modeling HIV dynamic and antiviral response with consideration of time-varying drug exposures, adherence and phenotypic sensitivity. Math. Biosci. 184(2), 165–186 (2003) 14. Lv, C., Huang, L., Yuan, Z.: Global stability for an HIV-1 infection model with Beddington–DeAngelis incidence rate and CTL immune response. Commun. Nonlinear Sci. Numer. Simul. 19, 121–127 (2014) 15. Shu, H., Wang, L., Watmough, J.: Global stability of a nonlinear viral infection model with infinitely distributed intracellular delays and CTL immune responses. SIAM J. Appl. Math. 73(3), 1280–1302 (2013) 16. Hattaf, K., Yousfi, N., Tridane, A.: Stability analysis of a virus dynamics model with general incidence rate and two delays. Appl. Math. Comput. 221, 514–521 (2013) 17. Hattaf, K., Yousfi, N., Tridane, A.: Mathematical analysis of a virus dynamics model with general incidence rate and cure rate. Nonlinear Anal., Real World Appl. 13(4), 1866–1872 (2012) 18. Tian, X., Xu, R.: Global stability and Hopf bifurcation of an HIV-1 infection model with saturation incidence and delayed CTL immune response. Appl. Math. Comput. 237, 146–154 (2014) 19. Huang, G., Takeuchi, Y., Ma, W.: Lyapunov functionals for delay differential equations model of viral infections. SIAM J. Appl. Math. 70(7), 2693–2708 (2010) 20. Elaiw, A.M., Azoz, S.A.: Global properties of a class of HIV infection models with Beddington–DeAngelis functional response. Math. Methods Appl. Sci. 36, 383–394 (2013) 21. Elaiw, A.M.: Global properties of a class of HIV models. Nonlinear Anal., Real World Appl. 11, 2253–2263 (2010) 22. Huang, D., Zhang, X., Guo, Y., Wang, H.: Analysis of an HIV infection model with treatments and delayed immune response. Appl. Math. Model. 40(4), 3081–3089 (2016) 23. Monica, C., Pitchaimani, M.: Analysis of stability and Hopf bifurcation for HIV-1 dynamics with PI and three intracellular delays. Nonlinear Anal., Real World Appl. 27, 55–69 (2016) 24. Li, M.Y., Wang, L.: Backward bifurcation in a mathematical model for HIV infection in vivo with anti-retroviral treatment. Nonlinear Anal., Real World Appl. 17, 147–160 (2014) 25. Liu, S., Wang, L.: Global stability of an HIV-1 model with distributed intracellular delays and a combination therapy. Math. Biosci. Eng. 7(3), 675–685 (2010) 26. Elaiw, A.M., AlShamrani, N.H.: Stability of a general delay-distributed virus dynamics model with multi-staged infected progression and immune response. Math. Methods Appl. Sci. 40(3), 699–719 (2017) 27. Elaiw, A.M., Althiabi, A.M., Alghamdi, M.A., Bellomo, N.: Dynamical behavior of a general HIV-1 infection model with HAART and cellular reservoirs. J. Comput. Anal. Appl. 24(4), 728–743 (2018) 28. Elaiw, A.M., Almuallem, N.A.: Global dynamics of delay-distributed HIV infection models with differential drug efficacy in cocirculating target cells. Math. Methods Appl. Sci. 39, 4–31 (2016)

Elaiw et al. Advances in Difference Equations (2018) 2018:85

Page 36 of 36

29. Alshorman, A., Wang, X., Meyer, J., Rong, L.: Analysis of HIV models with two time delays. J. Biol. Dyn. 11(S1), 40–64 (2017) 30. Wang, X., Tang, S., Song, X., Rong, L.: Mathematical analysis of an HIV latent infection model including both virus-to-cell infection and cell-to-cell transmission. J. Biol. Dyn. 11(S2), 455–483 (2017) 31. Wang, X., Elaiw, A.M., Song, X.: Global properties of a delayed HIV infection model with CTL immune response. Appl. Math. Comput. 218, 9405–9414 (2012) 32. Callaway, D.S., Perelson, A.S.: HIV-1 infection and low steady state viral loads. Bull. Math. Biol. 64, 29–64 (2002) 33. Buonomo, B., Vargas-De-Le, C.: Global stability for an HIV-1 infection model including an eclipse stage of infected cells. J. Math. Anal. Appl. 385, 709–720 (2012) 34. Wang, H., Xu, R., Wang, Z., Chen, H.: Global dynamics of a class of HIV-1 infection models with latently infected cells. Nonlinear Anal., Model. Control 20(1), 21–37 (2012) 35. Korobeinikov, A.: Global properties of basic virus dynamics models. Bull. Math. Biol. 66, 879–883 (2004) 36. Pankavich, S.: The effects of latent infection on the dynamics of HIV. Differ. Equ. Dyn. Syst. (2015). https://doi.org/10.1007/s12591-014-0234-6 37. Elaiw, A.M., AlShamrani, N.H.: Global stability of humoral immunity virus dynamics models with nonlinear infection rate and removal. Nonlinear Anal., Real World Appl. 26, 161–190 (2015) 38. Hlavacek, W.S., Stilianakis, N.I., Perelson, A.S.: Influence of follicular dendritic cells on HIV dynamics. Philos. Trans. R. Soc. Lond. B, Biol. Sci. 355, 1051–1058 (2000) 39. Hale, J.K., Lunel, S.M.V.: Introduction to Functional Differential Equations. Springer, New York (1993) 40. Yang, X., Chen, L.S., Chen, J.F.: Permanence and positive periodic solution for the single-species nonautonomous delay diffusive models. Comput. Math. Appl. 32(4), 109–116 (1996)