An age-structured within-host HIV model with T-cell competition

Nonlinear Analysis: Real World Applications 38 (2017) 1–20

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications www.elsevier.com/locate/nonrwa

An age-structured within-host HIV model with T-cell competition Yan Wang a , Kaihui Liu b , Yijun Lou b,∗ a

College of Science, China University of Petroleum (East China), Qingdao, Shandong 266580, China Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong b

highlights

• • • • •

The effects of viral shedding and T-cell competition are studied via a hybrid model. The well-posedness of the model is rigorously verified. Local stability of the infection-free equilibrium in terms of R0 is investigated. Two special forms of viral production rates are investigated. Reduced DDE models based on two special forms exhibit backward bifurcation.

article

info

Article history: Received 1 October 2016 Accepted 5 April 2017

Keywords: HIV Competition Age-structured model Stability Backward bifurcation

abstract Recent studies demonstrate that resource competition is an essential component of T-cell proliferation in HIV progression, which can contribute instructively to the disease development. In this paper, we formulate an age-structured withinhost HIV model, in the form of a hyperbolic partial differential equation (PDE) for infected target cells coupled with two ordinary differential equations for uninfected T-cells and the virions, to explore the effects of both the T-cell competition and viral shedding variations on the viral dynamics. The basic reproduction number is derived for a general viral production rate which determines the local stability of the infection-free equilibrium. Two special forms of viral production rates, which are extensively investigated in previous literature, the delayed exponential distribution and a step function rate, are further investigated, where the original system can be reduced into systems of delay differential equations. It is confirmed that there exists a unique positive equilibrium for two special viral production rates when the basic reproduction number is greater than one. However, the model exhibits the phenomenon of backward bifurcation, where two positive steady states coexist with the infection-free equilibrium when the basic reproduction number is less than one. © 2017 Elsevier Ltd. All rights reserved.

∗ Corresponding author. E-mail address: [email protected] (Y. Lou).

http://dx.doi.org/10.1016/j.nonrwa.2017.04.002 1468-1218/© 2017 Elsevier Ltd. All rights reserved.

2

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

1. Introduction Human immunodeficiency virus (HIV) is a retrovirus that eventually induces the infected individuals into acquired immune deficiency syndrome (AIDS). It mainly infects lymphocytes with CD4 molecules, especially CD4+ T-cells, the main driver of the immune response. The within-host virus dynamics have been extensively studied, which greatly improves our understanding on disease progression and control. The pioneer within-host HIV mathematical models are based on a linear growth of healthy CD4+ T-cells [1–3]. However, recent medical reports indicate that the proliferation rate of T-cells is density-dependent on the number of T-cells [4,5], and therefore, several mathematical models have been proposed with a logistic growth of T-cells [6–11]. More recently, Borghans et al. demonstrated that incorporating T-cell competition fits the data significantly better than proliferation functions without competition, and the T-cell competition plays an important role in vitro [12]. In another study, Gillet et al. illustrated that there exists competition among infected T-cells within each host through experimental data [13]. Therefore, it becomes necessary to include T-cell competition within the host environment. Competition among individuals in population dynamics has attracted a growing body of studies, see for example [14–16]. Several with-host models have also been proposed with T-cell competition [17–19]. For example, Wang and Ellermeyer [17] constructed an HIV model with all T-cell competition (both healthy T-cells and infected T-cells) and concluded that the stability of the unique positive equilibrium depends on the magnitude of the competition intensity. A further study [18] investigated the effect of HIV intracellular delay on the stability of the positive equilibrium and examined the occurrence of Hopf bifurcation by using the intracellular delay as a bifurcation parameter. After that, Li and Wang [19] considered two different competition rates based on the literature [17], and discovered that the backward bifurcation can occur in a suitable parameter region. A general extension of the basic model with all T-cell competition can be described as follows:   T + qI dT = λ − dT + r1 T 1 − − kV T, dt Tmax T + qI dI (1) = kV T − δI − r2 I , dt Tmax dV = N δI − cV. dt In this model, T and I represent the numbers of healthy T-cells and infected T-cells respectively, and V denotes the HIV viral load. Healthy T-cells are produced from the precursors at a rate λ + r1 T . The natural death rates of the healthy and infected T-cells are dT and δI respectively. The infection rate of healthy T-cells infected by virions is kV T . The number of new viral particles generated by one infected T-cell (burst size) is a constant N , and the clearance rate of virions is cV . To describe the competition of T-cells, density-dependent terms r1 T (T + qI)/Tmax and r2 I(T + qI)/Tmax are used where Tmax , r1 , r2 and q are constants [17–19]. In these terms, the constant q can take values in the interval [0, 1]: q = 0 represents that the competition terms are only dependent on the healthy T-cells, while q ̸= 0 implies that the competition terms are also dependent on both the healthy and infected T-cells. In fact, this model has been investigated in several previous studies: when r1 = 0 and r2 = 0, it is the basic HIV model [1,2,7]; when r2 = 0 and q = 0, it becomes the model with a simplified logistic growth term [7]; when r2 = 0 and q = 1, it is the model with a full logistic growth term [6,8]; when r1 = r2 = r and q = 1, this model is the competition model in [17]; when r1 ̸= r2 and q = 1, it is the model in [19]. In this paper, we extend the T-cell competition model (1) by including the age structure of infected T-cells with the consideration of the infection age-dependent variations on the viral shedding rate. Age-structured HIV models have been developed in previous studies [20–25]. However, to the best of our knowledge, the current work is the first one to explore the joint effects of infection age-dependent variations and competition

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

3

among T-cells on virus dynamics. To do that, we first formulate an age-structured HIV model with T-cell competition consisting of a hyperbolic partial differential equation (PDE) for infected target cells coupled with two ordinary differential equations for uninfected T-cells and virions in the next section, for which the well-posedness of the model system is validated rigorously. In Section 3, we derive the basic reproduction number and analyze the local asymptotic stability of the infection-free equilibrium. However, further analysis seems challenging for the general viral shedding distribution. Therefore, we focus on two special distributions extensively used in previous literature, in the form of a delayed exponential distribution and a step function rate. In these two cases, the original system can be reduced into a delay differential system with distributed delay. The existence of the unique positive equilibrium is confirmed in Section 4 for these two typical viral production rates if the basic reproduction number is greater than one. However, backward bifurcation is possible where two positive equilibria and one infection-free equilibrium coexist when the basic reproduction number is smaller than one. Finally, we summarize our work in Section 5. 2. The model and its well-posedness Incorporating the infected cells structured by infection age, system (1) can be extended into the following form:   ∞ T (t) + q 0 i(a, t)da dT (t) = λ − dT (t) + r1 T (t) 1 − − kV (t)T (t), dt Tmax  ∞ dV (t) (2) = p(a)i(a, t)da − cV (t), dt 0 ∞ T (t) + q 0 i(a, t)da ∂i(a, t) ∂i(a, t) + = −δi(a, t) − r2 i(a, t) ∂t ∂a Tmax with the boundary and initial conditions for the density of infected cells i(0, t) = kV (t)T (t),

i(a, 0) = i0 (a).

To make sure the solution is continuous, we assume that i(0, t) = i(a, 0) holds when a = t = 0, that is, kV (0)T (0) = i0 (0). Instead of using I(t) for the number of infected CD4+ T-cells in system (1), here i(a, t) is used to denote the density of infected CD4+ T-cells with infection age a at time t. Its dynamics can be described in the form of a McKendrick von-Foerster equation [14] for an age-structured population growth. Then the total number of infected T-cells at time t is  ∞ I(t) = i(a, t)da. 0

Most parameters are same as model (1). The kernel function p(a) is the viral shedding rate of an infected T-cell with infection age a. Two typical forms in the literatures [26–29] are   p∗ (1 − e−θ(a−a1 ) ), a ≥ a1 , p∗ , a ≥ a1 , p(a) = and p(a) = 0, a < a1 , 0, a < a1 . Here, a1 is the intracellular delay for infected cells shedding the virus. We defer the investigation of these two special cases in Section 4. For the biological model (2), it is essential to validate the well-posedness of the system, that is the existence and uniqueness of a non-negative solution for a given non-negative initial value. This is provided in Theorems 2.1 and 2.2. Here, we follow the approach similar to literature [30] by reformulating system (2) as a perturbation of a linear abstract Cauchy problem and studying the existence and uniqueness of the solution of this abstract Cauchy problem in the state space X := R × R × L1 ([0, ∞); R) based on the theory of nonlinear semigroups.

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

4

Let ˜i(t) ∈ L1 ([0, ∞); R) (with the usual L1 norm) be defined as: ˜i(t) : a → i(a, t), [0,∞)

a ≥ 0.

R

Then, the initial condition for i(a, t) when t = 0 is given by ˜i(0)(a) = i0 (a). Let X := R × R × L1 ([0, ∞); R) be the state space of model (2) with the norm ∥ · ∥X defined by:  ∞ ∥x∥X := |x1 | + |x2 | + |x3 (a)|da, for all x = (x1 , x2 , x3 ) ∈ X. 0

Then, for all t ≥ 0, we can write x(t) = (T (t), V (t), ˜i(t)) ∈ X. System (2) can be formulated as a perturbation of a linear abstract Cauchy problem: d x(t) + Ax(t) = F (x(t)), dt x(0) = (T (0), V (0), i0 (·)),

t > 0,

(3)

where the linear operator A and the nonlinear operator F on X are defined by: A1 (x1 , x2 , x3 ) = −(r1 − d)x1 , ∞

A2 (x1 , x2 , x3 ) = cx2 −

p(a)x3 (a)da, 0

∂ (x3 (a)) + δx3 (a), A3 (x1 , x2 , x3 )(a) = ∂a and F1 (x1 , x2 , x3 ) = λ −

c1 x21

 − c2 x1

∞

 x3 (a)da − kx1 x2 ,

0

F2 (x1 , x2 , x3 ) = 0, 

∞

F3 (x1 , x2 , x3 )(a) = −c3 x1 x3 (a) − c4 x3 (a)

 x3 (a)da ,

0 r1

r1 q Tmax , c3

r2

r2 q Tmax .

= Tmax and c4 = with c1 = Tmax , c2 = The following theorem shows the existence and uniqueness of a solution for the abstract Cauchy problem (3). Theorem 2.1. System (3) has a unique solution x(t) on [0, ∞) with initial value x0 = (T (0), V (0), i0 (·)) ∈ X. Proof. To obtain the existence of a unique solution on [0, ∞), it follows from [31, Chapter 6, Theorem 1.5] that we have to verify the following two conditions: (I) −A is the infinitesimal generator of a C0 -semigroup on X; (II) The nonlinear operator F is continuously differentiable in X. In what follows, we show these two conditions hold. Since −A corresponds to the linear homogeneous part of system (3) and the domain of A is D(A) = {x : x ∈ R × R × C 1 ([0, ∞)) such that A3 (x) ∈ L1 ([0, ∞); R)}, it is easy to check that −A is the infinitesimal generator of a C0 -semigroup on X. Therefore, (I) holds. Now we verify condition (II) through two steps. Step 1: For all x ∈ X, w ∈ X, we need to show the existence of the operator DF (x; w), which is defined by DF (x; w) := lim

h→0

F (x + hw) − F (x) . h

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

5

It is obvious that DF2 (x; w) = 0, and the differentiations for F1 and F3 are computed as follows: F1 (x + hw) − F1 (x) h   = lim −c1 (x1 + hw1 )2 − c2 (x1 + hw1 )

DF1 (x; w) = lim

h→0

∞

 (x3 (a) + hw3 (a))da h→0 0    ∞ 2 x3 (a)da + kx1 x2 h − k(x1 + hw1 )(x2 + hw2 ) + c1 x1 + c2 x1 0    ∞  ∞ w3 (a)da − kw2 x1 − kw1 x2 , x3 (a)da − c2 x1 = −2c1 x1 w1 − c2 w1 0

0

and F3 (x + hw)(a) − F3 (x)(a) h  = lim −c3 (x1 + hw1 )(x3 (a) + hw3 (a)) − c4 (x3 (a) + hw3 (a)) h→0  ∞   ∞  (x3 (a) + hw3 (a))da + c3 x1 x3 (a) + c4 x3 (a) x3 (a)da h 0  ∞  0  ∞  = −c3 x1 w3 (a) − c3 w1 x3 (a) − c4 w3 (a) x3 (a)da − c4 x3 (a) w3 (a)da .

DF3 (x; w)(a) = lim

h→0

0

0

Thus, the existence of these limits implies that the operator DF (x; ·) exists. Step 2: We have to show that the operator DF (x; ·) is continuous in x, that is, lim ∥DF (x; ·) − DF (y; ·)∥OP = 0,

for ∥x − y∥X → 0,

where ∥ · ∥OP is the operator norm defined as sup ∥DF (x; w) − DF (y; w)∥X ∥w∥X ≤1

 =

sup

 |DF1 (x; w) − DF1 (y; w)| + |DF2 (x; w) − DF2 (y; w)| +

∥w∥X ≤1

∞

 |(DF3 (x; w) − DF3 (y; w))(a)|da .

0

Since sup ∥DF1 (x; w) − DF1 (y; w)∥X ∥w∥X ≤1

 ∞ ≤ 2c1 |x1 − y1 ||w1 | + c2 |w1 | |x3 (a) − y3 (a)|da 0  ∞  + c2 |x1 − y1 | |w3 (a)|da + k|w2 ||x1 − y1 | + k|w1 ||x2 − y2 | 0  ∞ ≤ 2c1 |x1 − y1 | + c2 |x3 (a) − y3 (a)|da + c2 |x1 − y1 | + k|x1 − y1 | + k|x2 − y2 | 0  ∞ |x3 (a) − y3 (a)|da = (2c1 + c2 + k)|x1 − y1 | + k|x2 − y2 | + c2 0

≤ (2c1 + c2 + k)∥x − y∥. It follows that sup ∥DF1 (x; ·) − DF1 (y; ·)∥X → 0, ∥w∥X ≤1

when ∥x − y∥X → 0.

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

6

The operator norm for DF3 can be estimated through the following process:  ∞ |(DF3 (x; w) − DF3 (y; w))(a)|da sup ∥w∥X ≤1 0  ∞

 ∞ c3 |x1 − y1 ||w3 (a)|da + c3 |w1 ||x3 (a) − y3 (a)|da 0 0    ∞  ∞  ∞  ∞ |w3 (s)|ds da |x3 (s) − y3 (s)|ds da + c4 |x3 (a) − y3 (a)| + c4 |w3 (a)| 0 0 0 0  ∞  ∞ = c3 |x1 − y1 | |w3 (a)|da + c3 |w1 | |x3 (a) − y3 (a)|da 0 0  ∞  ∞  ∞  ∞ |x3 (a) − y3 (a)|da |w3 (s)|ds |w3 (a)|da + c4 |x3 (s) − y3 (s)|ds + c4 0 0 0 0  ∞  ∞  ∞ ≤ c3 |x1 − y1 | + c3 |x3 (a) − y3 (a)|da + c4 |x3 (s) − y3 (s)|ds + c4 |x3 (a) − y3 (a)|da 0 0 0  ∞ ≤ c3 |x1 − y1 | + (c3 + 2c4 ) |x3 (a) − y3 (a)|da ≤

0

≤ (c3 + 2c4 )∥x − y∥. It is obvious that 

∞

|(DF3 (x; w) − DF3 (y; w))(a)|da → 0,

sup ∥w∥X ≤1

when ∥x − y∥X → 0.

0

Therefore, we have sup ∥DF (x; ·) − DF (y; ·)∥X → 0,

when ∥x − y∥X → 0,

∥w∥X ≤1

which completes the proof.



In order to show the non-negativeness of the solution through an assigned initial value, it is easier to check the integral form of i(a, t) instead of the partial differential equation. Integrating the third equation of system (2) with respect to a from 0 to ∞, we obtain ∞  ∞  ∞  ∞  T (t) + q 0 i(a, t)da ∞ ∂i(a, t) ∂i(a, t) da + da = −δ i(a, t)da − r2 i(a, t)da, ∂t ∂a Tmax 0 0 0 0 that is, ∂ ∂t



∞

 i(a, t)da + i(∞, t) − i(0, t) = −δ

0

∞

 i(a, t)da − r2

0

0

∞

∞ T (t) + q 0 i(a, t)da i(a, t)da . Tmax

(4)

Since no infected T-cells can survival forever, mathematically, we assume i(∞, t) = 0. Using I(t) = ∞ i(a, t)da and substituting i(0, t) = kV (t)T (t) into (4) yield 0 dI(t) T (t) + qI(t) = kV (t)T (t) − δI(t) − r2 I(t) . dt Tmax Integrating along the characteristic lines of a − t =constant [32], the following partial differential equation  T (t) + qI(t)  ∂i(a, t) ∂i(a, t) + = −δi(a, t) − r2 i(a, t), ∂t ∂a Tmax  i(0, t) = kV (t)T (t), i(a, 0) = i0 (a), can be solved with the explicit solution  i (a − t) K0 (a, T (t), I(t)) , a > t, 0 i(a, t) = K0 (a − t, T (t), I(t))  kV (t − a)T (t − a)K0 (a, T (t), I(t)), a ≤ t,

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

7

where   K0 (a, T (t), I(t)) = exp − 0

a



T (t − a + s) + qI(t − a + s) δ + r2 Tmax



 ds .

Then, we can rewrite system (2) into the following integral–differential system:   T (t) + qI(t) dT (t) = λ − dT (t) + r1 T (t) 1 − − kV (t)T (t), dt Tmax  ∞ dV (t) = p(a)i(a, t)da − cV (t), dt 0 dI(t) T (t) + qI(t) = kV (t)T (t) − δI(t) − r2 I(t) , dt Tmax  i (a − t) K0 (a, T (t), I(t)) , a > t, 0 i(a, t) = K0 (a − t, T (t), I(t))  kV (t − a)T (t − a)K0 (a, T (t), I(t)), a ≤ t.

(5)

(6)

Since any solution of (2) should satisfy system (6), we can verify the non-negativeness of solutions for the above system, which is stated in the next result. Before conducting the analysis, we make some biologically ˜, reasonable assumptions on the function p, that is, p ∈ L∞ + and there exists an essential upper bound p i.e. p˜ := ess supa∈[0,∞) p(a) < +∞, which means, the set Mp˜ = {a : p(a) > p˜} has a measure µ(Mp˜) = 0. Theorem 2.2. System (6) has a non-negative solution with non-negative initial values T (0) ≥ 0, V (0) ≥ ∞ 0, i0 (·) ≥ 0 and I(0) = 0 i0 (a)da ≥ 0. Furthermore, the solution is ultimately bounded.   Proof. From the first equation of system (6), we can get dTdt(t)  = λ > 0. By [33, Theorem 5.2.1], T (t)=0

T (t) ≥ 0 for all t ≥ 0. Suppose t0 is the first time instant such that V (t0 ) = 0. Since  ∞ dV (t) p(a)i(a, t)da − cV (t) = dt 0  t  ∞ = p(a)i(a, t)da + p(a)i(a, t)da − cV (t) 0 t  t = p(a)kV (t − a)T (t − a)K0 (a, T (t), I(t))da 0  ∞ K0 (a, T (t), I(t)) da − cV (t), + p(a)i0 (a − t) K 0 (a − t, T (t), I(t)) t we have  ∞ dV (t)  K0 (a, T (t0 ), I(t0 )) ≥ p(a)i0 (a − t0 ) da ≥ 0,  dt t=t0 K0 (a − t0 , T (t0 ), I(t0 )) t0 which implies the non-negativeness of V (t) for all t ≥ 0. Based on the non-negativeness of the variables V (t) and T (t) for all t ≥ 0, we have i(a, t) ≥ 0 and ∞ I(t) = 0 i(a, t)da ≥ 0 for all t ≥ 0. Next, we show the ultimate boundedness of solutions. From the first equation of (6), we have   T ˙ T ≤ λ − dT + r1 T 1 − , Tmax therefore,    Tmax 4r λ 1 lim sup T (t) ≤ T0 := r1 − d + (r1 − d)2 + . 2r1 Tmax t→+∞

8

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

Let G(t) = T (t) + I(t), the total number of target cells. Calculating the derivative of G(t) along the solutions of system (6), we obtain   dG(t) T (t) + qI(t) T (t) + qI(t) = λ − dT (t) + r1 T (t) 1 − − δI(t) − r2 I(t) dt Tmax Tmax ≤ −dT (t) − δI(t) + λ + r1 T (t) ≤ − min{d, δ}G(t) + λ + r1 T (t), and therefore, lim sup(T (t) + I(t)) ≤ t→+∞

λ + r1 T0 . min{d, δ}

It then follows from the above assumption for p that  ∞ dV (t) = p(a)i(a, t)da − cV (t) dt 0  ∞ ≤ p˜ i(a, t)da − cV (t) 0

= p˜I(t) − cV (t). Thus, V (t) is also ultimately bounded with lim sup V (t) ≤ p˜ t→+∞

This completes the proof.

λ + r1 T0 . min{d, δ}c



3. R0 and local stability of the infection-free equilibrium In this section, we are going to derive an important index, the basic reproduction number, which plays a pivotal role on characterizing the virus dynamics [1,23]. Furthermore, this index is shown to max predict the local stability of the infection-free equilibrium E0 (T0 , 0, 0), where T0 is given as T0 = T2r 1    1λ r1 − d + (r1 − d)2 + T4rmax . Let N be the total number of infectious virus particles (burst size) produced by one infected T-cell during its lifespan, then      ∞  ∞ r2 T0 N = p(a)σ(a)da = p(a) exp − δ + a da, (7) Tmax 0 0 with     r2 T0 σ(a) = exp − δ + a Tmax representing the probability for an infected T-cell surviving to age a. Following the idea in literature [23], we can define the basic reproduction number as     ∞ r2 T0 kT p(a) exp − δ + 0 0 Tmax a da kT0 N R0 = = . c c Biologically, R0 is the average number of newly infected T-cells produced by one infected T-cell during its lifespan when the human body is only with healthy T-cells. Suppose the maximum survival age of infected T-cells is Amax , then i(a, t) = 0 for all a > Amax and t ≥ 0. Since we are focusing on the long term dynamics, by time shift, we have i(a, t) = kV (t − a)

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

9

I(t − a)K0 (a, T (t), I(t)) for t > Amax ≥ a. Now, we investigate the local stability of the infection-free equilibrium E0 for the following system corresponding to system (6) for t > Amax :   T (t) + qI(t) dT (t) = λ − dT (t) + r1 T (t) 1 − − kV (t)T (t), dt Tmax  ∞ dV (t) (8) p(a)kV (t − a)T (t − a)K0 (a, T (t), I(t))da − cV (t), = dt 0 dI(t) T (t) + qI(t) = kV (t)T (t) − δI(t) − r2 I(t) , dt Tmax with K0 (a, T (t), I(t)) defined in (5). Then the corresponding linearized system at E0 (T0 , 0, 0) is   2r1 T0 r1 qT0 dT (t) = −d + r1 − T (t) − kT0 V (t) − I(t), dt T Tmax max  ∞ dV (t) kT0 V (t − a)p(a)σ(a)da − cV (t), = dt 0   dI(t) r2 T0 = kT0 V (t) − δ + I(t). dt Tmax

(9)

The next result indicates that the basic reproduction number determines the stability of the infection-free equilibrium E0 . Theorem 3.1. The infection-free equilibrium E0 of system (6) is locally asymptotically stable if R0 < 1, and it is unstable if R0 > 1. Proof. The characteristic equation of (9) at the infection-free equilibrium E0 is     2r T r1 qT0  −d + r1 − 1 0 − ξ −kT − 0   Tmax Tmax  ∞     −ξa 0 kT0 e p(a)σ(a)da − c − ξ 0  = 0.    0   r2 T0  − ξ  0 kT0 −δ −  Tmax By simplification, we have      ∞ 2r1 T0 r2 T0 −ξa ξ + d − r1 + ξ+δ+ ξ + c − kT0 e p(a)σ(a)da = 0. Tmax Tmax 0 Using the equality r1 − d =

r1 T0 Tmax

ξ1 = r1 − d −

−

λ T0 ,

we can easily obtain two eigenvalues

2r1 T0 λ r1 T0 =− − <0 Tmax T0 Tmax

and ξ2 = −δ −

r2 T0 < 0. Tmax

Then, the local stability of E0 is determined by the sign of the real part of  ∞ ξ = kT0 e−ξa p(a)σ(a)da − c. 0

Substituting the equality

kT0 c

= ∞ 0

R0 p(a)σ(a)da

into the above equation gives

ξ R0 +1 = ∞ c p(a)σ(a)da 0



∞

e−ξa p(a)σ(a)da.

0

When R0 < 1, suppose Reξ ≥ 0, then      ∞   ξ  R0    −ξa e p(a)σ(a)da =  + 1 ≥ 1.  ∞  0 p(a)σ(a)da 0  c

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

10

However,    ∞   R0   −ξa e p(a)σ(a)da  ∞   0 p(a)σ(a)da 0       ∞  R0    e−ξa p(a)σ(a)da =  ∞   0 p(a)σ(a)da  0  ∞ R0 p(a)σ(a)da = R0 < 1, ≤ ∞ p(a)σ(a)da 0 0 a contradiction. Therefore, Reξ < 0 and the infection-free equilibrium E0 is locally asymptotically stable if R0 < 1. In the case of R0 > 1, let ξ R0 ϕ(ξ) = + 1 −  ∞ c p(a)σ(a)da 0



∞

e−ξa p(a)σ(a)da.

0

Then, ϕ(0) = 1 − R0 < 0 and limξ→∞ ϕ(ξ) = ∞. Thus, there exists a ξ ∗ > 0 such that ϕ(ξ ∗ ) = 0, which further indicates that E0 is unstable.  Till now, we have established basic results on the well-posedness of the model, identified the basic reproduction number and illustrated its relation to the stability of a boundary equilibrium (the infectionfree equilibrium) for the general distribution function p(a). It becomes very challenging to perform further analysis, regrading to the existence and number of positive equilibria for this general kernel. In the next section, we focus on two special forms of p(a). 4. Special forms of p(a) In this section, we consider two typical viral production functions p(a) that are commonly used in previous literature [26–29], a delayed exponential function and a step function. 4.1. p(a) is a delayed exponential function Some previous studies [26,27,29] employed the following function to characterize the viral reproducing rate  p∗ (1 − e−θ(a−a1 ) ), a ≥ a1 , p(a) = 0, a < a1 . In this form, p∗ represents the maximum production rate, and θ indicates how quickly virions can be produced. Constant a1 represents an intracellular delay, that is, the time needed from the initial infection to releasing the first virus particle. Let B0 = δ + r2 T0 /Tmax , then the total number of virus particles (burst size) produced by all ages of an infected T-cell during its lifespan N in Eq. (7) can be computed as  ∞  ∞ N = p(a)σ(a)da = p∗ (1 − e−θ(a−a1 ) )e−B0 a da 0 a1   1 1 p∗ θ = e−B0 a1 . = p∗ e−B0 a1 − B0 θ + B0 B0 (θ + B0 ) In this case, the basic reproduction number becomes R0 =

kT0 N kT0 p∗ θ = · e−B0 a1 . c c B0 (θ + B0 )

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

11

Using the explicit expression for p(a), we also have  ∞  ∞  ∞ p(a)i(a, t)da = p∗ i(a, t)da − p∗ e−θ(a−a1 ) i(a, t)da 0 a1 a1  ∞  ∞ = p∗ i(m + a1 , t)dm − p∗ e−θm i(m + a1 , t)dm 0 0  ∞   ∞ ∗ −θm i(m, t − a1 )dm − e =p i(m, t − a1 )dm 0

0 a1

  × exp −

   T (t − a1 + s) + qI(t − a1 + s) δ + r2 ds , Tmax

0

(10)

where a1

  i(m + a1 , t) = i(m, t − a1 )exp − 0

   T (t − a1 + s) + qI(t − a1 + s) ds . δ + r2 Tmax

Now, we are in a position to reduce the original hybrid system with a partial differential equation into a delay differential equation, which is more mathematically tractable. For that purpose, we introduce a new variable  ∞ W (t) = e−θm i(m, t)dm. 0

Then,  W (t − a1 ) =

∞

e−θm i(m, t − a1 )dm.

0

Using the third equation of system (2) and the equality ∂(e−θm i(m, t − a1 )) ∂i(m, t − a1 ) = + θe−θm i(m, t − a1 ), ∂m ∂m

e−θm we obtain dW (t − a1 ) = dt

 0

∞

e−θm

∂i(m, t − a1 ) dm ∂t

 ∞ T (t − a1 ) + qI(t − a1 ) = −δ e−θm i(m, t − a1 )dm − r2 Tmax  0∞  ∞ ∂i(m, t − a1 ) × e−θm i(m, t − a1 )dm − e−θm dm ∂m 0 0  T (t − a1 ) + qI(t − a1 ) = − δ + r2 W (t − a1 ) Tmax   ∞ ∂(e−θm i(m, t − a1 )) − + θe−θm i(m, t − a1 ) dm ∂m 0  T (t − a1 ) + qI(t − a1 ) = − θ + δ + r2 W (t − a1 ) + i(0, t − a1 ) Tmax   T (t − a1 ) + qI(t − a1 ) = − θ + δ + r2 W (t − a1 ) + kV (t − a1 )T (t − a1 ). Tmax Therefore, we arrive at   dW (t) T (t) + qI(t) = − θ + δ + r2 W (t) + kV (t)T (t). dt Tmax

12

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

By identity (10) and the previous equation, system (6) can be reduced into the following differential system:   T (t) + qI(t) dT (t) = λ − dT (t) + r1 T (t) 1 − − kV (t)T (t), dt Tmax   a1  dV (t) T (t − a1 + s) + qI(t − a1 + s) ∗ = p exp − (δ + r2 )ds × [I(t − a1 ) − W (t − a1 )] − cV (t), dt Tmax 0 (11) dI(t) T (t) + qI(t) = kV (t)T (t) − δI(t) − r2 I(t) , dt Tmax   dW (t) T (t) + qI(t) = kV (t)T (t) − θ + δ + r2 W (t). dt Tmax Please note that this system involves both the discrete and continuous delay terms. Generally, the dynamical behavior becomes very complicated for all parameter values [14]. However, we are able to discuss the existence of positive equilibrium. Let E1 = (T1 , V1 , I1 , W1 ) be a positive equilibrium. Then, E1 can be calculated through solving the following equations together:   T1 + qI1 λ − dT1 + r1 T1 1 − − kV1 T1 = 0, Tmax   T1 + qI1 )a1 (I1 − W1 ) − cV1 = 0, p∗ exp −(δ + r2 Tmax   T1 + qI1 kV1 T1 − δ + r2 I1 = 0, Tmax   T1 + qI1 W1 = 0. kV1 T1 − θ + δ + r2 Tmax From the third and the last equations of the above system, we obtain   1 I1 δ + r2 TT1 +qI kV1 T1 max . V1 = , and W1 = kT1 θ + δ + r2 TT1 +qI1 max

Inserting these equalities back (I1 ̸= 0), we have two equations:     T1 + qI1 T1 + qI1 λ − dT1 + r1 T1 1 − − δ + r2 I1 = 0,   Tmax  Tmax    T1 + qI1 T1 + qI1 T1 + qI1 θ + δ + r2 exp δa1 + r2 a1 − kp∗ θT1 = 0. c δ + r2 Tmax Tmax Tmax Normally it is impossible to get the explicit relationship between T1 and I1 from the above equations. Therefore, we consider a simplified case when q = 0 in the subsequent subsection. 4.1.1. A simplified case: q = 0 In the case of system (11) with q = 0, we have   T1 λ − dT1 + r1 T1 1 − − kV1 T1 = 0, T   max r2 T1 p∗ exp − δ + a1 (I1 − W1 ) − cV1 = 0, T max   r2 T1 kV1 T1 − δ + I1 = 0, T max   r2 T1 kV1 T1 − θ + δ + W1 = 0. Tmax

(12)

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

13

T0 − Tλ0 , we have Adding the first and the last equations and noting r1 − d = Tr1max     2 T0 1) 1 λ − dT1 + r1 T1 1 − TTmax λ + Tr1max − Tλ0 T1 − r1T(T max W1 = = T1 T1 θ + δ + Tr2max θ + δ + Tr2max     r1 T1 r1 T1 r 1 T0 λ λ λ + − − + (T − T ) 0 1 T1 Tmax T0 Tmax T0 Tmax   = = > 0, r2 T1 r2 T1 θ + δ + θ + δ + Tmax /T1 Tmax

since T1 ∈ (0, T0 ). Similarly, I1 can be reformulated in terms of T1 through adding the first and the third equations of algebraic system (12):     T1 1 λ − dT1 + r1 T1 1 − TTmax (T0 − T1 ) Tλ0 + Tr1max I1 = = > 0. T1 T1 δ + Tr2max δ + Tr2max Then, substituting these expressions into the second equation of system (12) gives     T1 a1 (I1 − W1 ) p∗ exp − δ + Tr2max V1 = c       r 2 T1 T1 ∗ p θ exp − δ + Tmax a1 (T0 − T1 ) Tλ0 + Tr1max    = > 0. T1 T1 θ + δ + Tr2max c δ + Tr2max Since V1 , I1 and W1 can be represented in terms of T1 , it suffices to study the equation for T1 only:      r2 T1 r2 T1 r2 T1 c δ+ θ+δ+ exp δa1 + a1 − kp∗ θT1 = 0. Tmax Tmax Tmax For all x ≥ 0, let      r2 x r2 x r2 x θ+δ+ exp δa1 + a1 − kp∗ θx. f1 (x) = c δ + Tmax Tmax Tmax Then, f1′ (x)

       r2 x r2 x cr2 r2 x r2 x exp δa1 + a1 θ + 2 δ + + a1 δ + θ+δ+ − kp∗ θ, = Tmax Tmax Tmax Tmax Tmax

and f1′′ (x) = c



r2

2

Tmax

       r2 x r2 x r2 x r2 x exp δa1 + a1 2 + 2a1 θ + 2δ + 2 + a21 δ + θ+δ+ Tmax Tmax Tmax Tmax

> 0. Noting that R0 =

kT0 c

·

p∗ θ −B0 a1 , B0 (θ+B0 ) e

we have

f1 (T0 ) = cB0 eB0 a1 (θ + B0 ) − kp∗ θT0 = cB0 eB0 a1 (θ + B0 ) − cR0 B0 eB0 a1 (θ + B0 ) = cB0 eB0 a1 (1 − R0 )(θ + B0 ). It is obvious that f1 (0) = c(θ + δ)δeδa1 > 0. If R0 > 1, then f1 (T0 ) < 0. Noting that f1 (0) > 0 and f1′′ (x) > 0, there exists a unique positive solution ∗ T ∈ (0, T0 ) such that f1 (T ∗ ) = 0. Therefore, we have the following result. Theorem 4.1. If R0 > 1, system (11) with q = 0 has only one positive equilibrium E1 .

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

14

a

b

Fig. 1. The backward bifurcation of system (11) with q = 0 and a1 = 0.004. (a) Solutions either converge to a positive equilibrium or the infection-free equilibrium, dependent on initial values. (b) The backward bifurcation curve R0 versus I1 (the number of infected T-cells at the positive equilibrium). When R0 < Rc = 0.9388, there exists a unique stable infection-free equilibrium (solid blue line). If Rc < R0 < 1, there exist a stable disease-free equilibrium, a stable positive equilibrium, and an unstable positive equilibrium. When R0 > 1, there exists a unique positive equilibrium which is stable, and the infection-free equilibrium becomes unstable.

If R0 < 1, we have f1 (T0 ) > 0. Since f1 (0) > 0 and f1 (x) is a convex function (f1′′ (x) > 0), there are at most two positive solutions between 0 and T0 such that f1 (x) = 0. It follows from the above argument that two positive equilibria may coexist when R0 < 1. That is to say, system (11) with q = 0 may undergo the backward bifurcation. The following result shows the conditions for backward bifurcation. Theorem 4.2. If R0 < 1, system (11) with q = 0 exhibits backward bifurcation if there exists a T1 ∈ (0, T0 ) such that f1′ (T1 ) = 0 and f1 (T1 ) < 0. Theorem 4.2 presents only implicit sufficient conditions for the existence of backward bifurcation, which brings more challenges for further theoretical analysis. In the subsequent subsection, some numerical simulations will be performed to validate our theoretical results. 4.1.2. Numerical simulation: q = 0 and a1 ̸= 0 The numerical simulation is performed to show the feasibility of the existence of backward bifurcation of system (11) with a small delay (a1 = 0.004) when q = 0 as described in Theorem 4.2. The parameters are chosen based on references [9,11,23,29] as follows: λ = 1 mm−3 day−1 , r2 = 1.95 day

−1

δ = 0.95 day−1 ,

,

d = 0.02 day−1 , Tmax = 1500 mm

−3

r1 = 0.2 day−1 , ,

p∗ = 854 virions/cell,

a1 = 0.004 day,

k = 0.00065 mm day−1 , 3

c = 30 day−1 ,

θ = 0.32.

For this set of parameters, we have R0 = 0.9654 < 1 and T0 = 1355.5. It is easy to check that there exists a 875.6 < T1 < 875.7 < T0 such that f1′ (T1 ) = 0 and f1 (T1 ) < 0, which implies that the conditions in Theorem 4.2 are satisfied. We simulate solutions of system (11) with six different initial values when q = 0 and a1 = 0.004. It can be seen from Fig. 1(a) that solutions either converge to the infection-free equilibrium or a positive equilibrium. This indicates that there exist two stable equilibria, the infection-free equilibrium and another positive equilibrium. The backward bifurcation occurs when R0 ranges from 0.9388 to 1 by adjusting the parameter r2 from 1.8950 to 1.9944, the diagram of which is shown in Fig. 1(b), where three equilibria coexist when Rc = 0.9388 < R0 < 1. In particular, both the infection-free equilibrium and

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

15

the larger positive equilibrium are locally asymptotically stable, while the smaller positive equilibrium is unstable. It can be easily seen that f1 (x) is a transcendental function due to the existence of the delay term a1 , which brings a big challenge to obtain the explicit conditions. Thus, we study a more simpler case when there is no delay, that is a1 = 0 in the upcoming subsection. 4.1.3. A simpler case: q = 0 and a1 = 0 Applying similar calculations as in the previous subsection, we obtain a single equation to locate T1 for the positive equilibrium (T1 , V1 , I1 , W1 ):    r2 T1 r2 T1 θ+δ+ − kp∗ θT1 = 0. c δ+ Tmax Tmax Let    r2 x r2 x f2 (x) = c δ + θ+δ+ − kp∗ θx Tmax Tmax  2   r2 cr2 2 ∗ =c x + (θ + 2δ) − kp θ x + cδ(θ + δ), Tmax Tmax then we have f2 (0) = cδ(θ + δ) > 0 and f2 (T0 ) = cB0 (1 − R0 )(θ + B0 ), where B0 = δ + r2 T0 /Tmax as defined in the previous section. In the case of R0 < 1, to obtain the conditions for the occurrence of backward bifurcation, it is equivalent to investigate under what circumstances the parabolic function f2 (x) = 0 has two positive solutions between 0 and T0 . Based on the properties of parabolic curves, the following three compatible conditions in the next theorem can guarantee f2 (x) = 0 has two positive solutions in (0, T0 ). Theorem 4.3. If q = 0 and a1 = 0, system (11) exhibits backward bifurcation when R0 < 1 if and only if the following three conditions exist simultaneously: 2 (a) f2′ (0) = Tcr (θ + 2δ) − kp∗ θ < 0, max r2 ′ 2 )2 + Tcr (θ + 2δ) − kp∗ θ > 0, (b) f2 (T0 ) = 2cT0 ( Tmax  max 2 ∗ 2 2 (c) the discriminant ∆ = Tcr (θ + 2δ) − kp θ − 4δ(θ + δ)( Tcr )2 > 0. max max

4.1.4. Numerical simulation: q = 0 and a1 = 0 We perform numerical simulation to illustrate the feasibility of the existence of backward bifurcation of system (11) with q = 0 and a1 = 0 as shown in Theorem 4.3. The parameters set is merely the same as the one in the previous case except that λ = 7 mm−3 day−1 in this case. Thus, it follows that R0 = 0.97 < 1, f1′ (0) = −0.09105 < 0, f1′ (T0 ) = 0.04967 > 0 and ∆ = 0.00095 > 0, which make the conditions (a)–(c) in Theorem 4.3 hold. We feed the simulation with six different initial values, and find that solutions of system (11) with q = 0 and a1 = 0 either converge to the infection-free equilibrium or a positive equilibrium (see Fig. 2(a)). This implies that there exist two stable equilibria, the infection-free equilibrium and another positive equilibrium. As R0 changes when r2 is varying between 1.012 and 2.012, the backward bifurcation of model (11) occurs. Fig. 2 (b) indicates the bifurcation diagram, where three equilibria coexist when Rc < R0 < 1. Moreover, the infection-free equilibrium and the larger positive equilibrium are locally asymptotically stable, while the smaller positive equilibrium is unstable.

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

16

a

b

Fig. 2. The backward bifurcation of system (11) with q = 0 and a1 = 0. (a) Solutions either converge to a positive equilibrium or the infection-free equilibrium, dependent on initial values. (b) The backward bifurcation curve R0 versus I1 (the number of infected T-cells at the positive equilibrium). When R0 < Rc = 0.933, there exists a unique stable infection-free equilibrium (solid blue line). If Rc < R0 < 1, there exist a stable disease-free equilibrium, a stable positive equilibrium, and an unstable positive equilibrium. When R0 > 1, there exists a unique positive equilibrium which is stable, and the infection-free equilibrium becomes unstable.

4.2. p(a) is a step function Another well-accepted distribution for the viral shedding rate is the step function [28], taking the following form  p∗ , a ≥ a1 , p(a) = 0, a < a1 , where p∗ is the maximum production rate and a1 represents the time needed from the initial infection to shedding the first virus particle. In this case, the total number of virus particles (burst size) N produced by an infected T-cell during its lifespan can be computed as  N =

∞



∞

p(a)σ(a)da =

p∗ e−B0 a da =

a1

0

p∗ −B0 a1 e , B0

and the basic reproduction number is R0 =

kT0 N kT0 p∗ −B0 a1 = e , c cB0

with B0 = δ + r2 T0 /Tmax . Moreover, we have 

∞



0

∞

p∗ i(a, t)da =



∞

p∗ i(m + a1 , t)dm 0 a1   a1   ∞ T (t − a1 + s) + qI(t − a1 + s) ∗ =p i(m, t − a1 ) exp − (δ + r2 )ds dm Tmax 0 0   a1  ∞ T (t − a1 + s) + qI(t − a1 + s) (δ + r2 i(m, t − a1 )dm = p∗ exp − )ds Tmax 0 0   a1  T (t − a1 + s) + qI(t − a1 + s) = p∗ exp − (δ + r2 )ds I(t − a1 ), Tmax 0

p(a)i(a, t)da =

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

and therefore, system (6) can be reformulated into the following system:   dT (t) T (t) + qI(t) = λ − dT (t) + r1 T (t) 1 − − kV (t)T (t), dt Tmax   a1  T (t − a1 + s) + qI(t − a1 + s) dV (t) ∗ = p exp − (δ + r2 )ds I(t − a1 ) − cV (t), dt Tmax 0 dI(t) T (t) + qI(t) = kV (t)T (t) − δI(t) − r2 I(t) . dt Tmax

17

(13)

It is challenging to investigate the existence of positive equilibrium for system (13). Here, we focus on the special case q = 0 to check the possibility of coexistence of multiple positive equilibria. The positive equilibrium E2 = (T2 , V2 , I2 ) of system (13) with q = 0 can be calculated through the following equations:   T2 − kV2 T2 = 0, λ − dT2 + r1 T2 1 − Tmax     T2 p∗ exp − δ + r2 a1 I2 − cV2 = 0, Tmax T2 = 0. kV2 T2 − δI2 − r2 I2 Tmax T0 − Tλ0 : It is easy to get the following relationships by noting that r1 − d = Tr1max   T2     (T0 − T2 ) Tλ0 + Tr1max p∗ T2 I2 = > 0 and V2 = exp − δ + r2 a1 I2 > 0 T2 c Tmax δ + Tr2max

since T2 ∈ (0, T0 ). Then, it suffices to study the single equation for T2 since both I2 and V2 can be represented in terms of T2 written as   r2 T2 r2 T2 e(δ+ Tmax )a1 − kp∗ T2 = 0. c δ+ Tmax Let 

r2 y f3 (y) = c δ + Tmax



r2 y

e(δ+ Tmax )a1 − kp∗ y,

y ≥ 0.

Then we have f3 (0) = cδeδa1 > 0, and  2   r2 y r2 r2 y δ+ Tmax a1 ′′ ( ) f3 (y) = ca1 e 2+δ+ > 0. Tmax Tmax Noticing that R0 =

kT0 p∗ −B0 a1 , cB0 e

we get

  r2 T0 r2 T0 e(δ+ Tmax )a1 − kp∗ T0 f3 (T0 ) = c δ + Tmax = cB0 eB0 a1 − cR0 B0 eB0 a1 = cB0 eB0 a1 (1 − R0 ). If R0 > 1, then f3 (T0 ) < 0. Noting that f3 (0) > 0 and f3′′ (y) > 0, it can be concluded that there exists only one T2 ∈ (0, T0 ) such that f3 (T2 ) = 0, which further indicates the existence and uniqueness of the positive equilibrium for system (13) with q = 0. Therefore, we have the following results. Theorem 4.4. If R0 > 1, system (13) with q = 0 has only one positive equilibrium E2 .

18

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

In the case of R0 < 1, we have f3 (T0 ) > 0. Since f3 (0) > 0 and f3′′ (y) > 0, there are at most two positive solutions for f3 (y) = 0, which probably gives rise to backward bifurcation. Based on the properties of convex functions, the backward bifurcation may appear if there is a T2 ∈ (0, T0 ) such that f3′ (T2 ) = 0 and f3 (T2 ) < 0. Therefore, we have the following result. Theorem 4.5. System (13) with q = 0 exhibits backward bifurcation when R0 < 1 provided that the following conditions are satisfied: there exists a T2 ∈ (0, T0 ) such that f3′ (T2 ) = 0 and f3 (T2 ) < 0. 5. Conclusion and discussion We proposed an age-structured HIV model incorporating T-cell competition and viral production variation, which contains two ordinary differential equations and a hyperbolic equation to describe the age-structured infected target cells. The well-posedness of the model system was rigorously analyzed. Furthermore, the basic reproduction number was identified, and shown to determine the local stability of the infection-free equilibrium. However, it became challenging to perform further investigation for general viral reproducing function p(a) due to the involvement of the competition terms. Two specific forms of functions p(a) were further investigated: a delayed exponential form and a step function form. For the delayed exponential form, by introducing a new variable, we formulated a delay differential system, which is different from the previous work [28] where an ordinary differential system was formulated from an age-structured disease model with latency and relapse when the distribution function is exponential. Although the idea of this kind of reduction is similar, however, the involvement of a delay term requires an integral form evolving from a previous time instance. For these two special scenarios, the appearance of backward bifurcation was discussed. By the theoretical analysis and numerical simulations, we demonstrated that the model can reveal some reasonable and complicated dynamical behaviors. On one hand, we have confirmed that system has only one positive steady state when the basic reproduction number is greater than one. On the other hand, we have proven the existence of backward bifurcation when the basic reproduction number is less than one. This means that it is not enough to eliminate the HIV virus even though the basic reproduction number is below one. To describe variations in the viral production and the mortality of infected T-cells, age-structured models have been developed to study the within-host HIV dynamics [20–25]. These studies demonstrate forward bifurcations in terms of the basic reproduction number, however, backward bifurcations have been found in ordinary differential systems in the last two decades [19,34–39]. Recently, backward bifurcations were also illustrated in few age-structured models [26,40,41]. Comparing our results with those age-structured HIV models [21–25,29] without T-cell competition, we found that the competition rate between healthy T-cells and infected T-cells is one of the important reasons that cause backward bifurcation. Also, comparing our model with the competition HIV models [17–19] that do not include infection age to characterize the viral production variations, the dynamics are more complex and the analysis brings novel challenges. For ordinary differential systems, Castillo-Chavez and Song [36] have developed a rigorous approach for the occurrence of backward bifurcation by using center manifold method, and we may use the similar idea to verify the backward bifurcation for delay differential systems. In addition, for a delay differential system with T-cell competition, Hu et al. [18] have shown the existence of Hopf bifurcation by using the delay as a bifurcation parameter, and we further conjecture that the Hopf bifurcation occurs in our system. These two topics are worth further investigation. Acknowledgments YW acknowledges the financial support of RGC (PolyU 253004/14P) while visiting the Hong Kong Polytechnic University. YW is supported by NSFC (11301543, 11401589) and the Fundamental Research

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

19

Funds for the Central Universities (17CX02066). YL is supported by NSFC (11301442) and RGC (PolyU 253004/14P). The authors thank Dr. Hui Cao for helpful discussion. References [1] A.S. Perelson, P.W. Nelson, Mathematical analysis of HIV-1 dynamics in vivo, SIAM Rev. 41 (1) (1999) 3–44. [2] M.A. Nowak, R.M. May, Virus Dynamics: Mathematical Principles of Immunology and Virology, Oxford University Press, 2000. [3] W. Zhang, L.M. Wahl, P. Yu, Viral blips may not need a trigger: How transient viremia can arise in deterministic in-host models, SIAM Rev. 56 (1) (2014) 127–155. [4] D.D. Ho, A.U. Neumann, A.S. Perelson, W. Chen, J.M. Leonard, M. Markowitz, Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection, Nature 373 (1995) 123–126. [5] N. Sachsenberg, A.S. Perelson, S. Yerly, G.A. Schockmel, D. Leduc, B. Hirschel, L. Perrin, Turnover of CD4+ and CD8+ T lymphocytes in HIV-1 infection as measured by Ki-67 antigen, J. Exp. Med. 187 (1998) 1295–1303. [6] R.V. Culshaw, S. Ruan, A delay-differential equation model of HIV infection of CD4+ T-cells, Math. Biosci. 165 (2000) 27–39. [7] H.L. Smith, P. De Leenheer, Virus dynamics: A global analysis, SIAM J. Appl. Math. 63 (2003) 1313–1327. [8] L. Wang, M.Y. Li, Mathematical analysis of the global dynamics of a model for HIV infection of CD4+ T cells, Math. Biosci. 200 (2006) 44–57. [9] Y. Wang, Y. Zhou, J. Wu, J.M. Heffernan, Oscillatory viral dynamics in a delayed HIV pathogenesis model, Math. Biosci. 219 (2) (2009) 104–112. [10] M.Y. Li, H. Shu, Joint effects of mitosis and intracellular delay on viral dynamics: Two-parameter bifurcation analysis, J. Math. Biol. 64 (6) (2012) 1005–1020. [11] Y. Wang, Y. Zhou, F. Brauer, J.M. Heffernan, Viral dynamics model with CTL immune response incorporating antiretroviral therapy, J. Math. Biol. 67 (4) (2013) 901–934. [12] J.A. Borghans, L.S. Taams, M.H.M. Wauben, R.J. De Boer, Competition for antigenic sites during T cell proliferation: A mathematical interpretation of in vitro data, Proc. Natl. Acad. Sci. USA 96 (19) (1999) 10782–10787. [13] N.A. Gillet, N. Malani, A. Melamed, N. Gormley, R. Carter, D. Bentley, C. Berry, F.D. Bushman, G.P. Taylor, C.R.M. Bangham, The host genomic environment of the provirus determines the abundance of HTLV-1-infected T-cell clones, Blood 117 (11) (2011) 3113–3122. [14] J. Fang, S.A. Gourley, Y. Lou, Stage-structured models of intra-and inter-specific competition within age classes, J. Differential Equations 260 (2) (2016) 1918–1953. [15] L. Korobenko, E. Braverman, On logistic models with a carrying capacity dependent diffusion: Stability of equilibria and coexistence with a regularly diffusing population, Nonlinear Anal. RWA 13 (6) (2012) 2648–2658. [16] S. Liu, L. Chen, G. Luo, Y. Jiang, Asymptotic behaviors of competitive Lotka–Volterra system with stage structure, J. Math. Anal. Appl. 271 (1) (2002) 124–138. [17] L. Wang, S.F. Ellermeyer, HIV infection and CD4+ T cell dynamics, Discrete Contin. Dyn. Syst. Ser. B 6 (6) (2006) 1417–1430. [18] Z. Hu, X. Liu, H. Wang, W. Ma, Analysis of the dynamics of a delayed HIV pathogenesis model, J. Comput. Appl. Math. 234 (2) (2010) 461–476. [19] M.Y. Li, L. Wang, Backward bifurcation in a mathematical model for HIV infection in vivo with anti-retroviral treatment, Nonlinear Anal. RWA 17 (2014) 147–160. [20] A. Alshorman, C. Samarasinghe, W. Lu, L. Rong, An HIV model with age-structured latently infected cells, J. Biol. Dyn. (2016) 1–24. [21] C.J. Browne, S.S. Pilyugin, Global analysis of age-structured within-host virus model, Discrete Contin. Dyn. Syst. Ser. B 18 (8) (2013) 1999–2017. [22] G. Huang, X. Liu, Y. Takeuchi, Lyapunov functions and global stability for age-structured HIV infection model, SIAM J. Appl. Math. 72 (1) (2012) 25–38. [23] L. Rong, Z. Feng, A.S. Perelson, Mathematical analysis of age-structured HIV-1 dynamics with combination antiretroviral therapy, SIAM J. Appl. Math. 67 (3) (2007) 731–756. [24] J. Wang, J. Lang, X. Zou, Analysis of an age structured HIV infection model with virus-to-cell infection and cell-to-cell transmission, Nonlinear Anal. RWA 34 (2017) 75–96. [25] X. Wang, Y. Lou, X. Song, Age-structured within-host HIV dynamics with multiple target cells, Stud. Appl. Math. (2016) 1–34. [26] R. Qesmi, S. ElSaadany, J.M. Heffernan, J. Wu, A hepatitis B and C virus model with age since infection that exhibits backward bifurcation, SIAM J. Appl. Math. 71 (4) (2011) 1509–1530. [27] X. Lai, X. Zou, Dynamics of evolutionary competition between budding and lytic viral release strategies, Math. Biosci. Eng. 11 (5) (2014) 1091–1113. [28] P. Van Den Driessche, L. Wang, X. Zou, Modeling diseases with latency and relapse, Math. Biosci. Eng. 4 (2) (2007) 205. [29] P.W. Nelson, M.A. Gilchrist, D. Coombs, J.M. Hyman, A.S. Perelson, An age-structured model of HIV infection that allows for variations in the production rate of viral particles and the death rate of productively infected cells, Math. Biosci. Eng. 1 (2) (2004) 267–288. [30] M.V. Barbarossa, G. R¨ ost, Immuno-epidemiology of a population structured by immune status: a mathematical study of waning immunity and immune system boosting, J. Math. Biol. 71 (2015) 1737–1770.

20

Y. Wang et al. / Nonlinear Analysis: Real World Applications 38 (2017) 1–20

[31] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, Springer, New York, 1983. [32] L.J.S. Allen, F. Brauer, P. van den Driessche, J. Wu, Mathematical Epidemiology, Springer-Verlag, Berlin Heidelberg, 2008. [33] H.L. Smith, Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, American Mathematical Society, 1995. [34] K.P. Hadeler, P. van den Driessche, Backward bifurcation in epidemic control, Math. Biosci. 146 (1) (1997) 15–35. [35] P. van den Driessche, J. Watmough, A simple SIS epidemic model with a backward bifurcation, J. Math. Biol. 40 (6) (2000) 525–540. [36] C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. 1 (2) (2004) 361–404. [37] H. Gomez-Acevedo, M.Y. Li, Backward bifurcation in a model for HTLV-I infection of CD4+ T cells, Bull. Math. Biol. 67 (1) (2005) 101–114. [38] W. Wang, Backward bifurcation of an epidemic model with treatment, Math. Biosci. 201 (1) (2006) 58–71. [39] H. Cao, Y. Zhou, Z. Ma, Bifurcation analysis of a discrete SIS model with bilinear incidence depending on new infection, Math. Biosci. Eng. 10 (5–6) (2012) 1399–1417. [40] M. Martcheva, H.R. Thieme, Progression age enhanced backward bifurcation in an epidemic model with super-infection, J. Math. Biol. 46 (5) (2003) 385–424. [41] T.C. Reluga, J. Medlock, A.S. Perelson, Backward bifurcations and multiple equilibria in epidemic models with structured immunity, J. Theoret. Biol. 252 (1) (2008) 155–165.