Inverse treatment planning for an electronic brachytherapy system delivering anisotropic radiation therapy

The dose distributions produced at each of the 54 operating points were calculated by Monte Carlo simulations using the Geant4 toolkit, version 10.4p2 (Geant4 Collaboration 2017). These simulations have been extensively validated against ion chamber measurements for this device, as described in (Badali et al 2019). In the simulations, the medium was assumed to be spatially homogeneous and composed of water. The output of the simulations is a 3D voxelized volume, where each voxel stores the dose accumulation per electron history for that region of space. By scaling the simulation output by the electron current at a specific beam power and kVp, the 3D voxelized volume is converted to units of Gy s⁻¹. The treatment dose distribution in Gy for a treatment plan can then be calculated by simply adding the dose distribution for each operating point, weighted by the corresponding exposure time.

The complete inverse RTP algorithm is illustrated in figure 2. The most important step in the inverse RTP algorithm is the optimization step, which calculates the 'best' treatment plan, given the target isodose surface and optimization goals. Mathematically, the optimized treatment plan $\hat{{\boldsymbol{t}}}$ is calculated according to:

$$\begin{eqnarray}&&\hat{{\boldsymbol{t}}}=\mathop{{\rm{\arg }}\,{\rm{\min }}}\limits_{{\boldsymbol{t}}\succcurlyeq 0}\left\{f\left({\boldsymbol{t}}\right)+\mu R\left({\boldsymbol{t}}\right)\right\},\end{eqnarray} \tag{ 1 }$$

where $f\left({\boldsymbol{t}}\right)$ is the objective function which penalizes treatment plans that differ from the PTV surface, $R\left({\boldsymbol{t}}\right)$ is a regularization function with a strength μ that will be discussed in section 2.1.2.2, and $\succcurlyeq$ denotes element-wise non-negativity (since the exposure times must be either zero or a positive number). The vector t contains the exposure times (in seconds) for each operating point at each translation position. Finally, the quality of the treatment plan is quantified by the calculation of metrics discussed in section 2.2.

Zoom In Zoom Out Reset image size

2.1.1. Preprocessing

2.1.2. Optimization

2.1.2.1. Objective function

There are many possible objective functions which can be used to quantify the agreement between the PTV and the dose distribution of a given treatment plan (see, for example, De Boeck et al 2014). Although many RTP systems use a data fidelity term that does a voxel-to-voxel comparison of the dose (Webb 1991, Sloboda 1992, Oldham and Webb 1995, Xing et al 1999, Cotrutz and Xing 2003), such objective functions do not penalize dose delivered to the healthy tissue adjacent to the margins. Instead, we define the objective function as:

$$\begin{eqnarray}&&f\left({\boldsymbol{t}}\right)=1-CN\left({\boldsymbol{t}}\right),\end{eqnarray} \tag{ 2 }$$

where $CN\left({\boldsymbol{t}}\right)$ is the conformation number (CN) introduced by van't Riet et al (1997). It is defined as:

$$\begin{eqnarray}&&CN=\displaystyle \frac{PT{V}_{TP}}{PTV}\times \displaystyle \frac{PT{V}_{TP}}{{V}_{TP}},\end{eqnarray} \tag{ 3 }$$

where $PTV$ is the volume of the target PTV in voxels, ${V}_{TP}$ is the size of the treatment plan's isodose volume in voxels at the prescribed dose, and $PT{V}_{TP}$ is the volume of their overlap (that is, the part of the PTV covered by the treatment plan's isodose volume). The first term represents the proportion of the PTV that receives a dose equal to or greater than the prescribed dose (e.g. $PT{V}_{TP}/PTV=1$ means that the PTV is fully encompased by the treatment plan's isodose volume). The second term represents the proportion of treatment plan's isodose volume that overlaps with the PTV (e.g. $1-PT{V}_{TP}/{V}_{TP}$ is the fraction of the treatment plan's isodose volume that is delivering dose to healthy tissue, at or above the prescribed dose). The combination of these two terms simultaneously permit treatment of the target (first term) and protection of healthy tissue (second term). The CN ranges from 0 to 1. Accordingly, the objective function defined in equation (2) encourages treatment plans with CNs close to 1, which corresponds to perfect conformation between the target and treatment volumes (Feuvret et al 2006). The CN is also referred to in the literature as the conformal index (COIN) (Baltas et al 1998).

We note that objectives based on a clinical index such as the CN have been successfully used in other inverse RTP implementations (for a recent review see De Boeck et al 2014). To our knowledge, the CN itself has been used as the objective in at least one other study (Lahanas et al 1999). The advantage of such objectives is that they are intuitively easy to understand for the medical physicist (Lahanas et al 2001), since they are the same metrics used to clinically evaluate the quality of the plan.

Additionally, we want to mention that using the CN in the objective function only considers the shape of the isodose volume and does not consider the dose within volume. Although this results in treatment plans with considerable dose inhomogeneity within the PTV, because the intended use case is to irradiate surgical margins, such inhomogeneities are likely to have minimal effects on such small volumes.

2.1.2.2. Regularization

A regularization term was added to the objective in equation (2) to penalize treatment plans that require a long treatment time (e.g. by using many operating points). Because the device is intended to irradiate surgical margins, such plans are undesirable because they increase the operating room time and associated costs. The regularization was thus chosen to be:

$$\begin{eqnarray}&&R\left({\boldsymbol{t}}\right)={\parallel {\boldsymbol{t}}\parallel }_{1},\end{eqnarray} \tag{ 4 }$$

where ${\parallel \cdot \parallel }_{1}$ is the L1 norm. Since we force all the exposure times to be positive, the L1 norm is simply the total treatment time, and the regularizer term thus encourages short treatment times. The regularization strength $\mu \geqslant 0$ in equation (1) allows for tuning the degree of regularization, which allows balancing of fidelity to the treatment plan against shorter treatment times.

2.1.2.3. Simulated annealing

The optimal exposure times were calculated by solving equation (1) with the simulated annealing algorithm (see, for example, Lessard and Pouliot 2001). Simulated annealing is a generic, global optimization algorithm that uses a probabilistic 'hill-climbing' approach to escape from local minima. Starting with an initial guess for the exposure times, a perturbation is applied to generate a new treatment plan ${{\boldsymbol{x}}}_{{\rm{new}}}.$ The new treatment plan is then accepted with a probability $P$ given by:

$$\begin{eqnarray}P=\left\{\begin{array}{cc}1 & {\rm{if}}\,{\rm{\Delta }}f\leqslant 0\\ \exp \left[-{\rm{\Delta }}f/{T}_{k}\right] & \mathrm{if}\,{\rm{\Delta }}f\gt 0\end{array}\right.,\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Delta }}f=f\left({{\boldsymbol{x}}}_{{\rm{new}}}\right)-f\left({{\boldsymbol{x}}}_{{\rm{old}}}\right).$ This allows treatment plans that increase the objective function (i.e. decrease the CN) to be occasionally accepted, thus enabling 'hill climbing' out of local minima. Accepted new treatment plans then become the 'old' treatment plan and the algorithm proceeds iteratively. The parameter ${T}_{k}$ is the so-called 'temperature' at iteration $k,$ and decreases according to the following linear 'cooling schedule':

$$\begin{eqnarray}&&{T}_{k+1}=\,\alpha {T}_{k},\end{eqnarray} \tag{ 6 }$$

where $\alpha \lt 1$ is the cooling rate. As the temperature decreases, fewer solutions that increase the objective function are accepted, according to equation (5). Typically, the temperature is not changed until a desired number of candidate treatment plans have been accepted. Because simulated annealing is a probabilistic technique, it can be rerun several times and always arrive at slightly different approximations of the global optimum.

For inverse treatment planning, the form of the perturbation was chosen randomly selecting one of the operating points in the treatment plan and adding a normally distributed random variable $n$ with a standard deviation $\sigma$ to it. That is,

$$\begin{eqnarray}{t}_{i}\to \left\{\begin{array}{cc}{t}_{i}+n & {\rm{if}}\,{t}_{i}+n\gt 0\\ 0 & {\rm{otherwise}}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

where ${t}_{i}$ is a random chosen element of ${\boldsymbol{t}}$ and $n\sim N\left(0,{\sigma }^{2}\right).$ The second case enforces the positivity constraint in equation (1).

Convergence is quantified by calculating the standard deviation of the CNs of the last accepted treatment plans, and if the value is below a threshold then convergence is declared. The parameters of the simulated annealing algorithm implemented for the inverse RTP can be found in table 1. These parameters were chosen by manual trail-and-error to optimize the average CN of the datasets described in section 3. A production implementation of such an RTP algorithm could be tuned based on real clinical PTVs.

Table 1. List of the parameters of the Simulated Annealing algorithm.

Parameter	Value
Initial temperature	0.01
Final temperature	1E-8
Rate at which temperature is decreased	0.8
The standard deviation of the normal distribution used to randomly perturb the exposure time of an operating point of the treatment plan	5 s
Maximum number of accepted treatment plans allowed before the temperature is decreased	20
Maximum number of perturbations allowed before the temperature is decreased	100
Minimum number of temperatures before the algorithm is allowed to exit	5
Threshold of the standard deviation of the conformation number that is used to test for convergence	1E-3
Number of conformation numbers stored to test for convergence	200

Simulated annealing has several benefits that make it ideally suited for inverse RTP algorithms. First, because it is a global optimization algorithm, it is able to avoid getting stuck in the local minima of equation (1). Second, the positivity constraints in equation (1) are trivially implemented by rejecting any perturbations in the current state that result in negative exposure times (see equation (7)). This algorithm has been used extensively in other inverse RTP implementations (Webb 1991, Oldham and Webb 1995, Pouliot et al 1996, Cho et al 1998, Pouliot et al 1999, Lessard and Pouliot 2001, Ayotte et al 2008, Kubicky et al 2008, Beliën et al 2009, Diest and Gorissen 2016).

All treatment plans were evaluated using several metrics, which are listed in table 2.

Table 2. List of the metrics used to quantify the treatment plan calculated by the inverse radiation treatment planning algorithm.

Metric	Notes
Total treatment time	The treatment times reported are the times for which the radiation source is inside the patient, delivering the radiation dose. The additional time required for treatment planning, robot movement, or insertion of the x-ray source are not included. Note that when the treatment is actually delivered, ordering the operating points to minimize the number of deflections, translations, and kVp changes might reduce the total operating room time required
Conformation number	Discussed in section 2.1.2.1
Gamma index analysis	The gamma index $\gamma$ was introduced by Low et al (1998) and is a clinically accepted metric of radiation treatment planning (Nelms and Simon 2007). It quantifies the agreement between the target isodose surface and the dose distribution of the treatment plan relative to dose tolerance (X) and spatial tolerance (Y) criteria that are specified as X%/Y mm (e.g. 5%/5 mm). A value $\gamma$ is assigned to each voxel in the target isodose surface, with $\gamma \leqslant 1$ if there is any voxel in the treatment dose map that is within Y mm that has a dose that differs by less than X%. If the acceptance criteria are not met, then $\gamma \gt 1.$ The fraction of voxels on the target isodose surface meeting the acceptance criteria (called the 'gamma passing rate' and ranging between 0 and 1) is typically reported. In calculating γ, local normalization was used for the dose-difference and a low-dose threshold of 10% was used, as recommended by AAPM Task Group 119 (Ezzell et al 2009)
Optimization time	This is the calculation time taken by the inverse RTP to create the optimized treatment plan. It does not include the time taken to load the dose distributions of each operating point into memory, since this will only be done once upon initialization of the RTP system. Shortening the optimization time improves the clinical workflow and reduces the treatment costs

Because of the choice of the perturbation $n$ used for the simulation annealing (equation (7)), calculating the dose distribution of the perturbed treatment plan is relatively fast. To see this, let ${{\boldsymbol{t}}}^{{\prime} }$ be a new treatment plan that is obtained by perturbing a treatment plan ${\boldsymbol{t}}.$ If ${\boldsymbol{D}}$ is the dose distribution corresponding to ${\boldsymbol{t}},$ then the dose ${D}_{i}$ in the ${i}^{th}$ voxel of ${\boldsymbol{D}}$ is given by:

$$\begin{eqnarray}&&{D}_{i}=\displaystyle \sum _{j=1}^{N}{d}_{ij}{t}_{j},\end{eqnarray} \tag{ 8 }$$

where $N$ is the number of operating points (i.e. the number of elements in ${\boldsymbol{t}}$), ${t}_{j}$ is the ${j{\rm{t}}{\rm{h}}}^{}$ element of ${\boldsymbol{t}}$ (i.e. the exposure time of the $j{\rm{t}}{\rm{h}}$ operating point), and ${d}_{ij}$ is the simulated dose rate in the $i{\rm{t}}{\rm{h}}$ voxel due to the ${j{\rm{t}}{\rm{h}}}^{}$ operating point. The perturbation can be written as ${\boldsymbol{t}}^{\prime} ={\boldsymbol{t}}+{\rm{\Delta }}{\boldsymbol{t}},$ where the elements of ${\rm{\Delta }}{\boldsymbol{t}}$ are all zero except for one nonzero element. Then, the new dose distribution ${{\boldsymbol{D}}}^{{\prime} }$ can be calculated according to:

$$\begin{eqnarray}\begin{array}{lll}{D}_{i}^{{\prime} } & = & \displaystyle \sum _{j=1}^{N}{d}_{ij}{t}_{j}^{{\prime} }\\ & = & \displaystyle \sum _{j=1}^{N}{d}_{ij}\left({t}_{j}+{\rm{\Delta }}{t}_{j}\right)\\ & = & \displaystyle \sum _{j=1}^{N}{d}_{ij}{t}_{j}+\displaystyle \sum _{j=1}^{N}{d}_{ij}{\rm{\Delta }}{t}_{j}\\ & = & {D}_{i}+\displaystyle \sum _{j=1}^{N}{d}_{ij}{\rm{\Delta }}{t}_{j}.\end{array}\end{eqnarray} \tag{ 9 }$$

That is, only the dose distribution of a single operating point (corresponding to the nonzero element of ${\rm{\Delta }}{\boldsymbol{t}}$) needs to be added to the dose distribution of the treatment plan. Because of this, updating the treatment dose distribution at each iteration of the simulated annealing algorithm is computationally inexpensive.

Because of its stochastic nature, simulated annealing is relatively insensitive to the quality of the starting guess. As such, the starting guess was chosen to be an exposure time of 1 s at each of the operating points.

The RTP algorithm was implemented in MATLAB^®, Release 2017a (The MathWorks, Inc., Natick, Massachusetts) and all tests were performed on a PC with a 2.8 GHz Intel i7 CPU and 16 GB of RAM. The dose maps for each operating point are loaded into memory once during the initialization of the RTP. While this makes the subsequent calculation time relatively quick, it comes at the cost of significant memory overhead. For example, for typical parameters (i.e. 301 × 301 × 301 volumes storing the dose rate in 64 bit floating point precision) a single dose map requires ∼218 MB. Thus storing the dose maps for all 54 unique operating points needs about 11.8 GB of RAM.

Three categories of datasets were used to characterize the performance of the inverse RTP. The first included simple shapes, namely a sphere and a prolate ellipsoid. The second category was composed of four treatment volumes that were calculated by combining an arbitrary set of exposure times from the device's available operating points (these datasets are referred to as Combo 1–4). These datasets test the algorithm's ability to reproduce known input exposure times. Finally, the third category included four datasets from The Radiotherapy Optimization Test Set (TROTS) (Breedveld and Heijmen 2017). TROTS provides voxelized PTVs that were obtained from treatment planning of radiation therapy patients. The four datasets were chosen from radiation therapy patients with liver tumors. These 10 datasets were chosen to cover a wide range of PTV shapes, sizes, and symmetries. The PTVs of the 10 datasets are visualized in the appendix. The inverse treatment planning was performed with the same parameters for all datasets, as listed in tables 1 and 3. Only a single kVp was chosen because it was found (Badali et al 2019) that shape of the dose distribution at each operating point did not depend significantly on the kVp. Thus choosing 100 kVp (which produces the highest dose rate at a fixed power of all available kVps) results in the shortest treatment times.

As discussed in section 2.1.2.2, the regularizer penalizes treatment plans with long treatment times. The choice of the regularizer strength μ is crucial to balance the quality metrics (CN and gamma passing rate) and the treatment time.

Figure 5 shows the dependence of these parameters on the regularizer strength, which was scanned from 0 to 1. A value of $\mu =0$ corresponds to an objective that only minimizes the difference between the PTV surface and the treatment plan, and thus gives the best agreement to the target isodose surface, with a mean CN of 0.81 and a mean gamma passing rate of 0.93. We see that if μ = 1E-3 is used, there is only a minor reduction in the mean CN and gamma passing rate (2.5% and 3%, respectively), whereas the mean treatment time is reduced by a full minute (a 20% reduction). This indicates that regularization has the desired effect of significantly reducing the treatment time without compromising the quality of the treatment plan. Values of μ above about 1E-2 are too heavily biased towards satisfying the regularizer and force all exposure times to be exactly zero. We note that choosing the regularizer strength in this way produces the optimal treatment plan on average, although it might not be the best regularization for each individual dataset. It is likely that the medical physicist will want to tune the regularization strength μ for a given PTV in order to appropriately balance between the treatment time and plan quality.

Zoom In Zoom Out Reset image size

The optimization time scales with the total number of voxels used (see figure 6). Since we found that using 256³ voxels requires roughly 1 min of computation time to perform the optimization, we would expect ∼8 min for 512³ voxels and an hour for 1024³ voxels. In addition to the undesirably long computation times required, large volumes would place substantial demands on other system resources such as memory.

Zoom In Zoom Out Reset image size

Because it is expected that the PTVs do not have sharp features, we expect that little information will be lost if the PTV volume is blurred by downsampling prior to performing the optimization. This would significantly improve the calculation time (and memory requirements) of the inverse RTP algorithm. We investigated this approach by using a 320 × 320 × 320 volume with 0.4 mm voxels as the 'baseline,' and then proceeded to downsample the volume to 256³, 128³, and 64³ voxels, with voxel sizes of 0.5 mm, 1 mm, and 2 mm, respectively. These parameters were chosen to maintain the same volume dimensions (e.g. 320 × 0.4 mm = 64 × 2 mm). This allows us to solely investigate the impact of blurring the dose distribution on the treatment plan quality, without having to consider the effects of cropping the dose distribution.

The results are shown in figure 6. As expected, the optimization time is dramatically reduced as the number of voxels is decreased, starting from over 2 min for 320³ and decreasing to less than 1 s for 64³. While achieving this 150× speedup, the quality of the treatment plan (as indicated by the CN and the gamma passing rate) is unaffected. That is, a significant speedup of the optimization time can be achieved without compromising the quality of the treatment plan by downsampling the volume prior to optimization.