- Email: [email protected]
Interactively exploring optimized treatment plans
Int. J. Radiation Oncology Biol. Phys., Vol. 61, No. 2, pp. 570 –582, 2005 Copyright © 2005 Elsevier Inc. Printed in the USA. All rights reserved 0360-3016/05/$–see front matter
doi:10.1016/j.ijrobp.2004.09.022
PHYSICS CONTRIBUTION
INTERACTIVELY EXPLORING OPTIMIZED TREATMENT PLANS ISAAC ROSEN, PH.D.,* H. HELEN LIU, PH.D.,* NATHAN CHILDRESS, PH.D.,* AND ZHONGXING LIAO, M.D.† *Departments of Radiation Physics and †Radiation Oncology, The University of Texas M. D. Anderson Cancer Center, Houston, TX Purpose: A new paradigm for treatment planning is proposed that embodies the concept of interactively exploring the space of optimized plans. In this approach, treatment planning ignores the details of individual plans and instead presents the physician with clinical summaries of sets of solutions to well-defined clinical goals in which every solution has been optimized in advance by computer algorithms. Methods and Materials: Before interactive planning, sets of optimized plans are created for a variety of treatment delivery options and critical structure dose–volume constraints. Then, the dose–volume parameters of the optimized plans are fit to linear functions. These linear functions are used to show in real time how the target dose–volume histogram (DVH) changes as the DVHs of the critical structures are changed interactively. A bitmap of the space of optimized plans is used to restrict the feasible solutions. The physician selects the critical structure dose–volume constraints that give the desired dose to the planning target volume (PTV) and then those constraints are used to create the corresponding optimized plan. Results: The method is demonstrated using prototype software, Treatment Plan Explorer (TPEx), and a clinical example of a patient with a tumor in the right lung. For this example, the delivery options included 4 open beams, 12 open beams, 4 wedged beams, and 12 wedged beams. Beam directions and relative weights were optimized for a range of critical structure dose–volume constraints for the lungs and esophagus. Cord dose was restricted to 45 Gy. Using the interactive interface, the physician explored how the tumor dose changed as critical structure dose–volume constraints were tightened or relaxed and selected the best compromise for each delivery option. The corresponding treatment plans were calculated and compared with the linear parameterization presented to the physician in TPEx. The linear fits were best for the maximum PTV dose and worst for the minimum PTV dose. Based on the root-mean-square error between the fit values and their corresponding data values, the linear fit appears to be adequate, although higher order polynomials could give better results. Some of the variance in fit is due to the stochastic nature of the simulated annealing optimization algorithm, which does not reproduce the exact same results in repetitions of the same calculation. Using a directed search algorithm for plan optimization should produce better parameter fits and, therefore, better predictions of plan characteristics by TPEx. Conclusions: Using TPEx, the physician can easily select the optimum plan for a patient, with no imposed arbitrary definition of the “best” plan. More importantly, the physician can readily see what can be achieved for the patient with a given delivery technique. There is no more uncertainty about whether or not a better plan exists. By comparing the “best” plans for different delivery options (e.g., three-dimensional conformal radiotherapy versus intensity-modulated radiation therapy), the physician can gauge the clinical benefits of greater technical complexity. However, before the TPEx process can be clinical useful, faster computers and/or algorithms are needed and more studies are needed to better model the spaces of optimized solutions. © 2005 Elsevier Inc. Radiation therapy treatment planning, Inverse planning, Treatment plan optimization.
INTRODUCTION The goal of treatment planning, whether forward-based or inverse-based, is to find the best plan for the individual patient. A major difficulty in accomplishing that goal is the lack of universally accepted criteria for defining the “best” plan. Clearly, more dose to the tumor is better and less dose to the normal tissues is better. However, the ideal plan of
100% prescribed dose to the entire tumor volume and zero dose elsewhere is not attainable. In evaluating any given treatment plan, the radiation oncologist is faced with a variety of issues and choices. The primary question is whether a dose can be delivered to the tumor that will achieve a high probability of local control without serious treatment morbidity or sequelae. If that goal can be achieved, then the focus shifts to whether treatment
Reprint requests to: Isaac Rosen, Ph.D., Department of Radiation Physics, Unit 94, The University of Texas M. D. Anderson Cancer Center, 1515 Holcombe Blvd, Houston, TX 77030; Tel: (713) 5632635; Fax: (713) 563-2629; E-mail: [email protected]
Support in part by a sponsored research agreement from Philips Medical Systems. Received Jun 8, 2004, and in revised form Sep 15, 2004. Accepted for publication Sep 17, 2004. 570
Exploring optimized plans
morbidity can be further reduced (i.e., can normal tissue doses be lowered?). When there is a desire to escalate tumor doses because local control is not adequate, the question becomes one of finding the best compromise between potential complications and potential recurrence. The answers to these questions depend on the details of treatment delivery. In general, more complex treatments with more degrees of freedom in beam geometry and fluence produce better plans with greater separation between tumor dose and normal tissue doses. In designing the best plan for a patient, the physician must often compromise between the prescribed tumor dose, the uniformity of dose in the tumor volume, and the doses to critical structures. Further compromises may be required between the different adverse reactions of different critical structures. Unfortunately, the probabilities of complications for an individual patient are not binary functions and not currently predictable. The physician must try to anticipate how the patient will respond to the planned dosage and may need to compromise tumor dose to manage the toxicity. Because different physicians may choose different compromises based on their training and clinical experience and their understanding of the patient’s unique situation, the results of treatment planning vary with physician and with patient. If it were possible to mathematically define the characteristics of the best plan, then computer programs could unerringly find it for each patient. Attempts to accomplish this go back more than 35 years, and there is an extensive body of literature exploring the use of algorithms such as exhaustive search, linear programming, quadratic programming, mixed-integer programming, gradient search, simulated annealing, genetic algorithms, feasibility search, downhill simplex, and neural networks (e.g., 1–22). In theory, inverse planning (or treatment plan optimization, as it used to be called) could provide the unequivocal best treatment plan for a patient if the biologic responses of the individual patient’s tissues could be accurately modeled and the clinical compromises acceptable to the clinician were objectively and unambiguously known. Although mathematical models of tumor and tissue response have been developed, they are not yet accurate enough for prediction of individual patient responses and clinician decision criteria are not well simulated by computer algorithms. Consequently, inverse planning is still an iterative process. A major improvement of inverse planning over forward planning is the focus on clinical factors, such as normal tissue tolerances, rather than on the details of beam geometries. The existence of multiple competing objectives in the treatment planning problem has led to approaches that focus on the relationships between these objectives rather than attempting to define a mathematically “best” solution (e.g., 23–28). In general, these methods compute multiple plans with different optimization criteria to find the boundary of the space of optimized (Pareto) or feasible solutions. They may use sensitivity analysis to show how much tumor dose can be increased for a given increase in normal tissue dose
● I. ROSEN et al.
571
or may directly optimize the weighting factors assigned to the different objectives. The new treatment planning paradigm presented here follows this general approach and extends it by using the optimization criteria directly, rather than through surrogate preference factors, and by providing a real-time interactive interface into the space of optimized plans. It is worthwhile to ask why treatment planning is an iterative process. The clinician presumably knows what the ideal (or acceptable) dose distribution should look like. The physicist or dosimetrist generates a plan that he or she believes is the closest approximation to what the clinician’s ideal plan is. Why should there be any iterations? One reason may be that the planner has not correctly anticipated the physician’s goals and compromise preferences. For example, the planner might have chosen to spare equal portions of each kidney, whereas the physician’s preference might be to sacrifice one kidney and completely spare the other. Another reason for multiple iterations may be that the physician implicitly believes that a better plan exists and that through more effort (and perhaps experience) it can be found. In either case, the physician is faced with a process in which the characteristics of the best achievable plan are unknown. The dilemma facing the physician can be described with a simple example. Consider a patient with a lung tumor. The planner generates a conventional beam arrangement using large anteroposterior beams for an initial dose followed by a boost dose delivered with parallel-opposed oblique beams reduced in size to the primary target volume. The physician must balance the doses to the tumor, lungs, heart, esophagus, and spinal cord. Changing the angles of the oblique beams will change this balance, as will adding more beams or fluence modulation. Presumably, the planner has already optimized the directions of the oblique beams either through trial and error or previous experience. However, based on the single plan presented, the physician cannot know how much improvement can be achieved by adding more beams or adding fluence modulation. Therefore, the physician may ask for additional plans to assess the benefits of these changes versus the added complexity of treatment delivery. As fluence modulation or beams are added, the optimality of the beam directions is thrown into question. The physician no longer knows how good a plan can be achieved or what level of complexity is needed to achieve the clinical goals. The information that the physician needs to make the best possible treatment planning decision for a given patient can be encapsulated in the question, “What is the greatest tumor dose that can be delivered for specified normal tissue dose– volume constraints, regardless of the beam characteristics?” Given this information in an interactive format, the physician can readily select the plan that delivers the most tumor dose with acceptable normal tissue risk or the one that most reduces the dose to selected organs while maintaining a desired tumor dose level. The answer to this question requires that the best possible (optimized) dose–volume histogram (DVH) for the planning target volume (PTV) be
572
I. J. Radiation Oncology
● Biology ● Physics
presented to the physician for every combination of normal tissue DVHs that they might want to select. We are studying a new paradigm for treatment planning that embodies this concept of interactively exploring the space of optimized plans. As with some other inverse planning techniques, the details of treatment delivery (such as beam weights, angles, and doses) are hidden. Instead, the physician is presented with the relevant clinical summaries for the best possible solutions to well-defined clinical goals in which every solution has been optimized in advance by computer algorithms. Using an interactive interface, the physician selects a general treatment delivery method and then explores in real time how the maximum deliverable tumor dose changes as critical structure dose–volume constraints are tightened or relaxed. The physician can then directly select the compromises between doses to the various organs, dose to the tumor, and complexity of treatment delivery that he or she desires for the patient with no imposed arbitrary definition of the “best” plan. More importantly, the physician can see what is achievable for the patient. There is no more uncertainty about whether or not a better plan exists.
Volume 61, Number 2, 2005
Selection of plan parameters One goal of TPEx is to avoid dealing with the details of beam geometry. The first user screen (Fig. 2) asks for only three parameters to define the treatment delivery—the beam energy, the number of beams, and the degree of fluence modulation. Based on these selections, the program reads a set of previously computed optimized plans. It then presents the user with choices of critical structures to manipulate during exploration and the corresponding available decision doses. The decision doses are dose levels at which the physician would like to control the irradiated volume. Typically, these would be dose levels at which the irradiated volume is predictive of complications. The choices available to the user are limited to those for which optimized sets of plans were previously generated. Therefore, the planner should consult with the physician before beginning the optimization phase or have sufficient experience to anticipate the physician’s general approach to treatment for the given anatomic site. The planner should know what decision doses the physician uses in evaluating plans, what beam energies are preferred, and how much treatment complexity will probably be needed to accomplish the clinical goals. In that way, the planner can choose the best set of plans to be optimized in preparation for treatment plan exploration.
Interactive exploration METHODS AND MATERIALS A prototype system named Treatment Plan Explorer (TPEx) was developed using MATLAB (The MathWorks Inc., Natick, MA) to test this approach to treatment planning. It imports a set of optimized treatment plans, parameterizes their dose–volume characteristics, and then presents the results to the physician or planner in an interactive format. Before using TPEx, multiple sets of optimized plans are created, each of which covers the range of potentially desirable normal tissue DVH constraints. Each set corresponds to a different treatment delivery technique that the user may wish to explore. TPEx does not include external beam radiation treatment planning (RTP) software. The optimized plans used by TPEx are created by the RTP software. The output of TPEx is a set of critical structure DVH constraints selected by the planner, which are then used by the RTP software to create the final optimized plan for the patient. Treatment Plan Explorer has two user interface screens. In the first, the planner chooses the basic plan parameters and the critical structures to be considered. In the second, the planner interactively adjusts the dose–volume constraints on the critical structures and observes the effects on the DVH of the PTV. During exploration, TPEx continuously displays the values of the user parameters and the resulting PTV parameters. The use of TPEx is demonstrated using a clinical example. Figure 1 shows three sections through a patient with a tumor in the right lung treated using three coplanar 6-MV beams. The anterior beam included a 30° wedge and the posterior beam included a 45° wedge. The right posterior-oblique beam included a 15° wedge and partial blocking. The dose distribution delivered 63 Gy to the 95% isodose line. The DVH of the treatment is also shown. The maximum dose to the spinal cord was 45.5 Gy; the volume of right lung exceeding 20 Gy was 59%; the volume of left lung exceeding 20 Gy was 19%; and the volume of esophagus exceeding 50 Gy was 21%. Because most of the PTV was superior to the heart, only a small volume of heart (10%) received more than 45 Gy from the coplanar beam arrangement.
The interactive exploration of optimized solutions is shown in Fig. 3. The upper graph is a simplified DVH for the PTV showing minimum dose, 95% volume dose, and maximum dose. The bottom graphs are simple DVHs for the critical structures selected in the previous input screen. Each critical structure DVH is defined by one or two control points that are individually movable by the user. The maximum dose control points can only move laterally, changing dose but not volume. The decision dose control points can only move vertically, changing volume but not decision dose value. For each point, a red bar shows the range of permissible motion. As the planner changes the value of one of the control points (the “active” control point), two things happen in real time. First, the DVH of the PTV changes to reflect the best plan that corresponds to the current settings of all the control points. Second, the red bars showing the ranges of motion of the control points change dynamically, including the range for the active control point. The active control point is not allowed to move outside of its range. Because of the complexity of the space of optimized plans, there are strong interactions between the control points and there may be discontinuities in the space. It is possible that the range of motion of a nonactive control point may change in such a way that the corresponding control point is no longer within its range. That situation represents a nonfeasible solution and should not be selected. When that control point is activated, it is automatically relocated to a value within its range. If a situation arises in which the maximum dose control point is reduced to the decision dose value, then the volume corresponding to the decision dose is automatically reduced to 0% because the decision dose is now also the maximum dose. The planner should explore all of the control points in various ways and sequences to find the set of DVHs (PTV and critical structures) that is most desirable. Just as in conventional treatment planning, lower doses to the critical structures result in lower doses to the PTV. However, by observing the changes in the PTV DVH with changes in the control points, it is easy to see how much each critical structure limits the PTV dose and how much sparing of
Exploring optimized plans
● I. ROSEN et al.
573
Fig. 1. Transverse, sagittal, and coronal sections through a patient with a tumor in the right lung. The patient was treated with three coplanar 6-MV beams. The anterior beam included a 30° wedge and the posterior beam included a 45° degree wedge. The right posterior-oblique beam included a 15° wedge and partial blocking. The dose distribution delivered 63 Gy to the 95% isodose line.
each critical structure can be achieved for a given PTV dose. In addition, the simultaneous changes in the magnitude of the prescription dose and the degree of dose homogeneity in the PTV can be observed in response to alterations in critical structure constraint values. After a set of control point values is chosen, those critical structure constraints are used as input to the RTP software to create an optimized plan with the beam characteristics (energy, number, modulation) chosen by the user at the start of the process. Because the user is interrogating a parameterization of the original data set of optimized plans and also interpolating, the actual plan achieved will not exactly match the results of TPEx. How closely the results of TPEx will match the results of RTP software optimization with the desired control point constraints depends on several factors— the density of data in the space of optimized plans, the optimization algorithm, and the quality of the function fitting. These relationships have not yet been explored.
Creation of optimized plan sets The patient treatment was originally planned using Pinnacle3 (Philips Radiation Oncology Systems, Madison, WI). Normal tissue and tumor volume structures were extracted from Pinnacle3 and imported into a research planning system, Plan3D, which uses a modified Clarkson scatter-air ratio dose calculation (29) and simulated annealing for plan optimization (14). All cal-
Fig. 2. The first Treatment Plan Explorer screen asks for only three parameters to define the treatment delivery—the beam energy, the number of beams, and the degree of fluence modulation. Based on these selections, the program reads a set of previously optimized plans. It then presents the user with choices of critical structures to manipulate during exploration and the corresponding available decision doses.
574
I. J. Radiation Oncology
● Biology ● Physics
Volume 61, Number 2, 2005
Fig. 3. The interactive exploration of all the possible optimized solutions is shown. The upper graph is a simplified dose–volume histogram (DVH) for the planning target volume showing minimum dose, 95% volume dose, and maximum dose. The bottom graphs are simple DVHs for the critical structures selected in the previous input screen. Each critical structure DVH is defined by one or two control points that are individually movable by the user. The maximum dose control points can only move laterally, changing dose but not volume. The decision dose control points can only move vertically, changing volume but not decision dose value. For each point, a red bar shows the range of permissible motion.
culations included the full three-dimensional anatomy and bulk tissue heterogeneity effects (Batho power-law correction). DVHs were computed using quasi-random points rather than a Cartesian grid. The number of points used per structure was 1000 for the PTV, 600 for each lung, 100 for the esophagus, and 100 for the spinal cord. Eight optimized data sets using 6-MV photons were created, as shown in Table 1. For each type of modulation, open and wedge, plan sets were created with 4, 6, 8, and 12 coplanar beams. In each case, the relative beam weights and their directions were optimized. Beam direction optimization was accomplished by starting with 36 eligible beams equally spaced about the patient at 10° angular increments. Beam weights were initially optimized for all 36 beams. Then the number of eligible beams was reduced by eliminating those with low contributions to the plan. This process was repeated up to four times until the desired number of beams remained. For wedge plans, the wedge angles and directions were also optimized. Beam shapes were automatically generated to conform to the PTV projection with 0.5-cm margins, similar to the margins used in the clinical plan. However, unlike the clinical plan, no partial beam blocking was included. For each delivery method (modulation and beam number), a set of optimized plans was automatically generated using the various com-
binations of dose–volume constraints for the critical structures shown in Table 1. The optimization objective was to maximize the minimum PTV dose subject to strictly enforced dose–volume constraints (“hard” constraints). For each decision dose, an optimized plan was created using each volume limit. For this example case, there were seven different volume levels permitted to exceed 20 Gy for the right lung, seven different volume levels permitted to exceed 20 Gy for the left lung, and six different volume levels permitted to exceed 50 Gy for the esophagus. The spinal cord was limited to 45 Gy maximum dose, which was not varied. An optimized plan was generated for each combination of dose–volume constraints. Therefore, there were 294 plans in each set (7 ⫻ 7 ⫻ 6). No constraints were applied to the PTV so that there was always a feasible solution. The actual volumes exceeding the decision dose values were extracted from the DVHs for all the structures in each plan and recorded. Because of the optimization method, these values were always less than or equal to the constraint values and represented what could theoretically be achieved by the plan. The maximum dose to each critical structure was also recorded. These critical structure data were the control points for TPEx. For the PTV, the minimum dose, maximum dose, and prescription dose (95% volume) were recorded for each plan. An example of the output from one optimization is shown in Table 2. The constraint values used
Exploring optimized plans
● I. ROSEN et al.
575
Table 1. Parameter values for the sets of optimized plans Modulation
Number of beams
Critical structure
Maximum dose (Gy)
Decision dose (Gy)
Volume limit for decision dose (%)
4 6 8 12
Right lung Left lung Esophagus Spinal cord
75 75 75 45
20 20 50 45
20, 30, 40, 50, 60, 70, 80 10, 20, 30, 40, 50, 60, 70 10, 20, 30, 40, 50, 60 0
Open Wedge
For each decision dose, an optimized plan was created using each volume limit shown.
the number of independent variables is twice the number of critical structures.
as input to the optimization process and the resulting beam parameters (directions and doses) that produced the optimum plan were not used in TPEx.
Boundaries of the space of optimized plans
Parameterization of optimized plans
It is important that the user’s exploration be limited to the space of optimized solutions. Because of the parameterization, it is possible to compute PTV values for any set of control point values. However, not all combinations of control point values correspond to feasible solutions. Therefore, the motions of control points must be restricted by the boundaries of the space of optimized plans. This boundary will generally be complex and nonanalytic. The boundaries for any given independent variable (control point) will vary with the specific values of the other independent variables. The space of optimized plans was quantized using a bitmap of dimension Rm, where R is the desired resolution of the bitmap and m is the number of independent variables. Each dimension of the bitmap corresponds to one of the independent variables. The scale associated with dimension i was computed by
The PTV characteristics of the optimized plans were fit to linear functions of the control points, PTVmin ⫽ a0min ⫹
n
兺a
VDD i ⫹
min i
i⫽1
PTV95% ⫽ a095% ⫹
VDD i ⫹
95% i
i⫽1
PTVmax ⫽ a0max ⫹
n
兺a
min i
Dmax i
(1)
i⫽1
n
兺a
n
兺b n
兺b
95% i
Dmax i
(2)
max i
Dmax i
(3)
i⫽1
VDD i ⫹
max i
i⫽1
n
兺b i⫽1
where
⌬xi ⫽ (xmax ⫺ xmin) ⁄ (R ⫺ 1)
PTVmin ⫽ minimum dose in PTV PTV95% ⫽ minimum dose in 95% of PTV PTVmax ⫽ maximum dose in PTV VDD ⫽ volume exceeding the decision dose value for critical i organ i Dmax ⫽ maximum dose in critical organ i i n ⫽ number of critical organs ai, bi ⫽ fit parameters The control points (VDD and Dmax i i ) are the independent variables and the PTV doses are the dependent variables. For this example,
(4)
where xmax and xmin are the maximum and minimum values, respectively, of the independent variable xi, ⌬xi is the increment in values along the xi axis, and R is the dimension of the bitmap along each axis. Applying Eq. 4 to each independent variable produces a set of control point values for each point in the bitmap of the space of optimized plans. Note that the bitmap can be very large. For this example with eight independent variables, a resolution of eight values per variable was used, giving a bitmap with 16,777,216 elements.
Table 2. Sample output from the optimization process Optimization input
Critical structure Right lung Left lung Esophagus Spinal cord PTV minimum dose (100% volume) PTV 95% volume dose PTV maximum dose (0% volume)
Decision dose (Gy)
Volume limit for decision dose (%)
20 20 50 45
20 10 5 0
Optimization output Plan volume exceeding decision dose (%) 20.0 3.0 0.0 0.0
Maximum dose (Gy) 27.6 33.3 30.7 30.6
18.0 Gy 19.0 Gy 33.7 Gy
Abbreviations: PTV ⫽ planning target volume; TPEx ⫽ treatment plan explorer. The decision dose values, plan volumes exceeding the decision dose, maximum dose, and PTV doses were used for interactive exploration (highlighted). The constraint volume limits for the decision doses were used for optimization but not in TPEx.
576
I. J. Radiation Oncology
● Biology ● Physics
Volume 61, Number 2, 2005
Fig. 4. A sequence of steps that leads to the best plan selected for 12 wedged beams. For each plan in the optimized set, the corresponding point in the bitmap is set to a value of 1. Consequently, the boundary of the space of optimized plans is defined by the edges of the volume of 1s in the bitmap. For any set of values of m-1 variables, the range of acceptable values for the mth variable can easily be found by scanning the mth axis of the bitmap. The values corresponding to the limits of the bits with value of 1 are the range for the variable. Depending on the density of data in the space of optimized plans, there may be gaps in the bitmap from undersampling. Therefore, the extreme limits of the 1s
along an axis are taken as the range of the variable, even though the 1s may not be continuous.
Initial solution The interactive exploration can start with any plan from the optimized set. Ideally, there should be no isolated regions in the space of optimized plans, so that it should be possible to reach any solution from any starting point. However, a good starting solution can be beneficial for two reasons. First, finite sampling
Exploring optimized plans
of the space might, in fact, produce isolated regions that could trap the user and prevent reaching the most desirable one. Second, there is no reason not to expect discontinuities in the space of optimized solutions. In suggesting that there are discontinuities, we mean that every point in the space is not linearly connected to every other point. Although it should be possible to reach every point in the space from every other point, the sequences of moves along the multiple axes of the multidimensional space may not be commutative. Finding a particular solution will depend on the sequence of steps as well as the starting point. Therefore, a good starting solution will reduce the amount of interactive searching needed to find the most desirable one. We chose as our initial solution the one with the largest PTV maximum dose.
RESULTS One of us (Z.L.) used TPEx to select the best plans for treatment configurations using 4 open beams, 12 open beams, 4 wedged beams, and 12 wedged beams. For each beam configuration, the control point values (decision dose volumes and maximum dose values) were interactively adjusted until a solution was found that gave the best compromise in doses to the critical structures and PTV. Figure 4 shows a sequence of steps that led to the best plan selected for 12 wedged beams. The initial solution is shown in Fig. 4a. First, the volume of left lung exceeding 20 Gy (V20) was reduced from 43.3% to 26.7% (the bottom of the available range) (Fig. 4b). Then, the maximum dose to the esophagus was reduced from 75 Gy to 68 Gy (Fig. 4c). This reduction changed the range of available values for the decision dose for the left lung, allowing the V20 to be further reduced to 9.4% (Fig. 4d). Next, reducing the maximum dose to the right lung from 75 Gy to 71.8 Gy improved the uniformity of dose in the PTV by reducing the maximum PTV dose (Fig. 4e). Finally, increasing the maximum esophagus dose from 68 Gy to 70.6 Gy further improved the PTV dose uniformity by increasing the minimum dose (Fig. 4f). The first two steps (Figs. 4b and 4c) show that reducing the maximum dose to the esophagus lowers the range of possible doses to V20 of the left lung. It seems reasonable that the lower dose levels for V20 of the left lung should have been available in the initial solution. We attribute this deficiency to the coarse sampling of the space of optimized plans and the coarse quantization of the space by the bitmap. Both of these are consequences of limited computing resources. The actual sequence of steps originally used to produce this final plan was not recorded. The results of the linear fits (Eqs. 1–3) to the control point values are shown in Fig. 5 for each of the delivery techniques. The fit value is plotted against the corresponding plan value for each optimized plan and the root-mean-square error is shown. The trend lines are least-squares fits to the data. The fits were best for the maximum PTV dose and worst for the minimum PTV dose. In general, the linear model seems to be appropriate and higher order polynomials have not yet been tested. The critical structure DVH constraint values selected in
● I. ROSEN et al.
577
TPEx were entered into Plan3D and used to create corresponding optimized plans. Figure 6 shows the results for each of the four beam configurations. The TPEx panel is shown with the selected critical structure dose limits. The dose distribution and DVH from the resulting plan are also shown. On the plan DVH, the critical structure dose–volume limits and resulting PTV values from TPEx are indicated for comparison. It is clear that the plan DVHs did not exactly match the TPEx values. The differences arise from the parameterization of the optimized plan characteristics and, as previously mentioned, the stochastic nature of the simulated annealing algorithm. Based on the DVHs in TPEx, the plans with more beams and with wedges were generally better, as judged by the amount of sparing of the left lung. However, the differences in plans became progressively less as the complexity increased. Therefore, the physician could readily gauge the clinical benefit versus delivery complexity. These results could be compared with the clinical plan used to treat the patient. However, these calculations used a dose model that is less accurate than Pinnacle3, especially with heterogeneities, and a different methodology for computing DVHs. Furthermore, there was no attempt to vary beam margins or to use partial beam blocking as was done in the clinical plan. DISCUSSION The treatment planning example presented is intended to illustrate the concept of interactive exploration. It is straightforward to extend this treatment-planning paradigm to include intensity-modulated radiation therapy (IMRT), noncoplanar beams, and electrons. Although the example presented used coplanar non-IMRT beams, this limitation is totally independent of the TPEx process or software. Because the data set (the Pareto front) that underlies the TPEx interface requires computing hundreds (preferably thousands) of optimized treatment plans, the size of the space that can be explored depends on the available computing power. Including IMRT delivery will require more computing capability than three-dimensional conformal radiation therapy because it adds another degree of freedom to the treatment planning problem, namely the fluence pattern for each beam. The optimization of each plan can still be based on maximizing the tumor dose subject to dose–volume constraints for the critical structures. The use of only four critical structures and one decision dose for each is also a computer limitation and not inherent in the interactive exploration concept. On a positive note, it is easy and straightforward to harness the power of parallel processing to this problem. After the set of optimization problems has been defined, each problem can be run independently in parallel on a separate processor. The example presented does not explore the space of all possible solutions for a given delivery method and set of decision doses. It is limited to the space defined by the parameters in Table 1. To explore the entire space of optimized solutions would require computing optimized plans for the complete range of volumes (0 –100%) for each
578
I. J. Radiation Oncology
● Biology ● Physics
Volume 61, Number 2, 2005
Fig. 5. The results of the linear fits (Eqs. 1–3) to the control point values for each of the delivery techniques— 4 open beams (a), 12 open beams (b), 4 wedged beams (c), and 12 wedged beams (d). The fit value is plotted against the corresponding plan value for each optimized plan and the root-mean-square error is shown. The trend lines are least-squares fits to the data.
critical structure and each decision dose. We expect that computers and calculation algorithms will be fast enough in the not-too-distant future to do that. Then, the finite speed of computation will restrict how many delivery options and decision doses can be explored, but the interactive exploration as shown in Fig. 3 will extend over the entire space of optimized solutions. We hypothesize that the PTV DVH of the true optimum solutions vary smoothly as each of the independent variables is changed independently. However, solutions produced by simulated annealing optimization have an associated uncertainty because of the stochastic nature of the method. Unlike directed search methods, simulated
annealing optimization has the potential to find a global minimum in a space with multiple local minima and it can easily incorporate “hard” constraints. However, its stochastic nature means that it never finds the exact solution to the problem nor does it produce the same solution in repetitions of the same calculation. It asymptotically approaches the true solution in a probabilistic way that depends on the length of the computation. The analytic fits, therefore, replace the measured “noisy” data with smoothly varying functions that should be better approximations to the true optimum solutions. These analytic functions also make the response of the user interface smoother and perhaps faster. A directed search
Exploring optimized plans
● I. ROSEN et al.
579
Fig. 5. (Continued).
optimization method, such as a gradient search, would likely give better parameter fits and might entirely eliminate the need for the analytic function. The multidimensional linear function was chosen because it is the simplest analytic function for the purpose. Figure 5 shows that it is a good approximation. Some of the variance in fit is due to the stochastic nature of the simulated annealing optimization algorithm. Of course, higher order functions could be applied and might improve the results. The dose distributions shown in Fig. 6 demonstrate the well-known hazard of making decisions based on DVHs alone. With open beams, there are rarely hot spots outside of the target volume. However, such hot spots often occur with wedged beams and IMRT beams if appropriate inverse
planning goals/constraints are not employed. The same problem exists with the TPEx process. In generating the space of optimized plans, suitable constraints must be applied to normal structures to avoid such hot spots. In addition, TPEx must display the DVHs for all normal structures of concern. If important features of the treatment plan, such as hot spots, are not controlled during the optimization process or displayed during the interactive exploration, then the physician is attempting to make decisions based on flawed or incomplete data and the process is meaningless. There is a minimum level of computer power needed to produce a data set that is sufficiently complete so that the DVHs presented during interactive exploration properly represent the underlying three-dimensional dose distributions. The parameters for that level of computing capability
580
I. J. Radiation Oncology
● Biology ● Physics
Volume 61, Number 2, 2005
Fig. 6. The critical structure dose–volume histogram (DVH) constraint values selected in Treatment Plan Explorer (TPEx) for each of the delivery techniques and the resulting optimized treatment plan— 4 open beams (a), 12 open beams (b), 4 wedged beams (c), and 12 wedged beams (d). On the plan DVH, the critical structure dose–volume limits and resulting PTV values from TPEx are indicated for comparison.
and for the completeness of the data set have not been determined. As computing power increases and more efficient methods for computing optimum treatment geometries are developed, more plans can be generated before interactive exploration. In addition to faster computer hardware, parallel processing can easily be used because each
optimization problem (different constraint values) can be computed independently. In the current process of preparing optimized plans for TPEx, critical structure constraint values are varied in uniform step sizes over a range of values. Other more efficient methods of searching the range of dose–volume constraints are certainly possible. Faster plan optimization means that a greater
Exploring optimized plans
● I. ROSEN et al.
581
Fig. 6. (Continued).
range of treatment delivery options can be computed, incorporating more critical structures and more decision doses per structure. Then, at interactive treatment planning, the physician or physicist could explore the treatment possibilities at many more levels of complexity and make a more informed decision about the greater effort required for more complex treatments versus the clinical benefit to the patient. The general goal of TPEx is the same as that of inverse
planning—to eliminate the details of treatment delivery from the planning process. However, rather than attempting to define a mathematical formulation of the clinical decision process (objectives and constraints), TPEx analyzes the relationships between normal tissue doses and target volume doses under conditions of optimum delivery and allows the user to interactively explore those relationships. The planner is effectively traveling along the boundary of all optimal solutions. We believe that
582
I. J. Radiation Oncology
● Biology ● Physics
this interactive process can make treatment planning more efficient and result in treatment plans that are more
Volume 61, Number 2, 2005
tailored to the needs of the patient and the clinical judgment of the physician.
REFERENCES 1. Bedford JL, Webb S. Elimination of importance factors for clinically accurate selection of beam orientations, beam weights and wedge angles in conformal radiation therapy. Med Phys 2003;30:1788 –1804. 2. Bednarz G, Michalski D, Houser C, et al. The use of mixedinteger programming for inverse treatment planning with predefined field segments. Phys Med Biol 2002;47:2235–2245. 3. Cotrutz C, Xing L. Segment-based dose optimization using a genetic algorithm. Phys Med Biol 2003;48:2987–2998. 4. Deasy JO. Multiple local minima in radiotherapy optimization problems with dose-volume constraints. Med Phys 1997;24: 1157–1161. 5. Djajaputra D, Wu Q, Wu Y, et al. Algorithm and performance of a clinical IMRT beam-angle optimization system. Phys Med Biol 2003;48:3191–3212. 6. Gopal R, Starkschall G. Plan space: Representation of treatment plans in multidimensional space. Int J Radiat Oncol Biol Phys 2002;53:1328 –1336. 7. Langer M, Lee EK, Deasy JO, et al. Operations research applied to radiotherapy, an NCI-NSF-sponsored workshop February 7–9, 2002. Int J Radiat Oncol Biol Phys 2003;57: 762–768. 8. Langer M, Morrill S, Brown R, et al. A comparison of mixed integer programming and fast simulated annealing for optimizing beam weights in radiation therapy. Med Phys 1996; 23:957–964. 9. Langer M, Brown R, Morrill S, et al. A generic genetic algorithm for generating beam weights. Med Phys 1996; 23:965–971. 10. Michalski D, Xiao Y, Censor Y, et al. The dose-volume constraint satisfaction problem for inverse treatment planning with field segments. Phys Med Biol 2004;49:601– 616. 11. Morrill SM, Lane RG, Jacobson G, et al. Treatment plan optimization using constrained simulated annealing. Phys Med Biol 1991;36:1341–1361. 12. Pugachev A, Xing L. Incorporating prior knowledge into beam orientation optimization in IMRT. Int J Radiat Oncol Biol Phys 2002;54:1565–1574. 13. Rosen II, Lane RG, Morrill SM, et al. Treatment plan optimization using linear programming. Med Phys 1991;18:141– 152. 14. Rosen II, Lam KS, Lane RG, et al. Comparison of simulated annealing algorithms for conformal therapy treatment planning. Int J Radiat Oncol Biol Phys 1995;33:1091–1099. 15. Rowbottom CG, Khoo VS, Webb S. Simultaneous optimization of beam orientations and beam weights in conformal radiotherapy. Med Phys 2001;28:1696 –1702.
16. Spirou SV, Chui CS. A gradient inverse planning algorithm with dose-volume constraints. Med Phys 1998;25:321–333. 17. Starkschall G, Pollack A, Stevens CW. Treatment planning using a dose-volume feasibility search algorithm. Int J Radiat Oncol Biol Phys 2001;49:1419 –1427. 18. Stein J, Mohan R, Wang XH, et al. Number and orientations of beams in intensity-modulated radiation treatments. Med Phys 1997;24:149 –160. 19. Willoughby TR, Starkschall G, Janjan NA, et al. Evaluation and scoring of radiotherapy treatment plans using an artificial neural network. Int J Radiat Oncol Biol Phys 1996;34:923– 930. 20. Wu Q, Djajaputra D, Wu Y, et al. Intensity-modulated radiotherapy optimization with gEUD-guided dose-volume objectives. Phys Med Biol 2003;48:279 –291. 21. Yu Y, Schell MC, Zhang JB. Decision theoretic steering and genetic algorithm optimization: Application to stereotactic radiosurgery treatment planning. Med Phys 1997;24:1742– 1750. 22. Zhang X, Liu H, Wang X, et al. Speed and convergence properties of gradient algorithms for optimization of IMRT. Med Phys 2004;31:1141–1151. 23. Cotrutz C, Lahanas M, Kappas C, et al. A multiobjective gradient-based dose optimization algorithm for external beam conformal radiotherapy. Phys Med Biol 2001;46:2161–2175. 24. Liu H, Rosen I, Janjan N, et al. Treatment planning optimization based on response-surface modeling of cost function versus multiple constraints. CD ROM Proceedings of the World Congress on Medical Physics and Biomedical Engineering, 4 pp, July 23–28, 2000. 25. Meyer J, Phillips MH, Cho PS, et al. Application of influence diagrams to prostate intensity-modulated radiation therapy plan selection. Phys Med Biol 2004;49:1637–1653. 26. 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. 27. Xing L, Li JG, Donaldson S, et al. Optimization of importance factors in inverse planning. Phys Med Biol 1999;44:2525– 2536. 28. Alber M, Birkner M, Nusslin F. Tools for the analysis of dose optimization: II. Sensitivity analysis. Phys Med Biol 2002;47: N265–N270. 29. Rosen II, Loyd MD, Lane RG. Collimator scatter in modeling radiation beam profiles. Med Phys 1990;17:422– 428.
doi:10.1016/j.ijrobp.2004.09.022
PHYSICS CONTRIBUTION
INTERACTIVELY EXPLORING OPTIMIZED TREATMENT PLANS ISAAC ROSEN, PH.D.,* H. HELEN LIU, PH.D.,* NATHAN CHILDRESS, PH.D.,* AND ZHONGXING LIAO, M.D.† *Departments of Radiation Physics and †Radiation Oncology, The University of Texas M. D. Anderson Cancer Center, Houston, TX Purpose: A new paradigm for treatment planning is proposed that embodies the concept of interactively exploring the space of optimized plans. In this approach, treatment planning ignores the details of individual plans and instead presents the physician with clinical summaries of sets of solutions to well-defined clinical goals in which every solution has been optimized in advance by computer algorithms. Methods and Materials: Before interactive planning, sets of optimized plans are created for a variety of treatment delivery options and critical structure dose–volume constraints. Then, the dose–volume parameters of the optimized plans are fit to linear functions. These linear functions are used to show in real time how the target dose–volume histogram (DVH) changes as the DVHs of the critical structures are changed interactively. A bitmap of the space of optimized plans is used to restrict the feasible solutions. The physician selects the critical structure dose–volume constraints that give the desired dose to the planning target volume (PTV) and then those constraints are used to create the corresponding optimized plan. Results: The method is demonstrated using prototype software, Treatment Plan Explorer (TPEx), and a clinical example of a patient with a tumor in the right lung. For this example, the delivery options included 4 open beams, 12 open beams, 4 wedged beams, and 12 wedged beams. Beam directions and relative weights were optimized for a range of critical structure dose–volume constraints for the lungs and esophagus. Cord dose was restricted to 45 Gy. Using the interactive interface, the physician explored how the tumor dose changed as critical structure dose–volume constraints were tightened or relaxed and selected the best compromise for each delivery option. The corresponding treatment plans were calculated and compared with the linear parameterization presented to the physician in TPEx. The linear fits were best for the maximum PTV dose and worst for the minimum PTV dose. Based on the root-mean-square error between the fit values and their corresponding data values, the linear fit appears to be adequate, although higher order polynomials could give better results. Some of the variance in fit is due to the stochastic nature of the simulated annealing optimization algorithm, which does not reproduce the exact same results in repetitions of the same calculation. Using a directed search algorithm for plan optimization should produce better parameter fits and, therefore, better predictions of plan characteristics by TPEx. Conclusions: Using TPEx, the physician can easily select the optimum plan for a patient, with no imposed arbitrary definition of the “best” plan. More importantly, the physician can readily see what can be achieved for the patient with a given delivery technique. There is no more uncertainty about whether or not a better plan exists. By comparing the “best” plans for different delivery options (e.g., three-dimensional conformal radiotherapy versus intensity-modulated radiation therapy), the physician can gauge the clinical benefits of greater technical complexity. However, before the TPEx process can be clinical useful, faster computers and/or algorithms are needed and more studies are needed to better model the spaces of optimized solutions. © 2005 Elsevier Inc. Radiation therapy treatment planning, Inverse planning, Treatment plan optimization.
INTRODUCTION The goal of treatment planning, whether forward-based or inverse-based, is to find the best plan for the individual patient. A major difficulty in accomplishing that goal is the lack of universally accepted criteria for defining the “best” plan. Clearly, more dose to the tumor is better and less dose to the normal tissues is better. However, the ideal plan of
100% prescribed dose to the entire tumor volume and zero dose elsewhere is not attainable. In evaluating any given treatment plan, the radiation oncologist is faced with a variety of issues and choices. The primary question is whether a dose can be delivered to the tumor that will achieve a high probability of local control without serious treatment morbidity or sequelae. If that goal can be achieved, then the focus shifts to whether treatment
Reprint requests to: Isaac Rosen, Ph.D., Department of Radiation Physics, Unit 94, The University of Texas M. D. Anderson Cancer Center, 1515 Holcombe Blvd, Houston, TX 77030; Tel: (713) 5632635; Fax: (713) 563-2629; E-mail: [email protected]
Support in part by a sponsored research agreement from Philips Medical Systems. Received Jun 8, 2004, and in revised form Sep 15, 2004. Accepted for publication Sep 17, 2004. 570
Exploring optimized plans
morbidity can be further reduced (i.e., can normal tissue doses be lowered?). When there is a desire to escalate tumor doses because local control is not adequate, the question becomes one of finding the best compromise between potential complications and potential recurrence. The answers to these questions depend on the details of treatment delivery. In general, more complex treatments with more degrees of freedom in beam geometry and fluence produce better plans with greater separation between tumor dose and normal tissue doses. In designing the best plan for a patient, the physician must often compromise between the prescribed tumor dose, the uniformity of dose in the tumor volume, and the doses to critical structures. Further compromises may be required between the different adverse reactions of different critical structures. Unfortunately, the probabilities of complications for an individual patient are not binary functions and not currently predictable. The physician must try to anticipate how the patient will respond to the planned dosage and may need to compromise tumor dose to manage the toxicity. Because different physicians may choose different compromises based on their training and clinical experience and their understanding of the patient’s unique situation, the results of treatment planning vary with physician and with patient. If it were possible to mathematically define the characteristics of the best plan, then computer programs could unerringly find it for each patient. Attempts to accomplish this go back more than 35 years, and there is an extensive body of literature exploring the use of algorithms such as exhaustive search, linear programming, quadratic programming, mixed-integer programming, gradient search, simulated annealing, genetic algorithms, feasibility search, downhill simplex, and neural networks (e.g., 1–22). In theory, inverse planning (or treatment plan optimization, as it used to be called) could provide the unequivocal best treatment plan for a patient if the biologic responses of the individual patient’s tissues could be accurately modeled and the clinical compromises acceptable to the clinician were objectively and unambiguously known. Although mathematical models of tumor and tissue response have been developed, they are not yet accurate enough for prediction of individual patient responses and clinician decision criteria are not well simulated by computer algorithms. Consequently, inverse planning is still an iterative process. A major improvement of inverse planning over forward planning is the focus on clinical factors, such as normal tissue tolerances, rather than on the details of beam geometries. The existence of multiple competing objectives in the treatment planning problem has led to approaches that focus on the relationships between these objectives rather than attempting to define a mathematically “best” solution (e.g., 23–28). In general, these methods compute multiple plans with different optimization criteria to find the boundary of the space of optimized (Pareto) or feasible solutions. They may use sensitivity analysis to show how much tumor dose can be increased for a given increase in normal tissue dose
● I. ROSEN et al.
571
or may directly optimize the weighting factors assigned to the different objectives. The new treatment planning paradigm presented here follows this general approach and extends it by using the optimization criteria directly, rather than through surrogate preference factors, and by providing a real-time interactive interface into the space of optimized plans. It is worthwhile to ask why treatment planning is an iterative process. The clinician presumably knows what the ideal (or acceptable) dose distribution should look like. The physicist or dosimetrist generates a plan that he or she believes is the closest approximation to what the clinician’s ideal plan is. Why should there be any iterations? One reason may be that the planner has not correctly anticipated the physician’s goals and compromise preferences. For example, the planner might have chosen to spare equal portions of each kidney, whereas the physician’s preference might be to sacrifice one kidney and completely spare the other. Another reason for multiple iterations may be that the physician implicitly believes that a better plan exists and that through more effort (and perhaps experience) it can be found. In either case, the physician is faced with a process in which the characteristics of the best achievable plan are unknown. The dilemma facing the physician can be described with a simple example. Consider a patient with a lung tumor. The planner generates a conventional beam arrangement using large anteroposterior beams for an initial dose followed by a boost dose delivered with parallel-opposed oblique beams reduced in size to the primary target volume. The physician must balance the doses to the tumor, lungs, heart, esophagus, and spinal cord. Changing the angles of the oblique beams will change this balance, as will adding more beams or fluence modulation. Presumably, the planner has already optimized the directions of the oblique beams either through trial and error or previous experience. However, based on the single plan presented, the physician cannot know how much improvement can be achieved by adding more beams or adding fluence modulation. Therefore, the physician may ask for additional plans to assess the benefits of these changes versus the added complexity of treatment delivery. As fluence modulation or beams are added, the optimality of the beam directions is thrown into question. The physician no longer knows how good a plan can be achieved or what level of complexity is needed to achieve the clinical goals. The information that the physician needs to make the best possible treatment planning decision for a given patient can be encapsulated in the question, “What is the greatest tumor dose that can be delivered for specified normal tissue dose– volume constraints, regardless of the beam characteristics?” Given this information in an interactive format, the physician can readily select the plan that delivers the most tumor dose with acceptable normal tissue risk or the one that most reduces the dose to selected organs while maintaining a desired tumor dose level. The answer to this question requires that the best possible (optimized) dose–volume histogram (DVH) for the planning target volume (PTV) be
572
I. J. Radiation Oncology
● Biology ● Physics
presented to the physician for every combination of normal tissue DVHs that they might want to select. We are studying a new paradigm for treatment planning that embodies this concept of interactively exploring the space of optimized plans. As with some other inverse planning techniques, the details of treatment delivery (such as beam weights, angles, and doses) are hidden. Instead, the physician is presented with the relevant clinical summaries for the best possible solutions to well-defined clinical goals in which every solution has been optimized in advance by computer algorithms. Using an interactive interface, the physician selects a general treatment delivery method and then explores in real time how the maximum deliverable tumor dose changes as critical structure dose–volume constraints are tightened or relaxed. The physician can then directly select the compromises between doses to the various organs, dose to the tumor, and complexity of treatment delivery that he or she desires for the patient with no imposed arbitrary definition of the “best” plan. More importantly, the physician can see what is achievable for the patient. There is no more uncertainty about whether or not a better plan exists.
Volume 61, Number 2, 2005
Selection of plan parameters One goal of TPEx is to avoid dealing with the details of beam geometry. The first user screen (Fig. 2) asks for only three parameters to define the treatment delivery—the beam energy, the number of beams, and the degree of fluence modulation. Based on these selections, the program reads a set of previously computed optimized plans. It then presents the user with choices of critical structures to manipulate during exploration and the corresponding available decision doses. The decision doses are dose levels at which the physician would like to control the irradiated volume. Typically, these would be dose levels at which the irradiated volume is predictive of complications. The choices available to the user are limited to those for which optimized sets of plans were previously generated. Therefore, the planner should consult with the physician before beginning the optimization phase or have sufficient experience to anticipate the physician’s general approach to treatment for the given anatomic site. The planner should know what decision doses the physician uses in evaluating plans, what beam energies are preferred, and how much treatment complexity will probably be needed to accomplish the clinical goals. In that way, the planner can choose the best set of plans to be optimized in preparation for treatment plan exploration.
Interactive exploration METHODS AND MATERIALS A prototype system named Treatment Plan Explorer (TPEx) was developed using MATLAB (The MathWorks Inc., Natick, MA) to test this approach to treatment planning. It imports a set of optimized treatment plans, parameterizes their dose–volume characteristics, and then presents the results to the physician or planner in an interactive format. Before using TPEx, multiple sets of optimized plans are created, each of which covers the range of potentially desirable normal tissue DVH constraints. Each set corresponds to a different treatment delivery technique that the user may wish to explore. TPEx does not include external beam radiation treatment planning (RTP) software. The optimized plans used by TPEx are created by the RTP software. The output of TPEx is a set of critical structure DVH constraints selected by the planner, which are then used by the RTP software to create the final optimized plan for the patient. Treatment Plan Explorer has two user interface screens. In the first, the planner chooses the basic plan parameters and the critical structures to be considered. In the second, the planner interactively adjusts the dose–volume constraints on the critical structures and observes the effects on the DVH of the PTV. During exploration, TPEx continuously displays the values of the user parameters and the resulting PTV parameters. The use of TPEx is demonstrated using a clinical example. Figure 1 shows three sections through a patient with a tumor in the right lung treated using three coplanar 6-MV beams. The anterior beam included a 30° wedge and the posterior beam included a 45° wedge. The right posterior-oblique beam included a 15° wedge and partial blocking. The dose distribution delivered 63 Gy to the 95% isodose line. The DVH of the treatment is also shown. The maximum dose to the spinal cord was 45.5 Gy; the volume of right lung exceeding 20 Gy was 59%; the volume of left lung exceeding 20 Gy was 19%; and the volume of esophagus exceeding 50 Gy was 21%. Because most of the PTV was superior to the heart, only a small volume of heart (10%) received more than 45 Gy from the coplanar beam arrangement.
The interactive exploration of optimized solutions is shown in Fig. 3. The upper graph is a simplified DVH for the PTV showing minimum dose, 95% volume dose, and maximum dose. The bottom graphs are simple DVHs for the critical structures selected in the previous input screen. Each critical structure DVH is defined by one or two control points that are individually movable by the user. The maximum dose control points can only move laterally, changing dose but not volume. The decision dose control points can only move vertically, changing volume but not decision dose value. For each point, a red bar shows the range of permissible motion. As the planner changes the value of one of the control points (the “active” control point), two things happen in real time. First, the DVH of the PTV changes to reflect the best plan that corresponds to the current settings of all the control points. Second, the red bars showing the ranges of motion of the control points change dynamically, including the range for the active control point. The active control point is not allowed to move outside of its range. Because of the complexity of the space of optimized plans, there are strong interactions between the control points and there may be discontinuities in the space. It is possible that the range of motion of a nonactive control point may change in such a way that the corresponding control point is no longer within its range. That situation represents a nonfeasible solution and should not be selected. When that control point is activated, it is automatically relocated to a value within its range. If a situation arises in which the maximum dose control point is reduced to the decision dose value, then the volume corresponding to the decision dose is automatically reduced to 0% because the decision dose is now also the maximum dose. The planner should explore all of the control points in various ways and sequences to find the set of DVHs (PTV and critical structures) that is most desirable. Just as in conventional treatment planning, lower doses to the critical structures result in lower doses to the PTV. However, by observing the changes in the PTV DVH with changes in the control points, it is easy to see how much each critical structure limits the PTV dose and how much sparing of
Exploring optimized plans
● I. ROSEN et al.
573
Fig. 1. Transverse, sagittal, and coronal sections through a patient with a tumor in the right lung. The patient was treated with three coplanar 6-MV beams. The anterior beam included a 30° wedge and the posterior beam included a 45° degree wedge. The right posterior-oblique beam included a 15° wedge and partial blocking. The dose distribution delivered 63 Gy to the 95% isodose line.
each critical structure can be achieved for a given PTV dose. In addition, the simultaneous changes in the magnitude of the prescription dose and the degree of dose homogeneity in the PTV can be observed in response to alterations in critical structure constraint values. After a set of control point values is chosen, those critical structure constraints are used as input to the RTP software to create an optimized plan with the beam characteristics (energy, number, modulation) chosen by the user at the start of the process. Because the user is interrogating a parameterization of the original data set of optimized plans and also interpolating, the actual plan achieved will not exactly match the results of TPEx. How closely the results of TPEx will match the results of RTP software optimization with the desired control point constraints depends on several factors— the density of data in the space of optimized plans, the optimization algorithm, and the quality of the function fitting. These relationships have not yet been explored.
Creation of optimized plan sets The patient treatment was originally planned using Pinnacle3 (Philips Radiation Oncology Systems, Madison, WI). Normal tissue and tumor volume structures were extracted from Pinnacle3 and imported into a research planning system, Plan3D, which uses a modified Clarkson scatter-air ratio dose calculation (29) and simulated annealing for plan optimization (14). All cal-
Fig. 2. The first Treatment Plan Explorer screen asks for only three parameters to define the treatment delivery—the beam energy, the number of beams, and the degree of fluence modulation. Based on these selections, the program reads a set of previously optimized plans. It then presents the user with choices of critical structures to manipulate during exploration and the corresponding available decision doses.
574
I. J. Radiation Oncology
● Biology ● Physics
Volume 61, Number 2, 2005
Fig. 3. The interactive exploration of all the possible optimized solutions is shown. The upper graph is a simplified dose–volume histogram (DVH) for the planning target volume showing minimum dose, 95% volume dose, and maximum dose. The bottom graphs are simple DVHs for the critical structures selected in the previous input screen. Each critical structure DVH is defined by one or two control points that are individually movable by the user. The maximum dose control points can only move laterally, changing dose but not volume. The decision dose control points can only move vertically, changing volume but not decision dose value. For each point, a red bar shows the range of permissible motion.
culations included the full three-dimensional anatomy and bulk tissue heterogeneity effects (Batho power-law correction). DVHs were computed using quasi-random points rather than a Cartesian grid. The number of points used per structure was 1000 for the PTV, 600 for each lung, 100 for the esophagus, and 100 for the spinal cord. Eight optimized data sets using 6-MV photons were created, as shown in Table 1. For each type of modulation, open and wedge, plan sets were created with 4, 6, 8, and 12 coplanar beams. In each case, the relative beam weights and their directions were optimized. Beam direction optimization was accomplished by starting with 36 eligible beams equally spaced about the patient at 10° angular increments. Beam weights were initially optimized for all 36 beams. Then the number of eligible beams was reduced by eliminating those with low contributions to the plan. This process was repeated up to four times until the desired number of beams remained. For wedge plans, the wedge angles and directions were also optimized. Beam shapes were automatically generated to conform to the PTV projection with 0.5-cm margins, similar to the margins used in the clinical plan. However, unlike the clinical plan, no partial beam blocking was included. For each delivery method (modulation and beam number), a set of optimized plans was automatically generated using the various com-
binations of dose–volume constraints for the critical structures shown in Table 1. The optimization objective was to maximize the minimum PTV dose subject to strictly enforced dose–volume constraints (“hard” constraints). For each decision dose, an optimized plan was created using each volume limit. For this example case, there were seven different volume levels permitted to exceed 20 Gy for the right lung, seven different volume levels permitted to exceed 20 Gy for the left lung, and six different volume levels permitted to exceed 50 Gy for the esophagus. The spinal cord was limited to 45 Gy maximum dose, which was not varied. An optimized plan was generated for each combination of dose–volume constraints. Therefore, there were 294 plans in each set (7 ⫻ 7 ⫻ 6). No constraints were applied to the PTV so that there was always a feasible solution. The actual volumes exceeding the decision dose values were extracted from the DVHs for all the structures in each plan and recorded. Because of the optimization method, these values were always less than or equal to the constraint values and represented what could theoretically be achieved by the plan. The maximum dose to each critical structure was also recorded. These critical structure data were the control points for TPEx. For the PTV, the minimum dose, maximum dose, and prescription dose (95% volume) were recorded for each plan. An example of the output from one optimization is shown in Table 2. The constraint values used
Exploring optimized plans
● I. ROSEN et al.
575
Table 1. Parameter values for the sets of optimized plans Modulation
Number of beams
Critical structure
Maximum dose (Gy)
Decision dose (Gy)
Volume limit for decision dose (%)
4 6 8 12
Right lung Left lung Esophagus Spinal cord
75 75 75 45
20 20 50 45
20, 30, 40, 50, 60, 70, 80 10, 20, 30, 40, 50, 60, 70 10, 20, 30, 40, 50, 60 0
Open Wedge
For each decision dose, an optimized plan was created using each volume limit shown.
the number of independent variables is twice the number of critical structures.
as input to the optimization process and the resulting beam parameters (directions and doses) that produced the optimum plan were not used in TPEx.
Boundaries of the space of optimized plans
Parameterization of optimized plans
It is important that the user’s exploration be limited to the space of optimized solutions. Because of the parameterization, it is possible to compute PTV values for any set of control point values. However, not all combinations of control point values correspond to feasible solutions. Therefore, the motions of control points must be restricted by the boundaries of the space of optimized plans. This boundary will generally be complex and nonanalytic. The boundaries for any given independent variable (control point) will vary with the specific values of the other independent variables. The space of optimized plans was quantized using a bitmap of dimension Rm, where R is the desired resolution of the bitmap and m is the number of independent variables. Each dimension of the bitmap corresponds to one of the independent variables. The scale associated with dimension i was computed by
The PTV characteristics of the optimized plans were fit to linear functions of the control points, PTVmin ⫽ a0min ⫹
n
兺a
VDD i ⫹
min i
i⫽1
PTV95% ⫽ a095% ⫹
VDD i ⫹
95% i
i⫽1
PTVmax ⫽ a0max ⫹
n
兺a
min i
Dmax i
(1)
i⫽1
n
兺a
n
兺b n
兺b
95% i
Dmax i
(2)
max i
Dmax i
(3)
i⫽1
VDD i ⫹
max i
i⫽1
n
兺b i⫽1
where
⌬xi ⫽ (xmax ⫺ xmin) ⁄ (R ⫺ 1)
PTVmin ⫽ minimum dose in PTV PTV95% ⫽ minimum dose in 95% of PTV PTVmax ⫽ maximum dose in PTV VDD ⫽ volume exceeding the decision dose value for critical i organ i Dmax ⫽ maximum dose in critical organ i i n ⫽ number of critical organs ai, bi ⫽ fit parameters The control points (VDD and Dmax i i ) are the independent variables and the PTV doses are the dependent variables. For this example,
(4)
where xmax and xmin are the maximum and minimum values, respectively, of the independent variable xi, ⌬xi is the increment in values along the xi axis, and R is the dimension of the bitmap along each axis. Applying Eq. 4 to each independent variable produces a set of control point values for each point in the bitmap of the space of optimized plans. Note that the bitmap can be very large. For this example with eight independent variables, a resolution of eight values per variable was used, giving a bitmap with 16,777,216 elements.
Table 2. Sample output from the optimization process Optimization input
Critical structure Right lung Left lung Esophagus Spinal cord PTV minimum dose (100% volume) PTV 95% volume dose PTV maximum dose (0% volume)
Decision dose (Gy)
Volume limit for decision dose (%)
20 20 50 45
20 10 5 0
Optimization output Plan volume exceeding decision dose (%) 20.0 3.0 0.0 0.0
Maximum dose (Gy) 27.6 33.3 30.7 30.6
18.0 Gy 19.0 Gy 33.7 Gy
Abbreviations: PTV ⫽ planning target volume; TPEx ⫽ treatment plan explorer. The decision dose values, plan volumes exceeding the decision dose, maximum dose, and PTV doses were used for interactive exploration (highlighted). The constraint volume limits for the decision doses were used for optimization but not in TPEx.
576
I. J. Radiation Oncology
● Biology ● Physics
Volume 61, Number 2, 2005
Fig. 4. A sequence of steps that leads to the best plan selected for 12 wedged beams. For each plan in the optimized set, the corresponding point in the bitmap is set to a value of 1. Consequently, the boundary of the space of optimized plans is defined by the edges of the volume of 1s in the bitmap. For any set of values of m-1 variables, the range of acceptable values for the mth variable can easily be found by scanning the mth axis of the bitmap. The values corresponding to the limits of the bits with value of 1 are the range for the variable. Depending on the density of data in the space of optimized plans, there may be gaps in the bitmap from undersampling. Therefore, the extreme limits of the 1s
along an axis are taken as the range of the variable, even though the 1s may not be continuous.
Initial solution The interactive exploration can start with any plan from the optimized set. Ideally, there should be no isolated regions in the space of optimized plans, so that it should be possible to reach any solution from any starting point. However, a good starting solution can be beneficial for two reasons. First, finite sampling
Exploring optimized plans
of the space might, in fact, produce isolated regions that could trap the user and prevent reaching the most desirable one. Second, there is no reason not to expect discontinuities in the space of optimized solutions. In suggesting that there are discontinuities, we mean that every point in the space is not linearly connected to every other point. Although it should be possible to reach every point in the space from every other point, the sequences of moves along the multiple axes of the multidimensional space may not be commutative. Finding a particular solution will depend on the sequence of steps as well as the starting point. Therefore, a good starting solution will reduce the amount of interactive searching needed to find the most desirable one. We chose as our initial solution the one with the largest PTV maximum dose.
RESULTS One of us (Z.L.) used TPEx to select the best plans for treatment configurations using 4 open beams, 12 open beams, 4 wedged beams, and 12 wedged beams. For each beam configuration, the control point values (decision dose volumes and maximum dose values) were interactively adjusted until a solution was found that gave the best compromise in doses to the critical structures and PTV. Figure 4 shows a sequence of steps that led to the best plan selected for 12 wedged beams. The initial solution is shown in Fig. 4a. First, the volume of left lung exceeding 20 Gy (V20) was reduced from 43.3% to 26.7% (the bottom of the available range) (Fig. 4b). Then, the maximum dose to the esophagus was reduced from 75 Gy to 68 Gy (Fig. 4c). This reduction changed the range of available values for the decision dose for the left lung, allowing the V20 to be further reduced to 9.4% (Fig. 4d). Next, reducing the maximum dose to the right lung from 75 Gy to 71.8 Gy improved the uniformity of dose in the PTV by reducing the maximum PTV dose (Fig. 4e). Finally, increasing the maximum esophagus dose from 68 Gy to 70.6 Gy further improved the PTV dose uniformity by increasing the minimum dose (Fig. 4f). The first two steps (Figs. 4b and 4c) show that reducing the maximum dose to the esophagus lowers the range of possible doses to V20 of the left lung. It seems reasonable that the lower dose levels for V20 of the left lung should have been available in the initial solution. We attribute this deficiency to the coarse sampling of the space of optimized plans and the coarse quantization of the space by the bitmap. Both of these are consequences of limited computing resources. The actual sequence of steps originally used to produce this final plan was not recorded. The results of the linear fits (Eqs. 1–3) to the control point values are shown in Fig. 5 for each of the delivery techniques. The fit value is plotted against the corresponding plan value for each optimized plan and the root-mean-square error is shown. The trend lines are least-squares fits to the data. The fits were best for the maximum PTV dose and worst for the minimum PTV dose. In general, the linear model seems to be appropriate and higher order polynomials have not yet been tested. The critical structure DVH constraint values selected in
● I. ROSEN et al.
577
TPEx were entered into Plan3D and used to create corresponding optimized plans. Figure 6 shows the results for each of the four beam configurations. The TPEx panel is shown with the selected critical structure dose limits. The dose distribution and DVH from the resulting plan are also shown. On the plan DVH, the critical structure dose–volume limits and resulting PTV values from TPEx are indicated for comparison. It is clear that the plan DVHs did not exactly match the TPEx values. The differences arise from the parameterization of the optimized plan characteristics and, as previously mentioned, the stochastic nature of the simulated annealing algorithm. Based on the DVHs in TPEx, the plans with more beams and with wedges were generally better, as judged by the amount of sparing of the left lung. However, the differences in plans became progressively less as the complexity increased. Therefore, the physician could readily gauge the clinical benefit versus delivery complexity. These results could be compared with the clinical plan used to treat the patient. However, these calculations used a dose model that is less accurate than Pinnacle3, especially with heterogeneities, and a different methodology for computing DVHs. Furthermore, there was no attempt to vary beam margins or to use partial beam blocking as was done in the clinical plan. DISCUSSION The treatment planning example presented is intended to illustrate the concept of interactive exploration. It is straightforward to extend this treatment-planning paradigm to include intensity-modulated radiation therapy (IMRT), noncoplanar beams, and electrons. Although the example presented used coplanar non-IMRT beams, this limitation is totally independent of the TPEx process or software. Because the data set (the Pareto front) that underlies the TPEx interface requires computing hundreds (preferably thousands) of optimized treatment plans, the size of the space that can be explored depends on the available computing power. Including IMRT delivery will require more computing capability than three-dimensional conformal radiation therapy because it adds another degree of freedom to the treatment planning problem, namely the fluence pattern for each beam. The optimization of each plan can still be based on maximizing the tumor dose subject to dose–volume constraints for the critical structures. The use of only four critical structures and one decision dose for each is also a computer limitation and not inherent in the interactive exploration concept. On a positive note, it is easy and straightforward to harness the power of parallel processing to this problem. After the set of optimization problems has been defined, each problem can be run independently in parallel on a separate processor. The example presented does not explore the space of all possible solutions for a given delivery method and set of decision doses. It is limited to the space defined by the parameters in Table 1. To explore the entire space of optimized solutions would require computing optimized plans for the complete range of volumes (0 –100%) for each
578
I. J. Radiation Oncology
● Biology ● Physics
Volume 61, Number 2, 2005
Fig. 5. The results of the linear fits (Eqs. 1–3) to the control point values for each of the delivery techniques— 4 open beams (a), 12 open beams (b), 4 wedged beams (c), and 12 wedged beams (d). The fit value is plotted against the corresponding plan value for each optimized plan and the root-mean-square error is shown. The trend lines are least-squares fits to the data.
critical structure and each decision dose. We expect that computers and calculation algorithms will be fast enough in the not-too-distant future to do that. Then, the finite speed of computation will restrict how many delivery options and decision doses can be explored, but the interactive exploration as shown in Fig. 3 will extend over the entire space of optimized solutions. We hypothesize that the PTV DVH of the true optimum solutions vary smoothly as each of the independent variables is changed independently. However, solutions produced by simulated annealing optimization have an associated uncertainty because of the stochastic nature of the method. Unlike directed search methods, simulated
annealing optimization has the potential to find a global minimum in a space with multiple local minima and it can easily incorporate “hard” constraints. However, its stochastic nature means that it never finds the exact solution to the problem nor does it produce the same solution in repetitions of the same calculation. It asymptotically approaches the true solution in a probabilistic way that depends on the length of the computation. The analytic fits, therefore, replace the measured “noisy” data with smoothly varying functions that should be better approximations to the true optimum solutions. These analytic functions also make the response of the user interface smoother and perhaps faster. A directed search
Exploring optimized plans
● I. ROSEN et al.
579
Fig. 5. (Continued).
optimization method, such as a gradient search, would likely give better parameter fits and might entirely eliminate the need for the analytic function. The multidimensional linear function was chosen because it is the simplest analytic function for the purpose. Figure 5 shows that it is a good approximation. Some of the variance in fit is due to the stochastic nature of the simulated annealing optimization algorithm. Of course, higher order functions could be applied and might improve the results. The dose distributions shown in Fig. 6 demonstrate the well-known hazard of making decisions based on DVHs alone. With open beams, there are rarely hot spots outside of the target volume. However, such hot spots often occur with wedged beams and IMRT beams if appropriate inverse
planning goals/constraints are not employed. The same problem exists with the TPEx process. In generating the space of optimized plans, suitable constraints must be applied to normal structures to avoid such hot spots. In addition, TPEx must display the DVHs for all normal structures of concern. If important features of the treatment plan, such as hot spots, are not controlled during the optimization process or displayed during the interactive exploration, then the physician is attempting to make decisions based on flawed or incomplete data and the process is meaningless. There is a minimum level of computer power needed to produce a data set that is sufficiently complete so that the DVHs presented during interactive exploration properly represent the underlying three-dimensional dose distributions. The parameters for that level of computing capability
580
I. J. Radiation Oncology
● Biology ● Physics
Volume 61, Number 2, 2005
Fig. 6. The critical structure dose–volume histogram (DVH) constraint values selected in Treatment Plan Explorer (TPEx) for each of the delivery techniques and the resulting optimized treatment plan— 4 open beams (a), 12 open beams (b), 4 wedged beams (c), and 12 wedged beams (d). On the plan DVH, the critical structure dose–volume limits and resulting PTV values from TPEx are indicated for comparison.
and for the completeness of the data set have not been determined. As computing power increases and more efficient methods for computing optimum treatment geometries are developed, more plans can be generated before interactive exploration. In addition to faster computer hardware, parallel processing can easily be used because each
optimization problem (different constraint values) can be computed independently. In the current process of preparing optimized plans for TPEx, critical structure constraint values are varied in uniform step sizes over a range of values. Other more efficient methods of searching the range of dose–volume constraints are certainly possible. Faster plan optimization means that a greater
Exploring optimized plans
● I. ROSEN et al.
581
Fig. 6. (Continued).
range of treatment delivery options can be computed, incorporating more critical structures and more decision doses per structure. Then, at interactive treatment planning, the physician or physicist could explore the treatment possibilities at many more levels of complexity and make a more informed decision about the greater effort required for more complex treatments versus the clinical benefit to the patient. The general goal of TPEx is the same as that of inverse
planning—to eliminate the details of treatment delivery from the planning process. However, rather than attempting to define a mathematical formulation of the clinical decision process (objectives and constraints), TPEx analyzes the relationships between normal tissue doses and target volume doses under conditions of optimum delivery and allows the user to interactively explore those relationships. The planner is effectively traveling along the boundary of all optimal solutions. We believe that
582
I. J. Radiation Oncology
● Biology ● Physics
this interactive process can make treatment planning more efficient and result in treatment plans that are more
Volume 61, Number 2, 2005
tailored to the needs of the patient and the clinical judgment of the physician.
REFERENCES 1. Bedford JL, Webb S. Elimination of importance factors for clinically accurate selection of beam orientations, beam weights and wedge angles in conformal radiation therapy. Med Phys 2003;30:1788 –1804. 2. Bednarz G, Michalski D, Houser C, et al. The use of mixedinteger programming for inverse treatment planning with predefined field segments. Phys Med Biol 2002;47:2235–2245. 3. Cotrutz C, Xing L. Segment-based dose optimization using a genetic algorithm. Phys Med Biol 2003;48:2987–2998. 4. Deasy JO. Multiple local minima in radiotherapy optimization problems with dose-volume constraints. Med Phys 1997;24: 1157–1161. 5. Djajaputra D, Wu Q, Wu Y, et al. Algorithm and performance of a clinical IMRT beam-angle optimization system. Phys Med Biol 2003;48:3191–3212. 6. Gopal R, Starkschall G. Plan space: Representation of treatment plans in multidimensional space. Int J Radiat Oncol Biol Phys 2002;53:1328 –1336. 7. Langer M, Lee EK, Deasy JO, et al. Operations research applied to radiotherapy, an NCI-NSF-sponsored workshop February 7–9, 2002. Int J Radiat Oncol Biol Phys 2003;57: 762–768. 8. Langer M, Morrill S, Brown R, et al. A comparison of mixed integer programming and fast simulated annealing for optimizing beam weights in radiation therapy. Med Phys 1996; 23:957–964. 9. Langer M, Brown R, Morrill S, et al. A generic genetic algorithm for generating beam weights. Med Phys 1996; 23:965–971. 10. Michalski D, Xiao Y, Censor Y, et al. The dose-volume constraint satisfaction problem for inverse treatment planning with field segments. Phys Med Biol 2004;49:601– 616. 11. Morrill SM, Lane RG, Jacobson G, et al. Treatment plan optimization using constrained simulated annealing. Phys Med Biol 1991;36:1341–1361. 12. Pugachev A, Xing L. Incorporating prior knowledge into beam orientation optimization in IMRT. Int J Radiat Oncol Biol Phys 2002;54:1565–1574. 13. Rosen II, Lane RG, Morrill SM, et al. Treatment plan optimization using linear programming. Med Phys 1991;18:141– 152. 14. Rosen II, Lam KS, Lane RG, et al. Comparison of simulated annealing algorithms for conformal therapy treatment planning. Int J Radiat Oncol Biol Phys 1995;33:1091–1099. 15. Rowbottom CG, Khoo VS, Webb S. Simultaneous optimization of beam orientations and beam weights in conformal radiotherapy. Med Phys 2001;28:1696 –1702.
16. Spirou SV, Chui CS. A gradient inverse planning algorithm with dose-volume constraints. Med Phys 1998;25:321–333. 17. Starkschall G, Pollack A, Stevens CW. Treatment planning using a dose-volume feasibility search algorithm. Int J Radiat Oncol Biol Phys 2001;49:1419 –1427. 18. Stein J, Mohan R, Wang XH, et al. Number and orientations of beams in intensity-modulated radiation treatments. Med Phys 1997;24:149 –160. 19. Willoughby TR, Starkschall G, Janjan NA, et al. Evaluation and scoring of radiotherapy treatment plans using an artificial neural network. Int J Radiat Oncol Biol Phys 1996;34:923– 930. 20. Wu Q, Djajaputra D, Wu Y, et al. Intensity-modulated radiotherapy optimization with gEUD-guided dose-volume objectives. Phys Med Biol 2003;48:279 –291. 21. Yu Y, Schell MC, Zhang JB. Decision theoretic steering and genetic algorithm optimization: Application to stereotactic radiosurgery treatment planning. Med Phys 1997;24:1742– 1750. 22. Zhang X, Liu H, Wang X, et al. Speed and convergence properties of gradient algorithms for optimization of IMRT. Med Phys 2004;31:1141–1151. 23. Cotrutz C, Lahanas M, Kappas C, et al. A multiobjective gradient-based dose optimization algorithm for external beam conformal radiotherapy. Phys Med Biol 2001;46:2161–2175. 24. Liu H, Rosen I, Janjan N, et al. Treatment planning optimization based on response-surface modeling of cost function versus multiple constraints. CD ROM Proceedings of the World Congress on Medical Physics and Biomedical Engineering, 4 pp, July 23–28, 2000. 25. Meyer J, Phillips MH, Cho PS, et al. Application of influence diagrams to prostate intensity-modulated radiation therapy plan selection. Phys Med Biol 2004;49:1637–1653. 26. 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. 27. Xing L, Li JG, Donaldson S, et al. Optimization of importance factors in inverse planning. Phys Med Biol 1999;44:2525– 2536. 28. Alber M, Birkner M, Nusslin F. Tools for the analysis of dose optimization: II. Sensitivity analysis. Phys Med Biol 2002;47: N265–N270. 29. Rosen II, Loyd MD, Lane RG. Collimator scatter in modeling radiation beam profiles. Med Phys 1990;17:422– 428.