**3. Reproduction Number**

In this section, we compute the reproduction number of model (1) in the cases of no drug, *R*0, and of drug therapy, *R<sup>d</sup> <sup>c</sup>* , and the local stability of its disease-free equilibrium. The basic reproduction number is defined as the number of secondary CD4<sup>+</sup> T cells infections due to a single infected cell in a completely susceptible population. We start with *R*0. We use the next generation method [24].

The disease-free equilibrium of model (1) is given by:

$$P\_0 \quad = \quad \left(T\_{0\prime}I\_{0\prime}I\_{0\prime}V\_{0\prime}T\_{0\prime}, I\_{R\_0\prime}R^0\right) = \left(\frac{\lambda^a}{\mu^a}, 0, 0, 0, 0, 0, 0\right) \tag{2}$$

Using the notation in [24] on system (1), matrices for the new infection terms, *F*1, and the other terms, *V*1, are given as follows. The chosen variables of the model are *L*, *I* and *V* and the procedure is identical to [24].

$$F\_1 = \begin{pmatrix} 0 & \eta \beta\_2^a T\_0 & \eta \beta\_1^a T\_0 \\ 0 & (1 - \eta) \beta\_2^a T\_0 & (1 - \eta) \beta\_1^a T\_0 \\ 0 & 0 & 0 \end{pmatrix},$$

$$V\_1 = \begin{pmatrix} a\_L^a + \mu\_L^a & 0 & 0 \\ -a\_L^a & a^a & 0 \\ 0 & -k^a & c^a \end{pmatrix}$$

The associative basic reproduction number *R*<sup>0</sup> is written as:

$$R\_0 = \rho(F\_1 V\_1^{-1}) = \frac{T\_0 \left(\rho\_1^a k^a + \rho\_2^a c^a\right) \left[(1-\eta)\mu\_L^a + a\_L^a\right]}{a^a c^a \left(a\_L^a + \mu\_L^a\right)}\tag{3}$$

where *<sup>ρ</sup>* indicates the spectral radius of *<sup>F</sup>*1*V*−<sup>1</sup> <sup>1</sup> . The local stability of *<sup>P</sup>*<sup>0</sup> can be determined using Lemmas 1 and 2.

**Lemma 1.** *[25] The disease-free equilibrium P*<sup>0</sup> *of the system (1) is locally asymptotically stable iff all eigenvalues <sup>λ</sup><sup>i</sup> of the linearization matrix of model (1), satisfy* <sup>|</sup>*arg*(*λi*)<sup>|</sup> <sup>&</sup>gt; *<sup>α</sup> <sup>π</sup>* 2 *.*

**Lemma 2.** *The disease-free equilibrium P*<sup>0</sup> *of the system (1) is unstable if R*<sup>0</sup> > 1*.*

**Proof.** Let *M*<sup>1</sup> be given by:

$$M\_1 = \begin{pmatrix} -\mu^a & 0 & -\beta\_2^a T\_0 & -\beta\_1^a T\_0 \\ 0 & -(a\_L^a + \mu\_L^a) & \eta \beta\_2^a T\_0 & \eta \beta\_1^a T\_0 \\ 0 & a\_L^a & (1 - \eta)\beta\_2^a T\_0 - a^a & (1 - \eta)\beta\_1^a T\_0 \\ 0 & 0 & k^a & -c^a \end{pmatrix}$$

Expanding, det(*λ<sup>p</sup>*I4 <sup>−</sup> *<sup>M</sup>*1) <sup>=</sup> 0, where I4 is the 4 <sup>×</sup> 4 identity matrix, we have the following equation in terms of *λ*:

$$\begin{aligned} \left(\lambda^p + \mu^a\right) \left[\lambda^{3p} + \left(a\_L^a + \mu\_L^a + a^a - (1 - \eta)\beta\_2^a T\_0\right)\lambda^{2p} + \left(c^a(a\_L^a + \mu\_L^a + a^a) + (a\_L^a + \mu\_L^a)\right)a^a \right] \\ - (1 - \eta)\beta\_2^a T\_0(c^a + \mu\_L^a) - a\_L^a \beta\_2^a T\_0 - k^a (1 - \eta)\beta\_1^a T\_0\right) \lambda^p + \left(a\_L^a + \mu\_L^a\right) + c^a a^a \\ - \beta\_2^a T\_0 c^a \left(\mu\_L^a (1 - \eta) + a\_L^a\right) - \beta\_1^a T\_0 k^a \left(1 - \eta\right)\mu\_L^a + a\_L^a\right) \end{aligned} \tag{4}$$

Now, the arguments of the roots of the equation, *λ<sup>p</sup>* + *μ<sup>α</sup>* = 0, are given by:

$$\arg(\lambda\_j) = \frac{\pi}{p} + j\frac{2\pi}{p} > \frac{\pi}{M} > \frac{\pi}{2M}$$

where *j* = 0, 1, ..,(*p* − 1).

Using Descartes' rule of signs, we find that there is exactly one sign change of the equation:

$$\begin{aligned} \lambda^{3p} + \left(a\_L^a + \mu\_L^a + a^a + c^a - (1 - \eta)\beta\_2^a T\_0\right)\lambda^{2p} + \left(c^a(a\_L^a + \mu\_L^a + a^a) + (a\_L^a + \mu\_L^a)a^a \right. \\\\ -(1 - \eta)\beta\_2^a T\_0(c^a + \mu\_L^a) - a\_L^a \beta\_2^a T\_0 - k^a (1 - \eta)\beta\_1^a T\_0\right)\lambda^p + \left(a\_L^a + \mu\_L^a\right) + c^a a^a \left[1 - R\_0\right] = 0 \end{aligned} \tag{5}$$

for *R*<sup>0</sup> > 1. Thus, there is exactly one positive real root of the aforesaid equation for which the argument is less than *<sup>π</sup>* <sup>2</sup>*<sup>M</sup>* . As such, we conclude that if *R*<sup>0</sup> > 1 the disease-free equilibrium *P*<sup>0</sup> of the system (1) is unstable.

Now, we discuss the dynamics of system (1) with drugs. The disease-free equilibrium of model (1) with drugs is given by:

$$P\_{11} = \begin{pmatrix} T\_{1\prime}L\_{1\prime}I\_{1\prime}V\_{1\prime}T\_{R\_{1\prime}}I\_{R\_{1\prime}}R\_{1\prime} \end{pmatrix} = \begin{pmatrix} \frac{\lambda^a}{\mu^b + q^a R^\star}, 0, 0, 0, \frac{q^a T\_{1\prime}R^\star}{\overline{d}^a}, 0, R^\star \end{pmatrix} \tag{6}$$

Using the notation in [24] on system (1), matrices for the new infection terms, *F*2, and the other terms, *V*2, are given by:

$$F\_2 = \begin{pmatrix} 0 & \frac{\eta \beta\_L^x \Lambda^a}{\mu^x + q^x R^a} & \frac{\eta \beta\_L^x \Lambda^a}{\mu^x + q^x R^a} & 0\\ 0 & \frac{(1 - \eta) \beta\_L^x \Lambda^a}{\mu^x + q^x R^a} & \frac{(1 - \eta) \beta\_L^x \Lambda^a}{\mu^x + q^x R^a} & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{pmatrix}$$

$$V\_2 = \begin{pmatrix} a\_L^a + \mu\_L^a & 0 & 0 & 0\\ -a\_L^a & a^a + p^a R^a & 0 & 0\\ 0 & -k^a & c^a & 0\\ 0 & 0 & 0 & a^a \end{pmatrix}$$

In this case, the basic reproduction number *R<sup>d</sup> <sup>c</sup>* is computed to be:

$$R\_c^d = \rho(F\_2 V\_2^{-1}) = \frac{\lambda^a \left(\beta\_1^a k^a + \beta\_2^a c^a\right) \left[ (1 - \eta)\mu\_L^a + \mu\_L^a \right]}{(\mu^a + \eta^a R^\star)(a^a + p^a R^\star) c^a \left(\mu\_L^a + \mu\_L^a\right)} \tag{7}$$

where *<sup>ρ</sup>* indicates the spectral radius of *<sup>F</sup>*2*V*−<sup>1</sup> <sup>2</sup> . The stability of disease-free equilibrium in the case of the drug therapy, *P*1, can be determined using the following lemmas:

**Lemma 3.** *[25] The disease-free equilibrium P*<sup>1</sup> *of the system (1) is locally asymptotically stable iff all eigenvalues <sup>λ</sup><sup>i</sup> of the linearization matrix of model (1), satisfy* <sup>|</sup>*arg*(*λi*)<sup>|</sup> <sup>&</sup>gt; *<sup>α</sup> <sup>π</sup>* 2 *.*

**Lemma 4.** *The disease-free equilibrium P*<sup>1</sup> *of the system (1) is unstable if R<sup>d</sup> <sup>c</sup>* > 1*.*

**Proof.** Let *M*<sup>2</sup> be given by:

$$M\_{2} = \begin{pmatrix} -\mu^{a} - q^{a}R^{\star} & 0 & -\frac{\beta\_{1}^{2}\lambda^{4}}{\mu^{a} + q^{a}R^{\star}} & -\frac{\beta\_{1}^{2}\lambda^{4}}{\mu^{a} + q^{a}R^{\star}} & 0 & 0 & 0\\ 0 & -(a\_{L}^{a} + \mu\_{L}^{a}) & \frac{\beta\_{1}^{2}\lambda^{4}}{\mu^{a} + q^{a}R^{\star}} & \frac{\mu\_{1}^{2}\lambda^{4}}{\mu^{a} + q^{a}R^{\star}} & 0 & 0 & 0\\ 0 & a\_{L}^{a} & \frac{(1-\eta)\beta\_{2}^{2}\lambda^{4}}{\mu^{a} + q^{a}R^{\star}} - a^{a} - p^{a}R^{\star} & \frac{(1-\eta)\beta\_{1}^{2}\lambda^{4}}{\mu^{a} + q^{a}R^{\star}} & 0 & 0 & 0\\ 0 & 0 & k^{a} & -c^{a} & 0 & 0 & 0\\ q^{a}R^{\star} & 0 & 0 & 0 & -d^{a} & 0 & \frac{q^{a}\lambda^{a}}{\mu^{a} + q^{a}R^{\star}}\\ 0 & 0 & p^{a}R^{\star} & 0 & 0 & -d^{a} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & -g^{a} \end{pmatrix}$$

Expanding, det(*λ<sup>p</sup>*I7 <sup>−</sup> *<sup>M</sup>*2) <sup>=</sup> 0, where I7 is the 7 <sup>×</sup> 7 identity matrix, we have the following equation in terms of *λ*:

(*λ<sup>p</sup>* + *μ<sup>α</sup>* + *qαR*-) (*λ<sup>p</sup>* + *dα*) (*λ<sup>p</sup>* + *aα*) (*λ<sup>p</sup>* + *gα*) *λ*3*<sup>p</sup>* + *aα <sup>L</sup>* + *<sup>μ</sup><sup>α</sup> <sup>L</sup>* <sup>+</sup> *<sup>a</sup><sup>α</sup>* <sup>+</sup> *<sup>p</sup>αR*- <sup>+</sup> *<sup>c</sup><sup>α</sup>* <sup>−</sup> (<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)*β<sup>α</sup>* <sup>2</sup> *<sup>λ</sup><sup>α</sup> μα*+*qαR*- *λ*2*p*+ *cα*(*a<sup>α</sup> <sup>L</sup>* + *<sup>μ</sup><sup>α</sup> <sup>L</sup>* <sup>+</sup> *<sup>a</sup><sup>α</sup>* <sup>+</sup> *<sup>p</sup>αR*-)+(*a<sup>α</sup> <sup>L</sup>* + *<sup>μ</sup><sup>α</sup> <sup>L</sup>*)(*a<sup>α</sup>* <sup>+</sup> *<sup>p</sup>αR*-) <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)*β<sup>α</sup>* <sup>2</sup> *<sup>λ</sup><sup>α</sup> μα*+*qαR*- (*c<sup>α</sup>* + *μ<sup>α</sup> <sup>L</sup>*) <sup>−</sup> *<sup>a</sup><sup>α</sup> Lβ<sup>α</sup>* <sup>2</sup> *<sup>λ</sup><sup>α</sup> μα*+*qαR*- <sup>−</sup> *<sup>k</sup>α*(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)*β<sup>α</sup>* <sup>1</sup> *<sup>λ</sup><sup>α</sup> μα*+*qαR*- *λp* +(*a<sup>α</sup> <sup>L</sup>* + *<sup>μ</sup><sup>α</sup> <sup>L</sup>*) + *<sup>c</sup><sup>α</sup>* (*a<sup>α</sup>* <sup>+</sup> *<sup>p</sup>αR*-) <sup>−</sup>*β<sup>α</sup>* <sup>2</sup> *<sup>λ</sup><sup>α</sup> μα*+*qαR cα*(*μ<sup>α</sup> <sup>L</sup>*(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*) + *<sup>a</sup><sup>α</sup> <sup>L</sup>*) <sup>−</sup> *<sup>β</sup><sup>α</sup>* <sup>1</sup> *<sup>λ</sup><sup>α</sup> μα*+*qαR k<sup>α</sup>* (<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)*μ<sup>α</sup> <sup>L</sup>* + *<sup>a</sup><sup>α</sup> L* = 0 (8)

Now, the arguments of the roots of the equation, *λ<sup>p</sup>* + *μ<sup>α</sup>* + *qαR*- = 0, *λ<sup>p</sup>* + *d<sup>α</sup>* = 0, *λ<sup>p</sup>* + *a<sup>α</sup>* = 0, and *λ<sup>p</sup>* + *g<sup>α</sup>* = 0, are given by:

$$\arg(\lambda\_j) = \frac{\pi}{p} + j\frac{2\pi}{p} > \frac{\pi}{M} > \frac{\pi}{2M}$$

where *j* = 0, 1, ..,(*p* − 1).

Using Descartes' rule of signs, we find that there is exactly one sign change of the equation:

$$\begin{aligned} \lambda^{3p} &+ \left(a\_L^a + \mu\_L^a + a^a + p^a R^\star + c^a - (1 - \eta)\beta\_2^a \frac{\lambda^a}{\mu^a + q^a R^\star}\right) \lambda^{2p} \\ &+ \left(c^a \left(a\_L^a + \mu\_L^a + a^a + p^a R^\star\right) + \left(a\_L^a + \mu\_L^a\right) \left(a^a + p^a R^\star\right)\right) \\ &- (1 - \eta)\beta\_2^a \frac{\lambda^a}{\mu^a + q^a R^\star} \left(c^a + \mu\_L^a\right) - a\_L^a \beta\_2^a \frac{\lambda^a}{\mu^a + q^a R^\star} - k^a (1 - \eta) \beta\_1^a \frac{\lambda^a}{\mu^a + q^a R^\star}\right) \lambda^p \\ &+ \left(a\_L^a + \mu\_L^a\right) + c^a \left(a^a + p^a R^\star\right) \left[1 - R\_c^d\right] = 0 \end{aligned} \tag{9}$$

for *R<sup>d</sup> <sup>c</sup>* > 1. Thus, there is exactly one positive real root of the aforesaid equation for which the argument is less than *<sup>π</sup>* <sup>2</sup>*<sup>M</sup>* . As such, we conclude that, if *<sup>R</sup><sup>d</sup> <sup>c</sup>* > 1, the disease-free equilibrium *P*<sup>1</sup> of the system (1) is unstable.

#### **4. Global Stability of the Disease-Free Equilibrium**

In this section, we compute the global stability of the disease-free equilibrium *P*<sup>1</sup> of the model (1). Following Castillo & Chavéz [26], we rewrite model (1) as:

$$\begin{aligned} \frac{d^\mu X}{dt^\mu} &= F(X, Z) \\\\ \frac{d^\mu Z}{dt^\mu} &= G(X, Z), \qquad G(X, 0) = 0 \end{aligned} \tag{10}$$

where *<sup>X</sup>* = (*T*, *TR*, *<sup>R</sup>*) and *<sup>Z</sup>* = (*L*, *<sup>I</sup>*, *<sup>V</sup>*, *IR*), with *<sup>X</sup>* <sup>∈</sup> **<sup>R</sup>**<sup>3</sup> <sup>+</sup> being the number of uninfected and non-susceptible CD4<sup>+</sup> T cells and drugs, and *<sup>Z</sup>* <sup>∈</sup> **<sup>R</sup>**<sup>4</sup> <sup>+</sup> denoting the number of latent and infected CD4<sup>+</sup> T cells, virus, and non-infectious CD4<sup>+</sup> T cells.

The disease-free equilibrium is written as *U* = (*X*-, 0), where *X*- = *T*1, *TR*<sup>1</sup> , *R*<sup>1</sup> = *λα μα*+*qαR*- , *<sup>q</sup>αT*1*R*- *<sup>d</sup><sup>α</sup>* , *<sup>R</sup>*- .

The conditions (*H*1) and (*H*2) must be met to guarantee the global asymptotic stability of the disease-free equilibrium of the model (1):

$$\begin{aligned} \text{( $H\_1$ )}: \quad \text{For } \frac{d^\bullet X}{dt^\bullet} = F(X, 0), \quad X^\star \text{ is globally asymptotically stable} \\\\ \text{( $H\_2$ )}: \quad \text{G}(X, Z) = AZ - \hat{\text{G}}(X, Z), \text{ } \hat{\text{G}} \ge 0, \text{ for } (X, Z) \in \Upsilon\_1 \end{aligned} \tag{11}$$

where *A* = *DZG*(*X*-, 0) can be written in the form *A* = *M* − *D*, where *M* ≥ 0 (*mij* ≥ 0) and *D* > 0 is a diagonal matrix. Υ<sup>1</sup> is the region where the model makes biological sense. If the system (10) satisfies the conditions in (11) the following theorem holds [26].

**Theorem 1.** *The fixed point U* = (*X*-, 0) *is a globally asymptotically stable equilibrium of the system (10) provided that R<sup>d</sup> <sup>c</sup>* < 1 *and that the assumptions in (11) are satisfied.*

**Proof.** Let

$$F(X,0) = \begin{bmatrix} \lambda^a - \mu^a T - q^a T R \\ q^a T R - d^a T\_R \\ R\_k^a - g^a R \end{bmatrix} \tag{12}$$

and

$$A = \begin{pmatrix} -(a\_L^a + \mu\_L^a) & \eta \beta\_2^a T\_1 & \eta \beta\_1^a T\_1 & 0\\ a\_L^a & (1 - \eta) \beta\_2^a T\_1 - (a^a + p^a R^\star) & (1 - \eta) \beta\_1^a T\_1 & 0\\ 0 & k^a & -c^a & 0\\ 0 & p^a R^\star & 0 & -a^a \end{pmatrix} \tag{13}$$

and

$$\mathcal{G}(X,Z) = \begin{pmatrix} \hat{\mathcal{G}}\_1(X,Z) \\ \hat{\mathcal{G}}\_2(X,Z) \\ \hat{\mathcal{G}}\_3(X,Z) \\ \hat{\mathcal{G}}\_4(X,Z) \end{pmatrix} = \begin{pmatrix} \eta \beta\_1^a V T\_1 \left(1 - \frac{T}{T\_1}\right) + \eta \beta\_2^a I T\_1 \left(1 - \frac{T}{T\_1}\right) \\ (1 - \eta) \beta\_1^a V T\_1 \left(1 - \frac{T}{T\_1}\right) + (1 - \eta) \beta\_2^a I T\_1 \left(1 - \frac{T}{T\_1}\right) \\ 0 \\ 0 \end{pmatrix} \tag{14}$$

All conditions are satisfied, thus *U*<sup>0</sup> is globally asymptotically stable.

## **5. Numerical Results**

We simulate the model (1) for different values of the order of the fractional derivative, *α* and for clinically relevant parameters. We have applied the Predictor–Evaluator–Corrector–Evaluator PECE method of Adams–Bashford–Moulton type [27]. The parameters used in the simulations, based on [8,28], are: *λ* = 100 μL−<sup>1</sup> day−*α*, *μ* = 0.1 day−*α*, *a* = 0.3 day−*α*, *c* = 3 day−*α*, *k* = 210 day−*α*, *<sup>β</sup>*<sup>1</sup> <sup>=</sup> 1.5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> day−*α*, *<sup>β</sup>*<sup>2</sup> <sup>=</sup> 1.5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> day−*α*, *<sup>p</sup>* <sup>=</sup> 0.1 <sup>μ</sup>M−<sup>1</sup> day−*α*, *<sup>q</sup>* <sup>=</sup> 0.1 <sup>μ</sup>M−<sup>1</sup> day−*α*, *<sup>g</sup>* <sup>=</sup> 2.7726 day−*α*, *<sup>η</sup>* <sup>=</sup> 0.03, *aL* <sup>=</sup> 0.1 day−*α*, *<sup>μ</sup><sup>L</sup>* <sup>=</sup> <sup>4</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> day−*α*, *Rk* <sup>=</sup> 2.5, *<sup>τ</sup>* <sup>=</sup> 0.5 day*α*, and the initial conditions are: *T*(0) = 1000, *L*(0) = *I*(0) = *TR*(0) = *IR*(0) = 0, *V*(0) = 50 and *R*(0) = 2.5.

In Figures 1 and 2, we consider model (1) without PEP, for different values of the order of the fractional derivative, *α*. The concentration of CD4<sup>+</sup> T cells decreases over time and with *α*. This suggests that the infection is more severe as *α* is lowered. This pattern is supported by the graphs in Figure 2, where it is observed a ratio of healthy T cells to total T cells starting with 0.5 for *α* = 1.0, and decreasing for *α* = 0.9 and *α* = 0.7. Moreover, this ratio points to chronic infection, as the disease evolves, for all *α*.

**Figure 1.** Dynamics of the CD4<sup>+</sup> T cells, *T*, of system (1) without PEP for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text.

**Figure 2.** Ratio of healthy CD4<sup>+</sup> T cells, *T*, to total CD4<sup>+</sup> T cells, *T* + *L* + *I*, of system (1) without PEP for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text.

In Figures 3 and 4, we plot the dynamics of the drug *R* and of the basic reproduction number *R<sup>d</sup> c* , for different values of the order of the fractional derivative, *α*. These figures show that the dosage of the drug is important for controlling HIV infection, since *R<sup>d</sup> <sup>c</sup>* varies with *R*. As *R* increases, smaller values of *R<sup>d</sup> <sup>c</sup>* are observed, which indicate less infection. Moreover, the value of the fractional derivative, *α*, may also contribute to controlling the severity of the infection, since smaller values of *R<sup>d</sup> <sup>c</sup>* are observed with decreasing *α*.

**Figure 3.** Drug concentration in the plasma, *R*, given by system (1) for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text.

**Figure 4.** Basic reproduction number *R<sup>d</sup> <sup>c</sup>* of system (1) for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text.

Figure 5 depicts the HIV viral load for a dosage *Rk* = 5 and dosing interval *τ* = 0.5 day, for different values of the order of the fractional derivative, *α*. As it is shown, the dosage of the drug and the dosing interval are sufficient to control the infection, with the viral load going asymptotically to zero. Similar patterns are seen for all values of *α*, with higher initial viral load for smaller *α*, but faster velocity of convergence.

**Figure 5.** HIV viral-load, *V*, of the system (1) for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text, except *Rk* = 5.

Figure 6 shows the ratio of the infected CD4<sup>+</sup> T cells to total CD4<sup>+</sup> T cells in the presence and absence of PEP, with low drug dosage, for different values of the order of the fractional derivative, *α*. The ratio of infected to total CD4<sup>+</sup> T cells is always smaller when patients are under PEP, when compared to the case without treatment.

**Figure 6.** Ratio of infected CD4<sup>+</sup> T cells, *I* + *IR*, to total CD4<sup>+</sup> T cells, *T* + *L* + *I* + *TR* + *IR* of system (1) for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text, except *Rk* = 1.

In the next figures, we study the effect of different treatment strategies in the dynamics of HIV. We start in Figure 7 with two different treatment strategies: drug perfect adherence and drug therapy breaks. The last strategy consists of intervals (days) in which the therapy is stopped (Δ*Rk* = 0) followed by intervals where there is perfect drug adherence. Perfect adherence therapy consists of taking a dosage Δ*Rk* = *Rk* for all *t* = *tk*. In Figure 7, the drug therapy breaks consist of stopping drug application for two days, followed by drug perfect adherence strategy for five days. It is observed that the elimination of HIV from the body takes longer for drug therapy breaks. This is seen for all *α*. Moreover, despite a higher initial peak, the asymptotic HIV viral load is reached faster for smaller *α*.

**Figure 7.** HIV viral load of system (1) for two different therapy strategies and *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions are in the text, except *Rk* = 4.5. For more information, see text.

Figure 8 shows another example of the dynamics of the HIV, this time for three distinct treatment strategies: without treatment, drug perfect adherence, and drug therapy breaks. The intervals for the drug therapy breaks are in this case as follows. Five days of no drug administration, which are followed by five more days of perfect drug adherence strategy. The model provides oscillating solutions for the case of drug therapy breaks, as is seen in the figure.

In Figure 9 we show the dynamics of HIV for increasing values of the cell to cell transmission rate, *β*2, for three treatment strategies: no treatment, drug perfect adherence and drug therapy breaks, and for varying *α*. The drug therapy breaks strategy stops drug application for 15 days, followed by another 15 days of perfect drug adherence strategy. We observe higher peaks of the viral load and the corresponding curve, in the case of drug therapy breaks, is between the curves of no treatment and drug perfect adherence. This behaviour is repeated for all *α*.

The simulations of the model reveal that a combination of sufficient drug dosage and drug frequency may induce better efficacy of PEP. Drug perfect adherence strategy is always better than the other two. Nevertheless, one must think about the side effects of ART, though their toxicity has been reduced as medicine evolves and new treatment options appear.

We now proceed with the simulation of the effect of the latent reservoir in the dynamics of HIV infection, under the conditions of Figure 8. We consider three treatment strategies: without treatment, drug perfect adherence, and drug therapy breaks. The intervals for the drug therapy breaks consist of 10 days. In the first five days, the drug is halted, whereas for the last five days, the perfect drug

adherence strategy is applied. The difference from Figure 8 is in the value of *η*, which represents the proportion of latently infected CD4<sup>+</sup> T cells. The value of *η* is reduced from 0.03 to 0.01. Figure 10 shows slight higher peaks of HIV for *η* = 0.01, in particular, for smaller values of *α*.

**Figure 8.** HIV, viral load, *V*, of system (1) for three different therapy strategies and for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text, except *Rk* = 1. For more information, see text.

**Figure 9.** HIV, viral load, *V*, of system (1) for three different therapy strategies and for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text, except *Rk* = 1 and *β*<sup>2</sup> = 0.0015.

In Figure 11, we plot the viral load for two values of *η*, the fraction of latent infected cells. We note that the asymptotic value of the virus is the same for all *α*. Nevertheless, there are subtle changes in the dynamics of the virus. In the transient are observed smaller values of HIV viral load for *η* = 0.03, whereas in the asymptotic value there is a switch in this behaviour, and higher values of HIV are seen for *η* = 0.03. This may be explained as follows. When *η* = 0.03 > *η* = 0.01, there are more latently infected cells in the body. If these cells encounter an antigen or are exposed to specific cytokines or chemokines, they become actively infected by proviral transcription. The latter causes viral rebound if a patient stops ART. This happens earlier for smaller values of *α*.

**Figure 10.** HIV, viral load, *V*, of system (1) for different therapy strategies and for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text, except *Rk* = 1 and *η* = 0.01. For more information, see text.

**Figure 11.** HIV, viral load, *V*, of system (1) for two values of *η*, the proportion of latently infected cells, and for *α* = 1 (**top left**), *α* = 0.9 (**top right**) and *α* = 0.7 (**bottom**). Parameter values and initial conditions in the text, except for *Rk* = 2.5. For more information, see text.
