# Immune Therapeutic Strategies Using Optimal Controls with $L^{1}$ and
$L^{2}$ Type Objectives111Work partially supported by the NSF-DMS Award
1122462 and a Fullbright Scholarship.

O. Bara<EMAIL_ADDRESS>S.M. Djouadi<EMAIL_ADDRESS>J.D Day
<EMAIL_ADDRESS>S. Lenhart<EMAIL_ADDRESS>Department of Electrical
Engineering and Computer Science, University of Tennessee, Knoxville, TN
37996, United States  Department of Mathematics, University of Tennessee,
Knoxville, TN 37996, United States

###### Abstract

Therapeutic strategies to correct an excessive immune response to pathogenic
infection is investigated as an optimal control problem. The control problem
is formulated around a four dimensional mathematical model describing the
inflammatory response to a pathogenic insult with two therapeutic control
inputs which have either a direct pro- or anti-inflammatory effect in the
given system. We use Pontryagin’s maximum principle and discuss necessary
optimality conditions. We consider both an $L^{1}$ type objective functional
as well as an $L^{2}$ type objective. For the former, the presence of singular
control will be addressed. For each case, numerical simulations using a
nonlinear programming optimization solver to acquire different drug treatment
strategies are presented and discussed. The results provide insight for
possible treatment strategies and the methods could be a relevant tool for
future practice to assist in better prediction of clinical outcomes and
subsequently better treatment for patients.

###### keywords:

Optimal control, Biological system , Inflammatory response, Sepsis

††journal: Journal of Mathematical Bioscience

##  1 Introduction

The immune response is a complex mechanism triggered in an organism as a
result of biological or physical stress, such as the presence of a foreign
body (pathogenic infection) or trauma. The natural behavior of an organism to
respond to these phenomena is through stress adaptation by trying to eliminate
the invading pathogen threat in case of an infection while also promoting
tissue repair due to the self-harming effects of inflammation. However, an
excessive and dysregulated inflammatory response may lead to systemic
inflammatory response syndrome (SIRS) and sepsis, resulting in tissue damage,
organ dysfunction or even death. According to [1], the overall mortality is
approximately $30\%$, rising to $40\%$ in the elderly and is $50\%$ or greater
in patients with severe SIRS and sepsis.

With the withdrawal of previously approved drugs to treat sepsis ([2]) and the
acknowledgment that potential therapies for sepsis such as anti-tumor necrosis
factor-$\alpha$ appear to require accurate timing ([3], [4]), identifying
effective therapies and patient-specific protocols for curbing sepsis is
increasingly important. Using computational tools and methodologies is one
means by which to explore this and offer insight. In particular, applying
automatic control for biomedical therapeutic intervention has recently gained
interest in such areas as glucose control for diabetes or in critically ill
patients and for anesthesia depth control, as seen for example, in [5], [6],
[7], [8], and the references therein. The use of optimal control theory for
biomedical applications can be seen in [9], [10], and [11] for example; and,
more recently, the use of model predictive control (MPC) in [12] and [13]
which used the same mathematical model as used herein. MPC was also used in
[14] and [15] for diabetes treatment strategies and [16] for an application to
a multi-compartment respiratory system.

In the cases where MPC is used, an optimal control problem is essentially
solved for each specified time interval over the duration of the time span, as
though implementing the control measures in real time (i.e. online) as the
response is evolving. In this work, we investigate the _offline_ strategy that
solves the control problem for the entire time period of interest: from the
start of treatment till the end of the observation period. Within this
strategy, we explore differences in the form of the objective function used as
well as determine the effectiveness of the identified treatment protocol.

As in our prior work ([17], [18], [12], [13]), we utilize a phenomenological
ordinary differential equations model of the system inflammatory response to
pathogenic infection developed in [19] (see also [20]). This model provides a
dynamical system with rich behavior that is ideal for testing various
theoretical control strategies. A summary of its important features are as
follows:

  * 1.

This model is based on the early non-specific protective mechanism, namely,
the _innate immune response_ , in contrast to the mechanisms of _adaptive
immunity_. It is a lumped variable/lumped parameter model, representing an
abstract view of the dynamics of an acute inflammatory response and implicitly
considering affects of typical intervention strategies such as antibiotic use
or the administration of resuscitating fluids;

  * 2.

An anti-inflammatory mediator is included as a dynamic variable of the model.
This mediator plays an important role in directly mitigating inflammation as a
way to prevent excessive tissue damage caused by severe inflammation; and

  * 3.

The model’s biological relevance has been confirmed through its good
qualitative reproduction of severe systemic inflammatory states as observed in
a clinical setting.

In this paper, we explore the use of optimal control theory to derive an
immune control strategy for the chosen model consisting of manipulating the
pro- and/or anti-inflammatory mediator in order to reestablish the healthy
equilibrium of the virtual patient. Note that the model does not target a
specific pathogen and the control doses represent concentrations of
inflammatory mediators with idealized effect on the system. Even so, the
careful balance required for effective treatment is evident. For instance, a
large dose of pro-inflammatory therapy, while aiming to eliminate the
pathogen, would also pose a risk in causing possibly irreparable damage to
organs. Similarly, using a large amount of anti-inflammatory therapy could
suppress the negative effects of the inflammatory-damage feedback loop and
consequently preserve organ health; however, this intervention may also
suppress the positive effects of inflammation in controlling pathogen which
may then be uninhibited to rapidly grow. The dual objectives to eliminate
pathogen but also minimize damage via the administration of pro- and anti-
inflammatory mediators, respectively, necessitates the use of a control
algorithm in order to find a suitable balance that produces the desired
outcome: a healthy patient.

In [18] and [17] the necessary conditions of Pontryagin’s maximum principle
[21] were used in order to formulate a two point boundary value problem
(TPBVP) which was solved numerically with the forward-backward sweep method
[22] (indirect method). In the current work, we choose to use a direct
approach to the numerical solutions that relies on approximating the original
optimal control problem by first discretizing the state and control history
and then solving it using nonlinear programming (NL). Further details about
the solver used will be given in section 4. Constraints are handled more
easily in the direct method than in the indirect method where one must derive
adjoint equations along with the transversality and optimality conditions
before the TPBVP problem can be solved. Furthermore, in the indirect method,
the presence of state inequality constraints requires _a priori_ knowledge of
the number and duration of constrained arcs which is unknown beforehand [23].

We compare the results of using a quadratic or $L^{2}$ objective function
versus a linear or $L^{1}$ one. The use of a quadratic cost to derive the
optimal control can be explained by two important facts. First the
corresponding Hamiltonian will be strictly convex in the control with a unique
minimum, and second, the mathematical problem is more straightforward to
solve. Alternatively, using an $L^{1}$ or linear objective avoids distortion
due to the square used in the $L^{2}$ form which puts a small penalty on lower
doses and therefore, potentially biases the real effect of the control in the
cost function [24]. Also, the use of an $L^{1}$ objective often results in
“bang-bang” controls, meaning the control is applied at either its maximum
level or its minimum level, which provide simple dosing protocols whose
amounts do not vary greatly over time. In this respect, such controls are more
analogous to how medical treatments are administered in practice and are
considered natural choices as candidates for optimality [25].

The simulation results presented herein show some similarities when compared
to [18] and [17], in that the optimal control doses follow a certain strategy
whether considering a quadratic or a linear objective. As will be noted later,
the optimal inputs for a virtual patient whose inflammatory mediators are
elevated above some prescribed threshold are characterized by a bolus of pro-
inflammatory therapy followed just after by anti-inflammatory therapy,
consistent with earlier findings in [12].

The paper is organized as follows. Sections 2 and 3 respectively explain the
mathematical model and the optimal control problem. A brief background
regarding the numerical solver (PSOPT) is provided in Section 4 together with
the numerical results for both the quadratic and linear cases. Suggestions for
future extensions to this work are discussed in section 5.

##  2 Immune response model

We consider the mathematical model describing the inflammatory response for a
host-pathogen interaction as proposed in [19] and consider the model as
representing the immune dynamics of a virtual patient. There are four system
states of the model:

  * 1.

$P$ : the bacterial pathogen population that initiates the inflammatory
response;

  * 2.

$N$ : combined collection of early pro-inflammatory mediators such as
activated phagocytes (e.g. neutrophils) and the pro-inflammatory cytokines
they produce;

  * 3.

$D$ : a marker of global tissue damage which incites further inflammation and
which is useful for determining response outcomes; and

  * 4.

$C_{a}$: combined collection of anti-inflammatory mediators, such as cortisol,
Interleukin-10 (IL-10) and Transforming Growth Factor-beta (TGF-$\beta$).

In [12] it was first proposed to include non-negative, time-varying inputs to
account for therapeutic controls in this model. Here, we denote these as
follows:

  * 1.

$u_{p}(t)$: a pro-inflammatory enhancer that provides direct input into the
$N$ equation; and

  * 2.

$u_{a}(t)$: an anti-inflammatory enhancer that provides direct input into the
$C_{a}$ equation.

The dynamic equations describing the interactions of the four system
variables, $P$, $N$, $D$, and $C_{a}$, are governed by various rate constants
or parameters whose numerical values can be found in [19]. As was done in
[12], we consider different sets of parameter values to represent different
virtual patients’ inflammatory responses. We refer to the values given in [19]
as the reference set of parameter values.

Figure 1: Diagram of natural and enhanced immune dynamics in response to
pathogenic stimulus. $P$: replicating pathogen, $N$: early pro-inflammatory
immune mediators, $D$: marker of tissue damage/dysfunction caused by
inflammatory response, $C_{a}$: inhibitory anti-inflammatory mediators,
$u_{a}$ and $u_{p}$: time-varying input controls for the anti- and pro-
inflammatory therapy, respectively. Solid lines with arrow heads and dashed
lines with nodes/circular heads represent upregulation and inhibition,
respectively.

The diagram in Figure 1 characterizes the different interactions that exist
between the states of the inflammatory model. A solid line with an arrow head
indicates an up-regulation, whereas a dashed line with circular head indicates
inhibition or downregulation of a process. For instance, early pro-
inflammatory mediators, $N$, respond to the presence of pathogen, $P$, by
initiating self-recruitment of additional inflammatory mediators and $N$ is
therefore upregulated by the interaction with $P$ to attempt to efficiently
eliminate the pathogen. The self up-regulation that exists for $P$ is due to
replication. Furthermore, $N$ suppresses $P$ by eliminating it at some rate
when they interact; however, the inflammation caused by $N$ results in tissue
damage, $D$, which can provide a positive feedback into the early inflammatory
mediators depending on their intensity. To balance this, anti-inflammatory
mediators (e.g., cortisol, IL-10, TGF-$\beta$) can mitigate the inflammation
and its harmful effect by suppressing the response by $N$ and the effects of
$D$ in various ways. The ordinary differential equations modeling these
interactions is given by Equations (1) - (4):

$\displaystyle\dot{P}(t)$
$\displaystyle=k_{pg}P(t)\left(1-\frac{P(t)}{p_{\infty}}\right)-\frac{k_{pm}s_{m}P(t)}{\mu_{m}+k_{mp}P(t)}-k_{pn}F(N(t))P(t)$
(1) $\displaystyle\dot{N}(t)$
$\displaystyle=\frac{s_{nr}R(P(t),N(t),D(t))}{\mu_{nr}+R(P(t),N(t),D(t))}-\mu_{n}N(t)+u_{p}(t)$
(2) $\displaystyle\dot{D}(t)$
$\displaystyle=k_{dn}\frac{F(N(t))^{6}}{x_{dn}^{6}+F(N(t))^{6}}-\mu_{d}D(t)$
(3) $\displaystyle\dot{C_{a}}(t)$
$\displaystyle=s_{c}+k_{cn}\frac{F(N(t)+k_{cnd}D(t))}{1+F(N(t)+k_{cnd}D(t))}-\mu_{c}C_{a}(t)+u_{a}(t),$
(4)

where

$\displaystyle R(P,N,D)$ $\displaystyle=F(k_{np}P(t)+k_{nn}N(t)+k_{nd}D(t))\ \
\textrm{ and}$ $\displaystyle F(x)$
$\displaystyle=\frac{x}{1+\left(\frac{C_{a}(t)}{c_{\infty}}\right)^{2}}.$

As explained in [19], $F(x)$ represent a Hill function of order $6$ that
models the impact of activated phagocytes and their by-products ($N$) on the
creation of damaged tissue. Tissue damage ($D$) increases in a switch-like
sigmoidal fashion as $N$ increases, so that it takes sufficiently high levels
of $N$ to incite a moderate increase in damage and that the increase in damage
saturates with sufficiently elevated and sustained $N$ levels. The coefficient
(exponent) 6 was chosen, therefore, to model this aspect which ensured that
the healthy equilibrium had a reasonable basin of attraction for the $N/D$
system.

The model was formulated to represent a most abstract form of the complex
processes involved in the systemic inflammatory response. Hence, the variables
$N$ and $C_{a}$ represent multiple mediators with similar inflammatory
characteristics, and D is an abstract representation of tissue damage. This
abstraction reduces the description to four essential variables which also
allows for tractable mathematical analysis. Therefore the units of these
variables are in arbitrary units of $N$-units, $C_{a}$-units, $D$-units, since
they represent various types of cells and thus, they qualitatively rather than
quantitatively describe the response of the inflammatory mediators and by-
products each represents. Pathogen, $P$, units are more closely related to
numbers of pathogens or colony forming units (CFU), but abstract units
$P$-units are simply used as well and this population is scaled by
$10^{6}$/cc. More details about the model development can be found in [19].

The $C_{a}$ variable maintains a small positive background level when the
system is in steady state in the presence of no pathogen insult. Following the
setup used in [12], the reference parameter value for $C_{a}(0)$ is set to
$0.125$ and virtual patients have a value that is $\pm 25\%$ of the reference
value. In addition, five other parameters as well as the initial conditions
for $P$ and $C_{a}$ are made to have differing (positive) values from the
reference set to distinguish one virtual patient from another and from the
virtual patient having the reference set of parameters.

At the end of 168 hours222168 hours, which is equivalent to one week, is
considered to be a clinically relevant time period over which to observe a
patient experiencing a severe inflammatory response and when outcome would
most likely could be determined. Our results match this time frame, as for
example, from Figure 2 one can see that by 168h it is relatively clear to
determine to which equilibrium the time courses are converging, whether septic
death or the aseptic death region. , a healthy outcome of a virtual patient’s
inflammatory response given a particular initial pathogen is defined as the
ending state in which $D=0$, $P=0$, $N=0$, and $C_{a}$ is at its initial
background level. On the other hand, a patient is considered to have an
unhealthy outcome if the ending state has either $P=0$ with all other variable
values elevated above their background non-infection steady state values
(aseptic death outcome) or when $P$ is also elevated along with the other
variables (septic death outcome). The fact that a virtual patient may, or may
not, return to a healthy state depends on the parameter values and initial
condition. When inflammatory levels of a response are deemed excessive (i.e.
$N(t)>0.5$ in this work), then intervention becomes necessary to stabilize the
patient to his healthy equilibrium (i.e. to homeostasis). In the work here, we
consider two virtual patients whose inflammatory response outcomes are septic
death and aseptic death, respectively, in the absence of intervention
measures. Below, we summarize for each virtual patient the corresponding
parameter values used as well as the initial conditions which represents the
state of the patient’s inflammatory response some past time when the initial
infection occurred, and the inflammatory response has become elevated enough
to warrant therapeutic intervention.

Figure 2: Time courses of system states for natural (open loop) response for
virtual patient 1 (solid curves), resulting in a (_septic death outcome_) and
for virtual patient 2 (dashed curves), resulting in an (_aseptic death
outcome._)

  * 1.

Baseline Patient

    * (a)

Healthy equilibrium: $P(0)=0$, $N(0)=0$, $D(0)=0$, $C_{a}(0)=0.125$.

    * (b)

Parameter values of the reference patient: $k_{pg}=0.6$, $k_{cn}=0.04$,
$k_{nd}=0.02$, $k_{np}=0.1$, $k_{cnd}=48$, $k_{nn}=0.01$.

  * 2.

Patient 1 (septic death outcome)

    * (a)

Initial conditions: $P(0)=0.5360$, $N(0)=0.0660$, $D(0)=0.0477$,
$C_{a}(0)=0.1635$.

    * (b)

Parameter values differing from reference set: $k_{pg}=0.5846$,
$k_{cn}=0.0409$, $k_{nd}=0.0242$, $k_{np}=0.1211$, $k_{cnd}=49.1243$,
$k_{nn}=0.012$.

  * 3.

Patient 2 (aseptic death outcome)

    * (a)

Initial conditions: $P(0)=1.0017$, $N(0)=0.0711$, $D(0)=0.0732$,
$C_{a}(0)=0.1314$.

    * (b)

Parameter values differing from reference set: $k_{pg}=0.4746$,
$k_{cn}=0.0386$, $k_{nd}=0.0223$, $k_{np}=0.1116$, $k_{cnd}=46.3367$,
$k_{nn}=0.0112$.

The differential equations model using the parameter values and initial
conditions corresponding to each virtual patient is numerically solved and the
resulting time courses for each are shown in Figure 2. Solid curves display
the solution time courses for patient 1 and dashed curves are for patient 2.
As illustrated, in the absence of any control measures, each virtual patient’s
trajectory evolves toward either of the final unhealthy steady states. Next,
we demonstrate how interventions via positive inputs into the system can
redirect trajectories away from the unhealthy steady states and toward the
healthy steady state.

##  3 The Optimal Control Problems

We consider optimal control problems in which the controls enter linearly into
the system. After algebraically rearranging the expressions in equations
(1)-(4) and renaming the resulting parameter expressions from $a$ to $w$ as
shown in Table 1, the model system with controls can now be written as:

$\displaystyle\begin{pmatrix}\dot{P}\\\ \dot{N}\\\ \dot{D}\\\
\dot{C_{a}}\end{pmatrix}=f(x)+g(x)u=\begin{pmatrix}aP-
bP^{2}-\frac{cP}{d+eP}-\frac{z}{g+C_{a}^{2}}PN\vspace{2mm}\\\
\frac{hP+iN+jD}{kP+lN+mD+nC_{a}^{2}+o}-pN\vspace{2mm}\\\
\frac{qN^{6}}{(r+sC_{a}^{2})^{6}+tN^{6}}-\mu_{d}D\vspace{2mm}\\\
s_{c}+\frac{uN+vD}{g+C_{a}^{2}+gN+wD}-\mu_{c}C_{a}\end{pmatrix}+\begin{pmatrix}0&0\\\
1&0\\\ 0&0\\\ 0&1\end{pmatrix}\begin{pmatrix}u_{p}\\\ u_{a}\end{pmatrix}$ (5)

$a=k_{pg}$ | $e=k_{mp}$ | $i=c_{\infty}^{2}s_{nr}k_{nn}$ | $m=c^{2}_{\infty}k_{nd}$ | $q=c^{12}_{\infty}k_{dn}$
---|---|---|---|---
$b=\frac{k_{pg}}{P_{\infty}}$ | $z=k_{pn}c_{\infty}^{2}$ | $j=c^{2}_{\infty}s_{nr}k_{nd}$ | $s=k_{dn}$ | $n=k_{nr}$
$c=k_{pm}s_{m}$ | $g=c_{\infty}^{2}$ | $k=c^{2}_{\infty}k_{np}$ | $o=c^{2}_{\infty}k_{nr}$ | $t=c^{12}_{\infty}$
$d=k_{m}$ | $h=c_{\infty}^{2}s_{nr}k_{np}$ | $l=c^{2}_{\infty}k_{nn}$ | $p=k_{n}$ | $u=c^{2}_{\infty}k_{cn}$
$v=c^{2}_{\infty}k_{cn}k_{cnd}$ | $r=c^{2}_{\infty}k_{dn}$ | $w=c^{2}_{\infty}k_{cnd}$ |  |
Table 1: Parameters of the algebraically rearranged immune model with control
terms as given in (5).

where $f(x)$ and $g(x)$ represent a vector and a matrix, respectively, that
can be easily inferred from equation (5). The state vector is given by
$x=(P,N,D,C_{a})^{T}$ and the control vector by $u=(u_{p},u_{a})^{T}$. The
above alternate representation of the system, after algebraically rearranging
the immune model equations, provides a more useful form for the
differentiation that will be carried out below in solving the control problem
and allows a better understanding of the role of each state in the system. Two
case studies will be investigated: an optimal control problem with
$L^{1}$-type objective (linear cost term) versus one with an $L^{2}$-type
objective (quadratic cost term) in $J$.
For the quadratic case the objective functional to be minimized is therefore
given as

$\displaystyle
J(u_{p},u_{a})=\int_{0}^{t_{f}}a_{1}P(t)^{2}+a_{2}N(t)^{2}+a_{3}D(t)^{2}+a_{4}C_{a}(t)^{2}+B_{1}u_{p}(t)^{2}+B_{2}u_{a}(t)^{2}dt$
(6)

For the linear case, the objective functional is given according to

$\displaystyle
J(u_{p},u_{a})=\int_{0}^{t_{f}}a_{1}P(t)+a_{2}N(t)+a_{3}D(t)+a_{4}C_{a}(t)+B_{1}u_{p}(t)+B_{2}u_{a}(t)dt$
(7)

The goal here, for both functionals is to determine the optimal control pair
$u_{p}^{*}$ and $u_{a}^{*}$ such that

$\displaystyle
J(u_{p}^{*},u_{a}^{*})=\underset{(u_{p},u_{a})\in\mathbb{U}}{\min}J(u_{p},u_{a})$
(8)

where the basic control set is

$\mathbb{U}=\\{(u_{p},u_{a})\in(L^{\infty}(0,t_{f}))^{2}|\quad 0\leq
u_{p}(t)\leq M_{p},0\leq u_{a}(t)\leq M_{a},\quad a.e.\\}$

and other state constraints may be considered later.

Note that during our simulation we also added the final cost term to our
objective functional as
$\phi(x(t_{f}))=P(t_{f})+N(t_{f})+D(t_{f})+C_{a}(t_{f})$ for the linear case
and similarly
$\phi(x(t_{f}))=x(t_{f})^{T}x(t_{f})=P^{2}(t_{f})+N^{2}(t_{f})+D^{2}(t_{f})+C_{a}^{2}(t_{f})$
for the quadratic one.

The optimal state $x^{*}$, although satisfying the equality constraints
$\dot{x}=f(x(t),u(t))$, is further constrained by the initial and final
boundary condition, or point constraints. In our case, the initial state is
$x(0)=(P(0),N(0),D(0),C_{a}(0))=(0.536,\,0.066,\,0.0477,\,0.1635)$ for patient
1 and $(1.0017,0.0711,0.0732,0.1314)$ for patient 2 and the final state for
each is
$x(t_{f})=(P(t_{f}),N(t_{f}),D(t_{f}),C_{a}(t_{f}))=(0,\,0,\,0,\,0.125)$.

The approach known as the indirect method, used in our previous publications
[18, 17], attempts to preserve the infinite dimensional nature of the problem
[26] and uses Pontraygin’s maximum principle to solve a TPBVP. The
disadvantage of the indirect method lies in the fact that the boundary problem
with the transversality condition is often difficult to solve, and can become
further encumbered in the presence of path constraints. In this work, we take
an approach that finds a finite dimensional representation necessary to solve
the problem. Such approaches are known as direct methods and to achieve this,
we discretize the state and/or control variables and formulate a nonlinear
programming problem (NLP) to solve the optimization problem.

###  3.1 Existence of an optimal control

Practically, before attempting to numerically find an optimal control
solution, it is interesting to see if an optimal control indeed exists. The
existence of an optimal control for the $L^{1}$ and $L^{2}$ objective
functional can be obtained from standard results like in [27, 28] or by using
a weak convergence direct approach [29]. These results use the bounded range
of the controls and the states as well as the structure of the ODE system
being linear in the controls. Also the right hand sides of the ODEs are
continuous functions of the states.

######  Assumption 1.

If the generated inflammatory doses $u_{p}$ and $u_{a}$ and the initial states
are non-negative (i.e. $u_{p}(t)\geq 0,u_{a}(t)\geq 0\quad\forall t$ and
$P(0)\geq 0,\quad N(0)\geq 0,\quad D(0)\geq 0,\quad C_{a}(0)\geq 0$) then the
solution trajectories of the differential equations (1)-(4) remain in the
positive octant.

This is a reasonable assumption since the states represent biological entities
that are inherently positive.

##  4 Numerical solution of optimal control problem

Numerical solutions of our optimal control problem are found by applying
direct optimization via the use of PSOPT [30], an open source pseudospectral
optimization software package. PSOPT is written in C++ and uses direct
collocation methods, including Legendre or Chebyshev pseudospectral
discretization as well as local transcriptions such as the trapezoidal or
Hermite-Simpson integration methods. Sparse nonlinear programming is then used
to find local optimal solutions. PSOPT can be interfaced to, IPOPT, the open
source Interior Point Optimization solver ([31]) for solving large scale
optimization problems. This solver will be used for all simulations in this
work.

A good example that highlights the relationship between nonlinear programming
(NLP) and optimal control can be found in [32], where the Karush-Kuhn-Tucker
(KKT) necessary conditions approach the necessary condition of the maximum
principle as the number of nodes increases (i.e. as the step size
$h\rightarrow 0$. The NLP Lagrange multipliers can be interpreted as a
discrete approximation of the adjoint variables found using the maximum
principle, which motivates the following remark.

######  Remark 4.1.

Using the Lagrange multipliers provided by PSOPT, the numerical solutions can
be verified by showing that they satisfy the necessary conditions of
optimality with good accuracy.

We observed that the Hermite-Simpson integration method provided better
convergence results and accuracy $(\mathcal{O}(h^{4}))$ ([30]), when compared
to either trapezoidal methods provided in PSOPT which are known to have an
accuracy of $\mathcal{O}(h^{2})$. The Pseudospectral discretization method of
Legendre or Chebyshev proved to be slow to converge and not accurate although
these have been successfully applied elsewhere (e.g., [33], [34]). In
addition, since we found that the solution of the optimal control problem of
the previously mentioned two patients to be very similar, we will focus on
providing simulation results related to the patient with a septic outcome.

We note that the control objectives in this work are to minimize the levels of
pathogen and of damage such that $P=0$ and $D=0$, respectively, while at the
same time making sure to use the lowest possible dosing amounts. In addition,
we make the following assumption

######  Assumption 2.

If the generated inflammatory doses $u_{p}$ and $u_{a}$ and the initial states
are positive semidefinite (i.e. $u_{p}(t)\geq 0,u_{a}(t)\geq 0\quad\forall t$
and $P(0)\geq 0,\quad N(0)\geq 0,\quad D(0)\geq 0,\quad C_{a}(0)\geq 0$) then
the solution trajectories of the differential equations (1)-(4) remain in the
positive octant.

A number of grid points $N=1000$ is used together with a tolerance of
$10^{-6}$ for the NLP solver specifications. All simulations are performed for
a week ($168\,hours$) and, therefore, the sampling period, given the number of
grid points, is every $0.168$ hours. For each simulation we report the cost
function $J_{z}$, where $z$ specifies the corresponding subsection in which
the particular cost function is defined; the entry and exit times, $t_{1}$ and
$t_{2}$, respectively, of the singular arc, if one exists in the case where an
$L^{1}$ type objective function is used; and the numerical value of the
states, $P$, $N$, $D$, and $C_{a}$ at the final time, $t_{f}$. Note that in
all the control simulations it can be verified that by further integrating the
system from the final state $x(t_{f})$ the virtual patient will eventually
converge to its healthy equilibrium. In the next section we will formulate the
optimal control problem with an $L^{1}$ type objective and define the
necessary optimality condition of the maximum principle related to problem
given in (7).

###  4.1 Optimal Control with $L^{1}$ type Objective

We now formulate the optimal control problem with an $L^{1}$ type objective
and define the necessary optimality conditions of Pontraygin’s maximum
principle. The $L^{1}$-type objective function has the following form:

$\displaystyle
J_{L^{1}}={x(t_{f})}+\int_{0}^{t_{f}}\Big{[}{Ax(t)}+{Bu(t)}\Big{]}\,dt,$ (9)

where $x=[x_{1},x_{2},x_{3},x_{4}]^{T}$ and where
$A=[a_{1},a_{2},a_{3},a_{4}]$ together with $B=[b_{1},b_{2}]$ are row vectors
with positive constants. Also, $x(t)=[P(t),\,N(t),\,D(t),\,C_{a}]^{T}\geq 0$
and $u(t)=[u_{p},u_{a}]^{T}\geq 0$, where $\geq$ means componentwise.

To ensure that the lowest possible dosing amounts are used to achieve the
objective, we include the input doses, $u_{p}$ and $u_{a}$, in the cost
function $J_{L^{1}}$ above. The mediators $N$ and $C_{a}$ are also included in
the cost function since the steady state for these values in healthy patients
is assumed to be very low for $C_{a}$ and zero for $N$. This inclusion will
ensure to keep the level of inflammatory mediators at lower values. It was
also noticed that the optimal control solver converges faster in the presence
of $N$ and $C_{a}$ in the cost function.333The necessary inclusion here of $N$
and $C_{a}$ in the objective function differs from that used in our previous
studies using MPC ([12], [13]).

######  Remark 4.2.

Note that the existence of optimal controls with state constraints are more
difficult but can be obtained when variations can be done within the
admissible set of controls. For background on problems with state constraints
and the corresponding necessary conditions, see [35, 36]

The consequence of using an objective that is linear in the controls is that
the resulting time-varying controls will be “bang-bang”, meaning that
treatment protocols alternate between the maximum dose of pro- and/or anti-
inflammatory therapies and the minimum (i.e. zero) dose. Singular controls may
also arise and will be discussed later. Since we are minimizing the cost of
$J_{L^{1}}$, the standard Hamiltonian is given by:

$\displaystyle
H(x,u,\lambda)=Ax(t)+Bu(t)+\lambda(t)^{T}(f(x)+g(x)u(t)),\qquad\lambda\in\mathbb{R}^{4}.$
(10)

In what follows, we characterize the optimal control which is then illustrated
in several simulations where we explore various constraints on the dosing
controls as well as mixed state-control constraints.

######  Theorem 1.

(Characterization of the optimal control) Given two optimal controls
$u_{p}^{*}$, $u_{a}^{*}$ with the corresponding solutions $P^{*}$, $N^{*}$,
$D^{*}$, $C_{a}^{*}$ of the ordinary differential equation given in
Eq.(1)-Eq.(4), there exist adjoint equations associated with the system states
given by

$\displaystyle\dot{\lambda}_{1}$
$\displaystyle=-a_{1}-\lambda_{1}\Big{[}a-2bP-\frac{c}{d+eP}+\frac{ceP}{(d+eP)^{2}}$
$\displaystyle-\frac{zN}{g+C_{a}^{2}}\Big{]}-\lambda_{2}\Big{[}\frac{h(nC_{a}^{2}+o)}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}\Big{]};$
(11) $\displaystyle\dot{\lambda}_{2}$
$\displaystyle=-a_{2}-\lambda_{1}\left(\frac{zP}{g+C_{a}^{2}}\right)-\lambda_{2}\Big{[}\frac{i(nC_{a}^{2}+o)}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}$
$\displaystyle-p\Big{]}-\lambda_{3}\Big{[}\frac{6qtN^{11}}{((sC_{a}^{2}+r)^{6}+tN^{6})^{2}}\Big{]}-\lambda_{4}\Big{[}\frac{u(g+C_{a}^{2})}{(g+C_{a}^{2}+wD+gN)^{2}}\Big{]};$
$\displaystyle\dot{\lambda}_{3}$
$\displaystyle=-a_{3}-\lambda_{2}\Big{[}\frac{j(nC_{a}^{2}+o)}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}\Big{]}$
$\displaystyle+\lambda_{3}k_{d}-\lambda_{4}\Big{[}\frac{v(g+C_{a}^{4})}{(g+wD+N)^{2}}\Big{]};$
(12)
$\displaystyle\dot{\lambda}_{4}=-a_{4}-2\lambda_{1}\frac{zPNC_{a}}{(g+C_{a}^{2})^{2}}+\lambda_{2}\Big{[}\frac{2nC_{a}(hP+iN+jD)}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}\Big{]}$
$\displaystyle+\lambda_{3}\Big{[}\frac{12qsN^{6}(sC_{a}^{2}+r)^{5}}{((sC_{a}^{2}+r)^{6}+tN^{6})^{2}}\Big{]}-\lambda_{4}\Big{[}\frac{2(uN+vD)C_{a}}{(gN+wD+C_{a}^{2}+g)^{2}}-k_{c}\Big{]};$
(13)

and satisfying the following transversality condition

$\displaystyle\mathbf{\lambda}(t_{f})=\Big{\\{}\frac{\partial\phi[{x}(t_{f})]}{\partial{x}}\Big{\\}}=1.$
(14)

Furthermore,

$u_{p}^{*}=\left\\{\begin{array}[]{lll}0&if&b_{1}+\lambda_{2}>0\\\
N_{max}&if&b_{1}+\lambda_{2}<0\end{array}\right.$ (15)

$u_{a}^{*}=\left\\{\begin{array}[]{lll}0&if&b_{2}+\lambda_{4}>0\\\
C_{a_{max}}&if&b_{2}+\lambda_{4}<0\end{array}\right.$ (16)

###### Proof.

The adjoint equations and the transversality conditions follow from
Pontryagin’s maximum principle [21]. Differentiating the negative of the
Hamiltonian with respect to the states, $P$, $N$, $D$, and $C_{a}$, gives the
adjoint system in Eq.(11)-(13). According to the maximum principle, the
optimality equations are then given by

$\displaystyle\frac{\partial H}{\partial u_{p}}$
$\displaystyle=b_{1}+\lambda_{2}=0$ (17) $\displaystyle\frac{\partial
H}{\partial u_{a}}$ $\displaystyle=b_{2}+\lambda_{4}=0.$ (18)

These are known as switching functions and if we denote them by

$\displaystyle\phi_{1}(t)$ $\displaystyle=b_{1}+\lambda_{2}\textrm{ and}$ (19)
$\displaystyle\phi_{2}(t)$ $\displaystyle=b_{2}+\lambda_{4},$ (20)

then the optimal control that minimizes the Hamiltonian becomes

$u_{p}^{*}=\left\\{\begin{array}[]{lll}0&if&\phi_{1}(t)>0\\\
N_{max}&if&\phi_{1}(t)<0\\\ \end{array}\right.$ (21)
$u_{a}^{*}=\left\\{\begin{array}[]{lll}0&if&\phi_{2}(t)>0\\\
C_{a_{max}}&if&\phi_{2}(t)<0.\\\ \end{array}\right.$ (22)

∎

As was done in our previous studies, we define $N_{max}=0.5$ and
$C_{a_{max}}=0.62$ as the maximum allowed levels of the variables $N$ and
$C_{a}$, respectively. This restriction is to emulate a realistic bound
appropriate in a clinical setting to avoid toxicity effects of potential
overdosing. As will be seen, these maximum bounds are violated when state
constraints are not included in the problem statement. We next provide several
numerical solutions to the optimal control problem with the $L^{1}$ type
objective and discuss issues regarding when singular controls arise as well as
the case necessitating an explicit inclusion of state constraints. The
following weights appearing in the objective function 9 for the $L^{1}$ case
will be used for the various constraint scenarios explored:
$(b_{1},b_{2})=(1,50)$ and $(a_{1},a_{2},a_{3},a_{4})=(100,5,30,10)$. The
largest weights are assigned to the states $P$ and $D$ since eliminating the
pathogen and decreasing the level of damage is our goal.

####  4.1.1 Numerical results for $L^{1}$ objective with dosing constraints
$\,0\leq u_{p}\leq 0.5$ and $0\leq u_{a}\leq 0.62$

In this scenario we assume the upper bounds on the dosing input to be the same
as those used in [12]: $\,0\leq u_{p}\leq 0.5$ and $0\leq u_{a}\leq 0.62$.
Figure 3 displays the simulation results of this scenarios, showing the
optimal states of Pathogen ($P$), the pro-inflammatory state ($N$), the level
of Damage ($D$) and the anti-inflammatory state ($C_{a}$) in Figure 3 (a)-(d)
and the resulting control doses, $u_{p}$ and $u_{a}$, in panels (e)-(f). The
pro-inflammatory control $u_{p}$ is provided at maximum level for about half
an hour, while the anti-inflammatory dose $u_{a}$ is followed just after for
about one hour. Although not very evident in Figure 3, $u_{a}$ is applied with
one sample difference after $u_{p}$.

When a switching function is not zero there is “bang-bang” control and,
although the controls are “bang-bang” for most of the simulation, $u_{a}$
actually portrays a singular control between time
$(t_{1},t_{2})=(2.354,23.88)$. The scaled switching function $\phi_{2}$ shown
in Figure 3(h) is zero during this interval, unlike $\phi_{1}$ in panel (g)
which does not vanish in any finite open interval. We first demonstrate that
the singular control is minimizing and then demonstrate that an approximation
to the singular control gives similar results and provides a more practical
dosing protocol to implement.

Figure 3: (a)-(d) Time courses for optimal states, (e)-(f) control states, and
(g)-(h) the switching functions $\phi_{1}$ and $\phi_{2}$. Control constraints
in this simulation which uses an $L^{1}$ type objective are $0\leq u_{p}\leq
0.5,$ $0\leq u_{a}\leq 0.62$.

To determine if the singular control generated is minimizing, we derive a
solution corresponding to the singular control. However, the fact that there
is no explicit dependence on the control in either switching function (19) or
(20) makes it problematic to determine where they vanish and higher order
derivatives must be examined in order to derive a solution that will
correspond to the singular control. As mentioned, singular control occurs in
the anti-inflammatory control $u_{a}$ and thus, the switching function
corresponding to $\phi_{2}$ will be of interest. The first derivative of
$\phi_{2}$ is given in Eq.(4.1.1); however, since the control does not appear
in the expression, a second derivative is needed.444If $\phi_{2}$ vanishes on
$(t_{1}\,,t_{2})$ then we must have
$\phi_{2}=\dot{\phi}_{2}=\ddot{\phi}_{2}=0$ until the first control appears

$\displaystyle\dot{\phi_{2}}=-a_{4}-2\lambda_{1}\frac{zPNC_{a}}{(g+C_{a}^{2})^{2}}+\lambda_{2}\Big{[}\frac{2nC_{a}(hP+iN+jD)}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}\Big{]}$
$\displaystyle+\lambda_{3}\Big{[}\frac{12qsN^{6}(sC_{a}^{2}+r)^{5}}{((sC_{a}^{2}+r)^{6}+tN^{6})^{2}}\Big{]}-\lambda_{4}\Big{[}-\frac{2(uN+vD)C_{a}}{(gN+wD+C_{a}^{2}+g)^{2}}-k_{c}\Big{]}=0$
(23)

The second derivative of $\phi_{2}$ is given by Equation (24) where the
control now appears, and which $\phi_{2}=0$ can be solved to derive the
solution corresponding to the singular control.

$\displaystyle\ddot{\phi}_{2}$
$\displaystyle=\left(-\frac{2\lambda_{1}zNC_{a}}{(C_{a}^{2}+g)^{2}}-\frac{2b_{1}nC_{a}h}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}+\frac{4b_{1}nC_{a}(hP+iN+jD)k}{(nC_{a}^{2}+kP+lN+mD+o)^{3}}\right)\dot{P}$
$\displaystyle+\theta_{1}\left(\frac{(hP+iN+jD)}{(nC_{a}^{2}+kP+lN+mD+o)}-pN\right)+\left[\frac{-2b_{1}nC_{a}j}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}\right.$
$\displaystyle\left.+\frac{4b_{1}nC_{a}(hP+iN+jD)m}{(nC_{a}^{2}+kP+lN+mD+o)^{3}}+b_{2}\left(\frac{2vC_{a}}{(gN+wD+C_{a}^{2}+g)^{2}}-\frac{2(2uN+2vD)C_{a}w}{(gN+wD+C_{a}^{2}+g)^{3}}\right)\right]\dot{D}$
$\displaystyle+\theta_{2}\left(\frac{sc+(uN+vD)}{(gN+wD+C_{a}^{2}+g)}-k_{c}C_{a}\right)-2\dot{\lambda}_{1}\frac{zPNC_{a}}{(g+C_{a}^{2})^{2}}+\dot{\lambda}_{2}\Big{[}\frac{2nC_{a}(hP+iN+jD)}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}\Big{]}$
$\displaystyle+\theta_{1}u_{p}+\theta_{2}u_{a}=0,$ (24)

where

$\displaystyle\theta_{1}$
$\displaystyle=-2\frac{\lambda_{1}zPC_{a}}{(C_{a}^{2}+g)^{2}}-\frac{2b_{1}niC_{a}}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}+\frac{4b_{1}nC_{a}(hP+iN+jD)l}{(nC_{a}^{2}+kP+lN+mD+o)^{3}}$
$\displaystyle+\frac{72qsN^{5}(sC_{a}^{2}+r)^{5}}{((sC_{a}^{2}+r)^{6}+tN^{6})^{2}}\lambda_{3}-\frac{144\lambda_{3}qsN^{11}(sC_{a}^{2}+r)^{5}t}{((sC_{a}^{2}+r)^{6}+tN^{6})^{3}}$
$\displaystyle+b_{2}\left(\frac{2uC_{a}}{(gN+wD+C_{a}^{2}+g)^{2}}-\frac{2(2uN+2vD))C_{a}g}{(gN+wD+C_{a}^{2}+g)^{3}}\right)$
(25)

and

$\displaystyle\theta_{2}$
$\displaystyle=8\frac{\lambda_{1}zPNC_{a}^{2}}{(C_{a}^{2}+g)^{3}}-2\frac{\lambda_{1}fPN}{(C_{a}^{2}+g)^{2}}-8\frac{\lambda_{2}(hP+iN+jD)n^{2}C_{a}^{2}}{(nC_{a}^{2}+kP+lN+mD+o)^{3}}$
$\displaystyle+2\frac{\lambda_{2}(hP+iN+jD)n}{(nC_{a}^{2}+kP+lN+mD+o)^{2}}+120\lambda_{3}\frac{qs^{2}N^{6}(sC_{a}^{2}+r)^{4}C_{a}}{((sC_{a}^{2}+r)^{6}+tN^{6})^{2}}$
$\displaystyle-\lambda_{4}\left(\frac{-(2uN+2vD)}{(gN+wD+C_{a}^{2}+g)^{2}}+\frac{(4(2uN+2vD))C_{a}^{2}}{(gN+w*D+C_{a}^{2}+g)^{3}}\right)-288\lambda_{3}\frac{qs^{2}N^{6}(sC_{a}^{2}+r)^{10}C_{a}}{((sC_{a}^{2}+r)^{6}+tN^{6})^{3}}.$
(26)

For singular control to be minimizing, the generalized Legendre-Clebsch
condition shown in (27) needs to be satisfied [37].

$\displaystyle(-1)\frac{\partial}{\partial
u_{a}}\frac{d^{2}}{dt^{2}}\frac{\partial H}{\partial u_{a}}\geq 0$ (27)

That is, during the simulation, the expression in $-\theta_{2}$ (See Figure
3(h)) must be positive semi-definite. Although not reported here, a plot of
$-\theta_{2}$ (See Eq.(26)) shows that it is strictly positive for
$t\in(t_{1},t_{2})$, where $t_{1}=2.354$ and $t_{2}=23.88$, proving the
optimality of the singular arc. If we consider the pro-inflammatory input,
$u_{p}$, to be constant then the singular control part of $u_{a}$ can be
easily derived from Eq.(24). Note that in all simulations, $u_{p}$ is zero
during the interval where a singular control is present in $u_{a}$.

Figure 3(a)-(f) displays results of the simulation. The cost function
$J_{\ref{sec:L1dosing}}$ associated with Figure 3, the entry and exit times of
the singular control, and the final states corresponding to the simulation are
given as follows:

$\displaystyle J_{\ref{sec:L1dosing}}$ $\displaystyle=2104.78$ $\displaystyle
t_{1}$ $\displaystyle=2.354$ $\displaystyle t_{2}$ $\displaystyle=23.88$
$\displaystyle P(t_{f})$ $\displaystyle=0$ $\displaystyle N(t_{f})$
$\displaystyle=0.0061$ $\displaystyle D(t_{f})$ $\displaystyle=0.025$
$\displaystyle C_{a}(t_{f})$ $\displaystyle=0.29$

Finally, we consider an approximation to the “bang-bang-singular” control (or
what we term the basic optimal control), where we replace the time course for
the singular control with a constant input, namely zero. In other words, we do
not provide any dosing of $u_{a}$ during the period where the singular control
occurred and determine if this approximation is acceptable in terms of the
resulting outcomes of the states. The motivation of replacing the singular
control with a constant is to suggest a more practical implementation of the
dosing protocol.

Figure 4: Comparison of basic optimal control with an approximate control for
the $L^{1}$ objective case with dosing constraints of $\,0\leq u_{p}\leq 0.5$
and $0\leq u_{a}\leq 0.62$. (a) Generated optimal controls $u_{p}$ and $u_{a}$
and (b) corresponding optimal states. (c) Generated optimal control $u_{p}$
and approximation of optimal control $u_{a}$ and (d) corresponding states. The
time frame displayed in both plots (a) and (c) of the control doses is for
$t_{f}=40$ hours since afterward the values of both controls are zero.

Figure 4 provides the comparison of the simulation using the optimal controls
generated for $u_{p}$ and $u_{a}$ (panels (a)-(b) which are the same as in
Figure 3(a)-(f)) versus the simulation where an approximation of the $u_{a}$
control was made for the singular arc (panels (c)-(d)). Comparison of Figure 4
panel (d) with (b) shows that the states $C_{a}$ and $D$ are the ones mainly
affected by this approximate control case. In other words, changing the level
of anti-inflammatory dose $C_{a}$, namely reducing it to zero for that time
period, has a direct impact on the mitigation of damage ($D$). Thus, a lack of
anti-inflammatory mediator at critical moments results in an increase of the
level of damage, which is seen in Figure 4(d). In this particular case, the
long-term dynamics of the two simulations are the same; however, such
approximations would need to be done carefully. An improved solution might be
to approximate the singular arc by a piece-wise constant function of say, four
pieces which would better capture the gradual decrease of $u_{a}$ occurring in
the singular arc.

####  4.1.2 Numerical results for $L^{1}$ objective with dosing constraints
$0\leq u_{p}\leq 0.25$ and $0\leq u_{a}\leq 0.31$

We next explored how a reduction of the maximum allowable doses affected the
ability of the control problem with $L^{1}$ objective to drive the patient to
its healthy state. Thus the upper bounds on $u_{p}$ and $u_{a}$ were reduced
to $0\leq u_{p}\leq 0.25$ and $0\leq u_{a}\leq 0.31$, which is half of the
amount used in the previous section. The resulting simulation is given in
Figure 5 (a)-(b).

Figure 5: Comparison of basic optimal control with an approximate control for
the $L^{1}$ objective case with upper bounds on dosing constraints given as
$u_{pmax}=0.25$ and $u_{amax}=0.31$. (a) Generated optimal controls $u_{p}$
and $u_{a}$ and (b) corresponding optimal states. (c) Generated optimal
control $u_{p}$ and approximation of optimal control $u_{a}$ and (d)
corresponding states. The time frame displayed in both plots (a) and (c) of
the control doses is for $t_{f}=40$ hours since afterward the values of both
controls are zero.

The numerical information about the cost, exit and entry times for the
singular arc, as well as the final states corresponding to Figure 5 are given
by

$\displaystyle J_{\ref{sec:L1halvedcon}}$ $\displaystyle=2469.904$
$\displaystyle t_{1}$ $\displaystyle=3.363$ $\displaystyle t_{2}$
$\displaystyle=21.02$ $\displaystyle P(t_{f})$ $\displaystyle=0$
$\displaystyle N(t_{f})$ $\displaystyle=0.0074$ $\displaystyle D(t_{f})$
$\displaystyle=0.0326$ $\displaystyle C_{a}(t_{f})$ $\displaystyle=0.311$

The singular control verifies the Legendre-Clebsch condition in (27) for
optimality. The simulation results show an increase in the cost function
$J_{\ref{sec:L1halvedcon}}$ when compared to $J_{\ref{sec:L1dosing}}$ of the
previous section and this can be explained by the increase of the level of
Damage, $D$, as a result of decreasing the total dose of the anti-inflammatory
mediator, $C_{a}$. Similar to what was done in the previous section, Figure 5
(d) displays the outcome of the states when an approximation of the optimal
control input, $u_{a}$, is made as shown in panel (c).

In either sets of simulations for this and previous section and in either the
optimal or approximate cases, we observe that the level of damage ($D$)
remains such that $D\leq 1$ (and can be verified to approach zero eventually
in the absence of further control) and the level of pathogen ($P$) decreases
relatively quickly to zero. An interesting pattern we observe is that each
dose associated with an increase of the pro-inflammatory mediator is followed
immediately after by a dose of anti-inflammatory therapy. This behavior can be
explained by the fact that the immune system requires an initial boost of
activated phagocytes, $N$, to eliminate the pathogen threat, but the resulting
inflammation causes some self-harm seen in the relative increase of tissue
damage, $D$, which then needs the inhibitory effects of $C_{A}$ via an anti-
inflammatory dose, $u_{a}$, to ensure the decay of damage.

The anti-inflammatory dosing input is applied for a period of time clearly
longer than $u_{p}$. It is interesting to notice in Figure 5 that when the
maximum allowable doses were decreased by half, the control doses seems to be
applied for a longer duration when compared to the initial scenario where
$N_{max}=0.5$ and $C_{amax}=0.62$. As a consequence of lowering the intensity
of the dose, the optimal control compensates by allowing a longer application.
This behavior was also observed when solving the TPBVP for the $L^{1}$ problem
using the forward backward sweep method in [22].

It can also be seen in Figure 4 of the previous section that the anti-
inflammatory states do not remain below the upper bound mentioned previously
on the anti-inflammatory state $C_{a}$ (i.e. $C_{amax}=0.62$555The max value
for $C_{a}$ was derived from the bifurcation study in [19] of the $N/D$
subsystem. It was determined in that study that the system loses bi-stability
between the healthy and the death state if $C_{a}>0.62$ with only the healthy
state being stable. For this reason, to ensure the presence of a stable death
state the latter bound on $C_{a}$ is enforced. The max value of $N=0.5$ was
determined experimentally (through simulations; See [12])) and the same is
used for the simulations of this section given in Figure 5, even though the
maximum upper bound on the doses for both controls was reduced by half. Thus,
in the next section, we impose state constraints on $N$ and $C_{a}$ in the
problem formulation.

####  4.1.3 Numerical results for $L^{1}$ objective with state and dosing
constraints

In this section, we now impose the state constraints on $N$ and $C_{a}$ we
defined earlier of $0\leq N(t)\leq 0.5$ and $0\leq C_{a}\leq 0.62$, along with
dosing constraints of $0\leq u_{p}\leq 0.5$ and $0\leq u_{a}\leq 0.62$ which
are the same dosing constraints of Section 4.1.2. The analysis of the
numerical results obtained when using the state constraints together with the
control constraints would be subject, in most part, to the same discussion as
in the section before; the only difference would be the additional jump
conditions on the adjoints and the transversality condition. To be more
pertinent, we will refer the reader to the similar discussion presented later
in section 4.2.2. For the sake of brevity, the details are not provided for
this case.

The resulting simulation is given in Figure 6 and the numerical information
about the costs and final states corresponding to Figure 6 are given by

$\displaystyle J_{\ref{sec:L2bothcon}}$ $\displaystyle=2393.339$
$\displaystyle t_{1}$ $\displaystyle=2.186$ $\displaystyle t_{2}$
$\displaystyle=4.036$ $\displaystyle P(t_{f})$ $\displaystyle=0.0000$
$\displaystyle N(t_{f})$ $\displaystyle=0.0073$ $\displaystyle D(t_{f})$
$\displaystyle=0.0318$ $\displaystyle C_{a}(t_{f})$ $\displaystyle=0.3093$

Figure 6:  Optimal doses and corresponding state trajectories for an $L^{1}$
type objective with state constraints $0\leq N(t)\leq 0.5$, $0\leq C_{a}\leq
0.62$ together with dosing constraints $0\leq u_{p}\leq 0.5$ and $0\leq
u_{a}\leq 0.62$. (a) Dosing profiles for $u_{p}$ and $u_{a}$ controls that
lead to (b) a healthy outcome as all states decrease toward background levels.
The time frame displayed in (a) of the control doses is for $t_{f}=40$ hours
since afterward the values of both controls are zero.

The simulation in Figure 6(a) shows that a singular arc occurs between time
$t_{1}=2.186$ and $t_{2}=4.036$, which can be verified as before by checking
the resulting function $\phi_{2}$, as the one given in Figure 3(h). The
optimality of the singular arc can be derived by verifying the Legendre-Clebsh
condition of (27). Notice that 6(a) shows that the pro-inflammatory dose is
applied for a duration of almost 20 minutes followed by a rest period of 50
minutes before the anti-inflammatory dose is injected for a longer duration
corresponding to one hour.

Thus far, we have studied optimal control scenarios in which the objective is
linear in the control. In the remainder of this paper the focus will be on
deriving the optimal control doses subject to an objective that is quadratic
in the control as well as different constraint types. A discussion comparing
the results of the uses of the different objective functions will also be
given.

###  4.2 Optimal Control with $L^{2}$ type Objective

In this section we present numerical simulations of the optimal control
problem when the objective is quadratic in the control.

$\displaystyle
J_{L^{2}}=x^{T}(t_{f})x(t_{f})+\int_{0}^{t_{f}}\Big{[}x(t)^{T}Ax(t)+u(t)^{T}Bu(t)\Big{]}dt$
(28)

where A and B are matrices with diagonals $(a_{1},a_{2},a_{3},a_{4})$ and
$B=(b_{1},b_{2})$ respectively. Also,
$x(t)=[P(t),\,N(t),\,D(t),\,C_{a}]^{T}\geq 0$ and $u(t)=[u_{p},u_{a}]^{T}\geq
0$, where $\geq$ means componentwise.

We also once again define $N_{max}=0.5$ and $C_{a_{max}}=0.62$ as the maximum
allowed levels of the variables $N$ and $C_{a}$, respectively. We provide
several numerical solutions to the optimal control problem with the $L^{2}$
type objective under various constraint definitions. The following weights
appearing in the objective function (28) for the $L^{2}$ case will be used for
the various constraint scenarios explored in the subsections that follow:
$(b_{1},b_{2})=(1,20)$ and $(a_{1},a_{2},a_{3},a_{4})=(10,1,10,10)$. Note that
the objective function includes the cost at the final time,
$\phi(x(t_{f})^{2})$. All other aspects of the numerical settings are as given
in the beginning of Section 4.

The Hamiltonian for the $L^{2}$ type objective is given according to

$\displaystyle
H(x,u,\lambda)={x(t)}^{T}Ax(t)+{u(t)}^{T}Bu(t)+\lambda(t)^{T}(f(x)+g(x)u)$
(29)

where $\lambda(t)=(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4})^{T}$.
While The $L^{1}$ type objective leads in general to more intricate analysis
involving bang-bang and singular control, the $L^{2}$ type objective analysis
often results in continuous control, when no state constraints are included.
The case of mixed state-control inequality constraint will be introduced in
next section, followed by a section that looks at separate or _pure_ control
and state constraints. For all constraints used in the following subsections,
those associated with $N$ and/or $u_{p}$ have a maximum bound of $0.5$, while
those associated with $C_{a}$ and/or $u_{a}$ have a maximum bound of $0.62$.

####  4.2.1 Numerical results for $L^{2}$ objective with mixed control-state
inequality constraints

Consider the mixed control-state inequality constraints 666Note that we adhere
to the sign convention in [23] given according to

$\displaystyle C_{1}(N(t),u_{p}(t))$ $\displaystyle=u_{p}(t)+N(t)-0.5\leq 0$
(30) $\displaystyle C_{2}(C_{a}(t),u_{a}(t))$
$\displaystyle=u_{a}(t)+C_{a}(t)-0.62\leq 0.$ (31)

Imposing mixed control-sate constraints allows the derived optimal control to
account for the actual value of the states, and hence only the necessary
control dose will be administered. The simulation results for this OCP are
given in Figure 7 and correspond to the numerical information about the costs
and final states given as follows:

$\displaystyle J_{\ref{sec:L2mixed}}$ $\displaystyle=602.67$ $\displaystyle
P(t_{f})$ $\displaystyle=0$ $\displaystyle N(t_{f})$ $\displaystyle=0.0082$
$\displaystyle D(t_{f})$ $\displaystyle=0.0376$ $\displaystyle C_{a}(t_{f})$
$\displaystyle=0.3216$ (32)

Figure 7: (a) Optimal doses and (b) corresponding state trajectories for a
quadratic type objective with mixed control state inequality constraints
$0\leq u_{p}(t)\leq 0.5-N(t)$, $0\leq u_{a}(t)\leq 0.62-C_{a}(t)$. The time
frame displayed in (a) of the control doses is for $t_{f}=40$ hours since
afterward the values of both controls are zero.

Note that both controls are constrained to be positive semidefinite, i.e.,
$0\leq u_{p}$ and $0\leq u_{a}$. All the constraints are adjoined such that
the augmented Hamiltonian in equation (29) can be expressed according to

$\displaystyle\mathcal{H}(x,u,\lambda,\mu)$
$\displaystyle=H(x,u,\lambda)+\mu_{1}C_{1}(N,u_{p})+\mu_{2}C_{2}(C_{a},u_{a})-\mu_{3}u_{p}-\mu_{4}u_{a}.$
(33)

Necessary conditions for minimizing the Hamiltonian can then be derived. The
multiplier $\mu$ is given such that it verifies the complementary condition

$\displaystyle\mu_{i}=\begin{cases}0&if\quad C_{i}<0\\\ \mu_{i}\geq 0&if\quad
C_{i}=0\quad for\quad i=1,2.\end{cases}$ (34)

The adjoint equations are given in the following compact form

$\displaystyle\lambda^{T}=-\mathcal{H}_{x}=\begin{cases}-H_{x}(x,u,\lambda)&if\quad
C<0\\\ -H_{x}(x,u,\lambda)-\mu C_{x}&if\quad C=0;\end{cases}$ (35)

and more explicitly as follows:

$\displaystyle\dot{\lambda}_{1}$ $\displaystyle=-\frac{\partial H}{\partial
P}$ $\displaystyle\dot{\lambda}_{3}$ $\displaystyle=-\frac{\partial
H}{\partial D}$ (36) $\displaystyle\dot{\lambda}_{2}$
$\displaystyle=-\frac{\partial H}{\partial N}-\mu_{1}$
$\displaystyle\dot{\lambda}_{4}$ $\displaystyle=-\frac{\partial H}{\partial
C_{a}}-\mu_{2}.$ (37)

The controls $u_{p}$ or $u_{a}$ can be derived from the mixed constraints
along the constrained arc. In other words, whenever

$\displaystyle C(x,u)$ $\displaystyle=0$ $\displaystyle for\quad
t\in(t_{1},t_{2}).$ (38)

Hence, the boundary controls are given by

$\displaystyle u_{p_{b}}$ $\displaystyle=0.5-N(t)$ (39) $\displaystyle
u_{a_{b}}$ $\displaystyle=0.62-C_{a}(t).$ (40)

The free controls can be derived from the minimum condition such that

$\displaystyle u_{p}$ $\displaystyle=-\frac{\lambda_{2}}{2b_{1}}$ (41)
$\displaystyle u_{a}$ $\displaystyle=-\frac{\lambda_{4}}{2b_{2}}.$ (42)

Since the regularity condition is verified, i.e., $C_{u}(x,u)\neq 0$, the
multiplier $\mu$ follows from the local minimum condition
$\mu=\mu(x,\lambda)=-H_{u}(x,\lambda,u_{b}(x))/C_{u}(x,u_{b}(x))$. In other
words, by

$\displaystyle\mu_{1}$ $\displaystyle=-2b_{1}u_{p}-\lambda_{2}$ (43)
$\displaystyle\mu_{2}$ $\displaystyle=-2b_{2}u_{a}-\lambda_{4},$ (44)

whenever the lower constraints on the controls are not active (i.e. $u_{p}>0$
and $u_{a}>0$). When the mixed constraints are active, the multipliers
$\mu_{1}$ and $\mu_{2}$ are given by substituting the equation for $u_{p_{b}}$
given in (39) into (43) and similarly for $u_{a_{b}}$ of (40) into (44), such
that

$\displaystyle\mu_{1}$ $\displaystyle=-2b_{1}(0.5-N(t))-\lambda_{2}$ (45)
$\displaystyle\mu_{2}$ $\displaystyle=-2b_{2}(0.62-C_{a}(t))-\lambda_{4}.$
(46)

The optimal controls are given as follows

$\displaystyle
u_{p}(t)=\begin{cases}-\frac{1}{2b_{1}}(\lambda_{2}+\mu_{1})&if\quad 0\leq
t\leq t_{0}=1.009\\\ 0&if\quad t_{0}\leq t\leq t_{f}=168;\end{cases}$ (47)
$\displaystyle u_{a}(t)=\begin{cases}0,&if\quad 0\leq t\leq t_{1}=1.009\\\
\frac{-1}{2b_{2}}(\lambda_{4}+\mu_{2}),&if\quad t_{1}\leq t\leq t_{2}=19\\\
0,&if\quad t_{2}\leq t\leq t_{f}=168.\end{cases}$ (48)

Note that the mixed control-state constraints, $C_{1}$ and $C_{2}$ in
equations (30) and (31), respectively, are active whenever the controls are
not zero; that is, $\mu_{1}\geq 0$ and $\mu_{2}\geq 0$. For $0\leq t\leq
t_{0}$ and $t_{1}\leq t\leq t_{2}$ the switching structure is given by the
junction conditions

$\displaystyle\mu_{1}(t)$ $\displaystyle=-2b_{1}(0.5-N(t))-\lambda_{2}(t)=0$
(49) $\displaystyle\mu_{2}(t)$
$\displaystyle=-2b_{2}(0.62-C_{a}(t))-\lambda_{4}(t)=0.$ (50)

The next section presents a scenario where we have pure state constraints on
$N$ and $C_{a}$ as well as constraints on the upper bounds of the control
doses.

####  4.2.2 Numerical results for $L^{2}$ objective with pure state and
control constraints

Now consider imposing path constraints on the states $N$ and $C_{a}$ such that
$N(t)\leq 0.5$ and $C_{a}(t)\leq 0.62$ $\forall t$. In addition, we will
require as before that the following control constraints be met: $0\leq
u_{p}\leq 0.5$ and $0\leq u_{a}\leq 0.62$. The simulation results for this OCP
are given in Figure 8 and correspond to the numerical information about the
costs and final states given as follows:

$\displaystyle J_{\ref{sec:L2pure}}$ $\displaystyle=513.77$ $\displaystyle
P(t_{f})$ $\displaystyle=0$ $\displaystyle N(t_{f})$ $\displaystyle=0.0071$
$\displaystyle D(t_{f})$ $\displaystyle=0.0306$ $\displaystyle C_{a}(t_{f})$
$\displaystyle=0.3065$

Figure 8:  (a) Optimal doses and (b) corresponding state trajectories for a
quadratic type objective with state path constraints $0\leq N(t)\leq 0.5$,
$0\leq C_{a}\leq 0.62$. The time frame displayed in (a) of the control doses
is for $t_{f}=40$ hours since afterward the values of both controls are zero.

Consider the last inequality in the (OCP) problem $S(x(t))\leq 0$, given
specifically for this case as

$\displaystyle S_{1}(N(t))$ $\displaystyle=N(t)-0.5\leq 0$ (51) $\displaystyle
S_{2}(C_{a}(t))$ $\displaystyle=C_{a}(t)-0.62\leq 0.$ (52)

Whenever one of the constraints becomes active on a subinterval we have

$\displaystyle S_{i}(x(t))$ $\displaystyle=0$ for all $\displaystyle
t\in(t_{1},t_{2})\subset[0,168].$ (53)

Let $r$ be the _order_ of the state constraint; that is, the smallest non-
negative numbers necessary for the control to first appear [38]. It is
straightforward to see that $r=1$ and furthermore,

$\displaystyle\dot{S}_{2}(x,u_{a})=\dot{C}_{a}=s_{c}+\frac{uN+vD}{g+C_{a}^{2}+gN+wD}-\mu_{c}C_{a}+u_{a}(t).$
(54)

The state constraint satisfies the following regularity condition

$\displaystyle\frac{\partial}{\partial
u_{a}}\dot{S}_{2}(x,u_{a})=(0,1)\neq(0,0).$ (55)

######  Remark 4.3.

We can derive similar conditions for $S_{1}$. We will focus on the state
$C_{a}$ since that is the only state that becomes active.

The boundary control $u_{a_{b}}$ on the constrained arc can be derived from
(54) by also replacing $C_{a}(t)=0.62$ such that

$\displaystyle u_{a_{b}}=-s_{c}-\frac{uN+vD}{g+0.62^{2}+gN+wD}+0.62\mu_{c}.$
(56)

The augmented Hamiltonian is given according to

$\displaystyle\mathcal{H}(x,u,\lambda,\lambda_{s})$
$\displaystyle={x(t)}^{T}Ax(t)+{u(t)}^{T}Bu(t)+\lambda(t)^{T}(f(x)+g(x)u)+\lambda_{s}S(x),$
(57)
$\displaystyle+\mu_{1}(u_{p}(t)-0.5)+\mu_{2}(u_{a}(t)-0.62)-\mu_{3}u_{p}-\mu_{4}u_{a},$
(58)

where $\lambda_{s}=[\lambda_{sp},\,\lambda_{sa}]$ is a row vector and
$S(x)=[S_{1},\,S_{2}]^{T}$ is a column vector. The adjoint equation can be
formulated in a similar fashion as in (35). The multiplier $\lambda_{sa}$
needs to verify the complementary slackness condition [38] given as:

$\displaystyle\lambda_{sa}=\begin{cases}0&if\quad S_{2}<0\\\ \lambda_{sa}\geq
0&if\quad S_{2}=0.\end{cases}$ (59)

The adjoint equation corresponding to the dynamic equation $\dot{C_{a}}$ is
given by

$\displaystyle\dot{\lambda}_{4}(t)=-\frac{\partial\mathcal{H}}{\partial
C_{a}}-\lambda_{sa}.$ (60)

The minimum condition can be derived as in the previous section when
$0=\frac{\partial\mathcal{H}}{\partial u_{p}}=2b_{2}u_{a}+\lambda_{4}$. On the
boundary arc, $u_{a}$ is constant (See Figure 9), and we get
$\dot{\lambda}_{4}=0$ as can be seen in the same figure. The multipliers
$\lambda_{sa}$ can be derived from (60) according to

$\displaystyle\lambda_{sa}=-\frac{\partial\mathcal{H}}{\partial
C_{a}}\Bigg{|}_{C_{a}=0.62}.$ (61)

The optimal controls are given as follows

$\displaystyle u_{p}(t)=\begin{cases}0.5&if\quad 0\leq t\leq t_{0}=0.3363\\\
0&if\quad t_{0}\leq t\leq t_{f}=168;\end{cases}$ (62) $\displaystyle
u_{a}(t)=\begin{cases}0,&if\quad 0\leq t\leq t_{1}=1.009\\\ 0.62,&if\quad
t_{1}\leq t\leq t_{2}=2.186\\\ u_{a_{b}},&if\quad t_{2}\leq t\leq
t_{3}=17.32\\\ \frac{-1}{2b_{2}}(\lambda_{4}),&if\quad t_{3}\leq t\leq
t_{4}=17.66\\\ 0,&if\quad t_{4}\leq t\leq t_{f}=168.\end{cases}$ (63)

The jump conditions and transversality condition are given, respectively, as

$\displaystyle\lambda(t_{s^{-}})$
$\displaystyle=\lambda(t_{s^{+}})+\beta_{s}S_{x}(x(t_{s}))\quad\beta_{s}\geq
0\textrm{ and }$ (64) $\displaystyle\lambda(t_{f})$
$\displaystyle=(2P(t_{f}),2N(t_{f}),2D(t_{f}),2C_{a}(t_{f})).$ (65)

Figure 9: Numerical solutions of the state and control constrained problem
with $L^{2}$ objective. Top row: control $u_{a}(t)$ and constrained anti-
inflammatory state $C_{a}(t)$. Bottom row: The scaled adjoint $\lambda_{4}$
corresponding to the state $C_{a}(t)$ together with the multiplier
$\lambda_{sa}$ corresponding to the constraint $S(x(t))=C_{a}(t)-0.62$. The
time frame displayed in the figures is for $t_{f}=40$ hours since afterward
the values of the control inputs are zero.

Figure 9 displays the computed anti-inflammatory dose $u_{a}$ and the
evolution of the constrained state $C_{a}$ together with the adjoint variable
$\lambda_{4}$ and the multiplier $\lambda_{s}$ corresponding to the state
constraint $S(x(t))=C_{a}(t)-0.62$. The Figure outlines that the constraint on
$C_{a}$ becomes active in the interval $[2.186,17.32]$. One can also see in
the same figure that the multiplier $\lambda_{s}$ satisfies the complementary
condition (59). The jump discontinuity at $t_{1}=2.186$ on the control $u_{a}$
for this $L^{2}$ optimal control is nonintuitive since one would expect a
continuous control as in the scenario of the previous section; however, this
result is common when using pure state constraints (See Eq. (64)). The scaled
adjoint variable $\lambda_{4}$ in Figure 9, which is associated with the
differential equation $\dot{x}_{4}$, has a jump discontinuity at the time the
constraint becomes active.

######  Remark 4.4.

The objective functional value corresponding to section 4.2.1, which is
associated with the mixed control-state constraints is higher than the one
given in section 4.2.2, that is associated with pure control and state
constraints. In the case that imposes mixed control-state inequality
constraints, the maximum allowed doses are updated based on the the actual
states of $C_{a}$ and $N$. For instance, if at a sampling time, $t$,
$C_{a}(t)=0.62$ the upper bound on $u_{a}$ is zero and similarly, if
$N(t)=0.5$, then the upper bound on $u_{p}$ is zero. However, this constraint
formulation is more restrictive with the amount of doses given which results
in an increase in the damage response (Figure 7) compared with the pure state
and control constraints (Figure 8), which may require additional or longer
dosing to achieve a successful outcome.

##  5 Conclusion

In this work, we determined successful treatment protocols to regulate the
inflammatory response modeled with an ordinary differential equations system
through solving the optimal control problem using a direct approach. Two cases
were considered regarding the form of the objective functional: linear versus
quadratic in the control. The results of the direct numerical methods used are
in agreement with the characterization of the optimal control. A noticeable
difference was not observed in the way the controls are provided when
considering an $L^{1}$ versus $L^{2}$ type objective. In either case, it was
observed that state constraints needed to be included explicitly in order for
the states not to exceed the prescribed bounds set for realistic
considerations. The case where singular control appears is also addressed when
the objective is linear in the control, the optimality of the singular arcs
was demonstrated by verifying the Legendre-Clebsch condition.

The dynamics of inflammatory immune response is more complex and varied than
represented by the reduced inflammatory model considered in this study,
however the obtained computational results can provide insight into possible
treatment strategies and the methods could prove to be useful tools to
incorporate in future practice. In particular, one can see a prevailing
pattern in the way the doses are implemented in this study and is in agreement
with results of previous computational studies ([12], [13]). Whether
considering the solution for the $L^{2}$ or $L^{1}$ type objective, generally
speaking, a successful control strategy is marked by an initial bolus of a
pro-inflammatory control ($u_{p}$) administered first to quickly eliminate the
harmful effect of the pathogen, followed just after by a two phase application
of the anti-inflammatory mediator ($u_{a}$). In the first phase, a bolus of
anti-inflammation is applied, followed by a longer duration application of
$u_{a}$, but with a very small amount. The similar patterns in the generated
control doses $u_{p}$ and $u_{a}$ when solving for either type of objective
seems clinically reasonable when thinking that in the presence of a pathogen,
the pro-inflammatory mediators will naturally increase to try to eliminate the
pathogen, which in turn will cause a certain increase in the level of
collateral tissue damage which incites further inflammation. If not
controlled, this positive feedback may result in the death of the patient.
Fortunately, however, the response also includes a release of anti-
inflammatory mediators which helps alleviate the harmful effect of the pro-
inflammation, ultimately driving the patient to its healthy equilibrium in the
absence of further pathogen threat. So essentially, the dynamic necessary for
stabilizing a patient with sever inflammation is reflected exactly by
increasing the pro-inflammatory effects with input $u_{p}$ followed just after
by a dose of anti-inflammatory mediator, with input $u_{a}$. It is important
to notice in our study that we assume the dosages represent the actual
concentrations with idealized effect on patient; that is, we do not assume any
model for the pharmacokinetics or pharmacodynamics. The effects of these two
phenomena on optimal protocols could be included in further study. Future work
would also benefit from the use of an immune response model that displays
similar complex behavior as the reduced model but does not lump too many
variables together to see how the obtained results would then relate to what
has been presented in this work.

## References

  * [1] J. Cohen, The immunopathogenesis of sepsis, Nature 420 (6917) (2002) 885–891.
  * [2] U. Food, D. A. (FDA)”, Fda drug safety communication: Voluntary market withdrawal of xigris [drotrecogin alfa (activated)] due to failure to show a survival benefit (Oct. 2011).
URL http://www.fda.gov/Drugs/DrugSafety/ucm277114.htm

  * [3] G. E. Grau, D. N. Maennel, Tnf inhibition and sepsis sounding a cautionary note, Nature Medicine 3 (1997) 1193–1195.
  * [4] K. Reinhart, W. Karzai, Anti-tumor necrosis factor therapy in sepsis: update on clinical trials and lessons learned, Crit Care Med 29 (7 Supplement 1) (2001) S121–S125.
  * [5] R. S. Parker, G. Clermont, Systems engineering medicine: engineering the inflammation response to infectious and traumatic challenges, Journal of The Royal Society Interface (2010) rsif20090517.
  * [6] B. W. Bequette, Challenges and recent progress in the development of a closed-loop artificial pancreas, Annual reviews in control 36 (2) (2012) 255–266.
  * [7] F. J. Doyle, L. M. Huyett, J. B. Lee, H. C. Zisser, E. Dassau, Closed-loop artificial pancreas systems: engineering the algorithms, Diabetes care 37 (5) (2014) 1191–1197.
  * [8] W. M. Haddad, T. Hayakawa, J. M. Bailey, Adaptive control for non-negative and compartmental dynamical systems with applications to general anesthesia, Int. J. Adapt. Control Signal Process 17 (2003) 209–235.
  * [9] R. F. Stengel, R. Ghigliazza, Stochastic optimal therapy for enhanced immune response, Mathematical biosciences 191 (2) (2004) 123–142.
  * [10] D. Kirschner, S. Lenhart, S. Serbin, Optimal control of the chemotherapy of HIV, Journal of mathematical biology 35 (7) (1997) 775–792.
  * [11] R. F. Stengel, R. Ghigliazza, N. Kulkarni, O. Laplace, Optimal control of innate immune response, Optimal Control Applications and Methods 23 (2) (2002) 91–104.
  * [12] J. Day, J. Rubin, G. Clermont, Using nonlinear model predictive control to find optimal therapeutic strategies to modulate inflammation, Mathematical Biosciences and Engineering (MBE) 7 (4) (2010) 739–763.
  * [13] G. Zitelli, S. M. Djouadi, J. D. Day, Combining robust state estimation with nonlinear model predictive control to regulate the acute inflammatory response to pathogen., Mathematical biosciences and engineering: MBE 12 (5) (2015) 1127–1139.
  * [14] R. Parker, E. Gatzke, F. Doye, Advanced model predictive control (mpc) for type i diabetic patient blood glucose control, in: American Control Conference, 2000\. Proceedings of the 2000, Vol. 5, IEEE, 2000, pp. 3483–3487.
  * [15] B. Grosman, E. Dassau, H. C. Zisser, L. Jovanovič, F. J. Doyle, Zone model predictive control: a strategy to minimize hyper-and hypoglycemic events, Journal of diabetes science and technology 4 (4) (2010) 961–975.
  * [16] H. Li, W. M. Haddad, Model predictive control for a multicompartment respiratory system, IEEE Transactions on Control Systems Technology 21 (5) (2013) 1988–1995.
  * [17] O. Bara, J. Day, S. M. Djouadi, Optimal control of an inflammatory immune response model, in: ”Proceedings of the IEEE 54th Annual Conference on Decision and Control (CDC), IEEE, 2015, pp. 1283–1288.  doi:10.1109/CDC.2015.7402388.
  * [18] O. Bara, S. Djouadi, J. Day, Immune therapy using optimal control with $l^{1}$\- type objective, in: Proceedings of the Americain Control Conference (ACC), IEEE, 2016, pp. 4895–4900.  doi:10.1109/ACC.2016.7526128.
  * [19] A. Reynolds, J. Rubin, G. Clermont, J. Day, Y. Vodovotz, G. B. Ermentrout, A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of anti-inflammation, Journal of theoretical biology 242 (1) (2006) 220–236.
  * [20] J. Day, J. Rubin, Y. Vodovotz, C. C. Chow, A. Reynolds, G. Clermont, A reduced mathematical model of the acute inflammatory response ii. capturing scenarios of repeated endotoxin administration, Journal of theoretical biology 242 (1) (2006) 237–256.
  * [21] L. Pontryagin, Mathematical Theory of Optimal Processes, Classics of Soviet Mathematics, Taylor & Francis, 1987.
  * [22] S. Lenhart, J. T. Workman, Optimal control applied to biological models, CRC Press, 2007.
  * [23] A. E. Bryson, Applied optimal control: optimization, estimation and control, CRC Press, 1975.
  * [24] U. Ledzewicz, H. Schättler, Optimal controls for a model with pharmacokinetics maximizing bone marrow in cancer chemotherapy, Mathematical biosciences 206 (2) (2007) 320–342.
  * [25] A. Swierniak, U. Ledzewicz, H. Schattler, Optimal control for a class of compartmental models in cancer chemotherapy, International Journal of Applied Mathematics and Computer Science 13 (3) (2003) 357–368.
  * [26] S. Subchan, R. Zbikowski, Computational optimal control: tools and practice, John Wiley & Sons, 2009.
  * [27] W. Fleming, R. Rishel, Deterministic and stochastic optimal control, Springer, 1975\.
  * [28] K. R. Fister, S. Lenhart, J. S. McNally, Optimizing chemotherapy in an hiv model, Electronic Journal of Differential Equations 1998 (32) (1998) 1–12.
  * [29] M. R. Kelly Jr, J. H. Tien, M. C. Eisenberg, S. Lenhart, The impact of spatial arrangements on epidemic disease dynamics and intervention strategies, Journal of biological dynamics 10 (1) (2016) 222–249.
  * [30] V. M. Becerra, Solving complex optimal control problems at no cost with psopt, in: Computer-Aided Control System Design (CACSD), 2010 IEEE International Symposium on, IEEE, 2010, pp. 1391–1396.
  * [31] A. Wächter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical programming 106 (1) (2006) 25–57.
  * [32] J. T. Betts, Practical methods for optimal control and estimation using nonlinear programming, Vol. 19, Siam, 2010.
  * [33] N. S. Bedrossian, S. Bhatt, W. Kang, I. M. Ross, Zero-propellant maneuver guidance, Control Systems, IEEE 29 (5) (2009) 53–73.
  * [34] I. M. Ross, F. Fahroo, Pseudospectral knotting methods for solving nonsmooth optimal control problems, Journal of Guidance, Control, and Dynamics 27 (3) (2004) 397–405.
  * [35] F. Clarke, M. De Pinho, Optimal control problems with mixed constraints, SIAM Journal on Control and Optimization 48 (7) (2010) 4500–4524.
  * [36] R. T. Rockafellar, State constraints in convex control problems of bolza, SIAM journal on Control 10 (4) (1972) 691–715.
  * [37] A. J. Krener, The high order maximal principle and its application to singular extremals, SIAM Journal on Control and Optimization 15 (2) (1977) 256–293.
  * [38] S. P. Sethi, G. L. Thompson, Optimal control theory: applications to management science and economics, Springer, 2005.