Analysis of HIV latent infection model with multiple infection stages and different drug classes

Abstract

Latently infected CD$4^{+}$ T cells represent one of the major obstacles to HIV eradication even after receiving prolonged highly active anti-retroviral therapy (HAART). Long-term use of HAART causes the emergence of drug-resistant virus which is then involved in HIV transmission. In this paper, we develop mathematical HIV models with staged disease progression by incorporating entry inhibitor and latently infected cells. We find that entry inhibitor has the same effect as protease inhibitor on the model dynamics and therefore would benefit HIV patients who developed resistance to many of current anti-HIV medications. Numerical simulations illustrate the theoretical results and show that the virus and latently infected cells reach an infected steady state in the absence of treatment and are eliminated under treatment whereas the model including homeostatic proliferation of latently infected cells maintains the virus at low level during suppressive treatment. Therefore, complete cure of HIV needs complete eradication of latent reservoirs.

Keywords:

• HIV infection stages
• stability
• latently infected CD$4^{+}$ T cells
• homeostatic proliferation of latently infected cells
• entry inhibitors

1. Introduction

Human immunodeficiency virus (HIV) attacks the immune system by invading CD$4^{+}$ T cells, which are critical for energizing regulated immune responses to fight infections. As we see in Figure 1, the virus contains an enzyme essential for infectivity called protease. There is also a protein shell surrounding the nucleus of the virus called a capsid, which involves two enzymes necessary for HIV replication, reverse transcriptase and integrase, and two strands of RNA. The protease enzyme and the capsid are covered with the envelope, and this is covered with spikes of glycoproteins gp120 and gp41.

Figure 1. The HIV replication cycle [36].

The virus uses CD$4^{+}$ T-cell components to release new viruses and initiate the next rounds of replication. This replication cycle goes through four main infection stages as shown in Figure 1; viral entry into the cell, reverse transcription from viral RNA to DNA, integration of viral DNA into the host cell DNA, and transcription and translation to release infectious virus from the infected cell. During the viral entry stage, with the help of glycoprotein gp120, the virus first attaches to the primary receptor CD4, followed by the CCR5 or CXCR4 co-receptors on the surface of CD$4^{+}$ T cells. Then the virus merges its envelope with the cell wall by the glycoprotein gp41. This process of fusion allows the HIV capsid to enter the CD$4^{+}$ T cell. Inside the cell, the reverse transcription stage starts, and the virus uses its reverse transcriptase enzyme to reverse the order of the transcription process of generating messenger RNA from nuclear DNA to generate proviral DNA. The viral DNA in the integration stage enters the cell nucleus and binds to the cell DNA using HIV integrase enzyme. Once the viral DNA is integrated into the CD$4^{+}$ T-cell DNA, the provirus begins using the cell machinery to generate the viral RNA. This RNA provides instructions for making new viral proteins in long chains. The new HIV proteins and RNA move to the cell surface to release immature (noninfectious) viruses. The newly formed noninfectious virus releases protease enzyme which break up the long protein chains into smaller chains that assemble to create a fully mature (infectious) virus at the cell wall. The infectious virus buds from the cell wall to initiate the next round of replication.

In HIV replication cycle, anti-HIV drug classes act at different stages of infection progression. They can slow the virus progression, reduce the transmission rate and suppress the viral load to below the detection limit. Entry inhibitors block the processes of virus attachment and entry to the target CD$4^{+}$ T cell. Reverse transcriptase inhibitors stop the action of reverse transcriptase and the creation of proviral DNA. Integrase inhibitors prevent the viral DNA from being incorporated into the host CD$4^{+}$ T-cell DNA. Protease inhibitors block the production of infectious virus particles by preventing noninfectious virus from being fully matured. Even if the antiretroviral therapy can eliminate the virus, one of the latently infected CD$4^{+}$ T cells could be activated and rebound the virus spread. These latently infected CD$4^{+}$ T cells present a major barrier to HIV cure and play an important role in low viral load persistence in patients receiving prolonged antiretroviral therapy [8,23,67]. On the other hand, latently infected CD$4^{+}$ T cells may not be completely affected by current anti-HIV therapy and can be activated to rejoin the replication cycle and release new infectious virus particles. Studying and understanding the properties of these cells are therefore important to eradicate HIV. Many HIV models including latent infection with a single infection stage have been developed to investigate the dynamics of HIV [2,3,28,38,43,68,71–74,79,81,83,87]. Most of the models consist of one compartment of productively infected CD$4^{+}$ T cells. They assumed that the HIV replication cycle goes through a single infection stage to produce infectious viruses. In fact, HIV goes through multiple distinct stages to release new infectious virus particles: viral entry, reverse transcription, integration, and transcription and translation. Many mathematical models have been introduced to study the dynamics of HIV during multiple infection stages, and these models better predict the virus dynamics [24,31–33,39,45,57,85,86,88,92].

Many studies including mathematical models, cell culture experiments and patients data introduced the effect of different drug classes on the dynamics of HIV decay. These studies showed that the viral load decay depends on the inhibited stage in HIV replication cycle. Wang et al. and Sedaghat et al. [76,77,85] studied mathematical models and introduced that when drug classes are $100\mathrm{\%}$ operative, the viral load goes through a single phase of decay while when drug classes are not $100\mathrm{\%}$ effective or act at later stages of the viral replication cycle, the viral load decline can have two phases. Moreover, later in HIV replication cycle an inhibitor acts, the more rapid the viral load decay under a $100\mathrm{\%}$ effective treatment [21,24,49,57,76,77,85], whereas inhibitors acting later in the cycle may not result in a faster viral load decay than reverse transcription inhibitor if the drug efficacy is not a $100\mathrm{\%}$ [21,85].

In fact, viral entry stage is different from other stages in the replication cycle since during which the cell has not been HIV infected yet and to infect the cell, the proteins on HIV surface must bind to the proteins on the surface of CD$4^{+}$ T cell, triggering the fusion of the viral envelope with the cell wall. Entry inhibitors work by attaching themselves to proteins on the surface of CD$4^{+}$ T cells or proteins on the surface of HIV to prevent viral entry into the cells. Some entry inhibitors target the proteins on HIV surface (gp120 and gp41) and some others target the receptor (CD4) and co-receptors (CCR5 and CXCR4) on the cell surface [91]. Thus if entry inhibitors successfully block the first stage in HIV replication cycle, HIV would be unable to enter the cells. Furthermore, entry inhibitors help prevent resistance to other two drug classes (nucleoside and non-nucleoside reverse transcriptase inhibitors [NRTIs and NNRTIs and protease inhibitors PIs)10,12,25,26,47,50,55,65]. Patients who have tried and failed to respond to NRTIs, NNRTIs and/or PIs will likely benefit from entry inhibitor drugs because they are from different class [10,12,25,26,37,41,47,50,55,65]. Entry inhibitors like post-attachment inhibitor ibalizumab, CCR5 inhibitor maraviroc and fusion inhibitor enfuvirtide are FDA approved [17,54,64,78,94]. Therefore, the inclusion of entry inhibitors class into the models is beneficial for HIV patients and more useful in studying the influence of drug classes.

Many mathematical models have been developed to describe the effect of the entry inhibitors on HIV progression. Mathematical models in [11,14] suggest that a combination of fusion inhibitors (enfuvirtide) and protease inhibitors gives a better outcome than single-drug activity in suppressing the viral load and inhibiting viral entry into the host cell. Magombedze et al. [48] developed mathematical models of the immuno-pathogenesis of HIV infection by incorporating three different stages of the HIV replication cycle: viral entry, transcription from viral RNA to DNA, and production of HIV viral particles to study the effect of their inhibitors. They showed that any therapy that includes a combination of entry inhibitor and protease inhibitor is the most potent against HIV replication. Mathematical models in [46,80] demonstrated the importance of perfect patient compliance during HIV treatment and interruption in taking enfuvirtide can be worse than no therapy at all.

Wang et al. [86] developed a multi-stage latent infection model, based on the model presented in [85], with all drug classes except entry inhibitors. They incorporated latently infected cells which established in the stage of integration and a compartment of cells that have un-integrated viral DNA with two copies of the viral long-terminal repeat (2-LTR circles) to study the effect of treatment intensification with raltegravir on both viral load and 2-LTR dynamics in patients with sustained viral suppression. They found that treatment intensification induces a minor decrease in the viral load and a minor increase in 2-LRT. Adding additional raltegravir into treatment regimens did not help improve the treatment outcomes. However, the drug class of entry inhibitors has not been incorporated into the model and investigated.

This paper develops an HIV latent infection model with four drug classes: entry, reverse transcriptase, integrase and protease inhibitors. The model is an extension of the models introduced by Wang et al. [85,86] and Gilmore et al. [24], and it includes both entry inhibitors and latently infected CD$4^{+}$ T cells. We investigate the local stability of the steady states. The basic reproduction number $R_0$ of the model is derived to determine whether the infection will die out or persist. We show that chronic disease will be established if $R_0>1$ and will be eradicated if $R_0<1$. Numerical simulations of the model are performed to explain the influence of drug classes on the dynamics of the model. We find that the impact of entry inhibitor on the model's dynamics is the same as protease inhibitor. We also find that integrase inhibitor is the most effective drug in reducing the viral load and latently infected cells compared with the other inhibitors and uninfected CD$4^{+}$ T cells attain the highest level under the effect of entry inhibitor. The minimum effectiveness of each drug class to eliminate the virus is calculated. Moreover, we show that combination therapy of four drug classes doesn't give better outcomes. We also introduce the influence of the latency fraction on the dynamics of the latently infected CD$4^{+}$ T cells and viral load decay. The model with homeostatic proliferation of latently infected cells is shown to be able to produce the low viral load persistence observed in HIV patients receiving lengthy suppressive treatment. We also evaluate the effect of proliferation rate and carrying capacity of latent infected cells on the virus level.

2. HIV latent infection model with four infection stages and four drug classes

We develop an HIV latent infection model with four infection progression stages and four drug classes based on the models in refs. [24,85,86]. Wang et al. introduced the model without latently infected CD$4^{+}$ T cells in [85] to analyse the effect of different drug classes on the dynamics of HIV decay and in [86] they introduced the model without entry inhibitors including the population of cells that contain 2-LTR DNA circles to evaluate the influence of raltegravir intensification on the viral load dynamics. We incorporate both entry inhibitors and latently infected CD$4^{+}$ T cells into the model. Considering that HIV infection inside the CD$4^{+}$ T cells progresses in four major phases: viral entry, reverse transcription, integration and creation infectious virus, the model can be described by the following system of differential equations: (1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) In the model, $T(t)$ is the population of uninfected CD$4^{+}$ T cells at time t and $I_i$, i = 1, 2, 3 represent the population of infected CD$4^{+}$ T cells that have finished stages 1 up to i in the virus life cycle. $L(t)$ is the concentration of infected cells that have finished the process of integration and gone to the latency state. In the fourth phase of HIV life cycle, protease enzyme causes infected cells to produce infectious virus particles. $V_I$ is the concentration of infectious virus particles (that have not been affected by the enzyme inhibitors), whereas $V_{NI}$ is the concentration of non-infectious virus particles (that have been influenced by protease inhibitors). Parameters ${\eta }_E$, ${\eta }_{RT}$, ${\eta }_I$ and ${\eta }_P$ denote the efficacy of four drug classes: entry inhibitors, reverse transcriptase inhibitors, integrase inhibitors and protease inhibitors, respectively. Table 1 gives a brief definition and the values of the parameters used in model 1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) .

Table 1. Parameter values of model 1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) .

Parameter	Description	Value
s	Production rate of uninfected CD$4^{+}$ T cells	${10}^4$ cell ml${}^{-1}$ day${}^{-1}$ [3,61]
$d_T$	Death rate of uninfected CD$4^{+}$ T cells	0.01 day${}^{-1}$ [61]
β	Infection rate of uninfected CD$4^{+}$ T cells with virus particles	$3.6\times {10}^{-7}$ ml virion${}^{-1}$ day${}^{-1}$ [59]
${\beta }_1$	Transition rate of infected CD$4^{+}$ T cells from stage 1 to stage 2	0.1 day${}^{-1}$ [85]
${\beta }_2$	Transition rate of infected CD$4^{+}$ T cells from stage 2 to stage 3	0.1 day${}^{-1}$ [85]
f	Proportion of latently infected CD$4^{+}$ T cells	0.001 [71]
α	Activation rate of latently infected CD$4^{+}$ T cells	0.01 day${}^{-1}$ [71]
${\delta }_1$	Death rate of infected CD$4^{+}$ T cells that have finished stage 1 of infection	0.1 day${}^{-1}$ [85]
${\delta }_2$	Death rate of infected CD$4^{+}$ T cells that have finished infection stages 1 and 2	0.5 day${}^{-1}$ [85]
${\delta }_3$	Death rate of infected CD$4^{+}$ T cells that have finished stages 1, 2 and 3 of infection	1 day${}^{-1}$ [85]
${\delta }_L$	Death rate of latently infected CD$4^{+}$ T cells	$4\times {10}^{-3}$ day${}^{-1}$ [3,71]
N	The number of noninfectious and infectious virus produced by one infected cell during its life span	2000 virus cell${}^{-1}$ [75]
c	Clearance rate of noninfectious and infectious viruses	23 day${}^{-1}$ [66]

We derive the basic reproductive number $R_0$ using the next-generation method in [19,84], we consider the infection and viral production term in the model and obtain the nonnegative matrix F and the nonsingular matrix V as follows: (2) $\begin{array}{rl}F& ={\left(\begin{array}{cccccc}0& 0& 0& 0& 0& (1-{\eta }_E)\beta {\displaystyle \frac{s}{d_T}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)}_{6\times 6},\end{array}$(2) (3) $\begin{array}{rl}V& ={\left(\begin{array}{cccccc}{\delta }_1+(1-{\eta }_{RT}){\beta }_1& 0& 0& 0& 0& 0\\ -(1-{\eta }_{RT}){\beta }_1& {\delta }_2+(1-{\eta }_I){\beta }_2& 0& 0& 0& 0\\ 0& -f(1-{\eta }_I){\beta }_2& {\delta }_L+\alpha & 0& 0& 0\\ 0& -(1-f)(1-{\eta }_I){\beta }_2& -\alpha & {\delta }_3& 0& 0\\ 0& 0& 0& -{\eta }_{P}N{\delta }_3& c& 0\\ 0& 0& 0& -(1-{\eta }_P)N{\delta }_3& 0& c\end{array}\right)}_{6\times 6}.\end{array}$(3) The basic reproduction number can be calculated by the spectral radius of the next-generation matrix $F{V}^{-1}$: $R_0=\rho (F{V}^{-1})$. (4) $F{V}^{-1}={\left(\begin{array}{cccccc}{\displaystyle \frac{(1-{\eta }_{RT}){\beta }_1}{({\delta }_1+(1-{\eta }_{RT}){\beta }_1)}h}& h& 0& 0& 0& {\displaystyle \frac{s(1-{\eta }_E)\beta }{d_{T}c}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)}_{6\times 6},$(4) where (5) $h=\frac{s(1-{\eta }_P)N}{d_{T}c}\frac{(1-{\eta }_E)\beta (1-{\eta }_I){\beta }_2}{({\delta }_2+(1-{\eta }_I){\beta }_2)}[1-f+\frac{\alpha f}{{\delta }_L+\alpha }]$(5) Therefore, $R_0$ is given as follows: (6) $R_0=\frac{sN\beta {\beta }_1{\beta }_2(1-{\eta }_E)(1-{\eta }_{RT})(1-{\eta }_I)(1-{\eta }_P)(({\delta }_L+\alpha )(1-f)+\alpha f)}{d_{T}c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)({\delta }_L+\alpha )}.$(6) Model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) has two equilibria: disease-free equilibrium $E_0=(T_0,0,0,0,0,0)$ and chronic disease equilibrium $E^{\ast }=(T^{\ast },{I_1}^{\ast },{I_2}^{\ast },L^{\ast },{I_3}^{\ast },{V_{NI}}^{\ast },{V_I}^{\ast })$, where (7) $\begin{array}{rl}{T}_0& =\frac{s}{d_T},\\ T^{\ast }& =\frac{s}{d_T{R}_0},\\ {I_1}^{\ast }& =\frac{s(R_0-1)}{R_0({\delta }_1+(1-{\eta }_{RT}){\beta }_1)},\\ {I_2}^{\ast }& =\frac{s{\beta }_1(1-{\eta }_{RT})(R_0-1)}{R_0({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)},\\ L^{\ast }& =\frac{s{\beta }_1{\beta }_2(1-{\eta }_{RT})(1-{\eta }_I)f(R_0-1)}{R_0({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)({\delta }_L+\alpha )},\\ {I_3}^{\ast }& =\frac{d_{T}c(R_0-1)}{(1-{\eta }_E)(1-{\eta }_P)\beta N{\delta }_3},\\ {V_{NI}}^{\ast }& =\frac{d_T{\eta }_P(R_0-1)}{\beta (1-{\eta }_E)(1-{\eta }_P)},\\ {V_I}^{\ast }& =\frac{d_T(R_0-1)}{\beta (1-{\eta }_E)}.\end{array}$(7)

3. Local stability analysis of model (1)

In this section, we study the local stability of the disease-free equilibrium and chronic disease equilibrium.

Theorem 3.1

The disease-free equilibrium $E_0$ of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) is locally asymptotically stable if $R_0<1$ and unstable when $R_0>1$.

Proof.

Using Lemma 1 and Theorem 2 of the paper by van den Driessche and Watmough [84], we find that the disease-free steady state is locally asymptotically stable when $R_0<1$ and unstable when $R_0>1$.

Theorem 3.2

The chronic disease equilibrium $E^{\ast }$ of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) is locally asymptotically stable when it exists, i.e. $R_0>1$.

Proof.

We start by calculating the Jacobian matrix at the infected steady state, it follows the characteristic equation: (8) $\begin{array}{rl}0& =|\lambda I-J(E^{\ast })|\\ 0& =(\lambda +c{)}^2(\lambda +d_T{R}_0)(\lambda +{\delta }_1+(1-{\eta }_{RT}){\beta }_1)(\lambda +{\delta }_3)\\ & \times (\lambda +{\delta }_L+\alpha )(\lambda +{\delta }_2+(1-{\eta }_I){\beta }_2)\\ & +\frac{s}{d_T{R}_0}\beta {\beta }_1{\beta }_2(1-{\eta }_E)(1-{\eta }_{RT})(1-{\eta }_I)(\alpha f+(1-f)(\lambda +{\delta }_L+\alpha ))\\ & \ast (\lambda +d_T{R}_0-(1-{\eta }_P)d_T(R_0-1)(\lambda +c)N{\delta }_3),\end{array}$(8) where λ is the eigenvalue.

The characteristic equation can be rewritten as (9) $\begin{array}{rl}& (\lambda +{\delta }_L+\alpha )\\ & \times [\frac{(\lambda +c{)}^2(\lambda +d_T{R}_0)(\lambda +{\delta }_1+(1-{\eta }_{RT}){\beta }_1)(\lambda +{\delta }_3)(\lambda +{\delta }_2+(1-{\eta }_I){\beta }_2)}{-\lambda -d_T{R}_0+(1-{\eta }_P)d_T(R_0-1)(\lambda +c)N{\delta }_3}\\ & -g\frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{N(1-{\eta }_P)}]\\ & =(1-g)\left[\frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)({\delta }_L+\alpha )}{N(1-{\eta }_P)}\right]\end{array}$(9) where (10) $g=\frac{(1-f)({\delta }_L+\alpha )}{(1-f)({\delta }_L+\alpha )+\alpha f}$(10) It is obvious that 0<g<1. Suppose that the eigenvalue λ has a non-negative real part, using Cauchy–Schwartz inequality and the condition of $R_0>1$ we obtain (11) $\begin{array}{rl}& |\lambda +{\delta }_L+\alpha ||{\displaystyle \frac{(\lambda +c{)}^2(\lambda +d_T{R}_0)(\lambda +{\delta }_1+(1-{\eta }_{RT}){\beta }_1)(\lambda +{\delta }_3)(\lambda +{\delta }_2+(1-{\eta }_I){\beta }_2)}{-\lambda -d_T{R}_0+(1-{\eta }_P)d_T(R_0-1)(\lambda +c)N{\delta }_3}}\\ & -g{\displaystyle \frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{N(1-{\eta }_P)}|}\\ & >({\delta }_L+\alpha )||{\displaystyle \frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{(1-{\eta }_P)N}|}\\ & -|g{\displaystyle \frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{(1-{\eta }_P)N}||}\\ & =|({\delta }_L+\alpha ){\displaystyle \frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{(1-{\eta }_P)N}(1-g)|.}\end{array}$(11) The modulus of the left-hand side of Equation (11(11) $\begin{array}{rl}& |\lambda +{\delta }_L+\alpha ||{\displaystyle \frac{(\lambda +c{)}^2(\lambda +d_T{R}_0)(\lambda +{\delta }_1+(1-{\eta }_{RT}){\beta }_1)(\lambda +{\delta }_3)(\lambda +{\delta }_2+(1-{\eta }_I){\beta }_2)}{-\lambda -d_T{R}_0+(1-{\eta }_P)d_T(R_0-1)(\lambda +c)N{\delta }_3}}\\ & -g{\displaystyle \frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{N(1-{\eta }_P)}|}\\ & >({\delta }_L+\alpha )||{\displaystyle \frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{(1-{\eta }_P)N}|}\\ & -|g{\displaystyle \frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{(1-{\eta }_P)N}||}\\ & =|({\delta }_L+\alpha ){\displaystyle \frac{c({\delta }_1+(1-{\eta }_{RT}){\beta }_1)({\delta }_2+(1-{\eta }_I){\beta }_2)}{(1-{\eta }_P)N}(1-g)|.}\end{array}$(11) ) is greater than the modulus of the right-hand side. This is a contradiction and means that all eigenvalues have negative real parts. Thus the endemic steady state $E^{\ast }$ is locally asymptotically stable when $R_0>1$.

4. Numerical results

In this section we observe numerical simulations using the ode23s solver in MATLAB. In all simulations we plot the total viral concentration V of infectious and non-infectious virus particles, that is $V=V_I+V_{NI}$, because non-infectious virus can still be detected by the reverse transcriptase polymerase chain reaction test (RT-PCR) and counted in the total viral load.

4.1. Numerical simulations of model (1)

We first present numerical simulations of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) without treatment in Figure (2). We chose the initial value of the chronic disease steady state to be $(T(0),I_1(0),I_2(0),L(0),I_3(0),V(0))=({10}^6,0,0,0,0,{10}^{-3})$. Parameter values used are all chosen from Table (1). We assume that the death rate of infected cells increases as the infection stages progress ${\delta }_1=0.1$, ${\delta }_2=0.5$ and ${\delta }_3=1$. The basic reproduction number with these parameter values is 2.6>1. Numerical simulation in Figure 2 confirms our result obtained in Theorem 3.2 and shows that the solution approaches the chronic disease steady state $(383470,30806,5134.3,36.6998,513.2832,44633)$. Thus in the absence of treatment, the infection is predicted to persist.

Figure 2. Dynamics of model 1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) without treatment. Parameter values are chosen from Table 1.

In Figure 3, numerical simulations of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) are performed under treatment. The four drug classes are assumed to be $99\mathrm{\%}$ effective (i.e. ${\eta }_E={\eta }_{RT}={\eta }_I={\eta }_P=0.99$). The same parameters in Table 1 are used and the endemic steady state $(383470,30806,5134.3,36.6998,513.2832,44633)$ is chosen to be the initial value of the model under treatment. The basic reproduction number with these drug values is $6.1847\times {10}^{-8}<1$. Numerical simulations show that the uninfected cells (T) rebound to the pre-infection level ${10}^6$ and all other model compartments: infected cells $(I_1,I_2,$ and $I_3)$, latently infected cells (L) and viral load (V ) are predicted to decline to 0 under treatment. It agrees with the stability result in Theorem 3.1. Therefore, the chronic disease is predicted to die out after treatment. We compare the influence of each drug class administration on HIV eradication taking in consideration eliminating the other drug classes. We find that if entry/protease inhibitor is at least $62\mathrm{\%}$ effective, it can efficiently block the viral entry into the cell and therefore eradicate the virus, whereas the efficacy of reverse transcriptase inhibitor and integrase inhibitor should be greater than or equal to $77\mathrm{\%}$ and $66\mathrm{\%}$, respectively to eliminate the virus.

Figure 3. Dynamics of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) with treatment. The four drug classes ${\eta }_E$, ${\eta }_{RT}$, ${\eta }_I$ and ${\eta }_P$ are assumed to be $99\mathrm{\%}$ effective. Parameter values are chosen from Table 1.

We also examine the sensitivity of latent reservoir and viral load in model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) on the fraction of latency (f). The parameter values and the initial condition are the same as Figure (2), but we test the model for a variety of f. We find that the level of latently infected cells is very sensitive to f before treatment and Figure 4 shows that, as the fraction latency increases from 0.0001 to 0.01 [2], the latently infected cells level increases because a higher fraction of latency activates more latently infected cells while in Figure (5) there is no difference in the latent infection levels under treatment because $99\mathrm{\%}$ efficacy of integrase inhibitor reduces the inflow of new latently infected cells, that is $f{\beta }_2{I}_2$, by $99\mathrm{\%}$ to a small fraction $0.01f{\beta }_2{I}_2$. Since the proportion of infection that went to latency can be activated and rejoin the HIV replication cycle to release new virus particles, the fraction of latency (f) doesn't affect the dynamics of $I_3$ cells and therefore the viral load V as we see in Figures 4 and 5.

Figure 4. Dynamics of latently infected cells and viral load in model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) before treatment with different fractions of latency, that is, f = 0.01, 0.001, and 0.0001 [2]. Parameter values are chosen from Table 1.

Figure 5. Dynamics of latently infected cells and viremia in model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) with different fractions of latency f = 0.01, 0.001, and 0.0001. Drugs are assumed to be ${\eta }_E=0$, ${\eta }_{RT}=0$, ${\eta }_I=0.99$, and ${\eta }_P=0$. Parameter values are chosen from Table 1.

In Figure 6, we show that the dynamics of the model depend on the inhibited stages in the viral replication cycle (see Table 2). We assume the initial values are the same as those in Figure 3. Protease inhibitor and entry inhibitor act to block two consecutive stages in the HIV replication cycle (last and first stages). Protease inhibitor prevents infectious virus production and entry inhibitor hinders viral entry into the cell. Both works to protect uninfected CD$4^{+}$ T cells from viral entry. Therefore, there is almost no difference in the dynamics of model 1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) whether with protease inhibitor or entry inhibitor treatment (Figure 6). That means entry inhibitor is beneficial for HIV patients who have become resistant to PIS or even NRTIs or NNRTIs because it is from different class of drugs. We also find that the most effective treatment in reducing the viral load and latently infected cells is integrase inhibitor followed by reverse transcriptase inhibitor and then entry inhibitor (which has the same effect of protease inhibitor). For the compartments $I_i$, i = 1, 2, 3, inhibitor acting at a stage closer to its stage in the cycle leads to a faster decay, that is, the lists of inhibitors in order of the efficacy on $I_1$, $I_2$ and $I_3$ are $entry/protease\to integrase\to reverse$ transcriptase, reverse $transcriptase\to entry/protease\to integrase$, and $integrase\to reverse$ $transcriptase\to entry/protease$, respectively. The list of drug classes in order of effectiveness on target CD$4^{+}$ T cells is as same as those on $I_1$ cells and preventing earlier stages doesn't make a better effect than protease inhibitor because it blocks the viral production in the very last stage of HIV replication cycle. Numerical simulations in Figure (6) suggest that treatment intensification of the four drug classes (entry inhibitor, reverse transcriptase inhibitor, integrase inhibitor and protease inhibitor) doesn't have better results.

Figure 6. Dynamics of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ) with different drug classes. Parameter values are chosen from Table 1.

Table 2. The drugs in order of effectiveness on model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) ).

Model compartments	Inhibitors in order of effectiveness as shown in Figure (6)
T	entry/protease $\to$ integrase $\to$ reverse transcriptase
$I_1$	entry/protease $\to$ integrase $\to$ reverse transcriptase
$I_2$	reverse transcriptase $\to$ entry/protease $\to$ integrase
$I_3$	integrase $\to$ reverse transcriptase $\to$ entry/protease
L and V	integrase $\to$ reverse transcriptase $\to$ entry/protease

4.2. Numerical results of the model with homeostatic proliferation of latently infected CD$4^{+}$ T cells

In this section, we only consider integrase inhibitor in the treatment because it is shown in the previous section to be the most effective drug class on latently infected cells and viremia. Using parameter values in Table 1, we find that the virus would be eliminated if the efficacy of integrase inhibitor is $66\mathrm{\%}$ or more. In fact, the virus in clinical observations persists at a low level between 1 and 50 RNA copies/ml during lengthy HAART. Here we extend model 1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) to describe this viral persistence by incorporating homeostatic proliferation of latently infected cells. Homeostatic proliferation is described by the logistic expression $pL(1-\frac{L}{L_{max}})$ where p is the maximum proliferation rate and $L_{max}$ is the carrying capacity of latent cells. Equations of the basic model 1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) stay the same except the latent cells' equation becomes the following (12) $\frac{\mathrm{d}L}{\mathrm{d}t}=pL(1-\frac{L}{L_{max}})+f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L.$(12) In the equation, $f(1-{\eta }_I){\beta }_2{I}_2$ represents latently infected cells generation by the ongoing viral replication cycle while $pL(1-\frac{L}{L_{max}})$ is the proliferation by existing latent cells. Parameters p and $L_{max}$ are chosen to be 0.02 $da{y}^{-1}$ and 1 cell/ml [2], respectively. The other parameter values are the same as those in Table 1. The infected steady state $(377700,30000,5018,502.5,5.1925,0,43700)$ is taken to be the initial value for the model under $99\mathrm{\%}$ effective integrase inhibitor. Numerical simulations for 500 days are shown in Figure (7) where p = 0 represents the basic model and a proliferation rate of p = 0.02 per day is used for the new model. Figure 7(a and b) show that the model including homeostatic proliferation of latent cells can generate the persistence of low-level viremia during suppressive integrase inhibitor and drug intensification leads to a reduction in the level of persistent viremia without eradicating the virus. We evaluate the sensitivity of the new model with respect to the parameters used in the logistic term of homeostatic proliferation. As p or $L_{max}$ increases, the viremia level increases (see Figure 7 b and c). Homeostatic proliferation of latently infected cells was included in many mathematical models [2,60,71,86] and experimental study [13] to describe the virus persistence at low level below the detection limit of clinical inspections.

Figure 7. Effect of parameters ${\eta }_I$, p and $L_{max}$ on the viral decay dynamics using the model with homeostatic proliferation of latently infected cells. ($\mathbf{a}$) Viral load decay with different integrase inhibitor effectiveness, ${\eta }_I=1$, 0.99, 0.9 and 0.8. Homeostatic proliferation rate and carrying capacity of latently infected cells are fixed at 0.02 and 1, respectively. ($\mathbf{b}$) Viremia dynamics with different rates of proliferation, p = 0, 0.01, 0.02, and 0.03 day${}^{-1}$ [2]. The carrying capacity of latently infected cells is fixed at $L_{max}=1$ cell/ml and integrase inhibitor is $99\mathrm{\%}$ effective. ($\mathbf{c}$) Viral load with different carrying capacities of latently infected cells, $L_{max}=0.1$, 0.3, 0.5 and 1 cell/ml [2]. Homeostatic proliferation rate and integrase inhibitor are fixed, p = 0.02 and ${\eta }_I=0.99$. All other parameter values are as same as those listed in Table 1. The HIV detection limit is 50 RNA copies/ml.

5. Conclusion and discussion

HAART employs a combination of NRTIs, NNRTIs and/or PIs. Durable suppression of HIV can improve the immune system and HIV-patient health, prevent HIV-associated mortality and reduce the risk of transmission [4,6,7,16,27,69,70]. However, the virus can be suppressed by HAART rather than eradicated [15,23,93] due to the presence of latent infection [8,23,67]. Furthermore, long-term use of anti-HIV therapies causes the emergence of multi-drug resistant HIV strains [5,9,64,65]. The drug-resistant virus is then involved in HIV transmission and more than $25\mathrm{\%}$ of newly infected individuals carrying HIV that are resistant to at least one drug of HAART [29,89]. Drug class of entry inhibitors can be used to address this issue [10,12,25,26,47,50,55,65]. Moreover, entry inhibitor has mechanisms of action different from those of HAART. They act outside the cell targeting the extracellular viral infection steps before the cell is infected. The first peptide-based HIV entry inhibitor, enfuvirtide, was approved by the U.S. FDA in 2003 to treat HIV patients who have become resistant to HAART. Entry inhibitors are potentially more valuable among HIV medications and the search for new therapies for patients who have failed to respond to HAART is a very active area. There are many studies on developing new therapies to counter HIV entry into the cells [1,10,12,18,22,25,26,34,35,37,40,42,47,50,52,53,55,58,62,63,65,82,90]. In this paper, we develop an HIV model of infection progression stages and their inhibitors by including both entry inhibitor treatment and latently infected cells.

Model 1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}t}& =s-d_{T}T-(1-{\eta }_E)\beta V_{I}T,\\ \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =(1-{\eta }_E)\beta V_{I}T-{\delta }_1{I}_1-(1-{\eta }_{RT}){\beta }_1{I}_1,\\ \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =(1-{\eta }_{RT}){\beta }_1{I}_1-{\delta }_2{I}_2-(1-{\eta }_I){\beta }_2{I}_2,\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =f(1-{\eta }_I){\beta }_2{I}_2-{\delta }_{L}L-\alpha L,\\ \frac{\mathrm{d}{I}_3}{\mathrm{d}t}& =(1-f)(1-{\eta }_I){\beta }_2{I}_2-{\delta }_3{I}_3+\alpha L,\\ \frac{\mathrm{d}{V}_{NI}}{\mathrm{d}t}& ={\eta }_{P}N{\delta }_3{I}_3-c{V}_{NI},\\ \frac{\mathrm{d}{V}_I}{\mathrm{d}t}& =(1-{\eta }_P)N{\delta }_3{I}_3-c{V}_I.\end{array}$(1) without latent infection compartment has been introduced by Wang et al. [85] to study viral load decline dynamics under the effect of different drug classes and various drug efficacy. Infected CD$4^{+}$ T cells that have un-integrated viral DNA with two copies of the viral long-terminal repeat were incorporated into the model to explain the influence of raltegravir intensification on the viral load and 2-LTR in patients on suppressive antiretroviral therapy that doesn't include inhibitors of viral entry [86]. They found that adding additional raltegravir into treatment regimens would not help to improve the treatment outcomes. We incorporate both entry inhibitors and latently infected CD$4^{+}$ T cells into the model. Using sensitivity test, we find that viral load and latent infection are not sensitive to the latency fraction despite that the higher fraction of latency activates more latent cells in the absence of treatment. We compare the influence of the four drug classes on the dynamics of the model. We find that entry inhibitors have the same effect as protease inhibitors on the dynamics of the model (Figure 6) and would benefit HIV patients who have tried and failed many of current HIV medications. We show that the dynamics of the model depend on the drug used and they in order of effectiveness are summarized in Table 2. The class of protease/entry inhibitors shows the best effect on uninfected CD$4^{+}$ T cells and preventing earlier stages doesn't make better outcomes. Integrase inhibitor is shown to be the most effective drug in suppressing viral load and latently infected cells. This prediction agrees with the modelling results obtained in [24,76,77,85]. It is also consistent with clinical data [49,56] and experimental data in cell culture [21]. Numerical simulations also show that targeting more than one stage in the HIV replication cycle by adding new drug classes wouldn't increase the chance of eradicating the virus over integrase inhibitors. This is consistent with several clinical trials [20,30,44,51] and previous modelling results [2, 86].

The basic model analysis predicts that the infection persists in the absence of treatment but if treatment is very effective, the basic reproductive number becomes less than 1, and both latently infected cells and viremia will die out. When each drug class is taken alone in the basic model, we find that entry/protease inhibitor can efficiently prevent HIV from entering the cell and then eradicate the virus if it is at least $62\mathrm{\%}$ effective while the efficacy of reverse transcriptase inhibitor and integrase inhibitor should be greater than or equal to $77\mathrm{\%}$ and $66\mathrm{\%}$, respectively to eliminate the virus. Numerical simulations show that the model with homeostatic proliferation of latently infected cells can maintain the viral load at low level and treatment intensification can only lower the virus stability level. Moreover, increase the value of parameters used in the logistic term $pL(1-\frac{L}{L_{max}})$ leads to a higher level of viral load. This shows that eliminating latent cells in HIV patients is essential to achieve the goal of HIV eradication.

In summary, we develop mathematical models of HIV replication cycle including latency state to study the impact of different drug classes on the infection progression and explain the role of latently infected cells in HIV persistence during prolonged potent treatment. We find that entry inhibitor has the same effect as protease inhibitor and would benefit HIV patients who became drug-resistant due to the long-term use of treatment. Integrase inhibitor is shown to be the most potent drug against the replication of latent infected cells and virus particles. Using the basic model, we show that the virus and latently infected cells reach an infected steady state in the absence of treatment and are eliminated under treatment while the model including homeostatic proliferation of latent infected cells maintains the virus at a low level below 50 RNA copies/ml even under lengthy HAART. Therefore, HIV latency is a major source contributing to the viral persistence in patients receiving prolonged HAART and complete cure of HIV needs complete removal of latent reservoirs.

Disclosure statement

No potential conflict of interest was reported by the author(s).