GLOBAL STABILITY AND BACKWARD BIFURCATION OF A GENERAL VIRAL INFECTION MODEL WITH VIRUS-DRIVEN PROLIFERATION OF TARGET CELLS

In this paper, a general viral model with virus-driven proliferation of target cells is studied. Global stability results are established by employing the Lyapunov method and a geometric approach developed by Li and Muldowney. It is shown that under certain conditions, the model exhibits a global threshold dynamics, while if these conditions are not met, then backward bifurcation and bistability are possible. An example is presented to provide some insights on how the virus-driven proliferation of target cells influences the virus dynamics and the drug therapy strategies.


1.
Introduction. Many efforts have been recently made in mathematical modeling of in vivo dynamics of HIV-1, HBV, HCV and other viruses. Mathematical modeling plays an important role in understanding the viral infection process [25] as the resulting mathematical analysis and numerical simulations may shed light on developing new antiviral drugs and designing optimal combination of therapies for patients [10,22,28].
The standard in-host model is given by Nowak and May [24], which describes the interactions among healthy target cells (T cells), infected virus-producing T cells and free virus particles. The model takes the following form: where T (t) denotes the concentration of healthy T cells at time t, T i (t) denotes the concentration of infected virus-producing T cells at time t, and V (t) denotes the concentration of free virus particles at time t; the positive constant s is the recruitment rate of healthy T cells, and µ i > 0, i = 1, 2, 3 are, respectively, the mortality rates of healthy T cells, infected T cells and free virus particles; k 1 > 0 is the contact rate between healthy T cells and virus particles; N > 0 is the average number of virions budding out from a single infected target cell.
The standard model and its various modifications have been extensively studied in the literature (see [4,5,6,23,25,26,33,35,36] and the references therein). If the recruitment rate of healthy T cells is assumed to be a constant, then a threshold dynamics is usually observed: the viral load will be stabilized at the equilibrium level if the basic reproduction number is greater than unity, otherwise the virus will be cleared [5,12]. The viral load exhibits oscillations if a logistic growth function is adopted to model the production of new healthy T cells by antigen or mitogen [18,34]. The proliferation of healthy T cells in the presence of virus is often neglected in many in-host models. However, it is known that both CD8 + and CD4 + T cells specific to HIV can be directly stimulated, and that activated T cells can stimulate other CD8 + and CD4 + [11]. Taking this into account, we consider the following general model: , Here the function f (T, V ) denotes the intrinsic growth rate of the healthy T cells, which includes the source of new T cells from the thymus, the natural mortality of cells and the stimulation of T cells to proliferate in the presence of virus. The incidence of new infections of T cells occurs at a rate k 1 h(T )g(V ), which includes the rate of contact between free virions and healthy T cells as well as the probability of cell entry per contact. Infected cells die at the rate µ 2 p(T i ), and prior to death, each infected cell produces N new virions. The natural death rate of free virions is µ 3 q(V ) and the term k 2 h(T )g(V ) represents the loss of free infectious virions when they enter target cells, where k 2 is assumed to satisfy k 2 < N k 1 . This term is neglected (i.e., k 2 = 0) in [24,25] by arguing that it can be absorbed into the death rate µ 3 q(V ). However, as suggested in [5,9,10,27], this loss term may also play a role in HIV (HBV/HCV) dynamics. For our general model (2), we make the following assumptions: (H1): There exists a T > 0 such that f (T , 0) = 0 and f (T, 0)(T − T ) < 0 for T = T , T ≥ 0; (H2): f ∈ C 1 is bounded from above, f (0, V ) > 0 for all V ≥ 0, and ∂f (T, V ) ∂V < k 1 h(T )g (V ) for all T, V > 0; (3) Examples of functions satisfying the above assumptions include many widely used ones in the literature. For example, f (T, V ) = s − µT + rT (1 − T /T max ) with s, µ, r, T max > 0 [25]; f (T, V ) = s − µT + rT V /(C + V ) with s, µ, r, C > 0 [11], h(T ) = T m /(T m + A 1 ), and g(V ) = V n /(V n + A 2 ) with m, n, m, n > 0 and A 1 , A 2 ≥ 0 [5,17,21,25].
We organize the rest of this paper as follows. Section 2 provides some preliminary results concerning the well-posedness of Model (2) as well as an explicit formula for the basic reproduction number. Our global stability results are presented in Section 3. In Section 4, we show that backward bifurcation can be induced by the virusdriven proliferation of target cells. An example is presented in Section 5 to illustrate our analytical results and their biological implications. We conclude our work in Section 6.
2. Preliminary results. Proposition 1. Consider System (2) with the initial condition If the assumptions (H1)-(H3) are satisfied, then System (2) has a unique solution, which is nonnegative for t ≥ 0 and is ultimately bounded in R 3 + . Moreover, T (t) > 0 for all t ≥ 0.
Next we show the solution is ultimately bounded.
This, together with the first equation of (2), implies that lim sup t→∞ T (t) ≤ T . Set s = max T,V ≥0 f (T, V ) > 0 and σ = min{s/T , µ 2 k p /2, µ 3 k q }. It follows from (2) that Therefore, the solution (T (t), T i (t), V (t)) is ultimately bounded in R 3 + . The above analysis shows that the omega limit sets of System (2) are contained in the following positively invariant region:

HONGYING SHU AND LIN WANG
It is clear that System (2) admits a unique infection free equilibrium (IFE) E 0 = (T , 0, 0). An infection persistent equilibrium (IPE) E * = (T * , T * i , V * ) is a positive solution to the following equations Making use of the next generation method developed in [31], we obtain the basic reproduction number R 0 for System (2): By L Hôpital s Rule, R 0 can be written as (2). We first show that the IFE is locally asymptotically stable if R 0 < 1. The characteristic equation associated with the linearization of System (2) at the IFE E 0 is
Lemma 3.1. Assume that (H1)-(H3) are satisfied. If R 0 < 1, then E 0 is locally asymptotically stable; If R 0 > 1, then E 0 is a saddle point with a two-dimensional stable manifold and a one-dimensional unstable manifold.
In addition, by constructing a suitable Lyapunov function, we can show that E 0 is globally asymptotically stable in Γ if R 0 ≤ 1.
Theorem 3.2. Assume that (H1)-(H3) are satisfied. If R 0 ≤ 1, then the IFE E 0 of System (2) is globally asymptotically stable in Γ; If R 0 > 1, then E 0 is unstable, all solutions starting on the invariant T-axis approach E 0 along the T-axis, and System (2) is uniformly persistent in Γ\{T-axis}.
Proof. Define a Lyapunov function L as Then the derivative of L along the solutions of (2) is [8,13] implies that all solutions in Γ converge to E 0 . Thus E 0 is globally asymptotically stable.
Hence the solution leaves this neighborhood. Solutions on the positively invariant T-axis satisfy T (t) = f (T, 0). Thus T (t) → T as t → ∞, that is all solutions converge to E 0 along the T-axis. The uniform persistence result in [7] and a similar argument as in the proof of Proposition 3.2 in [14] lead to the instability of E 0 and the uniform persistence of System (2).
Note that the uniform persistence, together with uniform boundedness of solutions in • Γ, implies the existence of the IPE of System (2) ([1, Theorem 2.8.6]). To obtain the uniqueness of IPE E * = (T * , T * i , V * ), we need another assumption Proof. We only need to show that all solutions enter Ω = {(T, T i , V ) ∈ R 3 + : T ≤ T } in finite time. Suppose to the contrary that there exist a t 1 > 0 and a solution T (t) with T (t) > T for all t > t 1 . It then follows from (5) that which shows that T (t) → −∞ as t → ∞. This is a contradiction. Thus, there must exist a t 2 > t 1 such that T (t 2 ) ≤ T . Similar arguments can be used to show that all solutions (T, Next, we establish the uniqueness of the IPE. Define is strictly decreasing in T . Therefore, the equation F (T, V ) = 0 can be uniquely solved for T ∈ (0, T ) as a function of V . That is, there exists a function Note that ∂F (T, V )/∂V < 0 and ∂F (T, V )/∂T < 0 for all T, V > 0. Then which is strictly increasing since h is strictly increasing and q/g is nondecreasing.
By the monotonicity of φ and ϕ, and note that lim .

It then follows from
Next, we investigate the stability of the IPE E * . The characteristic equation associated with the linearization of System (2) at E * is The Routh-Hurwitz stability criterion implies that all roots of the characteristic equation (12) have negative real parts, and we have the following result.
The assumptions (H2) and (H3) imply that System (2) is competitive with respect to the cone {(x, y, z) ∈ R 3 : x, z ≥ 0, y ≤ 0} [30]. It then follows from the uniform persistence of System (2) (Theorem 3.2) that the omega limit set of a solution with initial condition satisfying (4) cannot contain a point on the Taxis if R 0 > 1. By Theorem 3.4, there is only one equilibrium E * which is not on the T -axis. The generalized Poincaré-Bendixson theorem for three dimensional competitive systems implies that the omega limit set of every solution with initial condition (4) either contains the unique IPE E * or is a nontrivial periodic orbit.
Next we adopt the geometric approach developed by Li and Muldowney [15,16] to establish the global stability of the IPE E * . To this end, we first briefly introduce this approach. Consider a system of ordinary differential equations where U is an open set in R n . Denote by x(t, x 0 ) the solution of (13) satisfying where ∂F ∂x [2] is the second additive compound matrix of the Jacobian matrix ∂F ∂x [20]. Let | · | denote a vector norm in R n and the induced matrix norm in R n×n .
where K denotes a compact absorbing set in U , Q F obtained by replacing each entry q ij of Q with its derivative in the direction of F. The following result from [15,16] will be employed to establish the global stability of E * when R 0 > 1.
Lemma 3.5. Let U ⊂ R n be a simply connected region. Assume that System (13) has a compact absorbing set K ∈ U admitting a unique equilibriumx ∈ U . Then the unique equilibriumx is globally asymptotically stable in U if there exist a function Q(x) and a Lozinskiȋ measure ξ such thatq 2 < 0. Clearly, • Γ is a simply connected region in R 3 . By Theorem 3.4, System (2) has (2), together with the boundedness of solutions given in Proposition 1, implies the existence of a compact absorbing set K ⊂ • Γ [2,32]. Therefore, by Lemma 3.5, we only need to show that there exist a matrix-valued function Q(x) and a Lozinskiȋ measure ξ such thatq 2 < 0.
Let (T (t), T i (t), V (t)) be any solution starting in the compact absorbing set K ⊂ • Γ and lett be sufficiently large such that (T (t), T i (t), V (t)) ∈ K for all t ≥t. Then along each solution (T (t), T i (t), V (t)) with (T (0), T i (0), V (0)) ∈ K, for t >t, we have Denote By (H4) and ∂f (T, V )/∂V ≥ 0, we further have Consequently, if then by the boundedness of T i and the uniform persistence of System (2), we obtain q 2 < 0. The above analysis and Lemma 3.3 immediately give the following result. Remark 1. System (2) includes the models considered in [5] and [29] as special cases: if f (T, V ) = f (T ) and h(x) = g(x) = p(x) = q(x) = x, then System (2) reduces to the model considered in [5]; if f (T, V ) = s − µ 1 T + rT V /(C + V ), and h(x) = g(x) = p(x) = q(x) = x, then System (2) reduces to the model considered in [29].
If k 2 = 0 in System (2), then the global stability of the IPE E * can be proved by constructing a Lyapunov function.
Proof. The existence and uniqueness of E * and its local stability follow from Theorem 3.4. By Lemma 3.3, it suffices to show that E * is globally attractive in • Γ. Define a Lyapunov function L in Γ as Calculating the time derivative of L along the positive solutions of System (2), we obtain .

Backward bifurcation induced by virus-driven proliferation of target cells.
In this section, we restrict our attention to System (2) with h(x) = g(x) = p(x) = q(x) = x. Then System (2) reduces to The basic reproduction number R 0 defined in (7) reads as If R 0 > 1, assumptions (H1)-(H2) ensure the existence and uniqueness of the IPE E * . Applying Theorems 3.2 and 3.6 to System (24), we have the following result.
Corollary 1. Assume that f (T, V ) satisfies (H1)-(H2). If R 0 ≤ 1, then the IFE E 0 of System (24) is globally asymptotically stable in R 3 + . If R 0 > 1, then E 0 is unstable, System (24) is uniformly persistent in R 3 + \{T-axis}, and there exists a unique IPE E * . In addition, if f (T, V ) satisfies (H4), ∂f (T, V )/∂V ≥ 0 for all T, V ≥ 0 and k 2 m 2 < µ 2 η, where η and m 2 are defined in (18) and (20), respectively, then E * is globally asymptotically stable in int(R 3 + ). The above result shows that if (3) holds, then the model exhibits the same threshold dynamics as the standard in-host model does. This implies that the virusdriven proliferation of target cells does not affect the viral dynamics structurally. In the sequel we assume that (3) is not satisfied. As we shall show that System (24) undergoes a backward bifurcation, and thus the virus-driven proliferation of target cells has an impact on the viral dynamics.
We first make the following assumptions on f (T, V ).
(B1): f ∈ C 1 is bounded from above and satisfies (H1) and (H4), We remark that (25) implies M ≥ T since f (T , 0) = 0 and f (T, V ) − k 1 T V is decreasing in T . Similar to the proofs of Proposition 1 and Lemma 3.3, we can show that the omega limit sets of System (24) are contained in the following positively invariant region If System (24) has an IPE E * = (T * , T * i , V * ), then we must have and T * i = k 1 T * V * /µ 2 . The existence of the IPE can be analyzed geometrically through examining the intersections of h 1 (V ) = f (T * , V ) and h 2 (V ) = k 1 T * V .
If R 0 > 1, that is T > T * , then f (T * , 0) > 0, and hence (26) has a unique positive root since f (T, V ) is concave down in V (Figure 1(v)). Summarizing the above argument, we have the following result concerning the existence of equilibria of System (24).
(ii) If R 0 < 1 and R 1 > 1, then there are three equilibria: the IFE E 0 and two IPEs, which we denote by If either R 0 < 1 and R 1 = 1 or R 0 ≥ 1, then there are two equilibria: the IFE E 0 and a unique IPE E * .
Next we study the stability of the equilibria of System (24). The characteristic equation of System (24) at E 0 is given as follows As a consequence of a straightforward analysis on the distribution of the eigenvalues, we have the following result.
Lemma 4.1. Assume that (B1) is satisfied. If R 0 < 1, then E 0 is locally asymptotically stable. If R 0 > 1, then E 0 is a saddle point with a two-dimensional stable manifold and a one-dimensional unstable manifold.
If R 0 < 1 and R 1 < 1, then it follows from (B2) and the definition of V that This, together with the first equation of (24) and (B1), implies that Similar as the proof of Theorem 3.2, we arrive at the following result. where Obviously, If R 0 < 1 and R 1 = 1, then there is a unique IPE E * = (T * , T * i , V * ) and ∂f (T * , V * )/∂V = k 1 T * . Thus, c 0 = 0, which implies that one eigenvalue is 0, the other two eigenvalues have negative real parts. Therefore, E * is locally stable, at which a saddle-node bifurcation occurs.
If R 0 ≥ 1, then there is a unique IPE E * = (T * , T * i , V * ) with ∂f (T * , V * )/∂V < k 1 T * . Thus, c 0 > 0, and hence E * is locally asymptotically stable. Moreover, the global stability of E * can be obtained from the same geometric approach as in the proof of Theorem 3.6.
Summarizing the discussion above, we have the following result. (i) If R 0 < 1 and R 1 < 1, then the IFE E 0 is globally asymptotically stable in R 3 + . (ii) If R 0 < 1 and R 1 = 1, then the IFE E 0 is locally asymptotically stable, and the unique IPE E * is locally stable, at which a saddle-node bifurcation occurs. (iii) If R 0 < 1 and R 1 > 1, then both the IFE E 0 and the IPE E * 2 are locally asymptotically stable, and the IPE E * 1 is a saddle point with a two-dimensional stable manifold and a one-dimensional unstable manifold. (iv) If R 0 = 1, then the IFE E 0 is locally stable, at which a transcritical bifurcation occurs, and the unique IPE E * is locally asymptotically stable. (v) If R 0 > 1, the IFE E 0 is unstable, and the unique IPE E * is locally asymptotically stable. In addition, all solutions in int(R 3 + ) with initial conditions satisfying (4) converge to E * provided that k 2 m 2 < µ 2 η, where η and m 2 are defined in (18) and (20), respectively.
5. An example: Data fitting, numerical simulations and drug therapies.
Then System (2) reduces to the following system: Here the term rT (t)V (t)/(c + V (t)) represents the virus-driven proliferation of T cells, r is the maximal proliferation rate satisfying 0 ≤ r < µ 1 , and c is the half saturation constant of the virus-driven proliferation process [11]. Clearly, System (31) admits a unique IFE E 0 = (s/µ 1 , 0, 0). If V * is a positive root of The basic reproduction number defined in (7) is given by Obviously, assumptions (H1) and (H3)-(H4) are satisfied with T = s/µ 1 . If r ≤ k 1 c, then (H2) holds. By Theorems 3.2 and 3.6, we have then System (31) admits a unique equilibrium, the IFE E 0 , which is globally asymptotically stable in R 3 + . (ii) If R 0 > 1, then E 0 becomes unstable, and there is a unique IPE E * , which is locally asymptotically stable. Moreover, if k 2 s ≤ µ 2 min{µ 1 − r, µ 2 } holds, then all solutions in int(R 3 + ) with initial conditions satisfying (4) converge to E * .
Next, we assume that r > k 1 c, then assumptions (B1)-(B2) hold with M = s/(µ 1 − r). If R 0 < 1, then we can easily calculate V = rc/k 1 − c, and obtain Note that it is easy to show that R 1 < 1 if and only if Applying Theorem 4.3 to System (31), we have the following result.
then E 0 is locally asymptotically stable, and the unique IPE E * is locally stable, at which a saddle-node bifurcation occurs.
Both E 0 and E * 2 are locally asymptotically stable, and E * 1 is a saddle point with a two-dimensional stable manifold and a one-dimensional unstable manifold. (iv) If R 0 = 1, then the IFE E 0 is locally stable, at which a transcritical bifurcation occurs, and there exists an IPE E * , which is locally asymptotically stable. (v) If R 0 > 1, then the IFE E 0 is unstable, and there exists a unique IPE E * , which is locally asymptotically stable. Moreover, if k 2 s ≤ µ 2 min{µ 1 − r, µ 2 } holds, then all solutions in int(R 3 + ) with initial conditions satisfying (4) converge to E * .
The above theorem implies that a backward bifurcation occurs at R 0 = 1 and bistability occurs when R 0 ∈ (1−a, 1). A bifurcation diagram is plotted in Figure 2.
To reduce the viral load, it is common to apply antiretroviral therapies using reverse transcriptase inhibitors (RTIs). Denote the efficacy of the RTIs by n rt , then the infection rate k 1 is reduced to k 1 = k 1 (1 − n rt ), and hence System (31) reads as As n rt increases, the effective infective rate k 1 decreases from k 1 , therefore, R 0 decreases, and there exists a critical value of n c1,rt at which R 0 = 1. In the meantime,   it is likely that the assumption (D2) holds and Theorem 5.2 applies. Thus, there exists another critical value of n c2,rt at which R 0 = 1 − a. Take parameter values as: s = 100 mm −3 · day −1 , µ 1 = 3.96 × 10 −2 day −1 , r = 0.03 day −1 , c = 10 ml −1 , k 1 = 3.0 × 10 −5 ml · day −1 , µ 2 = 0.3873 day −1 , N = 200 mm 3 · ml −1 , µ 3 = 4 day −1 , k 2 = 9.9933 × 10 −4 mm 3 · day −1 , then n c1,rt ≈ 0.5694 and n c2,rt = 0.75. As shown in Figure 4, we find R 0 > 1 for n rt ∈ [0, 0.5694) and 1 − a < R 0 < 1 for n rt ∈ (0.5694, 0.75) and R 0 < 1 − a for n rt ∈ (0.75, 1]. Due to the existence of a backward bifurcation and the resulting bistability for R 0 ∈ (1 − a, 1), the virus may not be completely eradicated even if R 0 is less than 1. Suppose that the drug therapy starts at day 200 at which the viral load of the  patient is near the equilibrium level, as indicated in Figure 5, the equilibrium viral load drops as the efficacy n rt increases, but the virus persists until the efficacy is increased to n c,rt = 0.68 at which the viral load quickly drops to zero. We notice that as n rt increases, the equilibrium viral load drops, and more damped oscillations are observed and as n rt reaches a critical value, n rt = n c,rt = 0.68, a complete clearance of virus is achieved. Note that such a critical value n c,rt is less than n c2,rt and the value is initial condition dependent and also depends on the time when the drug therapy starts. To completely clear the virus from the host, it is required that therapies either lower R 0 to be less than 1 − a or shift the viral load to the basin of attraction of the infection free equilibrium. 6. Summary and discussion. In this paper, we have incorporated the virusdriven proliferation of target cells into a general viral infection model (2). We have shown that this general model exhibits a threshold dynamics under the assumptions (H1)-(H3): if R 0 ≤ 1, then the infection free equilibrium E 0 is globally stable (Theorem 3.2), whereas if R 0 > 1, then E 0 becomes unstable and there exists a unique infection persistent equilibrium E * = (T * , T * i , V * ), which is locally stable provided that (H4) and ∂f (T * , V * )/∂V ≥ 0 are further satisfied (Theorem 3.4). If some additional conditions are satisfied, then E * attracts all positive solutions (Theorem 3.6). The more interesting dynamics occurs when the assumption (H2) is not satisfied. In that case, we have shown that a backward bifurcation occurs (Theorem 5.2). As such a backward bifurcation does not occur in the corresponding model where the virus-driven proliferation of target cells is ignored, we conclude that the backward bifurcation is induced by the virus-driven proliferation of target cells.
An example, System (31), has been given to demonstrate our analytical results. Essentially if r ≤ k 1 c and R 0 > 1, then the virus persists. Drug therapies such as the use of RTIs aiming to control the viral infection can reduce the value of the infection rate k 1 such that r > k 1 c and thus the result of Theorem 5.2 applies. As a consequence, a backward bifurcation occurs. The efficacy of the RTIs should be high enough such that either R 0 < 1 − a or the solution is shifted to the basin of attraction of the infection free equilibrium. Our analysis and numerical simulations suggest that simply reducing the basic reproduction number below unity may not be enough to successfully clear the viral infection and there exists an initial condition dependent critical value for the efficacy of drug therapy, above which the viral infection can be cleared. This implies that the virus-driven proliferation of target cells does play an important role and thus cannot be neglected.