\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 $a**0$ 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}
**