 Research
 Open Access
 Published:
A fast optimization approach for treatment planning of volumetric modulated arc therapy
Radiation Oncology volume 13, Article number: 101 (2018)
Abstract
Background
Volumetric modulated arc therapy (VMAT) is widely used in clinical practice. It not only significantly reduces treatment time, but also produces highquality treatment plans. Current optimization approaches heavily rely on stochastic algorithms which are timeconsuming and less repeatable. In this study, a novel approach is proposed to provide a highefficient optimization algorithm for VMAT treatment planning.
Methods
A progressive sampling strategy is employed for beam arrangement of VMAT planning. The initial beams with equalspace are added to the plan in a coarse sampling resolution. Fluencemap optimization and leafsequencing are performed for these beams. Then, the coefficients of fluencemaps 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 highquality VMAT plans efficiently and consistently. It presents a new way to accelerate current optimization process of VMAT planning.
Background
VMAT is widely used in cancer treatment in radiation oncology departments due to its high efficiency in treatment delivery [1]. Compared to conventional IMRT, VMAT delivers radiation dose while MLC leaf, dose rate and gantry move simultaneously [2]. It uses less treatment time and total monitor units (MU) compared to conventional IMRT technique [3]. The predecessor of modern VMAT is intensity modulated arc therapy (IMAT) which was first developed by Cedric Yu in 1995 [4, 5]. It was motivated from the idea of delivering plans with a large number of gantry positions. The fluence map of the beam is precalculated and decomposed to several apertures. These apertures is then delivered at a given gantry position by multiple arcs. IMAT requires more arcs which causes extended treatment time. Direct aperture optimization (DAO) was thus proposed to handle the complexity of VMAT optimization using stochastic approach. The stochastic approach is computationally intensive and timeconsuming [6,7,8]. Later, Otto presented an iterative algorithm for VMAT optimization [9, 10]. This algorithm employs progressive sampling strategy and aperturebased algorithm for VMAT optimization where highquality dose plan can be achieved by a single arc. This technique is successfully adopted in commercial treatment planning system and has been applied to a wide range of treatment sites [11].
Compared to IMRT, VMAT presents a complex optimization problem because of the significantly increased number of plan parameters such as gantry angle, MLC leaf position, dose rate, etc. [12, 13]. Various approaches were proposed attempting to solve this problem. Most of them are evolved from approaches originally developed for IMRT [14,15,16,17,18,19,20,21,22]. Currently, most of commercial treatment planning systems provide VMAT optimization function and heavily rely on DAO algorithm. RapidArc (Varian Medical System, Palo Alto, CA, USA) employs a progressive sampling strategy and simulated annealingbased DAO algorithm for VMAT planning [9]. Due to the nature of stochastic approach, the optimization process is timeconsuming and the optimization result is seldom repeatable. As reported the maximal DVH variation could be up to 2% for OARs using RapidArc [10]. SmartArc (Philips Healthcare, Inc., Thornton, CO, USA) utilizes an IMRT plan (consisting of equalspace fields) to initialize the arc plan. The fluence maps of IMRT plan are then segmented and redistributed to their neighboring angles. A local gradientbased optimization approach is used to finetune these apertures to meet the mechanical constraints for actual delivery [16, 17]. In this approach, The fluence maps of all beams are actually derived from few static beams and not fully determined by fluence map optimization (FMO) algorithm.
The purpose of this study is to develop a highefficient optimization approach for VMAT planning. The remainder of this paper is organized as follow. The progressive sampling strategy is first introduced and the optimization process is described. Then, two important components, fluencemap optimization algorithm and leafsequencing algorithm, are explained in brief. These algorithms were developed on an inhouse developed treatment planning system. Three typical clinical cases (head and neck, lung, prostate) were evaluated on both inhouse developed and commercial treatment planning systems. The dosimetric quality, delivery efficiency, and running time of VMAT plans were then assessed. Finally, the advantage and disadvantage of this approach are discussed.
Methods
Progressive sampling strategy
Most of VMAT optimization algorithms model Linac source as a series of static gantry positions. The MLC positions and MU setting is then determined at each gantry position. Provided the complexity of VMAT optimization, to limit the scale of problem a progressive sampling strategy is employed in this study. It starts with a coarse sampling of gantry positions, then move to fine sampling of gantry positions. An example is demonstrated as shown in Fig. 1. At the beginning of VMAT optimization a coarse sampling of the gantry positions is used to model the gantry rotation range. The initial 5 beams with equalspace are chosen for the first optimization stage. The MLC positions and MUs for the initial 5 beams are achieved by the fluencemap optimization and leafsequencing algorithms. Next, the sampling resolution is doubled. And another 5 beams are added to the plan and optimized. The sampling resolution is continuously increased and more beams are added to the plan until the total number of beams is arrived. As shown in Fig. 1, the number of new beams in 5 optimization stages is 5, 10, 20, 40, and 60, while the corresponding beam spacing is 72°, 36°, 18°, 9°, and 6°.
A VMAT optimization approach was developed based on the existing fluencemap optimization and leafsequencing algorithms provided by an inhouse developed treatment planning system [23]. The flowchart of VMAT planning proposed in this study is presented in Fig. 2. First, the basic plan parameters such as couch angle and arc range are set by planner. An initial plan consisting of few beams with equalspace is created and their fluence maps are generated by FMO algorithm. These resulting fluence maps are then processed by the leaf sequencing algorithm, and the optimal MLC positions and MUs are determined. Next, the coefficients of FMO algorithm are adjusted based on plan objective and known fluence maps of these beams. In the next round, the sampling resolution is increased and more beams are added to the plan. The optimization process continues until the maximal number of beams is arrived.
Fluencemap optimization
A gradientbased optimization approach, fast monotonic descent (FMD) algorithm is employed for VMAT optimization. Due to the nature of gradientbased algorithm, the global minimum of objective function can be found quickly. The goal is to find minimum of objective function subject to nonnegative fluences for all pencil beams. A nonsynchronous updating scheme is used which allows only one pencil beam adjusted at a time. The detail of classic FMD algorithm is described in references [23,24,25]. To accommodate the progressive sampling strategy, FMD algorithm was modified to account for the varied number of beams in multiple optimization stages. The formulation of the modified FMD algorithm is described in Appendix 1. In brief, new beams are added to the plan in each optimization stage and their fluence maps are optimized with the fluence maps of old beams known. Here, x^{new} denotes the fluence map of new beam which is unsolved in current stage and x^{old} denotes the fluence map of old beam which is solved in previous stage. For a given voxel v_{ ijk }, its dose is calculated by the sum of dose contributed from x^{new} and x^{old}. For solving x^{new} the coefficients of FMD algorithm is adjusted based the planning dose of target and the dose contribution from x^{old}. As fluence map of x^{new} known, it is processed by a leaf sequencing algorithm and a uniform fluence map for a single aperture is obtained. The uniform fluence map replaces the original fluence map of x^{new} and is used as the final value of x^{new}.
Leaf sequencing
Several leaf sequencing algorithms were developed in the past years such as powers of 2, linear programming, and graphics searching [26,27,28,29]. These algorithms decompose fluence map into several apertures with uniform fluences. Given the dynamic feature of beams in VMAT, the leaf sequencing algorithm has to deal with more mechanical constraints related to gantry speed, leaf speed, and dose rate than those of static beams. Mixed integer linear programming (MILP) was successfully applied in many industrial applications due to its capability in handling largescale optimization problems [30, 31]. In this study, a MILP model was developed using a commercial optimization software package (IBM ILOG CPLEX). It assumes that the fluence map of a beam at a given control point is a uniform fluence map for a single aperture. The aperture of the new beam is determined subjected to the apertures of the two old beams closet to it (one is at front side and one is at end side). As the original flence map of x^{new} is nonuniform, it is necessary to convert it to a uniform fluence by a leafsequencing algorithm with certain constraints. These constraints are mostly related to mechanical limitation of Linac and MLC and described in Appendix 2. Due to the nature of progressive sampling strategy, the beams solved in the earlier stages have larger apertures and weights than the beams solved in the latter stages. This is because the number of beams is less in the earlier stages and more in the latter stages. The MILP model is described in Appendix 2. The goal of this model is to search for an optimal solution, uniform fluence map, from a nonuniform fluence map under certain constraints. The MILP problem is solved using branchandbound technique provided by CPLEX software. As an example of fluence map consisting of 30×30 bixels, the total number of constraints is ~ 1800 and the total number of independent variable is ~ 2000. The processing time is approximately 0.5 s per beam.
Plan evaluation
Three clinical cases were selected for evaluation including head and neck, prostate, and lung. The patient CT images and contours were exported from Pinnacle (Philips Healthcare, Inc., Thornton, CO, USA). Then they are segmented into regions of interest (ROIs) for plan optimization using an inhouse developed treatment planning system. In this study, the simulated leaf width is 0.5 cm and the leaf step is 0.5 cm, which resulted in beamlet size of 0.5×0.5 cm. The range of leaf speed is 0.5–1 cm per degree because the leaf speed is 3–6 cm per second and gantry speed is 6 degree per second. The prescription dose for PTV and dose constraints for OARs are the same for both IMRT and VMAT plans for each case. A single dynamic arc plan consisting of 180 beams and 2° spacing was created and fluence maps were optimized by the proposed VMAT optimization algorithm on the inhouse developed system. The numbers of beams in five optimization stages are 15, 30, 60, 120, and 180, while the corresponding beam spacing is 24°, 12°, 6°, 3°, and 2°. A corresponding IMRT plan consisting of 9 equalspace fields was created and fluence maps were optimized by FMD algorithm on the inhouse developed system. The dosimetric quality of plans is quantified by conformity index (CI), homogeneity index (HI), quality score (QS) [32], and weighted root meansquare error [23], which are defined as below.
CI = V_{TV}/V_{95%}, where V_{TV} is the target volume and V_{95%} is the volume corresponding to 95% dose prescription. If V_{95%} = V_{TV}, CI =1 which is perfect.
HI = D_{5%}/D_{95%}, where D_{5%} and D_{95%} are doses received at least 5 and 95% target volume respectively. if D_{5%} = D_{95%}, HI =1 which is perfect.
Quality score (QS) quantitatively measures the difference between the dose plan and the dose constraint of all anatomical structures. \( \mathrm{QS}\kern0.5em =\kern0.5em \sum \limits_j\left\frac{\left({M}_j{C}_j\right)}{C_j}\right \) if objective is violated by plan dose, where C_{ j } is the jth dosevolume constraint, M_{ j } is the corresponding plan dose. If no objective is violated, QS = 0 which is perfect.
Weighted root meansquare error (WE) measures the difference between the plan dose and prescribed dose of all voxels and represents the value of objective function. \( \mathrm{WE}\kern0.5em =\kern0.5em {\left(\frac{\sum \limits_{i,j,k\in {V}_{TV}}{w}_{ijk}{\left({P}_{ijk}{D}_{ijk}\right)}^2\kern0.5em +\kern0.5em \sum \limits_{i,j,k\in {V}_{CO}}{w}_{ijk}{\left({P}_{ijk}{D}_{ijk}\right)}^2+\sum \limits_{i,j,k\in {V}_{NT}}{w}_{ijk}{\left({P}_{ijk}{D}_{ijk}\right)}^2}{N}\right)}^{\frac{1}{2}} \), where P_{ ijk }, D_{ ijk }, and W_{ ijk } are defined in Eq. 1, and N is the total number of voxels belonging to TV, CO and NT. For a plan which P_{ ijk } = D_{ ijk } for all voxels, WE = 0 which is perfect.
The dose distribution of plans is evaluated based on crosssectional dose distribution and dosevolume histogram. The total number of segments and monitor unit for both plans are recorded. Additionally, estimated arc delivery times and the computation time of plan optimization are recorded. For demonstration purpose, few important ROIs are used in plan optimization. For head and neck case, they are PTV, spinal cord, left and right parotids and mandible. For prostate case, they are PTV, left and right femur heads, bladder and rectum. For lung case, they are PTV, left and right lungs, spinal cord and heart. All tests are performed on a DELL Optiplex N9010 computer, equipped with Intel(R) i7–3770 CPU and 72GB RAM.
Corresponding to the VMAT plans made on inhouse developed system, three VMAT plans made on Pinnacle planning system were assessed. They were made by experienced planners and approved for clinical treatment. These plans are VMAT plans consisting of 2 full arcs with 4° spacing and 6 MV beams computed with convolution superposition algorithm. The dosimetric and delivery statistics of these clinically approved plans are compared to those of IMRT and VMAT plans made on the inhouse developed system. Note that the optimization algorithm, the dose calculation engine, the dosevolume constraints and many others of clinically approved plans are different from those of plans made on the inhouse developed system. The metric, WS, is not calculated for the clinically approved VMAT plans because the weight specification in Pinnacle system is different. For distinguishing clinically approved VMAT plans from the plans made on the inhouse developed system, VMAT plans made on pinnacle system are represented by VMAT_P.
Results
The fluence map and plan dose distribution achieved by progressive sampling strategy is demonstrated in Fig. 3. For each plan, the subplots at the top are the segments of all beams while the subplots at the bottom are the dose distribution in axial view. The grey level of segment indicates the magnitude of its weight. White represents the highest value (255) and black represents the lowest value (0). The corresponding beam arrangements are shown in Fig. 1. In the earlier stages, the values of fluence maps of new beams are larger as there are fewer beams. As more beams included, the values of fluence map of new beams are smaller. The plan dose quality improves quickly in the earlier stages (as shown in Fig. 3b and c), and gradually saturates in latter stages (as shown in Fig. 3d and e.
For the head and neck case, the objective is to maintain uniform dose (60 Gy) to the PTV, while minimize the dose to the parotid gland (V30Gy < 50%) and maximize the sparing of the brainstem (Dmax< 54 Gy), the larynx (Dmax<40Gy), and the mandible (Dmax< 60 Gy) as much as possible. Dose distributions for IMRT and VMAT plans are shown in Fig. 4. DVHs for IMRT and VMAT plans are compared as shown in Fig. 5. The VMAT plan is better than the corresponding IMRT plan. The coverage of PTV is similar for both VMAT and IMRT plans. The OAR sparing from VMAT plan is better for brainstem and mandible, but is similar for larynx and parotids (left and right). The max dose to brainstem is reduced by 10 Gy. Dose distribution of VMAT plan is more uniform with fewer hot and cold spots. The QS, WE, CI, and HI for both plans are shown in Table 1. The WE is reduced by 10% for VMAT plan, while QS, CI and HI are similar for both plans. The total number of segments and monitor units for both plans are shown in Table 1. The number of segment is similar for both plans. The total MU for VMAT plan is reduced by 40%. The computation time is 40 s for IMRT plan and 312 s for VMAT plan as shown in Table 1. The estimated delivery time is between 96 and 192 s based on a Varian Trilogy machine with variable dose rate (300–600 MUs per minute). VMAT_P plan show better plan quality and fewer MUs than VMAT plan, but the number of segments and optimization time are higher.
For the prostate case, the objective is to maintain uniform dose (76 Gy) to the PTV while keep the rectal dose at V50 Gy < 50% and maximize sparing of bladder (V50 Gy < 50%) and femoral head (V50 Gy < 5%). Dose distributions for IMRT and VMAT plans are shown in Fig. 6. DVHs for IMRT and VMAT plans are compared as shown in Fig. 7. The VMAT plan is comparable to the corresponding IMRT plan. The coverage of PTV is similar between VMAT and IMRT plans. For VMAT plan, OAR sparing is better for rectum, femoral heads (left and right) while similar for bladder. Dose distribution of VMAT plan is more uniform and has less high dose region in normal tissue. The QS, WE, CI, and HI for both plans are shown in Table 1. The WE is reduced by 12% for VMAT plan, while QS, CI, and HI are similar for both plans. The total number of segments and monitor units for both plans are shown in Table 1. The number of segment is similar for both plans. The total MU for VMAT plan is reduced by 16%. The computation time is 36 s for IMRT plan and 294 s for VMAT plan as shown in Table 1. The estimated delivery time is between 118 s and 236 s based on a Varian Trilogy machine with variable dose rate (300–600 MUs per minute). VMAT_P plan show better plan quality and fewer MUs than VMAT plan, but the number of segments and optimization time are higher.
For the lung case, the objective is to maintain uniform dose (60 Gy) to the PTV while keep the lung dose at V20 Gy < 25% and maximize the sparing of cord (Dmax< 45 Gy), heart (V30Gy < 40%, V40Gy < 30%) and esophagus (V50Gy < 50%). Dose distributions for IMRT and VMAT plans are shown in Fig. 8. DVHs for IMRT and VMAT plans are compared as shown in Fig. 9. The dose coverage of PTV is similar in both plans. For VMAT plan, dose to heart, esophagus, and cord are slightly higher than those of IMRT plan, while dose to lung is similar. The Dose distribution of the VMAT plan is more uniform and focused on PTV, however the dose distribution of IMRT plan is spread along anteriorposterior direction. The QS, WE, CI, and HI for both plans are shown in Table 1. The WE, QS, CI, and HI are similar for both plans. The total number of segments and monitor units for both plans are shown in Table 1. The number of segment for VMAT plan is reduced by 16%. The total MU is similar for both plans. The computation time is 32 s for IMRT plan and 216 s for VMAT plan as shown in Table 1. The estimated VMAT delivery times are between 116 s and 232 s based on a general linear accelerator with variable dose rate (300–600 MUs per minute). VMAT_P plan show better plan quality and few MUs than VMAT plan, but the optimization time is higher.
Discussions
The results of clinical case study show that plan quality of VMAT plans is similar to or better than that of traditional IMRT plans. The VMAT plans, in many cases, are able to achieve better OAR sparing compared to IMRT plans. VMAT plan has the potential to tailor high volume low dose regions because of its flexibility to spread out the dose at many beam angles. This character reduces the need of specific structures defined to remove hot spots, thus reduces overall planning effort. Plans consisting of more than 2 arcs are not reported in this work, but it is a viable option in clinical practice. We also tested optimizing all beams in multiple arcs which has slightly improved dosimetric quality. These data are not included in this paper. The dose quality of clinically approved plans made on pinnacle planning system is better. This may be caused by the differences of optimization algorithm, dose calculation engine, dosevolume constraints between pinnacle system and inhouse developed system. However, the optimization and delivery time of VMAT plans achieved by our approach is less. It is promising to implement this approach in commercial planning system to accelerate the current process of VMAT optimization.
The selection of beam orientations in initial stages of optimization would affect the final dose considerably. However it could be minimized with increasing number of beams. It was reported that optimizing beam orientations is most valuable for a small numbers of beams (≤5) and the gain diminishes rapidly for higher numbers of beams (≥15) [33]. In several pioneering studies of VMAT optimization it showed 10–15 beams with equalspace would be a better choice in initial stage of optimization [9, 10, 17]. We also tested our algorithm with the different initial beam orientations in the initial stage. It was found that for a higher number of beams (≥15) in the initial stage the final value of objective function is less affected. In the proposed algorithm the number of beams in the initial stage is set to 15, while this value is 10 in RapidArc of Eclipse (Varian Medical System, Palo Alto, CA, USA) and 15 in SmartArc of Pinnacle (Philips Healthcare, Inc., Thornton, CO, USA).
Plan quality can be further improved by increasing the number of beams. However, as more beams included, the computation time for plan optimization and dose calculation could be increased substantially. We speculate that only in certain cases, e.g., with complex target shape requiring substantial leaf motion, it would be beneficial to increase gantry position sampling. It is observed that the number of segments for VMAT plan is less than or similar to that of IMRT plan in three cases. The total MU for VMAT plan is also less than or similar to that of IMRT plan in three cases. Overall, the delivery time for VMAT plan is significantly reduced compared to IMRT plan. Plan optimization times vary case by case but all within 5 min. Compared with the selected commercial treatment planning systems, the running time of plan optimization is reduced by 3–4 times. Also, due to the nature of gradientbased optimization algorithm, the optimization solution is repeatable as long as the initial setting of plan parameters is the same.
The current objective function is quadratic function which can be solved efficiently by gradient decent method. However, for real clinical use it is desired to incorporate more complex objectives and constraints such as dosevolume histogram (DVH) and general equivalent uniform dose (gEUD). For DVH constraints, the current algorithm is applicable with minor modification. Those voxels whose doses exceed the given threshold will be identified when comparing actual DVH to expected DVH. The weights for those voxels will be adjusted automatically in order to minimize the DVH difference. For gEUD constraints, the form of objective function will be changed as described by Wu [34] and can be solved by gradient descent method which is similar to FMD algorithm employed in this study. The workflow shown in Fig. 2 is also feasible for gEUDbased optimization algorithm. It is also possible to use gEUDbased objective to constrain mean dose for normal tissue by setting parameter α to small and positive value. For more challenging clinical cases with complex geometry and spatial distribution of anatomical structures, the overlapped region between PTV and OARs should be handled respectively and supplementary structures would be used.
For the first time the new VMAT plan optimization approach was implemented on an inhouse developed planning platform. The progressive sampling strategy is combined with gradientdescentbased FMO algorithm for a highperformance VMAT planning system. This work demonstrates a new way to transform the existing optimization algorithms originally designed for IMRT to the new algorithm for VMAT planning. Currently, the dynamic arcs are planned with a single 360° with 2° angle spacing, varying dose rate, and constant gantry speed. Sensitivity to parameters such collimator angle, couch angle, arc length, and delivery time are not explored methodically. However these parameters may provide plan quality improvements for some cases. In addition, there may be a need for multiple arcs for cases not included in this study. Determining optimal use of these parameters and the potential automation of the parameter settings are subjects to ongoing investigation. It is also necessary to perform phantom verification in order to implement this algorithm for clinical application.
Conclusion
A new approach was developed which is based on fast gradient descent algorithm and mixed integer programming technique to provide a highperformance VMAT planning. Results from clinical case studies demonstrated that plan quality of VMAT plans is similar to or better than that of IMRT plans. The optimization time and the number of segments are reduced considerably. This work demonstrates a way to transform the existing optimization algorithms originally designed for IMRT to the new algorithm designed for VMAT planning.
Abbreviations
 CI:

Conformity index
 CO:

Critical organ
 DAO:

Direct aperture optimization
 DVH:

Dose volume histogram
 DVH:

Dosevolume histogram
 FMD:

Fast monotonic descent
 FMO:

Fluence map optimization
 gEUD:

general equivalent uniform dose
 HI:

Homogeneity index
 IMAT:

Intensity modulated arc therapy
 IMRT:

Intensity modulated radiation therapy
 MILP:

Mixed integer linear programming
 MLC:

Multiple leaf collimator
 MU:

Monitor unit
 NT:

Normal tissue
 OAR:

Organ at risk
 PTV:

Planning target volume
 QS:

Quality score
 ROI:

Region of interesting
 TV:

Target volume
 VMAT:

Volumetric modulated arc therapy
 WE:

Weighted error
References
 1.
Teoh M, Clark CH, Wood K, Whitaker S, Nisbet A. Volumetric modulated arc therapy: a review of current literature and clinical use in practice. Br J Radiol. 2011;84(1007):967–96.
 2.
Bedford JL. Treatment planning for volumetric modulated arc therapy. Med Phys. 2009;36(11):5128–38.
 3.
Rao M, Yang W, Chen F, Sheng K, Ye J, Mehta V, Shepard D, Cao D. Comparison of Elekta VMAT with helical tomotherapy and fixed field IMRT: plan quality, delivery efficiency and accuracy. Med Phys. 2010;37(3):1350–9.
 4.
Yu CX, Tang G. Intensitymodulated arc therapy: principles, technologies and clinical implementation. Phys Med Biol. 2011;56(5):R31–54.
 5.
Yu CX. Intensitymodulated arc therapy with dynamic multileaf collimation: an alternative to tomotherapy. Phys Med Biol. 1995;40(9):1435–49.
 6.
Shepard DM, Earl MA, Li XA, Naqvi S, Yu C. Direct aperture optimization: a turnkey solution for stepandshoot IMRT. Med Phys. 2002;29(6):1007–18.
 7.
Crooks SM, Wu X, Takita C, Watzich M, Xing L. Aperture modulated arc therapy. Phys Med Biol. 2003;48(10):1333–44.
 8.
Bortfeld T, Webb S. Singlearc IMRT? Phys Med Biol. 2009;54(1):N9–20.
 9.
Otto K. Volumetric modulated arc therapy: IMRT in a single gantry arc. Med Phys. 2008;35(1):310–7.
 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 intensitymodulated arc therapy vs. conventional IMRT in headandneck 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 anatomybased, fluencebased and aperturebased 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. Arcmodulated radiation therapy (AMRT): a singlearc form of intensitymodulated arc therapy. Phys Med Biol. 2008;53(22):6291–303.
 15.
Ulrich S, Nill S, Oelfke U. Development of an optimization concept for arcmodulated cone beam therapy. Phys Med Biol. 2007;52(14):4099–119.
 16.
Cameron C. Sweepingwindow arc therapy: an implementation of rotational IMRT with automatic beamweight calculation. Phys Med Biol. 2005;50(18):4317–36.
 17.
Bzdusek K, Friberger H, Eriksson K, Hardemark B, Robinson D, Kaus M. Development and evaluation of an efficient approach to volumetric arc therapy planning. Med Phys. 2009;36(6):2328–39.
 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 modulated arc therapy planning with sliding window delivery. Med Phys. 2014;41(1):011701.
 20.
Craft D, McQuaid D, Wala J, Chen W, Salari E, Bortfeld T. Multicriteria VMAT optimization. Med Phys. 2012;39(2):686–96.
 21.
Matuszak MM, Steers JM, Long T, McShan DL, Fraass BA, Romeijn HE, Ten Haken RK. FusionArc optimization: a hybrid volumetric modulated arc therapy (VMAT) and intensity modulated radiation therapy (IMRT) planning strategy. Med Phys. 2013;40(7):071713.
 22.
Nguyen D, Lyu Q, Ruan D, O'Connor D, Low DA, Sheng KA. Comprehensive formulation for volumetric modulated arc therapy planning. Med Phys. 2016;43(7):4263.
 23.
Yan H, Yin FF, Guan HQ, Kim JH. AIguided parameter optimization in inverse treatment planning. Phys Med Biol. 2003;48:3565–80.
 24.
Li RP, Yin FF. Optimization of inverse treatment planning using a fuzzy 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; 25(7 Pt 1):1171–1179.
 26.
Xia P, Hwang AB, Verhey LJ. A leaf sequencing algorithm to enlarge 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 monitor units needed to deliver IMRT using multileaf collimators. Med Phys. 2001;28(12):2450–8.
 28.
Chen Y, Hou Q, Galvin JM. A graphsearching method for MLC leaf sequencing under constraints. Med Phys. 2004;31(6):1504–11.
 29.
Dai J, Zhu Y. Minimizing the number of segments in a delivery sequence for intensitymodulated radiation therapy with a multileaf collimator. Med Phys. 2001;28(10):2113–20.
 30.
Yang R, Dai J, Yang Y, Hu Y. Beam orientation optimization for intensitymodulated radiation therapy using mixed integer programming. Phys Med Biol. 2006;51(15):3653–66.
 31.
Zhu X, Thongphiew D, McMahon R, Li T, Chankong V, Yin FF, Arcmodulated WQJ. Radiation therapy based on linear models. Phys Med Biol. 2010;55(13):3873–83.
 32.
Bohsung J, Gillis S, Arrans R, Bakai A, De Wagter C, Knöös T, Mijnheer BJ, Paiusco M, Perrin BA, Welleweerd H, Williams P. IMRT treatment planning: a comparative intersystem and interCentre planning exercise of the ESTRO QUASIMODO group. Radiother Oncol. 2005;76(3):354–61.
 33.
Stein J, Mohan R, Wang XH, Bortfeld T, Wu Q, Preiser K, Ling CC, Schlegel W. Number and orientations of beams in intensitymodulated radiation treatments. Med Phys. 1997;24(2):149–60.
 34.
Wu Q, Mohan R, Niemierko A, SchmidtUllrich R. Optimization of intensitymodulated radiotherapy plans based on the equivalent uniform dose. Int J Radiat Oncol Biol Phys. 2002;52(1):224–35.
Funding
This work was supported by National Key Projects of Research and Development of China (2016YFC0904600) and National Science Foundation of China (NSFC11275270).
Author information
Affiliations
Contributions
HY, JD, YL designed the study. All authors performed critical review of the manuscript and finally approved the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix 1
The objective function of classic FMD is defined as below.
N is the total number of beamlet and x_{ n } is the fluence of beamlet n. D_{ ijk } denotes the dose absorbed by voxel v_{ ijk }. A_{ n,ijk } is the dose deposition coefficient defining the dose contribution per unit fluence of beamlet n to voxel (i, j, k). P_{ ijk } is the constraint doses for target volume (TV), critical organs (CO), and normal tissue (NT). W_{ ijk } is the weighting factors assigned for TV, CO, and NT. P_{ ijk } and W_{ ijk } are predefined according to clinical requirement. The objective is to find minimum of f(x) subject to x_{ n } > 0 for all beamlet n∈N. A nonsynchronous updating scheme as shown below is used which allows only one beamlet adjusted at a time.
Here
\( \Delta {x}_m=\sum \limits_{n=1}^N{B}_{mn}{x}_n(l)+{C}_m \), \( {B}_{mn}=\frac{\sum \limits_i\sum \limits_j\sum \limits_k{W}_{ijk}{A}_{m, ijk}{A}_{n, ijk,}}{\sum \limits_i\sum \limits_j\sum \limits_k{W}_{ijk}{A}_{m, ijk}^2} \)\( {C}_m=\frac{\sum \limits_i\sum \limits_j\sum \limits_k{W}_{ijk}{A}_{m, ijk}{P}_{ijk}}{\sum \limits_i\sum \limits_j\sum \limits_k{W}_{ijk}{A}_{m, ijk}^2} \).
l is the iteration number and m is the index of beamlet adjusted in l–th iteration. Only one beamlet adjusted at one iteration. If all beamlets are adjusted one time, there are N iterations and it is called one cycle. Usually it takes few cycles (~ 10) to reach a satisfactory solution for one beam set.
To accommodate the progressive sampling strategy of VMAT optimization, FMD algorithm was modified to account for the increasing number of beams. The optimization process starts from the initial beam set consisting of few beams (usually 5), and then the fluence maps are achieved by FMO and processed by leaf sequencing. The beam set keeps expanding and new beams are added. The beams from previous beam set are called old beams and the beams newly added are called new beams. In a given beam set, the fluence maps of old beams are fixed and the fluence maps of new beams are optimized by FMO for several cycles. The beam set will continuously grow until the maximal number of beams reaches. For each beam set, there are many beamlets corresponding to the given number of beams. Assuming there are N^{T} beamlets in Tth beam set and N^{T} = N^{T}_{ new } + N^{T}_{ old } where N^{T}_{ new } and N^{T}_{ old } are the numbers of beamlets for new and old beams. x_{ r }^{new} denotes the fluence of rth beamlet of new beam while x_{ s }^{old} denotes the fluence of sth beamlet of old beam. For a given voxel v_{ ijk }, its dose is calculated based on the fluences of all x_{ r }^{new} and x_{ s }^{old}.
Here r and s are the index of beamlets of new and old beams. Accordingly, the objective function (1) can be rewritten as below.
Alternatively, it can be expressed as.
Then the expression (5) can be extended as.
To minimize the objective function in (2) with respect to the fluence of mth beamlet, its gradient at l + 1 iteration should be zero.
The gradient of objective function at l + 1 iteration is expressed as.
For m∈, the update formula (2) is.
Put (8) and (9) together we have
Replace \( {x}_r^{new}\left(l+1\right) \) in (11) with [\( {x}_r^{new}(l) \)+Δx] in (10), we have
Equation (12) can be expressed as
With proper arrangement of both sides of Eq. (13), it is
Then Δx_{ m }is calculated as
(15a) can also be expressed as \( {B}_m^{new}=\sum \limits_{r=1}^{N_{new}^T}{x}_r^{new}(l){B}_{mr}^{new} \),where \( {B}_{mr}^{new}=\frac{\sum \limits_i\sum \limits_j\sum \limits_k{w}_{ijk}{A}_{m, ijk}{A}_{r, ijk}}{\sum \limits_i\sum \limits_j\sum \limits_k{w}_{ijk}{A}_{m, ijk}^2},r\in {N}_{new}^T \).
(15b) can also be expressed as \( {B}_m^{old}=\sum \limits_{s=1}^{N_{old}^T}{x}_s^{old}{B}_{ms}^{old} \), where \( {B}_{ms}^{old}=\frac{\sum \limits_i\sum \limits_j\sum \limits_k{w}_{ijk}{A}_{m, ijk}{A}_{s, ijk}}{\sum \limits_i\sum \limits_j\sum \limits_k{w}_{ijk}{A}_{m, ijk}^2},s\in {N}_{old}^T \).
Finally, the updating formula is below.
Note that \( {x}_s^{old} \) is solved in 1th to T1th beam set and constant in Tth beam set, while \( {x}_r^{new} \) is variable to be solved in Tth beam set. C_{ m }is calculated based on w_{ ijk },A_{m, ijk}and P_{ ijk } which can be precomputed prior to the optimization of 1th beam set. \( {B}_m^{old} \) is calculated based on constantsw_{ ijk },A_{m, ijk},P_{ ijk } and \( {x}_s^{old} \). \( {B}_m^{old} \) can be precomputed prior to the optimization of Tth beam set. \( {B}_m^{new} \) is calculated based on constantsw_{ ijk },A_{m, ijk},P_{ ijk } and variable \( {x}_r^{new} \). \( {B}_m^{old} \) can’t be precomputed and should be calculated online.
Appendix 2
The objective function of the MILP model is defined as below.
Here \( {I}_{ij}^c \) is the original fluence of bixel (beam pixel) after FMO, and \( {A}_{ij}^c \)is the unknown fluence of bixel for the single aperture. N and M are the dimensions of fluence map indexed with ij. For achieving the best uniform fluence mapA^{c}in approximating nonuniform fluence mapI^{c}, the following constraints are required.
The superscripts c, p, and f represent current beam, previous beam and following beam.
The constants in the above constraints are defined as:
d_{max}: Maximal dose rate.
v_{max}: Maximal leaf speed.
s^{ϕ}: Gantry angular speed.
Δϕ: Angular spacing between two beams.
\( {I}_{ij}^c \): Contains the original fluence of bixel ij resulted by FMO.
A^{c}: Maximal fluence of current beam and is less than \( \frac{v_{\mathrm{max}}}{s^{\phi }}{d}_{\mathrm{max}} \).
The variables in the above constraints are defined as:
\( {l}_{ij}^c \): Binary variable indicates whether bixel ij of current beam is blocked by left leaf or not (1: close; 0: open).
\( {r}_{ij}^c \): Binary variable indicates whether bixel ij of current beam is blocked by right leaf or not (1: close; 0: open).
\( {b}_{ij}^c \): Binary variable indicates whether bixel ij of current beam is open or not (1: open; 0: close).
\( {A}_{ij}^c \): Float variable contains the fluence of bixel ij.
\( {l}_{ij}^p \): Binary variable indicates whether bixel ij of pervious beam is blocked by left leaf or not (1: close; 0: open).
\( {l}_{ij}^f \): Binary variable indicates whether bixel ij of following beam is blocked by left leaf or not (1: close; 0: open).
\( {r}_{ij}^p \): Binary variable indicates whether bixel ij of pervious beam is blocked by right leaf or not (1: close; 0: open).
\( {r}_{ij}^f \): Binary variable indicates whether bixel ij of following beam is blocked by right leaf or not (1: close; 0: open).
The meanings of constraints are explained as below.
Constraint (2): requires bixel ij either blocked by leaf or open.
Constraint (3): requires mechanical continuity of left leaf.
Constraint (4): requires mechanical continuity of right leaf.
Constraint (5): requires fluence of beam either zero or fluence A^{c}.
Constraint (6): requires fluence of beam is less than fluence I_{ ij }^{c}.
Constraint (7): requires the distance difference between left leaf positions in current and previous beams less than \( \frac{v_{\mathrm{max}}}{s^{\phi }}\Delta \phi \), which is the maximal distance that a leaf can travel between two neighboring control points or angles.
Constraint (8): requires the distance difference between left leaf positions in current and following beams less than \( \frac{v_{\mathrm{max}}}{s^{\phi }}\Delta \phi \) which is the maximal distance that a leaf can travel between two neighboring control points or angles.
Constraint (9): requires the distance difference between right leaf positions in current and previous beams less than \( \frac{v_{\mathrm{max}}}{s^{\phi }}\Delta \phi \), which is the maximal distance that a leaf can travel between two neighboring control points or angles.
Constraint (10): requires the distance difference between right leaf positions in current and following beams less than \( \frac{v_{\mathrm{max}}}{s^{\phi }}\Delta \phi \), which is the maximal distance that a leaf can travel between two neighboring control points or angles.
Rights and permissions
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.
About this article
Cite this article
Yan, H., Dai, J. & Li, Y. A fast optimization approach for treatment planning of volumetric modulated arc therapy. Radiat Oncol 13, 101 (2018). https://doi.org/10.1186/s130140181050x
Received:
Accepted:
Published:
Keywords
 Volumetric modulated arc therapy
 Fluencemap optimization
 Leafsequencing