Threshold dynamics and threshold analysis of HIV infection model with treatment

In this paper, a human immunodeficiency virus (HIV) infection model that includes a protease inhibitor (PI), two intracellular delays, and a general incidence function is derived from biologically natural assumptions. The global dynamical behavior of the model in terms of the basic reproduction number R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}_{0}$\end{document} is investigated by the methods of Lyapunov functional and limiting system. The infection-free equilibrium is globally asymptotically stable if R0≤1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}_{0}\leq 1$\end{document}. If R0>1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}_{0}>1$\end{document}, then the positive equilibrium is globally asymptotically stable. Finally, numerical simulations are performed to illustrate the main results and to analyze thre effects of time delays and the efficacy of the PI on R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}_{0}$\end{document}.


Introduction
Human immunodeficiency virus (HIV) is one of the most dangerous viruses that continues to be a major contributor to the global burden of disease [13]. The interactions between the human immune system and HIV have many highly complex characteristics, most of which are still not understood. Mathematical models have provided powerful tools for studying the dynamics of HIV infection, quantifying the effectiveness of drug therapies, and helping guide treatment strategies. A system of autonomous ordinary differential equations was developed to investigate the dynamics between HIV pathogens and the targeted CD4 + T-cells; see, e.g., [20][21][22] and the references therein. These models mainly explored the dynamics of uninfected target cells, infected target cells, and free virus particles. In [21,22], the effects of antiretroviral drugs including protease inhibitors (PIs) were considered, and Perelson and Nelson [22] showed that therapy using a single drug was doomed to fail because of drug resistance.
In 1996, Perelson et al. [23] assumed that there were two types of delays that occurred between the administration of drug and the observed decline in viral load: a pharmacological delay that occurred between the ingestion of drug and its appearance within cells and an intracellular delay that was between initial infection of a cell by HIV and the release of new virions. In the same year a discrete intracellular time delay defined as the time between infection of a cell and production of new virus particles was first introduced into an HIV infection model by Herz et al. [8]. Since then delayed HIV models have been studied by many authors, and we will not attempt to list all related references. Authors in [16] proposed a model for the interaction of HIV-1 with target cells including a gamma distribution of the intracellular time delay and assuming the density of non-infected target cells to be constant and studied the effects of protease inhibitor (PI) therapy. A model that allowed for less than perfect drug effects of a PI and which included a fixed delay in the initiation of virus production was considered in [18], where detailed analysis of the delay differential equation model was presented and the results to a model without delay were compared.
Culshaw and Ruan [3] developed a delay-differential equation model of HIV infection of CD4 + T-cells with discrete intracellular delay and bilinear incidence rate and investigated the effect of the time delay on local stability of the endemically infected equilibrium. Nelson and Perelson [19] generalized the basic ordinary differential model of HIV-1 infection with a PI by allowing the intracellular delay that varied according to a probability distribution and obtained the local stability results of two steady states. Based on the model built in [19], paper [12] further established global stability of steady states by constructing Lyapunov functionals. Li and Ma [9] constructed an HIV-1 infection model with the discrete intracellular delay and Holling type-II functional response term and gave sufficient criteria for local asymptotic stability of the infected equilibrium and global asymptotic stability of the viral free equilibrium. Wang et al. [28] considered a delayed HIV-1 infection model with Beddington-DeAngelis functional response and discussed the global stability of the infection-free equilibrium and the local stability of the chronic-infection equilibrium. The results of the global stability of the positive equilibrium obtained by Cai et al. [2] enriched and improved the corresponding results in [28]. In paper [2], it was shown that the global dynamics of the system restudied was completely determined by the basic reproduction number. Xu [30] investigated an HIV-1 infection model with the discrete intracellular delay and a saturation infection rate and established both the local stability and the global stability of the infection-free equilibrium and the infected equilibrium of the model. Pitchaimani and Divya [24] discussed local asymptotic stability of an HIV infection delay model in which the effects of a PI were considered and the values of three time delays varied according to the corresponding probability distribution. Bairagi and Adak [1] modified a basic ordinary differential HIV model with Hill type infection rate and the discrete intracellular time delay, studied local stability and global stability of both the delayed and non-delayed systems, and showed that multi-blockers drug therapy is more appropriate in the treatment of HIV patients in comparison to any mono-blocker drug therapy by numerical simulating. An epidemic model with free-living pathogens was constructed by Xing et al. [29] to understand the impact of free-living pathogens on the epidemics.
More general biological models take into account also the diffusion laws of the involved densities, and the related mathematical formulation is expressed by PDEs or systems of coupled PDEs [10,27].
Recently, Guo and Ma [6] first proposed a class of delay differential equations model of HIV infection dynamics with the intracellular delay, apoptosis induced by infected cells, and a general incidence function, and then analyzed the global properties of the model.
Considering the effects of a PI on the dynamics of HIV-1 infection, authors in [17] formu-lated and studied a model with a PI therapy and three delays as follows:  [4], and τ 2 is the delay representing the time necessary for a newly infected virus to become mature and then infectious. The factor e -mτ is the probability that an infected cell survives the interval τ , where 1/m is the average lifetime of infected cells before they become productive, and so the term e -mτ 1 [17]. The probability of survival of immature virions is given by e -vτ 2 , where 1/v is the average lifetime of an immature virus.
All the parameters in system (1.1) are positive constants.
In system (1.1), the time delays τ and τ 2 are necessary. But the delay τ 1 associated with the loss of target cells by infection should be removed. According to [4], when the PI alone is used, every infected cell still is destined to produce virus, and so the delay τ 1 associated with the loss of target cells by infection vanishes. Hence no delay is associated with the loss of target cells by infection in system (1.1) in fact because PI action does not prevent infected cells from becoming productive. In addition, the mass action infection law has been challenged in various ways since it has some unrealistic features such as the number of newly infected cells produced by a single virus depends on T and becomes very high when T is large [1]. Therefore, in order to improve understanding of the mechanisms of HIV infection model with a PI, we derive a refined model with more general infection processes: Here, a PI mono-therapy is considered; the PI is assumed to have variable efficacy and the drug is less than completely efficacious, i.e., p ∈ (0, 1); the T cells are allowed to vary; the function f (T(t), V I (t)) under some prescribed condition indicates the rate of uninfected target cells becoming infected by the infectious HIV viral particles. τ IP and τ IM correspond to, respectively, the time delay between the time a cell becomes infected and the time the infected cell develops into infectious, and the time required for newly infected virus particles to mature. All the parameters and variables appearing in system (1.2) have completely the same biological meanings as those of system (1.1).
In the following, we assume that the incidence function f (T, V I ) is continuously differential in R 2 + and satisfies the following hypotheses [6,14]: . These assumptions are biologically motivated. One could see that (H1) and (H2) are reasonable. (H3) follows the biological fact: if the total number of virus is constant, the bigger the amount of cells is, the higher the average number of cells which are infected by each virus in the unit time will be. The biological significance of (H4) is: if the total number of cells is constant, the bigger the amount of virus is, the higher the number of cells which are infected in the unit time will be. And we note that the final assumption ensures that ≤ 0 for all positive V I . This has the biological interpretation that the per capita dependence of new infections on the number of infectious virus particles V I is a decreasing function of V I . It is possible to check that the incidence rate f (T, V I ) generalizes many common forms such as βTV I , βTV I T+V I , βTV I 1+αV I , βTV I 1+γ T+αV I , and βTV I 1+γ T+αV I +αγ TV I , with α, β, γ > 0. In the present article, our primary goal is to establish some threshold dynamical results of system (1.2) in terms of the basic reproduction ratio. The basic reproduction number R 0 for the viral infection is easily obtained. However, Lyapunov functionals constructed for the whole model system of four differential equations are always imbalance of energy in and out. So we first use the method of Lyapunov functional and LaSalle's invariance principle for the first three equations of the model to get the global asymptotic stability of the equilibria for the subsystem; then we use the method of limiting system for the last delay differential equation; and finally we obtain the global asymptotic stability of the equilibria for the entire model. The infection-free equilibrium is globally asymptotically stable, and the HIV viruses are cleared if R 0 ≤ 1. If R 0 > 1, then the positive equilibrium is globally asymptotically stable and the infection persists. Furthermore, some interesting numerical simulations are performed to illustrate the main results on one side, and to explore the basic reproduction number R 0 on the other side. Within the latter, more specifically, considering the time lag τ IP and the efficacy p of the PI as parameters, we evaluate the lowest value of p which can be such that R 0 is less than one by sensitivity analysis of the basic reproduction number, which show that HIV can be permanently controlled by a PI mono-therapy only if the PI can keep nearly perfect all the time.
The paper is organized as follows. In Sect. 2, we present some preliminary results such as the global existence, nonnegativity, and boundedness of the solutions of system (1.2), the threshold value R 0 , the existence of equilibria, and the local asymptotic stabilities of the equilibria. In Sect. 3, by constructing suitable Lyapunov functionals and using established theories on asymptotically autonomous systems, we obtain the global threshold dynamics of the model. In Sect. 4, numerical simulations are carried out to demonstrate the main results and to show some sensitivity analysis on R 0 . Finally, some biological implications of our results can be found in the conclusions section.

Preliminary results
The preliminary results section includes three subsections. The properties about the wellposedness of system (1.2) are established in Sect. 2.1. The threshold value R 0 and the existence of equilibria are discussed in Sect. 2.2, and the local asymptotic stabilities of the infection-free equilibrium and the positive equilibrium are studied in Sect. 2.3.
Proof Obviously, Theorem 3.1 and Remark 3.3 in [26] imply that system (1.2) with the initial conditions above has a unique solution (T(t), T * (t), V I (t), V NI (t)) defined on some interval [0, σ 0 ) where σ 0 > 0. Let us first verify that this solution has nonnegative components. If Calculating the time derivative of A(t) along with the solution of (1.2) for t ∈ [0, σ 0 ), we get where g = min{d T , δ}. By the nonnegativity of T(t) and T * (t), and applying the well-known comparison principle to the first equation of (1.2) and (2.2), we know that T(t) and From the third equation and the fourth equation of (1.2), it follows that Hence, one can know that V I (t) and V NI (t) are both bounded on [0, σ 0 ). Therefore, by using Theorem 3.2 and Remark 3.3 of [26], the existence and uniqueness of the solution (T(t), T * (t), V I (t), V NI (t)) of (1.2) on [0, σ 0 ) can be extended to [0, +∞). Similarly, we can obtain that the solution (T(t), T * (t), V I (t), V NI (t)) of (1.2) with the initial conditions above is also nonnegative on [0, +∞). Now we show that (2.1) holds. Firstly, from the first equation of (1.2), one could verify that lim sup t→+∞ T(t) ≤ s d T . Secondly, we have from (2.2) that lim sup t→+∞ A(t) ≤ s ge mτ IP . From the nonnegativity of T(t) and T * (t), it follows that lim sup t→+∞ T * (t) ≤ s ge mτ IP . Lastly, again from the third equation and the fourth equation of (1.2), we get lim sup t→+∞ V I (t) ≤ If T(0) ≤ s d T , from the first equation of system (1.2), we have T(t) ≤ s d T for all t > 0. One can see that is a positively invariant region with respect to system (1.2).

Equilibria and basic reproductive number
We can compute the threshold value of model (1.2) as follows: This means that in order to have T ≥ 0 and T * > 0 at an equilibrium, we must have T * ∈ (0, s δe mτ IP ]. From (2.5) we obtain Positive equilibria of (1.2) are given by zeros of P in the interval (0, s/δe mτ IP ]. It is clear from R 0 > 1 that lim T * →0 + P T * = δe mτ IP (R 0 -1) > 0 and P s δe mτ IP = -δe mτ IP < 0.
Therefore, it follows from the continuity of the function P(T * ) on (0, s/δe mτ IP ] that there exists at least one T * ∈ (0, s/δe mτ IP ] such that P(T * ) = 0. Consequently, by (2.5), it follows that T > 0, V I > 0, and V NI > 0. So there exists at least one positive equilibrium E * . Taking derivative of P(T * ) yields According to (2.5), and notice that V I and from ∂f (T,V I )

Local stability
In this subsection, we explore local and global asymptotic stability of equilibria E 0 , E * . To study the local asymptotic stability, we analyze the roots of the characteristic equations associated with the linearization of system (1.2) at E 0 and E * .
First, we study the local stability of equilibria of system (1.2). For simplicity, τ IP + τ IM is denoted by τ .
Then we investigate the local stability of the positive equilibrium E * and have the following results.

Theorem 2.3
If R 0 > 1, then the positive equilibrium E * of system (1.2) is absolutely stable; that is, E * is locally asymptotically stable for all τ ≥ 0.
Proof The characteristic equation at E * reads (2.8) One of the characteristic roots of equation (2.8) is λ 1 = -c, and the equilibrium (T, T * , and Notice that Using the Routh-Hurwitz criterion, we know that all roots of equation ( Squaring and adding the two equations, (2.12) and (2.13), we can obtain (2.14) Since

Global threshold dynamics
Motivated by the work of Gourley and Lou [5], we establish the global asymptotic stability of E 0 and E * by constructing suitable Lyapunov functionals and using established theories on asymptotically autonomous systems in [15]. First, we show that when R 0 ≤ 1, E 0 is globally asymptotically stable. Note that the T(t) equation, the T * (t) equation, and the V I (t) equation are independent of the variable of non-infectious virus V NI (t). Thus we can first study the globally asymptotic stability of the infection-free equilibrium E 0 = (s/d T , 0, 0) of the following model: (3.1) Then, using established theories on asymptotically autonomous systems, we can obtain the globally asymptotical stability of equilibrium E 0 of system (1.2). Correspondingly, the positively invariant region for model (3.1) is Calculating the time derivative of U(t) along a solution of (3.1), we obtain Since T(t) ≤ s/d T and ∂f (T,V I ) is decreasing for V I . Therefore, Thus, we have dU(t) dt ≤ 0 whenever R 0 ≤ 1. Now we show that if dU(t) dt = 0, then either V I (t) = 0, or R 0 = 1 and T(t) = s/d T . It follows that the largest invariant set M 0 ⊆ M = {(T(t), T * (t), V I (t))| dU(t) dt = 0.} is the singleton {E 0 }. By the Lyapunov-LaSalle invariance principle [7], the infection-free equilibrium E 0 of (3.1) is globally asymptotically stable when R 0 ≤ 1. Then the V NI (t) equation of (1.2) can be considered as an asymptotically autonomous equation with the following limit equation: the solution of which tends to V NI . Thus, equilibrium E * of (1.2) is globally asymptotically stable if R 0 > 1. Now we complete the proof.

Numerical simulations
In this section we apply the main results obtained in the last section to three special infection functions of the generalized model (1.2), i.e., the bilinear incidence function investigated by Monica and Pitchaimani [17], the saturation response incidence rate studied by Xu [30], and the Beddington-DeAngelis function explored by Wang et al. [28]. We also present the results of sensitivity analysis on R 0 by numerical simulating, especially the estimated value of p which can such that R 0 < 1.

Simulation of the basic results
To verify the above analytic results, we introduce three specific infection functions to perform: bilinear incidence function βTV I which is proposed in [17], saturation response incidence rate βTV I 1+αV I in [30], and Beddington-DeAngelis function βTV I 1+γ T+αV I in [28]. In all of these papers [17,30], and [28], an HIV infection model with the intracellular delay was studied. Following [24], we fix the following coefficients: s = 5 mm -3 day -1 , d T = 0.03 day -1 , β = 0.002 mm 3 day -1 , α = 1 × 10 -6 mm 3 , γ = 0.5 × 10 -6 mm 3 , m = 0.03 day -1 , δ = 0.32 day -1 , N = 480, v = 0.03 day -1 , c = 3 day -1 . We will regard τ IP , τ IM , and p as parameters to verify the analytic results in Theorem 3.1 and Theorem 3.2. Notice that τ IP is expected to last between a few hours and three days and τ IM ≈ 0.15τ IP as mentioned in [4]. In this subsection, each graph has red dotted line, green solid line, and blue dashed line. They represent three solution curves with different initial values, respectively.
First, we consider the case of bilinear response incidence rate. The corresponding basic reproduction number is given by the infection-free equilibrium E 0 1 = (s/d T , 0, 0, 0), and the positive equilibrium Let p = 0.98, τ IP = 3, τ IM = 0.45, then R 01 = 0.9618 < 1 and the virus eventually dies out. Numerical simulation demonstrates the theoretical result of Theorem 3.1 that the infection-free equilibrium E 0 1 is globally asymptotically stable if R 01 < 1, as shown in Fig. 1. With the decrease of the parameter p , the basic reproduction number R 01 increases. So, we choose p = 0.9792, τ IP = 3, τ IM = 0.45, then R 01 = 1. According to Theorem 3.1, the infection-free equilibrium E 0 1 is globally asymptotically stable. Figure 2 shows that the infection-free equilibrium E 0 1 is globally asymptotically stable when R 01 = 1. Furthermore, we can verify the analytic result in Theorem 3.2. Direct calculation shows that R 01 = 18.9154 > 1 when p = 0.62, τ IP = 2, and τ IM = 0.3. Therefore, by Theorem 3.2, the positive equilibrium E * 1 is globally asymptotically stable. Figure 3 illustrates this fact.
Last, we consider the case of Beddington-DeAngelis functional response. The basic reproduction number is given by

Simulation of the basic reproduction number
Recall the expression of (2.3). Then it follows from the formula of R 0 that the basic reproduction number is decreasing with respect to p , τ IP , and τ IM . But we want to see how R 0 varies with the change of p , τ IP , and τ IM in quantity. R 01 (R 02 ) and R 03 are also selected to simulate. R 0 decreases linearly with increasing p and exponential decay with respect to τ IP . But as both m and v are very small parameters, the effect of τ IP changes on the value of R 0 is much less significant than that of p . From Fig. 10 we can get this conclusion. It is clear to see that, from Fig. 11, only a very small domain of the τ IP p -plane can be such that the basic reproduction numbers are less than unity. The domain is a very small particular area and is especially strict with p .
Furthermore, in order to find out exact minimum values of p which can be such that R 01 , R 02 , and R 03 are less than one, some numerical simulations are also performed. From the numerical simulation results, Fig. 12, it can be concluded that the lowest value we need is about 0.98.

Conclusions
Improved understanding of the mechanisms of latent infection and the importance of reservoirs of infection might eventually lead to a cure [13]. In this paper, we have carried out complete analysis for a delayed HIV infection model with a protease inhibitor (PI) mono-therapy and a general incidence function, system (1.2). The transmission dynamics of the system are studied and the global dynamics of system (1.2) are established. Based on the obtained theoretical results and simulation results, we conclude that the threshold (R 0 ) can provide much insight into the biological events underlying the disease process. If R 0 is less than or equal to unity, then the virus is cleared from the T-cell population and CD4 + T-cell count returns to normal, i.e., the virus is controlled and even cured and the disease dies out. If R 0 is greater than unity, then the infection persists in the host and the viral concentration remains at a constant level. The global stability result rules out the possibility of periodic oscillations and Hopf bifurcations. Moreover, the numerical simulation results demonstrate that the value of the basic reproductive number R 0 is dominated by the efficacy of the PI drug and is very weakly dependent on the intracellular delays. More importantly, due to the simulation results, we note that HIV can be permanently controlled by a PI mono-therapy only if the efficacy of the PI can keep reaching very high levels (at least 98%), namely the PI are nearly perfect all the time during therapy. However, it is nearly impossible to make the effectiveness of the drug always close to 100% because of the ability of the virus to mutate into a drug-resistant form and HIV-infected patients' dynamic health conditions. Therefore, the numerical simulation results of the discrete delay models for HIV infection with a PI drug suggest that therapy using a single PI drug is also doomed to fail compared with the analysis results of the ordinary differential equation models with a single drug in [22]. In order to eradicate HIV from an infected patient or help the body control the infection during therapy, as experience has borne out, changing different anti-HIV drugs such as PIs in time for HIV-positive individual patient and combining highly effective drugs for patients are needed.