\documentclass[reqno]{amsart} \usepackage{graphicx} \usepackage{hyperref} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2005(2005), No. 52, pp. 1--22.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2005 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2005/52\hfil Optimal control of combined therapy] {Optimal control of combined therapy in a single strain HIV-1 model} \author[ W. Garira, S. D. Musekwa, T. Shiri\hfil EJDE-2005/52\hfilneg] {Winston Garira, Senelani D. Musekwa, Tinevimbo Shiri} % in alphabetical order \address{Winston Garira \hfill\break Department of Applied Mathematics, National University of Science and Technology, P O Box AC939 Ascot, Bulawayo, Zimbabwe} \email{wgarira@nust.ac.zw} \address{Senelani D. Musekwa \hfill\break Department of Applied Mathematics, National University of Science and Technology, P O Box AC939 Ascot, Bulawayo, Zimbabwe} \email{sdmusekwa@nust.ac.zw} \address{Tinevimbo Shiri \hfill\break Department of Applied Mathematics, National University of Science and Technology, P O Box AC939 Ascot, Bulawayo, Zimbabwe} \email{tshiri@nust.ac.zw} \date{} \thanks{Submitted April 16, 2005. Published May 17, 2005.} \subjclass[2000]{93C15, 92D30, 49K15, 49J15, 34B15} \keywords{HIV; mathematical model; optimal control, combined therapy} \begin{abstract} Highly active antiretroviral therapy (HAART) is administered to symptomatic human immunodeficiency virus (HIV) infected individuals to improve their health. Various administration schemes are used to improve patients' lives and at the same time suppressing development of drug resistance, reduce evolution of new viral strains, minimize serious side effects, improve patient adherence and also reduce the costs of drugs. We deduce an optimal drug administration scheme useful in improving patients' health especially in poor resourced settings. In this paper we use the Pontryagin's Maximum Principle to derive optimal drug dosages based on a mathematical dynamical model. We use methods of optimal control to determine optimal controls analytically, and then use the Runge-Kutta scheme of order four to numerically simulate different therapy effects. We simulate the different effects of a drug regimen composed of a protease inhibitor and a nucleoside reverse transcriptase inhibitor. Our results indicate that for highly toxic drugs, small dosage sizes and allowing drug holidays make a profound impact in both improving the quality of life and reducing economic costs of therapy. The results show that for drugs with less toxicity, continuous therapy is beneficial. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \allowdisplaybreaks \section{Introduction} Recently, there has been a rollout of antiretroviral (ARV) therapies in many countries around the world, but availability of ARVs in poor resourced settings is a major concern. The cost of these drugs is beyond reach of many infected patients, hence there is need to come up with a comprehensive drug administration scheme that makes a significant impact in conferring clinical benefits and cost effectiveness. Clinical benefits of drug therapy for HIV infected individuals include restoration of CD4+ T cells levels, suppressing viral levels below detection limits and minimizing detrimental side effects such as risk of cardiovascular, acute retroviral syndrome, fat loss, lactic acidosis, abnormal fat distribution and mitochondrial damage \cite{b1}. There are more than twenty anti-HIV-1 drugs available and these are administered in many different combinations of three or four drugs. The drugs fall into three main categories, that is, reverse transcriptase inhibitors (RTIs) (nucleoside, nucleotide and nonnucleoside), protease inhibitors (PIs) and fusion inhibitors (FIs). RTIs prevent new HIV-1 infections by disrupting the conversion of viral RNA into DNA that can be incorporated into the host cell's genome. PIs function by preventing the assembly of key viral proteins after they have been mistakenly produced by infected host cells \cite{k7}. FIs function by preventing the fusion of the virus and the host cells. HAART consists of combined drug regimens that includes two or three nucleoside agents alone or two nucleoside agents combined with a protease inhibitor or a nonnucleoside reverse trancriptase inhibitor \cite{b1}. Examples of such regimen combinations include EFV (Efavirenz) + (3TC (Lamivudine) or FTC (Emtricitabine))+ (AZT (Zidovudine) or TDF (Tenofovir Disoproxil Fumarate)), a combination of a nonnucleoside reverse transcriptase inhibitor (EFV) and two nucleoside reverse transcriptase inhibitors (3TC or FTC and AZT) and LPV/r (Lopinavir) + (3TC or FTC) + AZT, a combination of a protease inhibitor (LPV/r) and two nucleoside reverse trancriptase inhibitors (3TC or FTC and AZT) and other options that are selected by government agencies, although these options are limited by generic formulations \cite{g1}. In this paper we explore the effects of a combination of a protease inhibitor and a nucleoside reverse transcriptase inhibitor, that is, we only look at effects of two types of drugs that are used in a HAART regimen. Suppression of viremia to less than detection limits or maintenance of even partial viremic suppression by selection of an optimal regimen remains the goal of therapy. The ultimate goal is to prevent further immune deterioration. The new chemotherapies offer added dosing convenience and improved safety profiles. Various chemotherapies for patients with HIV-1 are being examined to determine the optimal scheme for treatment \cite{f1}. The primary attention of this paper is to establish when and how treatment should be initiated, dosage size and means to continue clinical benefit in the face of challenges like antiretroviral drug failure and antiretroviral resistances. The optimal controls in this paper represent percentage effects chemotherapies have on the interaction of the CD4+ T cells with the virus (infection of CD4+ T cells) and the virions produced by infected cells (burst size). Chemotherapy has side effects if administered in high dosage sizes or continuously, therefore the length of treatment is a limited time frame. The interval of treatment is necessary since a plausible assumption is made that chemotherapy only has a certain designated time for allowable treatment \cite{k2}; \cite{f1}. After some finite time frame, HIV-1 is able to build up resistance to the treatment due to its mutation ability. Therefore, in this paper we fix the length of treatment. In this paper we need to determine optimal methodology for administering anti-viral medication therapies to fight HIV-1 infection. The main reasons for such an optimal therapy are minimization of drug toxicity or systemic cost, maximization of CD4+ T cell count and minimize cost of drugs. Optimal control methods have been applied to the derivation of optimal therapies for HIV infection. Butler {\it et al.} \cite{b2} and Fister {\it et al.} \cite{k3} explored an optimal chemotherapy strategy using Pontryagin's Maximum Principle, with a single control that represents the percentage effect it has on viral infectivity (simulating a drug such as AZT (zidovudine)) using dynamical HIV models. Kirschner {\it et al.} \cite{k2} used an existing model which describes the interaction of the immune system with HIV. In Kirschner {\it et al.} \cite{k2} the authors used a single control representing the percentage effect chemotherapy has on viral production (simulating effects of a protease inhibitor). Kutch and Gufil \cite{k7} investigated the reasons underlying the development of drug-insensitive HIV-strains, and demonstrated that optimal drug administration may be useful in increasing patient health by delaying the emergence of drug-resistant mutant viral strains. Kutch and Gufil \cite{k7} used two controls representing the percentage effect chemotherapy has on CD4+ T cells infection and viral production and also incorporated drug efficacy. In the study by Kutch and Gufil \cite{k7}, an alternative approach to Pontryagin's Maximum Principle was adopted, that involves converting the standard optimal control problem into a parameter optimization problem by discretizing the control input vector. An HIV immune dynamics model with three viral strains was used. Joshi \cite{j1} explored optimal control of an ordinary differential equation model taken from \cite{k3}. In the paper \cite{j1}, Joshi considered two controls, one boosting the immune system and the other delaying HIV progression. The novel part of our work is that we explore optimal control of chemotherapy using an HIV dynamical model that incorporates explicit cellular immune response (lytic mechanism and two non-lytic mechanisms). We use two controls, one simulating effect of RTIs and the other control simulating effect of PIs, incorporating drug efficacy. The paper is structured as follows: Firstly in section 2 we formulate a model of HIV immune dynamics, with explicit immune response (lytic and non-lytic components). The model mimics virus and CD4+ T cells dynamics in an infected individual. We modify the model to capture the effects of combined therapy and derive an optimal control problem with an objective functional that maximizes CD4+ T cells and minimizes systemic costs. In section 3 we prove the existence of an optimal control pair and characterize the control pair in section 4. In section 5 we state the optimality system, which is the state system coupled with the adjoint system. In section 6, we prove the uniqueness of the optimality system and we present numerical illustrations for the optimality system in section 7. We make some concluding remarks in section 8. \section{The Model} Let $T$ denote the population density of uninfected CD4+ T cells, $T^*$ the density of infected CD4+ T cells, $V$ the density of free viral particles and $C$ the density of HIV-1 specific cytotoxic T lymphocytes (CTLs). The rate of change of each of these is governed by a first order differential equation. $T$ cell dynamics are governed by proliferation due to virus presence, apoptosis, natural death and thymus supply and viral infectivity inhibited by CTL chemokines. For $T$ the equation is \begin{equation} \frac{dT(t)}{dt} = s_{1} + \frac{rT(t)V(t)}{B_{V}+ V(t)} - e^{-a_{0}C(t)}\beta V(t)T(t)- \mu_{T}T(t) - \frac{k V(t)T(t)}{B_{T} + T(t)}. \label{e1} \end{equation} Here the first term on the right-hand side, $s_{1}$, represents the source of new CD4+ T cells from the thymus \cite{k1}. This is followed by the proliferation term of CD4+ T cells in the presence of the virus: $r$ is the proliferation rate and $B_{V}$ is a parameter that determines the amount of antigen needed to generate half maximal stimulation \cite{k1}. The third term describes the infection of CD4+ T cells by the virus. The presence of CTLs that release chemokines, such as $\beta$- chemokines that block the entry of certain virions into target cells \cite{p1}; \cite{k4}, prevent infection of new cells by a factor $e^{-a_{0}C}$ (effectiveness of CTLs), where $a_{0}$ is the efficiency of each CTL in reducing CD4+ T cells infection. The hypothesis is that reduction of infection of CD4+ T cells is enhanced by the number of HIV-specific CTLs available. The idea goes as follows: as $C\rightarrow\infty$, $e^{-a_{0}C}\rightarrow 0$ meaning that the availability of large quantities of CTLs reduce the rate of infection of CD4+ T cells. The extent of reduction depends on the effectiveness of CTLs ($e^{-a_{0}C}$). Conversely as $C\rightarrow 0$, $e^{-a_{0}C}\rightarrow 1$ meaning that for low CTL count or zero CTL, the infection rate of CD4+ T cells by virus is slightly reduced or not reduced at all. The effectiveness value of CTLs ranges from 0 to 1. We assume that reduction in infection rate has an exponential effect. Here $\beta$ is the rate of infection of CD4+ T cells by the virus. The fourth term is a natural death term, since cells have a finite life span. On average the life span is $1/\mu_{T}$. The last term represents the destruction of CD4+ T cells by the influence of toxic viral proteins. The idea is as follows: The parameter $k$ is the rate of apoptosis. There is a limit to the rate of T cell mortality due to the induction of apoptosis. The limit is a function of variables such as presentation of HIV-1 Env gp120/gp41, receptors involved (especially chemokines CCR5 and CXCR4) and the complexity of target cell contact \cite{a1}. In other words, there is a saturation effect in which the virus can only present itself to so many T cells even when the CD4+ T cell population is low. Conversely, there is an increase in the effect of apoptosis at low CD4+ T cell densities. If T cell density is low, there are more virions per cell and this could lead to higher engagement of apoptosis receptors. On the other hand, if the T cell density is high, there are less virions per cell therefore the chances of virus presentation decreases. Thus presentation exhibits this switching phenomenon and it is this behaviour which is represented by the Hollings Type II function \cite{k5}. The importance of the parameter $B_{T}$, is that it determines the scale at which engagement of apoptosis receptors begins to take effect. The rate of change of the infected CD4+ T cells is governed by the equation \begin{equation} \frac{dT^{*}(t)}{dt} =e^{-a_{0}C(t)}\beta V(t)T(t) - \alpha T^{*}(t) - hT^{*}(t)C(t). \label{e2} \end{equation} The first term on the right-hand side is a gain term for infected cells. The third term is a direct killing of virus infected cells through perforin-granzyme and Fas-FasL pathways. Infected cells are lysed by CTLs at a rate $h$ \cite{k6}. Infected cells are also lost by cytopathic effect of virus and natural death such that they have a finite life span that averages $1/\alpha$. The third equation of the system \begin{equation} \frac{dV(t)}{dt} = N\alpha T^{*}(t)e^{-a_{1}C(t)} - \mu_{V}V(t), \label{e3} \end{equation} describes the rate of change of viral load. The first term on the right-hand side explains the source of the virus. Virions are released by a burst of infected cells \cite{k1}, where an average of $N$ viral particles are released per infected cell. $N \alpha$ is the average rate of virus production per productively infected cell. CTLs release cytokines such as interferon-$\gamma$ (INF-$\gamma$) that can suppress the rate of virus production by virus infected cells \cite{a2}; \cite{s1}. Therefore, they reduce viral burst by a factor of $e^{-a_{1}C}$, where $a_{1}$ is the rate at which each CTL suppresses virus production. The last term describes natural loss of viral particles. The fourth equation \begin{equation} \frac{dC(t)}{dt} = s_{2} + p_{0}T(t)V(t)C(t) - \mu_{C}C(t), \label{e4} \end{equation} describes the dynamics of CTLs during HIV-1 infection. Naive CD8+ T cells differentiate into CTLs when stimulated by helper cells (CD4+ T Cells). HIV-1 specific CTLs decline with increased disease and decreased CD4+ T cell numbers, which means that the CTL population proliferation depends on the stimulation of CD4+ T cells. High numbers of CTLs are associated with low virus titers at equilibrium and loss of CTLs results in an increase in viral load. The first term on the right-hand side, $s_{2}$ models the production rate of HIV specific CD8+ T cells from pre-cursors \cite{k6} and the second term accounts for the differentiation of naive CD8+ T cells into CTLs in response to HIV. Differentiation of CD8+ T cells depends on the help of CD4+ T cells present where $p_{0}$ is the rate of the process. Wodarz and Nowak \cite{w1} used a similar term to model the proliferation of HIV specific CTLs. CTLs are cleared at a rate $\mu_{C}$, a blanket term for death (natural and apoptotic). The model of HIV immune dynamics given by equations (\ref{e1}), (\ref{e2}), (\ref{e3}) and (\ref{e4}) has two steady states in the presence of immune response. The first steady state is the uninfected state given by $$ \Big(\bar T_{un}=\frac{s_{1}}{\mu_{T}},\mbox{ } \bar T^{*}_{un}=0, \mbox{ } \bar V_{un} =0, \mbox{ } \bar C_{un}=\frac{s_{2}}{\mu_{C}}\Big). $$ If infection persists the system converges to a second steady state, an immune controlled equilibrium given by: $$ \bar T_{in}=\frac{\mu_{V}(\alpha+h\bar C_{in})e^{(a_{0}+a_{1}) \bar C_{in}}}{N\alpha \beta},\quad \bar T_{in}^{*}=\frac{\mu_{V}}{N\alpha}e^{a_{1}\bar C_{in}}\bar V_{in},\quad \bar C_{in}=\frac{s_{2}}{\mu_{C}-p\bar T_{in} \bar V_{in}}, $$ and \begin{align*} \bar V_{in} &= \frac{\left(s_{1}+(r-\mu_{T})\bar T_{in}-\beta B_{V} \bar T_{in} e^{-a_{0}\bar C_{in}}-\frac{kB_{V}\bar T_{in}}{B_{T} +\bar T_{in}}\right)}{2\left(\frac{k\bar T_{in}}{B_{T}+\bar T_{in}} +\beta \bar T_{in}e^{-a_{0}\bar C_{in}}\right)} \\ &\quad +\Bigg(\Big(\beta B_{V}\bar T_{in}e^{-a_{0}\bar C_{in}} +\frac{kB_{V}\bar T_{in}}{B_{T}+\bar T_{in}}+(\mu_{T}-r)\bar T_{in}-s_{1} \Big)^{2} \\ &\quad-4\big(\frac{k\bar T_{in}}{B_{T}+\bar T_{in}} +\beta \bar T_{in}e^{-a_{0}\bar C_{in}}\big) \big(\mu_{T}B_{V}\bar T_{in}-s_{1}B_{V}\big)\Bigg)^{1/2}\\ &\quad\div \Big( 2\big(\frac{k\bar T_{in}}{B_{T}+\bar T_{in}}+\beta \bar T_{in}e^{-a_{0} \bar C_{in}}\big)\Big). \end{align*} The virus reproductive number, $R_{0}$ which is the number of newly infected cells that arise from any one infected cell when almost all cells are uninfected, is given by $$ R_{0}=\frac{N\beta \alpha s_{1}e^{-(a_{0}+a_{1})\bar C_{un}}}{\mu_{V} \mu_{T}(\alpha+h\bar C_{un})} $$ where $\bar C_{un}=\frac{s_{2}}{\mu_{C}}$. The reproductive number is governed by several factors including the efficiency with which HIV infects CD4+ T cells, $\beta$ (infectivity constant), number of virions produced by one infected cell (burst size, $N$), rate of virion clearance from the body, $\mu_{V}$, death rate of uninfected CD4+ T cells, $\mu_{T}$, CD4+ T cells production rate, $s_{1}$, effectiveness of CTLs in reducing infection and reducing burst size ($e^{-(a_{0}+a_{1})\bar C_{un}}$), the effect of CTLs in killing virally infected cells, $h\bar C_{un}$ and the the cytopathic effect of the virus, $\alpha$. Determination of stability of equilibrium states give us the following results: if $R_{0}<1$, uninfected equilibrium state is asymptotically stable, that is, infection is abortive. If $R_{0}>1$, the uninfected state is unstable and it converges to an immune controlled equilibrium state that is locally asymptotically stable. The virus will spread after infection and the abundance of uninfected cells, infected cells, free viruses and CTLs is given by equations in $\bar T_{in}$, $\bar T^{*}_{in}$, $\bar V_{in}$ and $\bar C_{in}$ respectively. If $R_{0}=1$ the uninfected state and the infected state coincide. If $R_{0}>1$ infection persists, then it will eventually leads to the acquired immune deficiency syndrome (AIDS) stage, associated with a weakened immune system which has difficulty fighting off opportunistic infections \cite{s2}. It is at this stage when therapy is initiated to boost the health of infected individuals. After initiation of combined chemotherapy, combination of RTIs and PIs, infection rate of CD4+ T cells is reduced and the number of viral particles produced by an actively infected CD4+ T cell is reduced. If we let $u_{RTI}(t)$ represent the normalized RTI dosage as a function of time, then $\beta$ will be modified to become $(1-\frac{1}{2}u_{RTI}(t))\beta$ where $\frac{1}{2}$ models drug efficacy \cite{k7}) and it is meant to take into account the effectiveness of the delivery. If we also let $u_{PI}(t)$ be the normalized PI dosage, then the parameter $N$ will be modified to become $(1-\frac{1}{2}u_{PI}(t))N$ \cite{k7}. Hence the state system becomes \begin{equation} \label{e5} \begin{gathered} \begin{aligned} \frac{dT(t)}{dt} &= s_{1} + \frac{rT(t)V(t)}{B_{V}+ V(t)} - (1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}V(t)T(t) \\ &\quad - \mu_{T}T(t) - \frac{k V(t)T(t)}{B_{T} + T(t)} \end{aligned} \\ \frac{dT^{*}(t)}{dt}=(1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}V(t)T(t) - \alpha T^{*}(t) - hT^{*}(t)C(t) \\ \frac{dV(t)}{dt}=(1-\frac{1}{2}u_{PI}(t))Ne^{-a_{1}C(t)}\alpha T^{*}(t) - \mu_{V}V(t) \\ \frac{dC(t)}{dt} = s_{2}+ p_{0}T(t)V(t)C(t) - \mu_{C}C(t). \end{gathered} \end{equation} The controls $u_{RTI}(t)$ and $u_{PI}(t)$ represent the action of RTI (viral infectivity reduction) and PI (viral replication suppression) drugs respectively.\\ The objective functional is defined as, \begin{equation} J(u_{RTI}, u_{PI}) = \int ^{T_{f}}_{0} \big[T(t)-\big(\frac{A_{1}}{2}u^{2}_{RTI}(t) + \frac{A_{2}}{2}u^{2}_{PI}(t)\big)\big]dt \label{e6} \end{equation} where $T(t)$ is the benefit based on CD4+ T cells and the other terms are systemic costs of the drug treatments. The benefit of treatment is based on an increase of CD4+T cells and systemic costs of drugs are minimized. The positive constants $A_{1}$ and $A_{2}$ represent desired weight on the benefit and cost, and $u^{2}_{RTI}$, $u^{2}_{PI}$ reflect the severity of the side effects of the drugs \cite{j1}. The cost function is assumed to be nonlinear, basing on the fact that there is no linear relationship between the effects of treatment on CD4+ T cells or viral load hence the choice of a quadratic cost function \cite{k2}. We impose a condition for treatment time, $t\in [0, T_{f}]$, limited treatment window \cite{b2}, that monitors global effects of these phenomena; treatment lasts for a given period of time because HIV can mutate and develop resistance to treatment after some finite time frame and in addition treatment has potentially harmful side effects, and these side effects increase with duration of treatment. The time $t=0$ is the time when treatment is initiated and time $t=T_{f}$ is the time when treatment is stopped. The main objective is to maximize the benefit based on the CD4+ T cell count (increase in quality of life) and the systemic cost based on the percentage effect of the chemotherapy given (RTIs and PIs) is being minimized (toxic side effects being avoided as much as possible and not causing patient death). We seek an optimal control pair, $u^{*}_{RTI}$, $u^{*}_{PI}$ such that \begin{equation} J(u^{*}_{RTI}, u^{*}_{PI}) = \max \left\{J(u_{RTI}, u_{PI})|(u_{RTI}, u_{PI})\in U \right\} \label{e7} \end{equation} where \begin{align*} U = \big\{&(u_{RTI}, u_{PI}), u_{RTI}, u_{PI} \mbox{ measurable, }0 \leq a_{11}\leq u_{RTI}\leq b_{11}\leq 1 \\ &\mbox{ and } 0 \leq a_{22}\leq u_{PI} \leq b_{22}\leq 1\big\} \end{align*} is the control set where $t \in [0, T_{f}]$. The basic framework of this problem is to characterize the optimal control and prove the existence of the optimal control and uniqueness of the optimality system. \section{Existence of an Optimal Control Pair} The existence of the optimal control pair can be obtained using a result by Joshi \cite{j1}, Fister {\it et al.} \cite{f1}, and other references quoted therein. \begin{theorem} Given the objective functional $$ J(u_{RTI}, u_{PI}) = \int ^{T_{f}}_{0}\left[T(t) -\left(\frac{A_{1}}{2}u^{2}_{RTI}(t) + \frac{A_{2}}{2}u^{2}_{PI}(t) \right)\right]dt\,, $$ where $U = \{(u_{RTI}(t), u_{PI}(t))$ , piecewise continuous such that $0< a_{11}\leq u_{RTI}(t)\leq b_{11} < 1$, $0 < a_{22} \leq u_{PI}(t) \leq b_{22} < 1\}$ for all $t \in [0, T_{f}]$ subject to equations of system \eqref{e5} with $T(0)=T_{0}$, $T^{*}(0)=T^{*}_{0}$, $V(0)=V_{0}$ and $C(0)=C_{0}$, then there exists an optimal control pair $u^{*}_{RTI}$, $u^{*}_{PI}$ such that $$\max\{J(u_{RTI},u_{PI})|(u_{RTI},u_{PI}) \in U\} = J(u^{*}_{RTI},u^{*}_{PI})$$ if the following conditions are met: \begin{enumerate} \item The class of all initial conditions with an optimal control pair $u_{RTI}$, $u_{PI}$ in the admissible control set along with each state equation being satisfied is not empty. \item The admissible control set $U$ is closed and convex. \item Each right hand side of equations of system (\ref{e5}) is continuous, is bounded above by a sum of the bounded control and the state, and can be written as a linear function of an optimal control pair $u_{RTI}$, $u_{PI}$ with coefficients depending on time and the state. \item The integrand $J(u_{RTI},u_{PI})$ is concave. \item The integrand $J(u_{RTI},u_{PI})$ is bounded above by $C_{2}- C_{1}(|u_{RTI}|^{2} + |u_{PI}|^{2})$ with $C_{1} >0$. \end{enumerate} \end{theorem} \begin{proof} Our definition of the control set satisfies conditions 1 and 2. For the model to be realistic, we impose the restrictions that CD4+ T cells and CD8+ T cells do not grow unbounded, so we use $T(t)0$, $T_{\rm max}>0$ and $0 < e^{-a_{0}C_{\rm max}} <1$. $$ \frac{d\bar V}{dt}=N\alpha e^{-a_{1}C_{\rm max}}\bar T^{*}, \quad \bar V(0)=\bar T_{0}, $$ where $N>0$, $\alpha >0$ and $0 < e^{-a_{1}C_{\rm max}} < 1$. Since this system is linear in finite time with bounded coefficients, then the supersolutions $\bar T^{*}$ and $\bar V$ are uniformly bounded. Since our state system is bilinear in $u_{RTI}$ and $u_{PI}$, the right hand side of equations of system (\ref{e5}) satisfies condition 3. The right hand side of system (\ref{e5}) is continuous and it can be written as: $$ \mathbf{f}(t,\mathbf{T},\mathbf{u})=\boldsymbol{\alpha}(t,\mathbf{T}) +\boldsymbol{\gamma}(t,\mathbf{T})\mathbf{u} $$ and the boundedness of solutions gives $$ |\mathbf{f}(t,\mathbf{T},\mathbf{u})|\leq C_{1}(1+|\mathbf{T}|+|\mathbf{u}|) $$ for $0\leq t \leq T_{f}$ where $\mathbf{T}\in \Re^{4}$, $\mathbf{u}\in \Re^{2}$ where $\mathbf{T}=(T, T^{*}, V, C)$ and $\mathbf{u}=(u_{RTI},u_{PI}))$ and $C_{1}$ depends on the coefficients of the system. The vectors $\boldsymbol{\alpha}$ and $\boldsymbol{\gamma}$ are vector-valued functions of $\mathbf{T}$. In order to verify the convexity of the integrand of our objective functional, \textsl{J} we show that \begin{equation} \textsl{J}(t,\mathbf{T}, (1-\epsilon)\mathbf{u}+\epsilon \mathbf{v})\geq (1-\epsilon)\textsl{J}(t, \mathbf{T}, \mathbf{u}) + \epsilon \textsl{J}(t, \mathbf{T}, \mathbf{v}) \label{e8} \end{equation} for $0 < \epsilon < 1$ and $\textsl{J}(t, \mathbf{T}, \mathbf{u}) = T-(\frac{A_{1}}{2}u^{2}_{RTI} + \frac{A_{2}}{2}u^{2}_{PI})$. \begin{align*} &\textsl{J}(t,\mathbf{T}, (1-\epsilon)\mathbf{u}+\epsilon \mathbf{v})\\ &=\left[T-\frac{A_{1}}{2}\left((1-\epsilon)u_{RTI} + \epsilon v_{RTI}\right)^{2} -\frac{A_{2}}{2}\left((1-\epsilon)u_{PI}+ \epsilon v_{PI}\right)^{2}\right] \\ &=T-\frac{A_{1}}{2}\left[u^{2}_{RTI}-2\epsilon u^{2}_{RTI} \epsilon^{2} u^{2}_{RTI}+2(1-\epsilon)u_{RTI}\epsilon v_{RTI} + \epsilon^{2}v^{2}_{RTI}\right] \\ &\quad -\frac{A_{2}}{2}\left[u_{PI}^{2}-2\epsilon u^{2}_{PI} + \epsilon ^{2} u_{PI}^{2} + 2(1-\epsilon)u_{PI}\epsilon v_{PI} + \epsilon^{2}v^{2}_{PI}\right] \\ &= T-(\frac{A_{1}}{2}u_{RTI}^{2}+\frac{A_{2}}{2}u_{PI}^{2}) \\ &\quad -\frac{A_{1}}{2}[(\epsilon^{2}-2\epsilon)u_{RTI}^{2}+\epsilon^{2}v_{RTI}^{2}+2\epsilon(1-\epsilon)u_{RTI}v_{RTI}] \nonumber \\ &\quad -\frac{A_{2}}{2}[(\epsilon^{2}-2\epsilon)u_{PI}^{2}+\epsilon^{2} v_{PI}^{2}+2\epsilon(1-\epsilon)u_{PI}v_{PI}]. \end{align*} \begin{equation} \begin{aligned} &(1-\epsilon)\textsl{J}(t, \mathbf{T}, \mathbf{u}) +\epsilon \textsl{J}(t, \mathbf{T}, \mathbf{v}) \\ &=(1-\epsilon)[T-(\frac{A_{1}}{2}u_{RTI}^{2}+\frac{A_{2}}{2}u_{PI}^{2})] + \epsilon [T-(\frac{A_{1}}{2}v_{RTI}^{2}+\frac{A_{2}}{2}v_{PI}^{2})] \\ &=T-(\frac{A_{1}}{2}u_{RTI}^{2}+\frac{A_{2}}{2}u_{PI}^{2}) -\epsilon[T-(\frac{A_{1}}{2}u_{RTI}^{2}+ \frac{A_{2}}{2}u_{PI}^{2})]\\ &\quad +\epsilon [T-(\frac{A_{1}}{2}v_{RTI}^{2}+\frac{A_{2}}{2}v_{PI}^{2})]\\ &=T-(\frac{A_{1}}{2}u_{RTI}^{2}+\frac{A_{2}}{2}u_{PI}^{2}) -\frac{\epsilon}{2}(-A_{1}u_{RTI}^{2}-A_{2}u_{PI}^{2}+A_{1}v_{RTI}^{2} +A_{2}v_{PI}^{2}). \end{aligned}\label{e9} \end{equation} Thus to show that $\textsl{J}(t, \mathbf{T},.)$ is concave in $U$, we note that the following inequality holds \begin{equation} \begin{aligned} &\frac{A_{1}}{2}[(\epsilon^{2}-2\epsilon)u_{RTI}^{2}+\epsilon^{2}v_{RTI}^{2} +2\epsilon(1-\epsilon)u_{RTI}v_{RTI}]\\ &+\frac{A_{2}}{2}[(\epsilon^{2} -2\epsilon)u_{PI}^{2}+\epsilon^{2}v_{PI}^{2}+2\epsilon(1-\epsilon)u_{PI}v_{PI}] \\ &\leq \frac{\epsilon}{2}(-A_{1}u_{RTI}^{2}-A_{2}u_{PI}^{2}+A_{1}v_{RTI}^{2} +A_{2}v_{PI}^{2}). \end{aligned}\label{e10} \end{equation} This implies \begin{align*} & \frac{A_{1}}{2}\epsilon^{2}u_{RTI}^{2}-A_{1}\epsilon u_{RTI}^{2} +\frac{A_{1}}{2}\epsilon^{2}v_{RTI}^{2}+A_{1} \epsilon(1-\epsilon)u_{RTI}v_{PI}+\frac{A_{2}}{2}\epsilon^{2}u_{PI}^{2} -A_{2}\epsilon u_{PI}^{2}\\ &+\frac{A_{1}}{2}\epsilon^{2}v_{PI}^{2} + A_{2}\epsilon(1-\epsilon)u_{PI}v_{PI}+\frac{A_{1}}{2} \epsilon u_{RTI}^{2}+ \frac{A_{2}}{2}\epsilon u_{PI}^{2} -\frac{A_{1}}{2}\epsilon v_{RTI}^{2}-\frac{A_{2}}{2}\epsilon v_{PI}^{2}\leq 0. \end{align*} Finally this gives \begin{align*} &\frac{A_{1}}{2}(\epsilon^{2}-\epsilon)(u_{RTI}^{2}+v_{RTI}^{2}) +\frac{A_{2}}{2}(\epsilon^{2}-\epsilon)(u_{PI}^{2}+v_{PI}^{2})\\ &+ \epsilon(1-\epsilon)(A_{1}u_{RTI}v_{RTI}+A_{2}u_{PI}v_{PI})\leq 0, \end{align*} which is equivalent to \begin{align*} &\frac{A_{1}}{2}(\epsilon^{2}-\epsilon)(u_{RTI}^{2}+v_{RTI}^{2}) +(\epsilon-\epsilon^{2})A_{1}u_{RTI}v_{RTI}\\ &+\frac{A_{2}}{2}(\epsilon^{2}-\epsilon)(u_{PI}^{2}+v_{PI}^{2}) +(\epsilon-\epsilon^{2})A_{2}u_{PI}v_{PI} \leq 0 \end{align*} which can be written as \begin{equation} \begin{aligned} &-\frac{A_{1}}{2}\left(\sqrt{\epsilon(1-\epsilon)}u_{RTI} -\sqrt{\epsilon(1-\epsilon)}v_{RTI}\right)^{2}\\ &-\frac{A_{2}}{2}\left(\sqrt{\epsilon(1-\epsilon)}u_{PI} -\sqrt{\epsilon(1-\epsilon)}v_{PI}\right)^{2}\leq 0. \end{aligned}\label{e11} \end{equation} This holds since $A_{1}$, $A_{2} >0$, hence equation \ref{e8} holds. Finally we need to show that $\textsl{J}(t, \mathbf{T}, \mathbf{u}) \leq C_{2}-C_{1}|\mathbf{u}|^{\beta} $, where $C_{1} > 0$ and $\beta > 1$. For our case $$ \textsl{J}(t, \mathbf{T}, \mathbf{u})=T -\big(\frac{A_{1}}{2}u_{RTI}^{2}+ \frac{A_{2}}{2}u_{PI}^{2}\big) \leq C_{2}-C_{1}|\mathbf{u}|^{2} $$ where $C_{2}$ depends on the upper bound on CD4+ T cells, $T$, and $C_{1}>0$ since $A_{1}$, $A_{2}>0$. We conclude that there exists an optimal control pair. \end{proof} \section{Characterization} Since there exists an optimal control pair for maximizing the functional, equation (\ref{e6}), subject to system (\ref{e5}) we derive necessary conditions on the optimal control pair \cite{f1}. We discuss the theorem that relates to the characterization of the optimal control. In order to derive the necessary conditions for this optimal control pair, we use Pontryagin's Maximum Principle \cite{k5}. The Lagrangian is defined as \begin{equation} \begin{aligned} L =& T(t)-\left(\frac{A_{1}}{2}u^{2}_{RTI}(t)+ \frac{A_{2}}{2}u^{2}_{PI}(t)\right)+ \lambda_{1}\Big[s_{1} + \frac{rT(t)V(t)}{B_{V}+ V(t)} \nonumber \\ &- (1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}V(t)T(t)-\mu_{T}T(t) - \frac{k V(t)T(t)}{B_{T} + T(t)} \Big] \nonumber \\ &+ \lambda_{2}\left[ (1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}V(t)T(t) - \alpha T^{*}(t) - hT^{*}(t)C(t) \right] \nonumber \\ &+ \lambda_{3}\left[(1-\frac{1}{2}u_{PI}(t))Ne^{-a_{1}C(t)}\alpha T^{*}(t) - \mu_{V}V(t)\right] \nonumber \\ &+ \lambda_{4}\left[ s_{2}+p_{0}T(t)V(t)C(t) - \mu_{C}C(t)\right] \nonumber \\ &+ w_{11}(t)(b_{11} - u_{RTI}(t)) + w_{12}(t)(u_{RTI}(t) - a_{11}) \nonumber \\ &+ w_{21}(t)(b_{22} - u_{PI}(t)) + w_{22}(t)(u_{PI}(t) - a_{22})\,, \end{aligned}\label{e12} \end{equation} where $w_{11}(t)\geq 0$, $w_{12}(t)\geq 0$, $w_{21}(t)\geq 0$, $w_{22}(t)\geq 0$ are penalty multipliers satisfying $w_{11}(t)(b_{11}-u_{RTI}(t))=0$, $w_{12}(t)(u_{RTI}(t)-a_{11})=0$ at the optimal $u^{*}_{RTI}$, and $w_{21}(t)(b_{22}-u_{PI}(t))=0$, $w_{22}(t)(u_{PI}(t)-a_{22})=0$ at the optimal $u^{*}_{PI}$. \begin{theorem} Given a pair of optimal controls $u^{*}_{RTI}$, $u^{*}_{PI}$ and solutions $T, T^{*}, V, C$ of the corresponding state system \eqref{e5}, there exists adjoint variables $\lambda_{i}$ for $i=1,2,3,4 $ satisfying the following canonical equations \begin{equation} \begin{gathered} \begin{aligned} \frac{d \lambda_{1}}{dt} &=-\frac{\partial L}{\partial T}\\ &=-\Big[1+\lambda_{1}\Big(\frac{rV(t)}{B_{V}+V(t)}-(1-\frac{1}{2}u_{RTI}(t)) \beta e^{-a_{0}C(t)}V(t)-\mu_{T}-\frac{kV(t)B_{T}}{(B_{T}+T(t))^{2}}\Big)\Big] \\ &\quad - \Big[\lambda_{2}((1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}V(t)) + \lambda_{4}p_{0}V(t)C(t)\Big] \end{aligned} \\ \frac{d \lambda_{2}}{dt} =-\frac{\partial L}{\partial T^{*}} =-\Big[ \lambda_{2}(-\alpha -hC(t))+ \lambda_{3}((1-\frac{1}{2}u_{PI}(t)) Ne^{-a_{1}C(t)}\alpha)\Big]\\ \begin{aligned} \frac{d \lambda_{3}}{dt} &=-\frac{\partial L}{\partial V}\\ &=-\Big[\lambda_{1}\left(\frac{rT(t)B_{V}}{(B_{V}+V(t))^{2}} -(1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}T(t) -\frac{kT(t)}{B_{T}+T(t)}\right)\Big] \\ &\quad -\Big[\lambda_{2}((1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}T(t)) - \lambda_{3}\mu_{V} \Big] \end{aligned} \\ \begin{aligned} \frac{d \lambda_{4}}{dt} &=-\frac{\partial L}{\partial C}\\ &=-\left[\lambda_{1}(a_{0}(1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}V(t)T(t))\right] \nonumber \\ &\quad+ \Big[\lambda_{2}(a_{0}(1-\frac{1}{2}u_{RTI}(t))\beta e^{-a_{0}C(t)}V(t)T(t) +hT^{*}(t))\Big] \\ &\quad+\Big[\lambda_{3}(a_{1}(1-\frac{1}{2}u_{PI}(t))N e^{-a_{1}C(t)} \alpha T^{*}(t)) - \lambda_{4}(p_{0}T(t)V(t)-\mu_{C})\Big] \end{aligned} \end{gathered} \label{e13} \end{equation} with transversality conditions $\lambda_{i}(T_{f})=0$ for $i=1,2,3,4$. Further, the following characterization holds: \begin{gather*} u^{*}_{RTI}(t)=\min \big\{\max\big\{a_{11},\frac{1}{2A_{1}} (\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)\big\}, b_{11}\big\}, \\ u^{*}_{PI}(t)=\min \big\{\max \big\{a_{22},\frac{-\lambda_{3}}{2A_{2}}N e^{-a_{1}C(t)}\alpha T^{*}(t)\big\},b_{22}\big\}. \end{gather*} \end{theorem} \begin{proof} The form of the adjoint equations and transversality conditions are standard results from Pontryagin's Maximum Principle \cite{j1}; therefore, solutions to the adjoint system exists and are bounded. To determine the interior maximum of our Lagrangian, we take the partial derivatives of $L$ with respect to $u_{RTI}$ and $u_{PI}$ and set it equal to zero. Thus \begin{align*} \frac{\partial L}{\partial u_{RTI}} &=-A_{1}u^{*}_{RTI}(t) + \frac{\lambda_{1}}{2}\beta e^{-a_{0}C(t)}V(t)T(t) -\frac{\lambda_{2}}{2}\beta e^{-a_{0}C(t)}V(t)T(t)\\ &\quad -w_{11}(t)+w_{12}(t)=0 \quad \mbox{at } u^{*}_{RTI}. \end{align*} \[ \frac{\partial L}{\partial u_{PI}}=-A_{2}u^{*}_{PI}(t) -\frac{\lambda_{3}}{2}N e^{-a_{1}C(t)}\alpha T^{*}(t)-w_{21}(t)+w_{22}(t)=0 \quad \mbox{at } u^{*}_{PI}. \] Hence upon simplification, we obtain \begin{gather} u^{*}_{RTI}(t)=\frac{\frac{(\lambda_{1}-\lambda_{2})}{2}\beta e^{-a_{0}C(t)} V(t)T(t)-w_{11}(t)+w_{12}(t)}{A_{1}}\label{e14} \\ u^{*}_{PI}(t)=\frac{\frac{-\lambda_{3}}{2}N e^{-a_{1}C(t)}\alpha T^{*}(t) -w_{21}(t)+ w_{22}(t)}{A_{2}}\label{e15} \end{gather} \subsection{Case $u^{*}_{RTI}$} \begin{enumerate} \item On the set $\{t|a_{11}< u^{*}_{RTI}(t)< b_{11}\}$, $w_{11}(t)=w_{12}(t)=0$. From (\ref{e14}) we have $$ u^{*}_{RTI}(t)=\frac{(\lambda_{1}-\lambda_{2}) \beta e^{-a_{0}C(t)}V(t)T(t)}{2A_{1}} $$ \item On the set $\{t|u^{*}_{RTI}(t)=a_{11}\}$, $w_{11}(t)=0$. Consequently, $$ u^{*}_{RTI}(t)=a_{11}=\frac{(\lambda_{1}-\lambda_{2}) \beta e^{-a_{0}C(t)}V(t)T(t)}{2A_{1}} + \frac{w_{12}(t)}{A_{1}} $$ or $$ \frac{(\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)}{2A_{1}} \leq a_{11},\quad \mbox{since } w_{12}(t)\geq 0. $$ \item On the set $\{t|u^{*}_{RTI}(t)=b_{11}\}$, $w_{12}(t)=0$. Consequently, $$ u^{*}_{RTI}(t)=b_{11}=\frac{(\lambda_{1}-\lambda_{2}) \beta e^{-a_{0}C(t)}V(t)T(t)}{2A_{1}} - \frac{w_{11}(t)}{A_{1}} $$ or $$ \frac{(\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)}{2A_{1}} \geq b_{11},\quad \mbox{since } w_{11}(t)\geq 0. $$ Combining all the three cases in a compact form gives $$ u^{*}_{RTI}(t)=\min \big\{\max \big\{a_{11},\frac{1}{2A_{1}} (\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)\big\},b_{11}\big\}. $$ \end{enumerate} \subsection{Case $u^{*}_{PI}$} \begin{enumerate} \item On the set $\{t|a_{22}< u^{*}_{PI}(t)< b_{22}\}$, $w_{21}(t)=w_{22}(t)=0$. From (\ref{e15}) we have $$ u^{*}_{PI}(t)=-\frac{\lambda_{3}Ne^{-a_{1}C(t)}\alpha T^{*}(t)}{2A_{2}}. $$ \item On the set $\{t|u^{*}_{PI}(t)=a_{22}\}$, $w_{21}(t)=0$. Consequently, $$ u^{*}_{PI}(t)=a_{22}=-\frac{\lambda_{3}Ne^{-a_{1}C(t)}\alpha T^{*}(t)}{2A_{2}} + \frac{w_{22}(t)}{A_{2}} $$ or $$ -\frac{\lambda_{3}Ne^{-a_{1}C(t)}\alpha T^{*}(t)}{2A_{2}} \leq a_{22}, \quad \mbox{since } w_{22}(t)\geq 0. $$ \item On the set $\{t|u^{*}_{PI}(t)=b_{22}\}$, $w_{22}(t)=0$. Consequently, $$ u^{*}_{PI}(t)=b_{22}=-\frac{\lambda_{3}Ne^{-a_{1}C(t)}\alpha T^{*}(t)}{2A_{2}} - \frac{w_{21}(t)}{A_{2}} $$ or $$-\frac{\lambda_{3}Ne^{-a_{1}C(t)}\alpha T^{*}(t)}{2A_{2}} \geq b_{22}, \quad \mbox{since } w_{21}(t)\geq 0. $$ Combining all the three cases in compact form gives $$ u^{*}_{PI}(t)=\min \left\{\max \left\{a_{22}, \frac{-\lambda_{3}}{2A_{2}}N e^{-a_{1}C(t)}\alpha T^{*}(t)\right\},b_{22}\right\}. $$ \end{enumerate} \end{proof} \section{Optimality System} Incorporating the presentation of the optimal treatment controls, we have the state system coupled with the adjoint system. \begin{gather*} \begin{aligned} \frac{dT(t)}{dt} &=s_{1}+\frac{rT(t)V(t)}{(B_{V}+V(t))}-\mu_{T}T(t) -\frac{kV(t)T(t)}{(B_{T}+T(t))} \\ &\quad -\Big(1-\frac{1}{2}\min \{\max \{a_{11},\frac{1}{2A_{1}} (\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)\},b_{11}\}\Big)\\ &\quad\times \beta e^{-a_{0}C(t)}V(t)T(t) \end{aligned} \\ \begin{aligned} \frac{dT^{*}(t)}{dt} &= \Big(1-\frac{1}{2}\min \{\max \{a_{11},\frac{1}{2A_{1}}(\lambda_{1} -\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)\},b_{11}\}\Big)\\ &\quad \times \beta e^{-a_{0}C(t)}V(t)T(t) - \alpha T^{*}(t)-hT(t)C(t) \end{aligned} \\ \begin{aligned} \frac{dV(t)}{dt}&= \Big(1-\frac{1}{2}\min \{\max \{a_{22}, \frac{-\lambda_{3}}{2A_{2}}N e^{-a_{1}C(t)}\alpha T^{*}(t)\},b_{22}\}\Big)\\ &\quad \times Ne^{-a_{1}C(t)}\alpha T^{*}(t)-\mu_{V}V(t) \end{aligned} \\ \frac{dC(t)}{dt}=s_{2}+p_{0}T(t)V(t)C(t)-\mu_{C}C(t) \\ \begin{aligned} \frac{d\lambda_{1}}{dt} &=-1-\lambda_{1}\left(\frac{rV(t)}{(B_{V}+V(t))} -\mu_{T}-\frac{kV(t)B_{T}}{(B_{T}+T(t))^{2}}\right) -\lambda_{4}p_{0}C(t)V(t) \\ &\quad +\lambda_{1}\Big(1-\frac{1}{2}\min \{\max \{a_{11},\frac{1}{2A_{1}} (\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)\},b_{11}\}\Big)\\ &\quad \times \beta e^{-a_{0}C(t)}V(t) \\ &\quad -\lambda_{2}\Big(\Big(1-\frac{1}{2}\min \{\max \{a_{11}, \frac{1}{2A_{1}}(\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)\}, b_{11}\}\Big)\\ &\quad\times \beta e^{-a_{0}C(t)}V(t)\Big) \end{aligned} \\ \begin{aligned} \frac{d \lambda_{2}}{dt} &=\lambda_{2}(\alpha+hC(t)) \\ &\quad -\lambda_{3}\Big(1-\frac{1}{2}\min \{\max \{a_{22},\frac{-\lambda_{3}} {2A_{2}}N e^{-a_{1}C(t)}\alpha T^{*}(t)\},b_{22}\}\Big)Ne^{-a_{1}C(t)}\alpha \end{aligned} \\ \begin{aligned} \frac{d\lambda_{3}}{dt} &=-\lambda_{1}\Big(\frac{rT(t)B_{V}}{(B_{V}+V(t))^{2}} -\frac{kT(t)}{B_{T}+T(t)}\Big)+ \lambda_{3}\mu_{V} \\ &\quad +\lambda_{1}\Big(1-\frac{1}{2}\min \{\max \{a_{11},\frac{1}{2A_{1}}(\lambda_{1}-\lambda_{2})\\ &\quad \times \beta e^{-a_{0}C(t)}V(t)T(t)\},b_{11}\}\beta e^{-a_{0}C(t)}T(t)\Big) \\ &\quad-\lambda_{2}\Big(\Big(1-\frac{1}{2}\min \{\max\{a_{11},\frac{1}{2A_{1}} (\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)\},b_{11}\}\Big)\\ &\quad\times \beta e^{-a_{0}C(t)}T(t)\Big) -\lambda_{4}p_{0}T(t)C(t) \end{aligned} \end{gather*} \begin{equation} \begin{aligned} \frac{d \lambda_{4}}{dt} &=-\lambda_{1}\Big(a_{0}\Big(1-\frac{1}{2}\min \{\max\{a_{11},\frac{1}{2A_{1}} (\lambda_{1}-\lambda_{2})\beta e^{-a_{0}C(t)}V(t)T(t)\},b_{11}\}\Big)\\ &\quad\times \beta e^{-a_{0}C(t)}V(t)T(t)\Big) \\ &\quad +\lambda_{2}\Big(a_{0}\Big(1-\frac{1}{2}\min \{\max \{a_{11},\frac{1}{2A_{1}}(\lambda_{1}-\lambda_{2}) \beta e^{-a_{0}C(t)}V(t)T(t)\},b_{11}\}\Big)\\ &\quad\times \beta e^{-a_{0}C(t)}V(t)T(t)\Big) \\ &\quad + \lambda_{3}\Big(a_{1}\Big(1-\frac{1}{2}\min \{\max \{a_{22},\frac{-\lambda_{3}}{2A_{2}}N e^{-a_{1}C(t)}\alpha T^{*}(t)\}, b_{22}\}\Big)\\ &\quad\times Ne^{-a_{1}C(t)}\alpha T^{*}(t)\Big) -\lambda_{4}(p_{0}T(t)V(t)-\mu_{C}) +\lambda_{2}hT^{*}(t) \end{aligned} \label{e16} \end{equation} with $T(0)=T_{0}$, $T^{*}(0)=T^{*}_{0}$, $V(0)=V_{0}$, $C(0)=C_{0}$, $\lambda_{i}(T_{f})=0$ for $i=1,2,3,4$. \section{Uniqueness of the Optimality System} Since the state system moves forward in time and the adjoint system moves backward in time, we have a challenge with uniqueness. To prove uniqueness of solutions of the optimality system for the small time interval, we use the following theorems \cite{j1}. \begin{theorem} The function $u^{*}(c)=\min (\max (c,a),b)$ is Lipschitz continuous in $c$, where $a0$ is chosen. Further we let \begin{gather*} u^{*}_{RTI}(t)=\min \{\max \{a_{11},\frac{1}{2A_{1}}(w-z)\beta e^{-a_{0} e^{mt}x}pq\},b_{11}\}, \\ u^{*}_{PI}(t)=\min \{\max \{a_{22},\frac{-\alpha N}{2A_{2}} e^{-a_{1}e^{-mt}x}vp^{*}\},b_{22}\} \end{gather*} and \begin{gather*} \bar u^{*}_{RTI}(t)=\min \{\max \{a_{11},\frac{1}{2A_{1}}(\bar w-\bar z) \beta e^{-a_{0}e^{mt}\bar x}\bar p\bar q\},b_{11}\}, \\ \bar u^{*}_{PI}(t)=\min \{\max \{a_{22},\frac{-\alpha N}{2A_{2}} e^{-a_{1}e^{-mt}\bar x}\bar v\bar p^{*}\},b_{22}\}. \end{gather*} For the first equation of system (\ref{e16}) we substitute $T=e^{mt}p$ and get \begin{align*} e^{mt}(\dot p + mp) &=s_{1} + \frac{re^{2mt}pq}{B_{V}+e^{mt}q}-\beta e^{-ae^{mt}x}e^{2mt}pq + \frac{1}{2}\beta e^{-a_{0}e^{mt}x}e^{2mt}pqu^{*}_{RTI}\\ &\quad -\mu_{T}e^{mt}p-\frac{ke^{2mt}pq}{B_{T}+e^{mt}p} \end{align*} and for $\bar T=e^{mt}\bar p$ we have \begin{align*} e^{mt}(\dot {\bar p} + m\bar p) &=s_{1} + \frac{re^{2mt}\bar p\bar q}{B_{V}+e^{mt}\bar q} -\beta e^{-a_{0}e^{mt}\bar x}e^{2mt}\bar p\bar q + \frac{1}{2}\beta e^{-a_{0}e^{mt}\bar x}e^{2mt}\bar p\bar q\bar u^{*}_{RTI}\\ &\quad -\mu_{T}e^{mt}\bar p-\frac{ke^{2mt}\bar p\bar q}{B_{T}+e^{mt}\bar p}. \end{align*} Subtracting the expression for $\bar T$ from the expression for $T$ we have \begin{align*} &(\dot p-\dot {\bar p}) + m(p-\bar p)\\ &=re^{mt}\Big(\frac{pq}{B_{V}+e^{mt}q}-\frac{\bar p\bar q}{B_{V} +e^{mt}\bar q}\Big)-\beta e^{mt}\left(e^{-a_{0}e^{mt}x}pq-e^{-a_{0}e^{mt}\bar x}\bar p\bar q\right) \nonumber \\ &\quad +\frac{1}{2}\beta e^{mt}\left(e^{-a_{0}e^{mt}x}u^{*}_{RTI}pq-e^{-a_{0} e^{mt}\bar x}\bar u^{*}_{RTI}\bar p\bar q\right) \\ &\quad -\mu_{T}(p-\bar p)-ke^{mt} \Big(\frac{pq}{B_{T}+e^{mt}p}-\frac{\bar p\bar q}{B_{T}+e^{mt}\bar p}\Big). \end{align*} Multiplying by $(p-\bar p)$ and integrating from $t=0$ to $t=T_{f}$ we have \begin{equation} \begin{aligned} &\frac{1}{2}(p-\bar p)^{2}(T_{f})+m\int_{0}^{T_{f}}(p-\bar p)^{2}dt\\ &= r\int_{0}^{T_{f}}e^{mt} \Big(\frac{pq}{B_{V}+e^{mt}q}-\frac{\bar p\bar q}{B_{V}+e^{mt}\bar q}\Big) (p-\bar p)dt-\mu_{T}\int_{0}^{T_{f}}(p-\bar p)^{2}dt \\ &\quad -\beta \int_{0}^{T_{f}}e^{mt} \left(e^{-a_{0}e^{mt}x}pq-e^{-a_{0}e^{mt} \bar x}\bar p\bar q\right)(p-\bar p)dt \\ &\quad - k\int_{0}^{T_{f}}e^{mt}\left(\frac{pq}{B_{T}+e^{mt}p} -\frac{\bar p\bar q}{B_{T}+e^{mt}\bar p}\right)(p-\bar p)dt \\ &\quad + \frac{\beta}{2}\int_{0}^{T_{f}}e^{mt} \left(e^{-a_{0}e^{mt}x}pqu^{*}_{RTI}-e^{-a_{0}e^{mt}\bar x}\bar p\bar q \bar u^{*}_{RTI}\right)(p-\bar p)dt. \end{aligned} \label{e17} \end{equation} Similarly for $\lambda_{1}=e^{-mt}w$ and $\bar \lambda_{1}=e^{-mt}\bar w$ we have \begin{align*} -\dot w +mw &=e^{mt} + \frac{rwqe^{mt}}{B_{V}+e^{mt}q} -w\beta qe^{-a_{0}e^{mt}x}e^{mt}+\frac{1}{2}\beta we^{-a_{0}e^{mt}x} e^{mt}qu^{*}_{RTI} \\ &\quad -\mu_{T}w -\frac{kwqB_{T}e^{mt}}{(B_{T}+e^{mt}p)^{2}} \\ &\quad + \beta zqe^{-a_{0}e^{mt}x}e^{mt}-\frac{1}{2}\beta z e^{-a_{0}e^{mt}x}e^{mt}qu_{RTI}-yp_{0}x^{2}e^{mt}q \end{align*} and \begin{align*} &-\dot {\bar w} +m\bar w\\ &=e^{mt} + \frac{r\bar w\bar qe^{mt}}{B_{V}+e^{mt}\bar q} -\beta \bar w\bar qe^{-a_{0}e^{mt}\bar x}e^{mt}+\frac{1}{2}\beta \bar we^{-a_{0}e^{mt}\bar x}e^{mt}\bar q\bar u^{*}_{RTI}-\mu_{T}\bar w \\ &\quad -\frac{kB_{T}e^{mt}\bar w\bar q}{(B_{T}+e^{mt}\bar p)^{2}} + \beta \bar z\bar qe^{-a_{0}e^{mt}\bar x}e^{mt} -\frac{1}{2}\beta \bar ze^{-a_{0}e^{mt}\bar x}e^{mt}\bar q\bar u_{RTI} -p_{0}\bar y\bar x^{2}e^{mt}\bar q \end{align*} respectively. Subtracting the expression for $\bar \lambda_{1}$ from the expression for $\lambda_{1}$ and multiplying by $(w-\bar w)$ and integrating from $t=0$ to $t=T_{f}$ we have \begin{align*} &\frac{1}{2}(w-\bar w)^{2}(0)+m\int_{0}^{T_{f}}(w-\bar w)^{2}dt \\ &= r\int_{0}^{T_{f}}e^{mt}\left(\frac{wq}{B_{V}+e^{mt}q} -\frac{\bar w\bar q}{B_{V}+e^{mt}\bar q}\right)(w-\bar w)dt \\ &\quad-\beta \int_{0}^{T_{f}}e^{mt}\left(e^{-a_{0}e^{mt}x}wq-e^{-a_{0}e^{mt} \bar x}\bar w\bar q\right)(w-\bar w)dt \\ &+ \frac{\beta}{2}\int_{0}^{T_{f}}e^{mt}\left(e^{-a_{0}e^{mt}x}wqu^{*}_{RTI} -e^{-a_{0}e^{mt}\bar x}\bar w\bar q\bar u^{*}_{RTI}\right)(w-\bar w)dt \\ &\quad +\beta \int_{0}^{T_{f}}e^{mt}\left(e^{-a_{0}e^{mt}x}zq-e^{-a_{0}e^{mt} \bar x}\bar z\bar q\right)(w-\bar w)dt-\mu_{T}\int_{0}^{T_{f}}(w-\bar w)^{2}dt\\ &\quad -\frac{\beta}{2} \int_{0}^{T_{f}}e^{mt}\left(e^{-a_{0} e^{mt}x}zqu^{*}_{RTI}-e^{-a_{0}e^{mt}\bar x}\bar z\bar q\bar u^{*}_{RTI}\right) (w-\bar w)dt \\ &\quad -p_{0}\int_{0}^{T_{f}}e^{mt}(yxq-\bar y\bar x\bar q)(w-\bar w)dt \\ &\quad -kB_{T}\int_{0}^{T_{f}}e^{mt}\Big(\frac{wq}{(B_{T}+e^{mt}p)^{2}} -\frac{\bar w \bar q}{(B_{T}+e^{mt}\bar p)^{2}}\Big)(w-\bar w)^{2}dt. \end{align*} Similarly, the equations for $T^{*}$ and $\bar T^{*}$, $V$ and $\bar V$, $C$ and $\bar C$, $\lambda_{2}$ and $\bar \lambda_{2}$, $\lambda_{3}$ and $\bar \lambda_{3}$, $\lambda_{4}$ and $\bar \lambda_{4}$ are subtracted, then each expression is multiplied by an appropriate function and integrated from $t=0$ to $t=T_{f}$. We obtain eight integral equations and we use estimates to obtain the result. Several terms are estimated in these eight equations. For example the third term on the right-hand side of equation \ref{e17}, \begin{align*} &k\int_{0}^{T_{f}}e^{mt}\Big(\frac{pq}{B_{T}+e^{mt}p} -\frac{\bar p\bar q}{B_{T}+e^{mt}\bar p}\Big)(p-\bar p)dt\\ &\leq C_{1}e^{mt}\int_{0}^{T_{f}}((p-\bar p)^{2}+(q-\bar q)^{2})dt, \end{align*} utilizes upper bounds on the solutions. Other estimates can be presented by utilizing upper bounds on solutions. They involve separating terms that involve squares, powers, several multiplied terms, and quotients. Also using Theorem 6.1 we have \begin{align*} &|u^{*}_{RTI}(t)-\bar u^{*}_{RTI}(t)|\\ &\leq \frac{\beta}{2A_{1}}\left|e^{-a_{0}e^{mt}x}pq(w-z) -e^{-a_{0}e^{mt}\bar x}\bar p\bar q(\bar w-\bar z)\right| \\ &\leq \frac{\beta}{2A_{1}}\left|\big(e^{-a_{0}e^{mt}x}pqw-e^{-a_{0}e^{mt} \bar x}\bar p\bar q\bar w\big) -\big(e^{-a_{0}e^{mt}x}pqz-e^{-a_{0}e^{mt}\bar x}\bar p\bar q\bar z\big)\right| \end{align*} and $$ |u^{*}_{PI}(t)-\bar u^{*}_{PI}(t)|\leq \frac{\alpha N}{2A_{2}} \left|e^{-a_{1}e^{mt}x}vp^{*}-e^{-a_{1}e^{mt}\bar x}\bar v\bar p^{*}\right|. $$ To show uniqueness, the integral equations are combined. Adding all the eight estimates gives \begin{align*} &\frac{1}{2}(p-\bar p)^{2}(T_{f})+\frac{1}{2}(p^{*}-\bar p^{*})^{2}(T_{f}) +\frac{1}{2}(q-\bar q)^{2}(T_{f})+\frac{1}{2}(x-\bar x)^{2}(T_{f})\\ &+\frac{1}{2}(w-\bar w)^{2}(0)+\frac{1}{2}(z-\bar z)^{2}(0) +\frac{1}{2}(v-\bar v)^{2}(0)+\frac{1}{2}(y-\bar y)^{2}(0) \\ &+ m\int_{0}^{T_{f}}[(p-\bar p)^{2}+(p^{*}-\bar p^{*})^{2}+(q-\bar q)^{2} +(x-\bar x)^{2}+(w-\bar w)^{2}\\ &+(z-\bar z)^{2}+(v-\bar v)^{2} +(y-\bar y)^{2}]dt \\ &\leq (\tilde{C_{1}}+\tilde{C_{2}}e^{3mT_{f}})\int_{0}^{T_{f}}[(p-\bar p)^{2} +(p^{*}-\bar p^{*})^{2}+(q-\bar q)^{2}+(x-\bar x)^{2}]dt \\ &\quad +\int_{0}^{T_{f}}[(w-\bar w)^{2}+(z-\bar z)^{2}+(v-\bar v)^{2} +(y-\bar y)^{2}]dt. \end{align*} Thus from the above expression, using the non-negativity of the variable expressions evaluated at the initial and the final time and simplifying, the inequality is reduced to \begin{align*} &(m-\tilde{C_{1}}-\tilde{C_{2}}e^{3mT_{f}})\int_{0}^{T_{f}}[(p-\bar p)^{2} +(p^{*}-\bar p^{*})^{2}+(q-\bar q)^{2}+(x-\bar x)^{2}+(w-\bar w)^{2}\\ &+(z-\bar z)^{2}+(v-\bar v)^{2}+(y-\bar y)^{2}]dt \leq 0 \end{align*} where $\tilde{C_{1}}$ and $\tilde{C_{2}}$ depend on the coefficients \cite{f1}; \cite{j1} and the bounds on all solution variables $p, p^{*}, q, x, w, z, v, y$. If we choose $m$ such that $m -\tilde{C_{1}}-\tilde{C_{2}}e^{3mT_{f}}>0$, the above inequality holds if the integrand is identically zero. Since the natural logarithm is an increasing function, then $\ln\big(\frac{m-\tilde{C_{1}}}{\tilde{C_{2}}}\big)>3mT_{f}$ if $m > \tilde{C_{1}}+\tilde{C_{2}}$. This gives that $T_{f} < \frac{1}{3m}\ln\big(\frac{m-\tilde{C_{1}}}{\tilde{C_{2}}}\big)$, then $p=\bar p$, $p^{*}=\bar p^{*}$, $q=\bar q$, $x=\bar x$, $w=\bar w$, $z=\bar z$, $v=\bar v$, $y=\bar y$. Hence the solution is unique for small time. \end{proof} \section{Numerical Simulations} The optimality system in section 5 is solved using an iterative method with Runge-Kutta of order four scheme. The optimality system is a two-point boundary value problem, where initial conditions are specified for the state system and terminal conditions are specified for the adjoint system. The method of obtaining the optimal control is as follows \cite{j1}: \begin{enumerate} \item Take a guess for the two controls. \item Solve the state system forward using those controls and using a Runge-kutta method of order four algorithm. Use state variables initial conditions. \item Using the new state values, solve the adjoint system backwards using the final time zero boundary conditions and Runge-Kutta of order four scheme. \item Calculate the new control values from the characterization. \item Go to steps 2, 3 again with new control from step 4. \item Calculate other new control values from step 5. Compare controls from last iteration to new iteration and compare states also. Keep repeating control updates and forward and backward solving until the iterates converge. \end{enumerate} In the virtual simulations in this section we chose a set of parameters and initial values yielding approximately realistic population numbers. Our initial conditions resemble clinical data for HIV infected individuals during symptomatic phase. Since therapy is initiated when patients are symptomatic, we consider cases when CD4+ T cell count is less than or equal 250 cells $\mu l^{-1}$. The following parameter values have been used to generate the solutions in this section: $s_{1}=20$, $\mu_{T}=0.02$, $r=0.01$, $N=1000$ \cite{k1}, $s_{2}=10$, $\beta=2\times 10^{-4}$, $B_{T}=350$, $B_{V}=400$, $h=2\times 10^{-3}$, $\mu_{C}=1.5$, $\mu_{V}=0.95$ \cite{b1}, $k=2\times 10^{-3}$, $a_{0}=0.01$, $a_{1}=0.075$, $p_{0}=1\times 10^{-5}$ and $\alpha=0.25$ \cite{p1}. The bounds for controls are $a_{11}=0.0$, $a_{22}=0.0$, $b_{11}=0.002$ and $b_{22}=0.9$. Figure \ref{fig:figure1} shows the graph of the solution to the optimality system, showing propagation of CD4+ T cells, infected CD4+ T cells, reverse transcriptase inhibitor and a protease inhibitor when treatment is administered for 50 days. Here we initiate treatment when CD4+ T cell count, T(0) is 250 cells $\mu l^{-1}$, viral load, V(0) is 4000 copies $ml^{-1}$, infected CD4+ T cells, $T^{*}(0)$ of 200 cells $\mu l^{-1}$ and a CTL count, C(0) of 100 $\mu l^{-1}$. The value for the first weight factor is given by $A_{1}=250 000$ and the second weight factor $A_{2}=250$ \cite{j1}. Figure \ref{fig:figure1}(a) shows the propagation of uninfected CD4+ T cells after initiation of therapy. Initial decline of CD4+ T cells for a day is due to pharmacokinetics and pharmacology delay, and thereafter T cell count start to increase for about 25 days, nearly an increase of 100\% and thereafter CD4+ T cells start to gradually decline. Within the first week of drug administration, viral load drops to zero, figure \ref{fig:figure1}(c), followed by a sudden increase and oscillation at around 400 for another week and then stabilizes for the next 25 days and starts to slowly increase at day 40. The CD4+ T cell and viral kinetics are produced if the nucleoside RTI (figure \ref{fig:figure1}(b)) is administered in full scale for one day (normalised dosage size of 0.002) after 4 days, and the drug is stopped for 2-3 days. Drug administration resume with almost 10\% of the initial dosage size and gradually increased up to day 30 and then tapered off up to day 50. Therapy start with high doses (normalised dosage size of 0.9) of a protease inhibitor (figure \ref{fig:figure1}(d)) for the first day, then stopped for 5 days. The drug is administered at full scale (dosage size of 0.9) for a day and then stopped for one day or two days. Finally, 40\% of the initial dosage is administered for 1-2 days and stopped for 30-35 days and resumed with 2\%-5\% of the initial therapy for the last 3 days. The effect of the regimen, in a short term, managed to increase the CD4+ T cells to nearly 450 cells in 25 days and level of viremia is suppressed to low levels (below 500 copies of RNA $ml^{-1}$) which is beneficial to the infected individual's health. \begin{figure}[ht] \begin{center} \includegraphics[width=0.48\textwidth]{tcell1} \includegraphics[width=0.48\textwidth]{urti1}\\ (a) \hspace{6cm}(b)\\ \includegraphics[width=0.48\textwidth]{viral1} \includegraphics[width=0.48\textwidth]{upi1} \\ (c) \hspace{6cm}(d) \end{center} \caption{Graph of the numerical solution to the optimality system, showing propagation of CD4+ T cells, infected CD4+ T cells, reverse transcriptase inhibitor and a protease inhibitor when treatment is administered for 50 days. Here we initiate treatment when CD4+ T cell count, T(0) is 250 cells $\mu l^{-1}$, viral load, V(0) is 4000 copies $ml^{-1}$, infected CD4+ T cells, $T^{*}(0)$ of 200 cells $\mu l^{-1}$ and a CTL count, C(0) of 100 $\mu l^{-1}$. The value for the first weight factor is given by $A_{1}=250 000$ and the second weight factor $A_{2}=250$. (a) CD4+ T cell kinetics (b) Reverse Transcriptase Inhibitor dosage sizes (c) Viral Load (d) Protease Inhibitor dosage sizes.} \label{fig:figure1} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.48\textwidth]{tcell2} \includegraphics[width=0.48\textwidth]{urti2}\\ (a) \hspace{6cm}(b)\\ \includegraphics[width=0.48\textwidth]{viral2} \includegraphics[width=0.48\textwidth]{upi2} \\ (c) \hspace{6cm}(d) \end{center} \caption{Graph of the numerical solution to the optimality system, showing propagation of CD4+ T cells, infected CD4+ T cells, reverse transcriptase inhibitor and a protease inhibitor when treatment is administered for 50 days. Here we initiate treatment when CD4+ T cell count, T(0) is 250 cells $\mu l^{-1}$, viral load, V(0) is 4000 copies $ml^{-1}$, infected CD4+ T cells, $T^{*}(0)$ of 200 cells $\mu l^{-1}$ and a CTL count, C(0) of 100 $\mu l^{-1}$. The value for the first weight factor is given by $A_{1}=250 000$ and the second weight factor $A_{2}=250$. The activity CTL in reducing burst size and reducing infection has been decreased to $a_{1}=0.075$ and $a_{0}=0.001$ respectively. Here $a_{0}$ and $a_{1}$ are different from figure \ref{fig:figure1}. (a) CD4+ T cell kinetics (b) Reverse Transcriptase Inhibitor dosage sizes (c) Viral Load (d) Protease Inhibitor dosage sizes.} \label{fig:figure2} \end{figure} Figure \ref{fig:figure2} shows the graph of the solution to the optimality system, showing propagation of CD4+ T cells, infected CD4+ T cells, reverse transcriptase inhibitor and a protease inhibitor where CTL activity is decreased. Here we initiate treatment when CD4+ T cell count, T(0) is 250 cells $\mu l^{-1}$, viral load, V(0) is 4000 copies $ml^{-1}$, infected CD4+ T cells, $T^{*}(0)$ of 200 cells $\mu l^{-1}$ and a CTL count, C(0) of 100 $\mu l^{-1}$, and treatment is administered for 50 days. The value for the first weight factor is given by $A_{1}=250 000$ and the second weight factor $A_{2}=250$. A decrease in the effect of CTL activity, that is decrease in $a_{0}=0.001$ and $a_{1}=0.075$ with the same initial conditions leads to CD4+ T cells decline for 1-2 days due to pharmacology and pharmacokinetics (similar to Dixit and Perelson \cite{d1} results) just after initiation of therapy (figure \ref{fig:figure2}(a)) and build up for nearly 10 days before a decreasing tendency up to day 50. Viral load sharply decreases to very low levels (figure \ref{fig:figure2}(c)) during the first 5 days of therapy initiation, sharply increase and sharply decrease before it starts to gradually increase. The kinetics of the CD4+ T cells and viral load are given if the nucleoside reverse transcriptase inhibitor (normalised dosage size of 0.002) is administered for 1 day and then stopped for 4 days (figure \ref{fig:figure2}(b)). The drug is given in small dosage sizes, increased daily to almost 100\% of initial dosage at day 10 for 2 days. The dosage size is decreased to nearly 40\%, then increased at day 15 and then tapered off up to day 50. Protease inhibitor drug schedule is started at day 6 (figure \ref{fig:figure2}(d)), increased to nearly 0.9 at day 9 and lowered to nearly zero at day 12. Small quantities are given, nearly 10\% of the initial dosage size and maintained up to day 50. The scheme is effective during the first 10 days and thereafter the immune system is heavily compromised due to little activity of CTLs. Figure \ref{fig:figure3} shows the graph of the solution to the optimality system, showing propagation of CD4+ T cells, infected CD4+ T cells, reverse transcriptase inhibitor and a protease inhibitor when treatment is administered for 50 days. Here we initiate treatment when CD4+ T cell count T(0) is 250 cells $\mu l^{-1}$, viral load V(0) is 4000 copies $ml^{-1}$, infected CD4+ T cells, $T^{*}(0)$ is 200 cells $\mu l^{-1}$ and a CTL count, C(0) of 100 $\mu l^{-1}$. The value for the first weight factor is given by $A_{1}=250$ and the second weight factor $A_{2}=100$. Exploring the effects of drug toxicity, that is reducing the weight factors $A_{1}=250$ and $A_{2}=100$ and all the other parameters remain as in figure (\ref{fig:figure2}). We have the following numerical results, given by figure (\ref{fig:figure3}): Effect of reducing the weight factors simulate the effect of a decrease in drug toxicity, and therefore we observe that dosage sizes for both drug types are increased. Since the efficacy of drugs is not changed, CD4+ T cells and viral kinetics (figure \ref{fig:figure3}(a) and (c)) respectively) do not change. The nucleoside reverse transcriptase inhibitor (figure \ref{fig:figure3}(b)) is administered at full scale for 2 days and stopped for 1 to 2 days. Therapy is given again at full scale up to day 50. Protease inhibitor schedule (figure \ref{fig:figure3}(d)) is given at day 7, increased to full scale for 4-5 days and decreased to nearly 11\% up to day 50. \begin{figure}[ht] \begin{center} \includegraphics[width=0.48\textwidth]{tcell3} \includegraphics[width=0.48\textwidth]{urti3}\\ (a) \hspace{6cm}(b)\\ \includegraphics[width=0.48\textwidth]{viral3} \includegraphics[width=0.48\textwidth]{upi3} \\ (c) \hspace{6cm}(d) \end{center} \caption{Graph of the numerical solution to the optimality system, showing propagation of CD4+ T cells, infected CD4+ T cells, reverse transcriptase inhibitor and a protease inhibitor when treatment is administered for 50 days. Here we initiate treatment when CD4+ T cell count, T(0) is 250 cells $\mu l^{-1}$, viral load, $V(0)$ is 4000 copies $ml^{-1}$, infected CD4+ T cells, $T^{*}(0)$ of 200 cells $\mu l^{-1}$ and a CTL count, C(0) of 100 $\mu l^{-1}$. The value for the first weight factor is given by $A_{1}=250$ and the second weight factor $A_{2}=100$. (a) CD4+ T cell kinetics (b) Reverse Transcriptase Inhibitor dosage sizes (c) Viral Load (d) Protease Inhibitor dosage sizes.} \label{fig:figure3} \end{figure} \section{Discussion} The virtual combined therapy simulations in this paper are designed to provide insights of drug scheduling in short term therapy performance. Similar to Bajaria {\it et al.} \cite{b1}, reduced dosage sizes and drug holidays can achieve goals of antiretroviral therapy (increase of CD4+ T cells and suppression of viral load). In poor resourced settings, an effective schedule should improve the patient's life, be affordable (small dosages sizes) and allowing drug holidays should relax the concept of strict patient adherence to treatment and compliance. The scheme should increase or restore the immunological function (CD4+ T cell increases) through less exposure to drugs. We also observe that the weight factors in the objective functional model drug toxicity, therefore the more toxic the drugs are the less dosage sizes to be administered to reduce systemic costs. Figure \ref{fig:figure1} shows that if an individual's immune response is strong, virus can be effectively controlled during therapy whilst weak immune responses, figure \ref{fig:figure2} and figure \ref{fig:figure3}, lead to a short term control of the virus. We can conclude that effectiveness of therapy largely depends on the immune response of an individual, that is if the immune response is better, virus can be controlled effectively. Also an effective immune response leads to administration of small dosage sizes which are cost effective. We observe from numerical illustrations that if therapy can induce CTL activity, control of viremia is feasible and this can facilitate the implementation of drug holidays. During drug holiday periods, CTLs will be controlling viremia. Due to drug toxicity, allowing drug holidays can be beneficial in the short term implementation of HAART. If therapy has less toxic effects, continuous therapy is beneficial as there will be less harmful side effects. Our dynamical model did not take into account the effects of viral population mutation over time in response to drug therapy. These effects can become significant in the case of long term anti-HIV therapy. \subsection*{Acknowledgments:} The authors would like to acknowledge with thanks financial support from Eagle Insurance Company Limited, Zimbabwe. T. Shiri would like to thank National University of Science and Technology for the SDF scholarship. We thank Professor Suzanne Lenhart for the optimal control numerical algorithm and also Professor Les Jennings for critical and helpful comments on the numerical illustrations. \begin{thebibliography}{00} \bibitem{a1}Ahr, B., Robert-Hebmann, V., Devaux, C. and Biard-Piechaczyk, M.; Apoptosis of uninfected cells induced by HIV envelope glycoproteins. {\it Retrovirology.} (2004), 1:12. \bibitem{a2}Altes, H. K., Price, D. A., Jansen, V. A. A.; Effector cytotoxic T lymphocyte numbers induced by vaccination should exceed levels in chronic infection for protection from HIV. {\it Science.} {\bf 20} (2001), 3-6. \bibitem{b1}Bajaria, H. S., Webb, G., Kirschner, E.D.; Predicting differential responses to structured treatment interruptions during HAART. {\it Bulletin of Mathematical Biology.} (2003), 1-26. \bibitem{b2}Butler, S., Kirschner, D., Lenhart, S.; Optimal control of chemotherapy affecting the infectivity of HIV. {\it In: Arino O, Axelrod D, Kimmel M, Langlais M (eds): Advances in Mathematical Population Dynamics: Molecules, Cells, Man. World Scientific Publishing,} (1997) 104-120. \bibitem{d1}Dixit, M. N., Perelson, S. A.; Complex patterns of viral load decay under antiretroviral therapy: influence of pharmacokinetics and intracellular delay. {\it Journal of Theoretical Biology.} {\bf 226} (2004), 95-109. \bibitem{f1}Fister, R., Lenhart S., McNally, J. S.; Optimizing chemotherapy in an HIV model. {\it J.of Diff Eqn.} {\bf 32} (1998), 1-12. \bibitem{g1}Gallant, E. J.; Antiretroviral therapy update from the 44th ICAAC. {\it The Johns Hopkins University AIDS Service.} {\bf 17} (2005), 1-8. \bibitem{j1}Joshi, H. R.; Optimal Control of an HIV Immunology Model. {\it Optim. Contr.Appl. Math} {\bf 23} (2002), 199-213. \bibitem{k1}Kirschner, D.; Using mathematics to understand HIV immune dynamics. {\it Notices Amer Math Soc.} {\bf 43} (1996), 191-202. \bibitem{k2}Kirschner, D., Lenhart, S., Serbin, S.; Optimal control of the chemotherapy of HIV. {\it Journal of Mathematical Biology.} {\bf 35} (1997), 775-792. \bibitem{k3}Kirschner, E. D., Webb, F. G.; Immunotherapy of HIV-1 infection. {\it Journal of Biological Systems.} {\bf Vol. 6, No. 1} (1998), 71-83. \bibitem{k4}Klenerman, P., Wu, Y., Philips, R.; HIV: current opinion on escapology. {\it Current Opinion in Microbiology.} {\bf 5} (2002), 408-413. \bibitem{k5}Kot, M.; {\it Elements of Mathematical Ecology.} Cambridge University Press, 2001. \bibitem{k6}Krakauer, D. and Nowak, M. A.; T cell induced pathogenesis in HIV:bystander effects and latent infection. {\it Proc. Royal Soc. London B.} {\bf 266} (1999), 1069-1075. \bibitem{k7}Kutch, J.J., Gufil, P.; Optimal control of HIV infection with a continuously mutating viral population. {\it Proceedings of the American Control Conference}, (2002), 4033-4038. \bibitem{p1}Perelson, S. A.; Modelling viral and immune system dynamics. {\it MacMillan Magazines Ltd.} {\bf 2} (2002), 28-36. \bibitem{p2}Perelson, S. A., Kirschner, E. D., De Boer, R.; Dynamics of HIV-infection of CD4+ T cells. {\it Mathematical Biosciences.} {\bf 114:} (1993), 81-125. \bibitem{s1}Sewell, A. K., Price, D. A., Oxenius, A., Kelleher, A. D., Phillips, R. E.; Cytotoxic T lymphocyte responses to human immunodeficiency virus: control and escape. {\it Stem Cells.} {\bf 18} (2000), 230-244. \bibitem{s2}Shim, H., Han, C. C., Chung, S.W., Nam, W. S., Seo, H. J.; Optimal scheduling of drug treatment for HIV infection: continuous dose control and receding horizon control. {\it International Journal of Control, Automation, and Systems.} {\bf 1,} (2003), 401-407. \bibitem{w1}Wodarz, D., Nowak, A. M.; Immune responses and viral phenotype: do replication rate and cytopathogenecity influence virus load? {\it Journal of Theoretical Medicine.} {\bf2} (2000), 113-127. \end{thebibliography} \end{document}