A fast robust optimizer for intensity modulated proton therapy using GPU

Abstract Robust optimization has been shown to be effective for stabilizing treatment planning in intensity modulated proton therapy (IMPT), but existing algorithms for the optimization process is time‐consuming. This paper describes a fast robust optimization tool that takes advantage of the GPU parallel computing technologies. The new robust optimization model is based on nine boundary dose distributions — two for ±range uncertainties, six for ±set‐up uncertainties along anteroposterior (A‐P), lateral (R‐L) and superior‐inferior (S‐I) directions, and one for nominal situation. The nine boundary influence matrices were calculated using an in‐house finite size pencil beam dose engine, while the conjugate gradient method was applied to minimize the objective function. The proton dose calculation algorithm and the conjugate gradient method were tuned for heterogeneous platforms involving the CPU host and GPU device. Three clinical cases — one head and neck cancer case, one lung cancer case, and one prostate cancer case — were investigated to demonstrate the clinical feasibility of the proposed robust optimizer. Compared with results from Varian Eclipse (version 13.3), the proposed method is found to be conducive to robust treatment planning that is less sensitive to range and setup uncertainties. The three tested cases show that targets can achieve high dose uniformity while organs at risks (OARs) are in better protection against setup and range errors. Based on the CPU + GPU heterogeneous platform, the execution times of the head and neck cancer case and the prostate cancer case are much less than half of Eclipse, while the run time of the lung cancer case is similar to that of Eclipse. The fast robust optimizer developed in this study can improve the reliability of traditional proton treatment planning in a much faster speed, thus making it possible for clinical utility.

conversions from computed tomography (CT) numbers to stopping power values, inhomogeneity, setup errors, anatomical changes, and so on. [3][4][5][6][7][8] As a result, the actual delivered dose distribution of IMPT plan will be affected by these uncertainties. Scanning spots or scanning lines are preset and optimized to meet the dose objectives and constraints in the inverse treatment plan. However, due to the steep falloff of Bragg peaks, IMPT is extremely sensitive to proton range and patient setup uncertainties. For photon therapy, clinical target volume (CTV) is expanded to the planning target volume (PTV) by a given margin to compensate for plan uncertainties, but the PTV expanding technique is not as effective in IMPT. 9 The solution of the above problems is to consider these uncertainties during the inverse treatment plan optimizationa process known as the robust optimization. 10 Several related reports for robust optimization of IMPT have shown better results than common PTV-based IMPT optimization techniques. Unkelbach et al. 6 proposed a robust linear programing method that was found to yield treatment plans less sensitive to range variations. A worst case optimization method was developed by Pflugfelder et al. 10 to account for uncertainties based on several possible realizations of the uncertainties A minimax robust optimization method was utilized by Fredriksson et al. 11 to minimize the penalty of the worst scenario which yielded better results for the tested clinical cases. Generally, in worst-case-based robust algorithms, the nominal dose distribution and perturbed dose distributions are calculated. When computing the target function, upper bound constraints are applied to the maximum dose while lower bound constraints are applied to the minimum dose. [12][13][14][15][16][17] A selective robust optimization method was derived by Li et al. 18 from worst case optimization methods and its objective function was selectively computed from either the worst-case dose or the nominal dose. Li et al. 19 reported that robust optimization in IMPT of lung cancer can reduce the dose variation caused by setup uncertainty and anatomical changes during treatment compared with PTV-based planning. Although three worst-case-based robust methods (composite worst case, voxelwise worst case and objectivewise worst case) have different behaviors, and no particular method was superior to the others under all circumstances. 20 Compared with minmax robust optimization approaches, worst-case dose approaches were less sensitive to uncertainties for the prostate and skull base cancer patients, whereas the minmax approach was superior for the head and neck cancer patients. 21 A 4D robust optimization was developed by Liu et al. 22 and it produced more robust plans for targets and normal tissues, compared to 3D robust optimization. In order to accelerate the robust optimization, a constraint generation solution method was developed by Mahmoudzadeh et al, 23 which reduced the optimization time to about 12 min. A chance-constrained optimization method was proposed in IMPT planning to hedge against the influence of uncertainties, 24 but several hours were needed to finish the optimization. The chance-constrained optimization method explicitly controlled the tradeoff between plan quality and plan robustness, and it used linear programing without parallelization, making the method very slow. Jiasen Ma et al developed an all scenario and MC-based IMPT optimizer and employed GPU to reduce the computational time. 25 In this research, we propose the uncertainty model which contains nine boundary dose distributions corresponding to different setups and range perturbations, two for ±range uncertainty, six for ±set-up uncertainties along anteroposterior (A-P), lateral (R-L) and superior-inferior (S-I) directions, and one for the nominal situation.
All the nine dose distributions consider the target function and optimized using the same optimization objective in each optimization iteration until all the nine dose distributions approach the dose constraints as much as possible. In order to reduce the optimization time, the proposed method is implemented in a CPU-GPU parallel platform. Three clinical cases are used to test the robust optimization effectiveness.

| The uncertainty optimization model
For the traditional PTV-based IMPT optimization model, the standard quadratic objective function can be used to represent optimization objectives and constraints similar to IMRT. 26 Mathematically speaking, treatment plan optimization is the minimization of the objective function and is given below by: where i is the index of voxel i, j is the proton beamlet j and m is the For dose volume constraints, H is determined by sorting the doses, finding violating voxels. For example, D q <pcGy means no more than q% of the ROI may receive a dose greater than p. The voxels of the ROI are sorted in ascending order of dose received and only those voxels which cause q to be exceeded are included in the objective function (H = 1). 27 As mentioned earlier, φ i;j is perturbed by range and setup uncertainties impacting on the finally dose contributions of the PTV and OARs. For the IMPT plan delivery, the final dose distribution is affected by many factors including complex inhomogeneity, setup errors, anatomical changes, and so on. [3][4][5][6][7][8] All these interference factors will result finally in two kinds of deviations relative to the primary plan, proton range uncertainties, and patient setup uncertainties. Compared to the traditional plan optimization, robust plan optimization takes the two kinds of deviations into consideration beforehand. The challenge is that, for IMPT robust optimization, the computational complexity increases proportionally to the number of potential of considered uncertainties. In theory, more uncertainties can lead to better robust plan during robust optimization. Considering the large computation cost even for one possible situation, however, we have to balance between the uncertainties taken into consideration and the actual performance of robust optimization. As reported previously 12,13 , the simplification for IMPT robust optimization is usually handled in two processes. The first simplification process is discretizing the proton range uncertainties and patient set-up uncertainties before selecting the extreme uncertainties as a boundary. In this study, we use the same discretizing method to simplify the robust optimization task, as to be shown next. The second simplification is that only the worst dose distributions are included in the target function. Thus, we can only get a relatively conservative result due to the robust optimization simplification. In this research, we improved the target function by considering all distributions in the total target function so that not merely the worst distribution but all the dose distributions approached the dose constrains. The CPU-GPU parallel platform was further adopted to increase computational efficiency, as to be illustrated next.
The proposed robust optimization model adopted in this research is shown in Eq. (4) in which nine boundary dose distributions are considered, two for ±range uncertainties, six for ±set-up uncertainties along anteroposterior (A-P), lateral (R-L) and superior-inferior (S-I) directions, and one for nominal situation.
where R is the total nine boundary dose distributions considering the range and setup uncertainties, and k represents one possible boundary dose distributions of R. ρ k CTV and ρ k OARs are the weights of CTV and OARs in the kth uncertainty situation, usually equal to 1.

D k
i is the dose of voxel i in the kth uncertainty, and φ k i;j is the dose contribution of beamlet j to voxel i in the kth uncertainty situation.
The objective function is minimized by optimizing ω j . It is impossible to forecast which uncertainty will occur in the treatment, so it is necessary to ensure that the treatment planning is optimal on the premise of the existence of the nine boundary conditions.
As the same with the reported researches, 12 we made the assumptions: (1) all the beamlets were affected by the range uncertainties and changed synchronously; (2) the max range uncertainties and setup uncertainties value kept the same during the treatment course.

| The robust optimization method
The CG method is one of the most popular techniques for solving large scale unconstrained optimization problems owing to small memory footprint and simplicity. Figure 1 demonstrates the flow chart of the CG method. x 0 is the initial solution, and k is the current iteration number. g k is the gradient of x k , ɛ is the tolerance to stop the iterative process and N is the maximum number of iterations. d is the searching direction, and λ is the step length.
In the first iteration, the searching direction d 0 is equal to the negative gradient, that is rf x 0 ð Þ. For other iterations, d k is defined by.
where β k is a scalar. It guarantees that all elements of x kþ1 are greater than 0.
In this paper, the DYHS method [Eq. (8)] is adopted to calculate the CG parameter β.
The DYHS method is a combination of the DY method and the HS method. [28][29][30] The line search subroutine yields a step length λ k that satisfies the standard Wolfe conditions [Eqs. (11) and (12)], where g k ¼ rf x k ð Þ, δ ¼ 10 À4 , and θ ¼ 0:9. The first inequation is to guarantee that the decrease of the target function is proportional to the tangential decline at least. The second inequation ensures the slope at λ k is not strongly negative.

| Heterogeneous platform with CPU and GPU
A GPU dose engine is adopted to calculate the dose contribution matrix φ i;j and then the matrix is converted to the most memory efficient sparse matrix format, that is compressed sparse row (CSR) XU ET AL.
| 125 format. 31 The CSR format uses three arrays to store the nonzero     Figure 3 shows the DVH bands of the CTV, brain stem, and spinal cord for the head and neck cancer case;