- Email: [email protected]
An Approach for Practical Multiobjective IMRT Treatment Planning
Int. J. Radiation Oncology Biol. Phys., Vol. 69, No. 5, pp. 1600–1607, 2007 Copyright Ó 2007 Elsevier Inc. Printed in the USA. All rights reserved 0360-3016/07/$–see front matter
doi:10.1016/j.ijrobp.2007.08.019
PHYSICS CONTRIBUTION
AN APPROACH FOR PRACTICAL MULTIOBJECTIVE IMRT TREATMENT PLANNING DAVID CRAFT, PH.D. TAREK HALABI, PH.D., HELEN A. SHIH, M.D., AND THOMAS BORTFELD, PH.D. Department of Radiation Oncology, Massachusetts General Hospital and Harvard Medical School, Boston, MA Purpose: To introduce and demonstrate a practical multiobjective treatment planning procedure for intensitymodulated radiation therapy (IMRT) planning. Methods and Materials: The creation of a database of Pareto optimal treatment plans proceeds in two steps. The first step solves an optimization problem that finds a single treatment plan which is close to a set of clinical aspirations. This plan provides an example of what is feasible, and is then used to determine mutually satisfiable hard constraints for the subsequent generation of the plan database. All optimizations are done using linear programming. Results: The two-step procedure is applied to a brain, a prostate, and a lung case. The plan databases created allow for the selection of a final treatment plan based on the observed tradeoffs between the various organs involved. Conclusions: The proposed method reduces the human iteration time common in IMRT treatment planning. Additionally, the database of plans, when properly viewed, allows the decision maker to make an informed final plan selection. Ó 2007 Elsevier Inc. Multiobjective, Inverse planning, Optimization, IMRT, Linear programming, EUD, Multicriteria, Pareto.
such plans forms the Pareto surface, and an MCO plan database consists of a sampling of Pareto optimal plans (7). With each structure (targets and organs at risk) in a plan, we associate an indicator function that maps the dose distribution in that structure to a scalar value. From the perspective of a treatment planner, after such functions are defined, the easiest procedure is to form a Pareto surface of plans with each of the indicator functions as objective functions. This is problematic in practice, because it results in a large Pareto surface with many plans that are not of clinical interest. To overcome this, we would like to constrain some or all of the indicator functions. This can lead to a human iteration loop because it is not, in general, possible to determine a set of mutually satisfiable and useful constraint levels for structures. In this article, we address these issues with a two-step process, in which the first step attempts to find such a set of constraints and the second step generates the Pareto surface based on these constraints.
INTRODUCTION Several authors have discussed the use of multiobjective optimization, also called multicriteria optimization (MCO), for intensity-modulated radiation therapy (IMRT) treatment planning (1–6). There are typically two motivations for MCO: to reduce the time-consuming human iteration loop required to come up with a satisfactory plan and to provide tradeoff information (in the form of a database of diverse treatment plans) that allows the decision makers to make an informed final plan choice. To simultaneously satisfy these two goals, one needs to determine a method of creating a database of good plans in a way that does not require timeconsuming human iterations. This is the topic we address in this article. Given that the most suitable plan for a particular patient is impossible to articulate before seeing what is possible regarding treatment plans, we prefer the multiobjective paradigm of a posteriori preference articulation, in which the doctor makes a decision after seeing the possible plans. This requires the computation of a database of treatment plans. We use the terms database of treatment plans and Pareto surface somewhat interchangeably. A treatment plan is Pareto optimal for a given set of objectives and constraints if there is no other feasible treatment plan that is at least as good in all the objectives and strictly better in at least one objective. The set of all
METHODS AND MATERIALS Radiation therapy planning computed tomography scans and physician-defined target and nontarget volumes are read into computational environment for radiotherapy research (CERR) (8) using CERR’s Radiation Therapy Oncology Group import facility. Within CERR, beam angles are selected and the dose influence matrix D This work was supported by NCI Grant 1 R01 CA103904-01A1: Multi-criteria IMRT Optimization. Conflict of interest: none. Received April 3, 2007, and in revised form July 31, 2007. Accepted for publication Aug 1, 2007.
Reprint requests to: David Craft, Ph.D., Department of Radiation Oncology, Harvard Medical School, Massachusetts General Hospital, 0 Grove St., Boston, MA 02114. Tel: (617) 724-3658; Fax: (617) 724-0368; E-mail: [email protected] 1600
An approach for practical multiobjective IMRT treatment planning d D. CRAFT et al.
that maps beamlet intensities x to the dose vector d is calculated. This information is made available to the optimization software, AMPL/CPLEX, by way of text files. For the optimizations, we use linear programming because of its ability to solve problems to provable optimality and because it handles hard constraints exactly. Any optimization involving linear functions and piecewise linear convex objectives and constraints can be written as a linear program, a feature that has been used in a variety of ways in IMRT optimization (3, 9, 10). We use piecewise linear convex functions for both the target and healthy structure indicator functions. For target modeling, we use a piecewise linear convex function that linearly penalizes doses below the prescription dose. For healthy structures, we fit the generalized equivalent uniform dose model (EUD) (11) with piecewise linear convex approximations. Generalized EUD for an organ is given by P EUD ¼ ð1=n i di a Þ1=a , where n is the number of voxels in the structure, di is the dose to voxel i, and a is a parameter. For MCO, we can drop the ath root in that formula, as explained elsewhere (5) (essentially, minimizing x subject to a set of constraints is the same as minimizing xa for any positive a). We fit the resulting dose summation by fitting the individual terms dia over the appropriate dose range, and summing the resulting fits. The functions dia are fit using three-segment piecewise linear functions (Fig. 1). Table 1 contains a values (12, 13) and the best fit (minimal area between curves) locations and relative weights of the linear sections used for each healthy structure modeled. We also employ hard constraints for maximum doses in healthy structures and minimum and maximum doses for targets. Step 1 of the Pareto database creation consists of setting aspiration levels for each indicator function and finding a treatment plan that is near to satisfying all these aspirations. Letting Fk denote the kth indicator function (Fk is the piecewise linear EUD for healthy structures and the linear underdose function for targets) and Ak denote the maximum desired level for this function, step 1 consists of solving the following optimization problem:
1601
Table 1. Values for the optimization function parameters Organ Brain case Brain stem Optic nerve Optic chiasm Eye Prostate case Rectum Bladder Femoral heads Lung case Lung Esophagus Spinal cord
EUD ‘‘a’’ parameter
Pivot 1
Pivot 2
Weight 2/ Weight 1
6.25 4 4 3.3
37.7 25.7 25.7 5.37
53.8 45.1 45.1 10.67
2.85 2.22 2.22 1.92
6.3 3.8 4
52.5 36.9 19.3
74.6 66.7 33.8
2.86 2.14 2.22
1.2 5 10
2.24 30.5 39
29 47.5 48.7
0.25 2.25 3.36
EUD a values are from previous work (12, 13). We fit the function f(d) = da with a three-segment piecewise linear convex function over the appropriate range for that particular structure. An example fit is shown in Fig. 1. Abbreviation: EUD = equivalent uniform dose.
minimize
P ðFk ðdÞ Ak Þþ k
(1)
subject to d ¼ Dx x$0:
The ($)+ operator is shorthand for max(0, $). Because this is a piecewise linear convex function, this can be written as a linear program. In other words, Formulation 1 tries to find a solution x, which gets every function Fk at or below its aspiration level Ak. Notice that this program is always feasible. We occasionally put additional constraints—for example, maximum dose constraints on all voxels and minimum dose constraints on the targets—at the risk of having an infeasible program. The aspiration levels can be chosen based on what experience suggests is a mutally satisfiable, although possibly aggressive set of goals. In this report, to avoid seemingly ad hoc aspiration levels, we use aspiration levels derived from setting the normal tissue complication probability (NTCP) to 0.1%, using data from Burman (12), and use an aspiration of 0 underdosing for the targets. The solution to the Step 1 optimization problem is used to determine a constraint set for Step 2, which is the generation of the Pareto database. For example, if the optimal value of the Step 1 problem is 0, it means all the aspirations were met and can be used as constraints in Step 2. On the other hand, for any aspirations that is not met, we know how far we are from meeting it based on examining the Step 1 solution. This can be used to form a new bound that is Table 2. Optimization information for the three cases
Fig. 1. Fitting the dose power term of the equivalent uniform dose function with a piecewise convex linear function. For comparison, we also plot the Emami (20) data in this graph, which shows good agreement.
Case
No. voxels
No. beamlets
Brain Prostate Lung
50,708 77,287 114,033
680 550 567
Voxel size dx, dy, dz (cm)
Beamlet size (cm)
0.2, 0.2, 0.25 0.3, 0.3, 0.25 0.3, 0.3, 0.25
0.5 0.5 1.0 1.0 1.0 1.0
Critical structures and unclassified tissue are downsampled in each case by a factor of 2 or 4. Targets are not downsampled. Number of voxels reflects the number of voxels used in the optimization (with downsampling), but all results (e.g., DVHs) are reported on the full voxel grid.
I. J. Radiation Oncology d Biology d Physics
1602
Volume 69, Number 5, 2007
Table 3. Brain case Brain optimization (PTV prescription: 59.4 Gy) Step 1
Step 2
Organ
EUD aspiration level (Gy)
Amount by which aspiration was violated in Step 1 (Gy)
EUD constraint for Step 2 (Gy)
Traded off against PTV?
EUD range in database (Gy)
Brain stem Optic chiasm L optic nerve R optic nerve L orbit R orbit Unclassified tissue
36 36 36 36 3 3 15
0 0 0 0 0 0 4
NA 36 NA 36 3 3 19
Yes No Yes No No No No
32–37.5 NA 19.5–27.5 NA NA NA NA
The aspirations for all of the healthy structures except for the unclassified tissue were met, and held tight at those achievable levels in Step 2, the creation of the Pareto database. Abbreviations: PTV = planning target volume; EUD = equivalent uniform dose. attainable. Softening attainable constraints further will yield a more diverse Pareto surface. It is up to the decision makers which functions to trade off in Step 2 and how to employ hard constraints, which can be used to keep any structure within reasonable limits. We use the algorithm PGEN (Pareto Generation) (7) to compute the Pareto database. The plan databases are navigated by means of a Pareto surface navigator, compatible with CERR, which is displayed in Figure 2 for the lung case. This navigator does linear plan interpolation, which allows one to smoothly transition between plans in the database. Exact EUDs are computed and displayed by the navigator. During the navigation process, a planner may constrain any of the functions to control the part of the surface explored. For each slider, the navigator also reports a resitance value, which is measure of the cost (i.e., the detriment to the other objectives) of pulling that slider down from the current position.
RESULTS For each case, we use seven equispaced beams. Beamlet and voxel discretization information is shown in Table 2. Sampling by a factor of 2 or 4 is used for some healthy structures to speed up the optimizations, but all results are reported
on the original grids. Individual plan optimizations take anywhere from 3 min (brain and prostate) to 20 min (lung) and the databases computed contain approximately 30 plans each. To set aspirations for the normal structures in Step 1, we use the EUD value that corresponds to NTCP of 0.1% (12). These levels are shown in Tables 3–5. Tables 3–5 show our choices for objectives and constraints for three clinical cases of a brain tumor, prostate cancer, and lung cancer, respectively. In general, with familiarity for achievable EUD levels on a case-by-case basis, a treatment planner might use a more customized approach, but this recipe approach works well for the cases herein. We use a hard maximum dose constraint of 115% of the prescription dose for every voxel. A hard minimum of approximately 90% prescription is used for the target voxels. Step 1 results provide us with a set of EUD constraints that can be simultaneously satisfied, but for the prostate and lung cases, to create a wide range of plans in the Pareto database, we soften these achievable values by 10%. In the case of the lungs, for even more exploration freedom we let each lung EUD go as high as 30 Gy. We display the
Table 4. Prostate case Prostate optimization (PTV prescription: 79 Gy, seminal vesicles: 60 Gy) Step 1
Step 2
Organ
EUD aspiration level (Gy)
Amount by which aspiration was violated in Step 1 (Gy)
EUD constraint for Step 2 (Gy)
Traded off against PTV?
EUD range in database (Gy)
Bladder Rectum L femoral head R femoral head Unclassified tissue
52 42 40 40 15
0 12.4 0 0 15
57 60 30 30 33
Yes Yes No No Yes
46–57 53–60 NA NA 27–33
The rectum overlaps the PTV and as expected the aggressive aspiration for the rectum was not met in Step 1. Step 2 produces a four dimensional Pareto surface (bladder, rectum, unclassified tissue, and PTV coverage). Abbreviations: PTV = planning target volume; EUD = equivalent uniform dose.
An approach for practical multiobjective IMRT treatment planning d D. CRAFT et al.
1603
Table 5. Lung case Lung optimization (prescription: 66 Gy) Step 1
Step 2
Organ
EUD aspiration level (Gy)
Amount by which aspiration was violated in Step1 (Gy)
EUD constraint for Step 2 (Gy)
Traded off against PTV?
EUD range in database (Gy)
Spinal cord Esophagus R lung L lung Unclassified tissue
30 44 10 10 15
11 0 0 0 10
45 48 48 30 27
Yes Yes Yes Yes Yes
40–45 18–38 18–38 1–10.5 21–27
Abbreviations as in Table 4. To create a database with a large diversity of plans, the lungs were allowed to go up to an EUD of 30 Gy, even though the aspiration of 10 Gy was met in Step 1. Dose escalation to the tumor is of paramount concern for lung patients: relaxing the lung constraint thus allows the physician to see the effects of increasing the target dose.
resulting diversity of the plans by plotting the extreme dose– volume histograms (DVHs) in Figures 3–5 and by reporting the associated range of EUD values in Tables 3–5. Below the DVH range figures, we have selected a single well-balanced plan from each of the databases. In the following section, we discuss the patient details and present the tradeoff results for three clinical cases. Detailed
Fig. 2. Screenshot of the Pareto surface navigator used to interactively slide through the plans in the database. The navigator is Matlab based and operates within CERR, where dose displays are updated automatically after the user slides a slider. The navigator is available from the authors or as a free download at the MathWorks.com website.
tradeoff analyses (which are most useful with state-of-theart NTCP and tumor control probability models) are beyond the scope of this work, but will be reported for specific disease sites in future publications. The first case was a 59-year-old man with a hemangiopericytoma of the brain involving the left middle cranial fossa. After presentation with left facial numbness, imaging revealed a 4-cm enhancing tumor along the skull base. The tumor abutted the left cavernous sinus, invaginated the medial temporal lobe, and eroded the left petrous apex. After subtotal resection, he was recommended adjuvant radiation therapy and was treated with IMRT. Prescribed tumor dose was 59.4 Gy. Figure 3 depicts the tradeoffs between the left optic nerve, the brainstem, and the planning target volume (PTV). The constraints for the other structures were easily met in Step 1, and not in difficult planning areas, so were not included in the tradeoff. The second case was 73-year-old man with an early stage prostate cancer diagnosed after detection of a palpable nodule on routine rectal exam. After biopsy confirmation of prostate cancer, he chose to undergo definitive radiation therapy. He was treated with IMRT to a dose of 79 Gy. Figure 4 demonstrates the anterior and posterior rectum displayed as separate DVHs, but they are optimized together. By modeling the unclassified tissue with an EUD parameter of a = 3.8 (the value for the bladder) and including it in the tradeoff, we are able to control general dose to the unclassified tissue, and in particular ‘‘streaking’’ of high-dose values far into normal tissues, because the relatively high a parameter makes the function sensitive to small volumes of high dose. A hard minimum of 70 Gy is applied to the PTV. The DVHs plotted by CERR show some voxels with slightly less than 70 Gy, but this is an artifact of the CERR DHV computation, and every PTV voxel gets 70 Gy or more in every plan. The third case was a 72-year-old woman with a non–smallcell lung cancer. She presented with a 6.5 cm right upper lobe tumor abutting the posterior chest wall and underwent resection of the tumor and postoperative chemotherapy. Intraoperatively, her tumor was adherent to the chest wall and the surgical margin was focally positive. Given her high risk of
1604
I. J. Radiation Oncology d Biology d Physics
Volume 69, Number 5, 2007
Fig. 3. Along the top row, extreme dose–volume histograms (DVHs) for each structure are displayed to demonstrate the range of plans available in the database. The lower left DVH set is for a selected single plan from the database, with DVHs for all structures shown. A transversal dose wash for this plan is shown in the bottom right.
local recurrence, she was recommended adjuvant radiation therapy and was treated to the tumor bed to 66 Gy using IMRT. The six-dimensional tradeoff for her case is shown in Figure 5. There is considerable tradeoff between dosing the two lungs: dose can be brought to the tumor through either lung almost entirely. Reducing the dose to the right lung, where the tumor occurs, creates a more conformal plan, but the high-dose band that comes through the left lung in this case might be undesireable. For the spinal cord, we plan on a 3-mm expansion of the cord. This cord planning volume overlaps with the PTV, and the minimum dose of 55 Gy to the PTV results in the adjacent section of the spinal cord receiving that dose as well. If this were deemed unacceptable, the planner would reduce the minimum required dose to the PTV. DISCUSSION AND CONCLUSION The goal of MCO for IMRT optimization is to reduce treatment planning time while making treatment planning more intuitive. This is achieved by constructing a database of Pareto optimal plans that the treatment planner can navigate through interactively, thereby acquiring a feel for the tradeoffs involved before selecting a final plan. This navigation
depends on mixing plans in the database in real time, and is more manageable than the Pareto visualization aid presented by Aubry et al. (6), which displays and color codes the plans in the database but does not allow plan combinations. The primary challenge addressed by this work is how to efficiently create a database of Pareto optimal treatment plans that is both diverse and populated with clinically useful plans. These criteria are in conflict with another, because too much diversity means many plans that are far away from a plan a physician would give. This challenge has been addressed by breaking down the task into two steps. The first step establishes a balanced plan that attempts to satisfy a reasonable (meaning possibly achievable) set of constraints. The solution to this problem yields a set of simultaneously achievable constraints. The user can then set true dose constraints rather than provisional constraints often used in IMRT planning (14). The achievable constraints from Step 1 are then used in the Pareto database generation, Step 2. Without Step 1, one might either specify a set of constraints that leads to an infeasible problem, and have to iterate to remedy this, or specify a set of constraints that is too loose, resulting in a database with large percentage of plans outside the domain of what might be delivered. The technical details of Step 2 can be found elsewhere (7).
An approach for practical multiobjective IMRT treatment planning d D. CRAFT et al.
1605
Fig. 4. The top row shows the extreme dose–volume histograms (DVHs) in the database for the planning target volume, bladder, and rectum. The lower left DVH set is for a selected single plan from the database; DVHs for all structures are shown. Notice that the seminal vesicles, with a hard minimum constraint of 60 Gy, are dosed properly because the underlying technology is linear programming, which models constraints exactly. A transversal dose wash for this plan is shown in the bottom right panel.
Our approach is similar to the approach detailed by Kuefer et al. (15). They begin with an equibalanced solution that is similar to our Step 1 solution. Two key differences at this point are their choice of objective functions—maxmean EUD (16) for healthy structures, and tumor control probability on the tumor minimum dose voxel—and their definition of an equibalanced solution, which omits the ($)+ operator from our Step 1 formulation. The approaches also differ slightly in how the Pareto database is populated, but the main difference relevant to the IMRT community is that their work focuses on technical aspects of the solution process, including an adaptive voxel clustering scheme developed to greatly reduce the solution time, whereas in the present work, we focus on demonstrating the clinical feasibility of an automated approach to IMRT Pareto database generation. Although we have provided a systematic procedure for producing Pareto plan databases, there is flexibility in how Steps 1 and 2 are performed. An obvious example is the choice of objective functions. Regarding our choice of a linear underdose ramp for the targets, our experience shows that this function produces plans that range from rounded to sharp regarding the shoulder of the target DVH curves, and that this
function is highly correlated with tumor control probability. Alternately, one could choose as a tumor objective to maximize the minimum tumor dose. It is also possible to include an objective which controls the complexity of the IMRT fluence map (17). Other choices regarding the implementation of the two-step approach include what to set as aspiration levels, how much—if at all—to soften constraints in Step 2, which structures to trade off, and what to use for maximum and minimum voxel dose constraints. In the spectrum of approaches to multiobjective treatment planing, we view the technique of constructing and navigating a Pareto database of solutions as the most difficult, but also the most valuable. Other multiobjective approaches such as prioritized optimization and goal programming (18, 19) lead to a single treatment plan and therefore do not convey tradeoff information, which we believe should be used in deciding on the final plan. For example, if a small increase in bladder NTCP can lead to a significant increase in tumor control probability, the physician and patient might be inclined to accept this additional risk. A Pareto database allows this tradeoff to be quickly assessed, as well as tradeoffs amongst any other relevant organs. Although only briefly discussed in this article, the subject
1606
I. J. Radiation Oncology d Biology d Physics
Volume 69, Number 5, 2007
Fig. 5. The top two rows show the extreme dose–volume histograms (DVHs) computed in the lung plan database. In a high-dimensional tradeoff though, organs can be traded off with each other and with planning target volume coverage, an exploration that is made possible with a slider-based navigator as depicted in Figure 5. The lower left DVH shows results for a single selected plan, which is also dose washed in the lower right panel.
of how to best navigate a Pareto database is an area of active research by the present authors. In the absence of a single number that reliably evaluates a treatment plan, the
Pareto surface methodology may best leverage a physician’s intuition to synthesize a sensible, patient-specific treatment plan.
REFERENCES 1. Ku¨fer K-H, Hamacher H, Bortfeld T. A multicriteria optimization approach for inverse radiotherapy planning. In: Bortfeld TR, Schlegel W, editors. Proceedings of the XIIIth ICCR. Heidelberg, Germany: Springer; 2000. p. 26–29. 2. Ehrgott M, Burjony M. Radiation therapy planning by multicriteria optimization. Proceedings of the 36th Annual Conference of the Operational Research Society of New Zealand. Auckland: Operational Research Society of New Zealand; 2001. p. 244–253.
3. Craft D, Halabi T, Bortfeld T. Exploration of tradeoffs in intensitymodulated radiotherapy. Phys Med Biol 2005;50:5857–5868. 4. Schreibmann E, Lahanas M, Xing L, et al. Multiobjective evolutionary optimization of the number of beams, their orientations, and weights for intensity-modulated radiation therapy. Phys Med Biol 2004;49:747–770. 5. Romeijn H, Dempsey J, Li J. A unifying framework for multicriteria fluence map optimization models. Phys Med Biol 2004;49:1991–2013.
An approach for practical multiobjective IMRT treatment planning d D. CRAFT et al.
6. Aubry J, Beaulieu F, Sevigny C, et al. Multiobjective optimization with a modified simulated annealing algorithm for external beam radiotherapy treatment planning. Med Phys 2006;33: 4718–4729. 7. Craft D, Halabi T, Shih H, et al. Approximating convex Pareto surfaces in multi-objective radiotherapy planning. Med Phys 2006;33:3399–3407. 8. Deasy J, Blanco A, Clark V. CERR: A computational environment for radiotherapy research. Med Phys 2003;30: 979–985. 9. Romeijn H, Ahuja R, Dempsey J, et al. A novel linear programming approach to fluence map optimization for intensity modulated radiation therapy treatment planning. Phys Med Biol 2003; 48:3521–3542. 10. Halabi T, Craft D, Bortfeld T. Dose-volume objectives in multi-criteria optimization. Phys Med Biol 2006;51: 3809–3818. 11. Niemierko A. A generalized concept of equivalent uniform dose. Med Phys 1999;26:1100. 12. Burman A, Kutcher GJ, Emami B, et al. Fitting of normal tissue tolerance data to an analytic function. Int J Radiat Oncol Biol Phys 1991;21:123–135.
1607
13. Wu Q, Djajaputra D, Liu H, et al. Dose sculpting with generalized equivalent uniform dose. Med Phys 2005;32:1387–1396. 14. De Neve W, Wu Y, Ezzell G. Image-guided IMRT. Berlin: Springer; 2006. 15. Ku¨fer K-H, Scherrer A, Monz M, et al. Intensity modulated radiotherapy—a large scale multi-criteria programming problem. OR Spectrum 2003;25:223–249. 16. Thieke C, Bortfeld TR, Ku¨fer K-H. Characterization of dose distributions through the max and mean dose concept. Acta Oncol 2002;41:158–161. 17. Craft D, Su¨ss P, Bortfeld T. The tradeoff between treatment plan quality and required number of monitor units in imrt. Int J Radiat Oncol Biol Phys 2007;67:1596–1605. 18. Vineberg KA, Jee KW, Ben-Josef E, et al. Multicriteria decision making for individualization of imrt treatment plans. Int J Radiat Oncol Biol Phys 2005;63(Suppl. 1):S503. 19. Wilkens J, Alaly J, Zakarian K, et al. Imrt treatment planning based on prioritizing prescription goals. Phys Med Biol 2007; 52:1675–1692. 20. Emami B, Lyman J, Brown A, et al. Tolerance of normal tissue to therapeutic irradiation. Int J Radiat Oncol Biol Phys 1991;21: 109–122.
doi:10.1016/j.ijrobp.2007.08.019
PHYSICS CONTRIBUTION
AN APPROACH FOR PRACTICAL MULTIOBJECTIVE IMRT TREATMENT PLANNING DAVID CRAFT, PH.D. TAREK HALABI, PH.D., HELEN A. SHIH, M.D., AND THOMAS BORTFELD, PH.D. Department of Radiation Oncology, Massachusetts General Hospital and Harvard Medical School, Boston, MA Purpose: To introduce and demonstrate a practical multiobjective treatment planning procedure for intensitymodulated radiation therapy (IMRT) planning. Methods and Materials: The creation of a database of Pareto optimal treatment plans proceeds in two steps. The first step solves an optimization problem that finds a single treatment plan which is close to a set of clinical aspirations. This plan provides an example of what is feasible, and is then used to determine mutually satisfiable hard constraints for the subsequent generation of the plan database. All optimizations are done using linear programming. Results: The two-step procedure is applied to a brain, a prostate, and a lung case. The plan databases created allow for the selection of a final treatment plan based on the observed tradeoffs between the various organs involved. Conclusions: The proposed method reduces the human iteration time common in IMRT treatment planning. Additionally, the database of plans, when properly viewed, allows the decision maker to make an informed final plan selection. Ó 2007 Elsevier Inc. Multiobjective, Inverse planning, Optimization, IMRT, Linear programming, EUD, Multicriteria, Pareto.
such plans forms the Pareto surface, and an MCO plan database consists of a sampling of Pareto optimal plans (7). With each structure (targets and organs at risk) in a plan, we associate an indicator function that maps the dose distribution in that structure to a scalar value. From the perspective of a treatment planner, after such functions are defined, the easiest procedure is to form a Pareto surface of plans with each of the indicator functions as objective functions. This is problematic in practice, because it results in a large Pareto surface with many plans that are not of clinical interest. To overcome this, we would like to constrain some or all of the indicator functions. This can lead to a human iteration loop because it is not, in general, possible to determine a set of mutually satisfiable and useful constraint levels for structures. In this article, we address these issues with a two-step process, in which the first step attempts to find such a set of constraints and the second step generates the Pareto surface based on these constraints.
INTRODUCTION Several authors have discussed the use of multiobjective optimization, also called multicriteria optimization (MCO), for intensity-modulated radiation therapy (IMRT) treatment planning (1–6). There are typically two motivations for MCO: to reduce the time-consuming human iteration loop required to come up with a satisfactory plan and to provide tradeoff information (in the form of a database of diverse treatment plans) that allows the decision makers to make an informed final plan choice. To simultaneously satisfy these two goals, one needs to determine a method of creating a database of good plans in a way that does not require timeconsuming human iterations. This is the topic we address in this article. Given that the most suitable plan for a particular patient is impossible to articulate before seeing what is possible regarding treatment plans, we prefer the multiobjective paradigm of a posteriori preference articulation, in which the doctor makes a decision after seeing the possible plans. This requires the computation of a database of treatment plans. We use the terms database of treatment plans and Pareto surface somewhat interchangeably. A treatment plan is Pareto optimal for a given set of objectives and constraints if there is no other feasible treatment plan that is at least as good in all the objectives and strictly better in at least one objective. The set of all
METHODS AND MATERIALS Radiation therapy planning computed tomography scans and physician-defined target and nontarget volumes are read into computational environment for radiotherapy research (CERR) (8) using CERR’s Radiation Therapy Oncology Group import facility. Within CERR, beam angles are selected and the dose influence matrix D This work was supported by NCI Grant 1 R01 CA103904-01A1: Multi-criteria IMRT Optimization. Conflict of interest: none. Received April 3, 2007, and in revised form July 31, 2007. Accepted for publication Aug 1, 2007.
Reprint requests to: David Craft, Ph.D., Department of Radiation Oncology, Harvard Medical School, Massachusetts General Hospital, 0 Grove St., Boston, MA 02114. Tel: (617) 724-3658; Fax: (617) 724-0368; E-mail: [email protected] 1600
An approach for practical multiobjective IMRT treatment planning d D. CRAFT et al.
that maps beamlet intensities x to the dose vector d is calculated. This information is made available to the optimization software, AMPL/CPLEX, by way of text files. For the optimizations, we use linear programming because of its ability to solve problems to provable optimality and because it handles hard constraints exactly. Any optimization involving linear functions and piecewise linear convex objectives and constraints can be written as a linear program, a feature that has been used in a variety of ways in IMRT optimization (3, 9, 10). We use piecewise linear convex functions for both the target and healthy structure indicator functions. For target modeling, we use a piecewise linear convex function that linearly penalizes doses below the prescription dose. For healthy structures, we fit the generalized equivalent uniform dose model (EUD) (11) with piecewise linear convex approximations. Generalized EUD for an organ is given by P EUD ¼ ð1=n i di a Þ1=a , where n is the number of voxels in the structure, di is the dose to voxel i, and a is a parameter. For MCO, we can drop the ath root in that formula, as explained elsewhere (5) (essentially, minimizing x subject to a set of constraints is the same as minimizing xa for any positive a). We fit the resulting dose summation by fitting the individual terms dia over the appropriate dose range, and summing the resulting fits. The functions dia are fit using three-segment piecewise linear functions (Fig. 1). Table 1 contains a values (12, 13) and the best fit (minimal area between curves) locations and relative weights of the linear sections used for each healthy structure modeled. We also employ hard constraints for maximum doses in healthy structures and minimum and maximum doses for targets. Step 1 of the Pareto database creation consists of setting aspiration levels for each indicator function and finding a treatment plan that is near to satisfying all these aspirations. Letting Fk denote the kth indicator function (Fk is the piecewise linear EUD for healthy structures and the linear underdose function for targets) and Ak denote the maximum desired level for this function, step 1 consists of solving the following optimization problem:
1601
Table 1. Values for the optimization function parameters Organ Brain case Brain stem Optic nerve Optic chiasm Eye Prostate case Rectum Bladder Femoral heads Lung case Lung Esophagus Spinal cord
EUD ‘‘a’’ parameter
Pivot 1
Pivot 2
Weight 2/ Weight 1
6.25 4 4 3.3
37.7 25.7 25.7 5.37
53.8 45.1 45.1 10.67
2.85 2.22 2.22 1.92
6.3 3.8 4
52.5 36.9 19.3
74.6 66.7 33.8
2.86 2.14 2.22
1.2 5 10
2.24 30.5 39
29 47.5 48.7
0.25 2.25 3.36
EUD a values are from previous work (12, 13). We fit the function f(d) = da with a three-segment piecewise linear convex function over the appropriate range for that particular structure. An example fit is shown in Fig. 1. Abbreviation: EUD = equivalent uniform dose.
minimize
P ðFk ðdÞ Ak Þþ k
(1)
subject to d ¼ Dx x$0:
The ($)+ operator is shorthand for max(0, $). Because this is a piecewise linear convex function, this can be written as a linear program. In other words, Formulation 1 tries to find a solution x, which gets every function Fk at or below its aspiration level Ak. Notice that this program is always feasible. We occasionally put additional constraints—for example, maximum dose constraints on all voxels and minimum dose constraints on the targets—at the risk of having an infeasible program. The aspiration levels can be chosen based on what experience suggests is a mutally satisfiable, although possibly aggressive set of goals. In this report, to avoid seemingly ad hoc aspiration levels, we use aspiration levels derived from setting the normal tissue complication probability (NTCP) to 0.1%, using data from Burman (12), and use an aspiration of 0 underdosing for the targets. The solution to the Step 1 optimization problem is used to determine a constraint set for Step 2, which is the generation of the Pareto database. For example, if the optimal value of the Step 1 problem is 0, it means all the aspirations were met and can be used as constraints in Step 2. On the other hand, for any aspirations that is not met, we know how far we are from meeting it based on examining the Step 1 solution. This can be used to form a new bound that is Table 2. Optimization information for the three cases
Fig. 1. Fitting the dose power term of the equivalent uniform dose function with a piecewise convex linear function. For comparison, we also plot the Emami (20) data in this graph, which shows good agreement.
Case
No. voxels
No. beamlets
Brain Prostate Lung
50,708 77,287 114,033
680 550 567
Voxel size dx, dy, dz (cm)
Beamlet size (cm)
0.2, 0.2, 0.25 0.3, 0.3, 0.25 0.3, 0.3, 0.25
0.5 0.5 1.0 1.0 1.0 1.0
Critical structures and unclassified tissue are downsampled in each case by a factor of 2 or 4. Targets are not downsampled. Number of voxels reflects the number of voxels used in the optimization (with downsampling), but all results (e.g., DVHs) are reported on the full voxel grid.
I. J. Radiation Oncology d Biology d Physics
1602
Volume 69, Number 5, 2007
Table 3. Brain case Brain optimization (PTV prescription: 59.4 Gy) Step 1
Step 2
Organ
EUD aspiration level (Gy)
Amount by which aspiration was violated in Step 1 (Gy)
EUD constraint for Step 2 (Gy)
Traded off against PTV?
EUD range in database (Gy)
Brain stem Optic chiasm L optic nerve R optic nerve L orbit R orbit Unclassified tissue
36 36 36 36 3 3 15
0 0 0 0 0 0 4
NA 36 NA 36 3 3 19
Yes No Yes No No No No
32–37.5 NA 19.5–27.5 NA NA NA NA
The aspirations for all of the healthy structures except for the unclassified tissue were met, and held tight at those achievable levels in Step 2, the creation of the Pareto database. Abbreviations: PTV = planning target volume; EUD = equivalent uniform dose. attainable. Softening attainable constraints further will yield a more diverse Pareto surface. It is up to the decision makers which functions to trade off in Step 2 and how to employ hard constraints, which can be used to keep any structure within reasonable limits. We use the algorithm PGEN (Pareto Generation) (7) to compute the Pareto database. The plan databases are navigated by means of a Pareto surface navigator, compatible with CERR, which is displayed in Figure 2 for the lung case. This navigator does linear plan interpolation, which allows one to smoothly transition between plans in the database. Exact EUDs are computed and displayed by the navigator. During the navigation process, a planner may constrain any of the functions to control the part of the surface explored. For each slider, the navigator also reports a resitance value, which is measure of the cost (i.e., the detriment to the other objectives) of pulling that slider down from the current position.
RESULTS For each case, we use seven equispaced beams. Beamlet and voxel discretization information is shown in Table 2. Sampling by a factor of 2 or 4 is used for some healthy structures to speed up the optimizations, but all results are reported
on the original grids. Individual plan optimizations take anywhere from 3 min (brain and prostate) to 20 min (lung) and the databases computed contain approximately 30 plans each. To set aspirations for the normal structures in Step 1, we use the EUD value that corresponds to NTCP of 0.1% (12). These levels are shown in Tables 3–5. Tables 3–5 show our choices for objectives and constraints for three clinical cases of a brain tumor, prostate cancer, and lung cancer, respectively. In general, with familiarity for achievable EUD levels on a case-by-case basis, a treatment planner might use a more customized approach, but this recipe approach works well for the cases herein. We use a hard maximum dose constraint of 115% of the prescription dose for every voxel. A hard minimum of approximately 90% prescription is used for the target voxels. Step 1 results provide us with a set of EUD constraints that can be simultaneously satisfied, but for the prostate and lung cases, to create a wide range of plans in the Pareto database, we soften these achievable values by 10%. In the case of the lungs, for even more exploration freedom we let each lung EUD go as high as 30 Gy. We display the
Table 4. Prostate case Prostate optimization (PTV prescription: 79 Gy, seminal vesicles: 60 Gy) Step 1
Step 2
Organ
EUD aspiration level (Gy)
Amount by which aspiration was violated in Step 1 (Gy)
EUD constraint for Step 2 (Gy)
Traded off against PTV?
EUD range in database (Gy)
Bladder Rectum L femoral head R femoral head Unclassified tissue
52 42 40 40 15
0 12.4 0 0 15
57 60 30 30 33
Yes Yes No No Yes
46–57 53–60 NA NA 27–33
The rectum overlaps the PTV and as expected the aggressive aspiration for the rectum was not met in Step 1. Step 2 produces a four dimensional Pareto surface (bladder, rectum, unclassified tissue, and PTV coverage). Abbreviations: PTV = planning target volume; EUD = equivalent uniform dose.
An approach for practical multiobjective IMRT treatment planning d D. CRAFT et al.
1603
Table 5. Lung case Lung optimization (prescription: 66 Gy) Step 1
Step 2
Organ
EUD aspiration level (Gy)
Amount by which aspiration was violated in Step1 (Gy)
EUD constraint for Step 2 (Gy)
Traded off against PTV?
EUD range in database (Gy)
Spinal cord Esophagus R lung L lung Unclassified tissue
30 44 10 10 15
11 0 0 0 10
45 48 48 30 27
Yes Yes Yes Yes Yes
40–45 18–38 18–38 1–10.5 21–27
Abbreviations as in Table 4. To create a database with a large diversity of plans, the lungs were allowed to go up to an EUD of 30 Gy, even though the aspiration of 10 Gy was met in Step 1. Dose escalation to the tumor is of paramount concern for lung patients: relaxing the lung constraint thus allows the physician to see the effects of increasing the target dose.
resulting diversity of the plans by plotting the extreme dose– volume histograms (DVHs) in Figures 3–5 and by reporting the associated range of EUD values in Tables 3–5. Below the DVH range figures, we have selected a single well-balanced plan from each of the databases. In the following section, we discuss the patient details and present the tradeoff results for three clinical cases. Detailed
Fig. 2. Screenshot of the Pareto surface navigator used to interactively slide through the plans in the database. The navigator is Matlab based and operates within CERR, where dose displays are updated automatically after the user slides a slider. The navigator is available from the authors or as a free download at the MathWorks.com website.
tradeoff analyses (which are most useful with state-of-theart NTCP and tumor control probability models) are beyond the scope of this work, but will be reported for specific disease sites in future publications. The first case was a 59-year-old man with a hemangiopericytoma of the brain involving the left middle cranial fossa. After presentation with left facial numbness, imaging revealed a 4-cm enhancing tumor along the skull base. The tumor abutted the left cavernous sinus, invaginated the medial temporal lobe, and eroded the left petrous apex. After subtotal resection, he was recommended adjuvant radiation therapy and was treated with IMRT. Prescribed tumor dose was 59.4 Gy. Figure 3 depicts the tradeoffs between the left optic nerve, the brainstem, and the planning target volume (PTV). The constraints for the other structures were easily met in Step 1, and not in difficult planning areas, so were not included in the tradeoff. The second case was 73-year-old man with an early stage prostate cancer diagnosed after detection of a palpable nodule on routine rectal exam. After biopsy confirmation of prostate cancer, he chose to undergo definitive radiation therapy. He was treated with IMRT to a dose of 79 Gy. Figure 4 demonstrates the anterior and posterior rectum displayed as separate DVHs, but they are optimized together. By modeling the unclassified tissue with an EUD parameter of a = 3.8 (the value for the bladder) and including it in the tradeoff, we are able to control general dose to the unclassified tissue, and in particular ‘‘streaking’’ of high-dose values far into normal tissues, because the relatively high a parameter makes the function sensitive to small volumes of high dose. A hard minimum of 70 Gy is applied to the PTV. The DVHs plotted by CERR show some voxels with slightly less than 70 Gy, but this is an artifact of the CERR DHV computation, and every PTV voxel gets 70 Gy or more in every plan. The third case was a 72-year-old woman with a non–smallcell lung cancer. She presented with a 6.5 cm right upper lobe tumor abutting the posterior chest wall and underwent resection of the tumor and postoperative chemotherapy. Intraoperatively, her tumor was adherent to the chest wall and the surgical margin was focally positive. Given her high risk of
1604
I. J. Radiation Oncology d Biology d Physics
Volume 69, Number 5, 2007
Fig. 3. Along the top row, extreme dose–volume histograms (DVHs) for each structure are displayed to demonstrate the range of plans available in the database. The lower left DVH set is for a selected single plan from the database, with DVHs for all structures shown. A transversal dose wash for this plan is shown in the bottom right.
local recurrence, she was recommended adjuvant radiation therapy and was treated to the tumor bed to 66 Gy using IMRT. The six-dimensional tradeoff for her case is shown in Figure 5. There is considerable tradeoff between dosing the two lungs: dose can be brought to the tumor through either lung almost entirely. Reducing the dose to the right lung, where the tumor occurs, creates a more conformal plan, but the high-dose band that comes through the left lung in this case might be undesireable. For the spinal cord, we plan on a 3-mm expansion of the cord. This cord planning volume overlaps with the PTV, and the minimum dose of 55 Gy to the PTV results in the adjacent section of the spinal cord receiving that dose as well. If this were deemed unacceptable, the planner would reduce the minimum required dose to the PTV. DISCUSSION AND CONCLUSION The goal of MCO for IMRT optimization is to reduce treatment planning time while making treatment planning more intuitive. This is achieved by constructing a database of Pareto optimal plans that the treatment planner can navigate through interactively, thereby acquiring a feel for the tradeoffs involved before selecting a final plan. This navigation
depends on mixing plans in the database in real time, and is more manageable than the Pareto visualization aid presented by Aubry et al. (6), which displays and color codes the plans in the database but does not allow plan combinations. The primary challenge addressed by this work is how to efficiently create a database of Pareto optimal treatment plans that is both diverse and populated with clinically useful plans. These criteria are in conflict with another, because too much diversity means many plans that are far away from a plan a physician would give. This challenge has been addressed by breaking down the task into two steps. The first step establishes a balanced plan that attempts to satisfy a reasonable (meaning possibly achievable) set of constraints. The solution to this problem yields a set of simultaneously achievable constraints. The user can then set true dose constraints rather than provisional constraints often used in IMRT planning (14). The achievable constraints from Step 1 are then used in the Pareto database generation, Step 2. Without Step 1, one might either specify a set of constraints that leads to an infeasible problem, and have to iterate to remedy this, or specify a set of constraints that is too loose, resulting in a database with large percentage of plans outside the domain of what might be delivered. The technical details of Step 2 can be found elsewhere (7).
An approach for practical multiobjective IMRT treatment planning d D. CRAFT et al.
1605
Fig. 4. The top row shows the extreme dose–volume histograms (DVHs) in the database for the planning target volume, bladder, and rectum. The lower left DVH set is for a selected single plan from the database; DVHs for all structures are shown. Notice that the seminal vesicles, with a hard minimum constraint of 60 Gy, are dosed properly because the underlying technology is linear programming, which models constraints exactly. A transversal dose wash for this plan is shown in the bottom right panel.
Our approach is similar to the approach detailed by Kuefer et al. (15). They begin with an equibalanced solution that is similar to our Step 1 solution. Two key differences at this point are their choice of objective functions—maxmean EUD (16) for healthy structures, and tumor control probability on the tumor minimum dose voxel—and their definition of an equibalanced solution, which omits the ($)+ operator from our Step 1 formulation. The approaches also differ slightly in how the Pareto database is populated, but the main difference relevant to the IMRT community is that their work focuses on technical aspects of the solution process, including an adaptive voxel clustering scheme developed to greatly reduce the solution time, whereas in the present work, we focus on demonstrating the clinical feasibility of an automated approach to IMRT Pareto database generation. Although we have provided a systematic procedure for producing Pareto plan databases, there is flexibility in how Steps 1 and 2 are performed. An obvious example is the choice of objective functions. Regarding our choice of a linear underdose ramp for the targets, our experience shows that this function produces plans that range from rounded to sharp regarding the shoulder of the target DVH curves, and that this
function is highly correlated with tumor control probability. Alternately, one could choose as a tumor objective to maximize the minimum tumor dose. It is also possible to include an objective which controls the complexity of the IMRT fluence map (17). Other choices regarding the implementation of the two-step approach include what to set as aspiration levels, how much—if at all—to soften constraints in Step 2, which structures to trade off, and what to use for maximum and minimum voxel dose constraints. In the spectrum of approaches to multiobjective treatment planing, we view the technique of constructing and navigating a Pareto database of solutions as the most difficult, but also the most valuable. Other multiobjective approaches such as prioritized optimization and goal programming (18, 19) lead to a single treatment plan and therefore do not convey tradeoff information, which we believe should be used in deciding on the final plan. For example, if a small increase in bladder NTCP can lead to a significant increase in tumor control probability, the physician and patient might be inclined to accept this additional risk. A Pareto database allows this tradeoff to be quickly assessed, as well as tradeoffs amongst any other relevant organs. Although only briefly discussed in this article, the subject
1606
I. J. Radiation Oncology d Biology d Physics
Volume 69, Number 5, 2007
Fig. 5. The top two rows show the extreme dose–volume histograms (DVHs) computed in the lung plan database. In a high-dimensional tradeoff though, organs can be traded off with each other and with planning target volume coverage, an exploration that is made possible with a slider-based navigator as depicted in Figure 5. The lower left DVH shows results for a single selected plan, which is also dose washed in the lower right panel.
of how to best navigate a Pareto database is an area of active research by the present authors. In the absence of a single number that reliably evaluates a treatment plan, the
Pareto surface methodology may best leverage a physician’s intuition to synthesize a sensible, patient-specific treatment plan.
REFERENCES 1. Ku¨fer K-H, Hamacher H, Bortfeld T. A multicriteria optimization approach for inverse radiotherapy planning. In: Bortfeld TR, Schlegel W, editors. Proceedings of the XIIIth ICCR. Heidelberg, Germany: Springer; 2000. p. 26–29. 2. Ehrgott M, Burjony M. Radiation therapy planning by multicriteria optimization. Proceedings of the 36th Annual Conference of the Operational Research Society of New Zealand. Auckland: Operational Research Society of New Zealand; 2001. p. 244–253.
3. Craft D, Halabi T, Bortfeld T. Exploration of tradeoffs in intensitymodulated radiotherapy. Phys Med Biol 2005;50:5857–5868. 4. Schreibmann E, Lahanas M, Xing L, et al. Multiobjective evolutionary optimization of the number of beams, their orientations, and weights for intensity-modulated radiation therapy. Phys Med Biol 2004;49:747–770. 5. Romeijn H, Dempsey J, Li J. A unifying framework for multicriteria fluence map optimization models. Phys Med Biol 2004;49:1991–2013.
An approach for practical multiobjective IMRT treatment planning d D. CRAFT et al.
6. Aubry J, Beaulieu F, Sevigny C, et al. Multiobjective optimization with a modified simulated annealing algorithm for external beam radiotherapy treatment planning. Med Phys 2006;33: 4718–4729. 7. Craft D, Halabi T, Shih H, et al. Approximating convex Pareto surfaces in multi-objective radiotherapy planning. Med Phys 2006;33:3399–3407. 8. Deasy J, Blanco A, Clark V. CERR: A computational environment for radiotherapy research. Med Phys 2003;30: 979–985. 9. Romeijn H, Ahuja R, Dempsey J, et al. A novel linear programming approach to fluence map optimization for intensity modulated radiation therapy treatment planning. Phys Med Biol 2003; 48:3521–3542. 10. Halabi T, Craft D, Bortfeld T. Dose-volume objectives in multi-criteria optimization. Phys Med Biol 2006;51: 3809–3818. 11. Niemierko A. A generalized concept of equivalent uniform dose. Med Phys 1999;26:1100. 12. Burman A, Kutcher GJ, Emami B, et al. Fitting of normal tissue tolerance data to an analytic function. Int J Radiat Oncol Biol Phys 1991;21:123–135.
1607
13. Wu Q, Djajaputra D, Liu H, et al. Dose sculpting with generalized equivalent uniform dose. Med Phys 2005;32:1387–1396. 14. De Neve W, Wu Y, Ezzell G. Image-guided IMRT. Berlin: Springer; 2006. 15. Ku¨fer K-H, Scherrer A, Monz M, et al. Intensity modulated radiotherapy—a large scale multi-criteria programming problem. OR Spectrum 2003;25:223–249. 16. Thieke C, Bortfeld TR, Ku¨fer K-H. Characterization of dose distributions through the max and mean dose concept. Acta Oncol 2002;41:158–161. 17. Craft D, Su¨ss P, Bortfeld T. The tradeoff between treatment plan quality and required number of monitor units in imrt. Int J Radiat Oncol Biol Phys 2007;67:1596–1605. 18. Vineberg KA, Jee KW, Ben-Josef E, et al. Multicriteria decision making for individualization of imrt treatment plans. Int J Radiat Oncol Biol Phys 2005;63(Suppl. 1):S503. 19. Wilkens J, Alaly J, Zakarian K, et al. Imrt treatment planning based on prioritizing prescription goals. Phys Med Biol 2007; 52:1675–1692. 20. Emami B, Lyman J, Brown A, et al. Tolerance of normal tissue to therapeutic irradiation. Int J Radiat Oncol Biol Phys 1991;21: 109–122.