\documentclass[reqno]{amsart} \usepackage{longtable,amssymb} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2007(2007), No. 171, pp. 1--24.\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 2007 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2007/171\hfil Bang-bang and singular controls] {Seeking bang-bang solutions of mixed immuno-chemotherapy of tumors} \author[L. G. de Pillis, K. R. Fister, W. Gu, et al.\hfil EJDE-2007/171\hfilneg] {Lisette G. de Pillis, K. Renee Fister, Weiqing Gu, \\ Craig Collins, Michael Daub, David Gross, James Moore, Ben Preskill*} % in alphabetical order with professors first then students \address{Department of Mathematics, Harvey Mudd College, Claremont, CA 91711, USA} \email[L. G. de Pillis]{depillis@hmc.edu} \email[D. Gross]{david\_gross@hmc.edu} \email[W. Gu]{gu@math.hmc.edu} \email[J. Moore]{jmoore@hmc.edu} \email[B. Preskill]{bpreskill@hmc.edu} \address{Department of Mathematics and Statistics, Murray State University, Murray, KY 42071, USA} \email[C. Collins]{craig.collins@murraystate.edu} \email[K.R. Fister]{renee.fister@murraystate.edu} \address{Department of Mathematics and Statistics, Williams College, Williamstown, MA 01267, USA} \email[M. Daub]{Michael.W.Daub@williams.edu} \thanks{Submitted October 3, 2007. Published December 6, 2007.} \thanks{Supported by grant NSF-DMS-041-4011 from the National Science Foundation} \thanks{*The second line of authors are students who worked on the project.} \subjclass[2000]{37N25, 49J15, 49J30, 49N05, 62P10, 92C50} \keywords{Cancer modelling; mixed immuno-chemo-therapy; immunotherapy; \hfill\break\indent chemotherapy; linear optimal control} \begin{abstract} It is known that a beneficial cancer treatment approach for a single patient often involves the administration of more than one type of therapy. The question of how best to combine multiple cancer therapies, however, is still open. In this study, we investigate the theoretical interaction of three treatment types (two biological therapies and one chemotherapy) with a growing cancer, and present an analysis of an optimal control strategy for administering all three therapies in combination. In the situations with controls introduced linearly, we find that there are conditions on which the controls exist singularly. Although bang-bang controls (on-off) reflect the drug treatment approach that is often implemented clinically, we have demonstrated, in the context of our mathematical model, that there can exist regions on which this may not be the best strategy for minimizing a tumor burden. We characterize the controls in singular regions by taking time derivatives of the switching functions. We will examine these representations and the conditions necessary for the controls to be minimizing in the singular region. We begin by assuming only one of the controls is singular on a given interval. Then we analyze the conditions on which a pair and then all three controls are singular. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \allowdisplaybreaks \section{Introduction} The goal of applying optimal control theory to mathematical models representing the interaction between tumor, immune system, and chemotherapy is to determine the ideal mix of treatments that minimizes both tumor mass and negative effects upon the health of the patient. Recent research into the mixture of chemo- and immunotherapy regimens shows a great deal of potential for the success of such treatment schedules, c.f. \cite{Gogas2007}, \cite{Riker2007}, \cite{Rubinov2007}, \cite{Kofler2006}, \cite{Liu2006}, \cite{Neri2006}, \cite{Sinkovic2006}, \cite{Hermans2003}. The logic behind the development of a combination chemo-immunotherapy strategy is intuitive - use as little chemotherapy drug as possible to effectively kill tumor cells while utilizing immunotherapy to bolster the patient's immune system, thus strengthening the body's natural defenses against both the tumor cells and the dangerous side effects of the chemotherapy. The application of optimal control theory to mathematical models incorporating the interaction between tumors and treatments has provided valuable information in the past. Works by Kim \emph{et al.} \cite{Kim1977}, Swan and Vincent \cite{SwanVincent}, and Murray \cite{Mur} have successfully utilized control theory to maximize the effectiveness of chemotherapy against tumor cells and minimize the toxic effects of such treatments. Optimal control has also been effectively applied to immunotherapy models; for example, Swan \cite{Swan}, Kuznetsov and Knott \cite{Kuz}, and Kirschner \emph{et al.} \cite{Kirschner1997} have contributed research on applications to immunotherapy strategies for cancer and HIV. The more recent approach of combining chemo- and immunotherapy is being explored by several researchers; Kirschner and Panetta \cite{Kirschner1998} present a model detailing tumor-immune interaction with chemotherapy, while Burden \emph{et al.} \cite{Burden2004} apply quadratic control to that model. In addition, in the works by de Pillis and Radunskaya \cite{DePRad01}, \cite{DePRad02}, and de Pillis \emph{et al.} \cite{Wiseman2005}, \cite{dePillis2006} the authors explore various approaches to combining chemo- and immunotherapies through numerical simulation and the implementation of numerical linear controls.\\ \indent The model presented in this paper is an expansion of the model presented by de Pillis \emph{et al.} \cite{dePillis2006}. The modifications are introduced in response to newer research on the kinetics of IL-2 and immune cell populations. Another alteration to this model is the inclusion of constraints to limit both the tumor mass after treatment and the total concentration of lymphocytes during treatment. These additional constraints introduce a slight variation on the typical optimal control problem and must be dealt with by other means than the standard application of Pontryagin's Maximum/Minimum Principle \cite{Pontryagin1962}. We turn to works by Kirk \cite{Kirk1970}, Hartl \emph{et al.} \cite{HSV1995}, and Kamien and Schwartz \cite{KS1991} for the methods necessary to incorporate these constraints. The interest in applying control in a linear fashion to this model is somewhat pragmatic, given standard chemotherapy treatments. A widely used approach to cancer treatment is to give a maximum dosage of chemotherapy drug for some period of time, followed by a period of recuperation in which no drug is given. This type of treatment correlates theoretically with bang-bang control. However, when dealing with linear controls, we investigate the possibility of singular arcs and which conditions are necessary for those arcs to be optimal. In this paper, we use the Generalized Legendre-Clebsch conditions - as given by Krener \cite{Krener1977} - to generate the higher order necessary conditions for the optimality of the singular arcs. A more detailed discussion of this is given in Section 3. The outline of the paper is as follows. In Section 2, we present the model and discuss the model's components. Section 3 deals with the application of optimal control to the model, beginning with the description of the objective functional. The section continues with the proof of existence of an optimal control in this context and concludes with the conditions under which singular controls are optimal. Section 4 summarizes the conclusions reached from analysis of the data. In the appendices we provide full statements of certain theorems and detailed equations used to develop the work in this paper. \section{Three-Cell Three-Chemical Model}\label{Model} \subsection{Model Development and Cellular Dynamics} Medical advances have given doctors more options in cancer treatment. Traditional chemotherapy schedules may now be supplemented with various forms of immunotherapy. The existence of more options, however, can make it more challenging to find the best treatment. A mathematical model of all the processes involved can help to reveal improved treatment strategies. In this paper we analyze a model exploring the possible dynamics that can occur in different regions of parameter space. The hope is that the analysis will provide intuition into the interactions of the tumor, chemotherapy, and immunotherapy. The model tracks 6 quantities \begin{itemize} \item Tumor Cells, $T(t)$ (Units: Number of Cells) \item Natural Killer Cells, $N(t)$ (Units: Number of Cells per Liter) \item Circulating Lymphocytes, $C(t)$ (Units: Number of Cells per Liter) \item Tumor Specific $CD8^+$ $T$-Lymphocytes, $L(t)$ (Units: Number of Cells per Liter) \item Interleukin 2, $I(t)$ (Units: International Units (IUs) per Liter) \item Medicine, $M(t)$ (Units: Milligrams per Liter) \end{itemize} The circulating lymphocyte population represents all $B$ and $T$ lymphocytes in the bloodstream. Natural killer cells, which are also lymphocytes, are tracked as a separate population. The quantity $L$ represents the concentration of cells that have been activated by a tumor related antigen. \subsection{Tumor Equation ($T$)} Simple logistic growth of the tumor cell population, in the absence of medicine and immune interactions, is used. Death of tumor cells due to natural killer cells is given by a mass action term $cNT$, whereas death due to CTLs is given by a ratio dependent term, $D$. Note that unlike other cell populations, $T$ is measured as an absolute number and not a concentration. \subsection{Natural Killer Cell Equation ($N$)} A constant source term of Natural Killer cells differentiating from Circulating Lymphocytes and a linear death term are both assumed. We also assume that natural killer cells die when they have interacted with a tumor cell; so we include a mass action death term identical to that in the $T$ equation. For $N$, $L$, and $C$ the quantity given is cells per liter of blood. \subsection{Circulating Lymphocyte Equation} This population has the simplest dynamics: a constant source term and a linear death rate. \subsection{Tumor Specific T Cell (CD8+T CTL) Equation} We assume that these cells have a linear death rate, $-mL$, as well as a quadratic death rate, $-uL^2C$. The latter term represents the activity of regulatory T-cells, which are a subset of circulating lymphocytes. The CTLs may also die through interaction with the tumor and this is represented by a mass action term,$-qLT$. Interactions of the tumor with the larger lymphocyte populations, $N$ and $C$, stimulate CTL production. These stimulatory terms are represented by the two positive mass action terms. Additionally, tumor-related antigens stimulate the proliferation of CTLs through $j\frac{TL}{k+T}$. \subsection{Dynamics and Effects of Medicine} Once injected, medicine (chemotherapy) is assumed to have a linear decay rate. The medicine interacts with each of the four cell populations, $T$, $N$, $L$ and $C$ through a term of the form $K_X(1-e^{(-\delta_X M)})X$, \cite{Gardner}. For each cell population, this term represents cell death due to the medicine. \subsection{Dynamics and Effects of IL-2} Although naturally produced, the cytokine IL-2 is often used to treat cancer. This model assumes a linear decay rate and a constant source from circulating lymphocytes. Additionally, when a T-cell is stimulated by IL-2, the T-cell will secrete more IL-2 as represented by $\omega \frac{LI}{g_I+L}$. IL-2 also stimulates the proliferation of Natural Killer cells and CTLs. This stimulation is represented by the Michaelis-Menten terms, $p\frac{XI}{g+X}$, in each equation. Also, IL-2 inhibits the linear death term but stimulates the quadratic death term in the CTL equation. \subsection{The Equations} The system of differential equations describing the growth, death, and interactions of these populations with a che\-mo\-the\-rapy treatment is given by \begin{align} \frac{dT}{dt}& = aT(1-bT) - cNT - DT - K_T(1-e^{-\delta_T M})T \label{Teqn}\\ \frac{dN}{dt}& = f \big( \frac{e}{f}C -N\big) - pNT + \frac{p_NNI}{g_N+I} - K_N(1-e^{-\delta_N M})N \label{Neqn}\\ \frac{dL}{dt} &= -\frac{\theta mL}{\theta+I} - qLT + r_1NT + r_2CT + \frac{p_I L I}{g_I + I}-\frac{u_0L^2CI}{\kappa+I} \notag\\ &\quad + \frac{jTL}{k+T} -K_L(1-e^{-\delta_L M})L + \eta_1v_L(t)\label{Leqn}\\ \frac{dC}{dt}& = \beta\big( \frac{\alpha}{\beta} -C\big)- K_C(1-e^{-\delta_C M})C \label{Ceqn}\\ \frac{dM}{dt}& = -\gamma M + \eta_2v_M(t) \label{Meqn}\\ \frac{dI}{dt}& = -\mu_I I + \phi C +\frac{\omega LI}{\zeta+I} + \eta_3v_I(t) \label{Ieqn} \end{align} where $D=d \frac{(L/T)^{\ell} }{s/n^{\ell} +(L/T)^{\ell}}$, $T(0) = T_0$, $ N(0) = N_0$, $L(0) = L_0$, $I(0) = I_0$, $C(0) = C_0$, and $M(0) =M_0$. In Table 1, we have provided a summary of equation term descriptions, and in Table 2 we have a list of parameters with their units and biological interpretation. \begin{center} \begin{longtable}{|c|c|p{3.0in}|l|} \caption[Equation Descriptions]{Equation Descriptions}\label{eqntable} \\ \hline Eq. & Term & Description \\ \hline & $aT(1-bT)$ & Logistic tumor growth \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} $\frac{dT}{dt}$ & $-cNT$ & NK-induced tumor death \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} & $-DT$ & CD8+ T cell-induced tumor death \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} & $-K_T(1-e^{-\delta_T M})T$ & Chemotherapy-induced tumor death \\ \hline % & \cite{depillis:2006}, \cite{gardner:2000} \\ \hline & $eC$ & Production of NK cells from circulating lymphocytes \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} & $-fN$ & Natural killer breakdown \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} $\frac{dN}{dt}$ & $-pNT$ & Natural killer death by exhaustion of tumor-killing resources \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} & $\frac{p_NNI}{g_N + I}$ & Stimulatory effect of IL-2 on NK cells \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} & $-K_N(1-e^{\delta_N M})N$ & Death of NK cells due to medicine toxicity \\ \hline % & \cite{depillis:2006}, \cite{gardner:2000} \\ \hline & $-\frac{m\theta L}{\theta+I}$ & CD8+T cell breakdown \\ \cline{2-3} % & \cite{depillis:2006}, \cite{abbas:2005} \\ \cline{2-4} & $-qLT$ & CD8+T cell death by exhaustion of tumor-killing resources \\ \cline{2-3} % & \cite{depillis:2006}, \cite{kuznetsov:1994} \\ \cline{2-4} $\frac{dL}{dt}$ & $r_1NT$ & CD8+ T cell stimulation by NK-lysed tumor cell debris \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} & $r_2CT$ & Activation of naive CD8+T cells in the general lymphocyte population \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} & $\frac{p_I LI}{g_I + I}$ & Stimulatory effect of IL-2 on CD8+T cells \\ \cline{2-3} % & \cite{depillis:2006}, \cite{kirschner:1998} \\ \cline{2-4} & $-\frac{u_0L^2CI}{\kappa + I}$ & Breakdown of surplus CD8+T cells in the presence of IL-2 \\ \cline{2-3} % & \cite{depillis:2006}, \cite{abbas:2005} \\ \cline{2-4} & $\frac{jTL}{k + T}$ & CD8+ T cell stimulation by CD8+ T cell-lysed tumor cell debris \\ \cline{2-3} % & \cite{kuznetsov:1994} \\ \cline{2-4} & $-K_L(1-e^{-\delta_L M})L$ & Death of CD8+ T cells due to medicine toxicity \\ \cline{2-3} % & \cite{depillis:2006}, \cite{gardner:2000} \\ \hline & $\eta_1 v_L(t)$ & External TIL therapy, controllable \\ \hline % Added to table LdeP & $\alpha$ & Lymphocyte synthesis in bone marrow \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} $\frac{dC}{dt}$ & $-\beta C$ & Lymphocyte breakdown \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} & $-K_C(1-e^{-\delta_C M})C$ & Death of lymphocytes due to medicine toxicity \\ \hline % & \cite{depillis:2006}, \cite{gardner:2000} \\ \hline & $-\gamma M$ & Excretion and breakdown of medicine \\ \cline{2-3} % & \cite{depillis:2006} \\ \hline $\frac{dM}{dt}$ & $\eta_2 v_M(t)$ & External chemotherapy, controllable \\ \hline % Added to table LdeP & $-\mu_I I$ & IL-2 breakdown \\ \cline{2-3} % & \cite{depillis:2006} \\ \cline{2-4} $\frac{dI}{dt}$ & $\phi C$ & Production of IL-2 due to naive CD8+T cells and CD4+T cells \\ \cline{2-3} % & \cite{abbas:2005} \\ \cline{2-4} & $\frac{\omega LI}{\zeta + I}$ & Production of IL-2 from activated CD8+T cells \\ \cline{2-3} % & \cite{kirschner:1998} \\ \hline & $\eta_3 v_I(t)$ & External IL-2, controllable \\ \hline % Added to table LdeP \end{longtable} \end{center} \begin{center} \begin{longtable}{|c|c|p{68mm}|c|} \caption[Parameter Descriptions]{Parameter Descriptions}\label{paramdeftable} \\ \hline Eq. & Param. & Description & Units \\ \hline & $a$ & Growth rate of tumor & 1/day \\ \cline{2-4} & $b$ & Inverse of carrying capacity & 1/cells \\ \cline{2-4} $\frac{dT}{dt}$ & $c$ & Rate of NK-induced tumor death& liter/(cells$\cdot$day) \\ \cline{2-4} & $K_T$ & Rate of chemotherapy-induced tumor death & 1/day \\ \cline{2-4} & $\delta_T$ & Medicine efficacy coefficient & liter/mg \\ \hline & $e/f$ & Ratio of rate of NK cell creation with rate of cell death & unitless \\ \cline{2-4} & $f$ & Rate of NK cell death & 1/day \\ \cline{2-4} & $p$ & Rate of NK cell death due to tumor interaction & 1/(cells$\cdot$day) \\ \cline{2-4} $\frac{dN}{dt}$ & $p_N$ & Rate of IL-2 induced NK cell genesis & 1/day \\ \cline{2-4} & $g_N$ & Concentration of IL-2 for half-maximal NK cell genesis & IU/liter \\ \cline{2-4} & $K_N$ & Rate of NK depletion from medicine toxicity & 1/day \\ \cline{2-4} & $\delta_N$ & Medicine toxicity coefficient & liter/mg \\ \hline & $m$ & Natural decay rate of CD8+T cells & 1/day \\ \cline{2-4} & $\theta$ & Concentration of IL-2 to halve effectiveness of CD8+T self-regulation & IU/liter \\ \cline{2-4} & $q$ & Rate of CD8+T cell death due to tumor interaction & 1/(cells$\cdot$day) \\ \cline{2-4} & $r_1$ & Rate of NK-lysed tumor cell debris activation of CD8+T cells & 1/(cells$\cdot$day) \\ \cline{2-4} & $r_2$ & Rate of CD8+T cell production from circulating lymphocytes & 1/(cells$\cdot$day) \\ \cline{2-4} $\frac{dL}{dt}$ & $p_I$ & Rate of IL-2 induced CD8+T cell activation & 1/day \\ \cline{2-4} & $g_I$ & Concentration of IL-2 for half-maximal CD8+T cell activation & IU/liter \\ \cline{2-4} & $u_0$ & CD8 self-limitation feedback coefficient & liter$^2$/(cells$^2\cdot$day) \\ \cline{2-4} & $\kappa$ & Concentration of IL-2 for half-maximal IL-2-dependent CD8+T self-regulation& IU/liter \\ \cline{2-4} & $j$ & Rate of CD8+T-lysed tumor cell debris activation of CD8+T cells & 1/day \\ \cline{2-4} & $k$ & Tumor size for half-maximal CD8+T-lysed debris CD8+T activation & cells \\ \cline{2-4} & $K_L$ & Rate of CD8+T depletion from medicine toxicity & 1/day \\ \cline{2-4} & $\delta_L$ & Medicine toxicity coefficient & liter/mg \\ \cline{2-4} & $\eta_1$ & TIL therapy administration rate & 1/day \\ \cline{2-4} & $v_L(t)$ & Externally administered TIL CD8+T therapy concentration & cells/liter \\ \hline & $\alpha/\beta$ & Ratio of rate of circulating lymphocyte production to death rate & cells/liter \\ \cline{2-4} $\frac{dC}{dt}$ & $\beta$ & Rate of decay of circulating lymphocytes & 1/day \\ \cline{2-4} & $K_C$ & Circulating lymphocyte-toxicity of medicine & 1/day \\ \cline{2-4} & $\delta_C$ & Medicine toxicity coefficient & liter/mg \\ \hline & $\gamma$ & Rate of decay of medicine & 1/day \\ \cline{2-4} $\frac{dM}{dt}$ & $\eta_2$ & Chemotherapy administration rate & 1/day \\ \cline{2-4} & $v_M(t)$ & Externally administered chemotherapy concentration & mg/liter \\ \hline & $\mu_I$ & Rate of decay of IL-2 & 1/day \\ \cline{2-4} & $\phi$ & Rate of IL-2 production from circulating lymphocytes & IU/(cells$\cdot$day) \\ \cline{2-4} & $\omega$ & Rate of IL-2 production from CD8+T cells & IU/(cells$\cdot$day) \\ \cline{2-4} $\frac{dI}{dt}$ & $\zeta$ & Concentration of IL-2 for half-maximal CD8+T IL-2 production & IU/liter \\ \cline{2-4} & $\eta_3$ & Exogenous IL-2 administration rate & 1/day \\ \cline{2-4} & $v_I$ & Externally administered IL-2 conentration & IU/liter \\ \hline & $d$ & Immune system strength coefficient & liter/day \\ \cline{2-4} & $l$ & Immune strength scaling coefficient & unitless \\ \cline{2-4} $D$ & $s$ & Concentration of CD8+T cells in tumor for half-maximal tumor death & cells/liter \\ \cline{2-4} & $n$ & Number of CD8+T cells that infiltrate tumor & cells \\ \hline \end{longtable} \end{center} In addition to the system of differential equations, there are two constraints associated with this model. The first is a terminal condition in which the tumor population is required to be limited by an upper bound, \begin{equation} T(t_f)\le\Omega_T,\label{Tcon} \end{equation} where $\Omega_T$ is constant. The second constraint is a condition on the control $v_M$ in which the total drug administered is limited by a constant. The constant is divided by 2 for mathematical convenience. \begin{equation} \frac{\tau}{2}-\int_0^{t_f}v_M(t) \ge 0,\label{vMcon} \end{equation} where $\tau$ is constant. \section{Linear Control} The goal of implementing linear control on the model described in Section 2 is to determine theoretically the optimal treatment schedule for a cancer patient. We prove the existence of such a control using the Filippov-Cesari Theorem, as stated in Hartl, Sethi, and Vickson \cite{HSV1995}. Constraints \eqref{Tcon} and \eqref{vMcon} prevent the direct application of Pontryagin's Minimum Principle to find the first order necessary conditions for the control to be optimal. We therefore use conditions given by Kamien and Schwartz \cite{KS1991} and Hartl et al. \cite{HSV1995} to deal with the terminal inequality conditions generated by constraint \eqref{Tcon}. Constraint \eqref{vMcon} is incorporated into the state system using a method detailed by Kirk \cite{Kirk1970}. When implementing linear control, we must investigate the possibility of singular arcs. A singular arc is one for which one or more of the control variables $v_{\alpha}$ satisfies \begin{equation*} \frac{\partial H}{\partial v_{\alpha}} =0, \end{equation*} where $H$ is the Hamiltonian. In this situation, first order necessary conditions are inadequate; therefore, we use the generalized Legendre-Clebsch conditions (see Appendix A, \eqref{eqn:GLC}) as given by Krener \cite{Krener1977} to generate these higher order necessary conditions for optimality. \subsection{Objective Functional} Now, seeking a bang-bang solution, we wish to minimize the objective functional \begin{equation} J(v_L,v_M,v_I) = \int_0^{t_f} \left( T(t) + \epsilon_L v_L(t) + \epsilon_M v_M(t) + \epsilon_I v_I(t) \right)dt, \label{ObjFunc} \end{equation} which is linear in the three controls and where $\epsilon_L$, $\epsilon_M$ and $\epsilon_I$ are weight factors. \subsection{Existence} We first establish the existence of an optimal control building on the existence theorem of Filippov-Cesari, the full statement of which is provided in Appendix A, Theorem \eqref{eqn:FCT}. \begin{theorem}[Existence of a Linear Optimal Control] Given the objective functional \eqref{ObjFunc}, subject to system \eqref{Teqn}-\eqref{Ieqn} with $T(0)=T_0$, $N(0)=N_0$, $L(0)=L_0$, $C(0)=C_0$, $M(0)=M_0$, $I(0)=I_0$, $T(t_f)\le\Omega_T$, and $\frac{\tau}{2}-\int_0^{t_f}{v_M}dt\ge0$, and the control set \begin{equation*} U = \{v_L(t),v_M(t), v_I(t) \text{ piecewise cont. }: 0 \le v_L(t),v_M(t), v_I(t) \le 1, \forall t \in [0,t_f]\}, \end{equation*} the conditions $(1)-(4)$ of Theorem \eqref{eqn:FCT} are met, and therefore there exists an optimal control $\vec{V}^*(t)=(v_L^*(t), v_M^*(t), v_I^*(t))$ such that $$\min_{\vec{V}\in [0,1]} J(V)=J(V^*).$$ \end{theorem} Applying the notation of Theorem \ref{eqn:FCT} to the optimal control problem \eqref{Teqn}-\eqref{Ieqn}, \eqref{ObjFunc}, we have \begin{equation*} \mathbf{x}= \begin{pmatrix} T \\ N \\ L \\ C \\ M \\ I \end{pmatrix} \end{equation*} \begin{align*} &N(\mathbf{x},t) \\ &={\footnotesize \begin{pmatrix} T(t) + \epsilon_L v_L(t) +\epsilon_M v_M(t) + \epsilon_I v_I(t)\\ aT(1-bT) - cNT - DT - K_T(1-e^{-\delta_T M})T \\ f \left( \frac{e}{f}C -N\right) - pNT + \frac{p_NNI}{g_N+I} - K_N(1-e^{-\delta_N M})N \\ -\frac{\theta mL}{\theta+I} - qLT + r_1NT + r_2CT +\frac{p_I L I}{g_I + I}-\frac{u_0L^2CI}{\kappa+I} + \frac{jTL}{k+T} -K_L(1-e^{-\delta_L M})L + \eta_1v_L(t) \\ \beta\left( \frac{\alpha}{\beta} -C\right)- K_C(1-e^{-\delta_C M})C \\ -\gamma M + \eta_2v_M(t) \\ -\mu_I I + \phi C +\frac{\omega LI}{\zeta+I}+\eta_3v_I(t) \end{pmatrix} } \end{align*} where $\tilde{\gamma} \leq 0$, and $v_L,v_M,v_I \in \Omega(x,t)$. \begin{proof} We know there exists an admissible solution pair for the state and controls as seen in previous work, \cite{Immune}. For the second condition, we define $w_1=$ \begin{equation*} {\footnotesize \begin{pmatrix} T(t) + \epsilon_L v_{L_1}(t) +\epsilon_M v_{M_1}(t) + \epsilon_I v_{I_1}(t)\\ aT(1-bT) - cNT - DT - K_T(1-e^{-\delta_T M})T \\ f \left( \frac{e}{f}C -N\right) - pNT + \frac{p_NNI}{g_N+I} - K_N(1-e^{-\delta_N M})N \\ -\frac{\theta mL}{\theta+I} - qLT + r_1NT + r_2CT +\frac{p_I L I}{g_I + I}-\frac{u_0L^2CI}{\kappa+I}+ \frac{jTL}{k+T} -K_L(1-e^{-\delta_L M})L + \eta_1v_{L_1}(t)\\ \beta\left( \frac{\alpha}{\beta} -C\right)- K_C(1-e^{-\delta_C M})C \\ -\gamma M + \eta_2v_{M_1}(t) \\ -\mu_I I + \phi C +\frac{\omega LI}{\zeta+I}+\eta_3v_{I_1}(t) \end{pmatrix},} \end{equation*} and $w_2=$ \begin{equation*} {\footnotesize \begin{pmatrix} T(t) + \epsilon_L v_{L_2}(t) +\epsilon_M v_{M_2}(t) + \epsilon_I v_{I_2}(t)\\ aT(1-bT) - cNT - DT - K_T(1-e^{-\delta_T M})T \\ f \left( \frac{e}{f}C -N\right) - pNT + \frac{p_NNI}{g_N+I} - K_N(1-e^{-\delta_N M})N \\ -\frac{\theta mL}{\theta+I} - qLT + r_1NT + r_2CT +\frac{p_I L I}{g_I + I}-\frac{u_0L^2CI}{\kappa+I}+ \frac{jTL}{k+T} -K_L(1-e^{-\delta_L M})L + \eta_1v_{L_2}(t)\\ \beta\left( \frac{\alpha}{\beta} -C\right)- K_C(1-e^{-\delta_C M})C \\ -\gamma M + \eta_2v_{M_2}(t) \\ -\mu_I I + \phi C +\frac{\omega LI}{\zeta+I}+\eta_3v_{I_2}(t) \end{pmatrix} } \end{equation*} for some $\tilde{\gamma_1},\tilde{\gamma_2} \leq 0$, and $v_{L_j},v_{M_j},v_{I_j} \in\Omega(x,t)$, with $j=1,2$. We then let $w_3 = \lambda w_1 + (1-\lambda)w_2$ for $\lambda \in [0,1]$. To prove that $N(\mathbf{x},t)$is convex, we need to show that $w_3 \in N(\mathbf{x},t)$ However, since the controls appear linearly, we note that this occurs naturally. For the third condition we need to show that there exists a number $\delta$ such that $\| x \| \leq \delta, \;\forall t \in [0,t_f]$ and all admissible pairs $(\vec{\mathbf{x}},\vec{\mathbf{V}})$. We need to determine an upper bound on the right hand sides of the differential equations \eqref{Teqn}-\eqref{Ieqn}. To do this, we first consider the right hand side of \eqref{Teqn}, the tumor cell population. With $T_{\rm max}$ as an upper bound solution associated with $T$ and $T(t)\ge 0$, then $T_{\rm max}= T_0e^{at_f}$. By using the bound $T_{\rm max}$, we can form a set of supersolutions for system \eqref{Teqn}-\eqref{Ieqn}. This set of supersolutions $\bar{T}$, $\bar{N}$, $\bar{L}$, $\bar{C}$, $\bar{M}$, $\bar{I}$ of \begin{align} \frac{d\overline{T}}{dt} & = a\overline{T} \notag\\ \frac{d\overline{N}}{dt} & = p_N\overline{N} + e\overline{C}\notag\\ \frac{d\overline{L}}{dt} & = r_1\overline{N}T_{\rm max} + (p_I+j)\overline{L} + r_2\overline{C}T_{\rm max} + \eta_1v_L\notag\\ \frac{d\overline{C}}{dt} & = \alpha \notag\\ \frac{d\overline{M}}{dt} & = \eta_2v_M\notag\\ \frac{d\overline{I}}{dt} & = \omega\overline{L}+\phi\overline{C}+\eta_3v_I \notag \label{eq:superarrayfull} \end{align} is bounded on a finite time interval. We see that system \eqref{Teqn}-\eqref{Ieqn} can be written as \[ \begin{pmatrix} \bar{T}\\ \bar{N}\\ \bar{L}\\ \bar{C}\\ \bar{M}\\ \bar{I} \end{pmatrix}' = \begin{pmatrix} a & 0 & 0 & 0 & 0 & 0 \\ 0 & p_N & 0 & e & 0 & 0 \\ 0 & r_1T_{\rm max} & p_I+j & r_2T_{\rm max} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \omega & \phi & 0 & 0 \end{pmatrix} \begin{pmatrix} \bar{T}\\ \bar{N}\\ \bar{L}\\ \bar{C}\\ \bar{M}\\ \bar{I} \end{pmatrix} + \begin{pmatrix} 0\\ 0\\ \eta_1v_L\\ \alpha \\ \eta_2v_M\\ \eta_3v_I \end{pmatrix} \] Since this supersolution system involves only constants then it has a finite upper bound. Letting this upper bound be $\delta$ satisfies the third condition. The fourth condition is satisfied by definition, since $00\\ 1 & \text{if }\Phi_L <0\\ \text{singular }&\text{if }\Phi_L =0, \end{cases} \label{vLrep}\\ v_M(t) & = \begin{cases} 0 & \text{if }\Phi_M >0\\ 1 & \text{if }\Phi_M <0\\ \text{singular}&\text{if }\Phi_M =0, \end{cases}\label{vMrep}\\ v_I(t) & = \begin{cases} 0 & \text{if }\Phi_I >0\\ 1 & \text{if }\Phi_I <0\\ \text{singular}&\text{if }\Phi_I =0. \end{cases} \label{vIrep} \end{align} \end{theorem} \begin{proof} The Lagrangian associated with this problem is \begin{align*} \mathcal{L} &= H - W_1(t) v_L(t) - W_2(t) (1- v_L(t)) - W_3(t) v_M(t) - W_4(t) (1- v_M(t))\\ &\quad - W_5(t) v_I(t) - W_6(t) (1-v_I(t)), \end{align*} where $H$ is the Hamiltonian given by \begin{align*} H &= T(t) + \epsilon_L v_L(t) + \epsilon_M v_M(t) + \epsilon_I v_I(t) \\ &\quad + \lambda_1(aT(1-bT) - cNT- DT - K_T(1-e^{-\delta_T M})T) \\ &\quad + \lambda_2(f \big( \frac{e}{f}C -N\big) - pNT + \frac{p_NNI}{g_N+I} - K_N(1-e^{-\delta_N M})N) \\ &\quad + \lambda_3(-\frac{\theta mL}{\theta+I}- qLT + r_1NT + r_2CT + \frac{p_I L I}{g_I + I} -\frac{u_0L^2CI}{\kappa+I} + \frac{jT}{k+T}L \\ &\quad -K_L(1-e^{-\delta_L M})L + \eta_1v_L(t)) \\ &\quad +\lambda_4(\beta\big( \frac{\alpha}{\beta} -C\big) - K_C(1-e^{-\delta_C M})C)+\lambda_5(-\gamma M + \eta_2v_M(t)) \\ &\quad + \lambda_6(-\mu_I I + \phi C +\frac{\omega LI}{\zeta+I}+\eta_3v_I(t)) + \lambda_7v_M(t) \end{align*} The adjoint equations are formed from differentiating the Hamiltonian with respect to the corresponding state variables as $ \frac{d \lambda_1}{dt} =-\frac{\partial H}{\partial T}$. Since system \eqref{adjT}-\eqref{adjA} has bounded coefficients and the solutions are bounded on the finite time interval, we know that adjoint variables satisfying \eqref{adjT}-\eqref{adjA} exist from Lukes \cite{Lukes1982}, p.182. Also, the final conditions are free for all variables with the exception of $T(t)$ and $A(t)$; thus we have that $\lambda_i(t_f)=0$, for $i=2, \ldots, 6$. The additional conditions imposed by the terminal inequality constraints are taken from Kamien and Schwarz \cite{KS1991}. We find the representations of \eqref{vLrep}-\eqref{vIrep} of the optimal controls by looking at the coefficients of the controls in H; i.e. we consider the sign of $\Phi_L = \frac{\partial H}{\partial v_L} $ to determine the values of $v_L$. \end{proof} The characterization of the controls in singular regions is determined by taking time derivatives of the switching functions. We will examine these representations and the conditions necessary for the controls to be minimizing in the singular region provided that the representations of the singular controls satisfy the control bounds. Note that each interval of singularity is a subinterval of $[0,t_f]$. We will begin by assuming only one of the controls is singular on a given interval in each of Theorems 3.3-3.5. \begin{theorem} \label{thm3.3} If $v_L$ is singular on the interval $(s_L, t_L)$, the singularity is of degree two and the representation of the control is given by \begin{equation*} v_L = -\frac{P_L}{Q_L}, \quad Q_L\neq 0, \end{equation*} where $Q_L$ is given by \begin{equation} Q_L = \eta_1\big[\lambda_1 T\frac{\partial ^2D}{\partial{L^2}}- 2u_0C\frac{\epsilon_L}{\eta_1} \big(\frac{I}{\kappa+I}\big) \big]\label{Q_L} \end{equation} and $P_L$ is given in Appendix B. Furthermore, if $v_L$ is minimizing on this interval, it is necessary that \begin{equation*} \lambda_1T\frac{\partial^2 D}{\partial L^2}\le 2u_0C\frac{\epsilon_L}{\eta_1}\left(\frac{I}{\kappa+I}\right). \end{equation*} \end{theorem} \begin{proof} On the interval $(s_L,t_L)$, we know that time derivatives of $\Phi_L$ are identically $0;$ i.e. $\frac{d}{dt}\Phi_L =0$, $\frac{d^2}{dt^2}\Phi_L =0$, etc. In addition, we have that $\lambda_3 = -\frac{\epsilon_L}{\eta_1}$ in this region. We use this information to find $v_L$. To begin, notice that \begin{equation*} \frac{d}{dt}\Phi_L = \eta_1\frac{d}{dt}\lambda_3 = 0, \quad\text{or }\frac{d\lambda_3}{dt} =0. \end{equation*} Then, from the corresponding adjoint equation \eqref{adjL}, we have that \begin{align*} \frac{d\lambda_3}{dt} &=\lambda_1\left(\frac{\frac{s}{n^{\ell}}d\ell \big( \frac{L}{T} \big)^{\ell-1}}{\big[ \frac{s}{n^{\ell}} + \big( \frac{L}{T} \big)^{\ell}\big]^2}\right) - \lambda_6 \Big( \frac{\omega I}{\zeta+I} \Big) \\ &\quad+ \lambda_3 \left(\frac{\theta m}{\theta+I}+qT - \frac{p_I I}{g_I + I} + \frac{2u_0LCI}{\kappa+I} -\frac{jT}{k+T} + K_L (1-e^{-\delta_L M}) \right) \\ & = 0 \end{align*} Note that \begin{equation*} \frac{\partial D}{\partial L} T = \frac{\frac{s}{n^{\ell}}d\ell \big( \frac{L}{T} \big)^{\ell-1}}{\big[ \frac{s}{n^{\ell}} + \big( \frac{L}{T} \big)^{\ell}\big]^2} \end{equation*} This will simplify the notation in the following calculation. Taking a second time derivative yields \begin{align*} \frac{d^2\lambda_3}{dt^2} & = \frac{d\lambda_1}{dt}\left( \frac{\partial D}{\partial L}T\right) + \lambda_1 T \left[ \frac{\partial ^2D}{\partial L^2}\frac{dL}{dt} + \frac{\partial ^2 D}{\partial T\partial L}\frac{dT}{dt}\right]\\ &\quad + \lambda_1 \left(\frac{\partial D}{\partial L}\frac{dT}{dt}\right) - \frac{d\lambda_6}{dt}\left(\frac{\omega I}{\zeta +I}\right) -\lambda_6\left(\frac{\omega\zeta}{(\zeta +I)^2} \right)\frac{dI}{dt}\\ &\quad +\lambda_3\Big[-\frac{\theta m}{(\theta +I)^2}\frac{dI}{dt} + q\frac{dT}{dt} -\frac{p_Ig_I}{(g_I +I)^2}\frac{dI}{dt} + 2u_0LC\left(\frac{\kappa}{(\kappa +I)^2}\right)\frac{dI}{dt}\\ &\quad +2u_0\left(\frac{I}{\kappa +I}\right)\left[L\frac{dC}{dt}+C\frac{dL}{dt} \right]- \frac{jk}{(k+T)^2}\frac{dT}{dt} + \delta_LK_Le^{-\delta_LM}\frac{dM}{dt}\Big]\\ & = 0. \end{align*} Grouping like terms (by state derivative) and substituting equation \eqref{Leqn} results in \begin{align*} \frac{d^2\lambda_3}{dt^2} & = \frac{d\lambda_1}{dt}\frac{\partial D}{\partial L}T -\frac{d\lambda_6}{dt}\left(\frac{\omega I}{\zeta+I}\right) \\ &\quad + \left[ \lambda_1 T\frac{\partial^{2}D}{\partial T\partial L} + \lambda_1\frac{\partial D}{\partial L} +\lambda_3\left(q-\frac{jk}{(k+T)^2}\right)\right]\frac{dT}{dt}\\ &\quad +\left[ \lambda_3\left(-\frac{\theta m} {(\theta+I)^2} - \frac{p_Ig_I}{(g_I+I)^2} +2u_0LC\frac{\kappa}{(\kappa+I)^2}\right)-\lambda_6\frac{\omega\zeta}{(\zeta+I)^2}\right]\frac{dI}{dt}\\ &\quad +2 u_0\lambda_3L\left( \frac{I}{\kappa+I}\right)\frac{dC}{dt} +\lambda_3 \delta_LK_Le^{-\delta_LM}\frac{dM}{dt}\\ &\quad + \left[\lambda_1 T\frac{\partial ^2D}{\partial{L^2}}+ 2u_0\lambda_3C\left(\frac{I}{\kappa+I}\right) \right]\\ &\quad\times\Big[-\frac{\theta mL}{\theta+I} - qLT + r_1NT + r_2CT + \frac{p_I L I}{g_I + I}-\frac{u_0L^2CI}{\kappa+I} \\ & \quad + \frac{jTL}{k+T} -K_L(1-e^{-\delta_L M})L+\eta_1 v_L(t)\Big]\\ & = 0. \end{align*} Isolating the $v_L$ term and substituting $\lambda_3=-\frac{\epsilon_L}{\eta_1}$ gives \begin{align*} \frac{d^2\lambda_3}{dt^2} & = \eta_1\left[\lambda_1 T\frac{\partial ^2D}{\partial{L^2}}- 2u_0C\frac{\epsilon_L}{\eta_1}\left(\frac{I}{\kappa+I}\right)\right]v_L \\ &\quad+\frac{d\lambda_1}{dt}\frac{\partial D}{\partial L}T -\frac{d\lambda_6}{dt}\left(\frac{\omega I}{\zeta+I}\right) \\ &\quad + \left[ \lambda_1 T\frac{\partial^{2}D}{\partial T\partial L} + \lambda_1\frac{\partial D}{\partial L} -\frac{\epsilon_L}{\eta_1}\left(q-\frac{jk}{(k+T)^2}\right)\right]\frac{dT}{dt}\\ &\quad +\left[ \frac{\epsilon_L}{\eta_1}\left(\frac{\theta m} {(\theta+I)^2}+ \frac{p_Ig_I}{(g_I+I)^2} -2u_0LC\frac{\kappa}{(\kappa+I)^2}\right)-\lambda_6\frac{\omega\zeta}{(\zeta+I)^2}\right]\frac{dI}{dt}\\ &\quad -2 u_0L\frac{\epsilon_L}{\eta_1}\left( \frac{I}{\kappa+I}\right)\frac{dC}{dt} - \delta_LK_Le^{-\delta_LM}\frac{\epsilon_L}{\eta_1}\frac{dM}{dt}\\ &\quad + \left[\lambda_1 T\frac{\partial ^2D}{\partial{L^2}}- 2u_0C\frac{\epsilon_L}{\eta_1} \left(\frac{I}{\kappa+I}\right) \right]\\ &\quad\times\Big[-\frac{\theta mL}{\theta+I} - qLT + r_1NT + r_2CT + \frac{p_I L I}{g_I + I}-\frac{u_0L^2CI}{\kappa+I} \\ & \quad + \frac{jTL}{k+T} -K_L(1-e^{-\delta_L M})L\Big]\\ & = 0. \end{align*} Now, let $P_L$ denote the above sum \emph{excluding} the $v_L$ term. Notice that $P_L$ depends on $\lambda_1,\lambda_2,\lambda_6, T,N,L,C,M,I,v_M,\text{and }v_I$. Next, define $Q_L$ as \begin{equation*} Q_L = \eta_1\left[\lambda_1 T\frac{\partial ^2D}{\partial{L^2}}- 2u_0C\frac{\epsilon_L}{\eta_1} \big(\frac{I}{\kappa+I}\big) \right]. \end{equation*} Using this condensed notation, we see that \begin{equation*} v_L = -\frac{P_L}{Q_L}, \quad\text{if } Q_L\neq 0. \end{equation*} Here the degree of the singularity is 2.\\ Furthermore, \emph{if} $v_I$ and $v_M$ are assumed to be bang-bang on the interval $(s_L,t_L)$, then the Legendre-Clebsch condition for the control to be minimizing is \begin{equation*} (-1)\frac{\partial}{\partial v_L}\frac{d^2}{dt^2}\frac{\partial H}{\partial v_L}\ge0, \quad\text{or}\quad Q_L \le0; \end{equation*} i.e. \begin{equation*} \lambda_1T\frac{\partial^2 D}{\partial L^2}\le 2u_0C\frac{\epsilon_L}{\eta_1}\left(\frac{I}{\kappa+I}\right). \end{equation*} \end{proof} \begin{theorem} \label{thm3.4} If $v_M$ is singular on the interval $(s_M, t_M)$, the singularity is of degree two and the representation of the control is given by \begin{equation*} v_M = -\frac{P_M}{Q_M}, \quad Q_M\neq 0, \end{equation*} where $Q_M$ is given by \begin{align} Q_M &= -\eta_2\Big[\delta_TK_T\lambda_1T(\delta_Te^{-\delta_TM})+ \delta_NK_N\lambda_2N(-\delta_Ne^{\delta_NM}) \notag \\ &\quad+\delta_LK_L\lambda_3L(\delta_Le^{-\delta_LM})+ \delta_CK_C\lambda_4C(-\delta_Ce^{\delta_CM})\Big]\label{Q_M} \end{align} and $P_M$ is given in Appendix B.\\ Furthermore, if $v_M$ is minimizing on this interval, it is necessary that \begin{align*} &\Big[\delta_TK_T\lambda_1T(\delta_Te^{-\delta_TM}) +\delta_NK_N\lambda_2N(\delta_Ne^{-\delta_NM})\\ &+\delta_LK_L\lambda_3L(\delta_Le^{-\delta_LM}) +\delta_CK_C\lambda_4C(\delta_Ce^{-\delta_CM})\Big]\ge0. \end{align*} \end{theorem} \begin{proof} If $v_M$ is singular on some interval $(s_M,t_M)$, we know that time derivatives of $\Phi_M$ are identically zero in this region; in addition, we have that $\lambda_5 = -\epsilon_M/\eta_2$. As before, note that \begin{equation*} \frac{d}{dt}\Phi_M = \eta_2\frac{d}{dt}\lambda_5 + \frac{d}{dt}\lambda_7 =0, \quad\text{or } \frac{d\lambda_5}{dt} =0, \end{equation*} since $\frac{d\lambda_7}{dt}=0$. Coupling this information with equation \eqref{adjM} yields \begin{align*} \frac{d\lambda_5}{dt} &= \lambda_1\delta_TK_TTe^{-\delta_TM}+\lambda_2\delta_NK_NNe^{-\delta_NM} +\lambda_3\delta_LK_LLe^{-\delta_LM} \\ &\quad +\lambda_4\delta_CK_CCe^{-\delta_CM}+\gamma\lambda_5\\ &=0. \end{align*} Taking a time derivative gives the result \begin{align*} \frac{d^2\lambda_5}{dt^2} &= \delta_TK_T\left[\lambda_1T(-\delta_Te^{-\delta_TM})\frac{dM}{dt} + e^{-\delta_TM}\left(\lambda_1\frac{dT}{dt} + T\frac{d\lambda_1}{dt}\right)\right]\\ &\quad +\delta_NK_N\left[\lambda_2N(-\delta_Ne^{-\delta_NM})\frac{dM}{dt} + e^{-\delta_NM}\left(\lambda_2\frac{dN}{dt} + N\frac{d\lambda_2}{dt}\right)\right]\\ &\quad +\delta_LK_L\left[\lambda_3L(-\delta_Le^{-\delta_LM})\frac{dM}{dt} + e^{-\delta_LM}\left(\lambda_3\frac{dL}{dt} + L\frac{d\lambda_3}{dt}\right)\right]\\ &\quad +\delta_CK_C\left[\lambda_4C(-\delta_Ce^{-\delta_CM})\frac{dM}{dt} + e^{-\delta_CM}\left(\lambda_4\frac{dC}{dt} + C\frac{d\lambda_4}{dt}\right)\right]\\ &\quad +\gamma\frac{d\lambda_5}{dt}\\ &=0. \end{align*} Notice that the last term is zero, since $\frac{d\lambda_5}{dt} =0$. Grouping the $\frac{dM}{dt}$terms, substituting equation \eqref{Meqn}, and isolating the $v_M$ term in the above equation gives \begin{align*} \frac{d^2\lambda_5}{dt^2}&= \eta_2\Big[\delta_TK_T\lambda_1T(-\delta_Te^{-\delta_TM})+ \delta_NK_N\lambda_2N(-\delta_Ne^{-\delta_NM})\\ &\quad+\delta_LK_L\lambda_3L(-\delta_Le^{-\delta_LM})+ \delta_CK_C\lambda_4C(-\delta_Ce^{-\delta_CM})\Big]v_M\\ &\quad + \delta_TK_Te^{-\delta_TM}\left(\lambda_1\frac{dT}{dt}+T\frac{d\lambda_1}{dt}\right) +\delta_NK_Ne^{-\delta_NM}\left(\lambda_2\frac{dN}{dt}+N\frac{d\lambda_2}{dt}\right)\\ &\quad +\delta_LK_Le^{-\delta_LM}\left(\lambda_3\frac{dL}{dt}+L\frac{d\lambda_3}{dt}\right) +\delta_CK_Ce^{-\delta_CM}\left(\lambda_4\frac{dC}{dt}+C\frac{d\lambda_4}{dt}\right)\\ &\quad -\gamma M\Big[\delta_TK_T\lambda_1T(-\delta_Te^{-\delta_TM})+ \delta_NK_N\lambda_2N(-\delta_Ne^{-\delta_NM})\\ &\quad+\delta_LK_L\lambda_3L(-\delta_Le^{-\delta_LM})+ \delta_CK_C\lambda_4C(-\delta_Ce^{-\delta_CM})\Big]\\ &=0. \end{align*} Now, as before, let $P_M$ denote the above sum, excluding the $v_M$ term. Note that $P_M$ is dependent on all the state and adjoint variables, as well as the control $v_L$. Define $Q_M$ as \begin{align*} Q_M &= -\eta_2\Big[\delta_TK_T\lambda_1T(\delta_Te^{-\delta_TM})+ \delta_NK_N\lambda_2N(-\delta_Ne^{\delta_NM})\\ &\quad+\delta_LK_L\lambda_3L(\delta_Le^{-\delta_LM})+ \delta_CK_C\lambda_4C(-\delta_Ce^{\delta_CM})\Big]. \end{align*} Then we have \begin{equation*} v_M = -\frac{P_M}{Q_M}, \quad\text{if } Q_M\neq 0. \end{equation*} Here the degree of the singularity is two.\\ Furthermore, \emph{if} $v_L$ and $v_I$ are assumed to be bang-bang on the interval $(s_M,t_M)$, then the Legendre-Clebsch condition for the control to be minimizing is \begin{equation*} (-1)\frac{\partial}{\partial v_M}\frac{d^2}{dt^2}\frac{\partial H}{\partial v_M}\ge0, \end{equation*} or in this situation, \begin{equation*} Q_M \le0. \end{equation*} Thus, if $v_M$ is singular on some interval, we have the additional necessary condition \begin{align*} &\Big[\delta_TK_T\lambda_1T(\delta_Te^{-\delta_TM}) +\delta_NK_N\lambda_2N(\delta_Ne^{-\delta_NM})\\ &+\delta_LK_L\lambda_3L(\delta_Le^{-\delta_LM}) +\delta_CK_C\lambda_4C(\delta_Ce^{-\delta_CM})\Big]\ge0, \end{align*} since $\eta_2\ne0$, in general. \end{proof} \begin{theorem}\label{thm3.5} If $v_I$ is singular on the interval $(s_I, t_I)$, the singularity is of degree two and the representation of the control is given by \begin{equation*} v_I = -\frac{P_I}{Q_I}, \quad Q_I\neq 0, \end{equation*} where $Q_I$ is given by \begin{equation} Q_I=\eta_3\left[\frac{2\lambda_2p_Ng_NN}{(g_N+I)^3}+\frac{2\lambda_3\theta mL}{(\theta+I)^3}+\frac{2\lambda_3p_Ig_IL}{(g_I+I)^3}-\frac{2\lambda_3u_0\kappa L^2C}{(\kappa+I)^3}-\frac{2\epsilon_I\omega\zeta L}{\eta_3(\zeta+I)^3}\right] \label{Q_I} \end{equation} and $P_I$ is given in Appendix B. Furthermore, if $v_I$ is minimizing on this interval, it is necessary that \begin{equation*} \frac{\lambda_2p_Ng_NN}{(g_N+I)^3}+\frac{\lambda_3\theta mL}{(\theta+I)^3}+\frac{\lambda_3p_Ig_IL}{(g_I+I)^3}\le\frac{\lambda_3u_0\kappa L^2C}{(\kappa+I)^3}+\frac{\epsilon_I\omega\zeta L}{\eta_3(\zeta+I)^3}\text{ .} \end{equation*} \end{theorem} \begin{proof} Assume $v_I$ is singular on the interval $(s_I,t_I)$; then all derivatives with respect to time of $\Phi_I$ are equal to zero and $\lambda_6 = -\frac{\epsilon_I}{\eta_3}$. Then \begin{equation*} \frac{d}{dt}\Phi_I = \eta_3\frac{d}{dt}\lambda_6 =0,\quad\text{and so}\quad\frac{d\lambda_6}{dt} = 0. \end{equation*} Together with equation \eqref{adjI}, we have \begin{align*} \frac{d\lambda_6}{dt} &= -\lambda_2\left(\frac{p_Ng_NN}{(g_N+I)^2}\right) -\lambda_3\left(\frac{\theta mL}{(\theta+I)^2}+\frac{p_Ig_IL}{(g_I+I)^2}-\frac{u_0\kappa L^2C}{(\kappa+I)^2}\right)\\ &\quad +\lambda_6\left(\mu_I-\frac{\omega\zeta L}{(\zeta+I)^2}\right)\\ &=0. \end{align*} Differentiate with respect to $t$ for \begin{align} \frac{d^2\lambda_6}{dt^2} &= -\frac{d\lambda_2}{dt}\left(\frac{p_Ng_NN}{(g_N+I)^2}\right) -\lambda_2\Big[\frac{(g_N+I)^2p_Ng_N\frac{dN}{dt}-p_Ng_NN\cdot 2(g_N+I) \frac{dI}{dt}}{(g_N+I)^4}\Big]\notag\\ &\quad -\frac{d\lambda_3}{dt}\left[\frac{\theta mL}{(\theta+I)^2} +\frac{p_Ig_IL}{(g_I+I)^2}-\frac{u_0\kappa L^2C}{(\kappa+I)^2}\right]\notag\\ &\quad -\lambda_3\Bigg[\frac{(\theta+I)^2\theta m\frac{dL}{dt} -\theta mL\cdot 2(\theta+I)\frac{dI}{dt}}{(\theta+I)^4}\notag\\ &\quad +\frac{(g_I+I)^2p_Ig_I\frac{dL}{dt}-p_Ig_IL\cdot 2(g_I+I) \frac{dI}{dt}}{(g_I+I)^4}\notag\\ &\quad -\frac{u_0\kappa (\kappa+I)^2\left[L^2\frac{dC}{dt}+2LC\frac{dL}{dt} \right]-u_0\kappa L^2C\cdot 2(\kappa+I)\frac{dI}{dt}}{(\kappa+I)^4}\Bigg] \notag\\ &\quad +\frac{d\lambda_6}{dt}\left[\mu_I-\frac{\omega\zeta L}{(\zeta+I)^2} \right] -\lambda_6\Big[\frac{(\zeta+I)^2\omega\zeta\frac{dL}{dt} -\omega\zeta L\cdot 2(\zeta+I)\frac{dI}{dt}}{(\zeta+I)^4}\Big]\notag \\ &=0 . \label{vising} \end{align} We substitute $\lambda_6 = -\epsilon_I/\eta_3$, and note that the $\frac{d\lambda_6}{dt}$ term equals zero. Then we group terms, substitute equation \eqref{Ieqn}, and isolate the $v_I$ term to get \begin{align*} 0&=\frac{d^2\lambda_6}{dt^2} \\ &=\eta_3\left[\frac{2\lambda_2p_Ng_NN}{(g_N+I)^3}+\frac{2\lambda_3\theta mL}{(\theta+I)^3}+\frac{2\lambda_3p_Ig_IL}{(g_I+I)^3}-\frac{2\lambda_3u_0\kappa L^2C}{(\kappa+I)^3}-\frac{2\epsilon_I\omega\zeta L}{\eta_3(\zeta+I)^3}\right]v_I\\ &\quad -\frac{d\lambda_2}{dt}\left(\frac{p_Ng_NN}{(g_N+I)^2}\right)\\ &\quad -\frac{d\lambda_3}{dt}\left[\frac{\theta mL}{(\theta+I)^2}+\frac{p_Ig_IL}{(g_I+I)^2}-\frac{u_0\kappa L^2C}{(\kappa+I)^2}\right]\\ &\quad -\frac{dN}{dt}\left(\frac{\lambda_2p_Ng_N}{(g_N+I)^2}\right) -\frac{dC}{dt}\left(\frac{\lambda_3 u_0\kappa L^2}{(\kappa+I)^2}\right)\\ &\quad +\frac{dL}{dt}\left[-\frac{\lambda_3\theta m}{(\theta+I)^2}-\frac{\lambda_3p_Ig_I}{(g_I+I)^2}+\frac{2\lambda_3u_0\kappa LC}{(\kappa+I)^2}+\frac{\epsilon_I\omega\zeta}{\eta_3(\zeta+I)^2}\right]\\ &\quad +\left[\frac{2\lambda_2p_Ng_NN}{(g_N+I)^3} +\frac{2\lambda_3\theta mL}{(\theta+I)^3} +\frac{2\lambda_3p_Ig_IL}{(g_I+I)^3} -\frac{2\lambda_3u_0\kappa L^2C}{(\kappa+I)^3} -\frac{2\epsilon_I\omega\zeta L}{\eta_3(\zeta+I)^3}\right]\\ &\quad\times \Big(-\mu_I I + \phi C +\frac{\omega LI}{\zeta+I}\Big). \end{align*} Let $P_I$ denote the above sum excluding the $v_I$ term. Note that $P_I$ depends on all the state and adjoint variables, as well as the control $v_L$. Let $Q_I$ be defined by \[ Q_I=\eta_3\left[\frac{2\lambda_2p_Ng_NN}{(g_N+I)^3}+\frac{2\lambda_3\theta mL}{(\theta+I)^3}+\frac{2\lambda_3p_Ig_IL}{(g_I+I)^3}-\frac{2\lambda_3u_0\kappa L^2C}{(\kappa+I)^3}-\frac{2\epsilon_I\omega\zeta L}{\eta_3(\zeta+I)^3}\right]. \] Then we have \begin{equation*} v_I = -\frac{P_I}{Q_I}, \quad\text{if } Q_I\neq 0. \end{equation*} Here the degree of the singularity is 2. Furthermore, \emph{if} $v_L$ and $v_M$ are assumed to be bang-bang on the interval $(s_I,t_I)$, then the Legendre-Clebsch condition for the control to be minimizing is \begin{equation*} (-1)\frac{\partial}{\partial v_I}\frac{d^2}{dt^2}\frac{\partial H}{\partial v_I}\ge0, \end{equation*} or in this situation, $Q_I \le0$. Thus, if $v_I$ is singular on some interval, we have the additional necessary condition \begin{equation*} \frac{\lambda_2p_Ng_NN}{(g_N+I)^3}+\frac{\lambda_3\theta mL}{(\theta+I)^3}+\frac{\lambda_3p_Ig_IL}{(g_I+I)^3}\le\frac{\lambda_3u_0\kappa L^2C}{(\kappa+I)^3}+\frac{\epsilon_I\omega\zeta L}{\eta_3(\zeta+I)^3}, \end{equation*} since $\eta_3 \neq0$. \end{proof} Having found the representations for each control on a singular interval, we turn our attention to the necessary conditions generated when two of the controls are simultaneously singular on the same interval. Note that the degree of singularity of each control is two; both this fact and the representations found previously will be used to generate the necessary Legendre-Clebsch conditions. \begin{theorem} \label{thm3.6} If $v_L$ and $v_M$ are both singular on some interval $(s,t)$, then the following conditions must hold for $v_L$ and $v_M$ to be minimizing: \begin{gather*} Q_L \le0,\\ Q_LQ_M\ge \eta_1\eta_2(\lambda_3\delta_LK_Le^{-\delta_LM})^2, \end{gather*} where $Q_L$ and $Q_M$ are defined as in \eqref{Q_L} and \eqref{Q_M}, respectively. \end{theorem} \begin{proof} Assuming that $v_L$ and $v_M$ are both singular on the interval $(s,t)$, we use the generalized Legendre-Clebsch conditions to form the $2\times 2$ matrix \[ A_{ij} = (-1)^{\frac{h_j +1}{2}}\frac{\partial}{\partial v_i}\frac{d^{({\frac{h_i +h_j}{2} +1})}}{dt^{({\frac{h_i +h_j}{2}+1})}} \frac{\partial H}{\partial v_j} = (-1)\frac{\partial}{\partial v_L}\frac{d^{2}}{dt^{2}}\frac{\partial H}{\partial v_M}, \] which is given by \begin{equation*} A_{v_L,v_M}=\begin{pmatrix} -\eta_1Q_L & -\eta_1\eta_2\lambda_3\delta_LK_Le^{-\delta_LM}\\ -\eta_1\eta_2\lambda_3\delta_LK_Le^{-\delta_LM} & -\eta_2Q_M \end{pmatrix}, \end{equation*} Note that $Q_L$ and $Q_M$ are defined as in \eqref{Q_L} and \eqref{Q_M}, respectively. This matrix is clearly symmetric. In order for $v_L$ and $v_M$ to be minimizing on the interval $(s,t)$, the matrix $A_{v_L,v_M}$ must be positive definite. In other words, all upper left minors of $A_{v_L,v_M}$ must be positive; then the conditions \begin{gather*} Q_L \le0,\\ Q_LQ_M\ge \eta_1\eta_2(\lambda_3\delta_LK_Le^{-\delta_LM})^2. \end{gather*} must hold. \end{proof} In a similar fashion, conditions for the pairs $v_L$, $v_I$ and $v_M$, $v_I$ to be singular on a given interval are found as noted in the following theorems without proof. \begin{theorem} \label{thm3.7} If $v_L$ and $v_I$ are both singular on some interval $(s,t)$, then the following conditions must hold for $v_L$ and $v_I$ to be minimizing: \begin{gather*} Q_L \le0,\\ Q_LQ_I\ge \eta_1\eta_3(\circledast)^2. \end{gather*} where $Q_L$ and $Q_I$ are defined as in \eqref{Q_L} and \eqref{Q_I}, respectively, and \begin{equation} \circledast = \left(\frac{\lambda_3\theta{m}}{(\theta+I)^2}+\frac{\lambda_3p_Ig_I}{(g_I+I)^2}-\frac{2\lambda_3u_0\kappa LC}{(\kappa+I)^2}+\frac{\lambda_6\omega\zeta}{(\zeta+I)^2}\right).\label{v_LIduo} \end{equation} \end{theorem} \begin{theorem} \label{thm3.8} If $v_M$ and $v_I$ are both singular on some interval $(s,t)$, then the following conditions must hold for $v_M$ and $v_I$ to be minimizing: \begin{gather*} Q_M \le0,\\ Q_MQ_I\ge 0. \end{gather*} where $Q_M$ and $Q_I$ are defined as in \eqref{Q_M} and \eqref{Q_I}, respectively. \end{theorem} Finally, we will use the Legendre-Clebsch conditions to find the conditions necessary for all three controls to be minimizing on the same interval. \begin{theorem} \label{thm3.9} If $v_L$, $v_M$, and $v_I$ are simultaneously singular on some interval $(s,t)$, then the following conditions must hold for these controls to be minimizing: \begin{gather*} Q_L \le0,\\ Q_LQ_M\ge \eta_1\eta_2(\circleddash)^2,\\ Q_LQ_MQ_I \le \eta_1\eta_3Q_M\circledast + \eta_1\eta_2Q_I\circleddash, \end{gather*} where $Q_L$, $Q_M$, $Q_I$ and $\circledast$ are defined by \eqref{Q_L}, \eqref{Q_M}, \eqref{Q_I}, and \eqref{v_LIduo}, respectively, and \begin{equation} \circleddash = -\lambda_3\delta_LK_Le^{-\delta_LM}. \label{v_LMduo} \end{equation} \end{theorem} \begin{proof} Assuming $v_L$, $v_M$, and $v_I$ are all singular on the interval $(s,t)$, we use the generalized Legendre-Clebsch conditions to form the $3\times 3$ matrix whose entries are given by \begin{equation*} A_{ij} = (-1)^{\frac{h_j +1}{2}}\frac{\partial}{\partial v_i}\frac{d^{({\frac{h_i +h_j}{2} +1})}}{dt^{({\frac{h_i +h_j}{2} +1})}}\frac{\partial H}{\partial v_j} . \end{equation*} Here we have \begin{equation*} A_{v_L,v_M,v_I}=\begin{pmatrix} -\eta_1Q_L & \eta_1\eta_2\circleddash & \eta_1\eta_3\circledast \\ \eta_1\eta_2\circleddash & -\eta_2Q_M & 0\\ \eta_1\eta_3\circledast & 0 & -\eta_3Q_I \end{pmatrix}, \end{equation*} where $Q_L$, $Q_M$, $Q_I$, $\circledast$, and $\circleddash$ are defined previously in \eqref{v_LIduo} and \eqref{v_LMduo}. This matrix is symmetric; in order for $A_{v_L,v_M,v_I}$ to be positive semidefinite, we must have that \begin{gather*} Q_L \le0,\\ Q_LQ_M\ge \eta_1\eta_2(\circleddash)^2,\\ Q_LQ_MQ_I \le \eta_1\eta_3Q_M\circledast \, + \,\, \eta_1\eta_2Q_I\circleddash. \end{gather*} \end{proof} \subsection*{Conclusion} If the controls are singular either as a unit, in pairs, or as a triple, conditions are given for those controls to be minimizing. First, existence of the control triple is established. Within the characterization, the introduction of singular and/or bang-bang controls is given. It is worth noting that clinicians commonly give treatments in a bang-bang type scenario, i.e. give the drug combination for a period and then do not given any drug, c.f. \cite{Hermans2003}, \cite{Rosti1998}. However, in this work, conditions are given such that a singular control vector would minimize the proposed objective functional. This could change the dynamic in which the chemotherapy and immunotherapy treatments are administered. Due to the interdependence on the conditions, we note that the characterizations of the singular control on the selected intervals are not explicitly determined. In future work, numerical analyses will be investigated to aid in graphically characterizing these singular control expressions. \section{Appendix A} In this section, we will give precise statements of the theorems used for proving the existence and finding the characterizations of optimal quadratic controls. \begin{theorem}[Fillipov-Cesari Theorem \cite{HSV1995}]\label{eqn:FCT} Consider the following optimal control problem: \begin{align*} \min &\quad J = \int_0^{T} F(x(t),u(t),t)dt + S(x(T),T) \\ &\quad \dot{x}(t) = f(x(t),u(t),t), \quad x(0)=x_0 \\ &\quad g(x(t),u(t),t)\ge0\\ &\quad h(x(t),t)\ge0\\ &\quad a(x(T),T)\ge0\\ &\quad b(x(T),T)=0 \end{align*} where T is free on $[0,t_f]$. Assume that F,f,g,h,S,a, and b are continuous in all their arguments at all points $(x,u,t)$. Define the (state-dependent) control region \begin{equation*} \Omega (x,t) = \{u\in \mathbb{R}^{m}|\;g(x,u,t)\ge0\}\subset \mathbb{R}^{m} \end{equation*} and the set \begin{equation*} N(x,t)=\{(F(x,u,t)+\gamma,f(x,u,t))|\;\gamma\le0,u\in\Omega(x,t)\}\subset\mathbb{R}^{n+1} \end{equation*} where $m$ and $n$ are the number of control and state variables, respectively. Suppose that the following conditions hold: \begin{enumerate} \item There exits an admissible solution pair. \item $N(x,t)$ is convex for all $(x,t)\in\mathbb{R}^{n}\times[0,t_f]$. \item There exists $\delta > 0$ such that $\| x(t) \| <\delta$ for all admissible $\{x(t),u(t)\}$ and $t$. \item There exists $\delta_1 >0$ such that $\| u \| <\delta_1$ for all $u\in\Omega (x,t)$ with $\| x \| <\delta$. \end{enumerate} Then there exists an optimal triple $\{T^*,x^*,u^*\}$ with $u^{*}(\cdot)$ measurable. \end{theorem} \begin{theorem}[Pontraygin's Maximum/Minimum Principle \cite{KS1991}] \label{eqn:PMP} % \label{PMP} \quad \\ Let $\textbf{u}(t)=[u_1(t),\ldots,u_m(t)]$ be a piecewise continuous control vector and $\textbf{x}(t)=[x_1(t),\ldots,x_n(t)]$ be an associated continuous and piecewise differentiable state vector defined on the fixed time interval $[t_0,t_1]$ that minimizes \begin{equation*} \int_{t_0}^{t_1} f(t,\textbf{x}(t),\textbf{u}(t))dt \end{equation*} subject to the differential equations \begin{equation*} x_i(t)=g_i(t,\textbf{x}(t),\textbf{u}(t)),\quad i=1,\ldots,n, \end{equation*} initial conditions \begin{equation*} x_i(t_0)=x_{i0},\quad i=1,\ldots,n\quad (x_{i0}\text{ fixed}), \end{equation*} terminal conditions \begin{gather*} x_i(t_1) = x_{i1},\quad i=1,\ldots,p,\\ x_i(t_1) \ge x_{it},\quad i=p+1,\ldots,q\quad (x_{i1}, i=1,\ldots,q\text{ fixed}),\\ x_i(t_1) \text{free},\quad i=q+1,\ldots,n, \end{gather*} and control variable restriction \begin{equation*} \textbf{u}(t)\in U,\quad U\text{ a given set in }\mathbb{R}^m. \end{equation*} We assume that $f,g,\partial f / \partial x_j$, and $\partial g_i / \partial x_j$ are continuous functions of all their arguments, for all $i=1,\ldots,n$ and $j=1,\ldots,n$. Then there exists a constant $\lambda_0$ and continuous functions $\lambda(t)=(\lambda_1(t),\ldots,\lambda_n(t))$, where for all $t_0\le t\le t_1$ we have $(\lambda_0,\lambda(t))\neq (0,0)$ such that for every $t_0\le t\le t_1$, \begin{equation*} H(t,\textbf{x}^*(t),\textbf{u}(t),\lambda(t))\le H(t,\textbf{x}^*(t),\textbf{u}^*(t),\lambda(t)), \end{equation*} where the Hamiltonian function H is defined by \begin{equation*} H(t,\textbf{x},\textbf{u},\lambda)=\lambda_0 f(t,\textbf{x},\textbf{u}) + \sum_{i=1}^{n} \lambda_i g_i (t,\textbf{x},\textbf{u}). \end{equation*} Except at points of discontinuity of $\textbf{u}^*(t)$, \begin{equation*} \lambda ' (t)= -\partial H(t,\textbf{x}^*(t),\textbf{u}^*(t),\lambda(t)) / \partial x_i ,\quad i=1,\ldots,n. \end{equation*} Finally, the following transversality conditions are satisfied: \begin{gather*} \lambda_i(t_1) \quad \text{ no conditions},\quad i=1,\ldots,p,\\ \lambda_i(t_1) \ge0 \quad (=0 \text{ if } x_{i}^*(t_1)>x_{i1}))\quad i=p+1,\ldots,q,\\ \lambda_i(t_1) = 0, \quad i=q+1,\ldots,n.\\ \end{gather*} \end{theorem} In addition, the modifications to \eqref{eqn:PMP} generated by the terminal inequality are given in Kamien and Schwartz \cite{KS1991}, p. 160: \emph{If $K(x_q(t_1),\ldots,x_n(t_1))\ge0$ is required, then the transversality conditions \begin{align*} \lambda_i(t_1) &=\text{p } \partial K /\partial x_1,\quad i=q,\ldots,n,\\ p &\le0,\\ pK &=0 \end{align*} are necessary.} \begin{theorem}[Generalized Legendre-Clebsch Conditions \cite{Krener1977}]\label{eqn:GLC} Assume that $\mathbf{u}(t)$ and $\mathbf{x}(t)$ are defined for \eqref{Teqn}-\eqref{Ieqn} on $[0,t_f]$. Suppose $\mathbf{u}(t) \in \emph{interior }\Omega$ and each $\mathbf{u_i}$ is singular of degree $\frac{h_i+1}{2}$ on the subinterval $(t_0,t_1)$. If $\mathbf{u}(t)$ is minimal, then there exists a $\mathbf{\lambda}(t)$ satisfying the PMP on $[0,t_f]$ such that on the subinterval $(t_0,t_1)$, \begin{equation*} \frac{\partial}{\partial\mathbf{u_i}}\frac{d^k}{dt^k}\frac{\partial}{\partial\mathbf{u_j}}H(\mathbf{\lambda}(t),\mathbf{x}(t),\mathbf{u}(t))=0 \end{equation*} for $k=0,\dots, \frac{h_i+h_j}{2}, 1\le i, j\le l$ (where $l$ is the dimension of the control space). Moreover, if $h_i <\infty$ for $i = 1,\dots,k\le l$, then the $k\times k$ matrix whose $i,j$ entry is \begin{equation*} (-1)^{\frac{h_j+1}{2}}\frac{\partial}{\partial\mathbf{u_i}} \frac{d^{({\frac{h_i +h_j}{2} +1)}}}{dt^{({\frac{h_i +h_j}{2}+1)}}} \frac{\partial}{\partial\mathbf{u_j}}H(\mathbf{\lambda}(t), \mathbf{x}(t),\mathbf{u}(t)) \end{equation*} where $1\le i, j\le k$, must be symmetric and nonnegative definite. \end{theorem} Note that a given symmetric real matrix is nonnegative definite (i.e. positive semidefinite) if and only if all of its upper-left minors are positive; that is, for a given $n\times n$ matrix \begin{equation*} A_{n\times n} =\begin{pmatrix} A_{11} &\dots & A_{1n}\\ \vdots &\ddots & \vdots\\ A_{n1} &\dots & A_{nn} \end{pmatrix}, \end{equation*} we have that \[ A_{11}\ge0,\quad \left|\begin{matrix} A_{11} & A_{12}\\ A_{21} & A_{22} \end{matrix} \right| \ge0,\quad \left|\begin{matrix} A_{11} & A_{12} & A_{13}\\ A_{21} & A_{22} & A_{23}\\ A_{31} & A_{32} & A_{33} \end{matrix} \right| \ge0, \] etc. \section{Appendix B} In this section, we will state the equations $P_L, P_M$, and $P_I$ used in the characterizations of controls $v_L, v_M$, and $v_I$, respectively. \begin{align*} P_L &= \frac{d\lambda_1}{dt}\frac{\partial D}{\partial L}T -\frac{d\lambda_6}{dt}\left(\frac{\omega I}{\zeta+I}\right) \\ &\quad + \left[ \lambda_1 T\frac{\partial^{2}D}{\partial T\partial L} + \lambda_1\frac{\partial D}{\partial L} -\frac{\epsilon_L}{\eta_1}\left(q-\frac{jk}{(k+T)^2}\right)\right]\frac{dT}{dt}\\ &\quad +\left[ \frac{\epsilon_L}{\eta_1}\left(\frac{\theta m} {(\theta+I)^2}+ \frac{p_Ig_I}{(g_I+I)^2} -2u_0LC\frac{\kappa}{(\kappa+I)^2}\right)-\lambda_6\frac{\omega\zeta}{(\zeta+I)^2}\right]\frac{dI}{dt}\\ &\quad -2 u_0L\frac{\epsilon_L}{\eta_1}\left(\frac{I}{\kappa+I}\right)\frac{dC}{dt} - \delta_LK_Le^{-\delta_LM}\frac{\epsilon_L}{\eta_1}\frac{dM}{dt}\\ &\quad + \left[\lambda_1 T\frac{\partial ^2D}{\partial{L^2}}- 2u_0C\frac{\epsilon_L}{\eta_1} \left(\frac{I}{\kappa+I}\right) \right]\\ &\quad\times \Big[-\frac{\theta mL}{\theta+I} - qLT + r_1NT + r_2CT + \frac{p_I L I}{g_I + I}-\frac{u_0L^2CI}{\kappa+I} \\ &\quad + \frac{jTL}{k+T} -K_L(1-e^{-\delta_L M})L\Big]\,, \end{align*} \begin{align*} P_M &=\delta_TK_Te^{-\delta_TM}\left(\lambda_1\frac{dT}{dt}+T\frac{d\lambda_1}{dt}\right) +\delta_NK_Ne^{-\delta_NM}\left(\lambda_2\frac{dN}{dt}+N\frac{d\lambda_2}{dt}\right)\\ &\quad +\delta_LK_Le^{-\delta_LM}\left(\lambda_3\frac{dL}{dt}+L\frac{d\lambda_3}{dt}\right) +\delta_CK_Ce^{-\delta_CM}\left(\lambda_4\frac{dC}{dt}+C\frac{d\lambda_4}{dt}\right)\\ &\quad -\gamma M\Big[\delta_TK_T\lambda_1T(-\delta_Te^{-\delta_TM}) +\delta_NK_N\lambda_2N(-\delta_Ne^{-\delta_NM})\\ &\quad+\delta_LK_L\lambda_3L(-\delta_Le^{-\delta_LM})+ \delta_CK_C\lambda_4C(-\delta_Ce^{-\delta_CM})\Big] \end{align*}\,, \begin{align*} P_I &=-\frac{d\lambda_2}{dt}\left(\frac{p_Ng_NN}{(g_N+I)^2}\right) -\frac{d\lambda_3}{dt}\left[\frac{\theta mL}{(\theta+I)^2}+\frac{p_Ig_IL}{(g_I+I)^2}-\frac{u_0\kappa L^2C}{(\kappa+I)^2}\right]\\ &\quad -\frac{dN}{dt}\left(\frac{\lambda_2p_Ng_N}{(g_N+I)^2}\right) -\frac{dC}{dt}\left(\frac{\lambda_3 u_0\kappa L^2}{(\kappa+I)^2}\right)\\ &\quad +\frac{dL}{dt}\left[-\frac{\lambda_3\theta m}{(\theta+I)^2}-\frac{\lambda_3p_Ig_I}{(g_I+I)^2}+\frac{2\lambda_3u_0\kappa LC}{(\kappa+I)^2}+\frac{\epsilon_I\omega\zeta}{\eta_3(\zeta+I)^2}\right]\\ &\quad +\left[\frac{2\lambda_2p_Ng_NN}{(g_N+I)^3} +\frac{2\lambda_3\theta mL}{(\theta+I)^3} +\frac{2\lambda_3p_Ig_IL}{(g_I+I)^3} -\frac{2\lambda_3u_0\kappa L^2C}{(\kappa+I)^3} -\frac{2\epsilon_I\omega\zeta L}{\eta_3(\zeta+I)^3}\right]\\ &\quad \times\left(-\mu_I I + \phi C +\frac{\omega LI}{\zeta+I}\right)\,. \end{align*} \subsection*{Acknowledgements} This work was supported in part by generous funding from the W.M. Keck Foundation through the Harvey Mudd College Center for Quantitative Life Sciences as well as by the National Science Foundation Division of Mathematical Sciences, Grant Number DMS-0414011, ``Mathematical Modeling of the Chemotherapy, Immunotherapy, and Vaccine Therapy of Cancer." \begin{thebibliography}{20} \bibitem{BahKim1975} K. Bahrami and M. Kim, \newblock \emph{Optimal control of multiplicative control systems arising from cancer therapy} \newblock {IEEE Trans. Autom. Contr.}, (20), (1975),537--542. \bibitem{Burden2004} T. Burden, J. Ernstberger, and K. R. Fister, \newblock \emph{Optimal control applied to immunotherapy}, in \newblock {Discrete and Continuous Dynamical Systems-Series B}, (4) 1 , 2004, 135--146. \bibitem{dePGuFi} L. G. de Pillis, W. Gu, K. R. Fister, \emph{et al.}; \newblock \emph{Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls}, \newblock {Mathematical Biosciences}, (accepted work), 2006. \bibitem{dePillis2006} L. G. de Pillis, W. Gu, and A. E. Radunskaya; \newblock \emph{Mixed immunotherapy and chemotherapy of tumors: Modeling, applications and biological interpretations}, \newblock {Journal of Theoretical Biology}, (238) 4 (2006),841--862. \bibitem{DePRad01} L. G. de Pillis and A. E. Radunskaya; \newblock \emph{A mathematical tumor model with immune resistance and drug therapy: an optimal control approach}, \newblock {Journal of Theoretical Medicine}, (3) (2001), 79--100. \bibitem{DePRad02} L. G. de Pillis and A.E. Radunskaya; \newblock \emph{The dynamics of an optimally controlled tumor model: A case study}, \newblock {Mathematical and Computer Modelling}, (37)11 (2003),1221--1244. \bibitem{Wiseman2005} L. G. de Pillis, A. E. Radunskaya, and C.L. Wiseman; \newblock\emph{A validated mathematical model of cell-mediated immune response to tumor growth}, \newblock {Cancer Research}, (61)17 (2005), 7950--7958. \bibitem{Immune} K. R. Fister and J. H. Donnelly; \newblock\emph{Immunotherapy: An Optimal Control Theory Approach}, \newblock{SIAM Jounal on Applied Mathematics}, (2) (2005), 499-510. \bibitem{Fister2000} K. R. Fister and J. C. Panetta; \newblock \emph{Optimal control applied to cell-cycle-specific cancer chemotherapy}, in \newblock {SIAM J. Appl. Math}, (60) 3 (2000), 1059--1072. \bibitem{Fister2003} K. R. Fister and J. C. Panetta; \newblock \emph{Optimal control applied to competing chemotherapeutic cell-kill strategies}, in \newblock {SIAM J. Appl. Math}, (63) 6 (2003),1954--1971. \bibitem{Fleming1975} W. H. Fleming and R. W. Rishel; \newblock ``Deterministic and Stochastic Optimal Control'', \newblock Springer-Verlag, 1975. \bibitem{Gardner} S. N. Gardner; \newblock \emph{A mechanistic, predictive model of dose-response curves for cell cycle phase-specific and nonspecific drugs}, \newblock{ Cancer Research}, (60) (2000), 1417--1425. \bibitem{Gogas2007} H. J. Gogas, J. M. Kirkwood, V. K. Sondak; \newblock \emph{Chemotherapy for Metastatic Melanoma: Time for a Change}, \newblock \emph{Cancer}, (109) 3 (2007), 455--464. \bibitem{HSV1995} R. F. Hartl and S. P. Sethi and R. G. Vickson. \newblock \emph{A survey of the maximum principles for optimal control problems with state constraints}, \newblock {SIAM Review}, (37) 2 (1995), 181--218. \bibitem{Hermans2003} I. F. Hermans, T. W. Chong, M. J. Palmowski, A. L. Harris, V. Cerundolo; \newblock \emph{Synergistic Effect of Metronomic Dosing of Cyclophosphamide Combined with Specific Antitumor Immunotherapy in a Murine Melanoma Model}, \newblock {Cancer Research}, (63) (2003), 8408--8413. \bibitem{MISER3} L. S. Jennings, M. E. Fisher, K. L. Teo, C.J. Goh; \newblock \emph{{MISER3} {O}ptimal Control Software: {T}heory and User Maual}, \newblock Department of Mathematics, The University of Western Australia, Nedlands, WA 6907, Australia, 2004, \newblock Version 3. Available at {\tt http://www.cado.uwa.edu.au/miser/}. \bibitem{KS1991} M. I. Kamien, N. L. Schwartz; \newblock \emph{Dynamic Optimization: {T}he Calculus of Variations and Optimal Control in Economics and Management}, volume~31 in {``Advanced Textbooks in Economics''}, \newblock North-Holland, 2nd edition, 1991. \bibitem{Kim1977} M. Kim, S. Perry, K. B. Woo; \newblock \emph{Quantitative approach to the design of antitumor drug dosage schedule via cell cycle kinetics and systems theory}, \newblock Ann. Biomed. Engng., (5) (1977) 12--33. \bibitem{Kirk1970} D. E, Kirk \newblock ``Optimal Control Theory: An Introduction'', \newblock Dover Publications, Inc., 1970. \bibitem{Kirschner1997} D. Kirschner, S. Lenhart, S. Serbin; \newblock \emph{Optimal control of the chemotherapy of HIV}, \newblock {J. Math. Biol.}, (35) (1997), 775-792. \bibitem{Kirschner1998} D. Kirschner and J.C. Panetta; \newblock \emph{Modeling immunotherapy of the tumor - immune interaction}, \newblock { J. Math. Biol.}, (37) (1998), 235--252. \bibitem{Kofler2006} D. M. Kofler, C. Mayr, C. M. Wendtner; \newblock \emph{Current status of immunotherapy in B cell malignancies}, \newblock Current Drug Targets, (7) 10 (2006), 1371--1374. \bibitem{Krener1977} A. J. Krener, \newblock \emph{The high order maximal principle and its application to singular extremals}, \newblock {SIAM J. Control Optimization}, (15) 2 (1977) 256--293. \bibitem{Kuz} V.~Kuznetsov and G. D. Knott; \newblock \emph{Modeling tumor regrowth and immunotherapy} \newblock Math and Comp. Modelling, (33) (2001), 1275--1287. \bibitem{Liu2006} G. Liu, K. L. Black, J. S Yu; \newblock \emph{Sensitization of malignant glioma to chemotherapy through dendritic cell vaccination}, \newblock {Expert Review of Vaccines - Future Drugs}, (5) 2 (2006), 233--247. \bibitem{Neri2006} B. Neri, L. Vannozzi, C. Fulignati, P. Pantaleo, D. Pantalone, C. Paoletti, F. Perfetto, M. Turrini, R. Mazzanti; \newblock \emph{Long-term survival in metastatic melanoma patients treated with sequential biochemotherapy: report of a Phase II study} \newblock {Cancer Investigation}, (24) 5 (2006), 474--478. \bibitem{Urs2004} U. Ledzewicz, T. Brown, H. Schattler; \newblock \emph{Comparison of optimal controls for a model in cancer chemotherapy with $l_1$ and $l_2$ type objectives}, \newblock {Optimization Methods and Software}, (19) 3-4 (2004), 339--350. \bibitem{Urs2006} U. Ledzewicz and H. Schattler; \newblock \emph{Drug resistance in cancer chemotherapy as an optimal control problem}, \newblock {Discrete and Continuous Dynamical Systems - Series B}, (6) 1 (2006), 129--150. \bibitem{Lukes1982} D. L. Lukes, \newblock ``Differential equations: Classical to controlled'', volume 162, \newblock Academic Press, 1982. \bibitem{Mur} J. M. Murray, \newblock \emph{Some optimality control problems in cancer chemotherapy with a toxicity limit}, \newblock {Math. Biosci.}, (100) (1990), 49--67. \bibitem{Pontryagin1962} L. S. Pontryagin, V.G. Boltyanskii, R. V. Gamkrelidze, E. F. Mishchenko; \newblock ``The Mathematical Theory of Optimal Processes'', \newblock Gordon and Breach, 1962. \bibitem{Riker2007} A. I. Riker, S. Radfar, S. Liu, Y. Wang, H. T Khong; \newblock \emph{Immunotherapy of melanoma: a critical review of current concepts and future strategies}, \newblock {Expert Opinion on Biological Therapy}, (7) 3 (20070, 345-358. \bibitem{Rosti1998} G. Rosti, \newblock ``Dose and Timing: The Pillars of Successful Chemotherapy'' \newblock Elsevier Health Sciences, 1998. \bibitem{Rubinov2007} E. Gez, R. Rubinov, D. Galtini, S. Keretyk, L.A. Best, T. Mashiach, O. Native, A. Stein, A. Kuten; \newblock \emph{Immuno-chemotherapy in metastatic renal cell carcinoma: long-term results from the rambam and linn medical centers, Haifa, Israel}, \newblock {Journal of Chemotherapy}, (19) 1 (2007), 79--84. \bibitem{Sinkovic2006} J. G. Sinkovic and J. C. Horvath; \newblock \emph{Evidence accumulating in support of cancer vaccines combined with chemotherapy: a pragmatic review of past and present efforts}, \newblock {International Journal of Oncology}, (29) 4 (2006), 765--777. \bibitem{Swan} G. W. Swan, \newblock \emph{Optimal control applications in biomedical engineering-a survey}, \newblock { Opt. Control Appl. \& Meth.}, (2) (1981), 311--314. \bibitem{SwanVincent} G. W. Swan and T. L. Vincent; \newblock \emph{Optimal control analysis in the chemotherapy of IgG Multiple Myeloma}, \newblock {Bull. of Math Bio}, (39) (1977) , 317--337. \bibitem{Swierniak2003} A. Swierniak, U. Ledzewicz, H. Schattler; \newblock \emph{Optimal control for a class of compartmental models in cancer chemotherapy}, \newblock {Int. J. Appl. Math. Comput. Sci.}, (13) 3 (2003), 357--368. \end{thebibliography} \end{document}