\documentclass[letterpaper,10pt,conference]{cssconf}
\IEEEoverridecommandlockouts
\overrideIEEEmargins
\usepackage{amsmath}
%\documentclass[twocolumn]{article}
%\usepackage{verbatim}
%\usepackage{amsfonts}
\usepackage{graphicx}
%\usepackage{amsmath}
%\setlength{\textheight}{9.0in}
%\setlength{\textwidth}{6.5in}
%\voffset=-1.200in
%\renewcommand{\baselinestretch}{1.00}
\def\eq{\begin{eqnarray}}
\def\nq{\end{eqnarray}}
\def\equnn{\begin{eqnarray*}}
\def\nqunn{\end{eqnarray*}}
\def\ebf{\mathbf{e}}
\def\fbf{\mathbf{f}}
\def\Dbf{\mathbf{D}}
\def\Fbf{\mathbf{F}}
\def\Gbf{\mathbf{G}}
\def\jbf{\mathbf{j}}
\def\Nbf{\mathbf{N}}
\def\Ubf{\mathbf{U}}
\def\xbf{\mathbf{x}}
\def\ybf{\mathbf{y}}
\def\Ybf{\mathbf{Y}}
\def\Zbf{\mathbf{Z}}
\def\zerobf{\mathbf{0}}
\def\xibf{\mbox{\boldmath{$\xi$}}}
\def\etabf{\mbox{\boldmath{$\eta$}}}
\def\chibf{\mbox{\boldmath{$\chi$}}}
\def\Gambf{\mbox{\boldmath{$\Gamma$}}}
\def\Nbfhat{\widehat{\mathbf{N}}}
\def\Hc{\mathcal{H}}
\def\d{\delta}
\def\D{\Delta}
\def\gradx{\nabla_{x}}
\def\grady{\nabla_{Y}}
\def\And{ ~ \mathrm{and} ~ }
\def\For{ ~ \mathrm{for} ~ }
\def\On{ ~ \mathrm{on} ~ }
\def\Or{ ~ \mathrm{or} ~ }
\def\tol{\mathrm{tol}}
\def\half{\frac{1}{2}}
\def\DDot{\!\!:\!\!}
%
% Title *****
\title{\textbf Optimal control of drug delivery to brain tumors for a
distributed parameters model
\thanks{Work supported in part by the National Science Foundation
Grants DMS--99--73231 and DMS--02--07081. The first author was also
supported by the \textit{Provost's award for Graduate Research}
awarded by the Graduate College at the University of Illinois at
Chicago.  Submitted to \textit{Proceedings of American Control Conference
2005}.}}
%
% Authors *****
\author{
\parbox[c]{2.6in}{\centering Siddhartha P. Chakrabarty \\
{\small
Department of Mathematics, Statistics and Computer Science\\
University of Illinois at Chicago\\
851 S Morgan St.; m/c 249\\
Chicago, IL  60607--7045, USA
\\*[0.25ex]
E--Mail:  pratim@math.uic.edu
\\*[0.25ex]
URL: http://www.math.uic.edu/{$\sim$}pratim/}}~~ and ~~
\parbox[c]{2.6in}{\centering Floyd B. Hanson\\
{\small
Laboratory for Advanced Computing\\
University of Illinois at Chicago\\
851 S Morgan St.; m/c 249\\
Chicago, IL  60607--7045, USA
\\*[0.25ex]
E--Mail:  hanson@math.uic.edu
\\*[0.25ex]
URL: http://www.math.uic.edu/{$\sim$}hanson/\\}}}
%
% Document Begins *****
\begin{document}

\maketitle

% Abstract *****
\begin{abstract}
The growth and treatment of brain tumors is mathematically examined
using a distributed parameters model. The model is a system of three
coupled reaction diffusion equations involving the tumor cells,
normal tissue and the drug concentration. An optimal control problem
is designed, with the drug delivery rate as the control and solved
to obtain the state and co-state equations as well as the regular
control using a modified double shot forward-backward method. 
This gives rise to a
coupled system of equations with a forward state equation and a
backward co-state equation, which is solved using a double shot method.
%FBH 9/14:
A numerical procedure based upon the Crank-Nicolson method is
used to solve the coupled system of two three-dimension partial
differential equations.
%%
\end{abstract}
~\\*[-2.5em]

% Introduction *****
\section{Introduction}
%
The growth and control of brain tumors have been the subject of
medical and scientific scrutiny for a very long time. Simply
speaking a tumor, like most cancerous cells originate from a cell,
that proliferates and effects its neighboring normal tissues. As the
tumor cells become malignant they become more dangerous for the
host. Understanding the mechanism of tumor progression is necessary
for its diagnosis and treatment. As mentioned earlier, brain tumors
have been studied for centuries, dating back to as early as 2500 BC
\cite{Murray2}. The most common and deadly form of brain tumor are the
gliomas, which account for more than half of the brain tumor cases.
Gliomas are highly invasive and severely infiltrate the surrounding
tissues \cite{KRS99}. Despite improved diagnostic procedures such as
computerized tomography (CT) scan and magnetic resonance imaging
(MRI), the benefits of such modern accessories have been restricted
by the treatment options available. One major problem of
administering the drugs to the brain tumor site is the \textit{blood
brain barrier (BBB)} \cite {HHE00}, which exists in the human brain
as a protection for the brain cells and as a restriction on the
transport of water soluble substances between the blood and the
central nervous system. Another problem that arises is the resection
%FBH 9/14:
of a tumor, after the core mass of the tumor has been surgically
removed. One way to deliver the drug is to use the drug conjugated
with a \textit{non-toxic} polymer like the BCNU and place it in the
brain cavity for controlled release. A more traditional way would be
%FBH 9/12:
to deliver the drug by an optimal distribution
of the drug about the original tumor site. 
%
Wang et al. \cite{CHW98,CHW99} have studied the drug delivery
method to tumors in three dimension for drugs like IgG and BCNU.

%FBH 9/12:
The paper will mostly focus on the control for the  optimal distribution
of the drug about the original tumor site.  
%
Before that we need to look at the mechanisms behind
the growth of tumor cells as well as normal tissues and the drug
concentration in the tumor. Unlike a lot of other tumors, gliomas
can be highly diffusive \cite{Murray2}. Gatenby et al. \cite{RAG96}
and Mansuri \cite{JM02} study the mechanism of reaction
diffusion in the growth of tumors. They also take into account the
effects of competition for resources between the cancerous cells and
the healthy tissues. Westman et al. \cite{JJW02} look at the various
types of tumor growth, namely exponential, logistic and Gompertz.
%FBH 9/14:
Murray's books \cite{Murray1,Murray2} are excellent references for the study
of different types of growth mechanisms. Also, Woodward et al.
\cite{DEW96} study a model of glioma growth and the effects of
surgical resection. In this paper we set up a fairly generalized
% FBH 9/14:
distributed parameters model for the PDE driven system, 
define an objective functional to minimize drug delivery and tumor
burden costs, and use
a modified \textit{Lagrange multiplier method} \cite{MDG03} for 
including constraints in the
%%
control problem. Finally a \textit{double shot}, forward and backward, 
numerical method is
given to solve for the state, co-state and regular control. To solve
for the \textit{PDEs}  we use a Crank-Nicolson predictor-corrector
method developed by Hanson et al. \cite{FBH93} and Hanson
\cite{FBH96}, for stochastic dynamic programming.

% Mathematical Model *****
\section{Mathematical Model}
% FBH 9/10: must be density, not number! for spatial case.
Let $Y_{1}=n_{1}(\xbf,t)$ be the density of tumor cells, 
$Y_{2}=n_{2}(\xbf,t)$ be the density of  normal tissue and
$Y_{3}=c(\xbf,t)$ be the drug concentration at any vector position $\xbf$
and time $t$. This spatio-temporal model is a system of three coupled
reaction-diffusion equations.

\subsection{Tumor Cells}
%
It is assumed that the density of tumor cells, $n_1=n_1(\xbf,t)$,
satisfy a reaction-diffusion equation subject to competition with
the normal cells, $n_2=n_2(\xbf,t)$,
%
\eq\label{n1eq}
\frac{\partial n_1}{\partial t}\!\!=\!\!D_1\gradx^2[n_1]\!+\!a_1 f_1(n_1)\!-\!%
\left(\alpha_{1,2}n_2\!+\!b_1(c)\right)n_1,
\nq
%
where the tumor
diffusivity  is $D_1$. Let the term $a_{1}f_{1}(n_{1})$ be the
growth rate of the tumor cells, where $f_{1}$ could be exponential,
logistic or Gompertz growth ,
\[f_{1}(n_{1})=n_{1},~~n_{1}\left(1-\frac{n_{1}}{k_{1}}\right)~~
\textrm{or}~~n_{1}\ln \left(\frac{k_{1}}{n_{1}}\right),\]
respectively, where $a_1$ is the tumor cell intrinsic growth rate
and $k_1$ is the tumor cell carrying capacity. Let $\alpha_{1,2}$
denote the death rate of the tumor cells due to competition for
resources with the normal tissue. Let $b_1(c)$ be the death rate of
tumor cells due to drug treatment, which could be a function of the
localized drug concentration $c(\xbf,t)$ at the tumor site, e.g.,
purely linear function,
% FBH 9/10 reassigned kappas
% FBH 9/11-13 deleted ramp version to correct for page 3 avalanche and symmetry:
\equnn
%\begin{array}{c}
b_1(c(\xbf,t))\!=\!\kappa_{1,3}\!\cdot\!c(\xbf,t)  
%~  \Or 
%\\
%b_1(c(\xbf,t))\!=\!\kappa_{1,3}\!\cdot\!%
%\left(c(\xbf,t)\!-\!c_{\min}\right)\widetilde{\Hc}
%\left(c(\xbf,t)\!-\!c_{\min}\right), 
%\end{array}
\nqunn
%
respectively.  
%Here, $c_{\min}$ is
%the minimum drug concentration required to have effective treatment
%and $\widetilde{\Hc}$ is the Heavyside function.  In the current application,
%the pure linear form will be used.

\subsection{Normal Tissue}
%
Similar assumptions are made for the density of normal cells $ n_2 =
n_2(\xbf,t)$ with similar coefficients. Thus, the reaction-diffusion
equation for normal tissue evolution is as follows,
%
\eq\label{n2eq}
\frac{\partial n_{2}}{\partial t}\!\!=\!\!D_2\gradx^{2}[n_{2}]%
\!+\!a_{2}f_{2}(n_{2})\!-\!\left(\alpha_{2,1}n_{1}\!+\!b_{2}(c)\right)n_{2},
\nq
%
where $a_2$ is the normal cell intrinsic growth rate and the normal
cell growth function $f_{2}$ is either exponential, logistic or
Gompertz growth ,
\[f_{2}(n_{2})=n_{2},~~n_{2}\left(1-\frac{n_{2}}{k_{2}}\right)~~
\textrm{or}~~n_{2} \ln \left(\frac{k_{2}}{n_{2}}\right),\]
respectively, where $k_2$ is the normal tissue carrying capacity.
Note that the $b_{2}(c)$ term indicates that some normal tissues
could die as a result of the treatment and $b_{2}(c)$ must also
depend on the local drug concentration  $c(\xbf,t)$, e.g.,
$$b_2(c(\xbf,t))=\kappa_{2,3}\cdot c(\xbf,t).$$ % FBH 9/10 added

\subsection{Concentration}
%
It is assumed that the drug shows a diffusive behavior and that
there is a reabsorption at the rate $a_{3}$. Also let $u=u(\xbf,t)$ be
the rate at which the drug is being injected. The choice of the
letter $u$ indicates that we will use it as the control when dealing
with the optimal control system. The equation for drug concentration
at position $\xbf$ and time $t$ is,
% FBH 9/11: -u->+u
\eq\label{ceq}
\frac{\partial c}{\partial t}\!\!=\!\!D_3\gradx^{2}[c]\!+\!a_{3}f_{3}(c)\!+\!u,
% FBH 9/10 u0 belongs only in J:+\!(u_{0}\!-\!u(\xbf,t)),
\nq
%
where $f_{3}(c) = - c$ is the reabsorption function\cite{CHW99}.
% FBH 9/10 u0 belongs only in J: and $u_{0}$ is some drug rate reference value.

\subsection{Global State Vector}
%
Let the global state vector be
%
\eq\label{Ydefeq}
\Ybf(\xbf,t) = [Y_{i}(\xbf,t)]_{3\times{1}}
= \begin{bmatrix}n_{1}(\xbf,t)\\n_{2}(\xbf,t)\\c(\xbf,t)
\end{bmatrix},
\nq
%
at position $\xbf$  in the \textit{interior} $\Omega$ of the state domain and time $t$
on $[0,t_{f}]$.

\subsection{Initial and Boundary Conditions}
%
Let the initial conditions for the state be
%
\eq\label{Yic}
\Ybf(\xbf,0)=\begin{bmatrix}n_{1}(\xbf,0)\\n_{2}(\xbf,0)\\c(\xbf,0)
\end{bmatrix}=
\begin{bmatrix}{n_{1}}_{0}(\xbf)\\{n_{2}}_{0}(\xbf)\\{c}_{0}(\xbf)
\end{bmatrix} \equiv \Ybf_0(\xbf) ,
\nq
%
for $\xbf$ in $\Omega$.   
% FBH 9/11:
Murray  \cite{Murray2} recommends using Gaussian
distribution for the initial distributions of tumors.
% 
The \textit{no flux} boundary conditions are
%
\eq\label{Ybc}
-D(\Nbfhat\cdot\gradx)[\Ybf](\xbf,t)\!=\!\!\!%
\left[\begin{array}{c}\!\!\!-D_1(\Nbfhat\!\cdot\!\gradx)[n_{1}]
\!\!\!\\\!\!\!-D_2(\Nbfhat\!\cdot\!\gradx)[n_{2}]
\!\!\!\\\!\!\!-D_3(\Nbfhat\!\cdot\!\gradx)[c]
\end{array}\right]\!\!(\xbf,t)\!=\!\zerobf,
\nq
%
for $\xbf \in \Gambf\!=\!\partial \Omega$, i.e., on the boundary of the domain,
and for $t\in [0,t_{f}]$, 
assuming $D_{i}\neq 0$ or else the $D_{i}$ would not be used in the condition, 
where  $\Nbfhat(\xbf,t)$ is the normal to
the boundary, $-D_{i}\gradx[Y_{i}]$ is the flux of the $i$th
component and the diffusion matrix,
%
\[D(\xbf)
=\begin{bmatrix}D_1&0&0\\0&D_2&0\\0&0&D_3\end{bmatrix}
= [D_{i}\d_{i,j}]_{3\times{1}}\] is diagonal and 
inhomogeneous depending on the brain matter \cite{KRS99}, where
$\d_{i,j}$ is the Kronecker delta. Note that the no flux
condition at the boundary is motivated by the physical reality that
the brain is a finite and closed domain.

% Optimal Control Problem ******
\section{Optimal Control Problem}
%
%FBH 9/14:
A possible objective functional is the quadratic form of running and 
terminal costs,
%
\eq\label{Jdefeq}
J(u)\!\!\!&\!=\!\!&\!\!\half\int_{0}^{t_{f}}\!\!dt\int_{\Omega}\!\!d\xbf \left(
r n_1^2(\xbf,t)\!+\!s \left(u\!-\!u_0\right)^{2}(\xbf,t)
\right)
%FBH 9/10 add Terminal costs (also u0 may be distributed in x):
\nonumber\\\!\!&\!\!&\!\!
+\int_{\Omega}\!\!d\xbf\left(q_1n_1^2(\xbf,t_f)+q_3c^2(\xbf,t_f)\right)
.
\nq 
% FBH 9/10-11:
The goal is to minimize this functional with respect to
the drug input rate relative to some threshold rate $u_0$ and the terminal costs
%FBH 9/14:
at $t_f$, i.e., $\min_{u}\left[J(u)\right]$. Note that here $r>0$
is the tumor burden cost coefficient and $s>0$ is the drug delivery cost 
coefficient, while $q_1>0$ and  $q_3>0$ are the corresponding final costs. 
We are trying to minimize the density of tumor cells and
the drug delivery quadratic control term $(u(\xbf,t)-u_0(\xbf,t))^2$. 
%%
We could have chosen
a linear control which would have been less realistic, and also
would give rise to problems like singular control. In addition no
assumption is made about the control constraints, even though there
might be physical restriction on the amount of drugs that can be
administered.

% Vector Form *****
\section{Vector Form}
%
For the sake of \textit{brevity} we put the mathematical model in
vector form with vectors in boldface.

\subsection{Governing equations}
%
The vector state satisfies the PDE:
% FBH 9/11: -U --> +U
\eq\label{Ypde}
\frac{\partial\Ybf}{\partial t}=D\gradx^{2}[\Ybf]
+A\fbf(\Ybf)+B(\Ybf,t)\Ybf+\Ubf ,
%FBH 9/10 pure -\Ubf needed: +\Ubf ,
\nq
%
where
%
\[
A=\begin{bmatrix}a_{1}&0&0\\0&a_{2}&0\\0&0&a_{3}
\end{bmatrix} = [a_{i}\d_{i,j}]_{3\times{1}},
%\]
%
%\[
~~ 
\fbf(\Ybf)=\begin{bmatrix}f_{1}(n_{1})\\f_{2}(n_{2})\\ f_{3}(c) \end{bmatrix},
\]
%
\[
B(\Ybf,t)=-\begin{bmatrix}
\alpha_{1,2}n_{2}+\kappa_{1,3}c&0&0
\\0&\alpha_{2,1}n_{1}+\kappa_{2,3}c&0
\\0&0&0 \end{bmatrix},
\]
%
\eq\label{Ueq}
\Ubf(\xbf,t) =\begin{bmatrix} 0\\0\\
u(\xbf,t)
%FBH 9/10 replacing: u_0-u(\xbf,t) 
\end{bmatrix}
= U_{3}(\xbf,t)\ebf_{3} ,
\nq
% 
where $\ebf_{i}$ is the $i$th unit vector.

\subsection{Objective Functional}
%
The objective in vector form is
%
\eq\label{Jvecdef}
J[\Ybf,\Ubf]\!\!\!&\!\!=\!\!&\!\!\!\half\!\int_{0}^{t_{f}}\!\!\!dt\int_{\Omega}\!\!\!d\xbf%
\left(\Ybf^{\top}\!\!R\Ybf\!\!+\!\left(\Ubf\!\!-\!\!\Ubf_0\right)^{\top}\!\!\!S%
\left(\Ubf\!\!-\!\!\Ubf_0\right)\right)
\nonumber\\\!\!&&\!\!
+\half\int_{\Omega}\!\!d\xbf\left(\Ybf^{\top}\!\!Q\Ybf\right)\biggr|_{t=t_{f}}
,
\nq
%
where
%
$R
% FBH 9/10 delete:=\begin{bmatrix}r&0&0\\0&0&0\\0&0&0\\ \end{bmatrix} 
= r\ebf_{1}\ebf_{1}^{\top}$, 
$S
% FBH 9/10 delete:=\begin{bmatrix}0&0&0\\0&0&0\\0&0&s\\ \end{bmatrix} 
= s\ebf_{3}\ebf_{3}^{\top}$,
% FBH 9/10 add:
$Q = q_{1}\ebf_{1}\ebf_{1}^{\top}+q_{3}\ebf_{3}\ebf_{3}^{\top}$
and $\Ubf_0 = u_{0}(\xbf,t)\ebf_{3}$.

% Defining the pseudo-Hamiltonian
\section{Defining the Pseudo-Hamiltonian}
%
We have three vectors for the \textit{Lagrange multipliers}, two of
which are functions of space and time and one is independent of
time, needed to include the optimization constraints in the extended objective
for the state PDE (\ref{Ypde}), the boundary condition (\ref{Ybc})
and the initial condition (\ref{Yic}),
%
\eq\label{costate0}
\xibf(\xbf,t)=\begin{bmatrix} \xi_{1}\\\xi_{2}\\\xi_{3}
\end{bmatrix}\:,~
\etabf(\xbf,t)=\begin{bmatrix} \eta_{1}\\\eta_{2}\\\eta_{3}
\end{bmatrix}\:,~
\chibf(\xbf)=\begin{bmatrix} \chi_{1}\\\chi_{2}\\\chi_{3}
\end{bmatrix},
\nq
% 
i.e.,
$\xi_{i}=\xi_{i}(\xbf,t),\eta_{i}=\eta_{i}(\xbf,t)~\text{and}~\chi_{i}=\chi_{i}(\xbf)$,
for $i=1:3$.
Letting $\Zbf = (\Ybf, \Ubf, \xibf, \etabf,\chibf),$ define the \textit{pseudo-Hamiltonian} as,
%  CAUTION:  sign on Flux BC reversed 8/31:
\eq\label{Hameq}
\Hc(\Zbf)\!\!&\!\!=\!\!&\!\!\!%
\half\!\int_{0}^{t_{f}}\!\!\!dt\int_{\Omega}\!\!\!d\xbf%
\left(\Ybf^{\top}\!\!R\Ybf\!\!+\!\left(\Ubf\!\!-\!\!\Ubf_0\right)^{\top}\!\!\!S%
\left(\Ubf\!\!-\!\!\Ubf_0\right)\right)
\nonumber\\\!\!&&\!\!
+\half\int_{\Omega}\!\!d\xbf\left(\Ybf^{\top}\!\!Q\Ybf\right)\biggr|_{t=t_{f}}
\nonumber\\\nonumber
&&\!+\!\int_{0}^{t_{f}}\!\!\!dt\int_{\Omega}d\xbf~%
\xibf^{\top}\left( \frac{\partial\Ybf}{\partial t}\!-\!D\gradx^{2}[\Ybf]\!-\!A\fbf(\Ybf)\right.
\nonumber\\ %FBH 9/10-11 sign changeS on U:
&&\hspace*{1em}-B(\Ybf,t)\Ybf-\Ubf\biggr)
\nonumber\\
&&\!+\!\int_{0}^{t_{f}}\!\!dt\int_{\partial\Omega}d\Gambf~%\!
\etabf^{\top}\!\left(-D\left(\Nbfhat\!\cdot\!\gradx\right)[\Ybf]\right)
\nonumber\\
&&+\int_{\Omega}d\xbf~%
\chibf^{\top}\!\left(\Ybf\biggr|_{t=0}\!\!-\Ybf_{0}\right).%
\nq
%

% Lagrange multiplier method (?)
\section{Optimal Control Variational Formulation}
%
The \textit{calculus of variations} is used to find differential equation of 
optimal control for  the control, state and the co-state (adjoint or Lagrange multiplier) 
by seeking the functional critical point necessary conditions for the first variation 
\cite{MDG03,Kirk70} of the \textit{pseudo-Hamiltonian} $\Hc(\Zbf)$.    
% FBH 9/11: Not Necessary?
%Symbolically, these critical conditions are written
%$$
%\frac{\d\Hc}{\d\Ybf}=\zerobf, ~
%\frac{\d\Hc}{\d\xibf}=\zerobf, ~
%\frac{\d\Hc}{\d\etabf}=\zerobf, ~
%\frac{\d\Hc}{\d\chibf}=\zerobf, ~
%\frac{\d\Hc}{\d\Ubf}=\zerobf, ~
%$$
%where the $\d \Hc/\d \Zbf_{k}$ are the first variations with respect to 
%the five state/co-state/control vectors $\Zbf_{k}$ for $k=1\DDot{5}$ formulated in 
%functional derivative form.

\subsection{Pseudo-Hamiltonian First Variation}
%
Let the extended state vector be perturbed about the optimal trajectory $\Zbf^{*}$,
so that $\Zbf = \Zbf^{*}+\d\Zbf$, where $\d\Zbf$ is the perturbation.   Next expand
the pseudo-Hamiltonian 
$$
\Hc(\Zbf^{*}+\d\Zbf) = \Hc(\Zbf^{*}) + \d\Hc(\Zbf^{*},\d\Zbf) + O((\d\Zbf)^2).
$$
Neglecting the quadratic order terms, including the 2nd variation of $\Hc$, the
first variation is given by functional terms linear in $\d\Zbf$ using (\ref{Hameq}),
% 8/31 Change: %FBH 9/10 change U & FC costs:
\eq\label{H1stvar}
\d\Hc(\Zbf^{*},\d\Zbf)\!\!&\!\!=\!\!&\!\!\!%
\int_{0}^{t_{f}}\!\!\!\!dt\int_{\Omega}\!\!\!\!d\xbf
\left((\Ybf^{*})^{\top}\!R\d\Ybf\!%
+\!\left(\Ubf^{*}\!\!-\!\!\Ubf_0\right)^{\top}\!\!\!S\d\Ubf\right)
\nonumber\\\!\!&&\!\!
+\int_{\Omega}\!\!d\xbf\left((\Ybf^{*})^{\top}\!\!Q\d\Ybf\right)\biggr|_{t=t_{f}}
\nonumber\\\nonumber
&&\!+\!\int_{0}^{t_{f}}\!\!\!\!dt\int_{\Omega}\!\!\!\!d\xbf~%
\biggl((\xibf^{*})^{\top}\!\!\left(\d\Ybf_t\!-\!D\gradx^{2}[\d\Ybf]
\right.
\nonumber\\
&&\hspace*{1em}-A(\d\Ybf\!\cdot\!\grady)[\fbf](\Ybf^{*})-B(\Ybf^{*},t)\d\Ybf
\nonumber\\ %FBH 9/10 sign change on \d\Ubf; FBH 9/11 again:
&&\hspace*{1em}\left.-(\d\Ybf\!\cdot\!\grady)[B](\Ybf^{*},t)\Ybf^{*}-\d\Ubf\right)
\nonumber\\
&&\hspace*{1em}
+\d\xibf^{\top}\!\!\left(\Ybf^{*}_t\!-\!D\gradx^{2}[\Ybf^{*}]\!-\!A\fbf(\Ybf^{*})
\right.
\nonumber\\ %FBH 9/10 sign change on \Ubf; FBH 9/11 again::
&&\hspace*{1em}
\left.-B(\Ybf^{*},t)\Ybf^{*}-\Ubf^{*}\right)\biggr)%
\\
&&\!-\!\int_{0}^{t_{f}}\!\!\!\!\!dt\!\!\int_{\partial\Omega}\!\!\!\!\!d\Gambf%
\!\left((\etabf^{*})^{\top}\!\!D\left(\Nbfhat\!\cdot\!\gradx\right)[\d\Ybf]\right.
\nonumber\\
&&\hspace*{1em}
\left.+\d\etabf^{\top}\!D\left(\Nbfhat\!\cdot\!\gradx\right)[\Ybf^{*}]\right)%
\nonumber\\\nonumber
&&+\int_{\Omega}\!\!d\xbf
\left(\!(\chibf^{*})^{\top}\!\d\Ybf%
\!\!\!+\!\d\chibf^{\top}\!\!\left(\Ybf\biggr|_{t=0}\!\!-\Ybf_{0}\right)\!\!\right).
\nq
%
Before the critical conditions for first variation in (\ref{H1stvar}) can by used to 
obtain the extended state equations, the functional dependence of the higher 
derivatives in time and state of the extended state perturbations
must be eliminated on lower order terms by one or two integrations by parts,
i.e., by one,
%
\equnn
\int_{0}^{t_{f}}\!\!dt(\xibf^{*})^{\top}\!\d\Ybf_{t} = (\xibf^{*})^{\top}\!\d\Ybf\biggr|_{0}^{t_{f}}
\!-\!\int_{0}^{t_{f}}\!\!dt\d\Ybf^{\top}\xibf^{*}_{t}
\nqunn
%
and by two using the Green's formula \cite{Haberman83},
%
\equnn
\int_{\Omega}\!\!d\xbf(\xibf^{*})^{\top}\!\!D\gradx^{2}[\d\Ybf]%
\!\!&\!\!=\!\!&\!\!%
\int_{\Omega}\!\!d\xbf\d\Ybf^{\top}\!\!\gradx^{2}[D\xibf^{*}]
\\&&
+\int_{\partial\Omega}\!\!d\Gambf\biggl(
(\Nbfhat\!\cdot\!\gradx)[\d\Ybf^{\top}\!]D\xibf^{*}
\\&&
-\d\Ybf^{\top}(\Nbfhat\!\cdot\!\gradx)[D\xibf^{*}]\biggr).
\nqunn
%
Merging these identities with (\ref{H1stvar}), rearranging inner products and
collecting terms the extended state equations yields the following intermediate form:
%
% Optional: rearranged (\ref{H1stvar}) here.
%
% 8/31 Corrections:
\equnn%\label{H1stvar2}
\d\Hc(\Zbf^{*},\d\Zbf)\!\!&\!\!=\!\!&\!\!\!%
\int_{0}^{t_{f}}\!\!\!\!dt\int_{\Omega}\!\!d\xbf\:\d\Ybf^{\top}\!\!\left(
R\Ybf^{*}\!\!-\!\!\xibf^{*}_{t}\!\!-\!\!\gradx^{2}[D\xibf^{*}]
\right.
\\&&%
\!\!-\!\!\grady[\fbf](\Ybf^{*})A\xibf^{*}\!\!-\!\!B(\Ybf^{*}\!,t)\xibf^{*}
\\&&\left.
\!\!-\!\!\grady\![B](\Ybf^{*}\!,t))\DDot(\xibf^{*}\!(\Ybf^{*}\!)^{\top}%
\!\right)
\\&& %FBH 9/10 sign change on \d\Ubf with U_{0}; FBH 9/11 again for \xi:
\!+\!\int_{0}^{t_{f}}\!\!\!\!dt\int_{\Omega}\!\!d\xbf\:\d\Ubf^{\top}\!\!\left(
S\left(\Ubf^{*}\!\!-\!\!\Ubf_{0}\right)\!\!-\!\!\xibf^{*}\!\right)
\\
&&\!+\!\int_{0}^{t_{f}}\!\!\!\!dt\int_{\Omega}\!\!d\xbf\:\d\xibf^{\top}\!\!\left(
\Ybf^{*}_t\!\!-\!\!D\gradx^{2}[\Ybf^{*}]
\right.
\\&&\left. %FBH 9/10 sign change on \Ubf; FBH 9/11 again for \Ubf:
\!-\!A\fbf(\Ybf^{*}\!)\!\!-\!\!B(\Ybf^{*}\!,t)\Ybf^{*}-\Ubf ^{*}\!%
\right)
\\
&&\!-\!\int_{0}^{t_{f}}\!\!\!\!dt\int_{\partial\Omega}\!\!\!\!d\Gambf\:\d\etabf^{\top}
\!D\left(\Nbfhat\!\cdot\!\gradx\right)[\Ybf^{*}]
\\
&&\!+\!\int_{0}^{t_{f}}\!\!\!\!dt\int_{\partial\Omega}\!\!\!\!d\Gambf\:\d\Ybf^{\top}
\!\!\left(\Nbfhat\!\cdot\!\gradx\right)[D\xibf^{*}]
\\ % 8/31 reverse sign :
&&\!-\!\int_{0}^{t_{f}}\!\!\!\!dt\int_{\partial\Omega}\!\!\!\!d\Gambf
\!\left(\Nbfhat\!\cdot\!\gradx\!\right)\!\left[\d\Ybf^{\top}\!\right]\!D(\etabf^{*}\!\!\!+\!\!\xibf^{*})
\\
&&+\int_{\Omega}\!\!d\xbf\left(\d\chibf^{\top}\left(\Ybf^{*}\!\!-\!\!\Ybf_{0}(\xbf)\right)\!\right)\biggr|_{t=0}
\\
&&+\int_{\Omega}\!\!d\xbf\left(\d\Ybf^{\top}(\chibf^{*}\!\!-\!\!\xibf^{*})\right)\biggr|_{t=0}
\\
&&+\int_{\Omega}\!\!d\xbf\left(\d\Ybf^{\top}\!\!\left(\xibf^{*}\!\!+\!Q\Ybf\right)\right)\biggr|_{t=t_{f}},
\nqunn
%
%FBH 9/14:
where $A\DDot{B}$ denotes the trace of the matrix $AB$ or the double-dot product.


\subsection{State Equations}
%
The optimal state equation is recovered by setting the coefficient of $(\d\xibf^{*})^{\top}$
to zero:
%
\eq\label{Ystarpde}%FBH 9/10 sign change on \Ubf; FBH 9/11 again for \Ubf:
\frac{\partial \Ybf^{*}}{\partial t}=
D\gradx^{2}[\Ybf^{*}]+A\fbf(\Ybf^{*})+B(\Ybf^{*},t)\Ybf^{*}+\Ubf ^{*}
\nq
%
on  $\Omega \times (0,t_{f}]$,
with boundary conditions on $\partial\Omega \times [0,t_{f}]$ from the coefficient of $(\d\etabf^{*})^{\top}$, i.e.,
% 8/31 added minus for proper flux BC:
\eq\label{Ystarbc}
-D(\Nbfhat\!\cdot\!\gradx)[\Ybf^{*}](\xbf,t) = \zerobf,  
\nq
%
for $(\xbf,t) \in \partial\Omega\times[0,t_{f}]$
and with initial conditions on the interior $\Omega$ from the coefficient of $(\d\chibf^{*})^{\top}$,
i.e.,
%
\eq\label{Ystaric}
\Ybf^{*}(\xbf,0) =\Ybf_{0}(\xbf) 
\nq
%
for $\xbf \in \Omega$.  Due to the presence of the functions $\fbf(\Ybf)$ and $B(\Ybf,t)\Ybf$
the forward PDE (\ref{Ystarpde}) will be nonlinear.

\subsection{Degenerate Regular Optimal Control}
%
Since the control has been defined in (\ref{Ueq}) has been defined as only having 
on component, only the coefficient of $\d{U}_{3}$ is set to zero giving the regular,
but degenerate, control as
% 8/31 sign reversed: 
%FBH 9/10 More corrections for U_{0} and sign; FBH 9/11 again for \xi:
\eq\label{U3star}
U_{3}^{*}(\xbf,t) = U_{0,3}(\xbf,t)+\frac{1}{s}\xi_{3}^{*}(\xbf,t) ,
\nq
%
on  $\Omega \times [0,t_{f}]$,
provided $s \neq 0$.  Note that the control law with $\d{U}_{1}= 0 = \d{U}_{2}$
only requires solving for the 3rd component
of the first co-state vector $\xibf^{*}(\xbf,t)$.

\subsection{Co-State Equations}
%
Upon setting the functional coefficient of $(\d\Ybf^{*})^{\top}$ to zero yields
the primary co-state backward PDE:
% 8/31 -(\Nbfhat\!\cdot\!\gradx)[D\xibf^{*}] error deleted:
\eq\label{xistarbpde}
\zerobf\!\!&\!\!=\!\!&\!\!\frac{\partial\xibf^{*}}{\partial{t}}+\gradx^{2}[D\xibf^{*}]
+\grady[\fbf](\Ybf^{*})A\xibf^{*}
\nonumber\\&&
+B(\Ybf^{*},t)\xibf^{*}+\grady[B](\Ybf^{*},t))\DDot(\xibf^{*}(\Ybf^{*})^{\top})
\nonumber\\&&
-R\Ybf^{*},
\nq
%
for $(\xbf,t) \in \Omega\times[0,t_{f})$.
This PDE (\ref{xistarbpde}) is unidirectionally coupled to the state 
PDE (\ref{Ystarpde}), except that only the 3rd component $\xi^{*}_{3}(\xbf,t)$
is needed for the regular optimal control input $U^{*}_{3}(\xbf,t)$ from
(\ref{U3star}).
The boundary condition follows from setting the functional coefficient of 
$\d\Ybf(\xbf,t)$ for $\xbf$ on $\Gambf\!\!=\!\!\partial\Omega$ to zero, so
%
\eq\label{xistarbc}
(\Nbfhat\!\cdot\!\gradx)[D\xibf^{*}](\xbf,t) = \zerobf, ~~~ 
(\xbf,t) \in \partial\Omega\times[0,t_{f})
\nq
% 
and the final condition for this backward PDE follows from forcing the coefficient of
$\d\Ybf(\xbf,t_{f})$ to be zero on $\Omega$,
%
\eq\label{xistarfc}
\xibf^{*}(\xbf,t_{f}) = -Q\Ybf(\xbf,t_{f})
%FBH 9/10 not: \zerobf
.
\nq
%

The two other co-state vectors should not be needed, but satisfy rather simple
equations.   The 2nd co-state vector equation follows as the zero coefficient of 
$(\Nbfhat\!\cdot\!\gradx)[\d\Ybf^{\top}\!]$ on the state boundary 
$\Gambf\!\!=\!\!\partial\Omega$,
% FBH 9/10 Dless:
\equnn
\etabf^{*}(\xbf,t) = -\xibf^{*}(\xbf,t), ~~~ (\xbf,t) \in \partial\Omega\times[0,t_{f}].
\nqunn
%
The 3rd co-state vector equation follows as the zero coefficient of 
state initial condition $\d\Ybf(\xbf,0)$,
%
\equnn
\chibf^{*}(\xbf) = \xibf^{*}(\xbf,0), ~~~ \xbf \in \Omega.
\nqunn
%


% Algorithm for solving the coupled system
\section{Iterative Approximations for Coupled System}
%
We need to solve the system of equations developed in the previous
%FBH 9/14:
section, namely the \textit{state equations} (\ref{Ystarpde}) 
using the \textit{regular optimal control} (\ref{U3star})
and \textit{co-state equations} (\ref{xistarbpde}), with the understanding
that the state equations are forward equations while the co-state
equations are backward equations. 
The method a \textit{double shot} method, since model has two vector-valued
PDEs, so the method consists of one forward shot with (\ref{Ystarpde}) 
and and one backward shot with (\ref{xistarbpde}).
%%

The following steps can be used to
%FBH 9/14:
get approximate solutions:
%
\begin{enumerate}
%
\item
An initial guess for the first forward-backward shot iteration
is made for the control $U_{3}(\xbf,t) = U_{3}^{(1)}(\xbf,t)$ in (\ref{U3star}).  
Substituting it into the \textit{state
forward PDE} (\ref{Ystarpde}) solving for $\Ybf^{(1)}(\xbf,t)$, using the boundary condition
(\ref{Ystarbc}) and initial condition (\ref{Ystaric}).
%
\item
%FBH 9/14:
Next the approximate final condition  
$\xibf^{(1)}(\xbf,t_{f}) =-Q\Ybf^{(1)}(\xbf,t_{f})$ (\ref{xistarfc})
is used to start the backward \textit{co-state PDE}  (\ref{xistarbpde}) 
approximation using boundary conditions (\ref{xistarbc}) .
%
\item
Once  $\xibf^{(1)}(\xbf,t)$ is determined the \textit{regular optimal control
equation} is used to determine the updated value of the control,
%; FBH 9/11 sign change again for \xi:
$$U_{3}^{(2)}(\xbf,t) = U_{0,3}(\xbf,t)+\xi_{3}^{(1)}(\xbf,t)/s.$$
%
\item
This process is repeated for $\ell\!=\!2\DDot{L}$ double shot iterations
 until a
convergence criterion for sufficiently large $L$ is reached, e.g., the relative criterion
for the control,
% 8/31 edited:
\equnn
\left|\left|U_{3}^{(\ell)}(\xbf,t)\!\!-\!\!U_{3}^{(\ell\!-\!1)}(\xbf,t)\right|\right|
\!\!<\!\!\tol_{u}\!\left|\left|U_{3}^{(\ell\!-\!1)}(\xbf,t)\right|\right|,
\nqunn
% 
and, say,
% FBH 9/11:
\equnn
\left|\left|Y_{1}^{(\ell)}(\xbf,t)\!\!-\!\!Y_{1}^{(\ell\!-\!1)}(\xbf,t)\right|\right|
\!\!<\!\!\tol_{y}\!\left|\left|Y_{1}^{(\ell\!-\!1)}(\xbf,t)\right|\right|,
\nqunn
% 
%
for  $\ell=2\DDot{L}$ until satisfied, provided $||U_{3}^{(\ell-1)}(\xbf,t)|| \neq 0$ 
and $||Y_{1}^{(\ell-1)}(\xbf,t)|| \neq 0$,
where $\tol_{u} > 0$ and $\tol_{y} > 0$ are some prescribed tolerances.
\end{enumerate}

%A double shot method to solve the coupled system
\section{Double Shot Forward-Backward Computational Method}
%
In reality the problem is highly non-linear as are many problems in biology 
and we need numerical approximations of the solution. 
The main problem here is the fact that we have a forward state equation 
and a backward co-state equation.  
%FBH 9/14:
What we really do here is modify shooting methods \cite{Keller68} for initial-final-boundary value problems, where the starting aim is replaced by an estimate
of the full control law $U_{3}^{(\ell)}(\xbf,t)$ for the forward integration of the
state PDE (\ref{Ystarpde}) whose final approximation $\Ybf^{(\ell)}(\xbf,t_{f})$
serves as the backward aim (\ref{xistarfc}) for the backward integration of the  
co-state PDE  (\ref{xistarbpde}) producing an  approximation 
$\xibf^{(\ell)}(\xbf,t)$ whose third component is used to update (\ref{U3star}) 
the control law $U_{3}^{(\ell+1)}(\xbf,t)$. 
%%

Before we actually
go into the details of the shooting method, we will briefly discuss
the Crank-Nicolson, predictor-corrector central finite difference method
used to solve the nonlinear PDEs. The space is discretized as follows,
%
\[\xbf \rightarrow \xbf_{\jbf} =
[x_{{j_{i}},1}+(j_{i}-1)\cdot \Delta x_{i}]_{3\times 1}.\]
%
Here $\Delta x_{i}$ is the mesh size for state $i$ and
$\jbf=[j_{i}]_{3\times 1}$ where, $j_{i}=1\DDot{M}_{i}$ nodes per
node for states $i=1\DDot{3}$. For the forward state equation we have the
forward time discretization, $t \rightarrow t_{k} = k\Delta t$, for
$k=0\DDot{K}$ time steps where $\Delta t$ is the forward time step size,
$t_{0}=0$ and $t_{K}= t_{f}$.
Now we consider the vector state/co-state PDE system in general form,
% 8/31 Edited: 
\equnn
\Ybf^{*}_{t}\!\!&\!\!=\!\!&\!\!%
\Fbf(\xbf,t,\Ybf^{*}(\xbf,t),\Ubf^{*}(\xbf,t)),
\\
\zerobf\!\!&\!\!=\!\!&\xibf^{*}_{t} +
\Gbf(\xbf,t,\xibf^{*}(\xbf,t),\Ybf^{*}(\xbf,t)).
\nqunn 
%
In this equation, using the ordinate discretization
$\Ybf(\xbf_{\jbf},t_{k})\simeq \Ybf_{\jbf,k}$ and 
$\xibf(\xbf_{\jbf},t_{k})\simeq \xibf_{\jbf,k}$, the
relevant derivatives discretization are 
%
\equnn
\left(\Ybf^{*}_{t}\right)_ {\jbf,k+1/2 }\!\!\!&\!\!\!\simeq\!\!\!&\!\!\!%
\frac{\Ybf_{\jbf,k+1/2+1/2}-\Ybf_{\jbf,k+1/2-1/2}}{2(\Delta t/2)},
\\
\left(\xibf^{*}_{t}\right)_ {\jbf,k-1/2 }\!\!\!&\!\!\!\simeq\!\!\!&\!\!\!%
\frac{\xibf_{\jbf,k-1/2+1/2}-\xibf_{\jbf,k-1/2-1/2}}{2(\Delta t/2)}
\nqunn
%
and 
%
\equnn
\hspace*{-1ex}\left(\gradx^{2}[\Ybf^{*}\!]\right)_{\jbf,k+1/2}\!\!\!&\!\!\!\simeq\!\!\!&\!\!\!%
\frac{\Ybf_{\jbf+\ebf_{i},k+1/2}%
\!\!-\!\!2\Ybf_{\jbf,k+1/2}\!\!+\!\!\Ybf_{\jbf-\ebf_{i},k+1/2}}{(\Delta x_{i})^{2}},
\\
\hspace*{-1ex}\left(\gradx^{2}[\xibf^{*}\!]\right)_{\jbf,k-1/2}\!\!\!&\!\!\!\simeq\!\!\!&\!\!\!%
\frac{\xibf_{\jbf+\ebf_{i},k -1/2}%
\!\!-\!\!2\xibf_{\jbf,k-1/2}\!\!+\!\!\xibf_{\jbf-\ebf_{i},k-1/2}}{(\Delta x_{i})^{2}},
\nqunn
%
Consequently the forward and backward numerical schemes are given by
% 8/31 removed uneeded [ ] parens.:
\equnn
\Ybf_{\jbf,k+1}^{(\ell)}\!\!&\!\!=\!\!&\!\!\Ybf_{\jbf,k}^{(\ell)}+\Delta t
\Fbf_{\jbf,k+\frac{1}{2}}^{(\ell)} , \For k=0\DDot{K}\!\!-\!\!1,
\\
\xibf_{\jbf,k-1}^{(\ell)}\!\!&\!\!=\!\!&\!\!\xibf_{\jbf,k}^{(\ell)}+\Delta t
\Gbf_{\jbf,k-\frac{1}{2}}^{(\ell)} , \For k=\DDot{-1}\DDot{1},
\nqunn
% 
respectively, for $\ell = 1\DDot{L}$.   For each double shot $\ell$ the state starts from 
$\Ybf_{\jbf,0}^{(\ell)}= \Ybf_{0,\jbf,0} = \Ybf_{0}(\xbf_{\jbf},0)$ 
using 
%
\equnn
U_{3,\jbf,k}^{(\ell)} = U_{0,3,\jbf,k}+\xi_{3,\jbf,k}^{(\ell-1)}/s, \For k=0\DDot{K}\!\!-\!\!1,
\nqunn
% 
except when $\ell =1$ 
and the initial guess $U_{3,\jbf,k}^{(1)}$, e.g., $U_{0,3,\jbf,k}$, is used.  For each updated
forward state shot is completed, then the  backward co-state shot starts from
$\xibf_{\jbf,K}^{(\ell)}=Q\Ybf_{\jbf,K}^{(\ell)}$ 
using the whole state set $\Ybf_{\jbf,k}^{(\ell)}$ for 
$k=0\DDot{K}$.
The Crank-Nicolson temporal mid-point $\Fbf_{\jbf,k+\frac{1}{2}}^{(\ell)}$ state function is approximated by average with 
%
\equnn
\Ybf_{\jbf,k+\frac{1}{2}}^{(\ell)}=\frac{1}{2}\left(\Ybf_{\jbf,k+1}^{(\ell)}+\Ybf_{\jbf,k}^{(\ell)}
\right),
\nqunn
%
which can be used to construct finite differences for the derivatives, with a similar form
for $\Ubf_{\jbf,k+\frac{1}{2}}^{(\ell)}$, for $\ell=0\DDot{L}$.   Similarly, the backward
temporal mid-point $\Gbf_{\jbf,k-\frac{1}{2}}^{(\ell)}$ co-state function is approximated
with the average,
%
\equnn
\xibf_{\jbf,k-\frac{1}{2}}^{(\ell)}=\frac{1}{2}\left(\xibf_{\jbf,k}^{(\ell)}+\xibf_{\jbf,k-1}^{(\ell)}
\right),
\nqunn
%
and is used for the derivatives as well.

The no flux boundary conditions for both the state and co-state present some 
extra complexity, since the central differences of Crank-Nicolson are not
suitable at the  boundary if it is necessary to avoid using artificial external
points.   External points can be avoided by judicious use of forward and
backward differences of second order, matching the accuracy of the Crank-Nicolson
central differences.   In the simplest case of rectangular grids,  the discretized
no flux boundary conditions (\ref{Ystarbc},\ref{xistarbc}) with second order accuracy are 
% FBH 9/10 Dless:
\equnn
\zerobf\!\!&\!\!= \!\!&\!\!-\left((\Nbfhat\!\cdot\!\gradx)[\Ybf^{*}]\right)_{\jbf,k}^{(\ell)}
\\
\!\!&\!\!\simeq\!\!&\!\!-\frac{3\Ybf_{\jbf,k}^{(\ell)}-4\Ybf_{\jbf-\Nbf,k}^{(\ell)}
+\Ybf_{\jbf-2\Nbf,k}^{(\ell)}}{2|\Nbf\!\!\cdot\!\!\D\xbf|},
\\
\zerobf\!\!&\!\!= \!\!&\!\!\left((\Nbfhat\!\cdot\!\gradx)[(\xibf)^{*}]\right)_{\jbf,k}^{(\ell)}
\\
\!\!&\!\!\simeq\!\!&\!\!\frac{3(\xibf)_{\jbf,k}^{(\ell)}-4(\xibf)_{\jbf-\Nbf,k}^{(\ell)}
+(\xibf)_{\jbf-2\Nbf,k}^{(\ell)}}{2|\Nbf\!\!\cdot\!\!\D\xbf|},
\nqunn
%
respectively,
where $\Nbf = \Nbfhat_{\jbf,k}$, $\D\xbf=[\D{x}_{i}]_{3\times{1}} > \zerobf$, $D$ is not 
needed,
and, e.g., 
$\Ybf_{\jbf-\Nbf,k}^{(\ell)}\!\!=\!\!\Ybf^{(\ell)}(\xbf_{\jbf}-|\Nbf\!\cdot\!\D\xbf|\Nbf,t_{k})$.
For non-rectangular domains, interpolation would be needed to convert evaluations
to defined spatial nodes or else domain compatible grids should be used, e.g., for
%FBH 9/14: 
circular or spherical grid boundaries, $\Nbf = \ebf_{r}$, where here $r$ is the radius and
$\Nbf\cdot\D\xbf = \D{r}$.

%FBH 9/11 was not correct:
During each $\ell$th double shot, a prediction and corrections of the state and
co-state are used to account for the usual nonlinearities in the biological models \cite{FBH96}, stopping when the changes are sufficiently small.  

The above mentioned method is a sequential double shot method since
one shot is used to get $\Ybf_{\jbf,k}^{(\ell)}$ and a subsequent shot it used to get
$\xibf_{\jbf,k}^{(\ell)}$.    Alternately, a parallel two shot method could be used to
get an approximate solution by integrating for both $\Ybf_{\jbf,k}^{(\ell)}$ and
$\xibf_{\jbf,k}^{(\ell)}$ in the forward direction using a guess initial condition
for $\xibf_{\jbf,0}^{(\ell)}$ at $t_{0}=0$, with several genuine shooting method
shots until some $\ell^{*}$ shot where 
$||\xibf_{\jbf,K}^{(\ell^{*})}+Q\Ybf_{\jbf,K}^{(\ell^{*})}|| < \tol_{\xi}$, 
i.e., the final co-state value is small enough using some sufficiently small tolerance
$\tol_{\xi}$ to approximate the final condition (\ref{xistarfc}).  Interpolation should be used to accelerate the convergence of $\xibf_{\jbf,0}^{(\ell)}$
to some sufficiently small value.


% FBH 9/12:
\section{Test Results}
%
The double shot forward and backward algorithm outlined in the previous two sections 
has been tested on one space dimension, $x$, example with three state dimensions 
$\{Y_1^*\!\!=\!\!N_1^*,Y_2^*\!\!=\!\!N_2^*,Y_3^*\!\!=\!\!C^*\}$, plus
the drug input  control $U_3^*$.  The numerical parameter data come
from the BCNU drug  simulations for the brain of Wang et al. \cite{CHW99}
and the brain tumor modeling of Swanson \cite{KRS99} and Murray \cite{Murray2},
with some difficult to find parameters from Mansuri \cite{JM02} or from
reasonable estimates from other areas.  For example, diffusion diagonal vector
is $\Dbf = [$4.2e-3, 1.e-15, 2.2e-1$]$ cm$^{2}$ per day (normal tissue diffusion
is assumed to be insignificant), the quadratic cost coefficients are
$r=0.1 =s=q_1=q_3$, the net growth coefficients are $a=[$1.2e-2, 8.6e-7, 1.1e+1$]$
per day, the carrying capacities for tumor $k_1$ and normal $k_2$ tissues are scaled to 
one for the normal value and the interaction coefficients
$\{\alpha_{1,2}, \alpha_{2,1}, \kappa_{1,3}, \kappa_{2,3}\}$ are all given
the arbitrary value 1.0e-4.  The  initial states are given to be uniformly one for
the normal tissue, while the tumor density was assumed to be a spatial
Gaussian with spread 2.0e-2 about a mean of 0.0 with a weight of 1.e-3.
The initial drug concentration has a Gaussian spread of 1.0e-2 about a mean 0.0 with weight 1.e-4, while the threshold drug control $u_{0}(x,t)$ is similarly distributed,
but with weight 1.e-5.   

The results are given only for the tumor density $N_1^*= N_1^*(x,t)$ in 
Figure \ref{Fig1N1star} on the symmetric space interval $x \in [-1.0,+1.0]$ in 
centimeters over a $t_f = 5$ day treatment.   For this simple one space dimension
test example, we see that the optimal distribution of the tumor using
an optimal distribution of the drug delivery results  in the 16.2\% reduction
of the tumor density over this simulated five day drug treatment trial.
%
\begin{figure}[!hbtp]
\centering
{\includegraphics[width=3.25in]{acc05dim1n1.eps}}
\caption{%
Tumor density $N_1^*(x,t)$  versus the one-dimensional spatial coordinate $x$
with time $t$ at the rounded parameter values 
$\{0, t_{q1}\!\!=\!\!0.25t_f, t_{mid}\!=\!0.5t_f, t_{q3}\!=\!0.75t_f, t_f\}$, 
where $t_f= 5$ days.  The targeted tumor density rapidly decays in this simulated
5 day trial.
}%
\label{Fig1N1star}
\end{figure}
%


% Conclusion
\section{Conclusion}
%
The main interest of this paper was to provide the necessary
foundation to study the mechanism of drug delivery to the brain. We
have set up a fairly realistic distributed parameters model which
takes into account the spatial dependence of the state variables.
The main focus of the paper was to develop an algorithm to determine
%FBH 9/12:
the optimal drug delivery to brain tumors using an optimal distribution
of the drug about the original tumor site. 
%
This paper leaves room for a lot of new directions of work. 

\subsection{Future Directions}
One such
direction would be running a simulation to implement the algorithm,
using supercomputing tools.  Of course this would require more
realistic medical data.   Another important aspect that can be
examined is the effect of stochasticity, most notably the Gaussian
and Poisson type of noise.  The physical basis for such stochasticity
would be the phenomenon of metastasis, which gives rise to
additional tumor sites and also the side effects produced by the
drug.  We hope to examine all these aspects of the problem in future
works.

% Acknowledgments
\section{Acknowledgments}
% FBH 9/14:
The authors gratefully acknowledge, Dr. Herbert H. Engelhard of 
Neurological Surgery at the University of Illinois Medical Center 
for suggesting the problem and Professor William D. O'Neill of Bioengineering
at the University of Illinois at Chicago for putting us in contact.

% The Bibliography *****
\begin{thebibliography}{99}
%
\bibitem{HHE00}
H.H. Engelhard, ``Brain Tumors and the Blood-Brain Barrier,''
\textit{Neuro-Oncology, The Essentials}, pp. 49-53, Thieme Medical
Publishers, Inc., New York 2000.

\bibitem{RAG96}
R.A. Gatenby and E.T. Gawlinski, ``A Reaction-Diffusion Model of
Cancer Invasion,'' \textit{Cancer Research}, vol. 56, pp. 5745-5753,
December 15, 1996.

\bibitem{MDG03}
M.D. Gunzburger, \textit{Perspectives in Flow Control and Optimization},
SIAM Series in Design and Control, Philadelphia, 2003.
%%%%
\newpage
%%%%
\bibitem{Haberman83}
R. Haberman, 
\textit{Elementary Applied Partial Differential Equations with Fourier
Series and Boundary Value Problems},
Prentice-Hall, Englewood Cliffs, NJ, 1983.

\bibitem{FBH93}
F.B.Hanson and K.Naimipour, ``Convergence of Numerical Method for
Multistate Stochastic Dynamic Progamming,'' \textit{Proceedings of
International Federation of Automatic Control 12th World Congress},
vol. 9, pp.501-504, 1993.

\bibitem{FBH96}
F.B.Hanson, ``Techniques in Computational Stochastic Dynamic
Programming,'' \textit{Control and Dynamic Systems: Advances in Theory
and Applications}, vol. 76, pp.103-162, Academic Press, New York,
April 1996.

\bibitem{Keller68}
H. B. Keller, 
\textit{Numerical Methods for Two-Point Boundary Value Problems}, 
Gin Blaisdell, Waltham, MA, 1968.
%
\bibitem{Kirk70}
D. E. Kirk, 
\textit{Optimal Control Theory: An Introduction}, 
Prentice-Hall, Englewood Cliffs, NJ, 1970.

\bibitem{JM02}
J. Mansuri, \textit{The Modelling of Tumor Growth Using Reaction-Diffusion
Equations}, MSc Thesis, Oxford University, 2002.

\bibitem{Murray1}
J.D.Murray,  \textit{Mathematical Biology I:  An Introduction},
Springer, New York, 2002.

\bibitem{Murray2}
J.D.Murray,  \textit{Mathematical Biology  II:
Spatial Models and Biomedical Applications}, Springer, New York,
2003.

\bibitem{KRS99}
K.R. Swanson,  \textit{Mathematical Modeling of the Growth and Control of
Tumors}, PhD Thesis, University of Washington, 1999.

\bibitem{CHW98}
C-H. Wang and J. Li,  ``Three dimensional simulation of IgG delivery
to tumors,'' \textit{Chemical Engineering Science}, vol. 53, pp.
3579-3600, 1998.

\bibitem{CHW99}
C-H. Wang, J. Li, C.S. Teo and T. Lee, ``The delivery of BCNU to
brain tumors,''  \textit{Journal of Controlled Release}, vol. 61, pp.
21-41, 1999.

\bibitem{JJW02}
J. J. Westman, B. R. Fabijonas, D. L. Kern, and F. B. Hanson,
``Compartmental Model for Cancer Evolution: Chemotherapy and Drug
Resistance,'' \textit{Mathematical BioSciences}, pp. 1-20, May 2002,
under revision.

\bibitem{DEW96} D.E.Woodward, J.Cook, P.Tracqui, G.C.Cruywagen,
J. D Murray, E.C.Alvord, Jr., ``A mathematical model of glioma
growth: the effect of extent of surgical resection,'' \textit{Cell
Proliferation}, vol. 29, pp. 269-288, 1996.

\end{thebibliography}

\end{document}
