Optimal Control of an HIV Model with Gene Therapy and Latency Reversing Agents

Abstract

In this paper, we study the dynamics of HIV under gene therapy and latency reversing agents. While previous works modeled either the use of gene therapy or latency reversing agents, we consider the effects of a combination treatment strategy. For constant treatment controls, we establish global stability of the disease-free equilibrium and endemic equilibrium based on the value of ${\mathcal{R}}_0.$ We then consider time-dependent controls and formulate an associated optimal control problem that emphasizes reduction of the latent reservoir. Characterizations for the optimal control profiles are found using Pontryagin’s Maximum Principle. We perform numerical simulations of the optimal control model using the fourth-order Runge–Kutta forward-backward sweep method. We find that a combination treatment of gene therapy with latency reversing agents provides better remission times than gene therapy alone. We conclude with a discussion of our findings and future work.

1. Introduction

Uncovering a functional cure or new preventative methods for the Human Immunodeficiency Virus (HIV) remains a challenge and focus of medical research. Acquiring HIV infection leads to a weakened immune system that leaves the host more susceptible to other infections and diseases. If HIV goes untreated for a prolonged interval of time, the infection can progress to acquired immunodeficiency syndrome (AIDS).

Antiretroviral therapy (ART) remains the primary method of limiting the spread of the virus. Though ART is effective at reducing the viral load to undetectable levels, even short lapses in treatment grant the virus the capability to rapidly rebound within weeks [1]. This can be attributed to the presence of a minute number of latently infected cells that store the genetic information of HIV yet are not actively replicating virions, thus evading typical ART strategies. If an individual stops taking ART for a period of time, natural reactivation of this latent reservoir can quickly lead to an increase in viral load. In the search for a functional cure that would alleviate the need for lifelong ART medication, recent studies suggested the usage of other methods such as gene therapy and latency reversing agents to help overcome the limitations of ART treatments [1].

A motivating example for the usage of gene therapy is the Berlin Patient, an individual who received an allogeneic transplant from a donor in hopes of curing his leukemia. This transplant introduced $CCR5$-deficient bone marrow [1,2]. The CCR5 receptor is the main coreceptor in which HIV attaches itself to the target cell, and thus the introduction of this bone marrow transplant into the patient’s body reduced the number of susceptible target CD4${}^{+}$ T-cells by replacing the CCR5${}^{+}$ cells with the more resilient CCR5${}^{-}$ cells [3]. In the subsequent 10 years following the procedure, a rebound of the virus in the Berlin Patient was not observed [1].

The CRISPR-Cas9 gene therapy program seeks to provide a similar treatment to the bone marrow transplant that the Berlin Patient received by altering the DNA of the susceptible target CD4${}^{+}$ T-cells to negatively express the CCR5 coreceptor [4]. However, current levels of efficacy are insufficient to provide a reasonable level of reduction in viral load, and thus, combinations of strategies will likely be needed to increase the odds of reaching a functional cure [1]. To supplement the efficacy of current gene therapy, as well as to explicitly address the existence of the latent reservoir, we additionally consider latency-reversing agents (LRAs), also generally referred to as “shock and kill” therapy. LRAs activate latently infected cells so they can be detected and attacked by the immune system or by ART [5,6].

Mathematical models provide a useful framework for exploring disease dynamics and associated treatment programs [7]. Over the past several decades, systems of ordinary differential equations were used to better understand the complex interaction of HIV and the immune system, as well as explore potential treatment options in the search for a functional cure [8,9,10,11,12]. Additionally, recent HIV modeling approaches incorporated latent reservoirs to better understand viral rebounds [13,14,15,16,17]. In 2018, Ke et al. explored the dynamics of LRAs through a mathematical model with three compartments: latently infected cells, cells activated by LRA treatment, and refractory cells temporarily resistant to LRA treatment [6]. In 2019, Davenport et al. proposed a mathematical model for the use of gene therapy in HIV prevention [1]. The model describes the dynamics between a susceptible population and refractory (nonsusceptible) population that has received gene therapy, in addition to standard interactions with actively infected cells and free virions. The authors used the model to establish efficacy thresholds for gene therapy to achieve a functional cure for HIV for various ranges of basic reproductive ratios, and noted that combination strategies may serve to reduce the efficacy needed from any individual treatment.

The goal of this paper is to build upon the works of Ke et al. [6] and Davenport et al. [1] by constructing and analyzing a mathematical model that incorporates the combination of gene therapy and LRA treatment. We aim to explore possible synergies between these treatment methods for both constant and time-dependent descriptions of the treatment controls.

2. Model Development

In this paper, we use the following system of differential equations to model the dynamics of an HIV infection under the combined effects of CRISPR-Cas9 gene therapy, latency reversing agents, and antiretroviral therapy. For the purpose of our work, we are interested in the interactions between CCR5${}^{-}$ cells $T_u$, CCR5${}^{+}$ cells $T_s$, latently infected cells L, actively infected cells I, and free virions V.

$$\begin{array}{cc}\hfill & \frac{d{T}_u}{dt}=\lambda f-{\delta }_T{T}_u\hfill \\ \hfill & \frac{d{T}_s}{dt}=\lambda (1-f)-{\delta }_T{T}_s-\beta T_{s}V\hfill \\ \hfill & \frac{dL}{dt}=ϵ\beta T_{s}V-(\alpha +\eta +\delta )L+\gamma I\hfill \\ \hfill & \frac{dI}{dt}=(1-ϵ)\beta T_{s}V+(\alpha +\eta )L-(c+\delta +\gamma )I\hfill \\ \hfill & \frac{dV}{dt}=pI-{\delta }_{V}V\hfill \end{array}$$

(1)

In system (1), f represents the fraction of CD4${}^{+}$ T-cells that become nonsusceptible to HIV as a result of gene therapy. Susceptible $T_s$ cells become latently infected at an interaction rate of $ϵ\beta$ and actively infected at an interaction rate of $(1-ϵ)\beta .$ Latently infected cells are naturally activated at a rate of $\eta .$ Latency reversing agents are activated cells at a rate of $\alpha$ and cells stay activated for a duration of ${\displaystyle \frac{1}{\gamma }}$ days before returning to the latent stage. Actively infected cells are removed through antiretroviral therapy at a rate of c, and produce free virions at a rate of p. The target cell death rate, infected cell death rate, and virus clearance rate are given by ${\delta }_T$, $\delta$, and ${\delta }_V$, respectively. We summarize these parameter roles, along with baseline values gathered from literature, in

Table 1. Parameter Values.

Control	Description	Range	Source
f	Fraction of susceptible target cells that become refractory (nonsusceptible) as a result of gene therapy	0–1	[1]
$\alpha$	Cell activation as a result of latency reversing agents	0–1.8 day${}^{-1}$	[6]
c	Death of cells by antiretroviral treatment	0–1 day${}^{-1}$	[18]
Parameter	Description	Value	Source
$\lambda$	Rate of production of susceptible target cells	${10}^5$ cells day${}^{-1}$	[7]
${\delta }_T$	Death rate of target cells	$0.1$ day${}^{-1}$	[7]
$\beta$	The infectivity of free virus into target cells	$2.0\times {10}^{-7}$ cells${}^{-1}$ day${}^{-1}$	[19]
$ϵ$	Fraction of cells that become latently infected	${10}^{-4}$	[18]
$\eta$	Natural activation rate of latently infected cells	$1.63\times {10}^{-6}$ day${}^{-1}$	[1]
$\delta$	Death rate of infected cells	$0.5$ day${}^{-1}$	[7]
$\gamma$	Rate of transition from actively infected to latently infected	$0.5$ day${}^{-1}$	[6]
p	Production rate of free virus by infected cells	1000 virions cells${}^{-1}$ day${}^{-1}$	[18]
${\delta }_V$	Clearance rate of virions	23 day${}^{-1}$	[18]

.

3. Model Analysis

Assuming system (1) is subject to non-negative initial conditions, we note that the system is invariant in the non-negative orthant. Additionally, because the associated vector field is continuously differentiable, given a set of initial conditions there exists a unique solution to our system by the Picard–Lindelof Theorem [20].

Let N denote the total cell population so that $N\left(t\right)=T_u\left(t\right)+T_s\left(t\right)+L\left(t\right)+I\left(t\right)$ and $\Delta =min\{{\delta }_T,\delta \}.$ Then ${\displaystyle \underset{t\to \infty }{lim\; sup}N\left(t\right)\le \frac{\lambda }{\Delta }.}$ This also implies that ${\displaystyle \underset{t\to \infty }{lim\; sup}V\left(t\right)\le \frac{p\lambda }{\Delta {\delta }_V}.}$ Thus, solution trajectories of system (1) remain bounded over time.

3.1. Basic Reproductive Ratio: ${\mathcal{R}}_0$

We denote the disease-free equilibrium by ${\displaystyle E_0=\left(T_u^0,T_s^0,0,0,0\right)=\left(\frac{\lambda f}{{\delta }_T},\frac{\lambda (1-f)}{{\delta }_T},0,0,0\right)}$ and calculate ${\mathcal{R}}_0$ using the next generation matrix. To derive the next generation matrix, we first calculate the transmission matrix F and the transition matrix V evaluated at $E_0$. These are given by:

$$F=\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& -\frac{\beta \lambda (1-f)}{{\delta }_T}\\ 0& 0& 0& 0& \frac{ϵ\beta \lambda (1-f)}{{\delta }_T}\\ 0& 0& 0& 0& \frac{(1-ϵ)\beta \lambda (1-f)}{{\delta }_T}\\ 0& 0& 0& 0& 0\end{array}\right]$$

and

$$V=\left[\begin{array}{ccccc}{\delta }_T& 0& 0& 0& 0\\ 0& {\delta }_T& 0& 0& 0\\ 0& 0& \alpha +\eta +\delta & -\gamma & 0\\ 0& 0& -(\alpha +\eta )& c+\delta +\gamma & 0\\ 0& 0& 0& -p& {\delta }_V\end{array}\right].$$

${\mathcal{R}}_0$ is the spectral radius of the next generation matrix $F{V}^{-1},$ thus

$${\mathcal{R}}_0=\frac{\beta \lambda p(1-f)(\alpha +\eta +\delta (1-ϵ\left)\right)}{{\delta }_T{\delta }_V((\alpha +\eta +\delta )(c+\delta )+\delta \gamma )}.$$

3.2. Global Stability of the Disease Free Equilibrium

Theorem 1.

When ${\mathcal{R}}_0\le 1$, the disease free equilibrium $E_0$ is globally asymptotically stable.

Proof. 

Consider the Lyapunov function

$$\begin{array}{c}\mathcal{L}=\left(T_u-T_u^0-T_u^{0}ln\left(\frac{T_u}{T_u^0}\right)\right)+\left(T_s-T_s^0-T_s^{0}ln\left(\frac{T_s}{T_s^0}\right)\right)\hfill \\ +\frac{\alpha +\eta }{\alpha +\eta +(1-ϵ)\delta }L+\frac{\alpha +\delta +\eta }{\alpha +\eta +(1-ϵ)\delta }I+\frac{\beta \lambda (1-f)}{{\delta }_T{\delta }_V}V.\hfill \end{array}$$

Then

$$\begin{array}{cc}\hfill \frac{d\mathcal{L}}{dt}& =\left(1-\frac{T_u^0}{T_u}\right)(\lambda f-{\delta }_T{T}_u)+\left(1-\frac{T_s^0}{T_s}\right)(\lambda (1-f)-{\delta }_T{T}_s-\beta T_{s}V)+\hfill \\ & \frac{\alpha +\eta }{\alpha +\eta +(1-ϵ)\delta }(ϵ\beta T_{s}V-(\delta +\eta +\alpha )L+\gamma I)+\frac{\alpha +\delta +\eta }{\alpha +\eta +(1-ϵ)\delta }((1-ϵ)\beta T_{s}V+(\eta +\alpha )L-(\delta +c+\gamma )I)+\hfill \\ & \frac{\beta \lambda (1-f)}{{\delta }_T{\delta }_V}(pI-{\delta }_{V}V)\hfill \\ \hfill & =\lambda f\left(2-\frac{T_u}{T_u^0}-\frac{T_u^0}{T_u}\right)+\lambda (1-f)\left(2-\frac{T_s}{T_s^0}-\frac{T_s^0}{T_s}\right)+\left(\frac{(\alpha +\eta )\gamma }{\alpha +\eta +(1-ϵ)\delta }+\frac{\beta \lambda p(1-f)}{{\delta }_T{\delta }_V}-\frac{(\alpha +\delta +\eta )(\delta +c+\gamma )}{\alpha +\eta +(1-ϵ)\delta }\right)I\hfill \\ \hfill & =\lambda f\left(2-\frac{T_u}{T_u^0}-\frac{T_u^0}{T_u}\right)+\lambda (1-f)\left(2-\frac{T_s}{T_s^0}-\frac{T_s^0}{T_s}\right)+\frac{(\delta +c)\left(\right(\alpha +\delta +\eta )+\delta \gamma )}{\alpha +\eta +(1-ϵ)\delta }({\mathcal{R}}_0-1)I.\hfill \end{array}$$

Since the arithmetical mean is always greater than or equal to the geometric mean, it follows that ${\mathcal{R}}_0\le 1$ implies ${\displaystyle \frac{d\mathcal{L}}{dt}\le 0}$ for all $T_u,T_s,L,I,V\in {\mathbb{R}}_{\ge 0}$. Hence, it is a consequence of LaSalle’s Invariance Theorem [20] that the disease free equilibrium is globally asymptotically stable. □

3.3. Global Stability of the Endemic Equilibrium

Theorem 2.

When${\mathcal{R}}_0>1,$there exists a unique endemic equilibrium$(T_u^{*},T_s^{*},L^{*},I^{*},V^{*})$.

Proof. 

To establish equilibria of system (1), we must solve the following system of equations:

$$\begin{array}{cc}\hfill 0& =\lambda f-{\delta }_T{T}_u^{*}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill 0& =\lambda (1-f)-{\delta }_T{T}_s^{*}-\beta T_s^{*}{V}^{*}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill 0& =ϵ\beta T_s^{*}{V}^{*}-(\alpha +\eta +\delta )L^{*}+\gamma I^{*}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill 0& =(1-ϵ)\beta T_s^{*}{V}^{*}+(\alpha +\eta )L^{*}-(c+\delta +\gamma )I^{*}\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill 0& =p{I}^{*}-{\delta }_V{V}^{*}\hfill \end{array}$$

(6)

We see from Equation (2) that ${\displaystyle T_u^{*}=\frac{\lambda f}{{\delta }_T}.}$ If $I^{*}=0$, we recover our disease free equilibrium ${\displaystyle \left(\frac{\lambda f}{{\delta }_T},\frac{\lambda (1-f)}{{\delta }_T},0,0,0\right).}$ To solve for an endemic equilibrium, suppose $I^{*}\ne 0$ and use Equations (3), (4), and (6) to solve for $T_s^{*},L^{*}$ and $V^{*}$ in terms of $I^{*}$:

$$\begin{array}{cc}\hfill T_s^{*}& =\frac{\lambda (1-f)}{{\delta }_T+\beta V^{*}}\hfill \\ \hfill L^{*}& =\frac{ϵ\beta T_s^{*}{V}^{*}+\gamma I^{*}}{\alpha +\eta +\delta }\hfill \\ \hfill V^{*}& =\frac{p{I}^{*}}{{\delta }_V}.\hfill \end{array}$$

Substituting these expressions into Equation (5), we see that the number of internal equilibria is determined by the number of solutions to $f\left(I^{*}\right)=0,$ where

$$\begin{array}{cc}\hfill f\left(I^{*}\right)& =-(c+\gamma +\delta )+\frac{(1-f)p\beta (1-ϵ)\lambda }{{\delta }_T{\delta }_V+p\beta I^{*}}+\frac{(\alpha +\eta )(\gamma {\delta }_T{\delta }_V+p\beta (\gamma I^{*}+(1-f)ϵ\lambda )}{(\alpha +\delta +\eta )({\delta }_T{\delta }_V+p\beta I^{*})}.\hfill \end{array}$$

Note that

• ${\displaystyle f\left(0\right)=\frac{(c+\delta )(\alpha +\delta +\eta )+\delta \gamma }{\alpha +\delta +\eta }({\mathcal{R}}_0-1);}$
• ${\displaystyle f^{\prime }\left(I^{*}\right)=-\frac{(1-f)p^2{\beta }^2\lambda (\alpha +(1-ϵ)\delta +\eta )}{(\alpha +\delta +\eta ){(p\beta I^{*}+{\delta }_T{\delta }_V)}^2};}$ and
• ${\displaystyle \underset{I^{*}\to \infty }{lim}f\left(I^{*}\right)=-c-\frac{\delta (\alpha +\gamma +\delta +\eta )}{\alpha +\delta +\eta }.}$

Thus, when ${\mathcal{R}}_0>1$, $f\left(I^{*}\right)$ has a unique positive root. It follows that there exists a unique endemic equilibrium when ${\mathcal{R}}_0>1.$ □

Theorem 3.

When ${\mathcal{R}}_0>1$, the endemic equilibrium is globally asymptotically stable.

Proof. 

Notate the endemic equilibrium as $(T_u^{*},T_s^{*},L^{*},I^{*},V^{*})$ and consider the Lyapunov function

$$\begin{array}{cc}\hfill \mathcal{L}& =\left(T_u-T_u^{*}-T_u^{*}ln\left(\frac{T_u}{T_u^{*}}\right)\right)+\left(T_s-T_s^{*}-T_s^{*}ln\left(\frac{T_s}{T_s^{*}}\right)\right)+\left(\frac{\alpha +\eta }{\alpha +\eta +(1-ϵ)\delta }\right)\left(L-L^{*}-L^{*}ln\left(\frac{L}{L^{*}}\right)\right)\hfill \\ & +\left(\frac{\alpha +\eta +\delta }{\alpha +\eta +(1-ϵ)\delta }\right)\left(I-I^{*}-I^{*}ln\left(\frac{I}{I^{*}}\right)\right)+\left(\frac{\beta T_s^{*}}{{\delta }_V}\right)\left(V-V^{*}-V^{*}ln\left(\frac{V}{V^{*}}\right)\right).\hfill \end{array}$$

Then

$$\begin{array}{cc}\hfill \frac{d\mathcal{L}}{dt}& =\left(1-\frac{T_u^{*}}{T_u}\right)(\lambda f-{\delta }_T{T}_u)+\left(1-\frac{T_s^{*}}{T_s}\right)(\lambda (1-f)-{\delta }_T{T}_s-\beta T_{s}V)\hfill \\ & +\left(\frac{\alpha +\eta }{\alpha +\eta +(1-ϵ)\delta }\right)\left(1-\frac{L^{*}}{L}\right)(ϵ\beta T_{s}V-(\delta +\eta +\alpha )L+\gamma I)\hfill \\ & +\left(\frac{\alpha +\delta +\eta }{\alpha +\eta +(1-ϵ)\delta }\right)\left(1-\frac{I^{*}}{I}\right)((1-ϵ)\beta T_{s}V+(\eta +\alpha )L-(\delta +c+\gamma )I)+\left(\frac{\beta T_s^{*}}{{\delta }_V}\right)\left(1-\frac{V^{*}}{V}\right)(pI-{\delta }_{V}V)\hfill \\ \hfill & =\lambda f\left(2-\frac{T_u}{T_u^{*}}-\frac{T_u^{*}}{T_u}\right)+{\delta }_T{T}_s*\left(2-\frac{T_s}{T_s^{*}}-\frac{T_s^{*}}{T_s}\right)+2\beta T_s^{*}{V}^{*}-\beta T_s^{*}{V}^{*}\frac{T_s^{*}}{T_s}-\beta T_s^{*}{V}^{*}\frac{I{V}^{*}}{I^{*}V}\hfill \\ & +\left(\frac{(\alpha +\eta )ϵ\beta T_s^{*}{V}^{*}}{\alpha +\eta +(1-ϵ)\delta }\right)\left(1-\frac{L^{*}{T}_{s}V}{L{T}_s^{*}{V}^{*}}\right)+\left(\frac{(\alpha +\eta )\gamma I^{*}}{\alpha +\eta +(1-ϵ)\delta }\right)\left(1-\frac{I{L}^{*}}{I^{*}L}\right)\hfill \\ & +\left(\frac{(\alpha +\delta +\eta )(1-ϵ)\beta T_s^{*}{V}^{*}}{\alpha +\eta +(1-ϵ)\delta }\right)\left(1-\frac{I^{*}{T}_{s}V}{I{T}_s^{*}{V}^{*}}\right)+\left(\frac{(\alpha +\delta +\eta )(\alpha +\eta )L^{*}}{\alpha +\eta +(1-ϵ)\delta }\right)\left(1-\frac{I^{*}L}{I{L}^{*}}\right).\hfill \end{array}$$

Since

$$\beta T_s^{*}{V}^{*}=\frac{(\alpha +\eta )ϵ\beta T_s^{*}{V}^{*}}{\alpha +\eta +(1-ϵ)\delta }+\frac{(\alpha +\delta +\eta )(1-ϵ)\beta T_s^{*}{V}^{*}}{\alpha +\eta +(1-ϵ)\delta }$$

we have

$$\begin{array}{cc}\hfill \frac{d\mathcal{L}}{dt}& =\lambda f\left(2-\frac{T_u}{T_u^{*}}-\frac{T_u^{*}}{T_u}\right)+{\delta }_T{T}_s*\left(2-\frac{T_s}{T_s^{*}}-\frac{T_s^{*}}{T_s}\right)+\frac{(\alpha +\eta )\gamma I^{*}}{\alpha +\eta +(1-ϵ)\delta }\left(2-\frac{I{L}^{*}}{I^{*}L}-\frac{I^{*}L}{I{L}^{*}}\right)\hfill \\ & +\frac{(\alpha +\eta )ϵ\beta T_s^{*}{V}^{*}}{\alpha +\eta +(1-ϵ)\delta }\left(4-\frac{L^{*}{T}_{s}V}{L{T}_s^{*}{V}^{*}}-\frac{I^{*}L}{I{L}^{*}}-\frac{T_s^{*}}{T_s}-\frac{I{V}^{*}}{I^{*}V}\right)+\frac{(\alpha +\delta +\eta )(1-ϵ)\beta T_s^{*}{V}^{*}}{\alpha +\eta +(1-ϵ)\delta }\left(3-\frac{I^{*}{T}_{s}V}{I{T}_s^{*}{V}^{*}}-\frac{T_s^{*}}{T_s}-\frac{I{V}^{*}}{I^{*}V}\right).\hfill \end{array}$$

Thus, ${\displaystyle d\mathcal{L}/dt<0}$ for all $(T_u,T_s,L,I,V)\in {\mathbb{R}}_{\ge 0}^5,$$(T_u,T_s,L,I,V)\ne (T_u^{*},T_s^{*},L^{*},I^{*},V^{*})$, and so the endemic equilibrium is globally asymptotically stable when ${\mathcal{R}}_0>1.$ □

4. Model Predictions

Using the parameter values found in Table 1, in the absence of any treatment we find that system (1) has a baseline basic reproductive ratio of $R_0=8.6948.$ This aligns with typical HIV basal reproductive ratios as described in [1]. By allowing f to vary, we find that $R_0\left(f\right)<1$ when $f>0.8850$, i.e., gene therapy alone is capable of achieving a functional cure provided that at least 88.50% of susceptible target cells become refractory in response to gene therapy treatment. For a combination treatment strategy, we find that using upper limits for the LRA activation rate of $\alpha =1.8$ and ART removal rate of $c=1$ reduces the effectiveness of gene therapy required to achieve $R_0\left(f\right)<1$ to $f>0.8150.$

However, given current optimistic gene therapy efficacy rates of $f=0.3$ [1], we may consider instead the potential benefits of using gene therapy in combination with LRAs and ART therapy to achieve an elongated remission time. We assume the reactivation rate r of cells from the latently infected reservoir is directly proportional to the steady-state size of the latent reservoir $L^{*}$, i.e., $r=\eta L^{*}$, where $\eta$ is the natural activation rate of a latently infected cell. Using gene therapy alone at current efficacy levels of $f=0.3$, we find that the duration of remission $1/r$ is 10.49 days. However, using an LRA activation rate of $\alpha =1.8$ and ART removal rate of $c=1$ improves the duration of remission to 88.13 days. We summarize the added benefit of LRA+ART shock and kill treatment for other levels of gene therapy efficacy in Figure 1.

These results demonstrate that both gene therapy and shock and kill treatments can contribute to extending remission times and progressing toward a functional cure. To better understand optimal treatment protocols and allow for nonconstant treatment controls, we turn to an optimal control model in the sections to follow.

5. Optimal Control

Applying constant controls for arbitrarily long periods of time may be unsustainable, thus we next aim to find an optimal treatment strategy over a finite time horizon that minimizes the total burden of infection while also minimizing the cost of implementing such a strategy. We introduce into the model time-dependent controls for $f\left(t\right)$, the gene therapy treatment, $\alpha \left(t\right)$, the latent reservoir activation, and $c\left(t\right)$, the ART treatment. This yields the system

$$\begin{array}{cc}\hfill & \frac{d{T}_u}{dt}=\lambda f\left(t\right)-{\delta }_T{T}_u\hfill \\ \hfill & \frac{d{T}_s}{dt}=\lambda (1-f\left(t\right))-{\delta }_T{T}_s-\beta T_{s}V\hfill \\ \hfill & \frac{dL}{dt}=ϵ\beta T_{s}V-(\alpha \left(t\right)+\eta +\delta )L+\gamma I\hfill \\ \hfill & \frac{dI}{dt}=(1-ϵ)\beta T_{s}V+(\alpha \left(t\right)+\eta )L-(c\left(t\right)+\delta +\gamma )I\hfill \\ \hfill & \frac{dV}{dt}=pI-{\delta }_{V}V\hfill \end{array}$$

(7)

The objective functional we choose to minimize is given by

$$J(f,\alpha ,c)={\int }_0^{\tau }{A}_{1}L+A_{2}I+A_{3}V+A_4{f}^2+A_5{\alpha }^2+A_6{c}^{2}dt,$$

where $A_1$, $A_2$, $A_3$, $A_4$, $A_5$, and $A_6$ are positive weights and $\tau$ is the final time. The weights $A_4$, $A_5$, and $A_6$ are associated with quadratic costs of implementing the gene therapy, LRA, and ART treatments. The set of admissible control functions is given by

$$U=\{(f,\alpha ,c)\in L^{\infty }{(0,\tau )}^3|(f\left(t\right),\alpha \left(t\right),c\left(t\right))\in [0,1]\times [0,{\alpha }_{max}]\times [0,c_{max}]\forall t\in [0,\tau ]\}.$$

Thus, we seek an optimal control set $(f^{*},{\alpha }^{*},c^{*})$ such that

$$J(f^{*},{\alpha }^{*},c^{*})=min\left\{J(f,\alpha ,c)\right|(f,\alpha ,c)\in U\}.$$

We proceed by utilizing Pontryagin’s Maximum Principle [21] to exchange the problem of minimizing the objective functional for a problem of finding the pointwise minimum of the following Hamiltonian with respect to f, $\alpha$, and c:

$$\mathcal{H}=A_{1}L+A_{2}I+A_{3}V+A_4{f}^2+A_5{\alpha }^2+A_6{c}^2+{\lambda }_1(\lambda f-{\delta }_T{T}_u)+{\lambda }_2(\lambda (1-f)-{\delta }_T{T}_s-\beta T_{s}V)$$

$$+{\lambda }_3(ϵ\beta T_{s}V-L(\delta +\eta +\alpha ))+{\lambda }_4((1-ϵ)\beta T_{s}V+L(\eta +\alpha )-I(\gamma +\delta +c))+{\lambda }_5(pI-{\delta }_{V}V)$$

where ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$, ${\lambda }_4$, and ${\lambda }_5$ are the adjoint or costate variables associated with the state variables $T_u$, $T_s$, L, $I,$ and V. The following theorem establishes the existence of an optimal control set and characterizations for the optimal control functions.

Theorem 4.

There exists an optimal control set $(f^{*},{\alpha }^{*},c^{*})$ that minimizes $J(f,\alpha ,c)$ over U. Furthermore, there exist adjoint variables ${\lambda }_1,{\lambda }_2,{\lambda }_3$, ${\lambda }_4$, and ${\lambda }_5$ that satisfy the following system of differential equations:

$$\begin{array}{cc}\hfill \frac{d{\lambda }_1}{dt}& {\displaystyle =-\frac{\partial \mathcal{H}}{\partial T_u}={\delta }_T{\lambda }_1}\hfill \\ \hfill \frac{d{\lambda }_2}{dt}& {\displaystyle =-\frac{\partial \mathcal{H}}{\partial T_s}=({\delta }_T+\beta V){\lambda }_2-ϵ\beta V{\lambda }_3-(1-ϵ)\beta V{\lambda }_4}\hfill \\ \hfill \frac{d{\lambda }_3}{dt}& {\displaystyle =-\frac{\partial \mathcal{H}}{\partial L}=-A_1+(\delta +\eta +\alpha ){\lambda }_3-(\eta +\alpha ){\lambda }_4}\hfill \\ \hfill \frac{d{\lambda }_4}{dt}& {\displaystyle =-\frac{\partial \mathcal{H}}{\partial I}=-A_2-\gamma {\lambda }_3+(\delta +\gamma +c){\lambda }_4-p{\lambda }_5}\hfill \\ \hfill \frac{d{\lambda }_5}{dt}& {\displaystyle =-\frac{\partial \mathcal{H}}{\partial V}=-A_3+\beta T_S{\lambda }_2-ϵ\beta T_S{\lambda }_3-(1-ϵ)\beta T_S{\lambda }_4+{\delta }_V{\lambda }_5}\hfill \end{array}$$

(8)

with transversality conditions ${\lambda }_1\left(\tau \right)={\lambda }_2\left(\tau \right)={\lambda }_3\left(\tau \right)={\lambda }_4\left(\tau \right)={\lambda }_5\left(\tau \right)=0$. Also, the optimal control set has the characterization:

$$f^{*}=min\left\{1,max\left(0,\frac{\lambda ({\lambda }_2-{\lambda }_1)}{2A_4}\right)\right\}$$

$${\alpha }^{*}=min\left\{{\alpha }_{max},max\left(0,\frac{L{\lambda }_3+I{\lambda }_4}{2A_6}\right)\right\}$$

$$c^{*}=min\left\{c_{max},max\left(0,\frac{I{\lambda }_4}{2A_5}\right)\right\}$$

Proof. 

Using Corollary 4.1 of [22], the existence of an optimal control set follows from the Lipschitz property of system (1) with respect to the state variables, a priori boundedness of the state solutions established in Section 3, and the convexity of the integrand of J with respect to $(f,\alpha ,c)$.

By applying Pontryagin’s Maximum Principle, the adjoint variables must satisfy system (8) with zero final time conditions. To obtain the characterizations of the optimal controls, we solve ${\displaystyle \frac{\partial \mathcal{H}}{\partial f}=\frac{\partial \mathcal{H}}{\partial \alpha }=\frac{\partial \mathcal{H}}{\partial c}=0}$ on the interior of the control set U. Using the bounds on the controls, we obtain the desired characterization. □

We note that the optimality system consists of the state system (7) with initial time conditions, the adjoint system (8) with final time conditions, and the above characterizations of the optimal controls. Due to the opposite time orientations of the state and adjoint systems, the uniqueness of the optimal control solution is only guaranteed for small final time $\tau$ [23].

6. Numerical Results

We follow the method presented in [24] and utilize an iterative procedure to numerically solve the boundary value problem associated with the above optimality system. Using an initial guess for each of the three controls, we first solve the state system forward in time using a fourth-order Runge–Kutta procedure. Using the resulting state values and final time values, we next solve the adjoint system with transversality conditions backwards in time, again with a fourth-order Runge–Kutta scheme. The controls are then updated using a convex combination of the previous control values and the values resulting from the characterizations of each control. This iterative process then repeats until convergence, which occurs when the maximum relative error between each of the state variables, adjoint variables, and controls over successive iterations is less than a specified value $\delta$.

The optimality system was solved using parameter values found in Table 1, a final time of $t_f=10$ days, and the following initial conditions: $T_u\left(0\right)=0$, $T_S\left(0\right)={10}^6$, $L\left(0\right)=1$, $I\left(0\right)={10}^4$, and $V\left(0\right)={10}^6$. Furthermore, reasonable estimates for the weights $A_1$, $A_2$, $A_3$, $A_4$, $A_5$, and $A_6$ are necessary to ensure convergence is achieved and to obtain meaningful optimal control profiles. Using the extra degree of freedom to select $A_3=1$, we then solve $A_{1}L=A_{2}I=A_{3}V=A_4{\int }_0^{t_f}{f}^{2}dt=A_5{\int }_0^{t_f}{\alpha }^{2}dt=A_6{\int }_0^{t_f}{c}^{2}dt$ to balance the contributions of each weight, where we substitute the median values of each state variable and control: $L=16,810.8$, $I=36,405.3$, $V=1,585,277$, $f=0.5$, $\alpha =0.9$, $c=0.5$. We further adjust $A_1$ by a factor of 10 to prioritize reduction of the latent reservoir. This procedure yields weights of $A_1=943.01$, $A_2=43.54$, $A_3=1$, $A_4=634,111$, $A_5=195,713$, $A_6=634,111$.

The optimal control profiles $f^{*}$, ${\alpha }^{*}$, and $c^{*}$ are shown in Figure 2.

We see in each case the maximum allowable control being applied through approximately day 5–6 before steadily decreasing. This allows for optimal reduction in the total burden of infection in the early days of the treatment regimen before the associated quadratic costs of treatment become too large. We note that including shock and kill treatment allows for the optimal efficacy of gene therapy to decline sooner. The impact of the optimal controls on HIV state variables may be seen in Figure 3. In the absence of any controls, the latent reservoir reaches a peak of 269,051 cells, while using gene therapy alone reduces this maximum burden to 229,443 cells. A combination strategy with both gene therapy and shock and kill therapy is far more effective, reducing the peak burden to 66,486 cells and quickly suppressing the latent reservoir to negligible levels. We observe similar results with the actively infected cell population, beginning with a peak of 639,879 cells with no treatment while achieving a reduced maximum burden of 392,943 cells using an optimal combination treatment strategy. The free virion particles are also reduced from a peak of $2.769\times {10}^7$ virions to $1.699\times {10}^7$ virions under optimal combination treatment. To summarize, the optimal combination of gene therapy and LRA/ART therapy yields an overall decrease in the peak latent, actively infected, and free virion population levels of $75.33\%$, $38.59\%$, and $38.63\%$, respectively. As expected, the selection of weights above prioritized a larger reduction in the latent reservoir, which will be crucial for novel treatment strategies aiming to prolong remission times and progress towards a functional cure for HIV.

7. Conclusions

In this paper, we developed a mathematical model for HIV which includes the latent reservoir and captures the dynamics of combination treatment strategies using gene therapy, latency reversing agents, and antiretroviral therapy. We performed a complete global stability analysis of the model, calculating the basic reproductive ratio and establishing the existence and global stability of both disease free and endemic equilibria. We then computed threshold levels of treatment efficacy needed to achieve a functional cure, or alternatively, prolong remission times when reaching the disease-free equilibrium, which would require prohibitively high levels of treatment effectiveness. During this analysis, we showcased and quantified the added benefits of supporting gene therapy treatments with shock and kill methods designed to attack the latent reservoir.

Given that multiple treatment strategies provided tangible benefits to the long-term progression of HIV infection, we next formulated an optimal control problem to achieve the largest reduction in infected compartments using nonconstant treatment profiles. Our results yielded a substantial reduction in the latent reservoir when using the optimal combination of gene therapy, latency reversing agents, and antiretroviral medication.

The model presented here has some limitations which can inspire future directions of study. We do not capture the diminished efficacy of treatments over time due to acquired drug resistance or virus mutation. The mass action description of the interaction between susceptible cells and free virions could also be refined to consider other limited infection rates. Furthermore, additional data should be collected to establish more reliable parameters and weights for the optimal control system. This will potentially allow for extended time windows in the numerical simulation and improve the discovery of optimal combination treatment regimens.