Fast dose optimization for rotating shield brachytherapy

TL;DR

This paper introduces a rapid convex optimization method using POGS for treatment planning in rotating shield brachytherapy, significantly reducing computation time while maintaining plan quality, thus facilitating clinical adoption.

Contribution

The study presents a novel application of the POGS convex optimization solver for RSBT, achieving faster treatment plan calculations compared to traditional solvers like CPLEX.

Findings

POGS achieved 18 times faster optimization than CPLEX.

Treatment plan quality was comparable between POGS and CPLEX.

The method enables rapid RSBT planning suitable for clinical use.

Abstract

Purpose: To provide a fast computational method, based on the proximal graph solver (POGS) - a convex optimization solver using the alternating direction method of multipliers (ADMM), for calculating an optimal treatment plan in rotating shield brachytherapy (RSBT). RSBT treatment planning has more degrees of freedom than conventional high-dose-rate brachytherapy (HDR-BT) due to the addition of emission direction, and this necessitates a fast optimization technique to enable clinical usage. // Methods: The multi-helix RSBT (H-RSBT) delivery technique was considered with five representative cervical cancer patients. Treatment plans were generated for all patients using the POGS method and the previously considered commercial solver IBM CPLEX. The rectum, bladder, sigmoid, high-risk clinical target volume (HR-CTV), and HR-CTV boundary were the structures considered in our optimization problem, called the asymmetric dose-volume optimization with smoothness control. Dose calculation resolution was 1x1x3 mm^3 for all cases. The H-RSBT applicator has 6 helices, with 33.3 mm of translation along the applicator per helical rotation and 1.7 mm spacing between dwell positions, yielding 17.5 degree emission angle spacing per 5 mm along the applicator.// Results: For each patient, HR-CTV D90, HR-CTV D100, rectum D2cc, sigmoid D2cc, and bladder D2cc matched within 1% for CPLEX and POGS. Also, we obtained similar EQD2 figures between CPLEX and POGS. POGS was around 18 times faster than CPLEX. Over all patients, total optimization times were 32.1-65.4 seconds for CPLEX and 2.1-3.9 seconds for POGS. // Conclusions: POGS substantially reduced treatment plan optimization time around 18 times for RSBT with similar HR-CTV D90, OAR D2cc values, and EQD2 figure relative to CPLEX, which is significant progress toward clinical translation of RSBT. POGS is also applicable to conventional HDR-BT.

Tables

Table 1. Table 1: HR-CTV volumes and dimensions for all patients

Patient Num.	HR-CTV volume ($c{m}^3$)	HR-CTV maximum dimension ($cm$)
Case 1	42.2	6.3
Case 2	45.8	7.4
Case 3	78.0	8.6
Case 4	98.8	9.6
Case 5	75.0	7.5
Avg.	68.0	7.9
SDa	23.8	1.8

Table 2. Table 2: Parameter settings

Method	${\hat{d}}_{tumor}$	${\hat{d}}_{bladder}$	${\hat{d}}_{rectum}$	${\hat{d}}_{sigmoid}$	${\hat{d}}_{normal}$a	Tumor	Bladder	Rectum	Sigmoid	Normalb	$\beta$
						${\lambda }_i^{+}$/ ${\lambda }_i^{-}$	${\lambda }_i^{+}$/ ${\lambda }_i^{-}$	${\lambda }_i^{+}$/ ${\lambda }_i^{-}$	${\lambda }_i^{+}$/ ${\lambda }_i^{-}$	${\lambda }_i^{+}$/ ${\lambda }_i^{-}$
CPLEX	40	25	20	20	40	0/ 2	2/ 0	2/ 0	2/ 0	2/ 0	100
POGS	40	25	20	20	40	0/ 2	2/ 0	2/ 0	2/ 0	2/ 0	100

Table 3. Table 3: Comparison between POGS and CPLEX for 45 °azimuthal angle

Case	Method	HR-CTV	HR-CTV	Bladder	Rectum	Sigmoid	Execution
		$D_{90}$ (Gy)	$D_{100}$ (Gy)	$D_{2cc}$ (Gy)	$D_{2cc}$ (Gy)	$D_{2cc}$ (Gy)	time (sec.)
Case 1	CPLEX	110.8	54.0	89.9	62.4	75.0	32.1
	POGS	111.4	54.0	90.0	64.7	74.7	2.1
Case 2	CPLEX	111.5	44.3	90.0	72.2	54.4	37.0
	POGS	111.5	44.3	90.0	71.7	54.8	2.1
Case 3	CPLEX	96.0	44.3	85.9	57.3	75.0	65.4
	POGS	95.0	44.3	85.2	55.1	75.0	3.9
Case 4	CPLEX	107.0	55.3	90.0	69.9	54.0	39.4
	POGS	106.9	55.4	90.0	69.8	54.0	2.3
Case 5	CPLEX	112.7	44.3	90.0	68.1	59.2	65.4
	POGS	112.7	44.3	90.0	68.1	59.2	3.2
Average	CPLEX	107.6	48.4	89.2	66.0	63.5	47.9
	POGS	107.5	48.5	89.0	65.9	63.5	2.7

Table 4. Table 4: Dimension of D ∈ ℝ l × m ​ n 𝐷 superscript ℝ 𝑙 𝑚 𝑛 D\in\mathbb{R}^{l\times mn} in Eqn. ( II.3.1 )

Patient Num.	$l\times mn$
Case 1	$54693\times 144$
Case 2	$51109\times 126$
Case 3	$79065\times 222$
Case 4	$50680\times 144$
Case 5	$59220\times 228$
Avg.	$58953\times 173$

Equations

$$d_i=\sum _{j,b}{D}_i(j,b)t_{j,b},$$

$$\underset{t\in {\mathbb{R}}^{mn},d\in {\mathbb{R}}^l}{\mathrm{minimize}}\sum _{i\in {\mathcal{I}}_{VOIs}}h(d_i)+\sum _{j=1}^m\beta ||Lt[j]|{|}_1$$

$$\text{subject to}Dt=d,$$

$$t\ge 0,$$

$$h(d_i)\triangleq ({\lambda }_i^{+}H(d_i-\widehat{d_i})+{\lambda }_i^{-}H(\widehat{d_i}-d_i))(d_i-\widehat{d_i}{)}^2,$$

$$L\triangleq \left[\begin{array}{ccccccc}1& -1& 0& 0& ...& 0& 0\\ 0& 1& -1& 0& ...& 0& 0\\ 0& 0& 1& -1& ...& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 0& ...& 1& -1\\ 0& 0& 0& 0& ...& 0& 0\end{array}\right].$$

$$\underset{t,y\in {\mathbb{R}}^{mn},d\in {\mathbb{R}}^l}{\mathrm{minimize}}\sum _{i\in {\mathcal{I}}_{VOIs}}h(d_i)+\beta ||y|{|}_1+I(x\ge 0)$$

$$\text{subject to}Dt=d,$$

$$\bar{L}t=y,$$

$$A=\left(\begin{array}{c}D\\ \bar{L}\end{array}\right),z=\left(\begin{array}{c}d\\ y\end{array}\right),$$

$$\underset{t\in {\mathbb{R}}^{mn},z\in {\mathbb{R}}^{l+mn}}{\mathrm{minimize}}\sum _{i=1}^{l}h(z_i)+\beta ||z_{[l+1:l+mn]}|{|}_1+I(t\ge 0)$$

$$\text{subject to}At=z,$$

$$g(t)=I(t\ge 0),$$

$$f(z)=\sum _{i=1}^{l}h(z_i)+\beta ||z_{[l+1:l+mn]}|{|}_1.$$

$${\hat{t}}^{k+1}$$

$$=\underset{t}{\mathrm{argmin}}I(t\ge 0)+\frac{\rho }{2}||t-(t^k-{\tilde{t}}^k)|{|}^2,$$

$${\hat{z}}^{k+1}$$

$$=\underset{z}{\mathrm{argmin}}\sum _{i=1}^{l}h(z_i)+\beta ||z_{[l+1:l+mn]}|{|}_1$$

$$+\frac{\rho }{2}||z-(z^k-{\tilde{z}}^k)|{|}^2.$$

$${\hat{t}}^{k+1}=max(t^k-{\tilde{t}}^k,\text{0}),$$

$${\hat{z}}_i^{k+1}=\{\begin{array}{ll}& z_i^k-{\tilde{z}}_i^k-\frac{\beta }{\rho },\text{if}{z}_i^k-{\tilde{z}}_i^k>\frac{\beta }{\rho }\\ & z_i^k-{\tilde{z}}_i^k+\frac{\beta }{\rho },\text{if}{z}_i^k-{\tilde{z}}_i^k<-\frac{\beta }{\rho }\\ & 0,\text{otherwise}\end{array}.$$

$${\hat{z}}_i^{k+1}=\{\begin{array}{ll}& \frac{2{\lambda }_i^{+}{\hat{d}}_i+\rho (z_i^k-{\tilde{z}}_i^k)}{2{\lambda }_i^{+}+\rho },\text{if}{z}_i^k-{\tilde{z}}_i^k\ge {\hat{d}}_i,\\ & \frac{2{\lambda }_i^{-}{\hat{d}}_i+\rho (z_i^k-{\tilde{z}}_i^k)}{2{\lambda }_i^{-}+\rho },\text{if}{z}_i^k-{\tilde{z}}_i^k<{\hat{d}}_i.\end{array}.$$

$$\underset{t,z}{\mathrm{minimize}}\frac{1}{2}||t-({\hat{t}}^{k+1}+{\tilde{t}}^k)|{|}_2^2+\frac{1}{2}||z-({\hat{z}}^{k+1}+{\tilde{z}}^k)|{|}_2^2$$

$$\text{subject to}At=z$$

$$\left(\begin{array}{c}{t}^{k+1}\\ z^{k+1}\end{array}\right)={\left(\begin{array}{cc}I& A^T\\ A& -I\end{array}\right)}^{-1}\left(\begin{array}{c}{\hat{t}}^{k+1}+{\tilde{t}}^k+A^T({\hat{z}}^{k+1}+{\tilde{z}}^t)\\ 0\end{array}\right).$$

$${\tilde{t}}^{k+1}={\tilde{t}}^k+{\hat{t}}^{k+1}-t^{k+1},$$

$${\tilde{z}}^{k+1}={\tilde{z}}^k+{\hat{z}}^{k+1}-z^{k+1}.$$

$$||A{\hat{t}}^{k+1}-{\hat{z}}^{k+1}|{|}_2\le {ϵ}^{pri},$$

$$||A^T{\hat{v}}^{k+1}+{\hat{\mu }}^{k+1}|{|}_2\le {ϵ}^{dual},$$

$${ϵ}^{pri}={ϵ}^{abs}+{ϵ}^{rel}||{\hat{z}}^{k+1}|{|}_2,$$

$${ϵ}^{dual}={ϵ}^{abs}+{ϵ}^{rel}||{\hat{\mu }}^{k+1}|{|}_2,$$

Full text

Methods: The multi-helix RSBT (H-RSBT) delivery techniqueDadkhah et al. (2015) was considered with five representative cervical cancer patients. Treatment plans were generated for all patients using the POGS method and the previously considered commercial solver IBM ILOG CPLEXIBM . The rectum, bladder, sigmoid colon, high-risk clinical target volume (HR-CTV), and HR-CTV boundary were the structures considered in our optimization problem, which is called the asymmetric dose-volume optimization with smoothness control. Dose calculation resolution was $1\times 1\times 3$ $mm^{3}$ for all cases. The H-RSBT applicator has 6 helices, with $33.3\;mm$ of translation along the applicator per helical rotation and $1.7\;mm$ spacing between dwell positions, yielding 17.5°emission angle spacing per $5\;mm$ along the applicator.

Results: For each patient, HR-CTV $D_{90}$, HR-CTV $D_{100}$, rectum $D_{2cc}$, sigmoid $D_{2cc}$, and bladder $D_{2cc}$ matched within $1\%$ for CPLEX and POGS methods. Also, we obtained similar EQD2 figures between CPLEX and POGS methods. POGS was around 18 times faster than CPLEX. Over all patients, total optimization times were 32.1-65.4 seconds for CPLEX and 2.1-3.9 seconds for POGS.

Conclusions: POGS substantially reduced treatment plan optimization time around 18 times for RSBT with similar HR-CTV $D_{90}$, organ at risk (OAR) $D_{2cc}$ values, and EQD2 figure relative to CPLEX, which is significant progress toward clinical translation of RSBT. POGS is also applicable to conventional high-dose-rate brachytherapy.

Brachytherapy, rotating shield brachytherapy, gynecological cancer, cancer treatment planning, optimization

I INTRODUCTION

High-dose-rate brachytherapy (HDR-BT) involves placing a radiation source inside of or adjacent to a target organ, i.e., tumor. Conventional HDR-BT uses an unshielded brachytherapy source with a radially-symmetric dose distributionPȯtter et al. (2007); Kirisits et al. (2006), which limits the intensity modulation capacity of the approach. Rotating-shield brachytherapy (RSBT) has a rotating radiation-attenuating shield around a brachytherapy source. The RSBT concepts for single-catheter treatmentEbert (2002) and multi-catheter treatmentEbert (2006) were introduced by Ebert in 2002 and 2006 respectively.

In the multi-helix RSBT (H-RSBT) treatment, a radiation source travels inside a brachytherapy applicator having helical keyways. While moving along the applicator for a given keyway, the partial shield rotates around the radiation source simultaneously. In traveling along each keyway, the radiation source stops at designated locations called dwell positions. By adjusting the distance between adjacent dwell positions, the rotation angle of the partial shield is determined accordingly. Intensity modulated dose distributions can be delivered to the target with reduced dose exposure to non-target organs by controlling the treatment time in an optimal manner for each dwell position. Hence, it is reported that a radiation source with rotating shields can deliver more conformal dose distributions than an unshielded radiation source.Dadkhah et al. (2015)

HDR-BT treatment plans are often generated using inverse planning toolsLessard and Pouliot (2001); Akimoto et al. (2006); Baltas et al. (2009); Citrin et al. (2005); DeWitt et al. (2005); Manikandan et al. (2013); Jacob et al. (2008). Based on the given clinical prescription, various optimization problems were introduced previously ranging from minimizing treatment time under restrictionsKneschaurek, Schiessl, and Wehrmann (1999); Renner, O Connor, and Bermudez (1990) to minimizing dose errorLessard and Pouliot (2001); Manikandan et al. (2013); Jacob et al. (2008); Alterovitz et al. (2006). Inverse planning by simulated annealing (IPSA)Lessard and Pouliot (2001) is a well known method to optimize the dose volume histogram (DVH) directly with given constraints. The BrachyVision treatment planning system (Varian Medical System Inc., Palo Alto, CA) uses this type of DVH-based optimization algorithm.Manikandan et al. (2013)

Unlike the conventional HDR-BT optimization, the RSBT optimization problem has the additional optimization variables of radiation exposure time at each angle of the shield. Due to the increased degrees of freedom in RSBT, RSBT optimization is more difficult than that for the conventional HDR-BT. In addition, there is a compelling need to quickly obtain optimal treatment plans in RSBT to enable clinical usage. To achieve this, researchers have used the dose-surface optimization (DSO) methodLiu et al. (2013); Yang et al. (2013), which minimizes the total dose errors over only voxels on the HR-CTV surface. Instead of dealing with only voxels on the HR-CTV surface, Liu et al.Liu et al. (2014) defined the region of interest in tumor which includes the surface of HR-CTV as well as the inside voxels of tumor. Additionally, the authors used the total variation (TV) norm penalty in their optimization problem to make smooth changes in the emission times of adjacent beams in the treatment process to facilitate the efficient delivery of an RSBT plan. This optimization problem for RSBT is called asymmetric dose-volume optimization with smoothness control (ADOS).

In this paper, we consider the ADOS optimization problem. A fast computational method is proposed to solve the ADOS optimization problem for the optimal cancer treatment planning for RSBT. Liu et al. used a commercial optimization solver called CPLEXLiu et al. (2014). In order to efficiently solve the ADOS optimization problem, which is a large-scale RSBT optimization problem, we designed an optimization method based on the proximal graph solver (POGS)Fougner and Boyd (2015), which is a solver using the alternating direction method of multipliers (ADMM). For using POGS, we derived closed-form formulas for the proximal operators used in POGS. Further, we applied our method to the H-RSBT, which is a mechanically-feasible delivery technique for RSBT proposed by Dadkhah et al.Dadkhah et al. (2015). In the numerical experiments, we considered cervical cancer, even though our method is also applicable to other types of cancer such as breast cancer and prostate cancer.

II MATERIALS AND METHODS

II.1 Delivery method

In order to deliver the radiation dose to a target organ, we consider a mechanically-feasible delivery technique for RSBT, called the multi-helix RSBT (H-RSBT) Dadkhah et al. (2015). Fig. 1 shows the illustration of the H-RSBT method. The shield opening is represented by the azimuthal and zenith emission angles, denoted by $\Delta\varphi$ and $\Delta\theta$ respectively.

In H-RSBT, a Xoft Axxent radiation source, inside its cooling catheter, with a freely-rotating partial radiation shield is translated inside an applicator with six helical keyways carved out of the inner wall. The six keyways are evenly spaced on the applicator cross section, by $60\degree$, and each keyway has a helical pitch of one rotation per $33.3\;mm$ of translation. The partial shield has a protruding key that travels down a given keyway, and, due to the helical design of the keyways, the shield rotates about the radiation source as the source catheter is translated, and the emission angle of the shield is known for a given keyway and translational dwell position. As the H-RSBT applicator has 6 helices, with $33.3\;mm$ of translation along the applicator per helical rotation and $1.7\;mm$ spacing between dwell positions, it yields 17.5°of rotation for the shield per $5\;mm$ (standard dwell position spacing) of its translation along the applicator. The dose calculation resolution was $1\times 1\times 3\;mm^{3}$ for all cases. The transmission through the shield is 0.1% and approximated to be 0% for the dose calculation. For each patient considered, 45°azimuthal emission angle was used for treatment planning. The zenith emission angle of the modeled shields was an asymmetric 116°, which is consistent with previous work.Dadkhah et al. (2015)

II.2 Radiation source model and dose calculation

For H-RSBT, the delivery is parameterized by keyway number and dwell position number along the keyway. To quantitatively describe the structure of high-dose-regions formed by a partially shielded source, we introduce the notation of beamlet. A beamlet, denoted by $D_{i}(j,b)$, is defined as the dose rate at point $\vec{r}_{i}$ with the shielded source positioned at the $j$-th dwell position while the shield is aligned with the $b$-th keyway.

To calculate the beamlet, we use the TG-43 dose calculation model of Rivard et al.Rivard et al. (2006). The radiation source is assumed to be partially shielded 50kVp Xoft AxxentTM (Sunnyvale, CA). To be consistent with previous workLiu et al. (2013); Yang et al. (2013); Velpula et al. (2012), we consider that the dose to the points blocked by the shield is 0, since the transmission rate from 50kVp Xoft AxxentTM can be controlled to be less than 0.1% when using a $0.5\;mm$ tungsten shield.Yang et al. (2013); Liu et al. (2013) Then, we can quantify the radiation dose amount at the point $\vec{r}_{i}$, denoted by $d_{i}$, as a time-weighted sum of all beamlets as follows:

[TABLE]

where $t_{j,b}$ is the duration time of the beamlet $D_{i}(j,b)$.

In the next subsection, we introduce the RSBT optimization problem having asymmetric penalty parameters for HR-CTV and organs at risk (OARs), with TV regulation term for smoothness in the beamlet emission times.

II.3 Optimization problem for cancer treatment planning in RSBT

II.3.1 Problem formulation

Let $t\in\mathbb{R}^{mn\times 1}$ be the beamlet emission time vector for all keyways and all dwell positions, where $m$ and $n$ are the number of keyways, i.e, $m=6$, and the number of dwell positions along a keyway respectively. We can obtain $t$ by vectorizing $t_{j,b}$ in Eqn. (1); namely, the vector $t\in\mathbb{R}^{mn\times 1}$ is a concatenated vector, which is expressed as $t=[t[1]^{T},t[2]^{T},...,t[m]^{T}]^{T}$, where $t[j]\in\mathbb{R}^{n\times 1}$ is the dwell time vector for all the beamlets along a keyway, and the super-script $T$ represents the transpose. Let us define a concatenated dose rate matrix $D=[D[1],D[2],...,D[m]]\in\mathbb{R}^{l\times mn}$, where $D[j]\in\mathbb{R}^{l\times n}$, $j=1,...,m$, is the dose rate matrix for the $j$-th selected keyway, and $l$ is the number of voxels that we are interested in. We denote the whole index set for voxels of interest (VOIs) as $\mathcal{I}_{VOIs}$ and the index set for HR-CTV, bladder, rectum, sigmoid, and normal tissue around the HR-CTV as $\mathcal{I}_{tumor}$, $\mathcal{I}_{bladder}$, $\mathcal{I}_{rectum}$, $\mathcal{I}_{sigmoid}$, and $\mathcal{I}_{normal}$ respectively. The dose rate matrix $D$ has information about dose rate delivered to each tissue point in $\mathcal{I}_{VOIs}$ from each beamlet.

We consider the following RSBT optimization problem having a quadratic objective function with total variation (TV) regulation term for smooth beamlet emission times:

[TABLE]

where $D\in\mathbb{R}^{l\times mn}$ is a dose rate matrix, $t\geq 0$ is the element-wise non-negative emission time, $d_{i}$, which is the $i$-th element of $d$, is the dose amount at the $i$-th voxel as introduced in Eqn. (1). In Eqn. (II.3.1),

[TABLE]

where $H(x)$ is the unit step function, which is $H(a)=1$ if $a>0$, and $H(a)=0$ if $a\leq 0$, and $\lambda_{i}^{+}$ and $\lambda_{i}^{-}$ represent overdose and underdose penalty parameters for the $i$-th voxel respectively. $L\in\mathbb{R}^{n\times n}$ is the matrix which calculates TV norm of a vector; namely, $L$ is defined as follows:

[TABLE]

$\hat{d_{i}}$ is a prescribed dose amount for the $i$-th voxel. $\hat{d_{i}}$ can have a different value for each VOI. For example, $\hat{d_{i}}=\hat{d}_{tumor}$ if $i\in\mathcal{I}_{tumor}$, and $\hat{d_{i}}=\hat{d}_{badder}$ if $i\in\mathcal{I}_{bladder}$. We denote the prescribed dose amount for HR-CTV, bladder, rectum, sigmoid, and normal tissue around the HR-CTV as $\hat{d}_{tumor}$, $\hat{d}_{bladder}$, $\hat{d}_{rectum}$, $\hat{d}_{sigmoid}$, and $\hat{d}_{normal}$. For $\lambda_{i}^{+}$ and $\lambda_{i}^{-}$, we use different non-negative overdose and underdose parameter values including 0.

The TV norm alleviates the positioning uncertainty in the treatment process. If we have two dramatically different emission times $t_{j,b}$ and $t_{j+1,b}$ between two adjacent beamlets along the same keyway, a small error in the dwell positions may cause an unacceptable treatment result. By applying the smoothness term between two adjacent beamlets along the same keyway in H-RSBT, we can reduce the treatment error caused by the positioning uncertainty in the treatment process.

Since we use different penalty parameter values for the overdose and underdose of a voxel, we call Eqn. (II.3.1) as RSBT optimization problem having asymmetric penalty parameters or simply ADOS.

II.3.2 POGS implementation

In order to simplify the sum of the TV norms in Eqn. (II.3.1), let us introduce a matrix $\bar{L}\triangleq I_{m\times m}\otimes L$, where $\otimes$ is the Kronecker product, and $I_{m\times m}$ is an $m\times m$ identity matrix. By assigning $\bar{L}t=y$ and introducing the indicator function $I(\cdot)$, we restate Eqn. (II.3.1) as follows:

[TABLE]

where $D\in\mathbb{R}^{l\times mn}$, $\bar{L}\in\mathbb{R}^{mn\times mn}$, and $I(t\geq 0)$ is the element-wise indicator function; namely, $I(t_{i}\geq 0)=0$ if $t_{i}\geq 0$, and $I(t_{i}\geq 0)=\infty$ if $t_{i}<0$.

By letting

[TABLE]

we can further simplify Eqn. (II.3.2) into

[TABLE]

where $z_{[a:b]}$ is the partial vector of $z$ by taking vector $z$ from the $a$-th element to the $b$-th element. We define the following functions:

[TABLE]

Then we turn Eqn. (II.3.2) into a graph-form convex optimization problemFougner and Boyd (2015), where the constraint is $z=At$, and $A=[D^{T}\bar{L}^{T}]^{T}\in\mathbb{R}^{(l+mn)\times mn}$. We have derived the detailed updating rules for each optimization variable in the POGS solver for Eqn. (II.3.2). We introduce our derived results in detail for the proximal operators used in the POGS solver, updating steps, and stopping criteria as follows.

The POGS updates primal variables, conducts the projection onto the space $z=At$, and then, updates dual variables iteratively until the stopping criteria are satisfied or the maximum number of iterations, denoted by $MaxItr$, is reached. The primal variable and dual variable are updating variables to be used for optimality condition in the algorithm. For the primal variable, dual variable, and projection result, we used $(\hat{t},\hat{z})$, $(t,z)$, and $(\tilde{t},\tilde{z})$ respectively. We introduce each updating steps in detail for our optimization problem in Eqn. (II.3.2). We use the super-script $k$ to represent the $k$-th iteration.

Updating primal variables $\hat{t}^{k+1}$ and $\hat{z}^{k+1}$: In updating the primal variables, we use the following proximal operators with a penalty parameter $\rho$:

[TABLE]

The proximal operator is used to make a compromise between the solution at the $k$-th iteration and the function value with the solution at the $k+1$ iteration. We are able to explicitly derive closed-form formulas for the proximal operators. For $\hat{t}^{k+1}\in\mathbb{R}^{mn}$, we have

[TABLE]

where $max(a,b)$ provides the maximum value between $a$ and $b$ element-wise. For $\hat{z}^{k+1}_{i}$, $1\leq i\leq l$, we also derive

[TABLE]

For $\hat{z}^{k+1}_{i}$, $l+1\leq i\leq l+mn$, we obtain

[TABLE]

POGS uses the adaptive value for $\rho$ as default to further increase the convergence speed.

Projection onto $z=At$ from $(\hat{t}^{k+1}+\tilde{t}^{k},\hat{z}^{k+1}+\tilde{z}^{k})$: The projection operation is mapping the primal variables to the closest feasible solution. The projected variables onto $z=At$ from $(\hat{t}^{k+1}+\tilde{t}^{k},\hat{z}^{k+1}+\tilde{z}^{k})$, denoted as $t^{k+1}$ and $z^{k+1}$, are obtained by solving the following optimization:

[TABLE]

By solving this optimization and using Lagrange conditionsBoyd and Vandenberghe (2004), we have the following formulation:

[TABLE]

Updating dual variables $\tilde{t}^{k+1}$ and $\tilde{z}^{k+1}$: We obtain the dual variable at iteration $(k+1)$ by updating the dual variable at the $k$-th iteration as follows:

[TABLE]

We summarize the updating steps in Algorithm 1.

Stopping criteria: For the stopping criteria, we define the primal and dual residuals as follows:

[TABLE]

where $\hat{v}^{k+1}=-\rho(\hat{z}^{k+1}-z^{k}+\tilde{z}^{k})$, $\hat{\mu}^{k+1}=-\rho(\hat{t}^{k+1}-t^{k}+\tilde{t}^{k})$. Here, $\epsilon^{pri}$ and $\epsilon^{dual}$ are positive tolerances for primal and dual residuals respectively:

[TABLE]

where we used $(\epsilon^{abs},\epsilon^{rel})=(10^{-4},10^{-2})$ in the numerical experiments.

II.4 Treatment planning

Five patients with cervical cancer were considered, whose HR-CTV volumes range from $42.2$ to $98.8\;cm^{3}$. Table 1 shows the volume and maximum dimension of HR-CTV for five patients. All the HR-CTVs and OARs were manually contoured by physicians on T2 weighted $1\times 1\times 3$ $mm^{3}$ resolution MR images taken with a Siemens MAGNETOM 3T scanner (Siemens, Germany) at the beginning of the first fraction of brachytherapy. A titanium Fletcher-Suit-Delclos style tandem and ovoids (Varian Medical Systems, Palo Alto, CA) were used as the brachytherapy applicator. We used the same datasets as the previous research conducted by Liu et al.Liu et al. (2014) and Dadkhah et al.Dadkhah et al. (2015).

All the patients had external beam radiation treatment in 25 fractions at 1.8 Gy/fraction. We assume that the external beam radiotherapy dose was uniformly delivered to the HR-CTV and OARs for all the patient cases. The dose in each voxel was converted to equivalent dose in 2 Gy per fraction of external radiation therapy (EQD2) using the linear quadratic model,Joiner and Van der Kogel (2016) where the linear-quadratic parameter, $\nicefrac{{\alpha}}{{\beta}}$, is set to 3 Gy for OARs and 10 Gy for HR-CTV.

For VOIs, we define the voxels located at a distance between $3\;mm$ and $20\;mm$ to the radiation source path or those within $10\;mm$ inside and outside of the HR-CTV boundary surface.Liu et al. (2014) We deal with the HR-CTV, HR-CTV boundary, bladder, sigmoid, and rectum inside of VOIs in our optimization problem. The optimization parameter settings are shown in Table 2.

For all the brachytherapy treatment plans, we escalated the EQD2 of the HR-CTV without exceeding the $D_{2cc}$ tolerance of the bladder, rectum, and sigmoid colon. We used 90 Gy for bladder tolerance, and 75 Gy for rectum and sigmoid colon tolerances according to Groupe Européen de Curiethérapie, European Society for Therapeutic Radiology and Oncology (GEC ESTRO).Pȯtter et al. (2006); Haie-Meder et al. (2005)

II.5 Treatment plan evaluation

Optimized treatment plans were generated for all patients using the POGS method and the previously considered CPLEX method.Liu et al. (2014) The same objective function, with the same input parameters and beamlets, was minimized for each patient with both methods. A total variation term was included in the objective function as a regularization term, resulting in smoothly-varying emission times along each keyway. The regularization promotes robustness of the resulting overall dose distribution with respect to small errors (expected $\leq 1\;mm$) in source positioning. The rectum, bladder, sigmoid colon, HR-CTV, and HR-CTV boundary were the structures considered.

We compared our method with the previous research conducted by Liu et al.Liu et al. (2014) using CPLEXIBM for H-RSBT. We evaluated the quality of the delivery plans as well as the execution time to solve Eqn.  (II.3.1) with POGSFougner and Boyd (2015). Since Liu et al. compared their method based on CPLEX with other existing RSBT dose optimization methods ranging from DSO to IPSA in their previous researchLiu et al. (2014), we only compared POGS and CPLEX in this paper.

The comparison metrics for the quality of the delivery plans are the HR-CTV $D_{90}$, HR-CTV $D_{100}$, OARs $D_{2cc}$, DVH, and dose distributions. Since the goal of this research is achieving a fast solution to the RSBT dose optimization problem without compromising the plan quality, we compared the execution times to solve Eqn.  (II.3.1) for all five patient cases. We conducted our numerical experiments on HP Z220 CMT with an Intel Core i7-3770 dual core CPU @3.4GHz clock speed and 16GB DDR3 RAM, using Matlab (R2013b) on the Windows 7 operating system.

III RESULTS

Results for five cervical cancer patient cases are shown in Table 3. Table 4 shows the dimension of the dose matrix $D$ in Eqn.  (II.3.1). Fig. 3 shows the corresponding DVH in H-RSBT. With the same parameter settings as in Table 2, POGS can achieve an RSBT plan with almost the same $D_{90}$ value (less than 1% difference) as that achieved by CPLEX in each of the five patient cases. $D_{2cc}$ values for OARs obtained by POGS are also almost the same (less than 1% difference) as those from CPLEX. For the execution time, we achieved around 18 times faster speed to solve the ADOS problem for H-RSBT than the CPLEX based method on average. Over all patients, total optimization times were 32.1-65.4 seconds for CPLEX and 2.1-3.9 seconds for POGS.

Fig. 2 shows that the EQD2 figures were similar for each case between CPLEX and POGS.

IV DISCUSSION

Various treatment planning methods in radiation therapy have been studied. One of the well known methods, which is called IPSA, was introduced to directly optimize DVH with given constraints in heuristic way. Due to its heuristic nature, a global solution is not guaranteed. Instead of directly optimizing DVH, we consider the voxel-wise optimization problem for the RSBT treatment planning, called ADOS optimization problem, which can be expressed in a convex optimization problem. We can take advantage of the convexity to obtain a global solution.

In the ADOS optimization problem, we reduced the size of the ADOS optimization problem by defining VOIs. Instead of defining VOIs, the whole voxels can be considered in the optimization problem under the expectation of better treatment quality with heavy computation. Since parallel computing and GPU-based high performance computing can play an important role in solving extremely large-scale optimization problems, there is a rising question about the usability of POGS in parallel computing environment or GPU-based implantation. The implementation of POGS in such environment is another research work.

In addition, we used a partial shield with 45°azimuthal angle in H-RSBT for our numerical experiments. However, finding the optimal shield angle in H-RSBT is still an open problem. In order to determine the size of angle and the radiation exposure time at each angle of shield, considering both variables in the optimization problem is also a possible optimization problem.

Finally, POGS (and ADMM) was used in previous research on intensity-modulated radiation therapy (IMRT)Fougner and Boyd (2015); Zarepisheh, Xing, and Ye (2017), fluence map optimizationGao (2016), and external beam radiotherapy (EBRT) optimizationLiu et al. (2017). Right at the time of submitting our journal manuscript, we learned of the recently-appearing (published on April 12th, 2017) work Liu et al. (2017) which applied POGS algorithms to EBRT dose optimization. The work Liu et al. (2017) focused on EBRT, while our paper is the first work to use POGS in brachytherapy, including the mechanically-feasible delivery system called H-RSBT. In our paper, we use TV norm to promote smoothly-varying emission times along each keyway, such that the treatment plan is robust to positioning errors of dwell positions. By comparison, the TV norm is instead applied to promote the smoothness of the resulting fluence map in the researchLiu et al. (2017) and simplify the delivery. Our proposed method is applicable to conventional HDR-BT as well as dynamic modulated brachytherapyWebster et al. (2011) with simple modifications, since they share similar mechanisms as H-RSBT.

V CONCLUSIONS

POGS substantially reduced treatment plan optimization time around 18 times for RSBT with similar HR-CTV $D_{90}$, OAR $D_{2cc}$ values, and EQD2 figure comparing to CPLEX, which is significant progress toward clinical translation of RSBT. Over all cervical cancer patients, total optimization times were 32.1-65.4 seconds for CPLEX and 2.1-3.9 seconds for POGS. POGS is also applicable to conventional high-dose-rate brachytherapy.

VI ACKNOWLEDGEMENT

The authors would like to thank the generous support from NIH 1R01EB020665 and NSF DMS-1418737. The authors would also like to thank Jianfeng Cai, Qihang Lin, and Haiye Huo for discussions and advice on paper writing.