A new standard quadratic optimization approach to beam angle optimization for fixed-field intensity modulated radiation therapy

,

Department of Engineering and Applied Physics,University of Science and Technology of China,Hefei 230026,China

Abstract:Beam angle optimization (BAO)largely determines the performance of the fixed-field intensity modulated radiation therapy (IMRT),and it is usually considered as non-convex optimization and a non-deterministic polynomial(NP)hard problem.In this work,BAO is reformulated into a highly efficient framework of standard quadratic optimization.The maximum of beamlet intensities for each incident field as the surrogate variable indicates whether one radiation field has been selected.By converting the function of maximum value in the objective into a set of linear constraints,the problem is solved as standard quadratic optimization via reweighting l1-norm scheme.The performance of the proposed BAO has been verified on a digital phantom and two patients.And the conclusion is drawn:the proposed convex optimization framework is able to find an optimal set of beam angles,leading to improved dose sparing on organs-at-risk (OARs)in the fixed-field IMRT.

Keywords:beam angle optimization；intensity modulated radiation therapy；standard quadratic optimization

1 Introduction

Inverse treatment planning for intensity modulated radiation therapy (IMRT)aims to obtain a prescribed dose distribution on planning target volume (PTV)while sparing organs at risk (OARs).A fully optimized IMRT plan should consider all the system parameters of a clinical linear accelerators as control variables in the optimization process,including beam number,beam angle,multi-leaf collimator (MLC)leaf positions,and monitor unit (MU)for each segment.Convex formulation of such an optimization task,however,appears challenging since most control variables have a non-linear relationship with the delivered dose distribution.In this work,we improve fixed-field IMRT by including beam angled optimization into the inverse treatment planning process,via a newl1-norm minimization approach.

A large number of treatment beams prolongs dose delivery time and therefore increases potential dose errors due to patient motion.On the other hand,it is reported that the dose improvement of a treatment plan diminishes as beam number increases and less than 10 beam angles are often sufficient for IMRT[1].As a small beam number is used in current fixed-field IMRT,the selection of beam angles largely determines the treatment plan quality[2,3].Beam angle optimization (BAO)searches for an optimal set of beam orientations to obtain the best plan quality from all possible beam angle combinations,which is inherently an NP-hard combinatorial optimization problem with no efficient solutions yet.As such,BAO is not ubiquitously implemented in current clinical practice.Instead,beam number is first empirically determined,and beam angles are selected in a trial-and-error fashion.Due to the mathematical complexity of inverse planning in IMRT,empirical tuning of beam angle selection does not guarantee the optimality of treatment plan.For instance,the mathematically optimal beam configuration can be counterintuitive since the extra freedom of intensity modulation compensates for the visually sub-optimal beams[7].To shorten the treatment planning time of IMRT,equiangular beams are used in many radiation therapy scenarios,and the same beam angle setting is typically used for the same disease site on different patients,at the cost of reduced plan optimality[8,9].

BAO for IMRT has been an active research area for decades[6，7，10-14].Many existing BAO methods improve the empirical selection of beam angles by including dosimetric or geometric considerations[12，15-17],and they are not exactly optimization algorithms from a mathematical perspective.For example,some researchers use beam’s-eye-view dosimetry[18]to rank the possible beam orientations according to the quality of dose distribution inside the PTV when the tolerance for OARs is not exceeded.Another similar work[16]uses the ratio of OAR total dose to mean PTV dose as the quality metric for each incident field.The above strategies reduce the computation of BAO by analysing the contribution of individual beam to the overall quality of a treatment plan,which inevitably compromises the optimality of delivered dose distribution due to negligence of multiple-beam interplay[19].

Another category of BAO methods aims to find the optimal beam angles for IMRT using global optimization for a non-convex problem.Existing approaches include simulated annealing algorithms[7,13,19],genetic algorithms[20,21],particle swarm optimization method[22],and multi-objective optimization algorithms[23,24].As a weakness of non-convex optimization with a large solution pool in general,these methods typically require clinically unacceptable long computation and it is theoretically impossible to guarantee the global optimality of the solution due to the existence of multiple local minima[10,25].

Recent developments on optimization methods give rise to non-conventional treatment planning algorithms for IMRT.Sparse optimization was introduced to IMRT treatment planning by Zhu and Xing to obtain a satisfactory dose distribution with a simplified treatment plan[26,27].By minimizing a total-variation objective with quadratic constraints,the algorithm finds piece-wise constant fluence maps with sparse gradients,leading to a highly efficient treatment with a small number of segments.BAO searches for optimal sparse beams in the angular space,which can be formulated as a sparse recovery problem as well.The key challenge of solving BAO via sparse optimization is tond an appropriate control variable for the objective function to indicate the sparsity of beams while still preserving the convexity of the optimization problem.A probably first attempt of sparse optimization based BAO can be found in a recent literature[28].The authors find it difficult to formulate anl1-norm objective and propose a mixedl2,1-norm of beam intensities instead,which is inherently an adaptation of the group-lasso (also called group-sparsity)method[29]widely used in signal processing and statistical learning[30-34]for fearture selection.Theoretically,such a scheme compromises the optimality of BAO for the mixedl2,1-norm objective not only imposes the sparsity on the beam angle level,but also perform an additional smoothing penalty within each field.On the other hand,the authors also admit it that the use of a somewhat complicated group-lasso objective also complicates the computation since the proposedl2,1-norm minimization cannot be solved by either linear or quadratic programming.Some BAO algorithms are recently developed based on the group-lasso method and they both replace thel1-norm with the nonconvex function to better approximate thel0-norm while rataining thel2-normin each field as the control variable to indicate whether this field is selected or not[5,35].The quasil0-norm method (i.e.,l2,0-norm)eliminates thel2-norm penalization within each field during the BAO process,it still employs a complicated iteration scheme due to the nonideal group-sparsity objective.

In this work,by designing a new control variable in the sparse optimization framework,we propose an improved BAO algorithm with anl1-norm regularized quadratic objective.Since the algorithm is in a standard form of quadratic optimization,it accurately finds the global minimum with high computational efficiency.The method performance is demonstrated on one digital phantom,one prostate patient and one head-and-neck (HN)patient.

2 Methods

2.1 Inverse treatment planning of IMRT using l1-norm minimization

(1)

(2)

where the indexidenotes PTV or different OARs,Aiis the beamlet kernel for different structures,λiis the corresponding importance factor[37,38],anddiis the prescribed dose to each structure.The optimized beamlet intensity is finally converted to MLC leaf positions and MUs for different segments,using a leaf sequencing algorithm[39].

In current fixed-field IMRT,a small number of beam angles (typically 5-10)are pre-determined before the optimization of beamlet intensities.In this work,we aim to include a large number of beam angles from a full rotation into the beamlet optimization framework and use sparse optimization to automatically select the optimal beam combination,so the new optimization algorithm takes the following form ofl1-norm minimization:

(3)

(4)

(5)

(6)

(7.1)

(7.2)

(7.3)

whereIis an identity matrix with the size ofNbyN,andB0is expressed as

(7.4)

2.2 The proposed BAO with a reweighting scheme

Derived froml1-norm minimization,the optimization problem (6)sacrifices sparsity of the optimized solution for computational efficiency.At the cost of increased computation,the non-convexl0-norm minimization enhances the solution sparsity and therefore reduces the number of required beams.In this paper,we propose to balance the computational efficiency and the solution sparsity via a series of reweightedl1-norm minimization,a strategy commonly used in differentl1based sparse optimization problems[28,40,41].

The reweighting scheme approximatesl0-norm minimization by adaptively assigning large weights to the optimized vector elements with small values in the previous iteration ofl1-norm minimization[41].In each iteration,the optimization takes the following form:

(8)

(9)

AlgorithmBAO using quadratic optimization with reweighting

repeat

1.Solve the optimization problem (8);

untilNangis no longer larger thanNA.

2.3 Evaluation

We evaluate the proposed BAO method on a digital phantom,a prostate patient and a head-and-neck patient.For all evaluation studies,we consider 40 equiangular beams in a full rotation as the candidates of all available beam orientations.The PTV is centered at the axis of rotation,with a source-to-axis distance (SAD)of 100 cm.In the FMO for fixed-field IMRT,eacheld targets the center of PTV,and contains 20 by 20 beamlets,with a beamlet size of 5 mm by 5 mm at SAD,while during the BAO process,the beamlet size is downsampled to 1cm by 1cm to reduce the memory usage.To save computation in the Monte Carlo simulation of the dose kernel (i.e.,the matrixAi),the CT data are downsampled to a voxel size of 3.92 mm by 3.92 mm by 2 mm.All the algorithms are implemented in Matlab,using CVX,an open-source optimization software[42].On a 2.4 GHz workstation with 28 cores,the proposed BAO takes 3 min on the digital phantom,5 min on the prostate patient,and 6 min on the head and neck patient.

A theoretically optimal set of beam angles is difficult to derive on clinical cases,since it is dependent on the geometries of structures (i.e.,the dose kernelAi)as well as the parameters of treatment planning (i.e.,λi,diandβ).The study of digital phantom with a known optimal set of beam angles is designed to test the proposed BAO algorithm.We implement the conventional IMRT planning (i.e.,the optimization problem (2))with all beam angles included for comparison.In the patient studies,we investigate the dose performance of fixed-field IMRT using the proposed BAO and a set of equiangular beam angles.In addition to the final dose distributions,we compare the dose-volume-histogram (DVH)curves of OAR for different plans with a similar dose coverage on PTV.

A particular difficulty occurring in the design of patient studies is that,on the same patient,the parameters of IMRT planning,especially the importance factors (i.e.,λi),need to be fine-tuned for different sets of beam angles to achieve clinically acceptable treatment plans,leading to unfair comparisons of different algorithms.For a comprehensive evaluation of method performance,we consider fixed-field IMRT planning with the proposed BAO a multi-objective optimization problem,with the following objectives of minimization:

Figure 1.The validatio study o a simulatio phantom. (a)The simulatio phantom with PTV and OAR. The quantity of for each incident field without (b)and with (c)the proposed BAO.

PTV dose objective:

OAR dose objective:

beam number objective:

3 Results

3.1 The digital phantom study

Figure 1 (a)shows the digital phantom used in the simulation.This is a cylinder phantom composed of pure water to simulate the tissue and the diameters of the phantom and PTV are 20 cm and 2.5 cm respectively.The phantom is rotationally symmetric,and therefore the optimal beam angles are only dependent on the relative positions of PTV and OAR.Six passages at randomly selected angles of 0°,54°,81°,153°,216°,315° are designed,on which the radiation beams reach PTV without passing through OAR.As such,these six beam angles are considered the optimal orientations in this study.

Figure 2.Results of the prostate patient study. Pareto frontiers of the PTV dose objective (φ1 )and the OAR dose objective(φ2 )for the same number of beam angles using a equiangular pla and the proposed BAO.

3.2 The prostate cancer patient study

Figure 4.Compariso of DVH curves for the BAO plan (dashed)and the equiangular plan (solid). (a)shows the DVHs i PTV and overall OAR,and (b)shows the DVH i each critical structure.

Figures 2,3,and 4 show the results on the prostate patient.By tuning algorithm parameters,the proposed BAO is able to select different numbers of beam angles.With the same PTV dose coverage,Figure 2 compares the Pareto frontiers of fixed-field IMRT with five beam angles using an equiangular plan and the BAO method.It is seen that the proposed BAO substantially improves the dose performance over an equiangular plan with reduced dose objective values on both PTV and OARs.

The improved dose sparing on OARs achieved by the proposed BAO is better seen in the comparison of dose distributions in Figure 3 and DVHs in Figure 4.In this patient case,we find thatve beam angles (9°,36°,117°,234°,324°)obtained by the proposed BAO successfully achieve a clinically acceptable dose coverage on PTV.We then compare with the conventional IMRT planning usingve equiangular beams starting at 0°.The algorithm parameters are tuned such that both plans obtain the same dose performance on PTV.With the freedom of beam angle selection,the proposed BAO favors beam passages reaching PTV without intersecting OARs (see Figure 3)and therefore significantly improves dose sparing on OARs.The superior performance of the BAO plan over the equiangular plan is further seen in the DVH comparison of Figure 4.The proposed BAO reduces the overall dose exposure on OARs by 30.53%.Table 1 summarizes and compares the results of the BAO and equiangular treatment plans with the clinical acceptance criteria,which also demonstrates the validity of the proposed BAO.

Table 1.Prostate pla objectives and results.

3.3 The head-and-neck cancer patient study

Figure 5.Results of the head-and-neck patient study.Pareto frontiers of the PTV dose objective (φ1 )and the OAR dose objective(φ2 )for the same number of beam angles using a equiangular pla and the proposed BAO.

A similar performance of the proposed BAO is observed on the head-and-neck patient,as shown in Figures 5-8.Figure 5 compares the Pareto frontiers of the PTV and the OAR dose objectives for the IMRT plans using seven beams generated from the proposed BAO with seven and nine equiangular beams starting at 0°.It should be noticed that the seven beam angles selected by the BAO outperforms not only the equiangular seven beams but also the nine beams,which further confirms the necessity and essentiality of BAO for the complex and non-intuitive scenarios like head-and-neck cases.

The superior performance of the BAO plan for OAR avoidance is visually verified in the comparison of dose distributions in Figure 6.The algorithm parameters are tuned to obtain the same dose coverage on PTV in all three plans using seven beams selected by BAO (0°,45°,63°,99°,126°,261°,288°),equiangular seven and nine beams.We find that,compared with the prostate patient case,the results of BAO is more counterintuitive on the head-and-neck patient,mainly due to the geometric complexity of PTV and OARs.In this case,our algorithm selects the optimal beams distant from equiangular directions to better adapt the strip-shape of PTV as well as to avoid the OARs.The improved dose sparing on OARs in the BAO plan is seen in the DVH comparison in Figure 7.The proposed BAO reduces OAR dose by 16.1% and 12.3% from that of the equiangular plans using seven and nine beams,respectively.Table 2 compares the results of the BAO and equiangular plans with the acceptance criteria which indicates that the plans are all clinically acceptable and the superiority of the BAO plan.

Table 2.HN pla objectives and results.

4 Conclusions

Figure 6.Compariso of dose distributions of different slices o the HN patient with plans using equiangular seve beams (a,b),nine beams (c,d)，and seve beams using the proposed BAO(e,f),respectively.

In this paper,we propose a new BAO algorithm to improvexed-field IMRT.The problem of optimal angle selection is first formulated asl1-norm minimization based on the sparse optimization,and then converted into a highly efficient framework of standard quadratic optimization.On a digital phantom,the proposed BAO successfully finds the theoretically optimal set of beam angles from more than 3000000 possible combinations.Our algorithm reduces the delivered dose on OARs by 30.53% and 12.3% to 16.1% on a prostate patient and a head-and-neck patient,respectively,compared with that of equiangular IMRT plans with the same PTV dose coverage.

The optimal set of IMRT beam angles varies on different cancer patients.In the era of patient-specific radiation therapy,beam angle selection remains as one of very few procedures missing in the current clinical practice of fixed-field IMRT,mainly due to its high complexity of implementation.Compared with those of existing researches on non-convex or convex BAO algorithms,the main contribution of our work is to show that BAO can be accurately performed using a simple and efficient framework of standard quadratic optimization.As such,the proposed BAO method is practical for improving IMRT dose performance especially on patients with irregular shapes and/or positions of PTV and/or OARs (i.e.,head-and-neck patients).Larger dose benefits achieved by BAO are expected on non-conventional IMRT scenarios (e.g.,non-coplanar IMRT[17,35]),where beam angles have additional degrees of freedom.Our algorithm is therefore more attractive in these applications for its mathematical simplicity.

There remains a common issue in all existing group sparsity regularization methods.As mentioned in Section 2.3,the important factors for different OARs (λiin Eq.(2))need to be fine-tuned for different beam orientations selected,in other words,pre-fixed important factors can affect the final optimized beam angles.For example,if the weightingλifor OARiis very large,those beam angles directly passing through OARiwill be discarded in the optimized orientations.However,if weighting factors for all OARs are very large,the prescribed dose distribution inside the PTV may be impossible to be reached.So the best solution is to discard the dose objective during the BAO process,and the dose volume prescription should be implemented as constraints of the optimzation problem.In this way,the BAO method directly minimize the number of beam angles when the dose volume constraint is always satisfied.

Acknowledgments

This work is supported by the National Key R&D Program of China (2016YFC0101400),the National Natural Science Foundation of China (81671681),and the Fundamental Research Funds for the Central Universities (WK2030040089).

Conflictofinterest

The authors declare no conflict of interest.

Authorinformation

PENGJunbois currently a PhD student in Medical Physics at the George W.Woodruff School of Mechanical Engineering,Georgia Institute of Technology.His research mainly focuses on the sparsity optimization and its applications in medical imaging and radiation therapy.

ZHULei(corresponding author)received his PhD from the Department of Electrical Engineering at Stanford University.He is currently a professor at University of Science and Technology of China.His research interests include image reconstruction,image-guided therapy,and fluorescence imaging.