An HIV model with age-structured latently infected cells

ABSTRACT

HIV latency remains a major obstacle to viral elimination. The activation rate of latently infected cells may depend on the age of latent infection. In this paper, we develop a model of HIV infection including age-structured latently infected cells. We mathematically analyse the model and use numerical simulations with different activation functions to show that the model can explain the persistence of low-level viremia and the latent reservoir stability in patients on therapy. Sensitivity tests suggest that the model is robust to the changes of most parameters but is sensitive to the relative magnitude of the net generation rate and the long-term activation rate of latently infected cells. To reduce the sensitivity, we extend the model to include homeostatic proliferation of latently infected cells. The new model is robust in reproducing the long-term dynamics of the virus and latently infected cells observed in patients receiving prolonged combination therapy.

KEYWORDS:

• Antiretroviral therapy
• virus dynamics
• latent reservoir
• age structure
• proliferation
• cellactivation

AMS SUBJECT CLASSIFICATION:

• 92B05
• 92D25
• 92C37
• 35B35
• 37N25

1. Introduction

HIV-1 replication in many infected individuals can be controlled by antiretroviral therapy consisting of three to five drugs. Although treatment can suppress the plasma viral load to below the detection limit of standard assays, that is, 50 viral RNA copies/ml, it cannot cure HIV-infected patients. The plasma viral load in patients on effective antiretroviral therapy has a multiphasic decay to below the detection limit. The initial decay phase is rapid with a half life of approximately one day, followed by a slower second-phase decay with a half life of 1–4 weeks [33], in which the plasma HIV RNA level approaches the detection limit. The third phase of the viral load decay below the detection limit is extremely slow with a half life of months or years [37].

Long-term low viral load persistence in patients on antiretroviral therapy is due to the virus persisting in cellular compartments or reservoirs. The best known reservoir consists of latently infected CD4+ T cells [3]. These resting memory CD4+ T cells are neither interfered with antiretroviral therapy nor affected by immune responses but can produce virus when they are activated by relevant antigens and therefore keep the virus from being eliminated. Because of the lack of sensitive testing methods measuring very low levels of viral RNA and latently infected cells in patients, the decay characteristics of the plasma viremia and the latent reservoir in the third phase remain unclear. The estimate of the half life of the latent reservoir had a range from 4.6 months to 44 months (see review in ref. [37]). In addition to being an obstacle to viral eradication, the existence of this stable viral reservoir has been found to have important clinical implications such as limited treatment options due to the accumulation of drug-resistant viral variants residing in the latent reservoir [45].

Mathematical models have been developed to study the decay dynamics of the latent reservoir. Muller et al. [29] introduced a simple model that considers the heterogeneity of latent cell activation assuming 100% effective combination therapy. Kim and Perelson [25] developed an ordinary differential equation model that includes bystander proliferation and decreasing activation of latently infected cells. To explain the stability of the latent reservoir, low viral load persistence, and occasional appearance of transient viral load measurements above the detection limit (the so-called viral blips), asymmetric cell division [38], programmed expansion and contraction [36], and homeostatic proliferation [36] of latently infected cells upon activation have been proposed and tested by mathematical models. The hypothesis of homeostatic proliferation of latently infected cells was later experimentally confirmed [7]. Other models, including stochastic models that study the dynamics of the latent reservoir or viral blips, can be found in refs. [1,2,9,16,24,43,50,55] (see reviews in [31,37]).

The major latent reservoir consists of resting memory CD4+ T cells. Antigen specificity plays an important role in the activation of these latently infected cells. Cells specific to frequently encountered antigens can be activated quickly while cells specific to rarely encountered antigens need more time to be activated. Thus, the more time elapsed since the establishment of latent infection, the more likely the latently infected cell is specific to a rarer antigen, and the less likely it is to be activated. The population of latently infected cells will gradually shift towards cells with a small chance of activation, that is, cells specific to increasingly rare antigens [47,48]. In this paper, we develop a model including latently infected cell activation on the basis of the models in refs. [25,36]. We incorporate the time elapsed since latent infection (i.e. the age of latent infection) into a latent infection model that was introduced by Perelson et al. [33] to study the second phase viral load decline during combination therapy. The activation rate of latently infected cells is assumed to depend on the age of latent infection. The basic reproductive number is derived and stability results of the infection-free and infected steady states are obtained. Different activation rate functions are used in simulations. We perform sensitivity tests on a number of parameters to determine the effect of parameters on the dynamics of viral load and the latent reservoir. We also modify the model by including homeostatic proliferation of latently infected cells to explore whether the model can robustly explain the persistence of low viral load and the latent reservoir stability in patients receiving prolonged antiretroviral therapy.

2. Model description and preliminary results

We develop a model including latently infected cell activation, as given by the following system of differential equations  (1(1) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -\mathrm{d}T-kVT\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =fkVT,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)kVT-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}V\left(t\right)& =N\delta T^{\ast }-cV.\end{array}$(1) ). In the model, $T\left(t\right)$ is the concentration of uninfected CD4+ T cells at time t, $L(a,t)$ denotes the concentration of latently infected T cells with latency age a (i.e. the time elapsed since the establishment of latency is a) at time t, $T^{\ast }\left(t\right)$ is the concentration of productively infected cells, and $V\left(t\right)$ represents the concentration of virions in plasma at t. The parameter λ is the production rate of uninfected CD4+ T cells, d is the per capita death rate of uninfected cells, and k is the infection rate of the target cell by virus. We assume that a small fraction (f) of infected cells lead to latency and that the remaining become productively infected cells. Productively infected cells die at a rate constant δ per cell. N is the viral burst size, which is the total number of virions released by one infected cell in its lifespan, and c is the viral clearance rate. (1) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -\mathrm{d}T-kVT\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =fkVT,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)kVT-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}V\left(t\right)& =N\delta T^{\ast }-cV.\end{array}$(1)

In addition to the contribution from new infection (fkV T), latently infected cells can also be generated by (bystander) proliferation [25]. Instead of parameterizing the proliferation rate and the death rate of latently infected cells separately, we lump these two parameters into one parameter, that is, $\gamma =\rho -{\delta }_L$, where ρ is the proliferation rate and ${\delta }_L$ is the death rate. Therefore, γ represents the net generation rate of latently infected cells.

The pool size of latent reservoir decreases when latently infected cells are activated to become productively infected cells. As shown in Strain et al. [47,48], the latently infected cell population is very likely to be heterogeneous: cells specific to frequently encountered antigens may be preferentially activated and quickly removed from the reservoir, leaving latently infected cells specific to rarely encountered antigens. Thus, the activation rate depends on the time elapsed since latent infection and we assume an age-dependent activation rate $\alpha \left(a\right)$ in the model. The integral of $\alpha \left(a\right)L(a,t)$ over age gives the total number of productively infected cells gained per unit time from the activation of latently infected cells.

Both the proliferation rate ρ and the death rate ${\delta }_L$ of latently infected cells can be assumed to be age-dependent without introducing additional difficulties in the analysis and simulation of age-structured models [14,39,41,44,49]. However, for simplicity we assume that they are constant and only focus on the age-dependent activation rate $\alpha \left(a\right)$. During proliferation of latently infected cells, a mother cell divides and gives rise to two daughter cells. These daughter cells are assumed to have the same age of being latently infected as the mother cell (note that a is the time elapsed since latent infection; it is not necessarily the biological age of an infected cell). Moreover, the reason why we are interested in the latent infection age is that the time since infection may affect the activation rate. The two daughter cells are initially in the same place as the mother cell and have the same chance as the mother cell (if not divided) to encounter with their specific antigen.

Two classes of antiretroviral drugs are mainly used in highly active antiretroviral therapy. One is the reverse transcriptase (RT) inhibitor, which can block the reverse transcription from viral RNA to DNA and thus reduces viral infection. The other is the protease inhibitor, which prevents HIV protease from cleaving the viral polyprotein into functional units, causing infected cells to produce immature virions that are non-infectious. Let ${ϵ}_{\mathrm{RT}}$ and ${ϵ}_{\mathrm{PI}}$ denote the efficacy of RT inhibitor and protease inhibitor, respectively ($0\le {ϵ}_{\mathrm{RT}},{ϵ}_{\mathrm{PI}}\le 1$). The value 0 represents no effect and 1 means 100% effective. To study the effect of protease inhibitor, we divide the virus population into two classes: infectious virus $V_I$ and non-infectious virus $V_{NI}$. The model with combination therapy is given by the following system of differential equations  (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ). (2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2)

Notice that model (1(1) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -\mathrm{d}T-kVT\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =fkVT,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)kVT-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}V\left(t\right)& =N\delta T^{\ast }-cV.\end{array}$(1) ) is a special case of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) when no drug treatment is given (i.e. ${ϵ}_{\mathrm{RT}}={ϵ}_{\mathrm{PI}}=0$). System (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) can be reformulated as a system of Volterra integral equations. For the ease of simplicity, we introduce the following notations: (3) $\begin{array}{rl}{K}_0\left(a\right)& ={\mathrm{e}}^{-{\int }_0^a\left(\alpha \right(s)-\gamma )\mathrm{d}s}\\ K_1\left(a\right)& =\alpha \left(a\right)K_0\left(a\right)\\ {\mathcal{K}}_i& ={\int }_0^{\mathrm{\infty }}{K}_i\left(a\right)\mathrm{d}a,i=0,1.\end{array}$(3) Because the latent reservoir is relatively stable or decays very slowly during treatment, we assume that $\gamma \le \underset{a}{inf}\alpha \left(a\right)$. $K_0\left(a\right)$ is an exponential probability density function (not normalized), representing the probability of an infected cell staying in the latent state until age a.

For convenience we introduce a new variable $B\left(t\right)$ to describe the rate at which an uninfected CD4+ T cell becomes infected at time t in the absence of drug therapy. (4) $B\left(t\right)=k{V}_I\left(t\right)T\left(t\right).$(4) Integrating the $L(a,t)$ equation in system (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) along the characteristic lines, $t-a=$ constant, we get the following solution (5) $L(a,t)=\left\{\begin{array}{ll}f(1-{ϵ}_{\mathrm{RT}})B(t-a)K_0\left(a\right)& \text{for}a<t,\\ L_0(a-t)\frac{K_0\left(a\right)}{K_0(a-t)}& \text{for}a\ge t.\end{array}\right.$(5) Substituting Equation (5(5) $L(a,t)=\left\{\begin{array}{ll}f(1-{ϵ}_{\mathrm{RT}})B(t-a)K_0\left(a\right)& \text{for}a<t,\\ L_0(a-t)\frac{K_0\left(a\right)}{K_0(a-t)}& \text{for}a\ge t.\end{array}\right.$(5) ) into the $T^{\ast }$ equation in (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ), we have (6) $\frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)=(1-f)(1-{ϵ}_{\mathrm{RT}})B\left(t\right)-\delta T^{\ast }\left(t\right)+{\int }_0^{t}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)B(t-a)\mathrm{d}a+\tilde{F}\left(t\right),$(6) where (7) $\tilde{F}\left(t\right)={\int }_t^{\mathrm{\infty }}\alpha \left(a\right)L_0(a-t)\frac{K_0\left(a\right)}{K_0(a-t)}\mathrm{d}a.$(7) It is clear that $\tilde{F}\left(t\right)\to 0$ as $t\to \mathrm{\infty }$. Integrating the $T^{\ast }$ equation in (6(6) $\frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)=(1-f)(1-{ϵ}_{\mathrm{RT}})B\left(t\right)-\delta T^{\ast }\left(t\right)+{\int }_0^{t}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)B(t-a)\mathrm{d}a+\tilde{F}\left(t\right),$(6) ) and changing the order of integration, we have (8) $\begin{array}{rl}{T}^{\ast }\left(t\right)& =T_0^{\ast }{\mathrm{e}}^{-\delta t}+{\int }_0^t{\mathrm{e}}^{-\delta (t-u)}\left[(1-f)(1-{ϵ}_{\mathrm{RT}})B\left(u\right)\phantom{{\int }_0^u}\right.\\ & \left.+{\int }_0^{u}f(1-{ϵ}_{\mathrm{RT}})B(u-\tau )K_1\left(\tau \right)\mathrm{d}\tau +\tilde{F}\left(u\right)\right]\mathrm{d}u\\ & ={\int }_0^t\left[\right(1-f){\mathrm{e}}^{-\delta (t-u)}+fH(t-u\left)\right](1-{ϵ}_{\mathrm{RT}})B\left(u\right)\mathrm{d}u+F\left(t\right),\end{array}$(8) where (9) $\begin{array}{rl}H\left(t\right)& ={\mathrm{e}}^{-\delta t}{\int }_0^t{\mathrm{e}}^{\delta x}{K}_1\left(x\right)\mathrm{d}x\\ F\left(t\right)& =T_0^{\ast }{\mathrm{e}}^{-\delta t}+{\int }_0^t{\mathrm{e}}^{-\delta (t-u)}\tilde{F}\left(u\right)\mathrm{d}u.\end{array}$(9) For the derivation of Equation (8(8) $\begin{array}{rl}{T}^{\ast }\left(t\right)& =T_0^{\ast }{\mathrm{e}}^{-\delta t}+{\int }_0^t{\mathrm{e}}^{-\delta (t-u)}\left[(1-f)(1-{ϵ}_{\mathrm{RT}})B\left(u\right)\phantom{{\int }_0^u}\right.\\ & \left.+{\int }_0^{u}f(1-{ϵ}_{\mathrm{RT}})B(u-\tau )K_1\left(\tau \right)\mathrm{d}\tau +\tilde{F}\left(u\right)\right]\mathrm{d}u\\ & ={\int }_0^t\left[\right(1-f){\mathrm{e}}^{-\delta (t-u)}+fH(t-u\left)\right](1-{ϵ}_{\mathrm{RT}})B\left(u\right)\mathrm{d}u+F\left(t\right),\end{array}$(8) ) we used the following calculation $\begin{array}{rl}{\int }_0^t{\mathrm{e}}^{-\delta (t-u)}{\int }_0^{u}B(u-\tau )K_1\left(\tau \right)\mathrm{d}\tau \mathrm{d}u& ={\int }_0^t{\mathrm{e}}^{-\delta (t-u)}{\int }_0^{u}B\left(\xi \right)K_1(u-\xi )\mathrm{d}\xi \mathrm{d}u\\ & ={\int }_0^t{\int }_{\xi }^t{\mathrm{e}}^{-\delta (t-u)}{K}_1(u-\xi )B\left(\xi \right)\mathrm{d}u\mathrm{d}\xi \\ & ={\int }_0^{t}B\left(\xi \right){\mathrm{e}}^{-\delta (t-\xi )}{\int }_0^{(t-\xi )}{\mathrm{e}}^{\delta x}{K}_1\left(x\right)\mathrm{d}x\mathrm{d}\xi \\ & ={\int }_0^{t}B\left(\xi \right)H(t-\xi )\mathrm{d}\xi .\end{array}$

Similarly, by integrating the T and $V_I$ equations in (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) we obtain (10) $T\left(t\right)=T_0{\mathrm{e}}^{-dt}+{\int }_0^t{\mathrm{e}}^{-d(t-u)}[\lambda -(1-{ϵ}_{\mathrm{RT}}\left)B\right(u\left)\right]\mathrm{d}u,$(10) (11) $V_I\left(t\right)=V_{I0}{\mathrm{e}}^{-ct}+{\int }_0^t{\mathrm{e}}^{-c(t-u)}(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(u\right)\mathrm{d}u.$(11)

Equations (8(8) $\begin{array}{rl}{T}^{\ast }\left(t\right)& =T_0^{\ast }{\mathrm{e}}^{-\delta t}+{\int }_0^t{\mathrm{e}}^{-\delta (t-u)}\left[(1-f)(1-{ϵ}_{\mathrm{RT}})B\left(u\right)\phantom{{\int }_0^u}\right.\\ & \left.+{\int }_0^{u}f(1-{ϵ}_{\mathrm{RT}})B(u-\tau )K_1\left(\tau \right)\mathrm{d}\tau +\tilde{F}\left(u\right)\right]\mathrm{d}u\\ & ={\int }_0^t\left[\right(1-f){\mathrm{e}}^{-\delta (t-u)}+fH(t-u\left)\right](1-{ϵ}_{\mathrm{RT}})B\left(u\right)\mathrm{d}u+F\left(t\right),\end{array}$(8) ), (10(10) $T\left(t\right)=T_0{\mathrm{e}}^{-dt}+{\int }_0^t{\mathrm{e}}^{-d(t-u)}[\lambda -(1-{ϵ}_{\mathrm{RT}}\left)B\right(u\left)\right]\mathrm{d}u,$(10) ) and (11(11) $V_I\left(t\right)=V_{I0}{\mathrm{e}}^{-ct}+{\int }_0^t{\mathrm{e}}^{-c(t-u)}(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(u\right)\mathrm{d}u.$(11) ), with $B\left(t\right)$ replaced by $k{V}_I\left(t\right)T\left(t\right)$, form a system of Volterra integral equations, which is equivalent to the original system (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ). Hence, for determining the existence and uniqueness of the solution we only need to consider the following system (12) $\begin{array}{rl}T\left(t\right)& ={\int }_0^t{\mathrm{e}}^{-d(t-u)}[\lambda -(1-{ϵ}_{\mathrm{RT}}\left)k{V}_I\right(u\left)T\right(u\left)\right]\mathrm{d}u+T_0{\mathrm{e}}^{-dt},\\ T^{\ast }\left(t\right)& ={\int }_0^t\left[\right(1-f){\mathrm{e}}^{-\delta (t-u)}+fH(t-u\left)\right](1-{ϵ}_{\mathrm{RT}})k{V}_I\left(u\right)T\left(u\right)\mathrm{d}u+F\left(t\right),\\ V_I\left(t\right)& ={\int }_0^t{\mathrm{e}}^{-c(t-u)}(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(u\right)\mathrm{d}u+V_{I0}{\mathrm{e}}^{-ct},\end{array}$(12) where $H\left(t\right)$ and $F\left(t\right)$ are given in Equation (9(9) $\begin{array}{rl}H\left(t\right)& ={\mathrm{e}}^{-\delta t}{\int }_0^t{\mathrm{e}}^{\delta x}{K}_1\left(x\right)\mathrm{d}x\\ F\left(t\right)& =T_0^{\ast }{\mathrm{e}}^{-\delta t}+{\int }_0^t{\mathrm{e}}^{-\delta (t-u)}\tilde{F}\left(u\right)\mathrm{d}u.\end{array}$(9) ).

3. Analysis of the system (2)

In this section we provide analytic results on the existence and uniqueness of the positive solution, calculate all possible steady states, and show their stability for the system (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) or the equivalent system (12(11) $V_I\left(t\right)=V_{I0}{\mathrm{e}}^{-ct}+{\int }_0^t{\mathrm{e}}^{-c(t-u)}(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(u\right)\mathrm{d}u.$(11) ).

3.1. Existence of positive solutions

Let $\mathrm{\Phi }\left(t\right)=\left(T\right(t),T^{\ast }(t),V_I(t){)}^{\mathrm{T}}$, where $\mathrm{T}$ represents the transpose of a vector. System (12(11) $V_I\left(t\right)=V_{I0}{\mathrm{e}}^{-ct}+{\int }_0^t{\mathrm{e}}^{-c(t-u)}(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(u\right)\mathrm{d}u.$(11) ) can be written as $\mathrm{\Phi }\left(t\right)={\int }_0^t\kappa (t-u)g\left(\mathrm{\Phi }\right(u\left)\right)\mathrm{d}u+f\left(t\right),$ where $f\left(t\right)=(T_0{\mathrm{e}}^{-dt},F(t),V_{I0}{\mathrm{e}}^{-ct}{)}^{\mathrm{T}}$ is a continuous function from $[0,\mathrm{\infty })$ to $[0,\mathrm{\infty }{)}^3$, κ is the following $3\times 3$ matrix with entries being locally integrable functions on $[0,\mathrm{\infty })$, $\kappa \left(t\right)=\left(\begin{array}{ccc}\lambda {\mathrm{e}}^{-dt}& 0& -(1-{ϵ}_{\mathrm{RT}}){\mathrm{e}}^{-dt}\\ 0& 0& (1-{ϵ}_{\mathrm{RT}})\left[\right(1-f){\mathrm{e}}^{-\delta t}+fH(t\left)\right]\\ 0& (1-{ϵ}_{\mathrm{PI}})N\delta {\mathrm{e}}^{-ct}& 0\end{array}\right)$ and g is defined by $g\left(\mathrm{\Phi }\right)=(1,T^{\ast },k{V}_{I}T{)}^{\mathrm{T}}.$ It is clear that $f\in C\left(\right[0,\mathrm{\infty });{\mathbf{R}}^3)$, $g\in C({\mathbf{R}}^3,{\mathbf{R}}^3)$, and $\kappa \in L_{loc}^1\left(\right[0,\mathrm{\infty });{\mathbf{R}}^{3\times 3})$. From Theorem 1.1 in Gripenberg et al. [21], Section 12.1, there is a unique solution of the system given an initial condition and the solutions depend continuously on the initial conditions.

To see whether all solutions remain non-negative for positive initial value, we use the following system (see Equation (6(6) $\frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)=(1-f)(1-{ϵ}_{\mathrm{RT}})B\left(t\right)-\delta T^{\ast }\left(t\right)+{\int }_0^{t}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)B(t-a)\mathrm{d}a+\tilde{F}\left(t\right),$(6) )) that is also equivalent to system (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ). (13) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right),\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)-\delta T^{\ast }\left(t\right)\\ & +{\int }_0^{t}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)k{V}_I(t-a)T(t-a)\mathrm{d}a+\tilde{F}\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right),\end{array}$(13) where $\tilde{F}$ is given in Equation (7(7) $\tilde{F}\left(t\right)={\int }_t^{\mathrm{\infty }}\alpha \left(a\right)L_0(a-t)\frac{K_0\left(a\right)}{K_0(a-t)}\mathrm{d}a.$(7) ) and $\tilde{F}\left(t\right)>0$, $\underset{t\to \mathrm{\infty }}{lim}\tilde{F}\left(t\right)=0$.

We proceed to prove that $T\left(t\right)$, $T^{\ast }\left(t\right)$ and $V_I\left(t\right)$ are all non-negative for all t>0 and for all non-negative initial value. Suppose that there exists a $\overline{t}>0$ such that $T^{\ast }\left(\overline{t}\right)=0$, $T\left(\overline{t}\right)\ge 0$, $V_I\left(\overline{t}\right)\ge 0$, and $T\left(t\right),T^{\ast }\left(t\right),V_I\left(t\right)>0$ for $0<t<\overline{t}$. Then from the $T^{\ast }$ equation in (13(12) $\begin{array}{rl}T\left(t\right)& ={\int }_0^t{\mathrm{e}}^{-d(t-u)}[\lambda -(1-{ϵ}_{\mathrm{RT}}\left)k{V}_I\right(u\left)T\right(u\left)\right]\mathrm{d}u+T_0{\mathrm{e}}^{-dt},\\ T^{\ast }\left(t\right)& ={\int }_0^t\left[\right(1-f){\mathrm{e}}^{-\delta (t-u)}+fH(t-u\left)\right](1-{ϵ}_{\mathrm{RT}})k{V}_I\left(u\right)T\left(u\right)\mathrm{d}u+F\left(t\right),\\ V_I\left(t\right)& ={\int }_0^t{\mathrm{e}}^{-c(t-u)}(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(u\right)\mathrm{d}u+V_{I0}{\mathrm{e}}^{-ct},\end{array}$(12) ) we have $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(\overline{t}\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(\overline{t}\right)T\left(\overline{t}\right)\\ & +{\int }_0^{\overline{t}}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)k{V}_I(\overline{t}-a)T(\overline{t}-a)\mathrm{d}a+\tilde{F}\left(\overline{t}\right)\\ & >0.\end{array}$ This shows that $T^{\ast }\left(t\right)$ will never become negative for t>0. Similarly, we can prove that $T\left(t\right)\ge 0$ and $V_I\left(t\right)\ge 0$ for all t>0 and for all non-negative initial value.

3.2. Steady states and stability

We use the system (13(12) $\begin{array}{rl}T\left(t\right)& ={\int }_0^t{\mathrm{e}}^{-d(t-u)}[\lambda -(1-{ϵ}_{\mathrm{RT}}\left)k{V}_I\right(u\left)T\right(u\left)\right]\mathrm{d}u+T_0{\mathrm{e}}^{-dt},\\ T^{\ast }\left(t\right)& ={\int }_0^t\left[\right(1-f){\mathrm{e}}^{-\delta (t-u)}+fH(t-u\left)\right](1-{ϵ}_{\mathrm{RT}})k{V}_I\left(u\right)T\left(u\right)\mathrm{d}u+F\left(t\right),\\ V_I\left(t\right)& ={\int }_0^t{\mathrm{e}}^{-c(t-u)}(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(u\right)\mathrm{d}u+V_{I0}{\mathrm{e}}^{-ct},\end{array}$(12) ) for our stability analysis. According to Gripenberg et al. [21], Section 15.1, any equilibrium of system (13(12) $\begin{array}{rl}T\left(t\right)& ={\int }_0^t{\mathrm{e}}^{-d(t-u)}[\lambda -(1-{ϵ}_{\mathrm{RT}}\left)k{V}_I\right(u\left)T\right(u\left)\right]\mathrm{d}u+T_0{\mathrm{e}}^{-dt},\\ T^{\ast }\left(t\right)& ={\int }_0^t\left[\right(1-f){\mathrm{e}}^{-\delta (t-u)}+fH(t-u\left)\right](1-{ϵ}_{\mathrm{RT}})k{V}_I\left(u\right)T\left(u\right)\mathrm{d}u+F\left(t\right),\\ V_I\left(t\right)& ={\int }_0^t{\mathrm{e}}^{-c(t-u)}(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(u\right)\mathrm{d}u+V_{I0}{\mathrm{e}}^{-ct},\end{array}$(12) ), if it exists, must be a constant solution of the limiting system. (14) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)-\delta T^{\ast }\left(t\right)\\ & +{\int }_0^{\mathrm{\infty }}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)k{V}_I(t-a)T(t-a)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right).\end{array}$(14)

System (14(13) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right),\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)-\delta T^{\ast }\left(t\right)\\ & +{\int }_0^{t}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)k{V}_I(t-a)T(t-a)\mathrm{d}a+\tilde{F}\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right),\end{array}$(13) ) has two constant solutions: the infection-free steady state $\overline{E}=(\overline{T},{\overline{T}}^{\ast },{\overline{V}}_I)=\left(\frac{\lambda }{d},0,0\right)$ and the infected steady state $\tilde{E}=(\tilde{T},{\tilde{T}}^{\ast },{\tilde{V}}_I),$ where (15) $\begin{array}{rl}\tilde{T}& =\frac{c}{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN(1-f+f{\mathcal{K}}_1)}\\ {\tilde{T}}^{\ast }& =\frac{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\lambda (1-f+f{\mathcal{K}}_1)-dc}{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\delta }\\ {\tilde{V}}_I& =\frac{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\lambda (1-f+f{\mathcal{K}}_1)-dc}{(1-{ϵ}_{\mathrm{RT}})kc}\end{array}$(15) with ${\mathcal{K}}_1$ given in Equation (3(3) $\begin{array}{rl}{K}_0\left(a\right)& ={\mathrm{e}}^{-{\int }_0^a\left(\alpha \right(s)-\gamma )\mathrm{d}s}\\ K_1\left(a\right)& =\alpha \left(a\right)K_0\left(a\right)\\ {\mathcal{K}}_i& ={\int }_0^{\mathrm{\infty }}{K}_i\left(a\right)\mathrm{d}a,i=0,1.\end{array}$(3) ).

It is clear that ${\tilde{V}}_I>0$ if and only if $(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\lambda (1-f+f{\mathcal{K}}_1)-dc>0$ or $\mathcal{R}>1$, where (16) $\mathcal{R}=\frac{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\lambda (1-f+f{\mathcal{K}}_1)}{dc}.$(16) The infected steady state (15(14) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)-\delta T^{\ast }\left(t\right)\\ & +{\int }_0^{\mathrm{\infty }}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)k{V}_I(t-a)T(t-a)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right).\end{array}$(14) ) is feasible if and only if $\mathcal{R}>1$. Notice that $\lambda /d$ is the target cell population in the absence of infection, k and c are the cell infection and viral clearance rate, respectively, 1−f is the fraction of productive infection that leads to viral production, and $f{\mathcal{K}}_1$ represents the contribution to productively infected cells from activation of latently infected cells. Thus, $\mathcal{R}$ gives the effective reproductive ratio of the virus under drug treatment.

We now study the stability of steady states. Using the linearization process [20], we obtain the following linearized system of (14(13) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right),\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)-\delta T^{\ast }\left(t\right)\\ & +{\int }_0^{t}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)k{V}_I(t-a)T(t-a)\mathrm{d}a+\tilde{F}\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right),\end{array}$(13) ) about any equilibrium, denoted by $\breve{E}=(\breve{T},{\breve{T}}^{\ast },{\breve{V}}_I)$. (17) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =[-d-(1-{ϵ}_{\mathrm{RT}}\left)k{\breve{V}}_I\right]T\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k\breve{T}{V}_I\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-{ϵ}_{\mathrm{RT}})k{\breve{V}}_I\left[(1-f)T\left(t\right)+{\int }_0^{\mathrm{\infty }}f{K}_1\left(a\right)T(t-a)\mathrm{d}a\right]-\delta T^{\ast }\left(t\right)\\ & +(1-{ϵ}_{\mathrm{RT}})k\breve{T}\left[(1-f)V_I\left(t\right)+{\int }_0^{\mathrm{\infty }}f{K}_1\left(a\right)V_I(t-a)\mathrm{d}a\right]\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right).\end{array}$(17) The corresponding characteristic equation is (18) $\left|\begin{array}{ccc}-d-(1-{ϵ}_{\mathrm{RT}})k{\breve{V}}_I-\zeta & 0& -(1-{ϵ}_{\mathrm{RT}})k\breve{T}\\ (1-{ϵ}_{\mathrm{RT}})k{\breve{V}}_I[1-f+f{\hat{K}}_1(\zeta \left)\right]& -\delta -\zeta & (1-{ϵ}_{\mathrm{RT}})k\breve{T}[1-f+f{\hat{K}}_1(\zeta \left)\right]\\ 0& (1-{ϵ}_{\mathrm{PI}})N\delta & -c-\zeta \end{array}\right|=0,$(18) where ζ is an eigenvalue and ${\hat{K}}_1\left(\zeta \right)$ denotes the Laplace transform of $K_1\left(a\right)$, that is ${\hat{K}}_1\left(\zeta \right)={\int }_0^{\mathrm{\infty }}{K}_1\left(a\right){\mathrm{e}}^{-\zeta a}\mathrm{d}a.$ We will use this characteristic equation to study the stability of steady states.

Let us first consider the infection-free steady-state $\overline{E}$. The following result suggests that the virus is predicted to die out if the effective reproductive ratio is less than 1.

Theorem 3.1.

The infection-free steady-state $\overline{E}$ is locally asymptotically stable if $\mathcal{R}<1,$ and it is unstable if $\mathcal{R}>1$.

Proof.

Substituting the steady-state $\overline{E}$ for the equilibrium point $\breve{E}$ into the characteristic Equation (18(17) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =[-d-(1-{ϵ}_{\mathrm{RT}}\left)k{\breve{V}}_I\right]T\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k\breve{T}{V}_I\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-{ϵ}_{\mathrm{RT}})k{\breve{V}}_I\left[(1-f)T\left(t\right)+{\int }_0^{\mathrm{\infty }}f{K}_1\left(a\right)T(t-a)\mathrm{d}a\right]-\delta T^{\ast }\left(t\right)\\ & +(1-{ϵ}_{\mathrm{RT}})k\breve{T}\left[(1-f)V_I\left(t\right)+{\int }_0^{\mathrm{\infty }}f{K}_1\left(a\right)V_I(t-a)\mathrm{d}a\right]\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right).\end{array}$(17) ), we have (19) $\left|\begin{array}{ccc}-d-\zeta & 0& -(1-{ϵ}_{\mathrm{RT}})k\lambda /d\\ 0& -\delta -\zeta & (1-{ϵ}_{\mathrm{RT}})k\lambda [1-f+f{\hat{K}}_1(\zeta \left)\right]/d\\ 0& (1-{ϵ}_{\mathrm{PI}})N\delta & -c-\zeta \end{array}\right|=0$(19) that is, (20) $(\zeta +d)\left[\right(\zeta +c\left)\right(\zeta +\delta )-(1-{ϵ}_{\mathrm{RT}}\left)\right(1-{ϵ}_{\mathrm{PI}}\left)kN\lambda \delta \right(1-f+f{\hat{K}}_1\left(\zeta \right)\left)/d\right]=0.$(20) Using the effective reproductive ratio $\mathcal{R}$ defined in Equation (16(15) $\begin{array}{rl}\tilde{T}& =\frac{c}{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN(1-f+f{\mathcal{K}}_1)}\\ {\tilde{T}}^{\ast }& =\frac{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\lambda (1-f+f{\mathcal{K}}_1)-dc}{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\delta }\\ {\tilde{V}}_I& =\frac{(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\lambda (1-f+f{\mathcal{K}}_1)-dc}{(1-{ϵ}_{\mathrm{RT}})kc}\end{array}$(15) ), the above characteristic equation can be rewritten as (21) $(\zeta +d)\left[(\zeta +c)(\zeta +\delta )-\mathcal{R}c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}\right]=0.$(21) One negative root of Equation (21(20) $(\zeta +d)\left[\right(\zeta +c\left)\right(\zeta +\delta )-(1-{ϵ}_{\mathrm{RT}}\left)\right(1-{ϵ}_{\mathrm{PI}}\left)kN\lambda \delta \right(1-f+f{\hat{K}}_1\left(\zeta \right)\left)/d\right]=0.$(20) ) is $\zeta =-d$ and all the other roots are determined by the equation (22) $(\zeta +c)(\zeta +\delta )=\mathcal{R}c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}.$(22) Note that $|{\hat{K}}_1\left(\zeta \right)|\le {\mathcal{K}}_1$ for all complex roots ζ with non-negative real parts. It follows that $|1-f+f{\hat{K}}_1\left(\zeta \right)|\le 1-f+f{\mathcal{K}}_1$. Hence, the modulus of the right-hand side of Equation (22(21) $(\zeta +d)\left[(\zeta +c)(\zeta +\delta )-\mathcal{R}c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}\right]=0.$(21) ) is less than $c\delta$ when $\mathcal{R}<1$. Since the modulus of the left-hand side of Equation (22(21) $(\zeta +d)\left[(\zeta +c)(\zeta +\delta )-\mathcal{R}c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}\right]=0.$(21) ) is always greater than $c\delta$ when the real part of ζ is non-negative, we conclude that all roots of (22(21) $(\zeta +d)\left[(\zeta +c)(\zeta +\delta )-\mathcal{R}c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}\right]=0.$(21) ) have negative real parts if $\mathcal{R}<1$. It follows that $\overline{E}$ is locally asymptotically stable when $\mathcal{R}<1$.

In the case of $\mathcal{R}>1$, let $\varphi \left(\zeta \right)=(\zeta +c)(\zeta +\delta )-\mathcal{R}c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}$ Thus, any real root of $\varphi \left(\zeta \right)=0$ is also the root of Equation (22(21) $(\zeta +d)\left[(\zeta +c)(\zeta +\delta )-\mathcal{R}c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}\right]=0.$(21) ). Recognizing that $\varphi \left(0\right)=c\delta (1-\mathcal{R})<0\text{and}\underset{\zeta \to \mathrm{\infty }}{lim}\varphi \left(\zeta \right)=\mathrm{\infty }$ we know $\varphi \left(\zeta \right)=0$ has at least one positive root ${\zeta }^{\ast }>0$, which is a positive eigenvalue of the characteristic equation  (20(19) $\left|\begin{array}{ccc}-d-\zeta & 0& -(1-{ϵ}_{\mathrm{RT}})k\lambda /d\\ 0& -\delta -\zeta & (1-{ϵ}_{\mathrm{RT}})k\lambda [1-f+f{\hat{K}}_1(\zeta \left)\right]/d\\ 0& (1-{ϵ}_{\mathrm{PI}})N\delta & -c-\zeta \end{array}\right|=0$(19) ). This shows that the infection-free steady state is unstable when $\mathcal{R}>1$. This finishes the proof of the theorem.

Next, we consider the stability of the infected steady-state $\tilde{E}$. As noted earlier, this steady state exists if and only if $\mathcal{R}>1$. The following result suggests that the virus infection will be established if the effective reproductive ratio is greater than 1.

Theorem 3.2.

The infected steady-state $\tilde{E}$ is locally asymptotically stable if $\mathcal{R}>1$.

Proof.

Substituting the infected steady-state $\tilde{E}=(\tilde{T},{\tilde{T}}^{\ast },{\tilde{V}}_I)$ for the equilibrium point $\breve{E}$ into the characteristic equation  (18(17) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =[-d-(1-{ϵ}_{\mathrm{RT}}\left)k{\breve{V}}_I\right]T\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k\breve{T}{V}_I\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-{ϵ}_{\mathrm{RT}})k{\breve{V}}_I\left[(1-f)T\left(t\right)+{\int }_0^{\mathrm{\infty }}f{K}_1\left(a\right)T(t-a)\mathrm{d}a\right]-\delta T^{\ast }\left(t\right)\\ & +(1-{ϵ}_{\mathrm{RT}})k\breve{T}\left[(1-f)V_I\left(t\right)+{\int }_0^{\mathrm{\infty }}f{K}_1\left(a\right)V_I(t-a)\mathrm{d}a\right]\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right).\end{array}$(17) ), we have the following simplified equation after tedious calculation (23) $\begin{array}{rl}& (\zeta +c)(\zeta +\delta )[\zeta +d+(1-{ϵ}_{\mathrm{RT}}\left)k{\hat{V}}_I\right]\\ & =(\zeta +d)(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\delta \hat{T}(1-f+f{\hat{K}}_1(\zeta \left)\right).\end{array}$(23) Notice that the steady state of infected virus in Equation (15(14) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)-\delta T^{\ast }\left(t\right)\\ & +{\int }_0^{\mathrm{\infty }}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)k{V}_I(t-a)T(t-a)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right).\end{array}$(14) ), ${\tilde{V}}_I$, can be written as the following form in terms of the effective reproductive ratio $\mathcal{R}$ (24) ${\tilde{V}}_I=(\mathcal{R}-1)\frac{d}{(1-{ϵ}_{\mathrm{RT}})k}.$(24) Substituting ${\tilde{V}}_I$ in Equation (24(23) $\begin{array}{rl}& (\zeta +c)(\zeta +\delta )[\zeta +d+(1-{ϵ}_{\mathrm{RT}}\left)k{\hat{V}}_I\right]\\ & =(\zeta +d)(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})kN\delta \hat{T}(1-f+f{\hat{K}}_1(\zeta \left)\right).\end{array}$(23) ) and $\tilde{T}$ in Equation (15(14) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT\left(t\right)-(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_I\left(t\right)T\left(t\right)-\delta T^{\ast }\left(t\right)\\ & +{\int }_0^{\mathrm{\infty }}f(1-{ϵ}_{\mathrm{RT}})K_1\left(a\right)k{V}_I(t-a)T(t-a)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }\left(t\right)-c{V}_I\left(t\right).\end{array}$(14) ) into Equation (23(22) $(\zeta +c)(\zeta +\delta )=\mathcal{R}c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}.$(22) ), we get (25) $(\zeta +c)(\zeta +\delta )(\zeta +d\mathcal{R})=(\zeta +d)c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}.$(25)

Now we will exclude the possibility that (25(24) ${\tilde{V}}_I=(\mathcal{R}-1)\frac{d}{(1-{ϵ}_{\mathrm{RT}})k}.$(24) ) has a complex root λ with a non-negative real part. We prove this by contradiction. Suppose that ζ is a root with a non-negative real part. It is clear that $|(\zeta +c)(\zeta +\delta )|\ge c\delta$ and $|1-f+f{\hat{K}}_1\left(\zeta \right)|\le 1-f+f{\mathcal{K}}_1$. From the assumption $\mathcal{R}>1$, we also have $|\zeta +d\mathcal{R}|>|\zeta +d|$. Thus the modulus of the left-hand side of Equation (25(24) ${\tilde{V}}_I=(\mathcal{R}-1)\frac{d}{(1-{ϵ}_{\mathrm{RT}})k}.$(24) ) is greater than $|\zeta +d|c\delta$, while the modulus of the right-hand side of Equation (25(24) ${\tilde{V}}_I=(\mathcal{R}-1)\frac{d}{(1-{ϵ}_{\mathrm{RT}})k}.$(24) ) is less than or equal to $|\zeta +d|c\delta$. This leads to a contradiction. Therefore, the characteristic equation  (25(24) ${\tilde{V}}_I=(\mathcal{R}-1)\frac{d}{(1-{ϵ}_{\mathrm{RT}})k}.$(24) ) has no roots with non-negative real parts. The theorem is proved.

4. Numerical results

4.1. Parameter values

We fixed most of the parameters based on existing experimental data and our previous modelling studies for model simulations [36–38]. Before HIV infection, the CD4+ T cell level is approximately ${10}^3$ cells per µl [4]. We changed the unit to cells/ml and assumed the initial CD4+ T cell level to be ${10}^6$ cells/ml before infection. The death rate of uninfected CD4+ T cells was assumed to be d=0.01 day⁻¹ [32]. From the pre-infection steady state of CD4+ T cells, we obtained that the generation rate (λ) of CD4+ T cells is ${10}^6\left(0.01\right)={10}^4$ cells ml⁻¹ day⁻¹. The viral infection rate k was chosen to be $2.4\times {10}^{-8}$ ml virion⁻¹ day⁻¹ [40]. The death rate of infected T cells is $\delta =1$ day⁻¹ [27]. The total number of virions produced by one infected cell during its lifespan was assumed to be N=2000 virions cell⁻¹ [40]. The viral clearance rate was assumed to be c=23 day⁻¹ [34]. A very small fraction of new infection leads to latently infected cells. We assumed that the fraction of latency is f=0.0001 [37].

We used an overall drug efficacy to represent the effectiveness of combination therapy consisting of RT and protease inhibitors, that is, $ϵ=1-(1-{ϵ}_{\mathrm{RT}})(1-{ϵ}_{\mathrm{PI}})$ [38]. The overall drug efficacy was assumed to be $ϵ=0.99$ for most simulations but different treatment effectiveness have been used to evaluate the influence of suboptimal or 100% effective therapy on the dynamics of virus and the latent reservoir. Because non-infectious virions can still be detected by the RT polymerase chain reaction test and are counted in viral load measurements, we plotted the total viral load (i.e. $V=V_I+V_{NI}$) in all the simulations.

We assumed that the activation rate of latently infected cells (α) to be a function of the latent infection age, a. Because the population of latently infected cells shifts towards cells specific to increasingly rare antigens, the activation rate is eventually decreasing. Thus, we began the simulation with a simple exponential decay function. The activation rate may be low right after the latent infection, and experiences an increase when cells encounter specific antigens, followed by a decline to a lower level. A lognormal-like function will also be used for simulations later.

An example of the exponential decay function takes the following form (26) $\alpha \left(a\right)=a_{min}+(a_0-a_{min}){\mathrm{e}}^{-wa},$(26) where $a_{min}$ is the minimum activation rate of latently infected cells, $a_0$ is the initial activation rate, and w is the decay constant characterizing how fast the activation rate declines to the minimum value. The values of these parameters remain largely undetermined. We chose $a_0=0.05$ day⁻¹, $a_{min}=0.01$ day⁻¹, and w=0.01 day⁻¹ as default values in our simulation. Sensitivity tests of the solution of the model (Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) )) on these parameters have been performed. We also varied the net generation rate of latently infected cells (γ in Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) )) to explore when the model can generate low-level viral load persistence and a very stable latent reservoir in patients receiving long-term antiretroviral therapy. A summary of parameters and their default values are summarized in Table 1.

Table 1. Parameters and values used in the models and simulations.

Parameter	Description	Value and unit	Reference
λ	Generation rate of uninfected cells	${10}^4$ cell ml^−1 day^−1	See text
d	Death rate of uninfected cells	0.01 day^−1	[32]
k	Infection rate	$2.4\times {10}^{-8}$ ml virion^−1 day^−1	[40]
$\alpha \left(a\right)$	Activation rate function	–	–
$a_0$	Initial activation rate	0.05 day^−1	See text
$a_{min}$	Min activation rate	0.01 day^−1	See text
w	Rate of decay of activation	0.01 day^−1	See text
γ	Net generation rate of latently infected cells	varying	See text
p	Homeostatic proliferation rate of latently infected cells	varying	See text
$L_{max}$	Carrying capacity of latently infected cells	0.3 cell ml^−1	See text
f	Fraction of latency	0.0001	[37]
ε	Overall drug efficacy	0.99	See text
δ	Death rate of infected cells	1 day^−1	[27]
N	Viral burst size	2000 virion cell^−1	[40]
c	Viral clearance rate	23 day^−1	[34]

4.2. Numerical results of the model with linear proliferation of latently infected cells

We provided numerical simulations to illustrate the analytical results of the model with a linear proliferation rate (Equation (1(1) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -\mathrm{d}T-kVT\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =fkVT,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)kVT-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}V\left(t\right)& =N\delta T^{\ast }-cV.\end{array}$(1) ) before treatment and Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) under treatment). Backward Euler and the linearized finite difference method were used to discretize the ordinary and partial differential equation, respectively [14]. The integral over age of latent infection was evaluated using Simpson's rule. The net generation rate of latently infected cells γ was chosen to be a half of the minimum activation rate of latently infected cells, that is, $\gamma =0.5a_{min}$, where $a_{min}=0.01$ day⁻¹. The other parameter values were chosen to be the same as those listed in Table 1. Numerical simulation of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -\mathrm{d}T-kVT\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =fkVT,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)kVT-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}V\left(t\right)& =N\delta T^{\ast }-cV.\end{array}$(1) ) with initial CD4+ T cell count $T\left(0\right)={10}^6$ cells per ml, initial viral load $V\left(0\right)={10}^{-3}$ virions per ml, $T^{\ast }\left(0\right)=0$, and $L_0\left(a\right)=0$ for any a>0 yields an infected steady state in the absence of treatment. The dynamics of the latent reservoir and viral RNA level were showed in Figure 1 upper panel.

Figure 1. Dynamics of the latent reservoir and viral load predicted by model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ). Upper panel: The latent reservoir and viral RNA level reach an infected steady state in the absence of drug treatment. The generation rate of latently infected cells γ was chosen to be a half of the minimum activation rate of latently infected cells, that is, $\gamma =0.5a_{min}$, where $a_{min}=0.01$ day⁻¹. Middle panel: The latent reservoir and virus were predicted to die out under treatment. The overall drug efficacy was chosen to be 0.99 and $\gamma =0.1a_{min}$. Lower panel: The latent reservoir and virus were predicted to persist at a low level under treatment. The overall drug efficacy was chosen to be 0.99 and $\gamma =a_{min}$. The other parameter values were chosen to be the same as those listed in Table 1.

Using the steady state of the model in the absence of treatment as the initial condition of the model under antiretroviral therapy, we performed numerical simulation of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ). The overall drug efficacy was chosen to be 0.99 and the net generation rate of latently infected cells was assumed to be one tenth of the minimum activation rate (i.e. $\gamma =0.1a_{min}$). Both the latent reservoir and virus were predicted to die out in this scenario (Figure 1 middle panel). In the lower panel of Figure 1, we used the same overall drug efficacy but assumed that the net generation rate was equal to the minimum activation rate of latently infected cells. We found that both the latent reservoir and virus persist at a low level under prolonged treatment.

4.3. Viral and latent reservoir persistence depends on the relative magnitude of γ and $a_{min}$

Figure 1 suggested that the persistence of low-level latent reservoir and viral load depends on the value of the parameter γ, the net generation rate of latently infected cells. We further studied the effect of the relative magnitude of γ and $a_{min}$ (the minimum activation rate of latently infected cells) on the dynamics of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ). When γ was chosen to be greater than $a_{min}$ (e.g. $\gamma =1.2a_{min}$ in Figure 2), both the latent reservoir and viral load were predicted to rebound (dashed line). When $\gamma =a_{min}$, the latent reservoir and virus persist at a very stable low level (solid line). When γ is less than $a_{min}$ (e.g. $\gamma =0.8a_{min}$ in Figure 2), the latent reservoir and virus decline slowly (dotted line). The overall drug efficacy was fixed to 0.99 for all the three cases and the other parameter values were the same as those listed in Table 1.

Figure 2. The dynamics of the latent reservoir and virus predicted by model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) with different relative magnitudes of γ and $a_{min}$. When the net generation rate of latently infected cells was chosen to be 1.2 times of the minimum activation rate of latently infected cells, that is, $\gamma =1.2a_{min}$, both the latent reservoir and viral load were predicted to rebound (dashed line). When $\gamma =a_{min}$, the latent reservoir and virus persist (solid line). When $\gamma =0.8a_{min}$, the latent reservoir and virus decline slowly to the infection-free steady state (dotted line). The overall drug efficacy is 0.99 and the other parameter values are the same as those listed in Table 1.

The results in Figure 2 show that the relative magnitude of γ and $a_{min}$ determines whether the latent reservoir and viral load rebound, persist, or decline. This can also be explained from the equations in model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ). In the setting of effective antiretroviral therapy, the virus mainly comes from the activation of latently infected cells. When the net generation rate is greater than the long-term activation rate (i.e. the minimum activation value of the exponentially decay function) of latently infected cells, the size of the latent reservoir increases, which leads to more cell activation and viral production. When the net generation rate is equal to or slightly less than the long-term activation rate, the latent reservoir and the viral load persist at a stable level or decay very slowly.

4.4. Sensitivity tests of model (2) on parameters

In this section, we performed sensitivity tests of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) on a number of parameters. Figure 3 shows the sensitivity of the latent reservoir and viral load on the overall drug efficacy (ε) and the fraction of latency (f). We found that there is almost no difference in the level of latently infected cells and viral load when the treatment is able to suppress the viral load to below the detection limit (upper panel of Figure 3). This suggests that treatment intensification by adding additional antiretroviral drugs (leading to a higher overall drug efficacy) does not help to reduce the size of the latent reservoir or further lower the viral load. Similarly, the dynamics of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) under treatment (lower panel of Figure 3) are not sensitive to f because the influx of new latent infection, denoted by $f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T$, is small under effective treatment (i.e. ${ϵ}_{\mathrm{RT}}$ is close to 1). Therefore, although the fraction of latency plays an important role in determining the latent reservoir size during chronic infection before treatment, it does not affect the dynamics during effective combination therapy.

Figure 3. Sensitivity tests of the latent reservoir and viral load on the overall drug efficacy and the fraction of latency. Upper panel: The dynamics of latently infected cells and viral load predicted by model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) with different overall drug effectiveness, that is, $ϵ=0.9$, 0.95, and 1. Lower panel: The predicted dynamics by model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) with different fractions of infection that leads to latency, that is, f=0.01, 0.001, and 0.0001. The net generation rate (γ) was assumed to be the same as the minimum activation rate ($a_{min}$) of latently infected cells. The other parameter values are the same as those listed in Table 1.

In Figure 4, we examined the sensitivity of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) on the activation rate of latently infected cells. As the initial activation rate $a_0$ increases from 0.04 to 0.06 day⁻¹, both the latent reservoir and the viral load decrease because a higher initial activation rate activates or consumes more latently infected cells (upper panel of Figure 4). However, they still persist at a very low level (because the net generation rate γ was assumed to be the same as the minimum activation rate $a_{min}$ of latently infected cells). As the minimum activation rate $a_{min}$ increases from $a_0/4$ to $a_0/2$, it is faster to reach a balance between the net generation and activation of latently infected cells (note that $a_0$ was fixed and γ was assumed to be equal to $a_{min}$). Thus, both the latent reservoir and the viral load are maintained at higher levels (see lower panel of Figure 4).

Figure 4. Sensitivity tests of the latent reservoir and viral load on the activation rate of latently infected cells. Upper panel: The dynamics of latently infected cells and viral load predicted by model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) with different initial activation rates, that is, $a_0=0.06$, 0.05, and 0.04 day⁻¹. Lower panel: The predicted dynamics by model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) with different minimum activation rates of latently infected cells, that is, $a_{min}=a_0/4$, $a_{min}=a_0/3$, and $a_{min}=a_0/2$, where $a_0=0.05$ day⁻¹. The net generation rate (γ) was assumed to be the same as the minimum activation rate ($a_{min}$) of latently infected cells. The other parameter values are the same as those listed in Table 1.

4.5. Simulation with a lognormal-like activation rate function

In the previous simulations, we used an exponential decay function $\alpha \left(a\right)$ as the activation rate (Equation (26(25) $(\zeta +c)(\zeta +\delta )(\zeta +d\mathcal{R})=(\zeta +d)c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}.$(25) )). Latently infected cells may need some time to encounter specific antigens and thus the activation rate is likely to experience an increase from a low level to a peak, followed by a decline to a lower level. In this section, we chose a lognormal-like function to describe such an activation rate. An example of the lognormal-like function is (27) $\alpha \left(a\right)=\frac{1}{\sqrt{2\pi }\sigma a}{\mathrm{e}}^{-(\mathrm{lna}-\mu {)}^2/2{\sigma }^2}+a_{min},a>0,$(27) where the first term in the right-hand side is the probability density function of a lognormal distribution with the location parameter μ and the scale parameter σ, and the second term $a_{min}$ represents the minimum activation rate of latently infected cells, as used in the exponential decay function.

As the maximum activation rate may have different magnitudes and appear at different time, we fixed the location parameter $\mu =1.5$ but varied the scale parameter ($\sigma =0.5$ or 0.8) in the simulation. In the left column of Figure 5, we assumed the minimum activation rate to be $a_{min}=0.01$ day⁻¹, which is the same as the net generation rate γ of latently infected cells. Similar to the case of $\gamma =a_{min}$ in the previous simulations, the model is able to show the persistence of both the latent reservoir and the low viral load. For a larger scale parameter ($\sigma =0.8$), the activation rate increases faster to the maximum value, which induces a more substantial decrease in the latent reservoir. However, there is almost no difference in the decay of the viral load during treatment (Figure 5(a)). In the right column of Figure 5, we assumed $\sigma =0.5$ but chose a larger minimum activation rate ($a_{min}=0.015$ day⁻¹) than the net generation rate ($\gamma =0.01$ day⁻¹) of latently infected cells. Also similar to the case of $\gamma <a_{min}$ in the previous simulations, both the latent reservoir and the viral load are predicted to die out (Figure 5(b)). This further confirms that the stability of the latent reservoir and low viral load persistence depend on the relative magnitude of the net generation rate and the long-term activation rate of latently infected cells.

Figure 5. The dynamics of the latent reservoir and virus predicted by model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) with a lognormal-like activation rate. The lognormal-like function was given in Equation (27(26) $\alpha \left(a\right)=a_{min}+(a_0-a_{min}){\mathrm{e}}^{-wa},$(26) ). (a) The scale parameter of the lognormal function was set to $\sigma =0.5$ or 0.8. The minimum activation rate was assumed to be the same as the net generation rate of latently infected cells ($a_{min}=\gamma =0.01$ day⁻¹). (b) The scale parameter was fixed at $\sigma =0.5$ and the minimum activation rate was chosen to be larger than the net generation rate of latently infected cells ($a_{min}$=0.015 day⁻¹ and $\gamma =0.01$ day⁻¹). The overall drug efficacy is 0.99 and the other parameter values are the same as those listed in Table 1.

4.6. Numerical results of the model with homeostatic proliferation of latently infected cells

Numerical simulations of the model with linear proliferation of latently infected cells (Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) )) show that the modelling prediction is robust to most parameters but is sensitive to the relative magnitude of the net generation rate and the long-term activation rate of latently infected cells (Figure 2). In this section, we modified the model by including homeostatic proliferation of latently infected cells to reduce the sensitivity. Homeostatic proliferation was shown in a mathematical model [36] and confirmed in an experimental study [7] to be critical in maintaining HIV reservoir size and persistence. As used in ref. [36], a logistic term can describe homeostatic proliferation of latently infected cells. The second equation of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) ) becomes the following. The other equations in the system remain unchanged. (28) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)=pL\left(1-\frac{L}{L_{max}}\right)-\alpha \left(a\right)L.$(28) The parameter p is the maximum proliferation rate of latently infected cells. $L_{max}$ is the carrying capacity density of latently infected cells, beyond which proliferation does not occur. An estimate of 32 cells per ml was used for the maximum density of latently infected cells in ref. [36] assuming that there are ${10}^{12}$ total body lymphocytes and maximally $8\times {10}^6$ latently infected cells per patient under highly active antiretroviral therapy. In this paper, we let $L_{max}=0.3$ cells per ml for the age-specific maximum density of latently infected cells. Simulation will show that the predicted steady state of latently infected cells before treatment agrees with that obtained in ref. [36]. Different values of the maximum proliferation rate p and carrying capacity density $L_{max}$ have also been used to evaluate their effect on the dynamics of the latent reservoir and viral load during therapy. The same activation function (Equation (26(25) $(\zeta +c)(\zeta +\delta )(\zeta +d\mathcal{R})=(\zeta +d)c\delta \frac{1-f+f{\hat{K}}_1\left(\zeta \right)}{1-f+f{\mathcal{K}}_1}.$(25) )) was used in the simulation.

We again assumed that the initial CD4+ T cell count is $T\left(0\right)={10}^6$ cells per ml, initial viral load is $V\left(0\right)={10}^{-3}$ virions per ml, $T^{\ast }\left(0\right)=0$, and $L_0\left(a\right)=0$ for any a>0. The homeostatic proliferation rate of latently infected cells was chosen to be the same as the minimum activation rate of latently infected cells, that is, $p=a_{min}=0.01$ day⁻¹. The other parameter values were chosen to be the same as those listed in Table 1. Numerical simulation of the model with homeostatic proliferation of latently infected cells (Equation (28(27) $\alpha \left(a\right)=\frac{1}{\sqrt{2\pi }\sigma a}{\mathrm{e}}^{-(\mathrm{lna}-\mu {)}^2/2{\sigma }^2}+a_{min},a>0,$(27) )) in the absence of drug treatment generated chronic infection with an infected steady state. The dynamics of the latent reservoir and viral load change were shown in the upper panel of Figure 6. The predicted steady states of latently infected cells and viral load are approximately the same as those shown in the upper panel of Figure 1 using the model with linear proliferation of latently infected cells (Equation (1(1) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -\mathrm{d}T-kVT\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =fkVT,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)kVT-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}V\left(t\right)& =N\delta T^{\ast }-cV.\end{array}$(1) )). Using this steady state as the initial condition of the model under treatment and assuming that the overall drug efficacy is 0.99, the latent reservoir and virus were predicted to persist at a low level under long-term treatment (Figure 6 lower panel).

Figure 6. Dynamics of the latent reservoir and viral load predicted by the model including homeostatic proliferation of latently infected cells (Equation (28(27) $\alpha \left(a\right)=\frac{1}{\sqrt{2\pi }\sigma a}{\mathrm{e}}^{-(\mathrm{lna}-\mu {)}^2/2{\sigma }^2}+a_{min},a>0,$(27) )). Upper panel: The latent reservoir and viral RNA level reach an infected steady state in the absence of drug treatment. The homeostatic proliferation rate of latently infected cells p was chosen to be the same as the minimum activation rate of latently infected cells, that is, $p=a_{min}=0.01$ day⁻¹. The carrying capacity of latently infected cells was assumed to be 0.3 cells ml⁻¹. Lower panel: The latent reservoir and virus were predicted to persist at a low level under treatment. The overall drug efficacy was chosen to be 0.99. The other parameter values were chosen to be the same as those listed in Table 1.

We evaluated the sensitivity of the model including homeostatic proliferation (Equation (28(27) $\alpha \left(a\right)=\frac{1}{\sqrt{2\pi }\sigma a}{\mathrm{e}}^{-(\mathrm{lna}-\mu {)}^2/2{\sigma }^2}+a_{min},a>0,$(27) )) of latently infected cells with respect to several parameters. The model with a linear proliferation rate of latently infected cells is relatively sensitive to the proliferation rate γ (Figure 2). In the upper panel of Figure 7, we showed the predicted dynamics of latently infected cells and viral load using different homeostatic proliferation rates, that is, p=0.03, 0.02, and 0.01 day⁻¹. The minimum activation rate of latently infected cells was fixed at $a_{min}=0.01$ day⁻¹. The carrying capacity was fixed at $L_{max}=0.3$ cells ml⁻¹ and the overall drug efficacy was fixed at $ϵ=0.99$. When p is equal to (or slightly less than) the minimum activation rate $a_{min}$, the latent reservoir and viral load persist (or decay slowly) during treatment. This prediction is similar to that illustrated by the model (Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}T\left(t\right)& =\lambda -dT-(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T\\ \frac{\mathrm{\partial }}{\mathrm{\partial }t}L(a,t)+\frac{\mathrm{\partial }}{\mathrm{\partial }a}L(a,t)& =\gamma L(a,t)-\alpha \left(a\right)L(a,t)\\ L(0,t)& =f(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T,L(a,0)=L_0\left(a\right)\\ \frac{\mathrm{d}}{\mathrm{d}t}{T}^{\ast }\left(t\right)& =(1-f)(1-{ϵ}_{\mathrm{RT}})k{V}_{I}T-\delta T^{\ast }+{\int }_0^{\mathrm{\infty }}\alpha \left(a\right)L(a,t)\mathrm{d}a\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_I\left(t\right)& =(1-{ϵ}_{\mathrm{PI}})N\delta T^{\ast }-c{V}_I\\ \frac{\mathrm{d}}{\mathrm{d}t}{V}_{NI}\left(t\right)& ={ϵ}_{\mathrm{PI}}N\delta T^{\ast }-c{V}_{NI}.\end{array}$(2) )) with linear proliferation of latently infected cells (Figure 1 lower panel). When p is greater than $a_{min}$, there is an initial increase in the latent reservoir size and a minor initial increase in the viral load. As the latent reservoir size increases to above the carrying capacity, proliferation of latently infected cells stops and programmed cell death occurs due to the mechanism of homeostasis (i.e. the logistic term becomes negative). Both the latent reservoir and viral load approach a new steady state. Simulation in the first 700 days was shown in the upper panel of Figure 7 (solid and dashed lines). These simulations suggest that the model with homeostatic proliferation of latently infected cells is robust with respect to the relative magnitude of the proliferation rate and the long-term activation rate of latently infected cells in reproducing the persistence of low-level latent reservoir and plasma viremia.

Figure 7. Sensitivity tests of the model with homeostatic proliferation of latently infected cells (Equation (28(27) $\alpha \left(a\right)=\frac{1}{\sqrt{2\pi }\sigma a}{\mathrm{e}}^{-(\mathrm{lna}-\mu {)}^2/2{\sigma }^2}+a_{min},a>0,$(27) )) on parameters. Upper panel: The dynamics of latently infected cells and viral load predicted by model (28(27) $\alpha \left(a\right)=\frac{1}{\sqrt{2\pi }\sigma a}{\mathrm{e}}^{-(\mathrm{lna}-\mu {)}^2/2{\sigma }^2}+a_{min},a>0,$(27) ) with different homeostatic proliferation rates of latently infected cells, that is, p=0.03, 0.02, and 0.01 day⁻¹. The carrying capacity of latently infected cells was fixed at $L_{max}=0.3$ cells ml⁻¹ and the overall drug efficacy was fixed at $ϵ=0.99$. Middle panel: The predicted dynamics with different carrying capacities of latently infected cells, that is, $L_{max}=1$, 0.3, and 0.1 cells ml⁻¹. The homeostatic proliferation rate of latently infected cells was fixed at p=0.01 day⁻¹. Lower panel: The predicted dynamics by model (28(27) $\alpha \left(a\right)=\frac{1}{\sqrt{2\pi }\sigma a}{\mathrm{e}}^{-(\mathrm{lna}-\mu {)}^2/2{\sigma }^2}+a_{min},a>0,$(27) ) with different drug effectiveness, that is, $ϵ=0.8$, 0.9, and 1. $L_{max}$ was fixed at 0.3 cells ml⁻¹ and p was fixed at 0.01 day⁻¹. The other parameter values are the same as those listed in Table 1.

The sensitivity of the latent reservoir and viral load with respect to the carrying capacity of latently infected cells was shown in the middle panel of Figure 7. $L_{max}$ was chosen to be 1, 0.3, or 0.1 cells ml⁻¹. As $L_{max}$ decreases, the size of the latent reservoir decreases and the viral load undergoes a very minor decrease. We also evaluated the effect of different treatment effectiveness on the dynamics (lower panel of Figure 7). As the overall drug efficacy increases from 0.8 to 1, there is almost no difference in the size of the latent reservoir. There is a very small reduction in the viral load. This suggests that using the model with homeostatic proliferation of latently infected cells, treatment intensification does not further reduce the size of the latent reservoir or lower the viral load. This is consistent with the conclusion drawn from the model with linear proliferation of latently infected cells.

5. Discussion

Highly active antiretroviral therapy can control viral replication in many HIV patients but cannot eradicate the virus. The latent reservoir consisting of a population of latently infected CD4+T cells represents a major obstacle to viral elimination. Because of ongoing viral replication due to suboptimal treatment and some other mechanisms for the intrinsic stability of memory T cells [3], the size of the latent reservoir is very stable or decays extremely slowly. Latently infected cells can be activated by their relevant antigens and release new virions, rendering treatment termination or viral elimination impossible.

A basic model including the dynamics of latently infected cells has been introduced by Perelson et al. to study the biphasic viral load decline observed in patients receiving combination therapy of two RT inhibitors and one protease inhibitor [33]. The model was mathematically analysed [30]. Local and global stability results of steady states were obtained, which was shown to depend on the basic reproductive number. Because the steady states of the latent reservoir and viral load are very sensitive to the effectiveness of treatment, the model cannot be used to describe the persistence of low-level latently infected cells and viremia in patients on long-term combination therapy [37]. A few models have been shown to be able to reduce the sensitivity of viral load with respect to the drug efficacy [6]. One assumes the death rate of productively infected cells to be dependent on the infected cell density. The other model assumes that antiretroviral drugs have different effects in two distinct compartments, one of which is a drug sanctuary site and has a reduced drug efficacy. Some other models based on different biological mechanisms such as asymmetric cell division and proliferation of latently infected cells have also been developed to study the dynamics of the latent reservoir under treatment (see a review in ref. [37]).

The pool of latently infected CD4+ T cells is heterogeneous [47]. Cells specific to common antigens are activated soon and cells specific to rare antigens need more time to encounter with antigens and be activated. Thus, the activation rate depends on the time elapsed since the cell is latently infected, which is the age of latent infection. In this paper, we developed a mathematical model including age-structured latently infected cells. We assumed that the activation rate of latently infected cells depends on the age of latent infection. Using an exponentially decreasing activation rate and a lognormal-like function as examples, we found that the model is able to reproduce the dynamics of latently infected cells and low-level viremia during treatment. The relative magnitude of the net generation rate and the long-term activation rate (i.e. the minimum value of the activation function) of latently infected cells determines whether latently infected cells and viral load can persist at a low level. Only when the net generation rate is less than or equal to the long-term activation rate of latently infected cells can the model describe the persistence of latent reservoir and low viral load. If the net generation rate is slightly less than the long-term activation rate of latently infected cells, then the latent reservoir will be very stable or decay extremely slowly. If the difference between them becomes larger, then the decay of the latent reservoir will become more rapid. This can explain the divergent estimates (from 4.6 months to 44 months) of the half life of the decay of the latent reservoir in different clinical trials [8,15,35,46,54].

Using sensitivity test, we found the model is not sensitive to other parameters including the effectiveness of treatment. Even when treatment is assumed to be 100% effective, the viral load and latently infected cells can still persist. This suggests that treatment intensification adding more drugs to current treatment regimens may not help to reduce the size of the latent reservoir or eradicate the virus. This is consistent with a number of clinical trials in which addition of new drugs such as the integrase inhibitor raltegravir did not reduce persistent low-level viremia in patients with HIV suppression during combination antiretroviral therapy [5,18,22,23,28,53]. It also agrees with our recent results of using a multi-stage model to study the dynamics of 2-LTR (long terminal repeat) circles, a marker of recent viral infection [52].

Homeostatic proliferation of latently infected cells was included in a mathematical model [36] to explain the low viral load persistence and the stability of the latent reservoir, and was later confirmed by an experimental study [7]. We also included homeostatic proliferation of latent cells in the age-structured model and found the model is relatively robust in reproducing the dynamics of viral load and latently infected cells in HIV infected individuals under long-term effective treatment. This may explain why a ‘shock and kill’ approach involving activating the latent reservoir with latency-reversing agents has not achieved success up to now [37]. Novel treatment strategies, such as agents that can block the function of cytokines essential for the proliferation of memory T cells [19] or caspase-1 inhibitor that can inhibit pyroptosis (a programmed cell death which was shown to drive CD4+ T cell depletion [12,13,51]), may have the potential to reduce the size of the latent reservoir, maintain the CD4+ T cell population and eradicate the virus.

Viral blips are observed in many patients receiving prolonged antiretroviral therapy. Using an increased activation rate at some time when latently infected cells encounter with very relevant antigens, the model can generate a viral blip. When treatment stops, the viral load is predicted to increase rapidly, faster than the viral rebound observed in patients with sustained viral suppression who terminated highly active antiretroviral therapy [11]. One reason causing the rapid viral rebound may be that the immune response is not included in this model. Suppressive antiretroviral therapy was shown to be able to partially restore the immune response in HIV patients [17], particularly if therapy was initiated very early during primary infection. When the size of the latent reservoir is small, the immune response can control viral replication for a period of time after cessation of therapy. This provides an opportunity for post-treatment viral control [42]. A mathematical model including latently infected cell activation and immune response was used to explore when post-treatment control of HIV infection can occur [10]. Our age-structured model can be extended to include the immune response and evaluate whether treatment interruptions, such as structured treatment interruptions [26], can achieve optimal control of HIV infection.

In summary, we developed a model including age-structured latently infected cells to study the dynamics of the latent reservoir and the virus in HIV patients on combination therapy. The proliferation rate and the long-term activation rate of latently infected cells are important in governing the dynamics. The model with homeostatic proliferation of latently infected cells is very robust in reproducing the persistence of low-level viral load and latently infected cells. This work provides a modelling framework that can be used to evaluate treatment interruption strategies which may achieve ‘functional cure’ for HIV without administrating continuous and lifetime therapy.

Acknowledgments

We thank two reviewers for their detailed and constructive comments that improved the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work is supported by the NSF [grant DMS-1349939] (to LR).