High-dose-rate prostate brachytherapy inverse planning on dose-volume criteria by simulated annealing

TL;DR

This paper presents a simulated annealing algorithm for high-dose-rate prostate brachytherapy planning that improves speed and objective outcomes without relying on expensive solvers, using matrix operations and neighbor selection strategies.

Contribution

It introduces a novel simulated annealing approach for inverse planning in brachytherapy that outperforms existing methods in efficiency and solution quality.

Findings

Faster optimization compared to existing algorithms

Achieves better dose-volume objective values

Effective use of matrix multiplication and neighbor selection

Abstract

High-dose-rate brachytherapy is a tumor treatment method where a highly radioactive source is brought in close proximity to the tumor. In this paper we develop a simulated annealing algorithm to optimize the dwell times at preselected dwell positions to maximize tumor coverage under dose-volume constraints on the organs at risk. Compared to existing algorithms, our algorithm has advantages in terms of speed and objective value and does not require an expensive general purpose solver. Its success mainly depends on exploiting the efficiency of matrix multiplication and a careful selection of the neighboring states. In this paper we outline its details and make an in-depth comparison with existing methods using real patient data.

Tables

Table 1. Table 1: DVH protocol for a prescribed dose of 8.5 8.5 8.5 Gy for HDR prostate brachytherapy based on Hoskin et al. ( 2007 ) .

PTV	Rectum	Urethra
$V_{100\%}\ge 90\%$	$D_{10\%}\le 7.2\mathrm{Gy}$	$D_{10\%}\le 10\mathrm{Gy}$
	$D_{\max }\le 8\mathrm{Gy}$	$D_{\max }\le 10.6\mathrm{Gy}$

Table 2. Table 2: Radiation bounds and τ 𝜏 \tau for each structure.

	PTV	Rectum	Urethra	Shell
$L_i$	$8.5\mathrm{Gy}$	$7.2\mathrm{Gy}$	$10\mathrm{Gy}$	$8.5\mathrm{Gy}$
$U_i$	n/a	$8\mathrm{Gy}$	$10.6\mathrm{Gy}$	$8.5\mathrm{Gy}$
${\tau }_s$	n/a	$0.9$	$0.9$	$1$

Table 3. Table 3: The number of calculation points per structure.

Structure	Calculation Points
PTV	1750
Urethra	500
Rectum	250
Shell	250

Table 4. Table 4: The number of dwell positions and catheters for each patient.

Patient	Dwell Positions	Catheters
1	236	16
2	115	17
3	182	16

Equations

$${\left(\dot{D}t\right)}_i=\sum _{j\in J}{\dot{d}}_{ij}{t}_j,$$

$$\max$$

$${\left(\dot{D}t\right)}_i\ge L_i{v}_i$$

$${\left(\dot{D}t\right)}_i\le L_i+(U_i-L_i)(1-v_i)$$

$$\sum _{i\in I_s}{v}_i\ge {\tau }_s|I_s|$$

$$t_{j_1}\le \gamma t_{j_2}$$

$$v_i\in \left\{0,1\right\}$$

$$t_j\in \left[0,\infty \right)$$

$$\mathbb{P}\left(\text{accept tnew}\right)=\min \left\{1,\exp \frac{f\left(t^{new}\right)-f\left(t^{cur}\right)}{c_k}\right\}$$

$$c_k=c_0{\alpha }^{\lfloor \frac{k}{m}\rfloor },$$

$$t=t^{cur}+\lfloor r\rceil ,$$

$$|J_c|=\lceil \frac{V_{100\%}^k-V_{100\%}^{k-200,max}}{V_{100\%}^{k-200,max}}|J|\rceil ,$$

$$P\left(j\in J_c\right)=\frac{\sum _{i\in I_{PTV}:(\dot{D}t{)}_i<8.5}{\dot{d}}_{ij}}{\sum _{\tilde{j}\in J}\sum _{i\in I_{PTV}:(\dot{D}t{)}_i<8.5}{\dot{d}}_{i\tilde{j}}},$$

$$g=\lceil 40\left(1-\frac{V_{100\%}^k-V_{100\%}^{k-200,max}}{V_{100\%}^{k-200,max}}\right)\rceil .$$

$$\dot{D}{T}^{new}=\dot{D}{T}^{cur}+\dot{D}\left(T^{new}-T^{cur}\right).$$

$$y_i=\max \left\{0,\frac{8.5-{\left(\dot{D}t\right)}_i}{8.5}\right\}.$$

$$\underset{x\in [0,1],t\ge 0}{\max }\left\{\sum _{i\in I_{PTV}}{x}_i:(\dot{D}t{)}_i\ge 8.5x_i,(x,t)\in \mathcal{X}\right\},$$

$$\underset{y\in [0,1],t\ge 0}{\min }\left\{\sum _{i\in I_{PTV}}|y_i|:y_i\ge \frac{8.5-(\dot{D}t{)}_i}{8.5},(1-y,t)\in \mathcal{X}\right\}.$$

Full text

This is an author-created, un-copyedited version of an article published in Physics in Medicine and Biology DOI:10.1088/0031-9155/61/3/1155.

1 Introduction

High-dose-rate (HDR) brachytherapy is a form of radiation therapy in which the tumor is temporarily exposed to a highly radioactive source which dwells at different positions in or around the planning target volume (PTV). For prostate tumors, the dwell positions are offered by a temporary transperineal implant of catheters which run through the prostate. A remote afterloader, which is connected to all the catheters, then sends a radioactive source through the catheters one by one, stopping at several dwell positions inside the PTV. We assume that the catheter locations and the dwell positions are known. The goal of the treatment planner is to irradiate the PTV while sparing the surrounding organs at risk (OARs) by optimizing the dwell time at each dwell position. Traditionally, treatment planning was a forward process, where dwell times were adjusted until the dose distribution was satisfactory. In the last fifteen years, inverse planning has been developed as a computerized technique where the dwell times are optimized according to the treatment planner’s preferences on the dose distribution.

In the remainder of this introduction, we explain how dose distributions are evaluated, we provide a literature review on the ongoing progress of inverse planning, and we show how this work contributes to these developments.

DVH criteria. The dose distribution within the PTV or an OAR is described in a dose-volume histogram (DVH). Such a histogram displays the percentage of the structure receiving at least a certain dose. An example histogram is displayed in Figure 1. Points on this histogram are denoted by $D_{x}$ or $V_{y}$, e.g., $D_{10\%}$ is the minimum dose received in the hottest 10% of the structure under consideration and $V_{100\%}$ represents the percentage of the volume receiving 100% of the prescribed dose. In a slight deviation from the above notation, $D_{\max}$ is the maximal dose received in the structure. These statistics are currently the most important quantitative evaluation criteria (Hoskin et al., 2013). Besides these, the treatment planner inspects the isodose lines in order to avoid undesired hot spots and cold spots.

Table 1 displays clinically used DVH criteria from a local hospital. In order to achieve tumor control, $V_{100\%}$ of the PTV needs to be at least $90\%$, i.e., at least $90\%$ of the PTV’s volume need to receive at least $100\%$ of the prescribed dose. The limitation of complications is expressed in DVH constraints on the rectum and urethra. $D_{10\%}\leq 7.2\,\mathrm{Gy}$ for the rectum imposes that $90\%$ of the rectum’s tissue may not receive more than $7.2\,\mathrm{Gy}$. A hard upper bound for tissue in the rectum is set by $D_{\max}\leq 8\,\mathrm{Gy}$. For the urethra, $10\,\mathrm{Gy}$ is the upper bound for $90\%$ of the tissue. The remaining $10\%$ are allowed to reach radiation levels of at most $10.6\,\mathrm{Gy}$.

Dose measurement. In order to compute DVH statistics, the PTV and the OARs are discretized into a grid of calculation points. The dose rates $\dot{d}_{ij}$ describe the amount of radiation emitted from dwell position $j\in J$ towards calculation point $i\in I$ per second. The dose at a calculation point $i$ is equal to the sum of radiation received from each dwell position $j$:

[TABLE]

where $t_{j}$ is the dwell time at dwell position $j$. Let $\dot{D}$ be the $|I|\times|J|$ matrix with all dose rates $\dot{d}_{ij}$. The matrix-vector product $\dot{D}t$ yields a vector that denotes the dose per calculation point.

Literature review. Traditionally, inverse planning for HDR brachytherapy is dose-based: the treatment planner prescribes a lower bound $L_{i}$ and and an upper bound $U_{i}$ on the desired dose for each calculation point. Typically these bounds are the same for all calculation points within the same structure. At each calculation point $i$, the received dose $(\dot{D}t)_{i}$ is compared to its respective prescribed lower and upper bounds. If the dose lies outside the interval $[L_{i},U_{i}]$, a penalty is imposed that is linear or quadratic in the deviation. A treatment plan is optimal if it minimizes the sum of penalties over all calculation points. Linear penalty functions have been successfully implemented in the commercially available algorithms Inverse Planning by Simulated Annealing (IPSA) (Lessard and Pouliot, 2001) and Hybrid Inverse Treatment Planning and Optimization (HIPO) (Karabis et al., 2005). Quadratic penalty functions are used by Lahanas et al. (2003); Lahanas and Baltas (2003); Milickovic et al. (2002).

The disadvantage of dose-based penalty functions is that the resulting treatment plans may need a posteriori adjustments as they do not necessarily adhere to the DVH criteria. More recently, dose-volume based optimization methods have been developed that directly optimize on DVH criteria. Panchal (2008) has formulated a Harmony Search heuristic. Siauw et al. (2011) have presented a Mixed Integer Linear Programming (MILP) formulation and developed the heuristic Inverse Planning by Integer Program (IPIP) to solve it. IPIP determines a solution that is feasible and near-optimal for MILP by solving two Linear Programs (LPs). Gorissen et al. (2013) have devised a formulation in which they directly optimize an MILP using specific solver settings. Beliën et al. (2009) have introduced a hybrid approach based on simulated annealing and linear programming. Their objective is dose-based, while the delivered radiation to OARs is subject to dose-volume constraints. For the rectal constraint $D_{10\%}\leq 7.2\,\mathrm{Gy}$, the linear constraint $(\dot{D}t)_{i}\leq 7.2\,\mathrm{Gy}$ needs to hold for 90% of the calculation points. The authors utilize simulated annealing to determine which calculation points are in the group of 90%, and then use linear programming to determine the dwell times. According to Beliën et al. (2009), the algorithm yields results superior to running the integer program alone based on ten problem instances with $30$ minutes computation time each. The authors suggest the development of a pure SA algorithm without linear programming.

While our focus is on HDR brachytherapy, we briefly mention related research in Intensity Modulated Radiotherapy (IMRT). Instead of irradiating the PTV by bringing a radioactive source in close proximity to or inside the tumor, the patient is irradiated with radiation fields at different angles around the patient. For fixed beam angles, the optimization problem for determining the intensity map of each field is mathematically equivalent to the optimization problem we study. Historically, these research fields have not been connected, perhaps due to the different scale of the problems. For optimization in IMRT, we refer to Ehrgott et al. (2010) and the extensive literature review by Zaghian et al. (2014). Simulated annealing is also applied to IMRT treatment planning, e.g., Cho et al. (1998) optimize penalty functions to meet dose criteria.

Our contribution. This paper follows up on Beliën et al. (2009)’s suggestion to develop a pure simulated annealing (SA) heuristic to optimize on DVH statistics. We present DOPSA (DVH Optimization by Pure Simulated Annealing), a novel algorithm that solves the same MILP model as Gorissen et al. (2013) and Siauw et al. (2011). In this model the objective is to maximize $V_{100\%}$ for the PTV, i.e., the volume of the PTV receiving at least the prescribed dose, while meeting DVH constraints on the OARs. Our method exploits the speed of matrix multiplication and uses a specialized generator for neighboring solutions. The heuristic has been tuned by conducting many trial-and-error experiments. The parameters have been selected using a metamodel, and have been validated with an out-of-sample analysis.

The advantage of SA over the existing dose-volume based models lies in its simple implementation independent of costly LP or MILP solvers. Since optimization is often offered as a separate module for treatment planning systems, clinics may choose cheaper dose-based optimization modules or keep using them if the price difference is too large, leaving the potential of dose-volume based methods unused. Furthermore, the results in Section 3 show that our heuristic has clear advantages over existing methods in terms of speed and objective value.

A short optimization time is of practical importance. In the clinical workflow (Balvert et al., 2015), two treatment plans are made while the patient is anesthetized and continuously monitored by an anesthesia care team. Moreover, the radiation oncologist is on standby while the treatment planner makes the plan. Shortening the time needed to create a treatment plan improves the clinical workflow and reduces the treatment costs.

The implementation details of DOPSA will be explained in the upcoming section, followed by a discussion of the results and a comparison with existing algorithms in Section 3.

2 Methods

In this section we provide in-depth information about DOPSA. Its final implementation is the result of a long trial-and-error process where many ideas have been tested. Here, we restrict ourselves to a description of the final algorithm. For the selection procedure and parameter tuning based on trial-and-error and a metamodel we refer the reader to Deist (2013). There, values for the parameters (discussed below) have been determined by manual simulation and a linear metamodel has been fitted to fine-tune each parameter. However, the explanatory power of the model was low due to the presence of statistical noise. Note that only the data of Patient 1 has been employed in the parameter tuning.

First, we describe the optimization model (Section 2.1) and the main steps in the algorithm (Section 2.2). Then we focus on the two most important steps, which are the generation of neighborhood states (Section 2.3) and checking feasibility (Section 2.4).

2.1 Optimization model

We use the same model as Gorissen et al. (2013) and Siauw et al. (2011) to allow a comparison with their results. This model is based on the DVH protocol from Table 1 with exception for the constraint on $V_{100\%}$. Instead of formulating a constraint $V_{100\%}\geq 90\%$, the model is cast to maximize $V_{100\%}$ under $D_{10\%}$ and $D_{\max}$ constraints on the OARs. Additionally, to ensure that the dose distribution conforms to the PTV and to avoid high doses outside the PTV in regions where there is no OAR to limit the dose in that region, we also consider an artificial structure that surrounds the prostate at 2 mm distance. This collection of tissue is henceforth denoted as the ‘shell’, and a constraint is added to ensure that the maximum dose in this structure is below $8.5\,\mathrm{Gy}$. Furthermore, a dwell time modulation restriction (DTMR) is added to avoid excessively high dwell times. The DTMR prohibits the dwell times of two adjacent dwell positions within the same catheter to differ more than a factor two. A similar restriction is found in other algorithms (Baltas et al., 2009). The complete optimization model (Gorissen et al., 2013) is

[TABLE]

The binary variable $v_{i}$ indicates whether calculation point $i$ meets a certain dose criterion: for a calculation point in the PTV, i.e., $i\in I_{PTV}$, the corresponding $v_{i}$ is equal to one if the received dose at those points is at least the prescribed dose of $L_{i}=8.5\,\mathrm{Gy}$. These values are enforced by the objective (1a) and constraint (1b). For a calculation point in the OARs rectum (R), urethra (U), and shell (S), i.e., $i\in I_{R}\cup I_{U}\cup I_{S}$, constraint (1c) allows $v_{i}$ to be one if the received dose is at most the soft upper bound $L_{i}$. If the received dose is between $L_{i}$ and the hard upper bound $U_{i}$, $v_{i}=0$. Constraint (1d) enforces the soft upper bound on the dose for at least a fraction $\tau_{s}$ of the calculation points. The corresponding parameter values are displayed in Table 2. Note that for the shell, $L_{i}$ and $U_{i}$ coincide and $\tau_{s}=1$, since only a single dose constraint is imposed for this structure. Constraint (1e) enforces the DTMR, i.e., dwell times at neighboring dwell positions in the same catheter may not differ by more than a factor $\gamma$. For our implementation, $\gamma=2$. The set $\Gamma(j)$ contains all adjacent dwell positions of $j$.

2.2 Algorithm

The algorithm described below searches for a vector $t$ such that objective (1a) is maximal and all constraints (1b)-(1g) are satisfied. DOPSA has been implemented in MATLAB according to Algorithm 1, an adapted version of a standard simulated annealing process (Kirkpatrick et al., 1983). In simulated annealing, the feasible region is searched by selecting a new solution $t^{new}$ in the neighborhood of the current solution $t^{cur}$ and deciding whether to accept this new choice. A solution which is not decreasing the objective function is always accepted, whereas decreasing solutions are also accepted with a probability that decreases in the loss of objective function value and the duration of the search. This rule is also named the metropolis criterion. The acceptance probability is defined as

[TABLE]

with objective function $f$ and $c_{k}$ regulating the probability’s sensitivity to lower objective function values. The parameter $c_{k}$ is adjusted at each iteration $k$ according to the well-known exponential cooling schedule

[TABLE]

where $m\in(0,\infty)$ is a fixed number of iterations after which the initial temperature $c_{0}$ is rescaled by $\alpha\in(0,1)$.

DOPSA’s code is designed to generate and evaluate multiple states simultaneously. All time-consuming procedures within the heuristic thus become matrix-operations which decreases the required CPU time per state. Essentially it is faster to check the feasibility of $g$ states in one operation than to check the feasibility of the $g$ states individually, since computing $\dot{D}T$, where $T$ is a $|J|\times g$ matrix with the $g$ states as columns, takes less than $g$ times longer than computing $\dot{D}t$ for a single state $t$. This not only eliminates a for-loop in Matlab code, which is known to be relatively slow, but more importantly, uses the fact that an efficient implementation of a matrix multiplication procedure takes less operations than $\mathcal{O}(g|I||J|)$ when $g\geq 2$.

DOPSA is initialized with the dwell times set to zero because it should start in the feasible region and, a priori, any other feasible state can only be determined at relatively high computational effort. Then, one state is chosen from the set of feasible candidates in the neighborhood and is evaluated by the metropolis criterion as described in Section 2. Subsequently, the temperature is decreased in every iteration according to an exponential cooling schedule with parameters $\alpha=0.99$, $m=1$, and $c_{0}=1.5$. The search continues in the chosen state, if accepted, or restarts in the current state. After each 15,000 iterations, the state with the highest $V_{100\%}$ is selected and the search is restarted in the neighborhood of this state. This step is included to recover from an unsuccessful attempt to escape from a local optimum. Returning to the state with the highest $V_{100\%}$ focuses the search to areas with high objective function values and avoids unnecessary search efforts in low-potential regions.

The running time of DOPSA is restricted to a fixed limit (e.g., 180 seconds) to allow a comparison with existing methods and to test whether DOPSA can provide satisfactory results within a short time-span. In practice, an alternative stopping criterion can be used, such as the relative improvement of the objective value over the last $n$ iterations, or user intervention based on the coverage attained by the algorithm and the running time.

The performance of the algorithm was revealed to be insensitive to changes in the parameter values of $\alpha$ and $c_{0}$, suggesting that hill-climbing methods might be similarly successful.

2.3 Generating neighborhood states

DOPSA generates a neighborhood solution $t$ by perturbing the dwell times of the current state $t^{cur}$:

[TABLE]

where $r$ is a vector of the same size as $t$. The vector $r$ is chosen randomly according to $r_{j}\sim N(0,0.05)$ i.i.d. for all $j$ in $J_{c}$, the set of dwell positions that are changed, and $r_{j}=0$ for all $j$ in $J\backslash J_{c}$. The operator $\lfloor\rceil$ rounds the elements of $r$ to one decimal place, which is the input precision of the Flexitron remote afterloader (Nucletron BV, Veenendaal, the Netherlands). Negative dwell times are subsequently set to 0.

For the DTMR, the dwell times in each catheter are adjusted one by one to ensure the relative difference with the previous dwell time is less than a factor 2. This adjustment may increase the number of changed entries in the dwell time vector.

Number of changed entries. We observed that controlling the number of dwell time changes applied to a new state is fundamental to the performance of DOPSA. A lower number of changes is more likely to yield a feasible new state. However, lowering this number decreases the potential maximal increase in $V_{100\%}$ in the new state, which is most notable in early iterations. Since it is unlikely that big improvements can still be realized when the algorithm stalls, a negative correlation between the number of changed entries and the improvement over the past 200 iterations is introduced. Consequently, when the improvements get smaller, newly generated states get a larger probability of being feasible while the risk of missing out on potential large gains is minimal. In the first 200 iterations, all entries are subject to changes before the dynamic adjustment begins. The warmup period of 200 iterations only affects the initial phase of DOPSA since the total number of iterations is approximately 50000.

Let $J_{c}$ denote the index set of entries in $t^{cur}$ that are varied to construct a new vector $t$. The relation between the cardinality of $J_{c}$ and the improvement in $V_{100\%}$ is defined by

[TABLE]

where $V_{100\%}^{k}$ is the $V_{100\%}$ attained in the current state $t^{cur}$ at iteration $k$. $V_{100\%}^{k-200,max}$ is the overall highest $V_{100\%}$ attained until iteration $k-200$.

Note that the states chosen during the search can also exhibit $V_{100\%}$ lower than in preceding iterations due to the metropolis criterion, resulting in a negative $|J_{c}|$. For that reason and to maintain a natural upper bound, $|J_{c}|$ is restricted to the interval $[\lceil 0.02|J|\rceil,|J|]$. The lower bound is chosen to be at least 1 (since $|J|>50$ for practically all cases) since improvements in $V_{100\%}$ might only be attained with increases in dwell times for one dwell position while simultaneously decreasing the dwell time of another position. Taking $2\%$ of all entries as the lower bound ensures sufficient dwell time changes in generated states.

Empirical results indicate that a newly generated state has six times the probability of being feasible by using $J_{c}$ compared to perturbing all dwell times. The number of states that are dropped due to infeasibility decreases from 95% to 70%.

Dwell position selection. Once $|J_{c}|$, the number of dwell positions for which the dwell times are perturbed, is chosen, the next step is to select that many dwell positions. The intention is to alter dwell times such that the objective value is improved. This target is pursued by changing the dwell times of the dwell positions close to the calculation points in the PTV that do not yet receive the prescribed dose. Since the distance between dwell position $j$ and calculation point $i$ has an inverse squared relation to the dose rate $\dot{d}_{ij}$, the probability that dwell position $j$ is chosen in a single draw is therefore set to

[TABLE]

where $I_{PTV}$ is the set of calculation points in the PTV. The denominator normalizes the term to attain a probability distribution. Using probabilities rather than simply picking the closest dwell positions bears the advantage that the set of states that can be generated from the current state does not become too small. Occasionally also dwell positions further away are selected which diversifies the search.

This methods requires a random sample of dwell positions in each iteration of the search. Randomly sampling dwell positions without replacement in MATLAB is a computationally expensive procedure. Our experiments have shown that sampling $|J_{c}|$ items with replacement and removing duplicates yields better results, due to a larger number of iterations in the same running time. In a direct comparison, the average PTV coverage after a running time of 3 minutes over 50 replications was 0.12% higher when using sampling without replacement.

Number of generated states. In the early stages of the search process, improving states can be determined with few generated states per iteration since the improvement potential in the neighborhood of $t^{cur}$ is high. Later on, an intensification of the search in those neighborhoods is most important. Moreover, intensification becomes necessary only when the search approaches the boundary, where one expects to find optimal solutions. Therefore, DOPSA generates more states when the boundary is reached (as suggested by Hedar and Fukushima (2006)). Intensification can be controlled with an opposite relation between the number of states generated in a single iteration $g$ and $V_{100\%}$:

[TABLE]

$g$ is limited to the interval $[1,40]$ because trial runs indicated that excessive state generation does not further improve the DOPSA’s performance. In the first 200 iterations, $g=1$ is used.

2.4 Checking feasibility

The feasibility of each generated state needs to be verified, which is one of the computationally most expensive procedures next to computing $V_{100\%}$. Testing each constraint discussed in Section 1 requires a row subset of the matrix $\dot{D}$ corresponding to a specific OAR to be multiplied with the dwell time vector $t$. The computation time of checking a constraint increases with the number of calculation points in the corresponding OAR.

The sequence in which the constraints on the DVHs for the OARs are checked is chosen to minimize the required computation time: the most restrictive constraints with lowest computational cost are evaluated first. States that do not satisfy a constraint are immediately excluded in order to reduce the computation time for checking other constraints.

The calculation of constraint values and $V_{100\%}$ for newly generated states can be sped up by using the values for the current state $t^{cur}$ that have been computed in the preceding iteration. Let $T^{cur}$ denote the $|J|\times g$ matrix containing the current state as columns, and let $T^{new}$ denote the $|J|\times g$ matrix with the neighboring states as columns.

The dose computation for each generated state $\dot{D}T^{new}$ can be reformulated as follows:

[TABLE]

The matrix product $\dot{D}T^{cur}$ is known from the previous iteration, so only $\dot{D}(T^{new}-T^{cur})$ needs to be computed. $T^{new}-T^{cur}$ is a sparse matrix with at most $|J_{c}|$ non-zero entries per column (plus additional non-zero entries due to the DTMR adjustment), since states generated in the neighborhood of a current state only differ in approximately $|J_{c}|$-many entries. Moreover, the sparsity increases over the search process. Using MATLAB’s toolpack for multiplication of sparse matrices, an increase of the number of iterations over the total running time of 5% up to 21% could be observed, depending on the number of dwell positions.

3 Results

The performance of DOPSA is assessed by a direct comparison with two existing dose-volume-based methods: the MILP formulation by Gorissen et al. (2013) and the IPIP algorithm by Siauw et al. (2011). The running time has been fixed to 3 minutes on an Intel Xeon E5620 2.4 GHz with 6 GB of memory. Both MILP and IPIP are solved using the state-of-the-art solver CPLEX 12.4 (IBM).

3.1 Patient data

Tables 3 and 4 provide an overview of the three patient data sets. Each data set uses the same number of calculation points in the PTV and the OARs. The PTV contains the biggest share of calculation points since it is the largest organ. All dwell positions inside the PTV have been activated based on a step size of 2.5 mm, resulting in 115 to 236 dwell positions per patient. The number of inserted catheters is similar among patients with 17 for Patient 2 and 16 for the other two patients. The dose rates have been calculated using the TG-43 formalism.

The calculation points have been generated with a Sobol sequence, which is a quasi-random sampling technique that randomly distributes calculation points in such a way that the sampling density is approximately uniform. This avoids the possibility that small subregions are not covered with calculation points. Quasi-random sampling was shown to be excellent for evaluating treatment plans (Niemierko and Goitein, 1990), although this effect may disappear when the treatment plan is optimized with respect to the randomly sampled calculation points. To avoid excessively high doses in points close to dwell positions, distances smaller than 0.83 mm are set to 0.83 mm (in the same direction).

3.2 Running time and objective value

DOPSA has been tested on three patients. Since it is a random local search procedure, the result can be subject to random effects. Therefore, the heuristic has been run 250 times for each patient in order to provide a detailed analysis of the average treatment plan quality and the variation in quality. If a low variation can be ensured for this optimization algorithm, only a single optimization run per patient will be necessary in a clinical setting to obtain treatment plans with consistent quality.

Figures 2-4 show the results for all three patient data sets. The average and standard deviation in PTV coverage over 250 replications are displayed over the running time of 180 seconds. The positioning of catheters determines the maximally achievable level of $V_{100\%}$ in the PTV, which explains the difference in objective values between patients. The high dependence on catheter positioning is well observable for Patient 2 for whom only solutions with $V_{100\%}<50\%$ can be attained: inspecting the corresponding ultrasound scans reveals badly positioned catheters.

In all three cases, the PTV coverage increases rapidly within the first 30 seconds and further improvements on the treatment plan quality are achieved over the remaining running time. Every treatment plan determined by DOPSA is feasible, i.e., it always satisfies the DVH criteria imposed on the OARs and dwell time constraints (see constraints (1c)-(1g)).

The absolute standard deviations in PTV coverage after a running time of 3 minutes for each of the 250 replications are 0.09%, 0.33%, and 0.23%, respectively. Therefore, DOPSA delivers treatment plans with consistent quality across replications.

3.3 Comparison with existing DVH-based optimizers

For the comparison with existing DVH-based optimizers, we have run an MILP solver with the solver settings from Gorissen et al. (2013) and we have implemented IPIP (Siauw et al., 2011). The original IPIP has been extended by the DTMR to ensure a fair comparison. These two methods have been run on the same three patient data sets and on the same model. It is important to recognize that IPIP can only provide a single solution per data set, whereas MILP and DOPSA continue improving the solution. This is observable in Figures 5-7, where the line for IPIP remains flat after a single increase. In theory, MILP could stop when it proves that the solution is optimal, but this did not happen for any of the patients.

Figures 5 to 7 show the average PTV coverage obtained by DOPSA and the coverage obtained by the existing methods MILP and IPIP. The running time for Patient 3 has been extended to 5 minutes since MILP requires more than 3 minutes to find non-zero feasible points.

IPIP attains adequate PTV coverage within less than 6 seconds for each patient. All those solutions are initially superior to treatment plans found by DOPSA. Over the running time of 3 minutes, however, DOPSA consistently outperforms IPIP for all patient data sets. In each of the 250 replications for all three patients, the PTV coverage after 3 minutes is higher for DOPSA than for IPIP. The average PTV coverage attained by DOPSA exceeds the results by IPIP after only 18, 6, and 12 seconds, respectively.

MILP requires substantially more time to determine good treatment plans. For 2 out of 3 patients, MILP requires more than 1 minute to determine a treatment plan with non-zero PTV coverage. However, after 3 minutes for Patients 1-2 and 5 minutes for Patient 3 it has found the best plan among the three algorithms.

Unfortunately, the optimal objective values are not known. There is still a large gap between the best solution found by any of the algorithms and the upper bound provided by the MILP solver, and this gap does not decrease substantially after running the solver for a week. DOPSA differs greatly from the existing dose-volume based methods while it finds similar objective values. This could be an indication that much better solutions do not exist.

3.4 Model Discussion

The significance of the strong performance of DOPSA compared to the MILP by Gorissen et al. (2013) becomes apparent when closely investigating the structure of the MILP (1a)-(1g). Below, we show that the MILP is equivalent to the well-known $\ell_{0}$-norm minimization, and that an MILP solver implicitly uses a method that is known to be effective in approximating the $\ell_{0}$-norm. Therefore, we conclude that DOPSA displays a remarkable performance by outperforming MILP in terms of solution speed.

The objective (1a) is to maximize the number of calculation points that receive the prescribed dose of 8.5 Gy. Let us first define the vector $y$, indexed by $i$ in $I_{PTV}$ as follows:

[TABLE]

Since $y_{i}$ is zero if calculation point $i$ receives at least the prescribed dose, minimizing the number of non-zero elements in $y$ yields a treatment plan with maximal $V_{100\%}$. The number of non-zero elements is $||y||_{0}$ (recall that the $\ell_{0}$-norm counts the number of non-zero elements in a vector). $D_{10\%}$ and $D_{\max}$ can be rewritten similarly. Thus, the Lagrangean of the MILP (1a)-(1g) is equivalent to an $\ell_{0}$-norm minimization. Attaining sparse vectors by minimizing the $\ell_{0}$-norm has been shown to be NP-hard (Natarajan, 1995).

An MILP solver first solves the following LP relaxation:

[TABLE]

where $\mathcal{X}$ is a polyhedron. Substituting $y_{i}=1-x_{i}$ yields the following equivalent program:

[TABLE]

So, the LP relaxation (2) is equivalent to $\ell_{1}$-norm minimization. Literature, predominately in the field of image processing, suggests that an optimal solution for the $\ell_{1}$-norm minimization (linear relaxation of the MILP) provides a good approximation of an optimal solution for the $\ell_{0}$-norm minimization (MILP): Donoho (2006) shows that, in most cases, the $\ell_{1}$-approximation is also the sparsest solution for underdetermined systems of equations. Li (2010) discusses necessary and sufficient conditions for the equivalence of $\ell_{0}$-norm and $\ell_{1}$-norm solutions in general. These conditions do not hold in our case. Nevertheless, Candes et al. (2005) provide numerical evidence for the power of these $\ell_{1}$-approximations. So, when the MILP (1a)-(1g) is solved with an MILP solver, the solver implicitly uses an approximation that is known to be accurate. This explains the results by Gorissen et al. (2013), where an MILP solver quickly determines good treatment plans. However, our results in Section 3.3 show that DOPSA still outperforms an MILP solver in finding good solutions in a short time-span, which advocates the power of DOPSA.

Compared to MILP-based formulations, DOPSA has another advantage. When more DVH-constraints are added, the time for constraint verifications increases linearly, or sublinearly when solutions are removed in the process. Contrarily, both the number of variables and the number of constraints of the MILP formulation grow linearly in the number of DVH criteria, which slows down an MILP solver at least quadratically. For example, when the LP relaxation is solved with an interior point method, the time complexity of solving the KKT system increases at least quadratically in the number of variables. A dual solver cannot alleviate this slowdown since the number of constraints increases as well.

4 Conclusions

The existing dose-volume based optimization methods MILP and IPIP have not been compared before. Our results show that each method has its own advantages. IPIP is faster while MILP determines a better solution.

DOPSA, our new algorithm, is a viable alternative to both. Its main advantage is that it does not require a costly solver. The time it takes to determine a good solution is less than MILP but slightly more than IPIP. Since IPIP returns a single solution quickly but does not keep improving the solution, DOPSA eventually determines better treatment plans. This leads to two observations: first, advances in computing power will reduce the absolute time difference between both methods, making IPIP lose its edge. Second, it is possible to combine the strength of two methods and use the IPIP solution as a starting point for DOPSA. However, the additional costs for an LP solver required by IPIP and the added complexity is not justified by the slight improvement in solution time.

We have disclosed all information necessary to implement DOPSA. We hope it receives more testing by the community, and will eventually become available to treatment planners.

5 Acknowledgments

The authors would like to thank Elekta Brachytherapy for supplying anonymous patient data, T. Siauw (UCSF, San Francisco CA, USA) and anonymous referees for their constructive feedback.