Background: Volumetric modulated arc therapy (VMAT) is widely used in clinical practice. It not only significantly reduces treatment time, but also produces high-quality treatment plans. Current optimization approaches heavily rely on stochastic algorithms which are time-consuming and less repeatable. In this study, a novel approach is proposed to provide a high-efficient optimization algorithm for VMAT treatment planning. Methods: A progressive sampling strategy is employed for beam arrangement of VMAT planning. The initial beams with equal-space are added to the plan in a coarse sampling resolution. Fluence-map optimization and leaf- sequencing are performed for these beams. Then, the coefficients of fluence-maps optimization algorithm are adjusted according to the known fluence maps of these beams. In the next round the sampling resolution is doubled and more beams are added. This process continues until the total number of beams arrived. The performance of VMAT optimization algorithm was evaluated using three clinical cases and compared to those of a commercial planning system. Results: The dosimetric quality of VMAT plans is equal to or better than the corresponding IMRT plans for three clinical cases. The maximum dose to critical organs is reduced considerably for VMAT plans comparing to those of IMRT plans, especially in the head and neck case. The total number of segments and monitor units are reduced for VMAT plans. For three clinical cases, VMAT optimization takes < 5 min accomplished using proposed approach and is 3–4 times less than that of the commercial system. Conclusions: The proposed VMAT optimization algorithm is able to produce high-quality VMAT plans efficiently and consistently. It presents a new way to accelerate current optimization process of VMAT planning. Keywords: Volumetric modulated arc therapy, Fluence-map optimization, Leaf-sequencing Background from the idea of delivering plans with a large number of VMAT is widely used in cancer treatment in radiation on- gantry positions. The fluence map of the beam is cology departments due to its high efficiency in treatment pre-calculated and decomposed to several apertures. delivery . Compared to conventional IMRT, VMAT de- These apertures is then delivered at a given gantry pos- livers radiation dose while MLC leaf, dose rate and gantry ition by multiple arcs. IMAT requires more arcs which move simultaneously . It uses less treatment time and causes extended treatment time. Direct aperture total monitor units (MU) compared to conventional optimization (DAO) was thus proposed to handle the IMRT technique . The predecessor of modern VMAT complexity of VMAT optimization using stochastic ap- is intensity modulated arc therapy (IMAT) which was first proach. The stochastic approach is computationally inten- developed by Cedric Yu in 1995 [4, 5]. It was motivated sive and time-consuming [6–8]. Later, Otto presented an iterative algorithm for VMAT optimization [9, 10]. This algorithm employs progressive sampling strategy and * Correspondence: firstname.lastname@example.org; email@example.com aperture-based algorithm for VMAT optimization where Department of Radiation Oncology, National Cancer Center/Cancer Hospital, high-quality dose plan can be achieved by a single arc. Chinese Academy of Medical Sciences, Peking Union Medical College, Beijing 100021, China © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Yan et al. Radiation Oncology (2018) 13:101 Page 2 of 13 This technique is successfully adopted in commercial optimization a coarse sampling of the gantry positions is treatment planning system and has been applied to a wide used to model the gantry rotation range. The initial 5 range of treatment sites . beams with equal-space are chosen for the first Compared to IMRT, VMAT presents a complex optimization stage. The MLC positions and MUs for the optimization problem because of the significantly increased initial 5 beams are achieved by the fluence-map number of plan parameters such as gantry angle, MLC leaf optimization and leaf-sequencing algorithms. Next, the position, dose rate, etc. [12, 13]. Various approaches were sampling resolution is doubled. And another 5 beams proposed attempting to solve this problem. Most of them are added to the plan and optimized. The sampling reso- are evolved from approaches originally developed for IMRT lution is continuously increased and more beams are [14–22]. Currently, most of commercial treatment planning added to the plan until the total number of beams is ar- systems provide VMAT optimization function and heavily rived. As shown in Fig. 1, the number of new beams in 5 rely on DAO algorithm. RapidArc (Varian Medical System, optimization stages is 5, 10, 20, 40, and 60, while the Palo Alto, CA, USA) employs a progressive sampling strat- corresponding beam spacing is 72°, 36°, 18°, 9°, and 6°. egy and simulated annealing-based DAO algorithm for A VMAT optimization approach was developed based on VMAT planning . Due to the nature of stochastic ap- the existing fluence-map optimization and leaf-sequencing proach, the optimization process is time-consuming and algorithms provided by an in-house developed treatment the optimization result is seldom repeatable. As reported planning system . The flowchart of VMAT planning the maximal DVH variation could be up to 2% for OARs proposed in this study is presented in Fig. 2.First, the basic using RapidArc . SmartArc (Philips Healthcare, Inc., plan parameters such as couch angle and arc range are set Thornton, CO, USA) utilizes an IMRT plan (consisting of by planner. An initial plan consisting of few beams with equal-space fields) to initialize the arc plan. The fluence equal-space is created and their fluence maps are generated maps of IMRT plan are then segmented and redistributed by FMO algorithm. These resulting fluence maps are then to their neighboring angles. A local gradient-based processed by the leaf sequencing algorithm, and the optimal optimization approach is used to fine-tune these apertures MLC positions and MUs are determined. Next, the coeffi- to meet the mechanical constraints for actual delivery [16, cients of FMO algorithm are adjusted based on plan object- 17]. In this approach, The fluence maps of all beams are ac- ive and known fluence maps of these beams. In the next tually derived from few static beams and not fully deter- round, the sampling resolution is increased and more mined by fluence map optimization (FMO) algorithm. beams are added to the plan. The optimization process The purpose of this study is to develop a high-efficient continues until the maximal number of beams is arrived. optimization approach for VMAT planning. The remain- der of this paper is organized as follow. The progressive Fluence-map optimization sampling strategy is first introduced and the optimization A gradient-based optimization approach, fast monotonic process is described. Then, two important components, descent (FMD) algorithm is employed for VMAT fluence-map optimization algorithm and leaf-sequencing optimization. Due to the nature of gradient-based algo- algorithm, are explained in brief. These algorithms were rithm, the global minimum of objective function can be developed on an in-house developed treatment planning found quickly. The goal is to find minimum of objective system. Three typical clinical cases (head and neck, lung, function subject to non-negative fluences for all pencil prostate) were evaluated on both in-house developed and beams. A non-synchronous updating scheme is used commercial treatment planning systems. The dosimetric which allows only one pencil beam adjusted at a time. quality, delivery efficiency, and running time of VMAT The detail of classic FMD algorithm is described in ref- plans were then assessed. Finally, the advantage and disad- erences [23–25]. To accommodate the progressive sam- vantage of this approach are discussed. pling strategy, FMD algorithm was modified to account for the varied number of beams in multiple optimization Methods stages. The formulation of the modified FMD algorithm Progressive sampling strategy is described in Appendix 1. In brief, new beams are Most of VMAT optimization algorithms model Linac added to the plan in each optimization stage and their source as a series of static gantry positions. The MLC fluence maps are optimized with the fluence maps of old new positions and MU setting is then determined at each beams known. Here, x denotes the fluence map of old gantry position. Provided the complexity of VMAT new beam which is unsolved in current stage and x optimization, to limit the scale of problem a progressive denotes the fluence map of old beam which is solved in sampling strategy is employed in this study. It starts with previous stage. For a given voxel v , its dose is calcu- ijk new old a coarse sampling of gantry positions, then move to fine lated by the sum of dose contributed from x and x . new sampling of gantry positions. An example is demon- For solving x the coefficients of FMD algorithm is ad- strated as shown in Fig. 1. At the beginning of VMAT justed based the planning dose of target and the dose Yan et al. Radiation Oncology (2018) 13:101 Page 3 of 13 Fig. 1 The illustration of progressive sampling scheme. a The first beam set with 5 control points (beam spacing 72°). b The second beam set with additional 5 control points (beam spacing 36°). c The third beam set with additional 10 control points (beam spacing 18°). d The fourth beam set with additional 20 control points (beam spacing 9°). e The fifth beam set with additional 20 control points (beam spacing 6°) old new contribution from x . As fluence map of x known, it decompose fluence map into several apertures with uni- is processed by a leaf sequencing algorithm and a uni- form fluences. Given the dynamic feature of beams in form fluence map for a single aperture is obtained. The VMAT, the leaf sequencing algorithm has to deal with uniform fluence map replaces the original fluence map more mechanical constraints related to gantry speed, leaf new new of x and is used as the final value of x . speed, and dose rate than those of static beams. Mixed integer linear programming (MILP) was successfully ap- Leaf sequencing plied in many industrial applications due to its capability Several leaf sequencing algorithms were developed in in handling large-scale optimization problems [30, 31]. the past years such as powers of 2, linear programming, In this study, a MILP model was developed using a com- and graphics searching [26–29]. These algorithms mercial optimization software package (IBM ILOG Fig. 2 The flowchart of the proposed optimization algorithm for VMAT treatment planning Yan et al. Radiation Oncology (2018) 13:101 Page 4 of 13 CPLEX). It assumes that the fluence map of a beam at a CI = V /V , where V is the target volume and TV 95% TV given control point is a uniform fluence map for a single V is the volume corresponding to 95% dose prescrip- 95% aperture. The aperture of the new beam is determined tion. If V =V , CI =1 which is perfect. 95% TV subjected to the apertures of the two old beams closet to HI = D /D , where D and D are doses received 5% 95% 5% 95% it (one is at front side and one is at end side). As the ori- at least 5 and 95% target volume respectively. if D = 5% new ginal flence map of x is non-uniform, it is necessary D , HI =1 which is perfect. 95% to convert it to a uniform fluence by a leaf-sequencing Quality score (QS) quantitatively measures the differ- algorithm with certain constraints. These constraints are ence between the dose plan and the dose constraint of ðM −C Þ mostly related to mechanical limitation of Linac and j j all anatomical structures. QS ¼ j j if object- MLC and described in Appendix 2. Due to the nature of j progressive sampling strategy, the beams solved in the ive is violated by plan dose, where C is the j-th earlier stages have larger apertures and weights than the dose-volume constraint, M is the corresponding plan beams solved in the latter stages. This is because the dose. If no objective is violated, QS = 0 which is perfect. number of beams is less in the earlier stages and more Weighted root mean-square error (WE) measures the in the latter stages. The MILP model is described in difference between the plan dose and prescribed dose of Appendix 2. The goal of this model is to search for an all voxels and represents the value of objective function. X X X 1 2 2 2 optimal solution, uniform fluence map, from a w ðP −D Þ þ w ðP −D Þ þ w ðP −D Þ 2 ijk ijk ijk ijk ijk ijk ijk ijk ijk i; j;k∈V TV i; j;k∈V CO i; j;k∈V NT WE ¼ ð Þ , non-uniform fluence map under certain constraints. The where P , D , and W are defined in Eq. 1, and N is MILP problem is solved using branch-and-bound tech- ijk ijk ijk the total number of voxels belonging to TV, CO and nique provided by CPLEX software. As an example of NT. For a plan which P = D for all voxels, WE = 0 fluence map consisting of 30×30 bixels, the total number ijk ijk which is perfect. of constraints is ~ 1800 and the total number of inde- The dose distribution of plans is evaluated based on pendent variable is ~ 2000. The processing time is ap- cross-sectional dose distribution and dose-volume histo- proximately 0.5 s per beam. gram. The total number of segments and monitor unit for both plans are recorded. Additionally, estimated arc deliv- Plan evaluation ery times and the computation time of plan optimization Three clinical cases were selected for evaluation in- are recorded. For demonstration purpose, few important cluding head and neck, prostate, and lung. The pa- ROIs are used in plan optimization. For head and neck tient CT images and contours were exported from case, they are PTV, spinal cord, left and right parotids and Pinnacle (Philips Healthcare, Inc., Thornton, CO, mandible. For prostate case, they are PTV, left and right USA). Then they are segmented into regions of inter- femur heads, bladder and rectum. For lung case, they are est(ROIs)for plan optimization usinganin-house PTV, left and right lungs, spinal cord and heart. All tests developed treatment planning system. In this study, are performed on a DELL Optiplex N9010 computer, the simulated leaf width is 0.5 cm and the leaf step is equipped with Intel(R) i7–3770 CPU and 72GB RAM. 0.5 cm, which resulted in beamlet size of 0.5×0.5 cm. Corresponding to the VMAT plans made on in-house de- The range of leaf speed is 0.5–1cmper degreebe- veloped system, three VMAT plans made on Pinnacle plan- cause the leaf speed is 3–6 cm per second and gantry ning system were assessed. They were made by experienced speed is 6 degree per second. The prescription dose planners and approved for clinical treatment. These plans for PTV and dose constraints for OARs are the same are VMAT plans consisting of 2 full arcs with 4° spacing for both IMRT and VMAT plans for each case. A sin- and 6 MV beams computed with convolution superposition gle dynamic arc plan consisting of 180 beams and 2° algorithm. The dosimetric and delivery statistics of these spacing was created and fluence maps were optimized clinically approved plans are compared to those of IMRT by the proposed VMAT optimization algorithm on and VMAT plans made on the in-house developed system. the in-house developed system. The numbers of Note that the optimization algorithm, the dose calculation beams in five optimization stages are 15, 30, 60, 120, engine, the dose-volume constraints and many others of and 180, while the corresponding beam spacing is clinically approved plans are different from those of plans 24°, 12°, 6°, 3°, and 2°. A corresponding IMRT plan made on the in-house developed system. The metric, WS, is consisting of 9 equal-space fields was created and flu- not calculated for the clinically approved VMAT plans be- ence maps were optimized by FMD algorithm on the cause the weight specification in Pinnacle system is differ- in-house developed system. The dosimetric quality of ent. For distinguishing clinically approved VMAT plans plans is quantified by conformity index (CI), homo- from theplans madeon thein-housedeveloped system, geneity index (HI), quality score (QS) , and VMAT plans made on pinnacle system are represented by weighted root mean-square error , which are de- VMAT_P. fined as below. Yan et al. Radiation Oncology (2018) 13:101 Page 5 of 13 Results brainstem is reduced by 10 Gy. Dose distribution of The fluence map and plan dose distribution achieved by VMAT plan is more uniform with fewer hot and cold progressive sampling strategy is demonstrated in Fig. 3. spots. The QS, WE, CI, and HI for both plans are shown For each plan, the subplots at the top are the segments in Table 1. The WE is reduced by 10% for VMAT plan, of all beams while the subplots at the bottom are the while QS, CI and HI are similar for both plans. The total dose distribution in axial view. The grey level of segment number of segments and monitor units for both plans indicates the magnitude of its weight. White represents are shown in Table 1. The number of segment is similar the highest value (255) and black represents the lowest for both plans. The total MU for VMAT plan is reduced value (0). The corresponding beam arrangements are by 40%. The computation time is 40 s for IMRT plan shown in Fig. 1. In the earlier stages, the values of flu- and 312 s for VMAT plan as shown in Table 1. The esti- ence maps of new beams are larger as there are fewer mated delivery time is between 96 and 192 s based on a beams. As more beams included, the values of fluence Varian Trilogy machine with variable dose rate (300–600 map of new beams are smaller. The plan dose quality MUs per minute). VMAT_P plan show better plan qual- improves quickly in the earlier stages (as shown in Fig. ity and fewer MUs than VMAT plan, but the number of 3b and c), and gradually saturates in latter stages (as segments and optimization time are higher. shown in Fig. 3d and e. For the prostate case, the objective is to maintain uni- For the head and neck case, the objective is to main- form dose (76 Gy) to the PTV while keep the rectal dose tain uniform dose (60 Gy) to the PTV, while minimize at V50 Gy < 50% and maximize sparing of bladder the dose to the parotid gland (V30Gy < 50%) and (V50 Gy < 50%) and femoral head (V50 Gy < 5%). Dose maximize the sparing of the brainstem (Dmax< 54 Gy), distributions for IMRT and VMAT plans are shown in the larynx (Dmax<40Gy), and the mandible (Dmax< Fig. 6. DVHs for IMRT and VMAT plans are compared as 60 Gy) as much as possible. Dose distributions for IMRT shown in Fig. 7. The VMAT plan is comparable to the cor- and VMAT plans are shown in Fig. 4. DVHs for IMRT responding IMRT plan. The coverage of PTV is similar and VMAT plans are compared as shown in Fig. 5. The between VMAT and IMRT plans. For VMAT plan, OAR VMAT plan is better than the corresponding IMRT sparing is better for rectum, femoral heads (left and right) plan. The coverage of PTV is similar for both VMAT while similar for bladder. Dose distribution of VMAT plan and IMRT plans. The OAR sparing from VMAT plan is is more uniform and has less high dose region in normal better for brainstem and mandible, but is similar for lar- tissue. The QS, WE, CI, and HI for both plans are shown ynx and parotids (left and right). The max dose to in Table 1. The WE is reduced by 12% for VMAT plan, Fig. 3 The segments and dose distributions of VMAT plans for five beam sets with the increasing numbers of beams. a VMAT plan with 5 equal-space beams (72° spacing). b VMAT plan with 10 equal-space beams (36° spacing). c VMAT plan with 20 equal-space beams (36° spacing). d VMAT plan with 40 equal-space beams (9° spacing). e VMAT plan with 60 equal-space beams (6° spacing) Yan et al. Radiation Oncology (2018) 13:101 Page 6 of 13 Fig. 4 The transaxial plan dose distribution of head-and-neck case in (a) IMRT plan and (b) VMAT plan while QS, CI, and HI are similar for both plans. The total 45 Gy), heart (V30Gy < 40%, V40Gy < 30%) and esophagus number of segments and monitor units for both plans are (V50Gy < 50%). Dose distributions for IMRT and VMAT shown in Table 1. The number of segment is similar for plans are showninFig. 8. DVHs for IMRT and VMAT both plans. The total MU for VMAT plan is reduced by plans are compared as shown in Fig. 9.The dose coverage 16%. The computation time is 36 s for IMRT plan and of PTV is similar in both plans. For VMAT plan, dose to 294 s for VMAT plan as shown in Table 1. The estimated heart, esophagus, and cord are slightly higher than those of delivery time is between 118 s and 236 s based on a Varian IMRT plan, while dose to lung is similar. The Dose distri- Trilogy machine with variable dose rate (300–600 MUs bution of the VMAT plan is more uniform and focused on per minute). VMAT_P plan show better plan quality and PTV, however the dose distribution of IMRT plan is spread fewer MUs than VMAT plan, but the number of segments along anterior-posterior direction. The QS, WE, CI, and HI and optimization time are higher. for both plans are shown in Table 1.The WE,QS, CI, and For the lung case, the objective is to maintain uniform HI are similar for both plans. The total number of segments dose (60 Gy) to the PTV while keep the lung dose at and monitor units for both plans are shown in Table 1.The V20 Gy < 25% and maximize the sparing of cord (Dmax< number of segment for VMAT plan is reduced by 16%. Fig. 5 The DVH of head-and-neck case in IMRT plan and VMAT plan Yan et al. Radiation Oncology (2018) 13:101 Page 7 of 13 Table 1 Comparison of performance between VMAT plan and IMRT plan Treatment Plan Plan performance site type QS WE CI HI Number of segments MU Time (s) HN IMRT 0.92 79.91 1.83 1.26 129 1775 40 VMAT 0.84 71.63 1.82 1.23 108 960 312 Prostate VMAT_P 0.13 N/A 1.54 1.02 364 810 1080 IMRT 0.32 63.44 1.43 1.23 112 1418 36 VMAT 0.31 55.43 1.41 1.22 116 1180 294 VMAT_P 0.13 N/A 1.21 1.03 182 1076 1210 IMRT 0.13 15.9 1.83 1.14 166 1178 32 Lung VMAT 0.13 15.8 1.84 1.12 138 1162 216 VMAT_P 0.11 N/A 1.63 1.03 108 626 760 The total MU is similar for both plans. The computation optimization algorithm, dose calculation engine, time is 32 s for IMRT plan and 216 s for VMAT plan as dose-volume constraints between pinnacle system and showninTable 1. The estimated VMAT delivery times are in-house developed system. However, the optimization between 116 s and 232 s based on a general linear acceler- and delivery time of VMAT plans achieved by our ap- ator with variable dose rate (300–600 MUs per minute). proach is less. It is promising to implement this ap- VMAT_P plan show better plan quality and few MUs than proach in commercial planning system to accelerate the VMAT plan, but the optimization time is higher. current process of VMAT optimization. The selection of beam orientations in initial stages of Discussions optimization would affect the final dose considerably. The results of clinical case study show that plan quality However it could be minimized with increasing number of VMAT plans is similar to or better than that of trad- of beams. It was reported that optimizing beam orienta- itional IMRT plans. The VMAT plans, in many cases, tions is most valuable for a small numbers of beams (≤5) are able to achieve better OAR sparing compared to and the gain diminishes rapidly for higher numbers of IMRT plans. VMAT plan has the potential to tailor high beams (≥15) . In several pioneering studies of VMAT volume low dose regions because of its flexibility to optimization it showed 10–15 beams with equal-space spread out the dose at many beam angles. This character would be a better choice in initial stage of optimization reduces the need of specific structures defined to remove [9, 10, 17]. We also tested our algorithm with the differ- hot spots, thus reduces overall planning effort. Plans ent initial beam orientations in the initial stage. It was consisting of more than 2 arcs are not reported in this found that for a higher number of beams (≥15) in the work, but it is a viable option in clinical practice. We initial stage the final value of objective function is less also tested optimizing all beams in multiple arcs which affected. In the proposed algorithm the number of has slightly improved dosimetric quality. These data are beams in the initial stage is set to 15, while this value is not included in this paper. The dose quality of clinically 10 in RapidArc of Eclipse (Varian Medical System, Palo approved plans made on pinnacle planning system is Alto, CA, USA) and 15 in SmartArc of Pinnacle (Philips better. This may be caused by the differences of Healthcare, Inc., Thornton, CO, USA). Fig. 6 The transaxial plan dose distribution of prostate case in (a) IMRT plan and (b) VMAT plan Yan et al. Radiation Oncology (2018) 13:101 Page 8 of 13 Fig. 7 The DVH of prostate case in IMRT plan and VMAT plan Plan quality can be further improved by increasing the The current objective function is quadratic function number of beams. However, as more beams included, the which can be solved efficiently by gradient decent method. computation time for plan optimization and dose calcula- However, for real clinical use it is desired to incorporate tion could be increased substantially. We speculate that only more complex objectives and constraints such as in certain cases, e.g., with complex target shape requiring dose-volume histogram (DVH) and general equivalent uni- substantial leaf motion, it would be beneficial to increase form dose (gEUD). For DVH constraints, the current algo- gantry position sampling. It is observed that the number of rithm is applicable with minor modification. Those voxels segments for VMAT plan is less than or similar to that of whose doses exceed the given threshold will be identified IMRT plan in three cases. The total MU for VMAT plan is when comparing actual DVH to expected DVH. The also less than or similar to that of IMRT plan in three cases. weights for those voxels will be adjusted automatically in Overall, the delivery time for VMAT plan is significantly re- order to minimize the DVH difference. For gEUD con- duced compared to IMRT plan. Plan optimization times straints, the form of objective function will be changed as vary case by case but all within 5 min. Compared with the described by Wu andcanbesolved bygradientdes- selected commercial treatment planning systems, the run- cent method which is similar to FMD algorithm employed ning time of plan optimization is reduced by 3–4times. in this study. The workflow shown in Fig. 2 is also feasible Also, due to the nature of gradient-based optimization algo- for gEUD-based optimization algorithm. It is also possible rithm, the optimization solution is repeatable as long as the to use gEUD-based objective to constrain mean dose for initial setting of plan parameters is the same. normal tissue by setting parameter α to small and positive Fig. 8 The transaxial plan dose distribution of lung case in (a) IMRT plan and (b) VMAT plan Yan et al. Radiation Oncology (2018) 13:101 Page 9 of 13 Fig. 9 The DVH of lung case in IMRT plan and VMAT plan value. For more challenging clinical cases with complex better than that of IMRT plans. The optimization time geometry and spatial distribution of anatomical structures, and the number of segments are reduced considerably. the overlapped region between PTV and OARs should be This work demonstrates a way to transform the existing handled respectively and supplementary structures would optimization algorithms originally designed for IMRT to be used. the new algorithm designed for VMAT planning. For the first time the new VMAT plan optimization ap- proach was implemented on an in-house developed plan- Appendix 1 ning platform. The progressive sampling strategy is The objective function of classic FMD is defined as combined with gradient-descent-based FMO algorithm below. for a high-performance VMAT planning system. This X X X work demonstrates a new way to transform the existing min fðÞ x ¼ W D −P ð1Þ ijk ijk ijk optimization algorithms originally designed for IMRT to i j k the new algorithm for VMAT planning. Currently, the dy- namic arcs are planned with a single 360° with 2° angle where D ¼ A x ijk n;ijk n spacing, varying dose rate, and constant gantry speed. Sen- n¼1 sitivity to parameters such collimator angle, couch angle, arc length, and delivery time are not explored methodic- N is the total number of beamlet and x is the fluence ally. However these parameters may provide plan quality of beamlet n. D denotes the dose absorbed by voxel ijk improvements for some cases. In addition, there may be a v . A is the dose deposition coefficient defining the ijk n,ijk need for multiple arcs for cases not included in this study. dose contribution per unit fluence of beamlet n to voxel Determining optimal use of these parameters and the po- (i, j, k). P is the constraint doses for target volume ijk tential automation of the parameter settings are subjects (TV), critical organs (CO), and normal tissue (NT). W ijk to ongoing investigation. It is also necessary to perform is the weighting factors assigned for TV, CO, and NT. phantom verification in order to implement this algorithm P and W are predefined according to clinical require- ijk ijk for clinical application. ment. The objective is to find minimum of f(x) subject to x > 0 for all beamlet n∈N. A non-synchronous updat- Conclusion ing scheme as shown below is used which allows only A new approach was developed which is based on fast one beamlet adjusted at a time. gradient descent algorithm and mixed integer program- ming technique to provide a high-performance VMAT x ðÞ l þ 1 ¼ x ðÞ l þ Δx ; m∈N ð2Þ m m m planning. Results from clinical case studies demon- strated that plan quality of VMAT plans is similar to or Here Yan et al. Radiation Oncology (2018) 13:101 Page 10 of 13 X X X Δx ¼ B x ðlÞþ C , m mn n m n¼1 fxðÞ ¼ w ijk i j k X X X 20 1 3 N N new old − W A A X X ijk m;ijk n;ijk ; new old 4@ A 5 P − A x − A x i j k ijk r;ijk s;ijk r s X X X B ¼ mn r¼1 s¼1 W A ijk m;ijk i j k ð5Þ X X X W A P ijk m;ijk ijk Then the expression (5) can be extended as. i j k X X X C ¼ W A X X X ijk m;ijk fxðÞ ¼ i j i j l is the iteration number and m is the index of beamlet 0 1 0 1 T T adjusted in l–th iteration. Only one beamlet adjusted at one N N new new X X new new @ A @ A iteration. If all beamlets are adjusted one time, there are N ½w P − A x −2w P − A x ijk ijk r;ijk ijk ijk r;ijk r r iterations and it is called one cycle. Usually it takes few cy- r¼1 r¼1 0 1 cles (~ 10) to reach a satisfactory solution for one beam set. T T N N old old X X To accommodate the progressive sampling strategy of old old @ A A x þ w A x s;ijk ijk s;ijk s s VMAT optimization, FMD algorithm was modified to s¼1 s¼1 account for the increasing number of beams. The ð6Þ optimization process starts from the initial beam set consisting of few beams (usually 5), and then the fluence 0 1 new X X X X maps are achieved by FMO and processed by leaf se- new @ A fxðÞ ¼ w P − A x ijk ijk r;ijk quencing. The beam set keeps expanding and new i j k r¼1 beams are added. The beams from previous beam set are 20 1 3 N N called old beams and the beams newly added are called new old X X X X X new old 4@ A 5 − 2w P − A x A x new beams. In a given beam set, the fluence maps of old ijk ijk r;ijk s;ijk r s i j k r¼1 s¼1 beams are fixed and the fluence maps of new beams are 0 1 optimized by FMO for several cycles. The beam set will old X X X X old @ A continuously grow until the maximal number of beams þ w A x ijk s;ijk reaches. For each beam set, there are many beamlets i j s¼1 corresponding to the given number of beams. Assuming ð7Þ T T T there are N beamlets in T-th beam set and N = N new T T T + N where N and N are the numbers of old new old To minimize the objective function in (2) with respect new beamlets for new and old beams. x denotes the flu- to the fluence of m-th beamlet, its gradient at l+ 1 iter- old ence of r-th beamlet of new beam while x denotes the ation should be zero. fluence of s-th beamlet of old beam. For a given voxel v , its dose is calculated based on the fluences of all ∂fx½ ðÞ l þ 1 ijk new old ¼ 0 ð8Þ x and x . r s ∂x new old The gradient of objective function at l+ 1 iteration is X X new old D ¼ A x þ A x ð3Þ ijk r;ijk s;ijk expressed as. r s r¼1 s¼1 0 1 new X X X X Here r and s are the index of beamlets of new and old ∂fx½ ðÞ l þ 1 new @ A ¼ 2w P − A x ðÞ l þ 1 ijk ijk r;ijk beams. Accordingly, the objective function (1) can be re- ∂x i j k r¼1 2 3 written as below. old X X X X X X X old 4 5 −A − 2w A x −A m;ijk ijk s;ijk m;ijk fxðÞ ¼ w ijk i j s¼1 i j 2 0 13 ð9Þ new old X X new old 4 @ A5 P − A x þ A x ijk r;ijk s;ijk r s For m∈, the update formula (2) is. r¼1 s¼1 new new new new x ðÞ l þ 1 ¼ x ðÞ l þ Δx ; m∈N ð10Þ ð4Þ m m m T Alternatively, it can be expressed as. Put (8) and (9) together we have Yan et al. Radiation Oncology (2018) 13:101 Page 11 of 13 0 1 N new X X X X new X X X X w A x ðÞ l A ijk r;ijk r m;ijk new @ A 2w P − A x ðÞ l þ 1 ijk ijk r;ijk r i j k n ¼1 new X X X i j r¼1 B ¼ − 2 3 w A ijk T m;ijk old X X X X i j old X X X 4 5 −A − 2w A x −A ¼ 0: m;ijk ijk s;ijk m;ijk w A A T ijk m;ijk r;ijk s¼1 i j k new i j k new X X X ¼ − x ðÞ l : ð11Þ r w A ijk r¼1 m;ijk i j k new new Replace x ðl þ 1Þ in (11) with [x ðlÞ+Δx] in (10), r r ð15aÞ we have new new new (15a) can also be expressed as B ¼ x ðlÞ 8 2 39 m r¼1 r X X X < new = X X X X w A A new ijk m;ijk r;ijk 4 5 w A P − A x ðÞ l þ A Δx ijk m;ijk ijk r;ijk m;ijk m : ; i j k new new T i j k r¼1 X X X B ,where B ¼ − ; r∈N . mr mr new w A ijk m;ijk old X X X X i j old ¼ w A A x ijk m;ijk s;ijk i j s¼1 old X X X X old ð12Þ − w A x A ijk s;ijk m;ijk i j k s¼1 old X X X B ¼ Equation (12) can be expressed as w A ijk m;ijk i j k X X X X X X X X X − w A A T ijk m;ijk s;ijk w A P − w A ijk m;ijk ijk ijk m;ijk N old i j i j k i j k old X X X ¼ x ð15bÞ w A N ijk s¼1 new m;ijk X X X X new i j A x ðÞ l − w A A Δx r;ijk ijk m;ijk m;ijk m r¼1 i j k old old old old N (15b) can also be expressed as B ¼ x B , old X X X X m s¼1 s ms X X X old ¼ w A A x : ijk m;ijk s;ijk − w A A ijk m;ijk s;ijk i j k s¼1 i j k old T X X X where B ¼ ; s∈N . ð13Þ ms 2 old w A ijk m;ijk i j k With proper arrangement of both sides of Eq. (13), X X X w A P it is ijk m;ijk ijk i j k X X X C ¼ ð15cÞ X X X X X X m w A ijk w P A − w m;ijk ijk ijk m;ijk ijk i j k i j k i j k new X X X X Finally, the updating formula is below. new A A x ðÞ l − w m;ijk r;ijk ijk r¼1 i j k new new new x ðÞ l þ 1 ¼ x ðÞ l þ B x ðÞ l old X X X X m m mr r old 2 A A x ¼ w A Δx : m;ijk s;ijk ijk m r¼1 s m;ijk s¼1 i j k old old old T ð14Þ þ B x þ C ; m∈N ð16Þ ms s new s¼1 Then Δx is calculated as T T new old X X X X X X X X X X X new old w P A − w A A x ðÞ l − w A A x ijk ijk m;ijk ijk m;ijk r;ijk ijk m;ijk s;ijk r s i j k i j k r¼1 i j k s¼1 X X X Δx ¼ : w A ijk m;ijk i j k ð15Þ Yan et al. Radiation Oncology (2018) 13:101 Page 12 of 13 old Δϕ: Angular spacing between two beams. Note that x is solved in 1-th to T-1-th beam set and new I : Contains the original fluence of bixel ij resulted by constant in T-th beam set, while x is variable to be ij solved in T-th beam set. C is calculated based on w FMO. m ij- ,A and P which can be pre-computed prior to the A : Maximal fluence of current beam and is less than k m, ijk ijk old max d . optimization of 1-th beam set. B is calculated based ϕ max old old The variables in the above constraints are defined as: on constantsw ,A ,P and x . B can be ijk m, ijk ijk s m l : Binary variable indicates whether bixel ij of current pre-computed prior to the optimization of T-th beam ij new set. B is calculated based on constantsw ,A ,P beam is blocked by left leaf or not (1: close; 0: open). ijk m, ijk ijk new old r : Binary variable indicates whether bixel ij of current and variable x . B can’t be pre-computed and should ij r m be calculated online. beam is blocked by right leaf or not (1: close; 0: open). b : Binary variable indicates whether bixel ij of current ij Appendix 2 beam is open or not (1: open; 0: close). The objective function of the MILP model is defined as A : Float variable contains the fluence of bixel ij. ij below. l : Binary variable indicates whether bixel ij of pervious ij N M beam is blocked by left leaf or not (1: close; 0: open). XX c c min I −A ð1Þ ij ij l : Binary variable indicates whether bixel ij of follow- ij i¼1 j¼1 ing beam is blocked by left leaf or not (1: close; 0: open). r : Binary variable indicates whether bixel ij of pervious Here I is the original fluence of bixel (beam pixel) ij ij beam is blocked by right leaf or not (1: close; 0: open). after FMO, and A is the unknown fluence of bixel for ij r : Binary variable indicates whether bixel ij of following the single aperture. N and M are the dimensions of flu- ij ence map indexed with ij. For achieving the best uni- beam is blocked by right leaf or not (1: close; 0: open). form fluence mapA in approximating non-uniform The meanings of constraints are explained as below. fluence mapI , the following constraints are required. Constraint (2): requires bixel ij either blocked by leaf or open. c c c l þ r ¼ 1−b ð2Þ ij ij ij Constraint (3): requires mechanical continuity of left leaf. c c Constraint (4): requires mechanical continuity of right leaf. l ≤l ð3Þ ijþ1 ij Constraint (5): requires fluence of beam either zero or c c fluence A . r ≥r ð4Þ ijþ1 ij Constraint (6): requires fluence of beam is less than c c c A ¼ A b ð5Þ fluence I . ij ij ij Constraint (7): requires the distance difference between c c A ≤I ð6Þ ij ij left leaf positions in current and previous beams less than max Δϕ , which is the maximal distance that a leaf can M M X X max c p travel between two neighboring control points or angles. l − l ≤ Δϕ ð7Þ ij ij Constraint (8): requires the distance difference between j¼1 j¼1 left leaf positions in current and following beams less than M M max X X Δϕ which is the maximal distance that a leaf can travel f max c s l − l ≤ Δϕ ð8Þ ij ij between two neighboring control points or angles. j¼1 j¼1 Constraint (9): requires the distance difference between M M right leaf positions in current and previous beams less X X max max r − r ≤ Δϕ ð9Þ than Δϕ, which is the maximal distance that a leaf can ij ij ϕ j¼1 j¼1 travel between two neighboring control points or angles. Constraint (10): requires the distance difference between M M X X max right leaf positions in current and following beams less r − r ≤ Δϕ ð10Þ ij ij s max than Δϕ, which is the maximal distance that a leaf can j¼1 j¼1 travel between two neighboring control points or angles. The superscripts c, p, and f represent current beam, previous beam and following beam. Abbreviations The constants in the above constraints are defined as: CI: Conformity index; CO: Critical organ; DAO: Direct aperture optimization; DVH: Dose volume histogram; DVH: Dose-volume histogram; FMD: Fast d : Maximal dose rate. max monotonic descent; FMO: Fluence map optimization; gEUD: general v : Maximal leaf speed. max equivalent uniform dose; HI: Homogeneity index; IMAT: Intensity modulated s : Gantry angular speed. arc therapy; IMRT: Intensity modulated radiation therapy; MILP: Mixed integer Yan et al. Radiation Oncology (2018) 13:101 Page 13 of 13 linear programming; MLC: Multiple leaf collimator; MU: Monitor unit; 17. Bzdusek K, Friberger H, Eriksson K, Hardemark B, Robinson D, Kaus M. NT: Normal tissue; OAR: Organ at risk; PTV: Planning target volume; Development and evaluation of an efficient approach to volumetric arc QS: Quality score; ROI: Region of interesting; TV: Target volume; therapy planning. Med Phys. 2009;36(6):2328–39. VMAT: Volumetric modulated arc therapy; WE: Weighted error 18. Men C, Romeijn HE, Jia X, Jiang SB. Ultrafast treatment plan optimization for volumetric modulated arc therapy (VMAT). Med Phys. 2010;37:5787–91. 19. Papp D, Unkelbach J. Direct leaf trajectory optimization for volumetric Funding modulated arc therapy planning with sliding window delivery. Med Phys. This work was supported by National Key Projects of Research and Development 2014;41(1):011701. of China (2016YFC0904600) and National Science Foundation of China 20. Craft D, McQuaid D, Wala J, Chen W, Salari E, Bortfeld T. Multicriteria VMAT (NSFC11275270). optimization. Med Phys. 2012;39(2):686–96. 21. Matuszak MM, Steers JM, Long T, McShan DL, Fraass BA, Romeijn HE, Ten Authors’ contributions Haken RK. FusionArc optimization: a hybrid volumetric modulated arc HY, JD, YL designed the study. All authors performed critical review of the therapy (VMAT) and intensity modulated radiation therapy (IMRT) planning manuscript and finally approved the manuscript. strategy. Med Phys. 2013;40(7):071713. 22. Nguyen D, Lyu Q, Ruan D, O'Connor D, Low DA, Sheng KA. Comprehensive Ethics approval and consent to participate formulation for volumetric modulated arc therapy planning. Med Phys. Not applicable. 2016;43(7):4263. 23. Yan H, Yin FF, Guan HQ, Kim JH. AI-guided parameter optimization in inverse treatment planning. Phys Med Biol. 2003;48:3565–80. Competing interests 24. Li RP, Yin FF. Optimization of inverse treatment planning using a fuzzy The authors declare that they have no competing interests. weight function. Med Phys. 2000;27:691–700. 25. Rowbottom CG1, Webb S, Oldham M. Improvements in prostate radiotherapy from the customization of beam directions. Med Phys 1998; Publisher’sNote 25(7 Pt 1):1171–1179. Springer Nature remains neutral with regard to jurisdictional claims in 26. Xia P, Hwang AB, Verhey LJ. A leaf sequencing algorithm to enlarge published maps and institutional affiliations. treatment field length in IMRT. Med Phys. 2002 Jun;29(6):991–8. 27. Langer M, Thai V, Papiez L. Improved leaf sequencing reduces segments or Received: 18 October 2017 Accepted: 17 May 2018 monitor units needed to deliver IMRT using multileaf collimators. Med Phys. 2001;28(12):2450–8. 28. Chen Y, Hou Q, Galvin JM. A graph-searching method for MLC leaf References sequencing under constraints. Med Phys. 2004;31(6):1504–11. 1. Teoh M, Clark CH, Wood K, Whitaker S, Nisbet A. Volumetric modulated arc 29. Dai J, Zhu Y. Minimizing the number of segments in a delivery sequence therapy: a review of current literature and clinical use in practice. Br J Radiol. for intensity-modulated radiation therapy with a multileaf collimator. Med 2011;84(1007):967–96. Phys. 2001;28(10):2113–20. 2. Bedford JL. Treatment planning for volumetric modulated arc therapy. Med 30. Yang R, Dai J, Yang Y, Hu Y. Beam orientation optimization for intensity- Phys. 2009;36(11):5128–38. modulated radiation therapy using mixed integer programming. Phys Med 3. Rao M, Yang W, Chen F, Sheng K, Ye J, Mehta V, Shepard D, Cao D. Biol. 2006;51(15):3653–66. Comparison of Elekta VMAT with helical tomotherapy and fixed field IMRT: 31. Zhu X, Thongphiew D, McMahon R, Li T, Chankong V, Yin FF, Arc- plan quality, delivery efficiency and accuracy. Med Phys. 2010;37(3):1350–9. modulated WQJ. Radiation therapy based on linear models. Phys Med Biol. 4. Yu CX, Tang G. Intensity-modulated arc therapy: principles, technologies 2010;55(13):3873–83. and clinical implementation. Phys Med Biol. 2011;56(5):R31–54. 32. Bohsung J, Gillis S, Arrans R, Bakai A, De Wagter C, Knöös T, Mijnheer BJ, 5. Yu CX. Intensity-modulated arc therapy with dynamic multileaf collimation: Paiusco M, Perrin BA, Welleweerd H, Williams P. IMRT treatment planning:- a an alternative to tomotherapy. Phys Med Biol. 1995;40(9):1435–49. comparative inter-system and inter-Centre planning exercise of the ESTRO 6. Shepard DM, Earl MA, Li XA, Naqvi S, Yu C. Direct aperture optimization: a QUASIMODO group. Radiother Oncol. 2005;76(3):354–61. turnkey solution for step-and-shoot IMRT. Med Phys. 2002;29(6):1007–18. 33. Stein J, Mohan R, Wang XH, Bortfeld T, Wu Q, Preiser K, Ling CC, Schlegel W. 7. Crooks SM, Wu X, Takita C, Watzich M, Xing L. Aperture modulated arc Number and orientations of beams in intensity-modulated radiation therapy. Phys Med Biol. 2003;48(10):1333–44. treatments. Med Phys. 1997;24(2):149–60. 8. Bortfeld T, Webb S. Single-arc IMRT? Phys Med Biol. 2009;54(1):N9–20. 34. Wu Q, Mohan R, Niemierko A, Schmidt-Ullrich R. Optimization of intensity- 9. Otto K. Volumetric modulated arc therapy: IMRT in a single gantry arc. Med modulated radiotherapy plans based on the equivalent uniform dose. Int J Phys. 2008;35(1):310–7. Radiat Oncol Biol Phys. 2002;52(1):224–35. 10. Vanetti E, Nicolini G, Nord J, Peltola J, Clivio A, Fogliata A, Cozzi L. On the role of the optimization algorithm of RapidArc(®) volumetric modulated arc therapy on plan quality and efficiency. Med Phys. 2011;38(11):5844–56. 11. Verbakel WF, Cuijpers JP, Hoffmans D, Bieker M, Slotman BJ, Senan S. Volumetric intensity-modulated arc therapy vs. conventional IMRT in head- and-neck cancer: a comparative planning and dosimetric study. Int J Radiat Oncol Biol Phys. 2009;74(1):252–9. 12. Unkelbach J, Bortfeld T, Craft D, Alber M, Bangert M, Bokrantz R, Chen D, Li R, Xing L, Men C, Nill S, Papp D, Romeijn E, Salari E. Optimization approaches to volumetric modulated arc therapy planning. Med Phys. 2015; 42(3):1367–77. 13. Rao M, Cao D, Chen F, Ye J, Mehta V, Wong T, Shepard D. Comparison of anatomy-based, fluence-based and aperture-based treatment planning approaches for VMAT. Phys Med Biol. 2010;55(21):6475–90. 14. Wang C, Luan S, Tang G, Chen DZ, Earl MA, Yu CX. Arc-modulated radiation therapy (AMRT): a single-arc form of intensity-modulated arc therapy. Phys Med Biol. 2008;53(22):6291–303. 15. Ulrich S, Nill S, Oelfke U. Development of an optimization concept for arc- modulated cone beam therapy. Phys Med Biol. 2007;52(14):4099–119. 16. Cameron C. Sweeping-window arc therapy: an implementation of rotational IMRT with automatic beam-weight calculation. Phys Med Biol. 2005;50(18): 4317–36.
– Springer Journals
Published: May 30, 2018