\documentclass[a4paper]{article}
\usepackage[english]{babel}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{indentfirst}
\usepackage{wrapfig}
\usepackage{amsmath}
\usepackage{breakcites}
\usepackage[a4paper,top=2cm,bottom=2cm,left=2cm,right=2cm,marginparwidth=.5cm]{geometry}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage[colorinlistoftodos]{todonotes}
\usepackage[colorlinks=true, allcolors=blue]{hyperref}
\title{Optimal Therapeutic Window Model of IFN-$\alpha$ and Imatinib in Combination Therapy in CML Patients}
\author{Arianna Kalkandis, William McKenzie and Cynthia Sanchez}
\date{May 1, 2018}
\begin{document}
\maketitle
\begin{abstract}
Chronic Myeloid Leukemia is a type of cancer that starts in certain blood-forming cells of the bone marrow. Interferon-$\alpha$ was once the standard front-line treatment producing remission rates of only 28.3 percent in 1991 \cite{Tal}. After the highly effective drugs first became available in 2001, survival rates have increased immensely. According to the American Cancer Society, one large study of CML patients treated with a drug called imatinib found that about 90 percent of them were still alive 5 years after starting treatment. While imatinib has changed the way oncologists treat CML, remission is common after extended gaps in treatment.
In this paper, we will explore the long-term dynamics of CML under treatment through the use of use of theoretical and mathematical components. We closely base our methods upon the approach of Urszula Ledzewicz and Helen Moore \cite{Led}. We will introduce our unique model and explain component selections, while we move towards understanding the optimal interactions of imatinib and interferon-$\alpha$ against dormant and proliferating CML cells. Our future work involving optimal control dynamics will be briefly introduced and further solutions are currently ongoing.
\end{abstract}
\section{Introduction}
\linespread{1.5}
\large
Chronic myeloid leukemia (CML) is a type of cancer that begins in the bone marrow and later accumulates in the blood. This type of cancer can form slowly in its first chronic stage; as it passes through the second accelerated stage to the final blast phase, CML can become fatal-fast \cite{Led}. According to the Leukemia and Lymphoma Society, in 2017, 62,130 people were diagnosed with leukemia. Approximately every 9 minutes, someone in the US dies from a form of blood cancer.
CML is caused by a chromosomal translocation. Translocations occur when a gene on one chromosome is shifted onto another nonhomologous chromosome. This genetic abnormality can result in a wide array of lethal defects if the shift is not balanced. CML results from an unbalanced translocation where the ABL gene on chromosome 9 is shifted so that it sits next to the BCR gene on chromosome 22 forming what is called the "Philadelphia Chromosome" \cite{Prud}.
Recently, scientists have made vast strides in understanding the connection between changes in DNA and CML. Human genes play a crucial role in cell regulation. Specific genes called oncogenes are responsible for telling our cells when to grow and divide. Tumor suppressor genes signal to our cells to slow division or promote appropriate cell death \cite{Cong}. The newly formed BCR-ABL complex causes the gene to constantly be turned on, resulting in uncontrolled proliferation. Understanding the BCR-ABL gene plays a crucial role in CML detection and treatment.
Interferon is a protein released by immune cells when the body detects a threat. Interferons do not directly attack cancer cells; they slow down or stop cells from dividing by modifying the functions of the genes that regulate the secretion of several proteins that affect cell growth \cite{Will}. Interferons can also reduce the ability of the cancer cells to protect themselves from the immune system while simultaneously strengthening the immune system. Interferons can be replicated and used as a drug to fight cancer, such as interferon-alpha (IFN-$\alpha$).
IFN-$\alpha$ is administered to CML patients via injection once daily to stimulate the body’s own immune system. IFN-$\alpha$ was once the standard course of treatment for patients with CML; the majority of cases in the chronic stage could achieve stable remissions \cite{Tal}. Today, tyrosine kinase inhibitors (TKIs) have emerged as the forefront of successful targeted therapy. TKIs block the function of tyrosine kinase enzymes, preventing cell signaling and growth. This is significant for such enzymes are overly active resulting in uncontrolled cell proliferation and growth. The TKI we will be focusing on in our model is imatinib mesylate (known widely for its brand name Gleevec
\textsuperscript{\textregistered}).
Gleevec is generally administered to patients orally through a 400 mg tablet taken once daily. Imatinib is able to target and bind to the active site of the BCR-ABL kinase protein. This domain is a position reserved as an ATP-binding site. Now blocked by imatinib, the cell is unable to transfer a phosphate group to the tyrosine on the protein. This, prevents activation of the protein, blocks cell signaling to the nucleus and the leukemic cell is forced to induce apoptosis \cite{Toma}. The most innovative aspect of imatinib is its ability to target leukemic cells through detection of the specific BCR-ABL protein.
According to the National Cancer Institute, the 5-year survival rate of patients with CML in 1998 was only 39.6 percent. Highly effective drugs first became available in 2001 with the FDA approval of imatinib. In 2010, the survival rate was almost 70 percent. Today, someone with CML who is in remission after two years of imatinib treatment has the same life expectancy as someone who doesn’t have cancer.
It is important to note that while imatinib is an innovative form of targeted cancer therapy, the extremely harsh side effects and high toxicity of imatinib makes treatment difficult for patients to endure\cite{Tal}. Furthermore, the mechanisms in which IFN-$\alpha$ therapy works is not well understood \cite{Tal}. This highlights the importance in finding the optimal therapeutic window for using IFN-$\alpha$ and imatinib through combination therapy.
According to laboratory studies, quiescent CML stem cells are the most resistant to Imatinib-mesylate induced apoptosis \cite{Chu}. These quiescent cells are located in the bone marrow and remain dormant almost throughout life, only to awake in an instance of injury or blood loss. The possibility of specifically waking up these dormant stem cells opens up new prospects for cancer treatment\cite{elisa}, \cite{Anne}.
\section{A Mathematical Model for the Treatment of CML}
In our model, we will mimic the approach of Urszula Ledzewicz and Helen Moore. Our future goal is to analyze treatment as an optimal control problem. The CML cells will be divided into two separate populations: quiescent (Q) and proliferating (P). If Q cells are not needed to divide and to replenish tissue cells, they temporarily stop to progress through the cell cycle until further divisions are required. The two therapies we will be modeling are IFN-$\alpha$ and imatinib. The T-cell concentration (T) will demonstrate the effector T-cell's change in immune response as a result of both drugs. The use of theoretical and mathematical components is used to enlighten the reader on the long-term dynamics of CML under treatment \cite{Led}.
\subsection{A Brief Review of the Mathematical Model}
\centerline{\includegraphics[width=18cm,height=15.5cm]{FinalModel.JPG}}
\vskip 1.5 cm
\textbf{Figure 1.} "Diagram of the dynamical system. The green circular areas represent the 'populations' of the model. Solid long black arrows extending from or to the populations represents the changes in numbers, with inward-pointing arrows representing increases and outward-pointing arrows decreases. Dashed long black arrows indicate indirect effects on those increases and decreases. Bars represent inhibition of a production or an indirect effect, due to the represented treatment; arrows represent amplification of a rate.
The effects of IFN-$\alpha$ are shown in blue, and the effects of imatinib are shown in red" \cite{Led}.
\vskip .5 cm
\large
Our analysis will be dependent on fixing imatinib to a constant drug concentration, allowing us to analyze results with and without treatment. Representing the pharmacodynamic effects of the drugs using Michaelis-Menten terms results in the following equations:
\begin{align}
\begin{split}
\frac{dQ}{dt}&= r_Q Q\Big(1-\frac{U0_{max,2} u_0}{U0C_{50} + u_0}\Big) - \delta_Q Q \Bigg[\Big(1 + \frac{U0_{max,_1} u_0}{U0C_{50} + u_0}\Big)\Big(1 + \frac{E_{max,1} E}{EC_{50} + E})\Bigg]\\[15pt]
&\phantom{aa}-k_p\Big(\frac{K_{max} u_0}{K_n + u_0}\Big) Q
\end{split}\\[25pt]
\begin{split}
\frac{dP}{dt}&= \Bigg[\Big(1-\frac{U0_{max,2} u_0}{U0C_{50} + u_0}\Big)\Big(1-\frac{U1_{max,1} u_1}{U1C_{50} + u_1}\Big)\Big(r_{p}Pln\Big(\frac{P_{ss}}{P}\Big)\Bigg]+k_{p}Q\Big(1 + \frac{K_{max} u_0}{K_{n} + u_0}\Big)\Big(1-\frac{U1_{max,1} u_1}{U1C_{50} + u_1}\Big)\\
&\phantom{aa}-\delta_{P}P\Big(1 + \frac{U0_{max,_1} u_0}{U0C_{50} + u_0}\Big)\Big(1 + \frac{U1_{max,_2} u_1}{U1C_{50} + u_1}\Big)-\delta_{P}P\Big(1 + \frac{U0_{max,_1} u_0}{U0C_{50} + u_0}\Big)\Big(\frac{E_{max},_2 E}{EC_{50} + E}\Big)
\end{split}\\[25pt]
\frac{dE}{dt}&= s_E \Bigg[1 + \Big(1 + \frac{U0_{max,3} u_0}{U0C_{50} + u_0}\Big)
\Big(\frac{P_{max,1} P}{PC_{50} + P})\Bigg]E ln \Big(\frac{E_{ss}}{E}\Big) -
\delta_E E\Big(1 + \frac{P_{max,2} P}{PC_{50} + P}\Big)
\end{align}
\section{Model Components Explained}
{\normalsize \textit{In this section eq. (1), (2), and (3) are put into words, linking the biological and pharmacodynamic effects.}}
\vskip .2 cm
Figure 1, label (1) demonstrates inhibition against the natural reproduction of quiescent leukemic stem cells cause by IFN-$\alpha$. Through phosphorylation, cyclin-dependent kinases signal to the cell that it is ready to move on to the next stage of the cell cycle. By disrupting the signaling pathway, IFN-$\alpha$ inhibits the cell from growing properly, as well as reproducing at its natural rate \cite{Tal}.
\vskip .1 cm
Corresponding Term in $\frac{dQ}{dt}$:
$$r_Q Q\Big(1-\frac{U0_{max,2} u_0}{U0C_{50} + u_0}\Big)$$
\vskip .5 cm
Figure 1, label (2) represents the increased rate in which quiescent leukemic cells move to the proliferating leukemic cell population. Researchers have discovered that through combination therapy they are able to work around quiescent leukemic cell resistance to tyrosine kinase inhibitors. Researchers found that IFN-$\alpha$ increases the turnover rate in which quiescent cells become proliferating cells, maximizing the efficiency and capabilities of Imatinib \cite{Prud}. IFN-$\alpha$ allows imatinib to fight leukemic stem cells by breaking their dormancy and driving them to enter the proliferating population.
\vskip .1 cm
Corresponding Term in $\frac{dQ}{dt}$:
$$ -k_p\Big(\frac{K_{max} u_0}{k_n + u_0}\Big) Q$$
\\
Corresponding Term in $\frac{dP}{dt}$:
$$+k_{p}Q\Big(1 + \frac{K_{max} u_0}{K_{n} + u_0}\Big)\Big(1-\frac{U1_{max,1} u_1}{U1C_{50} + u_1}\Big) $$
\vskip .5 cm
Figure 1, label (3) displays the indirect resistance as a consequence of imatinib treatment. Chemokine receptor type 4 (CXCR4) is a molecule equipped with potent chemotactic activity for lymphocytes. It is expressed mostly on immature and mature hematopoietic cell types. The overexpression of BCR-ABL was reported to down-regulate CXCR4 expression, and this is associated
with the cell migration defects in CML. "It is proposed TKIs such as imatinib may restore CXCR4 expression and cause the migration of
CML cells to bone marrow microenvironment niches, which in turn results in acquisition of stroma-mediated chemoresistance of CML progenitor cells" \cite{Linh} A negative effect of imatinib is the unintended migration of CML cells. Therefore, we have included a dashed bar to demonstrate the inhibition that takes place.
Corresponding Term in $\frac{dP}{dt}$:
$$+k_{p}Q\Big(1 + \frac{K_{max} u_0}{K_{n} + u_0}\Big)\Big(1-\frac{U1_{max,1} u_1}{U1C_{50} + u_1}\Big) $$
\vskip .5 cm
Figure 1, label (4) demonstrates inhibition caused by IFN-$\alpha$. A cell must grow to a certain size before it can split and continue to proliferate; IFN-$\alpha$ inhibits cell growth. Therefore, it prevents population P from growing at it’s normal rate since it is able to interfere with the all phases of the cell cycle \cite{Tal}. Through such inhibition of growth, IFN-$\alpha$ is also able to induce apoptosis as a result of the cell cycle disruption.
\vskip .5 cm
Corresponding Term in $\frac{dP}{dt}$:
$$\Big(1-\frac{U0_{max,2} u_0}{U0C_{50} + u_0}\Big)\Big(1-\frac{U1_{max,1} u_1}{U1C_{50} + u_1}\Big)
\Big(r_{p}Pln\Big(\frac{P_{ss}}{P}\Big)\Big)$$
\vskip .5 cm
Figure 1, label (5) involves inhibition caused by imatinib. Transcription is when the DNA is translated into mRNA, and eventually becomes protein. BCR-ABL transcription levels will drop with treatment of Imatinib; this signifies that the Imatinib is successfully binding to the BCR-ABL protein and inhibiting all enzymatic activity. BCR-ABL is a key indicator of active CML cells \cite{Chu}. With this site blocked, the cell is unable to undergo cell signaling processes. Therefore, Imatinib is slowing the rate in which the proliferating population makes new CML cells by binding in the ATP site.
\vskip .5 cm
Corresponding Term in $\frac{dP}{dt}$:
$$ \Big(1-\frac{U0_{max,2} u_0}{U0C_{50} + u_0}\Big)\Big(1-\frac{U1_{max,1} u_1}{U1C_{50} + u_1}\Big)
\Big(r_{p}Pln\Big(\frac{P_{ss}}{P}\Big)\Big)$$
\vskip .5 cm
Figure 1, label (6) represents the amplification of the death rate caused by imatinib. Many tumor cells are dependent of the BCR-ABL gene. Through the specific targeted binding, imatinib is not only able to block the ATP site and prevent cell signaling, but it is also able to disrupt the nucleus causing the cell to be unable to preform its normal anti-apoptotic functions. This inevitably leads to increased cell death. \cite{vig}
Corresponding Term in $\frac{dP}{dt}$:
$$ -\delta_{P}P\Big(1 + \frac{U0_{max,_1} u_0}{U0C_{50} + u_0}\Big)\Big(1 + \frac{U1_{max,_2} u_1}{U1C_{50} + u_1}\Big) $$
\vskip .5 cm
Figure 1, label (7) shows how IFN-$\alpha$ is increasing the rate it which the proliferating cells are dying through forced apoptosis \cite{Bek}. IFN-$\alpha$ can initiate the apoptotic signal through the Janus kinase/signal transducers and activators of transcription (JAK/STAT) pathway \cite{Kot}. The JAK/STAT pathway is the primary signaling system for various cytokines and growth factors. JAK activation activates cell proliferation, differentiation, cell migration and apoptosis \cite{Raw}. This increases the rate of proliferating cell death.
Corresponding Term in $\frac{dP}{dt}$:
$$ -\delta_{P}P\Big(1 + \frac{U0_{max,_1} u_0}{U0C_{50} + u_0}\Big)\Big(1 + \frac{U1_{max,_2} u_1}{U1C_{50} + u_1}\Big)$$
\vskip .5 cm
Figure 1, label (8) represents the increased rate in which T-cells are activated as a result of the proliferating P population. Studies in mice have shown that IFN-$\alpha$ increases the activation rate of effector T-cells, thus speeding up the reproduction rate, expansion and long term survival of the cytotoxic T-cells \cite{Tou}.
\vskip .5 cm
Corresponding Term in $\frac{dE}{dt}$:
$$s_E \Bigg[1 + \Big(1 + \frac{U0_{max,3} u_0}{U0C_{50} + u_0}\Big)
\Big(\frac{P_{max,1} P}{PC_{50} + P}\Big)\Bigg]E ln \Big(\frac{E_{ss}}{E}\Big)$$
\vskip .5 cm
Figure 1, label (9) indicates the fighter T-cell antitumor immune response. Natural killer cells utilize their receptors to identify the major histocompatibility complex on the surface of tumor cells and induce apoptosis \cite{Danier}.
Corresponding Term in $\frac{dP}{dt}$:
$$ -\delta_{P}P\Big(1 + \frac{U0_{max,_1} u_0}{U0C_{50} + u_0}\Big)\Big(\frac{E_{max},_2 E}{EC_{50} + E}\Big)$$
\vskip .5 cm
Figure 1, label (10) represents the death of the fighter T-cells that are killed in the biological battle. NK cells are dysfunctional in chronic phase CML patients at diagnosis and NK cell numbers among lymphocytes are reduced, worsening with disease progression to advanced and blast crisis phase CML \cite{Hugh}. Therefore, it is important to include the loss of T-cells in our model.
\vskip .5 cm
Corresponding Term in $\frac{dE}{dt}$:
$$-\delta_E E\Big(1 + \frac{P_{max,2} P}{PC_{50} + P}\Big)$$
\vskip .5 cm
Figure 1, label (11) shows that IFN-$\alpha$ increases the rate in which T cells attack quiescent cells. IFN-$\alpha$ increases the proliferation and activation rate of T-cells, making the population larger and stronger against fighting CML cells, regardless of the cell's stage in the cell cycle \cite{Moor}.
\vskip .5 cm
Corresponding Term in $\frac{dQ}{dt}$:
$$- \delta_Q Q\Bigg[\Big(1 + \frac{U0_{max,_1} u_0}{U0C_{50} + u_0}\Big)\Big(1 + \frac{E_{max,1} E}{EC_{50} + E}\Big)\Bigg]$$
\vskip .5 cm
Figure 1, label (12) shows how IFN-$\alpha$ increases the death rate of quiescent cells. Treatment successfully induces apoptosis in quiescent CML progenitors resistant to elimination by imatinib alone. This eliminated of CML stem cells capable of engrafting was shown in mice \cite{zhang}. Therefore, the rate is amplified by IFN-$\alpha$.
\vskip .5 cm
Corresponding Term in $\frac{dQ}{dt}$:
$$- \delta_Q Q\Bigg[\Big(1 + \frac{U0_{max,_1} u_0}{U0C_{50} + u_0}\Big)\Big(1 + \frac{E_{max,1} E}{EC_{50} + E}\Big)\Bigg]$$
\vspace{7 cm}
\begin{center}
\begin{table}[ht]
\centering
\begin{tabular}{ |c|c| }
\hline
\textbf{Symbol} & \textbf{Interpretation} \\
\hline
Q & concentration of quiescent leukemic cells\\
\hline
P & concentration of proliferating leukemic cells \\
\hline
$P_{ss}$ & carrying capacity of proliferating leukemic cells\\
\hline
E & effector T cells\\
\hline
$E_{ss}$ & carrying capacity of effector T cells\\
\hline
$r_Q$ & replication rate constant of quiescent cells \\
\hline
$\delta_Q$ & natural death rate constant of quiescent cells \\
\hline
$k_P$ & rate constant for Q differentiating into P \\
\hline
$r_p$ & replication rate constant of proliferating leukemic cells \\
\hline
$\delta$ & natural death rate constant of proliferating leukemic cells \\
\hline
$s_E$ & growth rate constant for effector T cells \\
\hline
$\delta_E$ & natural death rate constant of effector T cells \\
\hline
$P_{max,1}$ & maximum stimulation effect of proliferating cells P on effector T cells E \\
\hline
$P_{max,2}$ & maximum death rate of T cells caused by P cell population\\
\hline
$PC_{50}$ & size of P with half the maximum effect \\
\hline
$E_{max,1}$ & maximum effect of effector T cells E on quiescent leukemic cells Q \\
\hline
$E_{max,2}$ & maximum effect of effector T cells E on proliferating leukemic cells P \\
\hline
$EC_{50}$ & size of E with half the maximum effect\\
\hline
\end{tabular}
\caption{States and parameters for the dynamical systems}
\label{Table1}
\end{table}
\end{center}
\begin{center}
\begin{table}
\centering
\begin{tabular}{ |c|c| }
\hline
\textbf{Symbol} & \textbf{Interpretation}\\
\hline
$u_0$ & normalized concentration of IFN-$\alpha$ \\
\hline
$u_{0max,1}$ & maximum possible effect of $u_0$ on death of Q and P cells \\
\hline
$u_{0max,2}$ & maximum possible effect of $u_0$ on inhibiting growth of Q and P cells \\
\hline
$u_{0max,3}$ & maximum possible effect of $u_0$ on stimulating proliferation of effector T cells \\
\hline
$K_{max}$ & maximum possible effect of $u_0$ on enhancing transfer from Q into P \\
\hline
$U0C_{50}$ & concentration of $u_0$ that gives half the maximum effect \\
\hline
$u_1$ & normalized concentration of imatinib \\
\hline
$u_{1max,1}$ & maximum possible effect of $u_1$ on slowing transfer or quiescent cells Q into P \\
\hline
$u_{1max,2}$ & maximum possible effect of $u_1$ on death of proliferating cells P\\
\hline
$U1C_{50}$ & concentration of $u_1$ that gives half the maximum effect\\
\hline
\end{tabular}
\caption{Controls and pharmacodynamic parameters}
\label{Table2}
\end{table}
\end{center}
\vskip 20 cm
\section{System Properties for Constant Concentrations- Dimensionless Model}
{\textit For a step-by-step visualization of the scaling and reduction of parameters, refer to the Appendix (8.1-8.5). After such changes, we are left with the following three equations:
\vskip .1 cm
\begin{align}
\begin{split}
\frac{d\tilde{Q}}{dt}&= r_Q \tilde{Q}\Big(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) - \delta_Q \tilde{Q}\Bigg[\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(1 + \frac{E_{max,1} \tilde{E}}{1 + \tilde{E}}\Big)\Bigg] -k_p\Big(\frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) \tilde{Q}\\[15pt]
\frac{d\tilde{P}}{dt}&= \hat{r_{p}}\Big(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
\tilde{P}ln\Big(\frac{\tilde{P_{ss}}}{\tilde{P}}\Big) +\hat{k_{p}}\tilde{Q}\Big(1 + \frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) -\hat{\delta_{P}}\tilde{P}\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) \\[10pt]
&\phantom{a0}-\hat{\delta_{P}}\tilde{P}\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(\frac{\hat{E_{max}},_2 \tilde{E}}{1 + \tilde{E}}\Big)\\[15pt]
\frac{d\tilde{E}}{dt}&= s_E \Bigg[1 + \Big(1 + \frac{U0_{max,3} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
\Big(\frac{P_{max,1} \tilde{P}}{1 + \tilde{P}}\Big)\Bigg]\tilde{E} ln \Big(\frac{\tilde{E}_{ss}}{\tilde{E}}\Big) -
\delta_E \tilde{E}\Big(1 + \frac{P_{max,2} \tilde{P}}{1 + \tilde{P}}\Big)
\end{split}
\label{TildeODESystem}
\end{align}
}
\section{Analysis}
In this section, we will show our analysis of various treatment windows and regimens. Parameter values used will be displayed in a corresponding table.
\vskip 5 cm
\begin{center}
\begin{table}[ht]
\begin{tabular}{ |c|c|c|c|c|c| }
\hline
Trial 1: & $u_0$ = 0.25 & Trial 2: & $u_0$ = 0. 35 & Trial 3: & $u_0$ = 0.45 \\
\hline
\hline
\textbf{Parameters} & \textbf{Values} & \textbf{Parameters} & \textbf{Values} & \textbf{Parameters} & \textbf{Values} \\
\hline
\hline
$Q_0$ & 40 & $r_Q$ & 0.4 & $\delta_Q$ & 0.25\\
\hline
$P_0$ & 50 & $r_P$ & 4 & $\delta_P$ & 0.75 \\
\hline
$E_0$ & 20 & $s_E$ & 0.5 & $\delta_E$ & 0.4\\
\hline
$u_1$ & 400 mg & & & & \\
\hline
\hline
\end{tabular}
\caption{States and parameters for the dynamical systems}
\label{Table6}
\end{table}
\end{center}
\subsection{Without Treatment}
\centerline{\includegraphics[width=18cm,height=15.5cm]{Figure_2.jpg}}
\vskip 1.5 mm
{\textbf{Figure 2.1}} In our model, we gather that without interferon treatment the quiescent and proliferating cell populations experience unlimited growth.
\vskip .5 cm
\subsection{With Treatment}
\subsubsection{Constant Treatment}
\centerline{\includegraphics[width=18cm,height=15.5cm]{Figure_8.jpg}}
\textbf{Figure 2.2} As expected with constant treatment of imatinib and interferon (400 mg) we are able to significantly reduce the quiescent and proliferating cell populations.
\subsubsection{Alternating Treatment 1 (Weekly)}
\centerline{\includegraphics[width=18cm,height=15.5cm]{Figure_1.jpg}}
\textbf{Figure 2.3} In three separate trials we used increasing interferon dosage regimens of 250 mg, 350 mg and 450 mg. Our model appropriately responded with a decreased population of proliferating leukemic cells.
\vskip 1 cm
\subsubsection{Alternating Treatment 2 (26 Weeks)}
\centerline{\includegraphics[width=18cm,height=15.5cm]{Figure_3.jpg}}
\textbf{Figure 2.4} In this trial we tested 400 mg of interferon for 26 weeks and then discontinued treatment. Similar to human trials done in
The interferon-alpha revival in CML \cite{Tal}, the cancer cell population is revived and then grows exponentially.
\vskip 2 cm
\subsubsection{Alternating Treatment 3}
\centerline{\includegraphics[width=18cm,height=15.5cm]{figure_5.jpg}}
\textbf{Figure 2.5} In this trial we alternated treatment by administering a maximum dosage of interferon (600 mg) for 1 week and discontinuing treatment for the following 2 weeks.
\vskip 2 cm
\subsubsection{Alternating Treatment 4}
\centerline{\includegraphics[width=18cm,height=15.5cm]{Figure_4.jpg}}
\textbf{Figure 2.6} In this trial we alternated treatment by administering a maximum dosage of interferon (600 mg) for 2 week and discontinuing treatment for the following 1 week.
\vskip 2 cm
\section{Discussion and Conclusions}
We considered the dynamical behavior of a mathematical model for CML that incorporated two types of therapies- imatinib and IFN-$\alpha$. We began our analysis of the long-term dynamical behavior of quiescent and proliferating leukemic cells and immune effects (represented by effector T cells.) Parameter values were estimated based on clinical applications, while some are still unknown at this current phase in our work. Based on our analysis, we were able to see overall expected trends in treatment. These trends are P populations decreasing after implementation of treatment. Naturally, we also saw P populations growing quickly after termination of treatment, showing the unbounded growth that can occur.
Although we have taken the first steps towards achieving optimal control solutions, there is still much future work ahead. Our future work efforts can be found in the Appendix of this paper where we introduce our Hamiltonian function.
{\bf Acknowledgments:}
We closely base our methods upon the approach of Urszula Ledzewicz and Helen Moore's paper \cite{Led}. Cynthia Sanchez provided resources and guidance throughout all components of this paper.
\section{Appendix}
Here you will find a step-by-step visualization of the scaling and reduction of parameters for all three populations introduced in section 4.
\subsection{Scaling of Parameters}
\large
\begin{align*}
&\tilde{Q} = \frac{Q}{Q_{ref}}, \qquad
\tilde{P} = \frac{P}{PC_{50}} \qquad\text{and}\qquad
\tilde{E} = \frac{E}{EC_{50}}\\[20pt]
&K_{n} = U0C_{50}, \qquad
\tilde{u_{0}} = \frac{u_{0}}{U0C_{50}} \qquad\text{and}\qquad
\tilde{u_{1}} = \frac{u_{1}}{U1C_{50}}
\end{align*}
\vskip .1 cm
Then we have,
\vskip .1 cm
\begingroup
\fontsize{16pt}{12pt}\selectfont
\centerline{$\frac{E}{EC_{50}+E} = \frac{\tilde{E}}{1 + \tilde{E}}$}
\endgroup
and analogously for the other terms.
\subsection{dQ/dt Analysis}
\begin{equation}
\frac{d\tilde{Q}}{dt}= \frac{1}{Q_{ref}}\frac{dQ}{dt}
\end{equation}
\begin{equation}
\frac{d\tilde{Q}}{dt}= \frac{1}{Q_{ref}}[r_Q Q(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}) - \delta_Q Q[(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1 + \frac{E_{max,1} \tilde{E}}{1 + \tilde{E}})] -k_p(\frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}) Q]
\end{equation}
\begin{equation}
\frac{d\tilde{Q}}{dt}= r_Q \tilde{Q}(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}) - \delta_Q \tilde{Q}[(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1 + \frac{E_{max,1} \tilde{E}}{1 + \tilde{E}})] -k_p(\frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}) \tilde{Q}
\end{equation}
\subsection{dP/dt Analysis}
\begin{equation}
\frac{d\tilde{P}}{dt}= \frac{1}{PC_{50}}\frac{dP}{dt}
\end{equation}
\vskip .1 cm
\begin{equation}
\begin{split}
\frac{d\tilde{P}}{dt}= \frac{1}{PC_{50}}[(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1-\frac{U1_{max,1} \tilde{u_{1}}}{1 + \tilde{u_{1}}})
(r_{p}Pln(\frac{P_{ss}}{P})) +k_{p}Q(1 + \frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1-\frac{U1_{max,1} \tilde{u_{1}}}{1 + \tilde{u_{1}}})\\
-\delta_{P}P(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1 + \frac{U1_{max,_2} \tilde{u_{1}}}{1 + \tilde{u_{1}}}) -\delta_{P}P(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(\frac{E_{max},_2 \tilde{E}}{1 + \tilde{E}})]
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\frac{d\tilde{P}}{dt}= (1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1-\frac{U1_{max,1} \tilde{u_{1}}}{1 + \tilde{u_{1}}})
(r_{p}\tilde{P}ln(\frac{P_{ss}}{\tilde{P}PC_{50}})) +k_{p}\frac{Q_{ref}}{PC_{50}}\tilde{Q}(1 + \frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1-\frac{U1_{max,1} \tilde{u_{1}}}{1 + \tilde{u_{1}}})\\
-\delta_{P}\tilde{P}(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1 + \frac{U1_{max,_2} \tilde{u_{1}}}{1 + \tilde{u_{1}}}) -\delta_{P}\tilde{P}(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(\frac{E_{max},_2 \tilde{E}}{1 + \tilde{E}})
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\frac{d\tilde{P}}{dt}= (1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1-\frac{U1_{max,1} \tilde{u_{1}}}{1 + \tilde{u_{1}}})
(r_{p}\tilde{P}ln(\frac{\tilde{P_{ss}}}{\tilde{P}})) +\tilde{k_p}\tilde{Q}(1 + \frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1-\frac{U1_{max,1} \tilde{u_{1}}}{1 + \tilde{u_{1}}})\\
-\delta_{P}\tilde{P}(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(1 + \frac{U1_{max,_2} \tilde{u_{1}}}{1 + \tilde{u_{1}}})
-\delta_{P}\tilde{P}(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}})(\frac{E_{max},_2 \tilde{E}}{1 + \tilde{E}})
\end{split}
\end{equation}
\subsection{dE/dt Analysis}
\begin{equation}
\frac{d\tilde{E}}{dt} = \frac{1}{EC_{50}}\frac{dE}{dt}
\end{equation}
\begin{equation}
\frac{d\tilde{E}}{dt}= \frac{1}{E_{ref}}\Bigg[s_E \Big[1 + (1 + \frac{U0_{max,3} \tilde{u_{0}}}{1 + \tilde{u_{0}}})
(\frac{P_{max,1} \tilde{P}}{1 + \tilde{P}})\Big]E \ln \Big(\frac{E_{ss}}{E}\Big) -
\delta_E E[(1 + \frac{P_{max,2} \tilde{P}}{1 + \tilde{P}})\Bigg]
\end{equation}
\begin{equation}
\frac{d\tilde{E}}{dt}= s_E [1 + (1 + \frac{U0_{max,3} \tilde{u_{0}}}{1 + \tilde{u_{0}}})
(\frac{P_{max,1} \tilde{P}}{1 + \tilde{P}})]\tilde{E} ln (\frac{\tilde{E}_{ss}}{\tilde{E}}) -
\delta_E \tilde{E}[(1 + \frac{P_{max,2} \tilde{P}}{1 + \tilde{P}})
\end{equation}
\vskip .3 cm
Thus, if we rescale
\vskip .3 cm
{$k_p$ as {$\tilde{k_p} = \frac{Q_{ref}}{PC{50}} k_p$}
and the steady-state values as
\begin{align*}
&&\tilde{P_{ss}} = \large\frac{P_{ss}}{PC_{50}} \qquad\text{and}\qquad
\tilde{E_{ss}} = \frac{E_{ss}}{EC_{50}}\\[20pt]
\end{align*}
"then formally the equations are the same as before with all "$C_50$" values in the Michaelis-Menten expressions normalized to 1. All other parameters remain unchanged and even their interpretation is the same as before. For the theoretical analysis and numerical computations this eliminated five parameters and introduces a favorable scaling to the variables. Naturally, the original parameters are still calculated for an interpretation of the results." \cite{Led}
\subsection{Reduction to Uncontrolled System and Basic Dynamical System Properties}
Keeping the "$C_{50}$" parameters in their original formulation in the controls, we define new drug-dependent parameters as
\begin{align*}
&\hat{r_{p}} = (1-\frac{U1_{max,1} \tilde{u_{1}}}{1 + \tilde{u_{1}}})r_{p}, \qquad
\hat{k_{p}} = (1-\frac{U1_{max,1} \tilde{u_{1}}}{1 + \tilde{u_{1}}}) \tilde{k_p}, \qquad
\end{align*}
\begin{align*}
\hat{\delta_{P}} = (1 + \frac{U1_{max,2} \tilde{u_{1}}}{1 + \tilde{u_{1}}}) \delta_{P}, \qquad
\hat{E}_{max,2}=\frac{E_{max,2}}{1+\frac{U_{1max,2}\tilde{u_1}}{1+\tilde{u_1}}}
\end{align*}
\vskip .1 cm
With these identifications, the dynamical system with constant controls is identical with the uncontrolled system and therefore, without loss of generality, the analysis can be done on the uncontrolled system. Returning to the original notation without the carets, we thus consider the following equations:
\begin{equation}
\frac{d\tilde{Q}}{dt}= r_Q \tilde{Q}\Big(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) - \delta_Q \tilde{Q}\Bigg[\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(1 + \frac{E_{max,1} \tilde{E}}{1 + \tilde{E}}\Big)\Bigg] -k_p\Big(\frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}) \tilde{Q}
\end{equation}
\begin{equation}
\begin{split}
\frac{d\tilde{P}}{dt}= \hat{r_{p}}\Big(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
\tilde{P}ln\Big(\frac{\tilde{P_{ss}}}{\tilde{P}}\Big) +\hat{k_{p}}\tilde{Q}\Big(1 + \frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) -\hat{\delta_{P}}\tilde{P}\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) \\
-\hat{\delta_{P}}\tilde{P}\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(\frac{E_{max},_2 \tilde{E}}{1 + \tilde{E}}\Big)
\end{split}
\end{equation}
\begin{equation}
\frac{d\tilde{E}}{dt}= s_E \Bigg[1 + \Big(1 + \frac{U0_{max,3} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
\Big(\frac{P_{max,1} \tilde{P}}{1 + \tilde{P}}\Big)\Bigg]\tilde{E} ln \Big(\frac{\tilde{E}_{ss}}{\tilde{E}}\Big) -
\delta_E \tilde{E}\Big(1 + \frac{P_{max,2} \tilde{P}}{1 + \tilde{P}}\Big)
\end{equation}
\subsection{Results 2}
\centerline{\includegraphics[width=18cm,height=15.5cm]{Figure_6.jpg}}
\begin{center}
\begin{table}[ht]
\begin{tabular}{ |c|c|c|c| }
\hline
\textbf{Parameter} & \textbf{Trial 1 Dosage} & \textbf{Trial 2 Dosage} & \textbf{Trial 3 Dosage} \\
\hline
\hline
$u_0$ & 400mg & $ 0 mg$ & 400 mg \\
\hline
$u_1$ & 0 mg & $ 200 mg $ & 200 mg \\
\hline
\hline
\end{tabular}
\caption{States and parameters for the dynamical systems}
\label{Table3}
\end{table}
\end{center}
\centerline{\includegraphics[width=18cm,height=15.5cm]{Figure_7.jpg}}
In this figure we used the same dosage in trial 3 but the regimen involves alternating interferon dosage weekly.
\vskip 2 cm
\subsection{Parameter Values Used In Graphs}
\begin{center}
\begin{table}
\begin{tabular}{ |c|c|c|c| }
\hline
\textbf{Parameter} & \textbf{Values used Results 1} & \textbf{Parameter} & \textbf{Values Results 1} \\
\hline
\hline
$Q_0$ & 40 & $P_{max,1}$ & 2 \\
\hline
$P_0$ & 50 & $P_{max,2}$ & 5 \\
\hline
$E_0$ & 20 & $E_{max,1}$ & 0 \\
\hline
$u_0$ & see graphachieving & $E_{max,2}$ & 1\\
\hline
$u_1$ & 0.4 & $K_{max}$ & l \\
\hline
$r_Q$ & 0.4 & $U0C_{50}$ & 1.64 \\
\hline
$r_P$ & 4 & $U1C_{50}$ & 1 \\
\hline
$s_E$ & 0.5 & $PC_{50}$ & 10e7 \\
\hline
$\delta_Q$ & 0.25 & $EC_{50}$ & 2000 \\
\hline
$\delta_P$ & 0.75 & $K_{n}$ & 1.64 \\
\hline
$\delta_E$ & 0.4 & $P_{ss}$ & 10 \\
\hline
$k_p$ & 0.5 & $E_{ss}$ & 1.75 \\
\hline
$U0_{max,1}$ & 1 & $U1_{max,1}$ & 0.8 \\
\hline
$U0_{max,2}$ & 2 & $U1_{max,2}$ & 10 \\
\hline
$U0_{max,3}$ & 1 & & \\
\hline
\hline
\hline
\end{tabular}
\caption{States and parameters for the dynamical systems}
\label{Table4}
\end{table}
\end{center}
\begin{center}
\begin{table}
\begin{tabular}{ |c|c|c|c| }
\hline
\textbf{Parameter} & \textbf{Values used Results 2} & \textbf{Parameter} & \textbf{Values Results 2} \\
\hline
\hline
$Q_0$ & 100 & $P_{max,1}$ & 2 \\
\hline
$P_0$ & 200 & $P_{max,2}$ & 5 \\
\hline
$E_0$ & 125 & $E_{max,1}$ & 1 \\
\hline
$u_0$ & see graph & $E_{max,2}$ & 1\\
\hline
$u_1$ & 0.2 & $K_{max}$ & l \\
\hline
$r_Q$ & 0.6 & $U0C_{50}$ & 1.64 \\
\hline
$r_P$ & 7 & $U1C_{50}$ & 1 \\
\hline
$s_E$ & 2 & $PC_{50}$ & 10e7 \\
\hline
$\delta_Q$ & 0.5 & $EC_{50}$ & 2000 \\
\hline
$\delta_P$ & 0.6 & $K_{n}$ & 1.64 \\
\hline
$\delta_E$ & 0.4 & $P_{ss}$ & 10 \\
\hline
$k_p$ & 0.7 & $E_{ss}$ & 1.75 \\
\hline
$U0_{max,1}$ & 1 & $U1_{max,1}$ & 0.8 \\
\hline
$U0_{max,2}$ & 2 & $U1_{max,2}$ & 10 \\
\hline
$U0_{max,3}$ & 1 & & \\
\hline
\hline
\hline
\end{tabular}
\caption{States and parameters for the dynamical systems}
\label{Table5}
\end{table}
\end{center}
\vskip 20cm
\subsection{Further Research: Hamiltonian}
\subsubsection{Objective Function}
Our goal is to minimize the objective function:
$$J[u] = rV(T) + \int_{0}^{T} V(t)dt + \int_{0}^T su_0 dt$$
where $$ V(T) = Q(T) + P(T)$$
therefore the Hamiltonian becomes:
\begin{equation}
H = P(t) + Q(t) + s_2\tilde{u_{0}}(t) + \lambda_{Q}\tilde{Q'}
+ \lambda_{P}\tilde{P'} + \lambda_{E}\tilde{E'}
\end{equation}
\\
\begin{equation}
H = [independent of u_0] + f(u_0)K(t)
\end{equation}
\begin{align}
\begin{split}
H = \lambda_Q \bigg[ r_Q \tilde{Q} - \delta_Q \tilde{Q}(1 +
\frac{E_{max,1}\tilde{E}}{1 + \tilde{E}})\bigg] \\
+ \lambda_P \bigg[\hat{r_{Q}}\tilde{P}ln\Big(\frac{\tilde{P_{ss}}}{\tilde{P}}\Big) +\hat{k_{p}}\tilde{Q} -\hat{\delta_{P}}\tilde{P} - \hat{\delta_{P}}\tilde{P}\Big(\frac{\hat{E_{max}},_2 \tilde{E}}{1 + \tilde{E}}\Big)\bigg] \\
+\lambda_E \bigg[s_E + s_E\Big(\frac{P_{max,1} \tilde{P}}{1 + \tilde{P}}\Big)\tilde{E} ln \Big(\frac{\tilde{E}_{ss}}{\tilde{E}}\Big) -
\delta_E \tilde{E}\Big(1 + \frac{P_{max,2} \tilde{P}}{1 + \tilde{P}}\Big) \bigg]
\\ + P(t) + Q(t) + s_2\tilde{u_0}(t) \\
+\lambda_{Q}\bigg[ (-1)r_Q \tilde{Q}\Big( \frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) - \delta_Q \tilde{Q}\Big(\frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(1 + \frac{E_{max,1} \tilde{E}}{1 + \tilde{E}}\Big)
-k_p\Big(\frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) \tilde{Q}
\bigg] \\
\lambda_P \bigg[ (-1)\hat{r_{p}}\Big(\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
\tilde{P}ln\Big(\frac{\tilde{P_{ss}}}{\tilde{P}}\Big) + +\hat{k_{p}}\tilde{Q}\Big(\frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) -\hat{\delta_{P}}\tilde{P}\Big(\frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big) \\
-\hat{\delta_{P}}\tilde{P}\Big(\frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(\frac{\hat{E_{max}},_2 \tilde{E}}{1 + \tilde{E}}\Big)
\bigg] \\
+\lambda_{E}\bigg[ s_E \Big(\frac{U0_{max,3} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
\Big(\frac{P_{max,1} \tilde{P}}{1 + \tilde{P}}\Big)\tilde{E} ln \Big(\frac{\tilde{E}_{ss}}{\tilde{E}}\Big)
\bigg]
\end{split}
\end{align}
\subsubsection{Adjoint Equations}
\begin{equation}
\dot{\lambda} = \begin{bmatrix}
\lambda_{Q}' \\ \lambda_{P}' \\ \lambda_{E}'
\end{bmatrix}
= \begin{bmatrix}
\frac{dH}{d\tilde{Q}} \\ \frac{dH}{d\tilde{P}} \\ \frac{dH}{d\tilde{E}}
\end{bmatrix}
\end{equation}
\\
\subsubsection{Optimality Condition}
\begin{align}
\begin{split}
-\frac{dH}{d\tilde{Q}} = (-1)\Bigg[
1 + \lambda_{Q}\Big[r_Q \big(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\big) - \delta_Q \big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\big)\big(1 + \frac{E_{max,1} \tilde{E}}{1 + \tilde{E}}\big) \\ -k_p\big(\frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\big)
\Big] + \lambda_{P}\Big[\hat{k_{p}}\big(1 + \frac{K_{max} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\big)
\Big]
\Bigg]
\end{split}
\end{align}
\begin{align}
\begin{split}
-\frac{dH}{d\tilde{P}} = (-1)\Bigg[1 + \lambda_{P}\bigg[
\hat{r_{p}}\Big(1-\frac{U0_{max,2} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
\big(ln\Big(\frac{\tilde{P_{ss}}}{\tilde{P}}\Big) - 1\Big)\\ -\hat{\delta_{P}}\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
-\hat{\delta_{P}}\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(\frac{\hat{E_{max}},_2 \tilde{E}}{1 + \tilde{E}}\Big)
\bigg] \\ + \lambda_{E}\Big[
s_E \Big(1 + \frac{U0_{max,3} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)
\big(\frac{P_{max,1} }{(1 + \tilde{P})^2}\big)\tilde{E} ln \Big(\frac{\tilde{E}_{ss}}{\tilde{E}}\Big) -
\delta_E \tilde{E} \frac{P_{max,2}}{(1 + \tilde{P})^2}
\Big]
\Bigg]
\end{split}
\end{align}
\\ \\
\begin{align}
\begin{split}
-\frac{dH}{d\tilde{E}} = (-1)\Bigg[
\lambda_{Q}\bigg[
(-1) \delta_Q \tilde{Q}\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(\frac{E_{max,1}}{(1 + \tilde{E})^2}\Big)
\bigg] \\
+ \lambda_{P} \bigg[
(-1)\hat{\delta_{P}}\tilde{P}\Big(1 + \frac{U0_{max,_1} \tilde{u_{0}}}{1 + \tilde{u_{0}}}\Big)\Big(\frac{\hat{E_{max}},_2}{(1 + \tilde{E})^2}\Big)
\bigg] \\
+\lambda_{E} \bigg[
s_E \big[1 + (1 + \frac{U0_{max,3} \tilde{u_{0}}}{1 + \tilde{u_{0}}})
(\frac{P_{max,1} \tilde{P}}{1 + \tilde{P}})\big]\big[ ln (\frac{\tilde{E}_{ss}}{\tilde{E}}) - 1 \big] -
\delta_E (1 + \frac{P_{max,2} \tilde{P}}{1 + \tilde{P}})
\bigg]
\Bigg]
\end{split}
\end{align} \\
\subsubsection{Stationary Condition}
\begin{align}
\begin{split}
\frac{dH}{d \tilde{u_0}} = s_2
+\lambda_{Q}\bigg[ (-1)r_Q \tilde{Q}\Big( \frac{U0_{max,2} }{(1 + \tilde{u_{0}})^2}\Big) - \delta_Q \tilde{Q}\Big(\frac{U0_{max,_1} }{(1 + \tilde{u_{0}})^2}\Big)\Big(1 + \frac{E_{max,1} \tilde{E}}{1 + \tilde{E}}\Big) \\
-k_p\Big(\frac{K_{max} }{(1 + \tilde{u_{0}})^2}\Big) \tilde{Q}
\bigg]
\lambda_P \bigg[ (-1)\hat{r_{p}}\Big(\frac{U0_{max,2} }{(1 + \tilde{u_{0}})^2}\Big)
\tilde{P}ln\Big(\frac{\tilde{P_{ss}}}{\tilde{P}}\Big) + +\hat{k_{p}}\tilde{Q}\Big(\frac{K_{max} }{(1 + \tilde{u_{0}})^2}\Big) \\-\hat{\delta_{P}}\tilde{P}\Big(\frac{U0_{max,_1} }{(1 + \tilde{u_{0}})^2}\Big)
-\hat{\delta_{P}}\tilde{P}\Big(\frac{U0_{max,_1} }{(1 + \tilde{u_{0}})^2}\Big)\Big(\frac{\hat{E_{max}},_2 \tilde{E}}{1 + \tilde{E}}\Big)
\bigg] \\
+\lambda_{E}\bigg[ s_E \Big(\frac{U0_{max,3} }{(1 + \tilde{u_{0}})^2}\Big)
\Big(\frac{P_{max,1} \tilde{P}}{1 + \tilde{P}}\Big)\tilde{E} ln \Big(\frac{\tilde{E}_{ss}}{\tilde{E}}\Big)
\bigg]
\end{split}
\end{align}
%% References with bibTeX database:
%\bibliographystyle{model1-num-names}
\bibliographystyle{plain}
\bibliography{sample.bib}
\end{document}
