# Operator Splitting for Adaptive Radiation Therapy with Nonlinear Health
Dynamics

Anqi Fu  Lei Xing  Stephen Boyd

###### Abstract

We present an optimization-based approach to radiation treatment planning over
time. Our approach formulates treatment planning as an optimal control problem
with nonlinear patient health dynamics derived from the standard linear-
quadratic cell survival model. As the formulation is nonconvex, we propose a
method for obtaining an approximate solution by solving a sequence of convex
optimization problems. This method is fast, efficient, and robust to model
error, adapting readily to changes in the patient’s health between treatment
sessions. Moreover, we show that it can be combined with the operator
splitting method ADMM to produce an algorithm that is highly scalable and can
handle large clinical cases. We introduce an open-source Python implementation
of our algorithm, AdaRad, and demonstrate its performance on several examples.

##  1 Introduction

In radiation therapy, beams of ionizing radiation are transmitted into a
patient, damaging both tumor cells and normal tissue. The goal of radiation
treatment planning is to deliver enough dose to the tumor so that diseased
cells are killed, while avoiding excessive injury to the normal tissue and
organs-at-risk (OARs). This is achieved by optimizing the beam intensity
profile, or fluence map, subject to constraints on the dose to certain parts
of the patient’s anatomy. The fluence map optimization problem is well-studied
[RAD+03, RDL04, AGRD10, Gao16], and technology like intensity-modulated
radiation therapy (IMRT) is now widespread in the clinic [ZFH+00, WCP+06,
GAJ+12, Web10].

Treatment in practice usually takes place over multiple sessions. A clinician
will divide up the total prescribed dose into smaller dose fractions, which
are delivered over the course of several weeks or months. This permits normal
tissue time to recover and repair cell damage, but also gives tumors an
opportunity to proliferate, especially when the treatment course is long. A
study of 4,338 prostate cancer patients showed that biochemical failure
increases by 6% for every 1 week increase in treatment time, with a dose
equivalent of proliferation of 0.24 Gy/day [TKL+10]. Thus, an important
question in treatment planning is how to choose the sequence of deliverable
doses such that they balance these temporal effects on a patient’s health.

###  1.1 Related work

Early clinical practitioners split the prescribed dose equally over a fixed
number of sessions. While convenient, this method does not account for errors
or uncertainty in the treatment process. For example, due to movements of the
patient’s anatomy, the expected dose may differ from the actual dose to an
anatomical structure. If the actual dose is observable, a common way to
compensate for this is to divide the residual dose (i.e., the difference
between the prescribed and cumulative actual dose) across the remaining
sessions. This then becomes the new per-session dose goal. In [dlZAX07], the
authors solve for the beam intensities by minimizing the sum-of-squared
difference between this goal dose and the expected dose. They compare the
results when errors are perfectly known, so the expected dose is equal to the
actual dose, with the results when errors are assumed to be zero. A similar
approach is taken in [FV04], except the errors are modeled explicitly as a
random shift in the surrounding voxels. Instead of the dose to each voxel,
[SEP12] work with the equivalent uniform dose (EUD), a value that captures the
biological effect of a dose distribution over a region. Their objective is to
minimize the sum of the EUD over all treatment criteria subject to bounds on
the EUD of the tumor and normal tissues. To solve this problem, they employ
methods from approximate dynamic programming coupled with a discrete
probabilistic model of the dose error.

The papers we have discussed so far only focus on the dose to the patient. By
contrast, [KGP09] introduces a Markov decision process model that includes
both the dose (action) and the patient’s health state. Each choice of dose
induces a transition to a particular health state with some probability.
Making this idea concrete, [MTD+12] define the health of a tumor (resp. OAR)
to be the radiation (resp. damage) effect of the delivered dose, as calculated
from the linear-quadratic (LQ) model of cell survival [Fow89]. They analyze a
simple example with one tumor and one OAR and find that the optimal
fractionation scheme is either a single session delivery of the full dose or
equal dose fractions, depending on the relationship between the LQ parameters.
The authors of [BRTU15] extend this analysis to incorporate accelerated tumor
repopulation and show that the dose per session increases over the treatment
course. Using simulated annealing, [YX05] solve a similar treatment planning
problem based on the LQR model, which captures all 4 Rs (repair of sublethal
damage, repopulation, redistribution, and reoxygenation) of cellular radiation
response [BHH+95].

In [CM13], the authors take a probabilistic approach to patient dynamics. They
model the patient’s breathing motion using a probability mass function (PMF)
over a finite set of states. They then solve a robust optimization problem
that enforces dose bounds over a set of PMFs, which represent uncertainty
during treatment. This uncertainty set is updated after the dose fraction is
delivered, and the problem is re-solved for the next session. A follow-up
paper [MC15] numerically studies the effects of adjusting the target dose
based on the dose delivered to date. It shows that this type of adjustment can
lead to a high degree of heterogeneity in the per-fraction dose distribution,
which is undesirable from a medical standpoint. The paper recommends an
alternative method of uncertainty set adaptation to mitigate these effects.

The above analyses provide insight into the tradeoffs between different
fractionation schemes in a simple setting. However, most clinical cases are
more complex, involving multiple tumors, OARs, and nonlinear constraints. For
instance, dose-volume (i.e., percentile) constraints are widely used to limit
the radiation exposure of a percentage of an anatomical structure, such as the
spine. These constraints are nonconvex, but can be approximated by a convex
restriction [HCB06, ZST+13, FUXB19]. In [SGK16], the authors consider a
dynamic setting with multiple OARs and dose-volume constraints. Starting from
a given set of beam intensities, they solve for the optimal number of sessions
and OAR sparing factors. They also derive sufficient conditions under which
the optimal treatment consists of equal dose fractions. In a follow-up paper
[SGK17], the authors integrate the spatial and temporal aspects of the
problem, treating both beam intensities and number of sessions as variables.
Restricting their attention to equal fractions, they propose a two-stage
solution algorithm: in the first stage, they solve for the optimal beams given
each potential fixed number of sessions, and in the second stage, they select
the number of sessions based on the optimal objectives from the first stage.
They show that their method achieves better tumor ablation than conventional
IMRT or the spatiotemporally separated method.

Perhaps the paper most similar to ours is [KGP12]. In it, the authors propose
a stochastic control formulation of the adaptive treatment planning problem
with multiple tumors and OARs. They estimate the radiation response of the
tumors with a log-linear cell kill model and the response of the OARs with the
standard LQ model. Their goal is to minimize the expected number of tumor
cells at the end of treatment subject to bounds on the radiobiological impact
on the OARs. Uncertainty arises in the cell model parameters, which may
fluctuate randomly between sessions, representing unpredictable changes in the
patient’s health status. The authors fix the number of sessions and focus on
optimizing with respect to the beam intensities. They show that their problem
is convex, so can be solved using a combination of standard stochastic control
methods and off-the-shelf convex solvers, and provide several examples
demonstrating the effectiveness of their approach.

###  1.2 Contribution

In this paper, we integrate the stochastic control approach with a distributed
optimization algorithm to produce a method for efficient large-scale adaptive
treatment planning. As clinical cases are quite complex, with tens of
thousands of beams and treatment that takes place over months, such methods
are necessary to construct plans in a timely fashion. (See [JZJ14] for a
review of previous work on high-performance computing in radiation therapy,
particularly treatment optimization). We formulate the adaptive treatment
planning problem as a finite-horizon nonconvex optimal control problem. To
solve it, we introduce an operator splitting algorithm, which is based on
solving a sequence of convex approximations. Our algorithm is naturally
parallelizable and can handle a large number of beams, sessions, and
anatomical targets or OARs. Moreover, it can be combined with model predictive
control to produce treatment plans that are robust to errors and uncertainty
about the patient’s health status. We illustrate our algorithm’s performance
on a synthetic case, as well as a large prostate cancer case, and provide an
implementation in the Python package AdaRad: https://github.com/anqif/adarad.

##  2 Problem formulation

In radiation treatment, beams of ionizing radiation are delivered to a patient
from an external source. The goal is to damage or kill diseased tissue, while
minimizing harm to surrounding healthy organs. A course of treatment is
generally divided into $T$ sessions. At the start of session $t$, the
clinician chooses the intensity levels of the $n$ beams, denoted by
$b_{t}\in{\mbox{\bf R}}_{+}^{n}$. Typically, $T\approx 20$ and $n$ is on the
order of $10^{3}$ to $10^{4}$. We are interested in determining the best
sequence of beam intensities $b=(b_{1},\ldots,b_{T})$, otherwise known as a
treatment plan, subject to upper bounds $B_{t}\in\bar{\mbox{\bf R}}_{+}^{n}$
on $b_{t}$ for $t=1,\ldots,T$.

#### Anatomy and doses.

The beams irradiate an area containing $K$ anatomical structures, labeled
$i\in\\{1,\ldots,K\\}$, where usually $K<10$. A subset
$\mathcal{T}\subset\\{1,\ldots,K\\}$ are targets/tumors and the rest are OARs.
The dose delivered to each structure is linear in the beam intensities. We
write the dose vector $d_{t}=A_{t}b_{t}$ with $A_{t}\in{\mbox{\bf
R}}_{+}^{K\times n}$ a known matrix that characterizes the physical effects
and define $d=(d_{1},\ldots,d_{T})$. Notice that since $b_{t}$ and $A_{t}$ are
nonnegative, $d_{t}\geq 0$.

In every session, we impose a penalty on $d_{t}$ via a dose penalty function
$\phi_{t}:{\mbox{\bf R}}^{K}\rightarrow{\mbox{\bf R}}\cup\\{\infty\\}$. A
common choice is

$\phi_{t}(d_{t})=\theta_{t}^{T}d_{t}+\xi_{t}^{T}d_{t}^{\circ 2},$

where $\theta_{t}\in{\mbox{\bf R}}^{K}$ and $\xi_{t}\in{\mbox{\bf R}}_{+}^{K}$
are constants. Here $d_{t}^{\circ 2}=d_{t}\odot d_{t}$ denotes the vector
$d_{t}$ with each element squared. The total dose penalty over all sessions is

$\phi(d)=\sum_{t=1}^{T}\phi_{t}(d_{t}).$

Additionally, we enforce upper bound constraints $d_{t}\leq D_{t}$, where
$D_{t}\in\bar{\mbox{\bf R}}_{+}^{K}$ is the maximum dose in session $t$.

#### Health dynamics.

To assess treatment progress, we examine the health status of each anatomical
structure and encode these statuses in a vector $h_{t}\in{\mbox{\bf R}}^{K}$.
For now, the details of this encoding do not matter. Typically, $h_{ti}$
represents an estimate of the total surviving cells in structure $i$. Hence if
$i\in\mathcal{T}$, a smaller $h_{t,i}$ is desirable (since the tumor is
shrinking), while if $i\notin\mathcal{T}$, a larger $h_{t,i}$ is desirable.

From an initial $h_{0}$, the health status evolves in response to the
radiation dose and various other biophysical factors that depend on the
patient’s anatomy, generating a health trajectory $h=(h_{1},\ldots,h_{T})$.
Here we represent its dynamics as

$h_{t}=f_{t}(h_{t-1},d_{t}),\quad t=1,\ldots,T,$ (1)

where $f_{t}:{\mbox{\bf R}}^{K}\times{\mbox{\bf R}}^{K}\rightarrow{\mbox{\bf
R}}^{K}$ is a known mapping function. In this paper, we focus on the linear-
quadratic (LQ) model in which

$f_{t,i}(h_{t-1},d_{t})=h_{t-1,i}-\alpha_{t,i}d_{t,i}-\beta_{t,i}d_{t,i}^{2}+\gamma_{t,i},\quad
i=1,\ldots,K,\quad t=1,\ldots,T$ (2)

with constants $\alpha_{t}\in{\mbox{\bf R}}^{K},\beta_{t}\in{\mbox{\bf
R}}_{+}^{K}$, and $\gamma_{t}\in{\mbox{\bf R}}^{K}$. This model is commonly
used to approximate cellular response to radiation [Fow89, TH87, Bre08].
Specifically, in the LQ + time framework [TT87], $h_{t,i}$ is the log of the
fraction of surviving cells in structure $i$ after a dose $d_{t,i}$, while
$\alpha_{t,i}/\beta_{t,i}$ and $\gamma_{t,i}$ are constants related to the
structure’s survival curve and repair/repopulation rate, respectively. Notice
that equation (2) implies that the health status of each structure evolves
independently of the others.

#### Health penalty and constraints.

In order to control the patient’s health, we introduce a health penalty
function $\psi_{t}:{\mbox{\bf R}}^{K}\rightarrow{\mbox{\bf
R}}\cup\\{\infty\\}$ that imposes a penalty on $h_{t}$. Moreover, we assume
that

$\psi_{t}(h_{t})=\psi_{t}(h_{t,1},\ldots,h_{t,K})\mbox{ is monotonically
}\begin{cases}\mbox{increasing in }h_{t,i}&i\in\mathcal{T}\\\ \mbox{decreasing
in }h_{t,i}&i\notin\mathcal{T}\end{cases}$ (3)

for $t=1,\ldots,T$. This means that for a target, the health penalty increases
as the health status increases, while for an organ-at-risk, the health penalty
decreases as the health status increases. The assumption is reasonable if, for
instance, $h_{t}$ is a measure of cell survival in session $t$, so a lower
(higher) status is desirable for a target (organ-at-risk). An example of a
penalty function that satisfies (3) is

$\psi_{t}(h_{t})=\overline{w}^{T}(h_{t}-h_{t}^{\text{goal}})_{+}+\underline{w}^{T}(h_{t}-h_{t}^{\text{goal}})_{-},$

where $h_{t}^{\text{goal}}\in{\mbox{\bf R}}^{K}$ is the desired health status
and $\underline{w}\in{\mbox{\bf R}}_{+}^{K}$ and $\bar{w}\in{\mbox{\bf
R}}_{+}^{K}$ are parameters with $\underline{w}_{i}=0$ for $i\in\mathcal{T}$
and $\overline{w}_{i}=0$ for $i\notin\mathcal{T}$. Here we define
$(x)_{+}=\max(x,0)$ and $(x)_{-}=-\min(x,0)$ applied elementwise to $x$. The
total health penalty is

$\psi(h)=\sum_{t=1}^{T}\psi_{t}(h_{t}).$

In addition, we enforce bounds $H_{t}\in\bar{\mbox{\bf R}}^{K}$ on the health
status such that $h_{t,i}\leq H_{t,i}$ for $i\in\mathcal{T}$ and $h_{t,i}\geq
H_{t,i}$ for $i\notin\mathcal{T}$.

#### Optimal control problem.

Given an initial health status $h_{0}$, we wish to select a treatment plan
that minimizes the total penalty across all sessions. Thus, our problem is

$\begin{array}[]{lll}\mbox{minimize}&\sum_{t=1}^{T}\phi_{t}(d_{t})+\sum_{t=1}^{T}\psi_{t}(h_{t})&\\\
\mbox{subject to}&h_{t}=f_{t}(h_{t-1},d_{t}),&t=1,\ldots,T,\\\ &h_{t,i}\leq
H_{t,i},~{}i\in\mathcal{T},\quad h_{t,i}\geq
H_{t,i},~{}i\notin\mathcal{T},&t=1,\ldots,T,\\\ &d_{t}=A_{t}b_{t},\quad 0\leq
d_{t}\leq D_{t},\quad 0\leq b_{t}\leq B_{t},&t=1,\ldots,T\end{array}$ (4)

with variables $(b_{1},\ldots,b_{T}),(d_{1},\ldots,d_{T})$, and
$(h_{1},\ldots,h_{T})$. This is a discrete-time optimal control problem. If
$\phi_{t}$ and $\psi_{t}$ are convex and $f_{t}$ is affine, e.g., $f_{t}$ is
given by (2) with quadratic dose effect $\beta_{t}=0$, it is also convex and
can be solved directly using standard convex solvers.

##  3 Lossless relaxation

For the remainder of this paper, we restrict our attention to a convex
objective function and linear-quadratic health dynamics (2). In this case,
condition (3) allows us to relax the health dynamics constraint so problem (4)
can be written equivalently as

$\begin{array}[]{lll}\mbox{minimize}&\sum_{t=1}^{T}\phi_{t}(d_{t})+\sum_{t=1}^{T}\psi_{t}(h_{t})&\\\
\mbox{subject to}&h_{t,i}\geq
f_{t,i}(h_{t-1},d_{t}),~{}i\in\mathcal{T},&t=1,\ldots,T,\\\ &h_{t,i}\leq
f_{t,i}(h_{t-1},d_{t}),~{}i\notin\mathcal{T},&t=1,\ldots,T,\\\ &h_{t,i}\leq
H_{t,i},~{}i\in\mathcal{T},\quad h_{t,i}\geq
H_{t,i},~{}i\notin\mathcal{T},&t=1,\ldots,T,\\\ &d_{t}=A_{t}b_{t},\quad 0\leq
d_{t}\leq D_{t},\quad 0\leq b_{t}\leq B_{t},&t=1,\ldots,T.\end{array}$ (5)

The equality constraint $h_{t}=f_{t}(h_{t-1},d_{t})$ has been replaced with
two inequality constraints: a lower bound for targets and an upper bound for
OARs. Notice that the first inequality is the only nonconvex constraint in
(5). Our relaxed problem has the same solution set as (4) because these two
inequalities are tight at the optimum.

######  Proposition 1.

Let $(b^{\star},d^{\star},h^{\star})$ be a solution to problem (5). If
conditions (2) and (3) hold,

$h_{t}^{\star}=f_{t}(h_{t-1}^{\star},d_{t}^{\star}),\quad t=1,\ldots,T.$

###### Proof.

Suppose there exist some $t\in\\{1,\ldots,T\\}$ and $i\in\mathcal{T}$ such
that $h_{t,i}^{\star}>f_{t,i}(h_{t-1}^{\star},d_{t}^{\star})$. Then, we can
choose an $\epsilon>0$ such that
$h_{t,i}^{\star}>h_{t,i}^{\star}-\epsilon>f_{t,i}(h_{t-1}^{\star},d_{t}^{\star})$.
Since $f_{s,i}(h_{s-1},d_{s})$ is nondecreasing in $h_{s-1,i}$ for all
$s\in\\{1,\ldots,T\\}$, the point $(b^{\star},d^{\star},\hat{h})$ with

$\hat{h}_{s,j}=\begin{cases}h_{s,j}^{\star}-\epsilon&s=t,~{}j=i\\\
h_{s,j}^{\star}&\mbox{otherwise}\end{cases}$

is feasible for problem (5) because

$\hat{h}_{t,i}>f_{t,i}(h_{t-1}^{\star},d_{t}^{\star})\geq
f_{t,i}(\hat{h}_{t-1},d_{t}^{\star}),\quad h_{t+1,i}^{\star}\geq
f_{t+1,i}(h_{t}^{\star},d_{t}^{\star})\geq
f_{t+1,i}(\hat{h}_{t},d_{t}^{\star}),$

and $\hat{h}_{t,i}<h_{t,i}^{\star}\leq H_{t,i}$. Moreover, by condition (3),
$\psi_{t}(\hat{h}_{t})<\psi_{t}(h_{t}^{\star})$ so
$(b^{\star},d^{\star},\hat{h})$ achieves a lower objective value than
$(b^{\star},d^{\star},h^{\star})$, contradicting our original assumption. An
analogous argument holds for $t\in\\{1,\ldots,T\\}$ and $i\notin\mathcal{T}$
such that $h_{t,i}^{\star}<f_{t,i}(h_{t-1}^{\star},d_{t}^{\star})$ with
$\hat{h}_{t,i}=h_{t,i}^{\star}+\epsilon$. ∎

##  4 Sequential convex optimization

###  4.1 Algorithm description

Problem (5) is in general nonconvex because the target’s health dynamics
constraint

$h_{t,i}\geq f_{t,i}(h_{t-1},d_{t}),~{}i\in\mathcal{T},\quad t=1,\ldots,T$ (6)

is nonconvex when any $\beta_{t}\neq 0$. However, we can derive an estimate of
its optimum by solving a sequence of convex approximations. Each approximation
is formed by linearizing the health dynamics function (2) around a fixed dose
point and replacing the right-hand side of (6) with this linearization minus a
slack variable. The slack allows for a degree of error in the approximation
and is penalized in the objective.

More precisely, let $d_{t}^{s}\in{\mbox{\bf R}}^{K}$ for $t=1,\ldots,T$.
Define the linearized dynamics function

$\hat{f}_{t,i}(h_{t-1},d_{t};d_{t}^{s})=h_{t-1,i}-\alpha_{t,i}d_{t,i}-\beta_{t,i}d_{t,i}^{s}(2d_{t,i}-d_{t,i}^{s})+\gamma_{t,i},\quad
i=1,\ldots,K.$ (7)

This function is an upper bound on the LQ function (2) because $\beta_{t}\geq
0$. We replace the nonconvex constraint (6) in problem (5) with the affine
constraint

$h_{t,i}=\hat{f}_{t,i}(h_{t-1},d_{t};d_{t}^{s})-\delta_{t,i},~{}i\in\mathcal{T},\quad
t=1,\ldots,T,$ (8)

where $\delta_{t}\in{\mbox{\bf R}}_{+}^{K}$ is a slack variable. (The
inequality can been tightened into an equality due to Proposition 1). Convex
approximation $s$ is then

$\begin{array}[]{lll}\mbox{minimize}&\sum_{t=1}^{T}\phi_{t}(d_{t})+\sum_{t=1}^{T}\psi_{t}(h_{t})+\lambda\sum_{t=1}^{T}\mathbf{1}^{T}\delta_{t}&\\\
\mbox{subject
to}&h_{t,i}=\hat{f}_{t,i}(h_{t-1},d_{t};d_{t}^{s})-\delta_{t,i},~{}i\in\mathcal{T},\quad\delta_{t}\geq
0&t=1,\ldots,T,\\\ &h_{t,i}\leq
f_{t,i}(h_{t-1},d_{t}),~{}i\notin\mathcal{T},&t=1,\ldots,T,\\\ &h_{t,i}\leq
H_{t,i},~{}i\in\mathcal{T},\quad h_{t,i}\geq
H_{t,i},~{}i\notin\mathcal{T},&t=1,\ldots,T,\\\ &d_{t}=A_{t}b_{t},\quad 0\leq
d_{t}\leq D_{t},\quad 0\leq b_{t}\leq B_{t},&t=1,\ldots,T\end{array}$ (9)

with variables
$(b_{1},\ldots,b_{T}),(d_{1},\ldots,d_{T}),(h_{1},\ldots,h_{T})$, and
$(\delta_{1},\ldots,\delta_{T})$ and slack penalty parameter $\lambda>0$. The
parameter $\lambda$ is typically determined empirically or by trial-and-error.
This problem is convex and can be solved using standard convex solvers. Given
a solution to (9), we set the next linearization point
$d^{s+1}=(d_{1}^{s+1},\ldots,d_{T}^{s+1})$ equal to the optimal dose.

  * Algorithm 4.1  _Sequential convex optimization._

input: initial point $d^{0}$, parameter $\lambda>0$.
---
for $s=0,1,\ldots$ do
|  1\.  _Linearize._ For $t=1,\ldots,T$, form the linearization (7) around
$d_{t}^{s}$.
|  2\.  _Solve._ Set $d^{s+1}$ equal to an optimal dose of problem (9).
until stopping criterion (10) is satisfied.

Algorithm 4.1 is a special case of the convex-concave procedure (CCP) [YR03,
LB16, SDGB16], which is itself a form of majorization-minimization [HL04,
SBP17]. CCP is a heuristic for finding a local optimum of a nonconvex
optimization problem. It is guaranteed to converge; indeed, when certain
differentiability conditions are met, it converges to a stationary point
[SL09]. As a descent algorithm, CCP is usually terminated when the change in
the objective falls below some user-specified threshold $\epsilon>0$, i.e.,

$p_{\text{opt}}^{s}-p_{\text{opt}}^{s+1}<\epsilon,$ (10)

where $p_{\text{opt}}^{s}$ is the optimal objective of problem (9). In our
simple experiments, we have found that an initial linearization point of
$d^{0}=0$ and threshold of $\epsilon=10^{-3}$ produce good results.

###  4.2 Illustrative example

#### Problem instance.

We consider an example with $n=1000$ beams divided into $50$ bundles of $20$
parallel beams each, positioned evenly around a half-circle. There are $K=5$
structures, a single target $\mathcal{T}=\\{1\\}$ and four OARs (including
generic body voxels) depicted in Figure 1. Treatment takes place over $T=20$
sessions, so the basic problem has $nT+2KT=20200$ variables.

Figure 1: Anatomical structures for Example 4.2. Red is the target ($i=1$),
while green ($i=2$), blue ($i=3$), and orange ($i=4$) are specific OARs. White
denotes the non-target body voxels ($i=5$).

The patient’s initial health status is $h_{0}=(1,0,0,0,0)$. His status evolves
according to equation (2) with

$\displaystyle\alpha_{t}$ $\displaystyle=(0.01,0.50,0.25,0.15,0.005),$
$\displaystyle\beta_{t}$ $\displaystyle=(0.001,0.05,0.025,0.015,0.0005),$
$\displaystyle\gamma_{t}$ $\displaystyle=(0.05,0,0,0,0)$

over all sessions $t=1,\ldots,T$.

We set the health penalty function to

$\psi_{t}(h_{t})=(h_{t,1})_{+}+\sum_{i=2}^{5}(h_{t,i})_{-},\quad
t=1,\ldots,T.$

This function penalizes positive statuses of the target and negative statuses
of the OARs. Moreover, we constrain the target’s health status to be
$h_{t,1}\leq 2.0$ for $t=1,\ldots,15$ and $h_{t,1}\leq 0.05$ for the remaining
sessions, and we enforce a bound on the other structures’ health statuses of
$(h_{t,2},h_{t,3},h_{t,4},h_{t,5})\geq(-1.0,-2.0,-2.0,-3.0)$. Thus,

$H_{t}=\begin{cases}(2.0,-1.0,-2.0,-2.0,-3.0)&t=1,\ldots,15\\\
(0.05,-1.0,-2.0,-2.0,-3.0)&t=16,\ldots,T.\end{cases}$

For the dose penalty function, we choose

$\phi_{t}(d_{t})=\sum_{i=1}^{4}d_{t,i}^{2}+0.25d_{t,5}^{2},\quad
t=1,\ldots,T.$

In addition, we restrict the dose and beam intensity to be no more than
$D_{t}=20$ and $B_{t}=1.0$, respectively, over all sessions $t$.

#### Computational details.

We implemented Algorithm 4.1 in Python using CVXPY [DB16] and solved problem
(9) with MOSEK [AA00]. From an initial $d^{0}=0$ and $\lambda=10^{4}$, the
algorithm converged in $11$ iterations to a threshold of $\epsilon=10^{-3}$.
Total runtime was approximately $17$ seconds on a 64-bit Ubuntu OS desktop
with $8$ 4-core Intel i7-4790k / $4.00$ GHz CPUs and $16$ GB of RAM.

#### Results and analysis.

The optimal treatment plan is depicted in Figure 2. Beams are densely
clustered diagonal from the vertical, striking the target while largely
sparing the OARs. As the sessions continue, the number of beams slowly
increases, damaging some of the less sensitive organs ($i=3$ and $4$). Then at
$t=16$, when the target’s health bound becomes more stringent, the beam
density drops precipitously so that only a narrow bundle remains focused on
the target, keeping its health status at the desired level.

Figure 3 shows the radiation dose and health status resulting from this plan.
The latter was computed by plugging the optimal dose into equation (2). Total
dose to the target ($i=1$) and body voxels ($i=5$) far exceed the dose to any
other structures. By the end of treatment, the target’s health status has
fallen to a steady $0.05$, while the health statuses of the OARs remain within
their respective lower limits.

Figure 2: Optimal beam intensities for Example 4.2.

(a) Dose Trajectories

(b) Health Trajectories

Figure 3: Optimal (a) radiation dose and (b) health status trajectories for
Example 4.2. Dashed lines represent the bounds on the dose and health status
($D_{t,i}$ and $H_{t,i}$, respectively).

##  5 Model predictive control

###  5.1 Algorithm description

So far, we have assumed that at the time of planning, $f_{t}$ perfectly
captures the health dynamics from $t=1,\ldots,T$. This is rarely true in
practice. A patient’s anatomy changes unpredictably between sessions,
affecting the dispersion of radiation beams and the course of their health
status. We can incorporate these changes into problem (4) using model
predictive control (MPC).

MPC is a powerful technique for automatic control of complex, nonlinear,
stochastic systems. It performs extremely well even when the dynamics are
approximated by a simple model, since the system’s state is updated regularly
and new information is incorporated into the solution. This is particularly
fitting for radiation treatment planning.

As is customary in MPC, we first convert the state variable constraints in the
original problem into soft constraints, i.e., we remove the inequality
constraints on $h$ in (4) and add a penalty for violating them to the
objective. Let $c_{\tau}:{\mbox{\bf R}}^{K}\rightarrow{\mbox{\bf R}}$ be the
corresponding health violation penalty function, defined as

$c_{\tau}(h_{\tau})=\sum_{i\in\mathcal{T}}(h_{\tau,i}-H_{\tau,i})_{+}+\sum_{i\notin\mathcal{T}}(H_{\tau,i}-h_{\tau,i})_{+},\quad\tau=1,\ldots,T.$

This penalty function allows us to accommodate new and unexpected changes in
the patient’s health, such as the metastasis of a tumor that renders it
impossible to control without exceeding the health damage limit of an OAR.

We are now ready to describe MPC for our model. At the beginning of each
session $t$, we observe $A_{t},f_{t}$, and the patient’s true health status,
$h_{t-1}$, then form the problem

$\begin{array}[]{lll}\mbox{minimize}&\sum_{\tau=t}^{T}\phi_{\tau}(d_{\tau})+\sum_{\tau=t}^{T}\psi_{\tau}(h_{\tau})+\eta\sum_{\tau=t}^{T}c_{\tau}(h_{\tau})&\\\
\mbox{subject to}&h_{\tau}=f_{t}(h_{\tau-1},d_{\tau}),&\tau=t,\ldots,T,\\\
&d_{\tau}=A_{t}b_{\tau},\quad 0\leq d_{\tau}\leq D_{\tau},\quad 0\leq
b_{\tau}\leq B_{\tau},&\tau=t,\ldots,T\end{array}$ (11)

with variables $(b_{t},\ldots,b_{T}),(d_{t},\ldots,d_{T})$, and
$(h_{t},\ldots,h_{T})$ and violation penalty parameter $\eta>0$. Since
$c_{\tau}$ is convex, problem (11) is convex and can be solved using a slight
variation on Algorithm 4.1. Let $\bar{b}=(\bar{b}_{t},\ldots,\bar{b}_{T})$ be
the optimal treatment plan. We carry out only the first treatment,
$\bar{b}_{t}$, and update our observations $A_{t+1},f_{t+1}$, and $h_{t}$
based on the patient’s response. This process repeats until all $T$ sessions
have been completed.

###  5.2 Illustrative example

#### Problem instance.

We return to the setting of Example 4.2, except now, the health dynamics are
modeled with some error. Specifically, let $h_{t-1}$ be the patient’s health
status at the beginning of session $t$ and $d_{t}$ the dose delivered during
session $t$. Our model predicts the status will become
$\hat{h}_{t}=f_{t}(h_{t-1},d_{t})$. In fact, at the beginning of the next
session, we observe the true health status to be

$(h_{t})_{i}=\begin{cases}\max(\hat{h}_{t}+\omega_{t},0)_{i}&i\in\mathcal{T}\\\
\min(\hat{h}_{t}+\omega_{t},0)_{i}&i\notin\mathcal{T}\end{cases},$

where $\omega_{t}\in{\mbox{\bf R}}^{K}$ is drawn from $N(\mu,\sigma^{2}I)$.
This random process continues for $t=1,\ldots,T$.

For this example, we choose $\mu=0$ and $\sigma=0.1$. The rest of the
functions and parameter values are identical to 4.2. In particular, we still
employ the LQ model (2) with constant $\alpha_{t},\beta_{t}$, and $\gamma_{t}$
even though the health status is now stochastic. We plan the treatment using
MPC with $\eta=10^{4}$ and compare the results to those generated by the naive
approach, which simply solves problem (4) once prior to session 1.

#### Computational details.

We solved problem (11) using Algorithm 4.1 with $\lambda=10^{4}$ and
$\epsilon=10^{-3}$. For the initial dose in session 1, we chose $d^{0}=0$. In
each subsequent session $t$, we set $d^{0}$ to be the (truncated) optimal dose
point from the previous session, $(d_{t}^{\star},\ldots,d_{T}^{\star})$. With
these parameters, the algorithm took an average of 7 iterations per session to
achieve convergence; most runs completed in only 3–4 iterations. The total
runtime was 116 seconds.

#### Results and analysis.

Figure 4 depicts the treatment plan output by MPC. Most beams are aimed
slightly diagonal from the vertical, similar to the naive plan (Figure 2) up
to session 14. Then, the bundles of beams start to grow sparser and fan out,
hitting more areas of the OARs. This sparse irradiation pattern continues
until the final session, when there is a brief spike in intensity to bring the
target’s health status into the desired range.

In Figure 5(a), we plot the dose trajectories of the MPC plan (green) and the
naive plan (blue). The MPC curves are more jagged with a large spike at the
end of treatment. However, in each structure, the area under the MPC and naive
dose curves remains on par. Thus, we conclude that the MPC plan delivers about
the same amount of radiation as the naive plan, only spread across a wider
range of beam angles/intensities so as to compensate for uncertainty in the
health dynamics model.

This strategy results in better patient health as shown in Figure 5(b). The
MPC plan reduces the target’s health status to $0.05$, while maintaining the
health status of the OARs at a high level. Indeed, the health of these organs
under the MPC plan exceeds their health under the naive plan by a significant
margin in all but structure 4, where the two are relatively equal up until the
last session.

Figure 4: Optimal beam intensities for Example 5.2 using MPC.

(a) Dose Trajectories

(b) Health Trajectories

Figure 5: Optimal (a) radiation dose and (b) health status trajectories for
Example 5.2 using MPC (green) and a naive planning approach (blue). Dashed
lines represent the bounds on the dose and health status ($D_{t,i}$ and
$H_{t,i}$, respectively). The MPC plan’s health trajectories all remain within
the desired bounds, despite the error in the health dynamics model.

##  6 Operator splitting

MPC enables us to robustly handle uncertainty over time. However, another
challenge in radiation treatment planning is the sheer size of problems, which
makes them computationally difficult to solve in practice. A typical case with
$K=15$ and $n=10^{4}$ requires approximately $10^{5}$ floating-point
operations for the beam-to-dose calculation alone. Over a month of sessions,
that comes out to 4.5 million operations on a single machine.

In this section, we propose a fast, efficient method for solving the radiation
treatment planning problem using operator splitting. Our method is distributed
and scales readily with the number of beams as well as the length of
treatment. It can be applied both to the original problem (4) and the soft
constrained MPC variant (11). Below, we describe the mathematical details for
the former; the latter is a straightforward extension.

###  6.1 Consensus form

We first rewrite problem (4) in an equivalent consensus form:

$\begin{array}[]{lll}\mbox{minimize}&\sum_{t=1}^{T}\phi_{t}(d_{t})+\sum_{t=1}^{T}\psi_{t}(h_{t})&\\\
\mbox{subject to}&h_{t}=f_{t}(h_{t-1},\tilde{d}_{t}),\quad
0\leq\tilde{d}_{t}\leq D_{t},&t=1,\ldots,T,\\\ &h_{t,i}\leq
H_{t,i},~{}i\in\mathcal{T},\quad h_{t,i}\geq
H_{t,i},~{}i\notin\mathcal{T},&t=1,\ldots,T,\\\ &d_{t}=A_{t}b_{t},\quad 0\leq
d_{t}\leq D_{t},\quad 0\leq b_{t}\leq B_{t},&t=1,\ldots,T,\\\
&d_{t}=\tilde{d}_{t},&t=1,\ldots,T\end{array}$ (12)

with additional variable $\tilde{d}=(\tilde{d}_{1},\ldots,\tilde{d}_{T})$.
This splits the problem into two parts, one that encapsulates the radiation
physics and the other that contains the health dynamics. The parts share no
variables. They are only linked by the consensus constraint,
$d_{t}=\tilde{d}_{t}$, which requires their doses be equal.

###  6.2 ADMM

We solve problem (12) using an iterative algorithm called the alternating
direction method of multipliers (ADMM) [BPC+10]. In ADMM, the beams and health
statuses are optimized separately, taking into account the difference between
their resulting dose values. This difference is associated with a dual
variable $u=(u_{1},\ldots,u_{T})$, where each $u_{t}\in{\mbox{\bf R}}^{K}$,
which is updated every iteration in order to promote consensus.

  * Algorithm 6.1  _ADMM algorithm._

input: initial point $(\tilde{d}^{0},u^{0})$, parameter $\rho>0$.
---
for $k=0,1,\ldots$ do
|  1\.  _Calculate beams._ For $t=1,\ldots,T$, set the value of
$(b_{t}^{k+1},d_{t}^{k+1})$ to a
| solution of the problem
|
$\begin{array}[]{ll}\mbox{minimize}&\phi_{t}(d_{t})+\frac{\rho}{2}\|d_{t}-\tilde{d}_{t}^{k}-u_{t}^{k}\|_{2}^{2}\\\
\mbox{subject to}&d_{t}=A_{t}b_{t},\quad 0\leq d_{t}\leq D_{t},\quad 0\leq
b_{t}\leq B_{t}.\end{array}$
|  2\.  _Calculate health trajectory._ Set the value of
$(h^{k+1},\tilde{d}^{k+1})$ to a solution
| of the problem
|
$\begin{array}[]{lll}\mbox{minimize}&\sum_{t=1}^{T}\psi_{t}(h_{t})+\frac{\rho}{2}\|\tilde{d}-d^{k+1}+u^{k}\|_{2}^{2}&\\\
\mbox{subject to}&h_{t}=f_{t}(h_{t-1},\tilde{d}_{t}),\quad
0\leq\tilde{d}_{t}\leq D_{t},&t=1,\ldots,T,\\\ &h_{t,i}\leq
H_{t,i},~{}i\in\mathcal{T},\quad h_{t,i}\geq
H_{t,i},~{}i\notin\mathcal{T},&t=1,\ldots,T.\end{array}$
|  3\.  _Update dual variables._ $u^{k+1}:=u^{k}+\tilde{d}^{k+1}-d^{k+1}$.
until stopping criterion (17) is satisfied.

Here $1/\rho>0$ may be interpreted as the step size. Notice that the first
step of Algorithm 6.2 can be parallelized across sessions. We impose the dose
bound constraint on both the beam and health subproblems because it produces
faster convergence in practice.

#### Initialization.

For complex problems, the initial dose point $\tilde{d}^{0}$ can have a
significant impact on the performance of Algorithm 6.2. Below, we describe one
heuristic that produces a good starting point by solving a series of simple
optimization problems. The idea is to first find the optimal beams in the
static setting, where the entire treatment is delivered in a single session,
then rescale these static beams each session to form a plan in the dynamic
setting. By limiting ourselves to scaling factors, we significantly reduce the
size of the dynamic problem.

We begin by solving the static treatment planning problem

$\begin{array}[]{ll}\mbox{minimize}&\phi_{1}(d_{1})+\psi_{1}(h_{1})+\mu\mathbf{1}^{T}\zeta\\\
\mbox{subject to}&h_{1}=f_{1}(h_{0},d_{1}),\quad\zeta\geq 0,\\\ &h_{1,i}\leq
H_{T,i},~{}i\in\mathcal{T},\quad h_{1,i}\geq
H_{T,i}-\zeta_{i},~{}i\notin\mathcal{T},\\\ &d_{1}=A_{1}b_{1},\quad 0\leq
d_{1}\leq\sum_{t=1}^{T}D_{t},\quad 0\leq
b_{1}\leq\sum_{t=1}^{T}B_{t}\end{array}$ (13)

with respect to $b_{1}\in{\mbox{\bf R}}^{n},d_{1}\in{\mbox{\bf
R}}^{K},h_{1}\in{\mbox{\bf R}}^{K}$, and $\zeta\in{\mbox{\bf R}}^{K}$, where
$\mu>0$ is a slack penalty parameter. A reasonable choice for
$\mu=\frac{1}{K-|\mathcal{T}|}$, assuming there is at least one non-target
structure. Problem (13) is convex and can be easily handled on a single
machine (e.g., via interior-point methods) for up to $10^{5}$ beams. Denote
the optimal beam intensities by $b^{\text{stat}}$.

Next, we consider the dynamic treatment planning problem in which the beams
for each session are restricted to be a scalar multiple of $b^{\text{stat}}$,

$\begin{array}[]{lll}\mbox{minimize}&\sum_{t=1}^{T}\phi_{t}(d_{t})+\sum_{t=1}^{T}\psi_{t}(h_{t})+\mu\sum_{t=1}^{T}\mathbf{1}^{T}\zeta_{t}&\\\
\mbox{subject to}&h_{t}=f_{t}(h_{t-1},d_{t}),\quad\zeta_{t}\geq
0,&t=1,\ldots,T,\\\ &h_{t,i}\leq H_{t,i},~{}i\in\mathcal{T},\quad h_{t,i}\geq
H_{t,i}-\zeta_{t,i},~{}i\notin\mathcal{T},&t=1,\ldots,T,\\\
&d_{t}=\nu_{t}b^{\text{stat}},\quad 0\leq d_{t}\leq D_{t},\quad\nu_{t}\geq
0,&t=1,\ldots,T\end{array}$ (14)

with variables
$(\nu_{1},\ldots,\nu_{T}),(d_{1},\ldots,d_{T}),(h_{1},\ldots,h_{T})$, and
$(\zeta_{1},\ldots,\zeta_{T})$, where each $\nu_{t}\in{\mbox{\bf R}}$ and
$\zeta_{t}\in{\mbox{\bf R}}^{K}$. This problem can be solved using a slight
variation on Algorithm 4.1. (For the initial CCP point, we may use the optimal
time-invariant $\nu_{t}=\nu$ when $\beta_{t}=0$; finding this value entails
solving a small convex problem). Since there are only $O(TK)$ variables,
convergence is generally quick, taking less than 5 iterations in our
experiments. We use the resulting doses as our initial dose point for ADMM,
i.e., $\tilde{d}_{t}^{0}=\nu_{t}^{\star}b^{\text{stat}}$ for $t=1,\ldots,T$.

Besides providing a good starting point, this initialization heuristic also
gives us a way to quickly tune problem parameters. If the health trajectory
from $\tilde{d}^{0}$ is poor, it is much faster to modify weights and re-solve
problems (13) and (14) than it is to re-run the full ADMM algorithm.

#### Stopping criterion.

If problem (12) is convex, then under mild conditions, ADMM converges to a
solution assuming one exists. Moreover, the primal and dual residuals

$\displaystyle r_{\text{prim}}^{k}$ $\displaystyle=d^{k}-\tilde{d}^{k}$ (15)
$\displaystyle r_{\text{dual}}^{k}$
$\displaystyle=\rho(\tilde{d}^{k}-\tilde{d}^{k-1})$ (16)

also converge to zero. Thus, a reasonable stopping criterion is

$\|r_{\text{prim}}^{k}\|_{2}\leq\epsilon_{\text{prim}}\quad\mbox{and}\quad\|r_{\text{dual}}^{k}\|_{2}\leq\epsilon_{\text{dual}},$
(17)

where $\epsilon_{\text{prim}}>0$ and $\epsilon_{\text{dual}}>0$ are tolerances
for primal and dual feasibility, respectively. Typically, these tolerances are
chosen with respect to absolute and relative cutoffs $\epsilon_{\text{abs}}>0$
and $\epsilon_{\text{rel}}>0$ using the relation

$\displaystyle\epsilon_{\text{prim}}$
$\displaystyle=\epsilon_{\text{abs}}\sqrt{TK}+\epsilon_{\text{rel}}\max(\|d^{k}\|_{2},\|\tilde{d}^{k}\|_{2})$
$\displaystyle\epsilon_{\text{dual}}$
$\displaystyle=\epsilon_{\text{abs}}\sqrt{TK}+\epsilon_{\text{rel}}\|u^{k}\|_{2}.$

A common choice for $\epsilon_{\text{rel}}=10^{-3}$, while the choice for
$\epsilon_{\text{abs}}$ depends on the scale of the treatment planning problem
[BPC+10, §3.3.1].

#### Convergence and choice of $\rho$.

When the problem is convex, i.e., the health dynamics function is affine,
Algorithm 6.2 converges to a solution for any $\rho>0$, although the value of
$\rho$ may have an impact on the practical convergence rate. When the problem
is nonconvex, ADMM is a heuristic and the final beam/dose plan can depend
directly on $\rho$ [BPC+10, §9]. The question of how to choose $\rho$ is still
unsettled; see [GTSJ15, XFG17, XLLY17] for further discussion on the topic. We
have found that for data on the order of one, values of $\rho$ between
$10^{-2}$ and $10^{2}$ work reasonably well.

###  6.3 Clinical example

#### Problem instance.

We test our method on a fluence map optimization of a prostate cancer IMRT
case with $n=34848$ beams and $K=7$ structures consisting of a single PTV
($i=1$), five OARs, and generic body voxels ($i=7$). Treatment is carried out
over $T=45$ sessions, so the planning problem has about 1.6 million variables.
The matrix $A_{t}$ remains constant over time and maps the beam intensities to
the average dose per structure, i.e., $d_{t,i}$ is the total dose to structure
$i$ divided by the number of voxels in $i$. Each beam’s intensity cannot
exceed $B_{t}=0.025$.

The LQ model parameters, initial health status, and dose and health status
bounds can be found in Table 1; these have been adapted from prior clinical
datasets [Keh05, GMH+10, MYJ+10, vLOC+18]. We choose the health and dose
penalty functions to be

$\psi_{t}(h_{t})=(h_{t,1})_{+}+\frac{1}{6}\sum_{i=2}^{7}(h_{t,i})_{-},\quad\phi_{t}(d_{t})=\sum_{i=1}^{6}d_{t,i}^{2}+0.25d_{t,7}^{2},\quad
t=1,\ldots,T.$

These penalties place greater importance on reducing the health status of the
PTV compared to sparing the OARs or generic body tissue.

Table 1: Prostate IMRT Problem Parameters

LQ model Health and dose $i$ Structure $\alpha_{t,i}$ $\beta_{t,i}$
$\gamma_{t,i}$ $h_{0,i}$ $H_{t,i}$ $D_{t,i}$ 1 Prostate 0.15 0.05
$\begin{cases}0&t\leq 28\\\ 0.0173&t>28\end{cases}$ 5.8579
$\begin{cases}5.8579&t\leq 14\\\ 4.4716&15\leq t\leq 31\\\ 0&t>31\end{cases}$
10 2 Urethra 1 0.2 0 0 -4.8 10 3 Bladder 1 0.2 0 0 -4.8 10 4 Rectum 1 0.2 0 0
-4.8 10 5 L. Femoral Head 1 0.25 0 0 -3.0 10 6 R. Femoral Head 1 0.25 0 0 -3.0
10 7 Body 1 0.3333 0 0 -6.0 10

#### Computational details.

The computational setup is the same as in Example 4.2. To solve the ADMM
subproblems, we used MOSEK and ran CCP ($\lambda=10^{4}$) on the health
trajectory subproblem. With $\rho=80$, ADMM converged in $82$ iterations to
cutoffs of $\epsilon_{\text{abs}}=10^{-2}$ and
$\epsilon_{\text{rel}}=10^{-3}$. The normed residuals,
$\|r_{\text{prim}}^{k}\|_{2}$ and $\|r_{\text{dual}}^{k}\|_{2}$, are shown in
Figure 6. Total runtime was about 43 minutes, with the bulk of that time spent
on the main ADMM loop (initialization took only 32 seconds). By contrast, a
straightforward application of Algorithm 4.1 to this problem required over an
hour.

#### Results and analysis.

Figure 7 depicts the dose trajectories resulting from the initial plan (green)
and the final plan output by ADMM (blue). The initial plan is essentially a
piecewise equal-dose fractionation scheme, reflected by the flat plateaus in
the corresponding dose trajectories. This already gives us a good
approximation of the final plan: both plans maintain a relatively high dose to
the PTV of about 0.9 Gy until session 31, then drop off sharply to the same
constant doses thereafter. However, during the high dose phase, the final plan
gradually increases the dosage over time to all structures except the bladder
($i=3$). By adapting dynamically to changes in the patient’s anatomy, it is
able to deliver more dose per session and thus achieve better tumor control,
while still respecting the limits on the OARs’ health statuses.

Indeed, we see in Figure 8 that the final plan exactly attains the desired PTV
health status of zero for $t>31$. It must sacrifice some OARs to do this,
reducing the health statuses of the urethra, rectum, and right femoral head
($i=2,4$, and $6$) to their lower bounds, but never violates those bounds. In
fact, by shifting radiation to other structures, the final plan actually
improves the health of the bladder over that from the initial plan, which
results in a $h_{3}(t)$ far below the limit of $-4.8$ for $t\geq 35$. Overall,
it is clear that the combination of a solid initialization heuristic and ADMM
produces a treatment plan that satisfies or even exceeds all of our clinical
goals.

Figure 6: Primal and dual residual $\ell_{2}$-norms for Example 6.3. Figure 7:
Optimal radiation dose trajectory for Example 6.3. The initial plan (green)
depicts the dose output by the initialization heuristic described in §6.2,
while the final plan (blue) depicts the dose output by the ADMM algorithm.
Figure 8: Optimal health trajectories resulting from the doses in Figure 7.
Dashed lines represent the health status bounds $H_{t,i}$.

##  7 Implementation

We provide an implementation of our adaptive radiation treatment planning
method in AdaRad, an open-source Python software package based on CVXPY
[DB16]. Our implementation is fully distributed, leveraging Python’s built-in
multiprocessing library to execute solves in parallel. Users can quickly
import patient data, define clinical goals, construct treatment plans, and
visualize the results. They can also rapidly modify and re-plan a case,
allowing for comparisons between different prescriptions and treatment
lengths. Moreover, since AdaRad is a Python library, it can be easily
integrated with other libraries (e.g., for image processing) used in radiation
therapy.

The code below imports some patient data and a prescription, solves for the
optimal treatment plan, and plots the resulting dose and health trajectories.

    import adarad, numpy
    from adarad import Case, CasePlotter

    # Construct the clinical case.
    > case = Case()
    > case.import_file("/examples/patient_01-case.yaml")
    > case.physics.dose_matrix = numpy.load("/examples/patient_01-dmat.npy")

    # Solve using ADMM algorithm.
    > status, result = case.plan(slack_weight = 50, max_iter = 100,
                                 solver = ECOS, use_admm = True)
    > print("Solve status: {}".format(status))
    > print("Solve time: {}".format(result.solver_stats.solve_time))
    > print("Iterations: {}".format(result.solver_stats.num_iters))

    # Plot the dose and health trajectories.
    > caseviz = CasePlotter(case)
    > caseviz.plot_treatment(result, stepsize = 10)
    > caseviz.plot_health(result, stepsize = 10)

In this example, the dose matrix $A_{t}$ is the same for all $t$ and stored in
a single *.npy file. AdaRad also supports other sparse data representations,
such as scipy.csc_matrix. To specify a time-varying dose matrix, the user
would input a list of matrices in order $[A_{1},\ldots,A_{T}]$.

We start by constructing a Case, which contains Anatomy, Physics, and
Prescription objects. The Anatomy and Physics must be defined prior to
planning, either by manually specifying them in the code or importing a case
description. A description is a YAML file that contains at minimum the keys
treatment_length and structures, where the latter is a list of anatomical
structures $i=1,\ldots,K$, each of which has a name, is_target boolean
indicator, and alpha, beta, and gamma values corresponding to the LQ model
parameters. The initial health status and health and dose bounds may also be
specified.

Once the Case is defined, we can solve for the optimal treatment plan. The
plan function implements Algorithms 4.1 and 6.2 (the latter with use_admm =
True). It takes as optional input d_init: the initial dose point, use_slack: a
boolean indicating whether to include slack variable $\delta$, slack_weight:
the slack penalty parameter $\lambda$, max_iter: the maximum number of
iterations, and solver: the convex solver to use for the beam and health
subproblems. In the above example, we call the solver ECOS [DCB13], one of
several free, open-source solvers packaged with CVXPY. If MOSEK is installed,
we can call it as well by passing solver = MOSEK into the planning function.

After the algorithm finishes, plan saves the results in case.current_plan and
returns the final solve status along with a RunRecord object that carries
solver performance data, such as the total runtime, and the optimal variable
values. To visualize the resulting plan, we instantiate a CasePlotter object
and call plot_treatment and plot_health on the RunRecord to display the dose
and health trajectories, respectively. We can also extract the optimal beams,
doses, and health statuses with, e.g., result.beams for further processing.

If we wish to explore alternate plans, we can easily modify the dose and
health status constraints of any structure and re-plan the case. Re-planning
is generally fast, since AdaRad uses the previously stored solution as a warm
start point. In a typical workflow, we may import a prescription formed from
general clinical guidelines, then repeatedly adjust the dose/health status
bounds until we obtain a treatment plan with our desired properties. The
case.current_plan will be updated with the new optimal values after each run.
To keep a history of plans for comparison, we can save our results in the Case
by calling save_plan before re-optimizing. The code below provides an example
of changing the upper dose bound on the PTV to $D_{t,i}=10$ Gy for all
sessions and plotting the dose and health trajectories under this new
constraint alongside the trajectories of the original plan.

    # Save previous treatment plan.
    > case.save_plan("Original Plan")

    # Constraint allows maximum of 10 Gy per session on the PTV.
    > case.prescription["PTV"].dose_upper = 10

    # Re-plan the case with new dose constraint.
    > status2, result2 = case.plan(slack_weight = 50, max_iter = 100,
                                   solver = ECOS, use_admm = True)
    > print("Solve status: {}".format(status2))

    # Compare original and new treatment plans.
    > caseviz.plot_treatment(result2, stepsize = 10, label = "New Plan",
                             plot_saved = True)
    > caseviz.plot_health(result2, stepsize = 10, label = "New Plan",
                          plot_saved = True)

For more details on AdaRad’s functions as well as additional examples, see the
documentation at https://github.com/anqif/adarad.

##  8 Conclusion

To achieve the best outcomes, radiation therapy must adapt to new information
about the patient’s health and anatomy during treatment. We have described one
method for adaptive radiation treatment planning using an operator splitting
algorithm. Our method is highly scalable, parallelizable, and can efficiently
handle a large number of beams and sessions. Moreover, it is robust to errors
in the patient’s health response model, as well as other sources of
uncertainty in the clinic. We demonstrated its effectiveness on a large
prostate cancer case and showed that the resulting plan improves markedly on a
standard equal-dose fractionation scheme.

Future work will focus on expanding our health response model to include
sublethal damage repair, redistribution, and reoxygenation effects. We will
also incorporate dose-volume constraints into the optimal control problem.
Finally, to increase our algorithm’s speed, we intend to release an
implementation that takes advantage of the parallel processing capabilities of
the GPU.

## Acknowledgments

We thank Peng Dong for providing the anonymized dataset for the prostate
cancer IMRT case. This research was supported by the Stanford Graduate
Fellowship.

## References

  * [AA00] E. D. Andersen and K. D. Andersen.  The MOSEK interior point optimizer for linear programming: An implementation of the homogeneous algorithm.  In High Performance Optimization, pages 197–232. Springer-Verlag, 2000.
  * [AGRD10] D. M. Aleman, D. Glaser, H. E. Romeijn, and J. F. Dempsey.  Interior point algorithms: Guaranteed optimality for fluence map optimization in IMRT.  Physics in Medicine and Biology, 55(18):5467–5482, 2010.
  * [BHH+95] D. J. Brenner, L. R. Hlatky, P. J. Hahnfeldt, E. J. Hall, and R. K. Sachs.  A convenient extension of the linear-quadratic model to include redistribution and reoxygenation.  International Journal of Radiation Oncology, Biology, Physics, 32(2):379–390, 1995.
  * [BPC+10] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends in Machine Learning, 3(1):1–122, 2010.
  * [Bre08] D. J. Brenner.  The linear-quadratic model is an appropriate methodology for determining isoeffective doses at large doses per fraction.  Seminars in Radiation Oncology, 18(4):234–239, 2008.
  * [BRTU15] T. Bortfeld, J. Ramakrishnan, J. N. Tsitsiklis, and J. Unkelbach.  Optimization of radiation therapy fractionation schedules in the presence of tumor repopulation.  INFORMS Journal on Computing, 27(4):788–803, 2015.
  * [BV04] S. Boyd and L. Vandenberghe.  Convex Optimization.  Cambridge University Press, 2004.
  * [CM13] T. Chan and V. Mišić.  Adaptive and robust radiation therapy optimization for lung cancer.  European Journal of Operational Research, 231(3):745–756, 2013\.
  * [DB16] S. Diamond and S. Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  Journal of Machine Learning Research, 17(83):1–5, 2016.
  * [DCB13] A. Domahidi, E. Chu, and S. Boyd.  ECOS: An SOCP solver for embedded systems.  In Proceedings of the European Control Conference, pages 3071–3076, 2013.
  * [dlZAX07] A. de la Zerda, B. Armbruster, and L. Xing.  Formulating adaptive radiation therapy (ART) treatment planning into a closed-loop control framework.  Physics in Medicine and Biology, 52(14):4137–4153, 2007.
  * [Fow89] J. F. Fowler.  The linear-quadratic formula and progress in fractionated radiotherapy.  The British Journal of Radiology, 62(740):679–694, 1989.
  * [FUXB19] A. Fu, B. Ungun, L. Xing, and S. Boyd.  A convex optimization approach to radiation treatment planning with dose constraints.  Optimization and Engineering, 20(1):277–300, March 2019.
  * [FV04] M. Ferris and M. Voelker.  Fractionation in radiation treatment planning.  Mathematical Programming, 101(2):387–413, 2004.
  * [GAJ+12] T. Gupta, J. Agarwal, S. Jain, R. Phurailatpam, S. Kannan, S. Ghosh-Laskar, V. Murthy, A. Budrukkar, K. Dinshaw, K. Prabhash, P. Chaturvedi, and A. D’Cruz.  Three-dimensional conformal radiotherapy (3D-CRT) versus intensity modulated radiation therapy (IMRT) in squamous cell carcinoma of the head and neck: A randomized controlled trial.  Radiotherapy and Oncology, 104(3):343–348, 2012.
  * [Gao16] H. Gao.  Robust fluence map optimization via alternating direction method of multipliers with empirical parameter optimization.  Physics in Medicine and Biology, 61(7):2838–2850, 2016.
  * [GMH+10] M. Gao, N. A. Mayr, Z . Huang, H. Zhang, and J. Z. Wang.  When tumor repopulation starts? The onset time of prostate cancer during radiation therapy.  Acta Oncologica, 49(8):1269–1275, 2010.
  * [GTSJ15] E. Ghadimi, A. Teixeira, I. Shames, and M. Johansson.  Optimal parameter selection for the alternating direction method of multipliers (ADMM): Quadratic problems.  IEEE Transactions on Automatic Control, 60(3):644–658, 2015.
  * [HCB06] T. Halabi, D. Craft, and T. Bortfeld.  Dose-volume objectives in multi-criteria optimization.  Physics in Medicine and Biology, 51:3809–3818, 2006.
  * [HL04] D. R. Hunter and K. Lange.  A tutorial on MM algorithms.  The American Statistician, 58(1):30–37, 2004.
  * [JZJ14] X. Jia, P. Ziegenhein, and S. B. Jiang.  GPU-based high performance computing for radiation therapy.  Physics in Medicine and Biology, 51(4):R151–R182, 2014.
  * [Keh05] T. Kehwar.  Analytical approach to estimate normal tissue complication probability using best fit of normal tissue tolerance doses into the NTCP equation of the linear quadratic model.  Journal of Cancer Research and Therapeautics, 1(3):168–179, 2005\.
  * [KGP09] M. Kim, A. Ghate, and M. H. Phillips.  A Markov decision process approach to temporal modulation of dose fractions in radiation therapy planning.  Physics in Medicine and Biology, 54(14):4455–4476, 2009.
  * [KGP12] M. Kim, A. Ghate, and M. H. Phillips.  A stochastic control formalism for dynamic biologically conformal radiation therapy.  European Journal of Operational Research, 219(3):541–556, 2012\.
  * [LB16] T. Lipp and S. Boyd.  Variations and extension of the convex-concave procedure.  Optimization and Engineering, 17(2):263–287, June 2016.
  * [MC15] V. Mišić and T. Chan.  The perils of adapting to dose errors in radiation therapy.  PLoS ONE, 10(5):e0125335, 2015.
  * [MTD+12] M. Mizuta, S. Takao, H. Date, N. Kishimoto, K. L. Sutherland, R. Onimaru, and H. Shirato.  A mathematical study to select fractionation regimen based on physical dose distribution and the linear-quadratic model.  International Journal of Radiation Oncology, Biology, Physics, 84(3):829–833, 2012.
  * [MYJ+10] L. B. Marks, E. D. Yorke, A. Jackson, R. K. Ten Haken, L. S. Constine, A. Eisbruch, S. M. Bentzen, J. Nam, and J. O. Deasy.  Use of normal tissue complication probability models in the clinic.  International Journal of Radiation Oncology Biology Physics, 76(3):S10–S19, 2010.
  * [RAD+03] H. E. Romeijn, R. K. Ahuja, J. F. Dempsey, A. Kumar, and J. G. Li.  A novel linear programming approach to fluence map optimization for intensity modulated radiation therapy treatment planning.  Physics in Medicine and Biology, 48(21):3521–3542, 2003.
  * [RDL04] H. E. Romeijn, J. F. Dempsey, and J. G. Li.  A unifying framework for multi-criteria fluence map optimization models.  Physics in Medicine and Biology, 49(10):1991–2013, 2004.
  * [SBP17] Y. Sun, P. Babu, and D. P. Palomar.  Majorization-minimization algorithms in signal processing, communications, and machine learning.  IEEE Transactions on Signal Processing, 65(3):794–816, February 2017.
  * [SDGB16] X. Shen, S. Diamond, Y. Gu, and S. Boyd.  Disciplined convex-concave programming.  In Proceedings of the 55th IEEE Conference on Decision and Control, pages 1009–1014, 2016.
  * [SEP12] M. Y. Sir, M. A. Epelman, and S. M. Pollock.  Stochastic programming for off-line adaptive radiotherapy.  Annals of Operation Research, 196(1):767–797, 2012.
  * [SGK16] F. Saberian, A. Ghate, and M. Kim.  Optimal fractionation in radiotherapy with multiple normal tissues.  Mathematical Medicine and Biology, 33(2):211–252, 2016.
  * [SGK17] F. Saberian, A. Ghate, and M. Kim.  Spatiotemporally optimal fractionation in radiotherapy.  INFORMS Journal on Computing, 29(3):422–437, 2017.
  * [SL09] B. K. Sriperumbudur and G. R. Lanckriet.  On the convergence of the concave-convex procedure.  Advances in Neural Information Processing Systems, pages 1759–1767, 2009.
  * [TH87] H. D. Thames and J. H. Hendry.  Fractionation in Radiotherapy.  Taylor and Francis, 1987.
  * [TKL+10] H. D. Thames, D. Kuban, L. B. Levy, E. M. Horwitz, P. Kupelian, A. Martinez, J. Michalski, T. Pisansky, H. Sandler, W. Shipley, M. Zelefsky, and A. Zietman.  The role of overall treatment time in the outcome of radiotherapy of prostate cancer: An analysis of biochemical failure in 4839 men treated between 1987 and 1995.  Radiotherapy and Oncology, 96(1):6–12, 2010.
  * [TT87] E. L. Travis and S. L. Tucker.  Isoeffect models and fractionated radiation therapy.  International Journal of Radiation Oncology, Biology, Physics, 13(2):283–287, 1987.
  * [vLOC+18] C. M. van Leeuwen, A. L. Oei, J. Crezee, A. Bel, N. A. P. Franken, L. J. A. Stalpers, and H. P. Kok.  The alfa and beta of tumours: A review of parameters of the linear-quadratic model, derived from clinical radiotherapy studies.  Radiation Oncology, 13(96):1–11, 2018.
  * [WCP+06] S. L. Wolden, W. C. Chen, D. G. Pfister, D. H. Kraus, S. L. Berry, and M. J. Zelefsky.  Intensity-modulated radiation therapy (IMRT) for nasopharynx cancer: Update of the Memorial Sloan-Kettering experience.  International Journal of Radiation Oncology Biology Physics, 64(1):57–62, 2006.
  * [Web10] S. Webb.  Contemporary IMRT: Developing Physics and Clinical Implementation.  CRC Press, 2010.
  * [XFG17] Z. Xu, M. Figueiredo, and T. Goldstein.  Adaptive ADMM with spectral penalty parameter selection.  In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, pages 718–727, 2017.
  * [XLLY17] Y. Xu, M. Liu, Q. Lin, and T. Yang.  ADMM without a fixed penalty parameter: Faster convergence with new adaptive penalization.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages 1267––1277, 2017.
  * [YR03] A. L. Yuille and A. Rangarajan.  The concave-convex procedure (CCCP).  Neural Computation, 15(4):915–936, April 2003.
  * [YX05] Y. Yang and L. Xing.  Optimization of radiotherapy dose-time fractionation with consideration of tumor specific biology.  Medical Physics, 32(12):3666–3677, 2005.
  * [ZFH+00] M. J. Zelefsky, Zvi Fuks, L. Happersett, H. J. Lee, C. C. Ling, C. M. Burman, M. Hunt, T. Wolfe, E. S. Venkatraman, A. Jackson, M. Skwarchuk, and S. A. Leibel.  Clinical experience with intensity modulated radiation therapy (IMRT) in prostate cancer.  Radiotherapy and Oncology, 55(3):241–249, 2000.
  * [ZST+13] M. Zarepisheh, M. Shakourifar, G. Trigila, P. S. Ghomi, S. Couzens, A. Abebe, L. Noreña, W. Shang, S. B. Jiang, and Y. Zinchenko.  A moment-based approach for DVH-guided radiotherapy treatment plan optimization.  Physics in Medicine and Biology, 58(6):1869–1887, 2013.