PAGE 1
NO VEL METHODS F OR INTENSITY MODULA TED RADIA TION THERAPY TREA TMENT PLANNING By AR VIND KUMAR A DISSER T A TION PRESENTED TO THE GRADUA TE SCHOOL OF THE UNIVERSITY OF FLORID A IN P AR TIAL FULFILLMENT OF THE REQUIREMENTS F OR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORID A 2005
PAGE 2
I dedicate this w ork to Sri Rani Sati Dadiji
PAGE 3
A CKNO WLEDGMENTS First of all, I w ould lik e to thank Dr. Ra vindra K. Ah uja, who has b een a great sup ervisor and men tor throughout m y four y ears at the Univ ersit y of Florida. His inspiring en th usiasm and energy ha v e b een con tagious, and his advice and constan t supp ort ha v e b een extremely helpful. I w ould also lik e to express m y gratitude to Dr. H. Edwin Romeijn and Dr. James F. Dempsey for their constan t guidance throughout m y do ctoral studies. I feel really fortunate to ha v e the opp ortunit y to w ork with p eople of their calib er. I w ould lik e to thank Dr. Siriphong La wphongpanic h, Dr. Stanisla v Ury asev, and Dr. Sanja y Rank a for serving as m y sup ervisory committee mem b ers and giving me though tful advice. My sp ecial ac kno wledgmen t go es to m y collab orators and colleagues at Univ ersit y of Florida and at other places, esp ecially , Dr. Jonathan G. Li, Dr. Jian Liu, Dr. Claudio Cunha, Dr Krishna C. Jha, G uv en c S ahin, Mic helle Hanna, On ur Seref, Bala V aidy anathan, Kamalesh Somani and Suc handan Guha. I am indebted to Srik an t Ranjan for pro ofreading m y dissertation and pro viding v aluable suggestions. My utmost appreciation go es to m y paren ts, the source of all m y strength, and to m y family for their encouragemen t and eternal supp ort. iii
PAGE 4
T ABLE OF CONTENTS page A CKNO WLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF T ABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii ABSTRA CT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x CHAPTER1 INTR ODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Dissertation Summary . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Fluence Map Optimization . . . . . . . . . . . . . . . . . . 3 1.1.2 Beam Orien tation Optimization . . . . . . . . . . . . . . . 4 1.1.3 T reatmen t Deliv ery . . . . . . . . . . . . . . . . . . . . . . 5 1.1.4 Ap erture Mo dulation . . . . . . . . . . . . . . . . . . . . . 7 1.2 Con tribution Summary . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.1 Fluence Map Optimization . . . . . . . . . . . . . . . . . . 8 1.2.2 Beam Orien tation Optimization . . . . . . . . . . . . . . . 8 1.2.3 T reatmen t Deliv ery . . . . . . . . . . . . . . . . . . . . . . 9 1.2.4 Ap erture Mo dulation . . . . . . . . . . . . . . . . . . . . . 9 2 FLUENCE MAP OPTIMIZA TION . . . . . . . . . . . . . . . . . . . . . 10 2.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Notation and Dose Calculation . . . . . . . . . . . . . . . . . . . . 13 2.3 Mo del F orm ulation . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 A Basic Linear Programming F orm ulation . . . . . . . . . . 15 2.3.2 The Dose-V olume Histogram . . . . . . . . . . . . . . . . . 21 2.3.3 Constrain ts on the Dose-V olume Histogram . . . . . . . . . 22 2.3.4 An Enhanced LP F orm ulation with Dose-V olume Histogram Constrain ts . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.5 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.5.1 Restoring F ractionation . . . . . . . . . . . . . . . . . . . . 35 2.5.2 Incorp orating Spatial Eects: The 2-Dimensional D VH . . 40 2.6 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . 44 iv
PAGE 5
3 BEAM ANGLE OPTIMIZA TION . . . . . . . . . . . . . . . . . . . . . . 46 3.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2 Metho dology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.2.1 A Mixed In teger Linear Programming F orm ulation for the BOO Problem . . . . . . . . . . . . . . . . . . . . . . . . 49 3.2.2 Searc hing the Solution Space for Equispaced Beams . . . . 51 3.3 Neigh b orho o d Searc h Algorithm . . . . . . . . . . . . . . . . . . . 52 3.3.1 Finding an Initial F easible Solution . . . . . . . . . . . . . 53 3.3.2 The Neigh b orho o d Searc h Algorithm . . . . . . . . . . . . 56 3.4 Computational Results and Analysis . . . . . . . . . . . . . . . . 60 3.4.1 Impro v emen t due to Equispaced BOO Metho dology . . . . 60 3.4.2 Results of the Neigh b orho o d Searc h Metho d . . . . . . . . 61 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4 DELIVER Y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.1.1 Problem Description . . . . . . . . . . . . . . . . . . . . . . 66 4.1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . 70 4.2 Net w ork-Based Mo del . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.2.1 Shap e Matrix Graph . . . . . . . . . . . . . . . . . . . . . 74 4.2.2 Inclusion of Connectedness Constrain ts . . . . . . . . . . . 77 4.2.3 Inclusion of T ongue and Gro o v e Constrain ts . . . . . . . . 78 4.3 Heuristics for Ecien t Deliv ery . . . . . . . . . . . . . . . . . . . 80 4.3.1 Extracting Shap e Matrices . . . . . . . . . . . . . . . . . . 80 4.3.2 Enforcing V arious Constrain ts in the Net w ork . . . . . . . 83 4.3.3 Sc hemes for Extraction . . . . . . . . . . . . . . . . . . . . 86 4.4 Maxim um Leaf Spread Constrain t . . . . . . . . . . . . . . . . . . 91 4.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . 94 4.5.1 Comparison of Metho ds for Unconstrained Case . . . . . . 94 4.5.2 Comparison of Metho ds With In terdigitation and T ongue and Gro o v e Constrain ts . . . . . . . . . . . . . . . . . . 95 4.5.3 Results for Maxim um Leaf Separation Constrain t . . . . . 96 4.6 Conclusions and F urther Researc h Directions . . . . . . . . . . . . 97 5 APER TURE MODULA TION . . . . . . . . . . . . . . . . . . . . . . . . 99 5.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.1.1 T raditional Tw o-Phase Approac h to T reatmen t Planning . 99 5.1.2 Ap erture Mo dulation Approac h to T reatmen t Planning . . 100 5.2 Mo del F orm ulation . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.3 A Column Generation Approac h . . . . . . . . . . . . . . . . . . . 105 5.3.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.3.2 The Pricing Problem . . . . . . . . . . . . . . . . . . . . . 106 5.3.3 Solving the Pricing Problem . . . . . . . . . . . . . . . . . 108 v
PAGE 6
5.4 Extending the Mo del with P artial V olume Constrain ts . . . . . . . 116 5.4.1 Mo deling P artial V olume Constrain ts . . . . . . . . . . . . 116 5.4.2 An Enhanced Problem F orm ulation . . . . . . . . . . . . . 118 5.4.3 The Pricing Problem . . . . . . . . . . . . . . . . . . . . . 119 5.5 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . 120 5.6 Concluding Remarks and F uture Researc h . . . . . . . . . . . . . 128 APPENDIXA TONGUE AND GR OO VE CONSTRAINT VIOLA TION . . . . . . . . . 129 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 BIOGRAPHICAL SKETCH . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 vi
PAGE 7
LIST OF T ABLES T able page 2{1 Mo del size and run times. . . . . . . . . . . . . . . . . . . . . . . . . . 32 2{2 V alues of the co ecien ts of the v o xel-based p enalt y functions. . . . . . 33 2{3 V alues of the co ecien ts for the CV aR constrain ts. . . . . . . . . . . 34 2{4 UF OR T automated treatmen t-planning results for ten head-and-nec k cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3{1 T able of results for the neigh b orho o d searc h algorithms. . . . . . . . . 62 4{1 Summary of results of the sequencing algorithms in unconstrained case. 94 4{2 Summary of results of the sequencing algorithms with in terdigitation. 95 4{3 Summary of results of the sequencing algorithms with in terdigitation and T&G correction. . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4{4 Summary of results for maxim um leaf separation. . . . . . . . . . . . . 96 5{1 Mo del sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5{2 V alues of the co ecien ts of the v o xel-based p enalt y functions. . . . . . 123 5{3 V alues of the co ecien ts for the CV aR constrain ts. . . . . . . . . . . 123 5{4 Num b er of ap ertures needed to satisfy eac h criterion. . . . . . . . . . . 125 5{5 Regression results of the relationship b et w een generated and used ap ertures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 vii
PAGE 8
LIST OF FIGURES Figure page 2{1 Irradiating a head-and-nec k cancer patien t using b eams from dieren t directions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2{2 T argets (PTV1, PTV2) and critical structures (RPG: righ t parotid gland, SC: spinal cord). . . . . . . . . . . . . . . . . . . . . . . . . 14 2{3 Examples of p enalt y functions for a target and a critical structure. . . 18 2{4 Examples of Dose-V olume Histograms for a target and a critical structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2{5 T raditional and new (lo w er) D VH constrain ts for a target. . . . . . . 26 2{6 Computation of (upp er) Conditional V alue-at-Risk. . . . . . . . . . . 29 2{7 (a) Cum ulativ e D VHs from the solution for case 5. (b) Ov erla y of axial CT iso dose curv es (45, 54, 73.8 Gy) and structures (PTV1, PTV2, Spinal cord, Righ t parotid gland) for the solution for case 5. 36 2{8 Cum ulativ e D VHs from the solutions obtained using our mo del with and without accoun ting for fractionation. . . . . . . . . . . . . . . . 41 2{9 An axial iso dose distribution from the original case 6, the spatially corrected case 6, and the dierence b et w een the t w o distributions. . 42 2{10 An o v erla y of D VH for case 6 with and without applying spatial p enalties. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3{1 Demonstration of degrees of freedom of an accelerator. . . . . . . . . 47 3{2 The neigh b orho o d of a 4-b eam plan. . . . . . . . . . . . . . . . . . . 57 3{3 Results of equi-spaced BOO for Case 01 for a 4-b eam plan. . . . . . 60 3{4 Results of neigh b orho o d searc h for Case 01 for a 3-b eam plan. . . . . 61 4{1 A m ulti-leaf collimator system . . . . . . . . . . . . . . . . . . . . . . 66 4{2 A single c hannel(ro w) of a m ultileaf collimator with left leaf at p osition l = 3 and righ t leaf at p osition r = n 2 : Radiation P asses through the bixels 3 ; : : : ; n 3 : . . . . . . . . . . . . . . . . . . . . 68 4{3 Realizing a ruence matrix using three shap e matrices . . . . . . . . . 68 viii
PAGE 9
4{4 An example of shap e matrix graph. . . . . . . . . . . . . . . . . . . . 75 4{5 Demonstrating in terdigitation constrain ts. . . . . . . . . . . . . . . . 75 4{6 Example of connectedness constrain t. . . . . . . . . . . . . . . . . . 77 4{7 Connectedness constrain t. . . . . . . . . . . . . . . . . . . . . . . . . 78 4{8 T ongue and gro o v e arrangemen t. . . . . . . . . . . . . . . . . . . . . 78 4{9 An example demonstrating length of the arcs. . . . . . . . . . . . . . 81 4{10 An example of the leaf setup . . . . . . . . . . . . . . . . . . . . . . 91 4{11 Example of a shap e matrix. . . . . . . . . . . . . . . . . . . . . . . . 92 4{12 Demonstrating p ossible w a ys in whic h the matrix can b e split. . . . . 92 4{13 An example of 3 split. . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5{1 Three t yp es of ap ertures. Ap erture (a) can b e deliv ered b y all three considered systems. Ap erture (b) can b e deliv ered b y systems that allo w disconnected ap ertures. Ap erture (c) can only b e deliv ered b y systems allo wing in terdigitation. . . . . . . . . . . . . . . . . . . 111 5{2 Net w ork for the pricing problem in case in terdigitation is not allo w ed. 114 5{3 (a) Ob jectiv e function v alue as a function of n um b er of ap ertures generated and used, and comparison with the optimal v alue to the b eamlet FMO problem. (b) Num b er of ap ertures used as a function of n um b er of ap ertures generated. . . . . . . . . . . . . . . . . . . . . 127 5{4 Co v erage of (a) PTV1 and (b) PTV2 as a function of n um b er of ap ertures used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5{5 Sparing of saliv a glands according to (a) D VH criterion (v olume > 26 Gy) and (b) mean dose as a function of n um b er of ap ertures used. . 127 5{6 Sparing of (a) spinal cord (v olume > 45 Gy) and brainstem (v olume > 54 Gy) and (b) unsp ecied tissue as a function of n um b er of ap ertures used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 A{1 An example of T&G constrain t. . . . . . . . . . . . . . . . . . . . . . 130 ix
PAGE 10
Abstract of Dissertation Presen ted to the Graduate Sc ho ol of the Univ ersit y of Florida in P artial F ulllmen t of the Requiremen ts for the Degree of Do ctor of Philosoph y NO VEL METHODS F OR INTENSITY MODULA TED RADIA TION THERAPY TREA TMENT PLANNING By Arvind Kumar August 2005 Chair: Ra vindra K. Ah uja Co c hair: Hilbrand E. Romeijn Ma jor Departmen t: Industrial and Systems Engineering W e ha v e considered the problem of radiation therap y treatmen t planning for cancer patien ts. Recen t tec hnological adv ancemen ts ha v e led to rapid dev elopmen t and widespread clinical implemen tation of an external-b eam radiation deliv ery tec hnique kno wn as in tensit y mo dulated radiation therap y (IMR T). IMR T allo ws for the creation of v ery complex non-uniform dose distributions that allo w the deliv ery of sucien tly high radiation doses to targets while limiting the radiation dose deliv ered to health y tissues. In IMR T, the patien t is irradiated from sev eral b eams, eac h of whic h is decomp osed in to h undreds of small b eamlets, the in tensities of whic h can b e con trolled individually . W e consider the problem of designing a treatmen t plan (determining in tensit y of eac h b eamlet, also referred as in tensit y mo dulation) for IMR T when the b eam orien tations are giv en. W e ha v e prop osed a new form ulation that incorp orates all asp ects whic h con trol the qualit y of a treatmen t plan that ha v e b een considered to date. Ho w ev er, in con trast with established mixed-in teger and global optimization form ulations, w e do so while retaining linearit y of the optimization problem, and thereb y ensuring that the x
PAGE 11
problem can b e solv ed ecien tly . W e also dev elop ed a linear programming based algorithm to in tegrate b eam orien tation optimization with in tensit y mo dulation. W e ha v e dev elop ed neigh b orho o d searc h metho ds based on fast sensitivit y analysis and greedy searc h based heuristics for solving this in tegrated problem. Deliv ering complex dose distributions sometimes requires a set of v ery complex b eam cross sections (ruence proles). The most widespread hardw are to ac hiev e this is the m ultileaf collimator (MLC), and the most common mo de of op eration for the MLC, called step-and-sho ot deliv ery , is based on a decomp osition of the ruence prole in to ap ertures with asso ciated in tensities. W e dev elop net w ork ro w based algorithms for solving this decomp osition problem. P olynomial-time algorithms are dev elop ed to solv e the decomp osition problem while accoun ting for the constrain ts that the ap ertures m ust satisfy in one or more of the curren tly a v ailable commercial IMR T equipmen t. W e also consider the problem of in tegrating the problem of in tensit y mo dulation with the nal decomp osition step. The problem is form ulated as a large-scale con v ex programming problem, and solv ed via a column generation approac h. xi
PAGE 12
CHAPTER 1 INTR ODUCTION Appro ximately 1.2 million U.S. citizens are newly diagnosed with cancer [ 19 ] ev ery y ear, and more than half of these cancer patien ts are treated b y some form of radiation therap y . Half of the patien ts treated with radiation (appro ximately 300,000 patien ts p er y ear) ma y signican tly b enet from conformal radiation therap y . Despite sophisticated treatmen t, man y patien ts who are initially considered curable do in fact die of their disease. Others ma y suer unin tended side eects from radiation therap y , sometimes sev erely reducing the qualit y of life. These eects are due to radiation therap y treatmen t plans that often deliv er to o little radiation dose to the targets, to o m uc h radiation dose to health y organs, or b oth. Th us, the preserv ation of health y or functional tissues and hence the qualit y of a patien t's life m ust b e balanced against the probabilit y of the eradication of the patien t's disease. In this thesis, w e ha v e prop osed sev eral new linear programming and net w ork ro w based approac hes to designing impro v ed radiation therap y treatmen t plans. During radiation therap y , b eams of radiation pass through a patien t, dep ositing energy along the path of the b eams. This radiation kills b oth cancerous and normal cells. Th us, the radiation therap y treatmen t m ust b e carefully planned, so that a clinically prescrib ed dose is deliv ered to cancerous cells while sparing normal cells in nearb y organs and tissues. T ypically , there are sev eral clinical targets that w e wish to irradiate, and there are sev eral nearb y organs, called critical structures, that w e wish to spare. Note that w e usually treat targets whic h con tain kno wn tumors, as w ell as regions whic h con tain the p ossibilit y of disease spread or accoun t for patien t motion. If w e w ere to treat a patien t with a single b eam of radiation, 1
PAGE 13
2 it migh t b e p ossible to kill all the cells in the targets. Ho w ev er, it w ould also risk damaging normal cells in critical structures lo cated along the path of the b eam. T o a v oid this, b eams are deliv ered from a n um b er of dieren t orien tations spaced around the patien t so that the in tersection of these b eams includes the targets, whic h th us receiv e the highest radiation dose, whereas the critical structures receiv e dose from some, but not all, b eams and can th us b e spared. Conformal radiation therap y seeks to conform the geometric shap e of the deliv ered radiation dose as closely as p ossible to that of the in tended targets. In con v en tional conformal radiation therap y , this usually means that from eac h b eam direction w e deliv er a single b eam with uniform in tensit y lev el whose shap e conforms to the b eam's-ey e view (t w o-dimensional pro jection) of the targets in the patien t as seen from that b eam. Recen t tec hnological adv ancemen ts ha v e led to rapid dev elopmen t and widespread clinical implemen tation of an external-b eam radiation deliv ery tec hnique kno wn as in tensit y mo dulated radiation therap y (IMR T). IMR T allo ws for the creation of v ery complex non uniform dose distributions that allo w the deliv ery of sucien tly high radiation doses to targets while limiting the radiation dose deliv ered to health y tissues. These dose distributions are obtained b y dynamically blo c king dieren t parts of the b eam. The application of optimization tec hniques is essen tial to enable ph ysicians, and thereb y patien ts, to fully b enet from the added rexibilit y that this tec hnique promises. 1.1 Dissertation Summary There are sev eral asp ects of IMR T treatmen t planning where optimization ma y b e applied. These asp ects include the selection of b eam n um b er, selection of b eam orien tation, selection of b eam qualit y , selection of the ruence distributions and their decomp osition in to deliv erable sequences. Ideally , w e w ould lik e to build one in tegrated mo del that w ould consider all the goals and constrain ts of the problem sim ultaneously and th us whose solution w ould pro vide a true
PAGE 14
3 globally optimal treatmen t plan. Though, w e migh t succeed in form ulating suc h a mo del, it has an astronomical dimensionalit y and cannot b e practically solv ed b y an y curren tly feasible optimization algorithm. Therefore, in practice, this optimization problem is t ypically brok en in to three subproblems that can b e solv ed or heuristically appro ximated. When solving the IMR T treatmen t planning problem, the optimizations or decisions to b e made are t ypically divided as suc h: (i) determine the n um b er and orien tations at whic h radiation b eams will b e deliv ered; (ii) determine the ruence map for eac h radiation b eam selected in (i); and (iii) determine a metho d of discretization of the ruence map pro duced in (ii) to pro duce a deliv erable treatmen t plan. 1.1.1 Fluence Map Optimization The problem of designing an optimal radiation in tensit y prole in the patien t is often referred to as the ruence map optimization (FMO) problem. The goal of this problem is to design a radiation therap y treatmen t plan that deliv ers a sp ecic lev el of radiation dose, a so-called prescription dose, to the targets. A t the same time to spare critical structures, it should b e ensured that the lev el of radiation dose receiv ed b y these structures do es not exceed some structure-sp ecic tolerance dose. These t w o goals are inheren tly con tradictory if the targets are lo cated in the vicinit y of critical structures. F or example, tumors in the head-and-nec k area are often lo cated v ery close to critical structures lik e spinal cord, brainstem, and saliv ary glands. A common approac h is to searc h for a radiation therap y treatmen t plan that satises the prescription and tolerance dose requiremen ts to the largest exten t p ossible. In this thesis, w e ha v e prop osed a new form ulation that incorp orates all asp ects that con trol the qualit y of a treatmen t plan that ha v e b een considered to date. Ho w ev er, unlik e mixed-in teger and global optimization form ulations in the literature, our form ulation is linear, thereb y ensuring that the problem
PAGE 15
4 can b e solv ed ecien tly . W e presen t a no v el linear programming (LP) based approac h for ecien tly solving the IMR T FMO problem to global optimalit y . Our mo del o v ercomes the apparen t limitations of a linear-programming approac h b y appro ximating an y con v ex ob jectiv e function b y a piecewise linear con v ex function. This approac h allo ws us to retain the rexibilit y oered b y general con v ex ob jectiv e functions, while allo wing us to form ulate the FMO problem as a LP problem. In addition, a no v el t yp e of partial-v olume constrain t that b ounds the tail a v erages of the dieren tial dose-v olume histograms of structures is imp osed while retaining linearit y as an alternativ e approac h to impro v e dose homogeneit y in the target v olumes, and to attempt to spare as man y critical structures as p ossible. The goal of this w ork is to dev elop a v ery rapid global optimization approac h that nds high qualit y dose distributions. Implemen tation of this mo del has demonstrated excellen t results. W e found globally optimal solutions for eigh t 7-b eam headand-nec k cases in less than 1 min ute of computational time on a single pro cessor p ersonal computer without the use of partial-v olume constrain ts. Adding suc h constrain ts increased the running times b y a factor of 2-3, but impro v ed the sparing of critical structures. All cases demonstrated excellen t target co v erage, target homogeneit y and organ sparing. 1.1.2 Beam Orien tation Optimization P atien ts receiving radiation therap y are t ypically treated on a clinical radiation deliv ery device called a linear accelerator. The radiation source is con tained in the accelerator head and can b e view ed as a p oin t source of high-energy photons. The patien t is immobilized with restrain ts on a couc h, and the accelerator head can rotate for a full 360 ab out a p oin t, called \iso cen ter", lo cated one meter from the source. In practice, it is common to use b et w een 3 and 9 accelerator b eam orien tations. Although additional degrees of freedom in the mo v emen t of the couc h
PAGE 16
5 as w ell as the accelerator head are a v ailable, only a single, xed couc h lo cation is used in practice. T ypical commercially a v ailable IMR T treatmen t planning systems select the b eam orien tations in an ad-ho c manner, and then optimization is done to determine ruence proles. Dev eloping an algorithmic approac h that sim ultaneously determines optimal b eam orien tations and ruences will b e of substan tial v alue to clinical IMR T practices. In this thesis, w e describ e a heuristic approac h that sim ultaneously addresses BOO and FMO. Our approac h starts with some b eam orien tations and rep eatedly impro v es it b y adding and dropping new b eams. A no v el feature of our approac h is that w e determine the optimal b eam in tensities b y solving a linear programming problem and use the sensitivit y analysis information of this problem to iden tify b eam orien tations to b e added. Our computational results indicate that the metho ds dev elop ed can appro ximate the b eha vior of the solution space v ery precisely and th us can guide the searc h in the correct direction. 1.1.3 T reatmen t Deliv ery IMR T allo ws for the creation of v ery complex non-uniform dose distributions that facilitates the deliv ery of sucien tly high radiation doses to targets while limiting the radiation dose deliv ered to health y tissues. Ho w ev er, IMR T often requires v ery complex dose proles. Because of v ery complex proles in IMR T, it is no longer p ossible to deliv er these doses using st yrofoam blo c ks (the con v en tional metho d) as it w ould mak e the pro cess v ery inecien t. Computer-con trolled Multileaf Collimators (MLCs) are th us in tro duced to deliv er IMR T ecien tly . The problem of con v erting the dose proles in to a set of leaf sequence les that con trol the mo v emen t of the MLC during deliv ery is referred to as the MLC sequencing problem. T o deliv er a desired dose prole, the ruence maps are decomp osed in to ap ertures (or shap e matrices) with asso ciated in tensities (or b eam-on time). The MLC sequencing problem is to nd a decomp osition whic h minimizes the total
PAGE 17
6 n um b er of ap ertures used and the sum of the corresp onding in tensities. The latter problem is also referred as the total minim um BOT problem. Because of ph ysical limitations, the MLCs cannot deliv er an y random pattern generated b y the sequencing algorithm. In our w ork, w e ha v e prop osed p olynomial time algorithms for solving the BOT problem while accoun ting for these constrain ts. In particular, w e pro v e that BOT is p olynomially solv able ev en with tongue and gro o v e constrain ts (a correction whic h is needed to mak e sure that the area b et w een lea v es in MLC is not seriously underdosed b ecause of mec hanical arrangemen t). The problem of determining minim um n um b er of ap ertures to deliv er a ruence prole has b een sho wn to b e NP-hard; hence, it is unlik ely to b e solv able in p olynomial time. Th us, heuristics ha v e b een prop osed for the sequencing problem whic h minimizes the n um b er of ap ertures used while also k eeping a tab on total BOT to o. The heuristics curren tly used in practice are based on iterativ ely extracting ap ertures from a giv en in tensit y matrix in an ad-ho c manner. As a result, the ap ertures extracted migh t not satisfy all the constrain ts. In that case, the ap ertures are c hanged iterativ ely to nd a solution whic h do es not violate an y constrain ts. Since the constrain t violation is considered after the ap erture extraction, these t w o steps are not in tegrated, and this ma y yield v ery inecien t shap es. W e presen t a net w ork-based form ulation to in tegrate the extraction pro cess with constrain t enforcemen t. W e also dev elop v ery ecien t metho ds to p erform the extraction pro cess on these net w orks. This form ulation can em ulate most of the existing metho ds to solv e the MLC sequencing problem and can also p oten tially impro v e them. Our form ulation can th us serv e as a unied platform to compare the dieren t form ulations in the literature. W e ha v e also prop osed a no v el metho d to extract the ap ertures. Our metho d outp erforms the b est rep orted metho d b oth in
PAGE 18
7 terms of BOT and n um b er of set-ups used. The run time of our algorithm is on the order of a fraction of a second. 1.1.4 Ap erture Mo dulation The MLC-based IMR T systems cannot directly deliv er the ruence maps determined b y the FMO routine. Th us, these maps need to b e discretized and subsequen tly decomp osed b efore deliv ery . Ho w ev er, the deterioration in the solution qualit y of the plan on discretization is inevitable. This m ulti-stage approac h is also rigid in the sense that a sound trade-o b et w een treatmen t time and treatmen t plan qualit y cannot b e made due to the decoupling of the decomp osition phase from the treatmen t plan optimization phase. Some empirical evidence suggests that one can obtain high-qualit y treatmen t plans with only a limited n um b er of ap ertures. Recen tly , heuristic lo cal searc h metho ds and mixedin teger linear programming mo dels for limited searc h space ha v e b een prop osed to address this issue. These approac hes ha v e b een applied to head-and-nec k and prostate cancer cases. Ho w ev er, they are all heuristic in nature, and therefore an optimal solution cannot b e guaran teed. The k ey to obtaining the b est ap ertures seems to b e to in tegrate the FMO and decomp osition problems formally . W e address the problem of designing a treatmen t plan for IMR T that, unlik e established t w o-phase approac hes, sim ultaneously determines an optimal set of the shap es (or ap ertures) and their corresp onding in tensities. The problem is form ulated as a large-scale con v ex programming problem, and a column generation approac h is dev elop ed to deal with the dimensionalit y of the problem. The asso ciated pricing problem determines, in eac h iteration, one or more ap ertures to b e added to our problem. Sev eral v arian ts of this pricing problem are discussed, eac h corresp onding to a particular set of constrain ts that the ap ertures m ust satisfy in one or more of the curren tly a v ailable commercial IMR T equipmen t. P olynomial-time algorithms for solving eac h of these v arian ts of the pricing problem
PAGE 19
8 to optimalit y are presen ted. The eectiv eness of our approac h is demonstrated on clinical data. 1.2 Con tribution Summary 1.2.1 Fluence Map Optimization 1. W e presen t a no v el linear programming based approac h for ecien tly solving the IMR T FMO problem to global optimalit y . The mo del o v ercomes the apparen t limitations of linear ob jectiv e function b y appro ximating con v ex functions b y piecewise linear functions. 2. W e ha v e also prop osed a set of no v el linear dose-v olume histogram constrain ts. Unlik e con v en tional dose-v olume constrain ts, these constrain ts also con trol the spread of dose distribution. These constrain ts are linear in nature and can b e included in the mo del without adding m uc h complexit y . 3. W e tested our metho dology on eigh t clinical cases of head-and-nec k cancer patien ts. Our results demonstrate that the linear mo del dev elop ed can translate the ph ysicist's ob jectiv es v ery precisely . Av erage solution times for the cases are under a min ute. 1.2.2 Beam Orien tation Optimization 1. W e suggest an exact approac h to solv e the limited BOO problem that exhaustiv ely en umerates all p ossible orien tations of equi-spaced b eams and solv es a linear programming problem for eac h orien tation. The b est solution obtained among these solutions is a global optimal solution for the case of equi-spaced b eam orien tations. 2. W e dev elop a neigh b orho o d searc h based algorithm for in tegrating FMO with BOO for the case of non-equispaced b eam orien tations. 3. W e presen t extensiv e computational results of our algorithm. The computational results indicate that the algorithms dev elop ed can appro ximate the BOO problem v ery w ell.
PAGE 20
9 1.2.3 T reatmen t Deliv ery 1. W e consider the problem of decomp osing a giv en ruence map in deliv erable sequences. Most of the practical deliv ery constrain ts are considered in the form ulation. 2. W e pro v e that minim um BOT problem with or without T ongue & Gro o v e constrain ts is p olynomially solv able. 3. W e dev elop a net w ork based algorithm to solv e the BOT problem ecien tly . 4. W e dev elop a net w ork based heuristic for the total sequencing problem (BOT minimization and n um b er of shap e minimization). Our results demonstrate that our algorithm outp erforms the b est rep orted metho d in literature. 1.2.4 Ap erture Mo dulation 1. W e consider the problem of in tegrating the ruence map generation problem with the problem of decomp osing the ruence map in deliv erable sequences. W e ha v e dev elop ed column generation based metho ds to tac kle this problem ecien tly . 2. W e dev elop net w ork based mo dels to incorp orate imp ortan t deliv ery constrain ts. 3. W e p erform an extensiv e computational analysis of our results. Our results demonstrate that n um b er of ap ertures used can b e signican tly reduced b y in tegrating the pro cess of map generation and decomp osition. Run time of the instances is of the order of 15 min utes.
PAGE 21
CHAPTER 2 FLUENCE MAP OPTIMIZA TION 2.1 In tro duction During radiation therap y , b eams of radiation pass through a patien t, dep ositing energy along the path of the b eams. This radiation kills b oth cancerous and normal cells. Th us, the radiation therap y treatmen t m ust b e carefully planned, so that a clinically prescrib ed dose is deliv ered to cancerous cells while sparing normal cells in nearb y organs and tissues. T ypically , there are sev eral clinical tar gets that w e wish to irradiate, and there are sev eral nearb y organs, called critic al structur es , that w e wish to spare. Note that w e usually treat targets whic h con tain kno wn tumors, as w ell as regions whic h con tain the p ossibilit y of disease spread or accoun t for patien t motion. If w e w ere to treat a patien t with a single b eam of radiation, it migh t b e p ossible to kill all the cells in the targets. Ho w ev er, it w ould also risk damaging normal cells in critical structures lo cated along the path of the b eam. T o a v oid this, b eams are deliv ered from a n um b er of dieren t orien tations spaced around the patien t so that the in tersection of these b eams includes the targets, whic h th us receiv e the highest radiation dose, whereas the critical structures receiv e dose from some, but not all, b eams and can th us b e spared (see Fig. 2{1 ). Conformal radiation therap y seeks to conform the geometric shap e of the deliv ered radiation dose as closely as p ossible to that of the in tended targets. In con v en tional conformal radiation therap y , this usually means that from eac h b eam direction w e deliv er a single b eam with uniform in tensit y lev el whose shap e conforms to the b eam's-ey e view of the targets in the patien t as seen from that b eam. Recen t tec hnological adv ancemen ts ha v e led to rapid dev elopmen t and widespread clinical implemen tation of an external-b eam radiation deliv ery 10
PAGE 22
11 Figure 2{1: Irradiating a head-and-nec k cancer patien t using b eams from dieren t directions.tec hnique kno wn as intensity mo dulate d r adiation ther apy (IMR T) (see W ebb [ 81 ]). IMR T allo ws for the creation of v ery complex non-uniform dose distributions that allo w the deliv ery of sucien tly high radiation doses to targets while limiting the radiation dose deliv ered to health y tissues. These dose distributions are obtained b y dynamic al ly blo c king dieren t parts of the b eam. The application of optimization tec hniques is essen tial to enable ph ysicians, and thereb y patien ts, to fully b enet from the added rexibilit y that this tec hnique promises. P atien ts receiving radiation therap y are t ypically treated on a clinical radiation deliv ery device called a line ar ac c eler ator . The radiation source is con tained in the accelerator head and can b e view ed as a p oin t source of high-energy photons. The patien t is immobilized with restrain ts on a couc h, and the accelerator head can rotate for a full 360 ab out a p oin t that is 1 meter from the source, called the iso c enter . An IMR T capable accelerator is equipp ed with a so-called m ulti-leaf collimator system, whic h is used to generate the desired dose distribution at a giv en angle. Although additional degrees of freedom in the mo v emen t of the couc h as w ell
PAGE 23
12 as the accelerator head are a v ailable, it is common in practice to use only a single, xed couc h lo cation, and consider a v ery limited n um b er (3-9) of lo cations of the accelerator head in the circle around the iso cen ter. Giv en the lo cations of the accelerator head, the problem of designing an optimal radiation densit y prole in the patien t is often referred to as the ruenc e map optimization (FMO) pr oblem . The goal of this problem is to design a radiation therap y treatmen t plan that deliv ers a sp ecic lev el of radiation dose, a so-called pr escription dose , to the targets, while on the other hand sparing critical structures b y ensuring that the lev el of radiation dose receiv ed b y these structures do es not exceed some structure-sp ecic toler anc e dose . These t w o goals are inheren tly con tradictory if the targets are lo cated in the vicinit y of critical structures. This is esp ecially problematic for certain cancers, suc h as tumors in the head-and-nec k area, whic h are often lo cated v ery close to for instance the spinal cord, brainstem, and saliv ary glands. A common approac h is to searc h for a radiation therap y treatmen t plan that satises the prescription and tolerance dose requiremen ts to the largest exten t p ossible. All commercially a v ailable IMR T treatmen t planning systems, as w ell as most of the researc h in the medical ph ysics literature to date, use lo cal searc h (suc h as the conjugate gradien t metho d; see, e.g., Shepard et al . [ 70 ], Xing and Chen [ 87 ], Holmes and Mac kie [ 34 ], Bortfeld et al . [ 12 ]) or sim ulated annealing tec hniques (see, e.g., Rosen et al . [ 62 ], Mageras and Mohan [ 43 ], W ebb [ 78 , 79 ]) to nd a satisfactory treatmen t plan. In a recen t comprehensiv e review of the radiation therap y literature, Shepard et al . [ 69 ] surv ey ed man y tec hniques previously used, including some simple linear and con v ex programming form ulations of the problem. More recen t op erations researc h approac hes are m ulti-criteria optimization (see, e.g., Hamac her and K ufer [ 32 ]; Bortfeld et al . [ 16 ]), and mixed-in teger linear programming (see, e.g., Lee et al . [ 40 , 41 ]).
PAGE 24
13 In this dissertation, w e prop ose a new approac h to the FMO problem. The strength of our approac h is that w e are able to enric h existing mo dels b y incorp orating common measures of treatmen t plan qualit y while retaining linearit y , and thereb y ecien t solv abilit y , of the problem. The eciency of linear programming, in turn, allo ws us to consider extended mo dels that incorp orate additional asp ects of the problem that ha v e b een heretofore ignored in the IMR T literature. This c hapter is organized as follo ws. In Section 2.2 , w e will formalize the radiation treatmen t planning problem and discuss v arious mo deling issues. Then, in Section 2.3 , w e will presen t our form ulation of the FMO problem, whic h balances the trade-o b et w een radiation of the targets and critical structures, as a linear programming problem. In Section 2.4 , w e discuss the results of our approac h on ten clinical data sets, and compare these results with an op erational planning system that is widely used. Extensions of the mo del are discussed in Section 2.5 , and w e end the c hapter with some remarks on the clinical imp ortance of our mo dels and results in Section 2.6 . 2.2 Notation and Dose Calculation The targets and critical structures (together simply referred to as structur es ) are irradiated using a predetermined set of b e ams , eac h corresp onding to a particular b eam angle. F or eac h b eam angle, w e ha v e a rectangular or square ap erture, t ypically of size 10 10 cm 2 to 20 20 cm 2 . Eac h of these ap ertures is decomp osed or discretized in to small b e amlets , t ypically of size 1 1 cm 2 , yielding on the order of 100-400 b eamlets p er b eam. A 2-dimensional matrix is used to represen t the p osition and in tensit y of all b eamlets in a b eam. Eac h elemen t in this matrix is called a bixel , and its v alue represen ts the intensity (or, more correctly , ruence; see W ebb and Lomax [ 82 ]) of the corresp onding b eamlet. W e denote the n um b er of bixels con tained in all b eams b y N . The core task in IMR T treatmen t planning is to nd radiation in tensities (sometimes also called w eigh ts) for all b eamlets. W e
PAGE 25
14 will denote the decision v ariable represen ting the w eigh t of b eamlet i b y x i , and the v ector of all b eamlet w eigh ts b y x 2 R N+ . W e denote the total n um b er of structures b y S , and the n um b er of targets b y T . F or con v enience, w e assume that the targets are indexed b y s = 1 ; : : : ; T , and the critical structures b y s = T + 1 ; : : : ; S . F urthermore, eac h of the structures s is discretized in to a nite n um b er v s of cub es, t ypically of size 4 4 4 mm 3 . In the con text of radiation therap y , these cub es are usually called voxels . The n um b er of v o xels required to accurately represen t the targets and surrounding critical structures is t ypically around 100,000. Note that, in general, structures ma y o v erlap. F or instance, if a target has in v aded a critical structure, some v o xels will b e in b oth the target and the critical structure (see Fig. 2{2 for an example). Ho w ev er, in our mo del w e asso ciate a unique (i.e., dominan t) structure with eac h v o xel, based on a priorit y list of all structures (where targets usually ha v e the highest priorities, follo w ed b y the critical structures). Note that this is not carried o v er to our treatmen t planning system, where the obtained treatmen t plans are alw a ys appropriately ev aluated using the original, o v erlapping, structures. As a Figure 2{2: T argets (PTV1, PTV2) and critical structures (RPG: righ t parotid gland, SC: spinal cord).
PAGE 26
15 b eamlet of radiation passes through a patien t, it dep osits an amoun t of energy along its path that decreases with the distance tra v eled. Eac h b eamlet therefore aects a large n um b er of v o xels in the patien t. Clearly , eac h v o xel receiv es dose from sev eral b eamlets. Although not en tirely accurate, it is common practice to represen t the relationship b et w een the b eamlet in tensities and the dose receiv ed b y v o xels in a linear w a y . Letting D ij s denote the dose receiv ed b y v o xel j in structure s from b eamlet i at unit in tensit y , w e can express the dose receiv ed b y eac h v o xel as a function of the b eamlet in tensities x as follo ws: D j s ( x ) = N X i =1 D ij s x i j = 1 ; : : : ; v s ; s = 1 ; : : : ; S: 2.3 Mo del F orm ulation The metho ds describ ed in this c hapter ha v e b een published in Romeijn et al . [ 59 ]. 2.3.1 A Basic Linear Programming F orm ulation The basic constrain ts in a FMO problem are hard b ounds on the dose receiv ed b y the structure of the follo wing form: (F 1 ) Deliv er a minim um prescription dose, sa y L , to all v o xels in the targets (where v o xels in dieren t targets ma y ha v e dieren t prescription doses). (F 2 ) Do not deliv er more than a maxim um tolerance dose, sa y U , to all v o xels in b oth targets and critical structures (where v o xels in dieren t structures ma y ha v e dieren t tolerance doses). Constrain ts (F 1 ) ensure a high lik eliho o d of eradicating the disease, while constrain ts (F 2 ) for the critical structures ensure a high lik eliho o d that the functionalit y of these structures is retained. In addition, constrain ts (F 2 ) for the targets represen t the undesirabilit y of deliv ering v ery high doses (or v ery non-uniform dose distributions) to targets. Constrain ts of this t yp e are often referred to as ful l volume c onstr aints , since they need to b e satised ev erywhere in a particular
PAGE 27
16 structure. Suc h constrain ts are particularly imp ortan t for targets, as w ell as for certain structures suc h as the spinal cord or brainstem, where ev en exp osing a v ery small v olume in the structure to a dose ab o v e the upp er b ound can cause a loss of the functionalit y of the en tire structure. The full v olume constrain ts that need to b e satised b y a treatmen t plan, c haracterized b y the b eamlet in tensities x , can easily b e expressed as follo ws D j s ( x ) L s j = 1 ; : : : ; v s ; s = 1 ; : : : ; T D j s ( x ) U s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S where L s and U s denote the lo w erb ound and upp erb ound, resp ectiv ely , on the dose receiv ed b y all v o xels in structure s . In practice, not all plans satisfying these constrain ts are equally desirable. Unfortunately , no satisfactory metho d or scien tic basis for distinguishing b et w een dieren t solutions has b een established. Th us, there is no ob vious fundamen tal or a priori form for the ob jectiv e function. W e will follo w a widely used approac h, in whic h w e dene, for eac h structure, a p enalt y function of the dose receiv ed b y the v o xels in that structure. Denoting the p enalt y function for structure s b y F s w e obtain the follo wing ob jectiv e function to b e minimized S X s =1 1 v s v s X j =1 F s ( z j s ) : The scaling term 1 =v s is included to ensure that the p enalt y functions are insensitiv e to the relativ e sizes of the structures found in dieren t problem instances. T ogether with the full v olume constrain ts, this leads to our basic mo del for FMO minimize S X s =1 1 v s v s X j =1 F s ( z j s )
PAGE 28
17 sub ject to (FMO 1 ) N X i =1 D ij s x i = z j s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S z j s L s j = 1 ; : : : ; v s ; s = 1 ; : : : ; T z j s U s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S x i 0 i = 1 ; : : : ; N z j s free j = 1 ; : : : ; v s ; s = 1 ; : : : ; S: T ypical p enalt y functions that ha v e b een prop osed in the literature ha v e a v alue prop ortional to the dose receiv ed b y all v o xels, or only b y v o xels in critical structures (see Rosen et al . [ 63 ]). A more sophisticated ob jectiv e function is obtained when the p enalt y is equal to the absolute deviation of the dose receiv ed from a prescription or tolerance dose (see Shepard et al . [ 69 ]). Since large deviations from the prescription or tolerance dose are considered to b e m uc h more imp ortan t than small deviations, a m uc h more frequen tly used alternativ e is to use a w eigh ted least squares ob jectiv e function. Since lo w er doses in critical structures are alw a ys preferred to higher doses, for v o xels in these structures the p enalt y function is generally one-sided, i.e., only surpluses o v er the tolerance dose are p enalized in the ob jectiv e function. W e prop ose to write the p enalt y function for structure s as the sum of t w o functions F s ( z ) = F U s ( z ) + F O s ( z ) : The former term accoun ts for the underdose p enalt y , and the latter term accoun ts for the o v erdose p enalt y . W e usually start b y c ho osing the t w o terms to b e onesided p olynomials with resp ect to an underdose threshold T U s and an o v erdose
PAGE 29
18 threshold T O s , resp ectiv ely F U s ( z ) = U s max (0 ; T U s z ) p Us F O s ( z ) = O s max(0 ; z T O s ) p Os where U s ; O s 0 and p Us ; p Os 1 to ensure con v exit y of the p enalt y function. T ypically , for critical structures w e will ha v e U s = 0, represen ting that only v o xels that are o v erdosed (with resp ect to the o v erdose threshold) are p enalized. Note that in our standard form ulation w e ha v e assumed that the p enalt y functions are the same for all v o xels in a particular structure. In general, w e ma y of course ev en relax this assumption b y allo wing for an individual p enalt y function for eac h v o xel (see, for example, our extended mo del in Section 2.5.2 ). In Fig. 2{3 , w e illustrate some examples of p enalt y functions. Figure 2{3: Examples of p enalt y functions for a target and a critical structure. T o allo w us to form ulate the FMO problem as a purely linear programming problem, w e appro ximate the resulting p enalt y functions b y piecewise linear con v ex functions (as illustrated in Fig. 2{3 ). With suc h p enalt y functions, w e can reform ulate (FMO 1 ) as a purely linear programming problem using one of t w o
PAGE 30
19 approac hes. In the rst approac h (see, e.g., Murt y [ 44 ]), w e asso ciate a decision v ariable, sa y z k j s , with the k th linear segmen t (out of a total of K s segmen ts) in the piecewise linear con v ex p enalt y function for eac h v o xel ( j; s ). Eac h of these v ariables has a zero lo w er b ound and an upp er b ound giv en b y the size of the segmen t. The dose to a v o xel is then giv en b y the sum of all segmen t v ariables for that v o xel. Letting a ks denote the slop e of the k th segmen t for structure s , and w k s the width of that segmen t, w e add the follo wing decision v ariables and constrain ts to (FMO 1 ) K s X k =1 z k j s = z j s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S z k j s w k s k = 1 ; : : : ; K s ; j = 1 ; : : : ; v s ; s = 1 ; : : : ; S z k j s 0 k = 1 ; : : : ; K s ; j = 1 ; : : : ; v s ; s = 1 ; : : : ; S and write the ob jectiv e as minimize S X s =1 1 v s v s X j =1 K s X k =1 a kj s z k j s : Note that the con v exit y of the p enalt y function implies that the slop es of the segmen ts are increasing as the dose increases, so that an y linear programming solution will use the segmen ts in order, from left to righ t, so that the p enalties are correctly determined. The second approac h is based on the observ ation that the piecewise linear con v ex p enalt y function ma y b e written as the maxim um of a set of linear p enalt y functions, sa y F s ( z j s ) = max k =1 ;::: ;K 0 s b ks z j s + c ks :
PAGE 31
20 W e then in tro duce a decision v ariable, sa y t j s , for eac h v o xel ( j; s ) that represen ts the p enalt y for that v o xel, whic h leads to a simple linear ob jectiv e minimize S X s =1 v s X j =1 t j s : W e then ensure that the correct p enalt y is applied to eac h v o xel b y adding constrain ts that b ound the p enalt y v ariables from ab o v e b y eac h of the linear functions. More formally , w e add the follo wing decision v ariables and constrain ts to the mo del b ks z j s + c ks t j s k = 1 ; : : : ; K 0 s ; j = 1 ; : : : ; v s ; s = 1 ; : : : ; S t j s 0 j = 1 ; : : : ; v s ; s = 1 ; : : : ; S: T o giv e an impression of the dimensions of this problem, w e study the n um b er of decision v ariables and constrain ts for a t ypical problem instance when emplo ying the rst approac h for implemen ting the piecewise linear ob jectiv e functions, as w e ha v e done in this study . If w e form ulate problem (FMO 1 ) using piecewise linear p enalt y functions with k O segmen ts for o v erdosing and k U segmen ts for underdosing, w e will ha v e a total of N + k U T X s =1 v s + k O S X s =1 v s in tensit y and dose v ariables. The mo del will ha v e T X s =1 v s + 2 S X s =1 v s constrain ts. In a t ypical problem instance, w e w ould ha v e ab out N = 1 ; 0002 ; 000 b eamlet in tensities, ab out P Ts =1 mv s = 10 ; 000 v o xels in targets, and ab out P Ss =1 v s = 100 ; 000 v o xels in all structures com bined. F urthermore, w e usually emplo y k U = 4 segmen ts for underdosing and k O = 2 segmen ts for o v erdosing in the piecewise linear ob jectiv e function. This w ould lead to a linear programming
PAGE 32
21 problem with appro ximately 240,000 decision v ariables and 220,000 constrain ts. Often, a reduction in the n um b er of v o xels used to represen t critical structures can signican tly reduce the problem dimensions without signican tly aecting the solutions. Ho w ev er, ev en when the n um b er of v o xels in all structures com bined is reduced to sa y P Ss =1 v s = 20 ; 000, w e will still obtain a problem with appro ximately 80,000 decision v ariables and 50,000 constrain ts. 2.3.2 The Dose-V olume Histogram Unfortunately , in most practical situations there do es not exist a feasible treatmen t plan as dened ab o v e, as it is usually imp ossible to satisfy the full v olume constrain ts for eac h v o xel in the target v olumes as w ell as all critical structures. There will b e some v o xels in the target v olumes that ha v e to b e underdosed and some v o xels in the critical structures that ha v e to b e o v erdosed. The task of the ph ysician is to someho w mak e a sound trade-o b et w een violations of conricting b ound constrain ts. The to ol that is most commonly emplo y ed b y ph ysicians to mak e suc h trade-os and judge the qualit y of a treatmen t plan is the so-called (cum ulativ e) Dose-V olume Histo gr am (D VH) . This histogram sp ecies, for a giv en target or critical structure, the fraction of its v olume that receiv es at least a certain amoun t of dose. More formally , for a giv en structure s , suc h a histogram is a nonincreasing function H s : R + ! [0 ; 1] where H s ( d ) is dened as the relativ e v olume of the part of a structure that receiv es d units of dose or more. Clearly , H s (0) = 1 and lim d !1 H s ( d ) = 0. Using the nite represen tation of all structures using v o xels, the D VH for structure s under b eamlet in tensities x can b e expressed as H s ( d ; x ) = jf j = 1 ; : : : ; v s : D j s ( x ) d gj v s
PAGE 33
22 i.e., the v alue of the D VH at dose d is precisely the fraction of v o xels in the structure who receiv e a dose of at least d units. F or instance, in Fig. 2{4 b elo w, the p oin t on the D VH for the target (indicated b y the righ t v ertical dotted line) indicates that appro ximately 95% of the target v olume receiv es at least 73.8 units of dose (usually measured in Gra y (Gy)). Similarly , the p oin t on the D VH for the critical structure (indicated b y the left v ertical dotted line) indicates that appro ximately 26% of the critical structure receiv es in at least 30 Gy . In the remainder of this c hapter, a generic D VH will b e denoted b y the function H ( d ; x ). Figure 2{4: Examples of Dose-V olume Histograms for a target and a critical structure. W e can no w enhance our mo del for FMO b y including constrain ts that con trol the shap e and lo cation of the D VHs for all structures. In the next section, w e will discuss ho w to mo del suc h constrain ts. 2.3.3 Constrain ts on the Dose-V olume Histogram Constrain ts on the shap e and lo cation of a D VH are often referred to as p artial volume c onstr aints . Suc h constrain ts can b e view ed as a t yp e of soft b ound constrain ts on the dose receiv ed b y a structure, since the corresp onding b ounds only need to b e satised b y a fraction of the v o xels in a particular structure.
PAGE 34
23 T raditional partial v olume constrain ts . T raditionally , optimization mo dels ha v e b een form ulated that incorp orate constrain ts on the shap e and lo cation of the D VH of a structure of the follo wing form: (P 1 ) A t least a fraction of the v o xels in a target receiv e a dose of at least L : H ( L ; x ) : (P 2 ) A t least a fraction of the v o xels in a structure receiv e a dose of at most U : H ( U ; x ) : Clearly , the full v olume constrain ts incorp orated in our basic mo del (FMO 1 ) can b e view ed as constrain ts on the D VHs of the structures b y c ho osing = 1. T o date, there ha v e b een t w o main approac hes to deal with traditional partial v olume constrain ts. The rst (heuristic) approac h is to deal with suc h constrain ts b y adding an appropriate nonlinear p enalt y term to the ob jectiv e function that p enalizes solutions that do not satisfy the partial v olume constrain ts. Suc h heuristic approac hes ha v e b een prop osed b y Bortfeld et al . [ 15 ], and ha v e b een incorp orated in to quadratic programming-based algorithms for the FMO problem (see, e.g., Oelfk e and Bortfeld [ 45 ], W u and Mohan [ 83 ], Cho et al . [ 24 ], Spirou and Ch ui [ 74 ], Deasy [ 30 ]), or a sim ulated annealing algorithms (see, e.g., Carol et al . [ 20 ], Cattaneo et al . [ 21 ], Cho et al . [ 23 ]). The resulting ob jectiv e function is noncon v ex, and the qualit y of the treatmen t plans obtained b y heuristics that ha v e incorp orated this approac h to handling partial v olume constrain ts has not b een established. The second (exact) approac h is to handle the nonlinearit y of the constrain ts b y in tro ducing a binary v ariable, sa y y j s , for eac h v o xel ( j; s ) and eac h partial v olume constrain t, whic h represen ts whether (v alue 1) or not (v alue 0) this v o xel receiv es at least L s (in the case of a lo w er partial v olume constrain t on a target)
PAGE 35
24 or at most U s units of dose (in the case of an upp er partial v olume constrain t on a target). This approac h results in a v ery large-scale mixed in teger programming problem (see, e.g., Langer et al . [ 37 ] and Lee et al . [ 40 , 41 ]). In particular, a single lo w er partial v olume constrain t for some target structure s can b e incorp orated in (FMO 1 ) b y adding z j s L s y j s j = 1 ; : : : ; v s v s X j =1 y j s d v s e y j s 2 f 0 ; 1 g j = 1 ; : : : ; v s for one or more suitable v alues of . Similarly , a single lo w er partial v olume constrain t for some structure s can b e incorp orated in (FMO 1 ) b y adding z j s U s y j s + M (1 y j s ) j = 1 ; : : : ; v s v s X j =1 y j s d v s e y j s 2 f 0 ; 1 g j = 1 ; : : : ; v s for one or more suitable v alues of . Note that eac h partial v olume constrain t imp osed on a structure requires the in tro duction of a set of binary v ariables and corresp onding constrain ts as giv en ab o v e. Both the heuristic and the exact approac h signican tly c hange the nature and the dicult y of the optimization problem, to a nonlinear glob al optimization problem or a v ery large-scale mixe d-inte ger programming problem, resp ectiv ely . Both these approac hes require large amoun ts of computation time for realistically sized problems, prohibiting a real-time solution of clinically relev an t mo dels. F or example, using the t ypical problem dimensions men tioned in Section 2.3.1 , w e w ould need 10,000 binary v ariables for e ach traditional partial v olume constrain t applied to the targets. As alternativ es for the partial v olume constrain ts, constrain ts that
PAGE 36
25 limit the v ariabilit y in b eamlet w eigh ts ha v e b een prop osed (see Shepard et al . [ 69 ]), as w ell as constrain ts that limit the v ariabilit y in the dose receiv ed b y v o xels within targets (see Hamac her and K ufer [ 32 ]). Ho w ev er, these t yp es of constrain ts do not oer sucien t rexibilit y to inruence the shap e of the D VHs. A new class of partial v olume constrain ts . The problem of form ulating constrain ts on the D VH resem bles, to a large exten t, a problem that has receiv ed m uc h atten tion in the nancial engineering literature, in particular in risk managemen t. Constrain ts on the D VH in that con text are constrain ts on the cum ulativ e distribution function of (usually in v estmen t) risks, and are notoriously dicult to solv e. Ho w ev er, recen tly an alternativ e approac h that form ulates risk managemen t constrain ts not in terms of the probabilit y distribution function of risk, but rather in terms of so-called tail aver ages has b een dev elop ed. It has b een sho wn that this approac h leads to a purely linear system of constrain ts. W e will apply this metho dology , whic h has b een dev elop ed in Ro c k afellar and Ury asev [ 56 , 57 ] and Ury asev [ 77 ], to the FMO problem. Instead of traditional D VH constrain ts, w e prop ose to incorp orate, in addition to the hard lo w er and upp er b ounds on the dose receiv ed b y eac h of the structures, alternativ e partial v olume constrain ts of the follo wing form: (P 01 ) The aver age dose receiv ed b y the subset of a target of relativ e v olume 1 receiving the lowest amoun t of dose m ust b e at least equal to L . (P 02 ) The aver age dose receiv ed b y the subset of a structure of relativ e v olume 1 receiving the highest amoun t of dose ma y b e no more than U . T o compare the traditional D VH constrain ts and our D VH constrain ts, it is useful to considering the follo wing e quivalent (although less in tuitiv e) form ulation of the traditional partial v olume constrain ts: (P 1 ) The maximum dose receiv ed b y the subset of a target of relativ e v olume 1 receiving the lowest amoun t of dose m ust b e at least equal to L .
PAGE 37
26 (P 2 ) The minimum dose receiv ed b y the subset of a structure of relativ e v olume 1 receiving the highest amoun t of dose ma y b e no more than U . Note that the dierence b et w een the traditional D VH constrain ts and our new D VH constrain ts lies in the fact that w e constrain the aver age dose in a tail of a giv en size as opp osed to the maximum or minimum dose in that tail. Figure 2{5 illustrates the t w o t yp es of D VH constrain ts for a particular dose densit y function. (Note that the D VH is obtained from the dose densit y function b y computing the upp er tail in tegrals, i.e., if the dose densit y function is denoted b y h ( d ; x ) then H ( d ; x ) = R 1 d h ( ; x ) d . In the medical ph ysics literature, the dose densit y function is often referred to as the dieren tial D VH.) Figure 2{5: T raditional and new (lo w er) D VH constrain ts for a target. Our approac h has t w o ma jor adv an tages. F rom a treatmen t qualit y p oin t of view, it not only con trols the exten t to whic h a soft b ound is violated, but also the actual magnitude of the violation. F rom a computational p oin t of view, our approac h yields a linear programming problem, while (as men tioned ab o v e) the traditional partial v olume constrain ts w ould either yield a v ery large-scale mixedin teger programming or a global optimization problem. Since v ery large-scale linear
PAGE 38
27 programming problems can b e solv ed v ery rapidly , our approac h allo ws a user to in v estigate man y more high qualit y candidate treatmen t plans, obtained with dieren t soft and hard b ound constrain ts. T o form ulate partial v olume constrain ts of the t yp e describ ed ab o v e, it will b e con v enien t to dene a loss function , to allo w us to treat the lo w er and upp er b ound constrain ts using the same metho dology . F or a target, lo w radiation doses are undesirable; w e therefore dene a loss function that is the negativ e of the dose calculation function L j s ( x ) = D j s ( x ) j = 1 ; : : : ; v s ; s = 1 ; : : : ; T F or all structures, v ery high dose lev els are undesirable; the corresp onding loss function is giv en b y the dose calculation function itself: L j s ( x ) = D j s ( x ) j = 1 ; : : : ; v s ; s = 1 ; : : : ; S: In the follo wing, w e will in tro duce and describ e the concepts of V alue-at-Risk and Conditional V alue-at-Risk in terms of a general loss function L j s ( x ). The V alue-at-R isk (V aR) at lev el (for short, the -V aR) is the smallest loss v alue with the prop ert y that no more than 100(1 )% of structure s exp eriences a larger loss. More formally , the -V aR is dened as s ( x ) = sup f 2 R : H s ( ; x ) g (2.1) when the loss function L is used, and s ( x ) = inf f 2 R : H s ( ; x ) 1 g (2.2) when the loss function L is used. The Conditional V alue-at-R isk (CV aR) at lev el (for short, the -CV aR) is then the aver age of all losses that exceed the -V aR.
PAGE 39
28 Ro c k afellar and Ury asev [ 56 ] sho w that the -CV aR can b e computed as s ( x ) = min 2 R ( + 1 (1 ) v s v s X j =1 ( L j s ( x ) ) + ) : (2.3) If the inequalit y in the appropriate denition of the -V ar in equation ( 2.2 ) or ( 2.1 ) is binding, the minim um in equation ( 2.3 ) is attained b y the -V aR, so that w e ha v e s ( x ) = s ( x ) + 1 (1 ) v s v s X j =1 ( L j s ( x ) s ( x )) + : (2.4) Otherwise, equation ( 2.3 ) suggests that the -V aR should b e redened sligh tly to appropriately accoun t for the discrepancy in the size of the tail. Equation ( 2.4 ) can b e explained b y recalling that the -CV aR is equal to the mean v alue of all loss v alues exceeding the -V aR. This mean is determined b y rst determining the -V aR v alue for a particular (see Fig. 2{6 (a) for the upp er 80%-V aR of a t ypical target dose densit y). Next, w e determine the normalized tail distribution of the losses exceeding the -V aR as w ell as its mean (see Fig. 2{6 (b)). Finally , the mean of the tail distribution is translated bac k b y adding the -V aR to it, thereb y undoing the earlier oset (see Fig. 2{6 (c)). F or eac h structure, w e ma y no w dene a set of soft b ounds on losses (i.e., lo w and high tail a v erages). F or targets, these will b e in the form of lo w er b ounds L s for 2 A s , where A s is a (nite) subset of (0 ; 1), and upp er b ounds U s for 2 A s , where A s is a (nite) subset of (0 ; 1) ( s = 1 ; : : : ; T ). F or the organs at risk, these will b e in the form of only upp er b ounds U s for 2 A s , where A s is a (nite) subset of (0 ; 1) ( s = 1 ; : : : ; S ). In case a target has in v aded a critical structure, an y partial v olume constrain ts that w e wish to apply to that critical structure should b e applied to the en tire critical structure rather than only the part of the critical structure outside the target. W e implemen t this b y form ulating suc h CV aR constrain ts to only apply
PAGE 40
29 (a) (b) (c) Figure 2{6: Computation of (upp er) Conditional V alue-at-Risk.
PAGE 41
30 to the part of the critical structure outside the target, b y suitably (re)dening the fraction . 2.3.4 An Enhanced LP F orm ulation with Dose-V olume Histogram Constrain ts W e can no w summarize our partial v olume constrain ts in terms of the in tensit y v ariables c haracterizing a treatmen t plan x : s ( x ) L s 2 A s ; s = 1 ; : : : ; T (2.5) s ( x ) U s 2 A s ; s = 1 ; : : : ; S (2.6) where and indicate whether the loss function L or L is b eing used, resp ectiv ely . W e can no w use the reform ulation of constrain ts ( 2.5 )-( 2.6 ) as describ ed in Ro c kafellar and Ury asev [ 56 ] and Ury asev [ 77 ], whic h leads to the follo wing enhanced form ulation of the FMO problem: minimize S X s =1 1 v s v s X j =1 F s ( z j s ) sub ject to (FMO 2 ) N X i =1 D ij s x i = z j s j = 1 ; : : : ; v s ; s 2 S z j s L s j = 1 ; : : : ; v s ; s = 1 ; : : : ; T z j s U s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S s 1 (1 ) v s v s X j =1 s z j s + L s 2 A s ; s = 1 ; : : : ; T (2.7) s + 1 (1 ) v s v s X j =1 z j s s + U s 2 A s ; s = 1 ; : : : ; S (2.8) x i 0 i = 1 ; : : : ; N z j s free j = 1 ; : : : ; v s ; s = 1 ; : : : ; S s free 2 A s ; s = 1 ; : : : ; T s free 2 A s ; s = 1 ; : : : ; S:
PAGE 42
31 Finally , observ e that equations ( 2.7 ) and ( 2.8 ) can easily b e reform ulated as purely linear ones b y in tro ducing articial v ariables, leading to a linear programming form ulation of the FMO problem. Note that CV aR partial v olume constrain t applied to a structure s w ould, after reform ulation using articial v ariables, add v s + 1 v ariables and constrain ts to the linear programming mo del. With the reduced t ypical problem dimensions as men tioned in Section 2.3.1 , w e w ould, for example, add ab out 10,000 v ariables and constrain ts to the mo del when a CV aR constrain t is applied to the targets. Due to the linearit y of the problem, suc h problem sizes remain manageable. It is in teresting to note that, when is c hosen to b e equal to 0, the corresp onding lo w er and upp er -CV aR v alues in fact coincide, and are equal to the mean dose in the en tire structure. The mean dose receiv ed b y all v o xels in a structure is sometimes also considered to b e an imp ortan t measure of the adv erse eect of radiation therap y on certain critical structures, suc h as saliv ary glands. The CV aR constrain ts can then b e simplied to so-called me an dose -constrain ts 1 v s v s X j =1 D j s ( x ) M s s = 1 ; : : : ; T 1 v s v s X j =1 D j s ( x ) M s s = 1 ; : : : ; S: 2.4 Results T o test our mo dels, w e ha v e studied cases of head-and-nec k cancer patien ts where the primary reason for using IMR T is to preserv e saliv ary function while obtaining adequate target co v erage. Three-dimensional treatmen t planning data for ten head-and-nec k cancer patien ts w ere exp orted from a commercial patien t imaging and anatom y segmen tation system (V o xelQ, Philips Medical Systems) at the Univ ersit y of Florida Departmen t Of Radiation Oncology . These data w ere then imp orted in to the Univ ersit y of Florida Optimized Radiation Therap y
PAGE 43
32 (UF OR T) treatmen t planning system and used to generate the data required b y the mo del describ ed in Section 2.3 . The UF OR T System generated a v o xel grid with a v o xel spacing of 4 mm, resulting in b et w een 58,665 and 195,506 v o xels. T o reduce the size of the optimization mo dels, w e increased the v o xel spacing to 8 mm for unsp ecied tissue. This resulted in an appro ximately 2-3 fold reduction in computation time. Ho w ev er, the plans w ere ev aluated using the original ner v o xel grid. T able 2{1 sho ws the dimensions and solution times for the ten cases. These results w ere obtained on a 2.8 GHz P en tium 4 laptop computer with 1 GB of RAM, and using CPLEX 8 as our linear programming solv er. T able 2{1: Mo del size and run times. case run time (sec.) # b eamlets # v o xels 1 48 1,017 17,471 2 52 1,195 19,543 3 104 1,745 36,136 4 164 1,804 38,724 5 52 1,182 17,466 6 32 1,015 14,147 7 56 1,178 22,300 8 62 1,090 19,417 9 29 921 15,416 10 155 2,073 30,664 As our planning criteria, w e ha v e adopted the criteria published in t w o studies of the Radiation Therap y Oncology Group [ 65 , 66 ] with minor mo dications. Summarizing the most imp ortan t criteria, the goal w as to treat t w o targets: Planning T ar get V olume 1 (PTV1) and Planning T ar get V olume 2 (PTV2) , whic h con tains PTV1. PTV1 consists of actual iden tied disease, while PTV2 is an expansion of PTV1 to accoun t for microscopic spread of disease as w ell as a margin to accoun t for the p ossibilit y of patien t and target mo v emen t and patien t setup uncertain t y . PTV1 should b e treated to a prescription dose lev el of D PTV1 R x = 73 : 8 Gy , and PTV2 should b e treated to a prescription dose lev el of D PTV2 R x = 54 Gy . In
PAGE 44
33 particular, this means that at least 95% of eac h target should receiv e the prescrib ed dose, no more than 1% of eac h target should b e underdosed b y more than 7%, and no more than 20% of the high-dose target should b e o v erdosed b y more than 10%. F or the four saliv ary glands (LPG: left parotid gland, RPG: righ t parotid gland, LSG: left submandibular gland, and RSG: righ t submandibular gland), the criterion for sparing is that either no more than 50% of the v o xels in the gland receiv e a dose exceeding 30 Gy , or the mean dose to the gland is no more than 26 Gy . No v o xels in the spinal cord (SC) should receiv e a dose exceeding 45 Gy , and no v o xels in the brainstem (BS) should receiv e a dose exceeding 54 Gy . T o prev en t undesirable lev els of o v erdosing in unsp ecied tissue to a v oid side eects of treatmen t, less than 1% of the corresp onding v o xels should receiv e more than 65 Gy . F or all ten cases, w e designed plans using sev en equispaced b eams at nominal settings, i.e., their orien tation w as not optimized for the individual patien ts. W e tuned the parameters of our mo del b y man ual adjustmen t using three of the patien t cases, resulting in the parameter v alues sho wn in T ables 2{2 and 2{3 . T o ensure feasibilit y , w e also included linear p enalizations of violations of the hard as w ell as the CV aR b ounds with v ery steep slop es (in the range 10 8 -10 12 ). Similar to Tsien et al . [ 76 ]), w e found that high p o w ers of dose dierence lead to excellen t results. T able 2{2: V alues of the co ecien ts of the v o xel-based p enalt y functions. s T C s = T H s L s C s p Cs U s H s p Hs PTV1 D PTV1 R x + 2 : 5 D PTV1 R x 0 : 5 7,500 12 D PTV1 R x + 5 : 5 7,500 6 PTV2 D PTV2 R x + 2 0 : 93 D PTV2 R x 75,000 12 D PTV1 R x + 5 : 5 75,000 6 Tissue 35 0 0 { D PTV1 R x + 5 : 5 8,000,000 5 LPG 0 0 0 { D PTV1 R x + 5 : 5 350,000 4 RPG 0 0 0 { D PTV1 R x + 5 : 5 350,000 4 LSG 0 0 0 { D PTV1 R x + 5 : 5 150,000 4 RSG 0 0 0 { D PTV1 R x + 5 : 5 150,000 4 SC 28 0 0 { 40 150 3 BS 33.6 0 0 { 48 200 3
PAGE 45
34 T able 2{3: V alues of the co ecien ts for the CV aR constrain ts. L ower CV aR-c onstr aints Upp er CV aR-c onstr aints s L s U s PTV1 0.90 0 : 97 D PTV1 R x 0.99 1 : 07 D PTV1 R x PTV2 0.90 0 : 96 D PTV2 R x 0.95 1 : 1 D PTV2 R x LPG { { 0.60 26 RPG { { 0.60 26 LSG { { 0.60 26 RSG { { 0.60 26 T o test robustness of the parameter v alues th us obtained, the tuned mo del w as then applied to all ten cases. W e c hose to appro ximate the piecewise p olynomial p enalt y function b y t w o PWL segmen ts for underdosing, four segmen ts for target o v erdosing, and t w o segmen ts for critical structure doses. The relev an t prop erties of the treatmen t plans are sho wn in T able 2{4 b elo w. Note that the criteria for all omitted structures (including spinal cord and brainstem) w ere satised in all cases. T able 2{4: UF OR T automated treatmen t-planning results for ten head-and-nec k cases. Case Structur e Criterion 1 2 3 4 5 6 7 8 9 10 All pass/fail pass pass pass pass pass pass pass fail pass pass PTV1 % D PTV1 R x 95 95 95 95 95 95 95 95 95 95 PTV1 % 0 : 93 D PTV1 R x 100 100 100 100 100 100 100 100 100 100 PTV1 % 1 : 1 D PTV1 R x 5 6 0 1 5 8 3 3 1 0 PTV2 % D PTV2 R x 96 97 98 98 98 96 96 96 96 96 PTV2 % 0 : 93 D PTV2 R x 99 100 100 99 99 99 99 99 99 100 Tissue dose @ 1% 60 57 55 57 60 59 56 57 58 58 LPG % 30 Gy 19 19 26 13 29 20 7 27 17 7 mean dose 19 19 17 14 25 21 11 24 18 14 RPG % 30 Gy 22 { 25 18 1 26 35 22 26 26 mean dose 19 { 17 17 11 23 28 22 26 26 LSG % 30 Gy 28 14 { 20 30 { 18 { 30 16 mean dose 26 20 { 24 25 { 20 { 21 19 RSG % 30 Gy { { { 35 25 { { 78 { { mean dose { { { 31 26 { { 41 { { SC % 45 Gy 0 0 0 0 0 0 0 0 0 0 BS % 54 Gy 0 0 0 0 0 0 0 0 0 0
PAGE 46
35 As can b e seen in T able 2{4 , nine out of ten cases passed all of the treatmen t planning criteria men tioned ab o v e. The only failure o ccurred in case 8, where the righ t submandibular gland could not b e spared. (Ho w ev er, note that in this case 61% of this gland is actually lo cated inside the target, and should th us receiv e a dose far exceeding the 30 Gy partial v olume constrain t.) T o ensure 95% co v erage of PTV1 at D PTV1 R x = 73 : 8 Gy , w e normalized all plans b efore ev aluating the other criteria. These results demonstrate that our mo del has the abilit y to pro duce high qualit y results that satised all or most of the published clinical criteria in a limited amoun t of time. In Fig. 2{7 , w e illustrate the solution to case 5 b y sho wing cum ulativ e D VHs for targets, saliv ary glands, spinal cord, and unsp ecied tissue for the solution found, as w ell as an example of iso dose curv e for an axial CT slice. These ndings ha v e b een published in Romeijn et al . [ 59 ]. 2.5 Extensions In this section, w e will discuss t w o imp ortan t issues in the determination of a go o d set of ruence maps (one for eac h b eam) that ha v e not b een addressed to date, but can b e incorp orated in our mo del while retaining linearit y of the problem. 2.5.1 Restoring F ractionation The v ast ma jorit y of IMR T treatmen t plans are not deliv ered to a patien t in a single session, but rather as a sequence of daily treatmen ts (often called fr actions ) o v er an extended p erio d of time (usually on the order of 4-8 w eeks), to tak e adv an tage of the fact that health y cells reco v er faster than cancer cells. The ma jorit y of con v en tional conformal external-b eam radiotherap y has b een dev elop ed using a narro w range of dose deliv ered p er fraction (see, e.g., Cheng and Das [ 27 ] and Bro c kstein and Masters [ 17 ]). This range has t ypically b een limited to 1.8-2 Gy p er daily fraction. Instead of a single ruence map set, w e therefore in fact need
PAGE 47
36 (a) (b) Figure 2{7: (a) Cum ulativ e D VHs from the solution for case 5. (b) Ov erla y of axial CT iso dose curv es (45, 54, 73.8 Gy) and structures (PTV1, PTV2, Spinal cord, Righ t parotid gland) for the solution for case 5.
PAGE 48
37 to determine m ultiple ruence map sets for eac h patien t. W e will next discuss the problems with the curren t practice of this issue. In t ypical clinical cases, a distinction is made b et w een t w o or more targets. F or example, w e ma y distinguish b et w een t w o Planning T arget V olumes, sa y PTV1 and PTV2. If the gross disease is lo cated in PTV1, this target will require a higher dose than PTV2. F or instance, adopting the dose lev els from the previous section, PTV1 ma y require 73.8 Gy , while PTV2 requires only 54 Gy . With con v en tional conformal radiation therap y , treatmen ts w ere t ypically carried out using sev eral dieren t conformal plans that w ere man ually designed for eac h target dose lev el. F or example, one conformal plan w ould then b e used to treat PTV2 to 54 Gy in 30 daily fractions of 1.8 Gy , while another conformal plan w ould treat PTV1 to an additional 19.8 Gy using 11 daily fractions of 1.8 Gy , ensuring a single dose p er fraction for all targets. With the adv en t of IMR T, it has b ecome a practical imp ossibilit y to repro duce this constan t dose-p er-fraction deliv ery tec hnique. This is due to the fact that it has b ecome a standard of practice in commercial planning systems to only solv e a single FMO problem. The single treatmen t plan is then simply divided in to a n um b er of daily fractions. Ho w ev er, if a single plan is optimized to deliv er dose to m ultiple target-dose lev els, then the dose p er fraction deliv ered to eac h target m ust c hange in the ratio of a giv en dose lev el to the maxim um dose lev el. There are essen tially t w o approac hes that are used in practice to deal with this issue. The rst is to set the highest target-dose lev el to the traditional dose p er fraction and allo w lo w er-dose targets to b e treated at lo w er doses p er fraction. In the example ab o v e, this w ould mean that 41 daily fractions are emplo y ed, eac h deliv ering 1.8 Gy to PTV1 for a total dose of 73.8 Gy , as desired. Ho w ev er, suc h a sc heme w ould then yield a dose of 54 = 41 = 1 : 32 Gy p er daily fraction to PTV2. The concern with this approac h is that the lo w er-dose targets ma y not
PAGE 49
38 receiv e sucien t cell kill, so the total prescription doses are often increased for these targets to yield a higher dose p er fraction. The second approac h is to set the lo w est target-dose lev el to the traditional dose p er fraction and allo w higher-dose targets to b e treated at higher doses p er fraction. In the example ab o v e, this w ould mean that 30 daily fractions are emplo y ed, eac h deliv ering 1.8 Gy to PTV2 for a total dose of 54 Gy , as desired. Ho w ev er, suc h a sc heme w ould then yield a dose of 73 : 8 = 30 = 2 : 46 Gy p er daily fraction to PTV2. The concern with this approac h is that the higher-dose targets receiv e m uc h higher daily doses that could result in an increase in undesirable side eects. T o a v oid ha ving dieren t doses deliv ered p er fraction to dieren t target-dose lev els, m ultiple ruence map sets m ust b e pro duced. The only w a y to ac hiev e this using curren tly a v ailable commercial systems is to determine these sets indep enden tly , whic h is clearly sub optimal. It is also p ossibly dangerous and requires extra care in reviewing the cum ulativ e plan, as critical structure sparing is not ensured for the cum ulativ e dose deliv ered. In this section, w e extend our mo del to allo w for the sim ultaneous optimization of m ultiple ruence map sets, one for eac h target-dose lev el. This will allo w IMR T practitioners to repro duce traditional dose fractionation sc hemes and tak e adv an tage of the clinical exp erience that has b een dev elop ed with these. In order to optimize m ultiple ruence map sets (one for eac h target or prescription-dose lev el) sim ultaneously , w e dene decision v ariables represen ting the in tensit y of b eamlet i in ruence-map set f as x if , and the decision v ector of all b eamlet in tensities in ruence-map set f b y x f . W e then express the dose receiv ed b y eac h v o xel in ruence-map set f as a linear function of the b eamlet in tensities as follo ws, analogous to the expression for D j s ( x ): D f j s ( x f ) = N X i =1 D ij s x if j = 1 ; : : : ; v s ; s = 1 ; : : : ; T ; f = 1 ; : : : ; T
PAGE 50
39 where w e ha v e assumed, for notational simplicit y , that all targets ha v e a distinct prescription dose lev el. F or v o xels in critical structures, the relev an t dose is the cum ulativ e dose receiv ed in all ruence-map sets. Ho w ev er, the situation is more complicated for targets. F or notational con v enience, w e will assume in this section that the targets are ordered in increasing order of prescription dose. This means that in ruence map set f , w e are in terested in treating all targets with a prescription dose larger than or equal to T f . Ho w ev er, an underdosing of v o xels in these targets in earlier ruence-map sets cannot b e comp ensated for in the curren t ruence-map set, as this w ould allo w the optimization to c hange the dose deliv ered p er fraction. W e therefore dene an articial cum ulativ e dose receiv ed b y target v o xels in the rst f ruence-map sets as the sum of the prescription dose for the preceding ruence-map set and the dose receiv ed in ruence-map set f . Denoting the prescription dose for ruence map set f b y P f , w e dene ~ D f j s ( x f ) = P f 1 + D f j s ( x f ) j = 1 ; : : : ; v s ; s = 1 ; : : : ; T ; f = 1 ; : : : ; T : In the optimization problem, w e will denote the articial doses b y the decision v ariables ~ z f j s . In the ob jectiv e function, w e p enalize the articial cum ulativ e doses for target v o xels in eac h ruence-map set according to the appropriate target-dep enden t p enalt y function. Eac h ruence-map set will only see the target v o xels that are included in its dose lev el. T o this end, w e dene v f = P Ts = f v s to b e the n um b er of target v o xels that are treated in fraction f ( f = 1 ; : : : ; T ). In addition, the true cum ulativ e dose receiv ed b y target v o xels will b e p enalized as unsp ecied tissue for eac h ruence-map set that do es not see it (where w e ha v e assumed that unsp ecied issue is the highest-n um b ered structure, S ). Finally , w e p enalize the true cum ulativ e dose receiv ed b y v o xels in all critical structures as in our single
PAGE 51
40 FMO mo del. Grouping the p enalt y terms b y fraction, this leads to the follo wing ob jectiv e function T X f =1 ( 1 v 1 v f + v S " f 1 X s =1 v s X j =1 F S P s + T X ` = s +1 z ` j s ! + v s X j =1 F S T X ` =1 z ` j S !# + 1 v f T X s = f v s X j =1 F f ( ~ z f j s ) + S 1 X s = T +1 1 v s v s X j =1 F s T X ` =1 z ` j s !) : Grouping the p enalt y terms b y structure leads to the equiv alen t ob jectiv e function T X s =1 v s X j =1 ( s X f =1 1 v f F f ( ~ z f j s ) + T X f = s +1 1 v 1 v f + v S ! F S P s + T X ` = s +1 z ` j s !) + S 1 X s = T +1 v s X j =1 T X f =1 1 v s F s f X ` =1 z ` j s ! + T X f =1 1 v 1 v f + v S ! v s X j =1 F S T X ` =1 z ` j S ! : It is clear that w e ma y incorp orate full and partial v olume constrain ts corresp onding to the targets for eac h individual ruence-map set, as w ell as full and partial v olume constrain ts corresp onding to all structures for the aggregate treatmen t plan. Using the same reform ulation as men tioned b efore, the mo del allo wing for fractionation can again b e form ulated as a linear program when the p enalt y functions are c hosen to b e piecewise linear and con v ex. Figure 2{8 illustrates the similarit y of the cum ulativ e D VHs for case 3 obtained using our mo dels for a single ruence map set (solid lines) and for t w o ruence map sets (dashed lines). The run time of the single-FMO mo del w as 104 seconds, while the run time of the m ultiple-FMO mo del w as 215 seconds. 2.5.2 Incorp orating Spatial Eects: The 2-Dimensional D VH The D VH that is often used in practice to ev aluate and compare treatmen t plans seems to con tain a ma jor shortcoming: it ignores spatial eects. In con v entional conformal radiation treatmen t planning, regions of (in particular) the target structures that receiv e less than some pre-sp ecied target or prescription dose are t ypically found on the p eriphery of the structure. Since the target usually includes
PAGE 52
41 Figure 2{8: Cum ulativ e D VHs from the solutions obtained using our mo del with and without accoun ting for fractionation. some safet y margin around the actual tumor, the consequences of these regions that are underdosed are usually minor. Ho w ev er, in IMR T treatmen t planning, it is v ery p ossible that underdosed regions actually app ear in the core of the target, whic h can ha v e p oten tially v ery serious eects. Since the D VH only measures the fr action of the structure that receiv es at least a certain dose, it is inheren tly incapable of distinguishing b et w een a treatmen t plan with underdosing near the b oundary of the target and a treatmen t plan with underdosing near the core of the target, when b oth underdosed regions are similar in size and amoun t of dose receiv ed. This problem is further amplied b y the fact that w e are dealing with 3-dimensional structures. One p ossible approac h ma y to dealing with this problem is to assign weights to v o xels, and mo dify the p enalt y function used in the mo del accordingly: S X s =1 P v s j =1 U j s F U s ( z ) + O j s F O s ( z ) P v s j =1 U j s + O j s where U j s and O j s are w eigh ting factors for v o xel j in structure s asso ciated with underdosing and o v erdosing, resp ectiv ely . F or example, w e ma y assign a higher w eigh t to v o xels that are deep inside a target (i.e., far a w a y from the b oundary of
PAGE 53
42 the target). F or critical structures, w e ma y assign a higher w eigh t to v o xels that are farther a w a y from the targets. T o demonstrate the application of spatial information to critical structures, w e applied a spatially dep enden t imp ortance to the unsp ecied tissue that increased linearly with distance from the surface of the union of all targets for case 6. Case 6 w as found to ha v e unacceptable o v erdoses at large distances from the targets and near the brainstem (see Fig. 2{9 ). W e applied the spatial p enalt y to attempt to solv e the former problem. W e found that o v erdosed regions in the plan could b e signican tly reduced at the cost of lo w ering the target co v erage for the lo w-dose target from 96% to 92% v olume at D R x . The dramatic impro v emen t in dose to tissue, with reductions of up to 20 Gy , is sho wn in Fig. 2{9 . The nal trade-o b et w een the t w o candidate plans should b e made b y a ph ysician. Fig. 2{10 sho ws the D VHs corresp onding to the original plan and the plan obtained using spatial correlations. Figure 2{9: An axial iso dose distribution from the original case 6, the spatially corrected case 6, and the dierence b et w een the t w o distributions. An alternativ e w ould b e to add an extra dimension to the D VH. Dempsey et al . [ 31 ] and Chao et al . [ 22 ] in tro duced a 2-dimensional generalization of the D VH that is a con v enien t to ol for recognizing the presence of regions of o v er or
PAGE 54
43 Figure 2{10: An o v erla y of D VH for case 6 with and without applying spatial p enalties. underdosing in unacceptable lo cations. The generalized D VH, whic h w e call a 2-dimensional dose-v olume-distance histogram (2D D VDH) has b oth dose as w ell as a distance measure as its argumen ts. The distance measure could, for instance, b e c hosen as the distance to the b oundary of a target, where outside the target the actual distance is used, and inside the target the negativ e of the actual distance is used. The D VH w ould then b e a function H : R + R ! [0 ; 1] where H ( d; r ) denotes the v olume of the part of a structure that receiv es d units of dose and has a distance measure of at most r , relativ e to the v olume of the part of the structure that has a distance measure of at most r . Denoting the set of v o xels in structure s that ha v e a distance measure of at most r b y V s ( r ), w e ha v e H s ( d; r ; x ) = jf j 2 V s ( r ) : D j s ( x ) d gj j V s ( r ) j :
PAGE 55
44 This clearly generalizes the usual concept of a D VH, since lim r !1 H s ( d; r ; x ) = H s ( d ; x ) : Note that the 2-dimensional D VH essen tially creates a D VH for eac h set V s ( r ). It is clear that, in a similar manner as for the original D VH, w e ma y incorp orate CV aR-constrain ts on the 2-dimensional D VH, thereb y more closely con trolling the dose distribution in eac h structure. F urther dev elopmen t of the concept and application of the 2D D VH is a topic for future researc h. 2.6 Concluding Remarks The linear programming mo dels that w e dev elop ed ha v e b een applied to realsized clinical problems and sho w great promise to adv ance the state-of-the-art of IMR T radiation-therap y treatmen t plan optimization. The mo dels are rexible and extendable, allo wing for the incorp oration of m ultiple ruence-map optimization to pro vide traditional dose fractionation and the incorp oration of spatial correlations in to the ob jectiv e function to allo w for spatial con trol o v er treatmen t plan features. Restoring traditional dose fractionation is extremely imp ortan t since a v ast b o dy of kno wledge and clinical exp erience exists with prescription dose lev els based on constan t daily dose fractions deliv ered to targets, while only limited kno wledge exists ab out the eects of treating dieren t targets using dieren t doses p er fraction. Incorp orating spatial correlations pro vides us with the p ossibilit y of con trolling the lo cation of o v erdosing and underdosing. F or example, o v erdosing tissue or a saliv a gland v ery near to a target, or underdosing a target near a spinal cord ma y b e acceptable, while underdosing cen tral regions in a target or o v erdosing critical structures a w a y from targets is not acceptable. Mo dels that only emplo y structure-dep enden t v o xel-based p enalt y functions and D VH constrain ts are blind to these dierences and will not b e able to distinguish b et w een certain clinically acceptable and unacceptable plans in a satisfactory manner.
PAGE 56
45 Our clinical exp erience sho ws that commercial treatmen t planning systems use heuristic metho ds, whic h o ccasionally terminate or get trapp ed in lo cal minima, corresp onding to inferior treatmen t plans. F urthermore, in curren t practice, generating a single IMR T treatmen t plan can require up to 20 min utes of computation time, and man y trials ma y b e required to pro duce a clinically acceptable treatmen t plan. This situation frustrates attempts b y treatmen t planners to ac hiev e go o d treatmen t plans b y logically adjusting the mo del input parameters, ev en if a go o d solution exists. F urthermore, these heuristics often ha v e no means to determine within a reasonable amoun t of computational time whether or not a giv en solution is indeed optimal. Our approac h has sho wn to consisten tly pro vide go o d treatmen t plans with a single set of problem parameters. In cases where trade-os m ust b e made, the computational eciency of our approac h mak es it feasible to quic kly generate alternativ e plans b y adjusting some of the problem parameters, emphasizing or de-emphasizing certain conricting goals based on the initial plan. Finally the com bination of con trol o v er solutions (through D VH constrain ts and spatial correlations) and the eciency of our mo del allo ws the treatmen t planner to ev aluate man y more plans than curren tly feasible, allo wing for the design of sup erior treatmen t plans for patien ts.
PAGE 57
CHAPTER 3 BEAM ANGLE OPTIMIZA TION 3.1 In tro duction When solving the IMR T treatmen t planning problem, the optimizations or decisions to b e made are t ypically divided as suc h: (i) determine the n um b er and orien tations at whic h radiation b eams will b e deliv ered; (ii) determine the ruence map for eac h radiation b eam selected in (i); and (iii) determine a metho d of decomp osition of the ruence map pro duced in (ii) to pro duce a deliv erable treatmen t plan. While subproblems (ii) and (iii) ha v e t ypically required optimization for a high qualit y result, b eam orien tations from subproblem (i) are often selected ad-ho c in concordance with previous conformal therap y practices. In all commercially a v ailable systems for treatmen t planning, these subproblems are considered separately as w ell. There is, ho w ev er, a gro wing realization that m uc h b etter solutions ma y b e obtained b y using a more formal approac h of in tegrating these optimization problems. Recen tly , the in tegration of pairs of these subproblems is b eing in v estigated, at least in some limited w a y . Examples include in tegration of (i) and (ii): Cro oks et al . [ 26 ], Das et al . [ 29 ], Pugac hev and Lei [ 51 ], Pugac hev and Xing [ 53 ], Pugac hev et al . [ 52 ], Pugac hev and Xing [ 54 ], Pugac hev, Bo y er et al . [ 50 ], W ebb [ 80 ]; (ii) and (iii) Bednarz et al . [ 9 ], Cho and Marks [ 25 ], Seco et al . [ 67 ], Shepard et al . [ 68 ], Sieb ers et al . [ 71 ]; and (i) and (iii) Otto and Clark [ 49 ] in to single mo dels. These in tegrated approac hes can pro vide the most desirable outcome as they allo w one to p erform a single optimization that considers all of the asp ects of b oth subproblems sim ultaneously . When subproblems are considered indep enden tly and sequen tially , the solutions are often not optimal with resp ect to all of the goals and constrain ts of b oth problems. This is particularly 46
PAGE 58
47 imp ortan t in ligh t of the recen t progress made b y Alb er et al . [ 4 ] in elucidating the degeneracy of the optimized solutions of the FMO problem (assuming of course that one has the optimal solution). In this w ork, it is noted that FMO problem solutions are often highly degenerate, p ossibly allo wing for subproblem in tegration with sub-problem (i) without signican t deterioration of the solution. Ho w ev er, the task of in tegrating requires v ery ecien t algorithms for solving the subproblems to com bat the increased dimensionalit y of the problem. Figure 3{1: Demonstration of degrees of freedom of an accelerator. [ 86 ] Beam-orien tation optimization (BOO) presen ts a signican t mathematical c hallenge as the in tro duction of b eam-orien tation degrees of freedom (see Fig. 3{1 for demonstration of rotational and translational degrees of freedom of a commercial accelerator) in to a FMO problem results in either a noncon v ex ob jectiv e function, ev en if the ob jectiv e function w as con v ex for the original FMO problem (Bortfeld et al . [ 13 ]) or a problem with disjoin t feasible solutions spaces (Lee et al . [ 41 ]). In either case, this leads to the existence of lo cal minima and a hard problem to solv e. Th us, in tegration of the BOO and FMO problems ha v e made use of heuristics based on con v en tional conformal radiation therap y . The approac hes that ha v e b een prop osed to date for optimizing the b eam angles can b e categorized
PAGE 59
48 in to t w o broad classes. In the rst approac h (Das and Marks [ 28 ], Pugac hev and Lei [ 51 ], Ro wb ottom et al . [ 64 ]), eac h candidate b eam angle is giv en a n umerical rating based on its eect on the targets and critical structures. Then, a presp ecied n um b er of b eams, sa y p , with the b est n umerical ratings are selected. The second approac h, studied b y Pugac hev and Xing [ 54 ], and Pugac hev et al . [ 50 ], is a lo cal searc h based approac h. Initially , p b eam p ositions are selected, and corresp onding b eamlet in tensities are determined. Then, one or more b eam p ositions are c hanged, new b eam or b eamlet in tensities are determined, and the impact of this c hange is assessed. Generally , if the c hange impro v es the treatmen t plan, the c hange is made; otherwise another c hange is considered. This pro cess is then rep eated iterativ ely . V arian ts of this approac h as part of a sim ulated annealing or genetic algorithm framew ork ha v e also b een implemen ted and studied (Langer et al . [ 38 ], W u and Mohan [ 83 ], W u and Zh u [ 84 ]). While man y of these studies ha v e found that optimizing b eam angles can signican tly impro v e the resulting treatmen t plan, the ab o v e approac hes are somewhat limited as it is w ell kno wn that the true inruences of a set of gan try angles on the nal dose distribution is not kno wn un til optimization of the b eam in tensit y prole is p erformed. Moreo v er, o wing to the heuristic nature of the algorithms, the curren t metho dologies cannot pro vide a fair comparison b et w een plans with dieren t n um b ers of b eams. In a t ypical commercial a v ailable IMR T treatmen t planning system, the b eam orien tations are selected in an ad-ho c manner, and then optimization is done to determine b eam ruences. Dev eloping an algorithmic approac h that in tegrates b eam orien tation optimization (BOO) with the ruence map optimization (FMO) and sim ultaneously determines b eam orien tations and b eam ruences will b e of substan tial v alue to clinical IMR T practices. In this w ork, w e describ e a heuristic approac h to in tegrate BOO with FMO. Our approac h starts with a set of b eam orien tations and rep eatedly impro v es it b y adding and dropping new b eams at
PAGE 60
49 dieren t orien tations. A no v el feature of our approac h is that w e determine the optimal b eam in tensities b y solving a linear programming problem and use the sensitivit y analysis information of this problem to iden tify b eam orien tations to b e added. Our computational results indicate that the metho ds dev elop ed can appro ximate the b eha vior of the solution space v ery precisely and th us can guide the searc h in correct direction. This c hapter is organized as follo ws. In Section 3.2 , w e discuss the metho dology used for form ulating the BOO problem. W e also presen t t w o exact approac hes for BOO problem in this section. In Section 3.3 , w e discuss neigh b orho o d searc h based heuristics based on sensitivit y and greedy searc h based routines. In Section 3.4 , w e presen t the results and nally conclude the ndings in Section 3.5 . 3.2 Metho dology Our algorithm to in tegrate BOO with FMO relies on the linear programming form ulation FMO 1 of the FMO problem giv en in Section 2.3 . In form ulation FMO 1 , hard upp er and lo w er b ounds are imp osed on the dose receiv ed b y a v o xel. In this c hapter, w e will discuss the form ulation with the implicit assumption that these b ounds are enforced as soft b ounds and th us, are part of the ob jectiv e function no w. The form ulation FMO 1 has b een further extended in Section 2.3 to incorp orate the dose-v olume histogram (D VH) constrain ts. F or simplicit y , in this c hapter, w e ha v e not dealt with these constrain ts. The concepts discussed in this c hapter can b e easily extended to accommo date these set of constrain ts without m uc h impact on the complexit y of the problem. 3.2.1 A Mixed In teger Linear Programming F orm ulation for the BOO Problem W e will no w use our form ulation FMO 1 to dev elop a mixed in teger linear programming (MILP) form ulation of the BOO problem. In this form ulation, w e
PAGE 61
50 dene an in teger v ariable y l for eac h l 2 B as follo ws: y l = 8><>: 1, if beam at or ientation l is used 0 ; other w ise W e refer to y l 's as b eam v ariables. In terms of these b eam v ariables and the b eamlet in tensit y v ariables ( x i 's), the BOO can b e form ulated as: minimize S X s =1 1 v s v s X j =1 F j s ( z j s ) (3.1) sub ject to X i 2 N D ij s x i = z j s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S (3.2) X i : i 2 N l x i M y l l 2 B (3.3) X l 2 B y l p (3.4) y l = f 0 ; 1 g l 2 B (3.5) x i 0 i 2 N (3.6) (3.7) and M is a sucien tly large n um b er. In this form ulation, the constrain ts ( 3.2 ) represen t the dose receiv ed b y the v o xels in terms of the b eamlet in tensities. The constrain ts ( 3.3 ) ensure that b eamlets coming from a b eam cannot tak e p ositiv e in tensit y v alues unless the corresp onding b eam is turned on. Finally , the constrain t ( 3.4 ) states that no more than p b eams are turned on. Let us no w in v estigate the size of this MIP form ulation. Let us consider 72 p ossible b eam orien tations discretized at 5 in terv als (suc h as, 0 , 5 , 10 , 15 ,. . . , 355 , whic h giv es 72 in teger v ariables y l 's denoting whic h b eams will b e turned on. F or eac h b eam v ariable y l , w e ha v e around 150 b eamlet v ariable, x i giving
PAGE 62
51 a total of 11,800 b eamlet v ariables. In addition, w e ha v e 100,000 v o xel v ariables z j s . The form ulation consists of 40,000 constrain ts. W e attempted to solv e this problem using CPLEX 8.0, the state-of-the-art commercial solv er whic h could not solv e it optimally or near-optimalit y in 3 da ys on a P en tium IV 2.8 Ghz computer. It app ears that MIP of this t yp e and size are v ery dicult to solv e b y general purp ose co des in acceptable running times. Th us, w e need sp ecial-purp ose heuristic metho ds that can solv e BOO eectiv ely and ecien tly . W e next describ e suc h an approac h. 3.2.2 Searc hing the Solution Space for Equispaced Beams The feasible solution space of the p ossible b eam orien tations is v ery large. Supp ose that w e allo w candidate b eam orien tations to b e m ultiples of 5 . Then, w e ha v e 360 /5 = 72 c hoices of b eam orien tations. Supp ose further that w e are in terested in a v e b eam plan. Then, there are 72 C 5 = 13,991,544 p ossible com binations of b eam orien tations. In principle, w e can solv e the BOO problem b y considering eac h p ossible set of b eam orien tations and solving a FMO problem for eac h of these orien tations. The set of b eam orien tations whic h yields the minim um ob jectiv e function v alue of the FMO problem giv es the global optimal solution of the BOO problem. Ho w ev er, this approac h en tails solving to o man y FMO problems and is to o time-consuming to b e of m uc h practical use. One p ossible approac h to limit the size of the searc h space of p ossible b eam orien tations is to consider only equi-spaced b eam orien tations. Though highly constrained, equi-spaced b eams sets ha v e b een demonstrated to yield high qualit y plans, as b eams cannot b e to o closely spaced without incurring signican t hot sp ots at shallo w b eamlet in tersections in the patien t. F or the equi-spaced b eam orien tations case, w e solv e the BOO problem in the follo wing manner. (W e illustrate our metho d for the v e b eam case and translation to dieren t n um b er of b eams is straigh tforw ard.) W e consider the equi-spaced b eams at angles 0 ,
PAGE 63
52 72 , 144 , 216 , and 288 and solv e a FMO problem to determine optimal b eamlet in tensities; let z 0 denote the ob jectiv e function v alue of this FMO problem. W e next rotate all the b eams b y 1 to obtain a new set of b eam orien tations 1 , 73 , 145 , 217 , and 289 , and solv e a FMO problem to determine optimal b eamlet in tensities with z 1 as the corresp onding ob jectiv e function v alue. W e again rotate all the b eams b y 1 , and rep eat this pro cess. After this pro cess has b een rep eated 72 times, w e obtain the starting set of b eam orien tations and the pro cess is terminated. No w, w e compare the ob jectiv e function v alues z 0 , z 1 ,. . . , z 71 , and the set of b eam orien tations corresp onding to the lo w est ob jectiv e function v alue is the optimal solution of the BOO problem. The preceding metho d requires solving 72 FMO problems for the v e b eam case. In general, for the case of L b eams, it will require solving b 360 =L c FMO problems. Though this approac h is far more practical than en umerating all nonequispaced b eam orien tations, it is still a time-consuming approac h. W e can sp eed up this pro cess further b y rotating the b eam orien tations b y a n um b er larger than 1 . F or example, if w e rotate the b eam orien tations b y 5 , then the b eam orien tations considered b y this approac h will b e f 0 , 72 , 144 , 216 , 288 g , f 5 , 77 , 149 , 221 , 293 g , . . . , f 65 , 137 , 209 , 281 , 353 g , and w e need to solv e b 72 = 5 c =14 FMO problems. Our linear programming problem (1) for the FMO problem can t ypically b e solv ed within a min ute, and hence the BOO problem can b e solv e within 15 min utes whic h is acceptable from the clinical p oin t of view. W e presen t computational results of this approac h in Section 3.4 . 3.3 Neigh b orho o d Searc h Algorithm W e will no w discuss a more adv anced searc h metho d. The BOO problem is a discrete optimization problem, where w e need to select p b eams out of a discrete set of 72 b eams. Neigh b orho o d searc h algorithms are the most p opular algorithms
PAGE 64
53 to solv e discrete optimization problems. A neigh b orho o d searc h algorithm t ypically starts with an initial feasible solution of an optimization problem, and rep eatedly replaces it b y an impro v ed solution in its neigh b orho o d un til it cannot b e impro v ed an y more. Let us represen t an y feasible solution feasible solution of the BOO problem b y set ( x , y ) where v ector x giv es the in tensities of b eamlets in the b eams that are turned on and the b o olean v ector y indicates if the b eam at dieren t orien tations are on or o. Giv en a feasible solution ( x , y ) of the BOO problem, the neigh b orho o d searc h algorithm denes a set of neigh b oring solutions that are obtained b y p erturbing the solution y to another solution in its neigh b orho o d. Let us denote the neigh b orho o d b y y neig hbor . W e will describ e later in detail ho w w e dene the neigh b orho o d of a solution. W e ev aluate eac h neigh b or in y neig hbor and try to iden tify a neigh b or y 0 , whose ob jectiv e function v alue z ( x 0 , y 0 ) is less than curren t ob jectiv e function v alue z ( x , y ). If w e nd an impro v ed neigh b or y 0 , w e replace y b y y 0 , dene the new neigh b orho o d with resp ect to the new solution and rep eat this pro cess. If w e are unable to nd an impro v ed neigh b or of y , w e declare y to b e a lo cal optimal solution and the algorithm terminates. Belo w is the summary of the steps of our neigh b orho o d searc h algorithm for the BOO problem. Step 0. Obtain an initial feasible solution ( x , y ) of BOO; Step 1. If there exists a neigh b or y 0 2 y neig hbor with z ( x 0 , y 0 ) < z ( x , y ), go to step 2. Else, go to step 3. Step 2. Replace y b y y 0 . Go to step 1. Step 3. Output ( x , y ) as a lo cal optimal solution of BOO. In the follo wing sections, w e will describ e the steps of the neigh b orho o d searc h algorithm in detail. 3.3.1 Finding an Initial F easible Solution In searc h algorithms, the con v ergence sp eed and the qualit y of the nal solution usually dep end on the qualit y of initial solution. So, it is imp ortan t to nd
PAGE 65
54 a go o d qualit y solution to start the searc h. W e discuss t w o metho ds to nd a go o d initial solution for the neigh b orho o d searc h algorithm. An equispaced b eam plan. W e can assume that the initial feasible solution b eams are equi-spaced. Recall from Section 3.2.2 that if w e consider 5 discretization of b eam orien tations, then there are 14 distinct equi-spaced b eam orien tations sets. W e exp erimen ted with t w o p ossibilities: (i) when w e select an y arbitrary set of equi-spaced b eam orien tations; or (ii) when w e consider eac h distinct equi-spaced b eam orien tation, solv e a FMO for eac h orien tation, and select the b est b eam orien tation. W e found in our computational in v estigations that though the second p ossibilit y ga v e the b etter initial solutions, the nal solution obtained when the neigh b orho o d searc h algorithm terminated, had comparable ob jectiv e function v alues. W e th us selected an arbitrary set of equi-spaced b eam orien tations as the initial solution since it can b e found without an y computational eort. A gradien t based approac h for generating go o d initial solution . In this subsection, w e will dev elop a metho dology to decide the go o dness of a b eam to b e selected using concepts similar to the sensitivit y analysis in linear programming. In this metho d, w e start with a 0-b eam plan and then b eams are added progressiv ely to the existing set. Eac h new b eam is c hosen in suc h a w a y that it aids the existing b eams in ac hieving b etter co v erage of tumor while sparing the health y organs. The metho d analyzes individual b eamlets and decides their usefulness b y their p oten tial to impro v e the ob jectiv e function. This information is then used to decide go o dness of the b eams. By replacing the dose expressions z j s in the ob jectiv e function in form ulation FMO 1 b y the dose calculation functions P i 2 N D ij s x i , the FMO problem can b e giv en as (Note that the hard upp er and lo w er b ounds imp osed on dose receiv ed b y a v o xel will b e enforced b y p enalizing these in the ob jectiv e function)
PAGE 66
55 minimiz e w ( x ) = S P s =1 v s P j =1 F j s P i 2 N D ij s x i sub ject to x i 0 i 2 N whic h is a con v ex optimization problem with just non-negativit y constrain ts. Note that in this section, the ob jectiv e functions are the con v ex dieren tiable functions men tioned in Section 2.3.1 . W e do not use the piecewise linear equiv alen t of the functions here. No w, let us dene some notations here. A t some bixel in tensit y solution x , dene gradien t +i = @ w ( x ) @ x i i = 8><>: @ w ( x ) @ x i , if x i > 0 0 if x i = 0 F or a small c hange in the in tensit y of bixel i , +i represen ts the impro v emen t in the ob jectiv e function p er unit increase in the in tensit y of bixel i . It is easy to observ e here that a bixel i with negativ e +i will impro v e the ob jectiv e function if its in tensit y is increased. It can also b e observ ed that lo w er (more negativ e) the v alue of +i is, the more promising the bixel is. Similarly , i represen ts the impro v emen t in the ob jectiv e function p er unit decrease in the in tensit y of bixel i . W e prop ose to decide the go o dness of a b eam based on the gradien t information of the bixels in the follo wing w a y: i) A t an y solution x , ev aluate the gradien t +i for all the bixel ( i 2 N b ) in b eam b . ii) Denote the set of useful bixels U b = f i : +i < 0 ; i 2 N b g . iii) Dene the go o dness of the b eam b as the sum of the gradien ts of the bixels in set U b , i.e., Go o dness(b) = P i 2 U b +i ! .
PAGE 67
56 Using the ste ep est gr adient approac h for deciding go o dness of b eams, w e can generate a go o d p -b eam plan in the follo wing w a y . W e ev aluate the gradien t information of the b eams and calculate the go o dness of the b eams based on this information. The b eam with maxim um go o dness is c hosen. Then, w e form ulate a FMO problem for b eams selected and solv e the LP form ulation to get the optimal solution. A t the optimal solution, the b eamlet gradien t information is calculated. Using this information, the go o dness of the b eams whic h are not curren tly activ e is ev aluated. The b eam with maxim um go o dness is c hosen to b e selected next. Rep eating the pro cess p times giv es us an initial p -b eam plan. The algorithm for this sc heme is giv en b elo w. Step 0. Set x 0i = 0 ; 8 i 2 N . Set y b = 0 8 b 2 B . Set t = 0. Step 1. A t solution x t , ev aluate the go o dness of the b eams not included in the v ector y . Find the b est b eam (sa y , b eam l ) from this information. Set y l = 1. Set t = t +1. Step 2. If t = k , stop. Else, form ulate the FMO problem for the b eams turned on in v ector y and solv e to get new FMO solution x t and go to step 1. In the results section, w e demonstrate that this sc heme do es b etter than sc heme 1. In most of the cases, the initial solutions are of b etter qualit y and th us result in lesser n um b er of iterations of neigh b orho o d searc h. 3.3.2 The Neigh b orho o d Searc h Algorithm Our searc h algorithms are drop-add algorithms. A t eac h step, w e drop one b eam and try to nd a b etter b eam in its vicinity . Sa y , b eam l (at angle ) w as dropp ed. W e will refer to the b eams whic h are not p erturb ed in this iteration as xe d b eams. W e dene the vicinit y of the b eam b y b eams with are no more than degrees apart from b eam l . The candidate b eams are selected at some discrete steps (sa y , 5 ) in the range ( ; + ). Then, for the iteration when b eam l is dropp ed and a candidate b eam is b eing searc hed for, the neigh b orho o d , y neig hbor is
PAGE 68
57 dened as, y neig hbor b = 1 ; if y b = 1 ; for all b 2 B ; b 6 = l y b =1, for a b eam c hosen in vicinit y of b eam l y b =0, otherwise Beams cannot b e to o closely spaced without incurring signican t hot sp ots at shallo w b eamlet in tersections in the patien t. So, the candidate b eams include only those b eams whic h are at least r (sa y 15 ) apart from the xe d b eams. Figure 3{2: The neigh b orho o d of a 4-b eam plan. Complete en umeration of the neigh b orho o d. W e ha v e dev elop ed v ery fast (t ypical run-time, less than a min ute) LP mo dels to solv e the FMO problems. W e can use these as the basic routines for the neigh b orho o d searc h. Sa y , for example, for ev ery set of b eams in y neig hbor , w e can form ulate an LP and solv e the FMO to nd the optimal ob jectiv e. After solving FMO for all the instances in the set, w e can compare their ob jectiv e v alues to nd the b est b eam set in the neigh b orho o d. Steep est descen t metho d . A complete en umeration helps us get exact information ab out the neigh b orho o d y neig hbor . It giv es the \optimal" solution resp ect to the ob jectiv e function dened for all the b eam sets in the neigh b orho o d. But, if w e ha v e to p erform large n um b er of iterations b efore terminating, w e will need to solv e
PAGE 69
58 lots of FMO problems as subproblems. Th us, this metho d migh t not b e acceptable b ecause of high run-time. In last section, w e sho w ed ho w w e can appro ximate the go o dness of a b eam resp ect to an y giv en solution x using the gradien t information. W e can use this metho d to searc h the neigh b orho o d. F or ev ery elemen t in y neig hbor , w e do the follo wing: Step 1 :F orm ulate and solv e FMO problem for the p -1 xe d b eams (i.e., excluding the b eam l ) and denote the solution b y x 0 . Then ev aluate +i with resp ect to the solution v ector x 0 for eac h bixel in the candidate b eam to b e added. Step 2 :Find the go o dness of the b eams using ste ep est desc ent metho d. Pic k the b eam with highest go o dness (sa y , l 0 ) as the b eam whic h will p oten tially replace the b eam l . The b ottlenec k op eration in the searc h metho d usually is to nd an impro v ed neigh b or. The ste ep est gr adient metho d is v ery fast. F or a neigh b or with 10{15 candidate b eams, the go o dness of the b eams can b e ev aluated in less than a second of CPU time. W e discuss the computational results for this metho d in Section 3.4 . Greedy algorithm . W e discuss an alternate metho d to nd the go o dness of the candidate b eam in this sub-section. W e form ulate and solv e the FMO problem for ( p -1) xe d b eams to get the optimal ruence map. A t this solution, w e p erform a greedy searc h with the in tensit y of the b eamlets in the candidate b eam as v ariables. The solution of lo cal searc h appro ximates the p oten tial of the candidate b eam to impro v e the plan. F or an y candidate b eam, the greedy algorithm starts with the FMO solution of ( p -1) b eam plan and increases in tensit y of the bixels in candidate b eam to impro v e the solution. Sa y , w e increase the in tensit y of bixels in small steps of size . Our algorithm uses the follo wing greedy principle to iden tify the b eamlet in tensit y to b e increased: for eac h b eamlet i , it considers c hanging its in tensit y from x i to x i + ,
PAGE 70
59 and mak es the c hange with the maxim um impro v emen t on the ob jectiv e function v alue FMO 1 . The algorithm main tains a term i for eac h b eamlet i , whic h denotes its impr ovement p otential. The i is computed with resp ect to a giv en b eamlet in tensit y v ector u in the follo wing manner: i = z ( x 1 ; x 2 ; : : : ; x i + ; : : : ; x n ) z ( x 1 ; x 2 ; : : : ; x i ; : : : ; x n ) A t an y iteration, the greedy algorithm selects a b eamlet i with the minim um v alue of i , and c hanges its in tensit y to x i + . Since impr ovement p otential i 's are functions of the b eamlet in tensities x , their v alues c hange as x c hanges. Hence, after up dating in tensit y of an y b eamlet, i 's need to b e up dated. The algorithm up dates i 's and rep eats the pro cess un til it do es not nd an y more impro v emen t. The b ottlenec k op eration in the greedy algorithm is the computation of i 's in ev ery iteration. In our algorithm, w e do not compute eac h i from scratc h whenev er a b eamlet in tensit y is increased, but up date only those i 's whic h c hange. Up dating i 's is substan tially faster than re-computing these and mak es the greedy heuristic v ery ecien t. T o p erform the up dates ecien tly , w e main tain the follo wing link ed lists: Set of v o xels aected b y b eamlet i , that is, V O X E LS ( i ) = f j 2 V : D ij > 0 g ; Also, for eac h v o xel j , w e store the set of b eamlet that aect v o xel j , that is, B E AM LE T S ( j ) = f i 2 N : D ij > 0 g . Consider an itreration where in tensit y of b eamlet k is c hanged. Then, only for the b eamlets in the set [ j 2 V O X E LS ( k ) B E AM LE T S ( j ), i v alues ha v e to b e up dated. All other i 's will remain unc hanged. This metho d giv es excellen t results with an acceptable run time of 10{15 min utes for our sample cases. The results are discussed in Section 3.4 .
PAGE 71
60 3.4 Computational Results and Analysis In head and nec k cases, b ecause of pro ximit y of critical organs to the tumors and the o ccasional o v erlap of these with the Planning T arget V olume (PTV), in tensit y mo dulation and b eam-orien tation pla y a critical role. By orien ting the b eams correctly , w e try to a v oid the b eams coming from the directions of the critical organs whic h are tough to spare. W e ha v e applied our approac hes to 3 cases of head-and-nec k cancer patien ts. The patien t data is generated as discussed in Section 2.4 . 3.4.1 Impro v emen t due to Equispaced BOO Metho dology W e observ ed that rotating the b eam orien tations b y 5 w as not noticeably w orse than rotating the b eam orien tations b y 1 in terms of the solution qualit y and, since it sp eeds up the running time b y a factor of 5. T o amplify the results, w e sho w the b est and w orst 4 b eam plans in Fig. 3{3 . (a) W orst Plan (b) Best Plan Figure 3{3: Results of equi-spaced BOO for Case 01 for a 4-b eam plan. In t w o separate cases studied, equispaced BOO allo w ed b etter critical structure sparing when sparing w as not ac hiev ed for the initial angle set. The spinal cord
PAGE 72
61 sparing is ac hiev ed in plan (b) while plan (a) fails to do so. The dose conformit y of the tumors is also impro v ed. 3.4.2 Results of the Neigh b orho o d Searc h Metho d W e compare the results of the b est plan ac hiev ed b y the equi-BOO metho d after applying our neigh b orho o d searc h sc heme on it. It can b e observ ed that the hotsp ots (the fraction of target receiving dose more than the prescription dose)in target PTV 1 are reduced and b etter sparing of the saliv a glands is ac hiev ed b y using the b eams c hosen b y the neigh b orho o d searc h metho d. Figure 3{4: Results of neigh b orho o d searc h for Case 01 for a 3-b eam plan. W e also dene a h ybrid sc heme whic h searc hes the neigh b orho o d in eac h iteration using b oth steep est gradien t sc heme and greedy searc h and pic ks the b etter c hoice from the options pro vided b y the t w o metho ds. In the table b elo w, w e ha v e summarized these results. In the table, w e compare the ob jectiv e function v alue of the plans. As is eviden t from the results in the table, the h ybrid sc heme
PAGE 73
62 p erforms b etter than b oth the individual sc hemes and the results are comparable to the en umeration metho d. T able 3{1: T able of results for the neigh b orho o d searc h algorithms. Initial Solution Final Solution Case No. ofBeams Metho d Used for InitialSolution InitialOb-jectiv e V alue ( 10 13 ) Metho ds and Ob jectiv e v alues ( 10 13 ) En umeration Steep est Gradien t GreedySearc h Hybrid 1 3 Equispaced 5.45 4.41 4.75 4.67 4.59 3 Steep est Gradien t 5.26 4.31 4.38 5.19 4.43 4 Equispaced 8.89 3.31 3.38 3.67 3.31 4 Steep est Gradien t 3.60 3.32 3.49 3.38 3.38 2 3 Equispaced 13.41 4.85 5.39 11.31 5.16 3 Steep est Gradien t 9.13 4.98 5.84 5.75 5.27 4 Equispaced 6.82 4.15 4.29 4.30 4.15 4 Steep est Gradien t 5.17 3.36 3.52 4.35 3.43 3 3 Equispaced 1.17 0.97 1.02 1.12 0.97 3 Steep est Gradien t 1.09 0.97 1.02 0.97 0.97 4 Equispaced 1.44 0.77 0.74 0.74 0.74 4 Steep est Gradien t 0.78 0.73 0.78 0.76 0.76 These results, though studied o v er a small sample of clinical cases, demonstrate that BOO optimization can signican tly help in impro ving plan qualit y . Lik e other researc hers, w e observ ed that the BOO optimization approac hes are useful only in plans with small n um b er of b eams. In a plan with large n um b er of
PAGE 74
63 b eams, a randomly equi-spaced b eam plan is almost as go o d as an optimized b eam set-up. Since in clinic, it is desirable to ha v e less n um b er of b eams, our study aids in ac hieving this goal. The equi-BOO metho dology dev elop ed in Section 3.2.2 pro vides excellen t plans. These metho d ha v e additional adv an tage of b eing pro v en globally optimal solution for the constrained space considered. Th us, it can serv e as a b enc hmark for other BOO optimization algorithms. In our study , w e ha v e used the ob jectiv e function minimization to b e our goal. Ho w ev er, as previously noted, no clear guidelines ha v e b een dev elop ed to decide the preference of a solution o v er other. Should w e need to adopt a dieren t metho dology (for example, comparing the Dose-V olume-Histograms based on some pre-dened criteria), the searc h based metho d dev elop ed in this c hapter can easily b e mo died to adopt this rubric easily . Our sensitivit y analysis based neigh b orho o d searc h are rst attempt of its kind for BOO optimization. Earlier metho ds ha v e b een based on either complex MIP form ulation or b eam-ey e-view heuristics. W e b eliev e that this approac h will pro vide a new direction to the ongoing researc h on this complex problem. W e ha v e only considered in tegration of BOO with FMO problem here. In Chapter 5, w e will consider in tegration of FMO with deliv ery . As a future researc h direction, in tegration of all the three subproblems together will b e an in teresting topic to in v estigate. 3.5 Conclusions In clinic, the usual practice is to use more b eams in a plan if it is not p ossible to generate go o d qualit y plan with less n um b er of b eams. Large n um b er of b eams, ho w ev er, results in higher v erication and treatmen t time and is also prone to deliv ering more normal tissue dose. Th us, it is desired to k eep the n um b er of b eams used to the minimal. Our results clearly demonstrate that b eam orien tation optimization can help in suc h cases. BOO is imp ortan t and can comp ensate for
PAGE 75
64 some observ ed deliv ery tec hnique dierences. In this c hapter, w e ha v e prop osed practical metho ds to address this problem. W e also prop osed a metho dology to b enc hmark heuristic approac hes for the b eam-angle-optimization problem.
PAGE 76
CHAPTER 4 DELIVER Y 4.1 In tro duction Conformal radiation therap y seeks to conform the geometric shap e of the deliv ered radiation dose as closely as p ossible to that of the in tended targets. In c onventional conformal radiation therap y , this is tried to ac hiev e b y deliv ering a radiation dose b y using a relativ ely small n um b er of dieren t b eam shap es, called ap ertur es (or shap e matrices), eac h with uniform in tensit y lev el. Suc h ap ertures are formed b y partially blo c king the radiation source using customized st yrofoam blo c ks. IMR T allo ws for the creation of v ery complex non-uniform dose distributions that allo w the deliv ery of sucien tly high radiation doses to targets while limiting the radiation dose deliv ered to health y tissues. Ho w ev er, IMR T also results in requiring the system to deliv er v ery complex dose proles. Because of v ery complex proles in IMR T, it is no longer p ossible to deliv er these doses using blo c ks as it w ould mak e the pro cess v ery inecien t. Computer-con trolled Multileaf Collimators (MLCs, see Fig. 4{1 ) are th us in tro duced to deliv er IMR T ecien tly . T o deliv er a desired dose pr ole (more appropriately , ruenc e map ), the dose proles are decomp osed in to ap ertures (or shap e matrices) with asso ciated in tensities. The MLC system is able to generate complex dose distributions b y dynamic al ly blo c king dieren t parts of the b eam to create dieren t ap ertures. The problem of con v erting the dose proles in to a set of le af se quenc e les that con trol the mo v emen t of the MLC during deliv ery is referred to as MLC se quencing pr oblem . 65
PAGE 77
66 Figure 4{1: A m ulti-leaf collimator system 4.1.1 Problem Description The b eam head at eac h gan try angle can b e conceiv ed of as a rectangle whic h is discretized in to an m n rectangular grid; eac h individual rectangle in this grid is called a bixel . A t eac h gan try angle, the radiation head deliv ers the radiation dose whic h is sp ecied b y an p q matrix I , called the ruenc e matrix . The ruence matrices are usually generated in separate step (see Chapter 2). The ruence maps th us generated need to b e discretized and decomp osed in to shap e matrices. F or example, for a 5 4 grid, consider the follo wing ruence matrix 4 . 1 r 3 . 9 r 3 . 2 r 0 . 1 r 0 . 9 r 6 . 1 r 2 . 9 r 0 r 3 . 2 r 3 . 8 r 0 . 8 r 0 r 4 . 1 r 4 . 2 r 3 . 1 r 0 r 3 . 1 r 5 . 8 r 3 . 8 r 2 . 9 r T o allo w a decomp osition of the ruence map in to a manageable n um b er of shap e matrices, the b eamlet in tensities in a giv en b eam are rst discretized using an y sc heme dened b y the user. Using a simple round-o sc heme for the in tensit y of the b eamlets to the nearest in teger, w e get
PAGE 78
67 4 r 4 r 3 r 0 r 1 r 6 r 3 r 0 r 3 r 4 r 1 r 0 r 4 r 4 r 3 r 0 r 3 r 6 r 4 r 3 r Fluence map discretization is a decision problem in its o wn (see Bortfeld et al . [ 14 ]). In our study , w e will discretize the ruence map using the follo wing standard sc heme. The user sp ecies a lev el of discretization l % and b eamlet in tensities need to b e adjusted for this lev el of discretization. Let u b = max f x i : i 2 N b g denote the maxim um in tensit y of an y b eamlet in an y b eam b . F urther, let b = ( l = 100) u b denote the discretization step for b eam b . Then, discretized in tensit y of eac h b eamlet in b eam b is in teger m ultiple of b closest to b eamlet in tensit y in FMO. F or example, if u b = 100 for some b eam b and the user sp ecies the lev el of discretization as 10%, then b = 10. Th us, the in tensit y of eac h b eam in b m ust b e an in teger m ultiple of 10. In our study , w e will use the discretization lev el of 10%. The radiation generated in the linear accelerator pro vides uniform ruence for all bixels. T o generate a non-uniform ruence matrix from a uniform ruence source, a MLC system is used. A MLC has m ro ws, called channels , and eac h ro w has a left le af and a right le af , whose p ositions can b e c hanged and the radiation can pass in b et w een the left and righ t lea v es. If I has n columns 1 ; 2 ; : : : ; n , then for eac h ro w i , there are n + 1 p ositions, 1 ; 2 ; : : : ; n; n + 1, at whic h the left and righ t lea v es can b e p ositioned (see Fig. 4{2 .). If the left leaf is at p osition l and the
PAGE 79
68 righ t leaf is at p osition r , then the radiation will pass through the bixels n um b ered l ; l + 1 ; : : : ; r 1. r nr -r 3r r 3r r 1r r 2r r nr -r 2r r nr -r 1r r nr r bixelsr r leafr r positionsr r 1r r 2r r 3r r nr -r 2r r nr -r 1r r nr r n+1r r Figure 4{2: A single c hannel(ro w) of a m ultileaf collimator with left leaf at p osition l = 3 and righ t leaf at p osition r = n 2 : Radiation P asses through the bixels 3 ; : : : ; n 3 : Eac h c hoice of the left and righ t lea v es in ev ery ro w can b e sp ecied b y a 0-1 matrix called a shap e matrix . A \0" in the shap e matrix indicates that the corresp onding bixel is blo c k ed b y the lea v es and do es not deliv er an y radiation, and a \1" indicates that the corresp onding bixel deliv ers the radiation. Eac h shap e matrix satises the consecutiv e 1's prop ert y that all the 1's in eac h ro w are consecutiv e. T o deliv er the giv en ruence matrix I , the linear accelerator sends a uniform b eam of radiation through MLCs with dieren t shap e matrices S 1 ; S 2 ; : : : ; S K for dieren t lengths of times x 1 ; x 2 ; : : : ; x K , called b e am-on times (BOT, whic h is measured in terms of monitor units, MUs) suc h that I = P Kk =1 S k x k : In Fig. 4{3 , w e sho w a realization of a giv en ruence matrix through three shap e matrices. r r 4r r 4r r 3r r 0r r r r r 1r r 1r r 0r r 0r r r r r 1r r 1r r 1r r 0r r r r r 0r r 0r r 1r r 0r r r r 1r r 6r r 3r r 0r r r r r 0r r 1r r 1r r 0r r r r r 1r r 1r r 0r r 0r r r r r 0r r 1r r 0r r 0r r r r 3r r 4r r 1r r 0r r r = 3 r r r 0r r 1r r 0r r 0r r r +1r r r 1r r 1r r 1r r 0r r r +2r r r 1r r 0r r 0r r 1r r r r 4r r 4r r 3r r 0r r r r r 1r r 1r r 1r r 0r r r r r 1r r 1r r 0r r 0r r r r r 0r r 0r r 0r r 0r r r r 3r r 6r r 4r r 3r r r r r 0r r 1r r 1r r 1r r r r r 1r r 1r r 1r r 0r r r r r 1r r 1r r 0r r 0r r r Figure 4{3: Realizing a ruence matrix using three shap e matrices The total deliv ery time is the sum of the b eam-on times plus the setup times needed to go from one shap e matrix to another shap e matrix. Let c( S k ; S k +1 )
PAGE 80
69 denote the setup time needed to go from the shap e matrix S k to the shap e matrix S k +1 . In the ab o v e example, w e need three setups and its b eam-on time is 6 units. A p ertaining optimization problem to determine the shap e matrices S 1 ; S 2 ; : : : ; S K and their corresp onding BOTs is denoted as minimize K X k =1 x k + K 1 X k =0 c ( S k ; S k +1 ) sub ject to (MLC 1 ) K X k =1 S k x k = I (4.1) S k is a valid shap e matrix and x k > 0 for al l k = 1 ; 2 ; : : : ; K : (4.2) where S 0 is the initial setup cost whic h can corresp ond to the MLC b eing completely closed. W e refer to this problem as the total delivery time pr oblem . Note that in this form ulation all v ariables are required to b e strictly p ositiv e since only for those v ariables, setup costs ha v e to b e considered. In particular, K , whic h denotes the n um b er of shap e matrices used in the decomp osition, is a part of the output of the problem. A sp ecial case of (1) arises where setup times are assumed to b e negligible compared to the BOTs and th us can b e ignored. This problem is referred to as the c ol limator se quencing pr oblem with zer o setup times , or alternativ ely , the BOT pr oblem . Since BOT is measured in terms of MUs, a plan with less BOT is also referred as a MU ecien t plan. In tuitiv ely , it seems that the ob jectiv e of minimizing BOT should facilitate minimizing the n um b er of setups (and vice v ersa). Ho w ev er, empirical evidence indicates that after a limit, one has to compromise on BOT to reduce the n um b er of setups (and vice v ersa). Hence, nding ecien t deliv ery (i.e. minimizing BOT and n um b er of setups sim ultaneously) is a multi-criteria pr oblem . W e w ould lik e to p oin t out here that the imp ortance of these ob jectiv es compared to eac h other
PAGE 81
70 is not clearly dened. Moreo v er, the imp ortance of these factors migh t v ary from man ufacturer to man ufacturer. 4.1.2 Literature Review The approac hes dev elop ed for solving the MLC sequencing problem in literature can b e broadly classied in 2 categories based on the ob jectiv e; 1) Minimizing T otal Beam-On Time Minimization of the total BOT has b een the sub ject of n umerous in v estigations. See, for example, Bortfeld et al . [ 14 ], Sio c hi [ 72 ], Langer et al . [ 39 ], Boland et al . [ 11 ], Ah uja and Hamac her [ 3 ], Kamath et al . [ 35 , 36 ]. Because of ph ysical limitations, the MLCs cannot deliv er an y random pattern generated b y the sequencing algorithm. F or example, b ecause of the leaf structure that is used to form the ap ertures, in eac h ro w of the b eamlet grid, the b eamlets that are exp osed should b e c onse cutive (ro w-con v exit y). There are sev eral system-sp ecic design constrain ts lik e in terdigitation, minim um leaf separation, maxim um leaf spread, and maxim um o v ercen ter limit. Based on some imp ortan t constrain ts accoun ted for, most MLC related articles can b e classied in follo wing groups: 1. Unconstrained Case: Bortfeld et at. [ 14 ], Ah uja and Hamac her [ 3 ], and Kamath et al . [ 35 ] discussed algorithms to solv e this problem optimally . It is in teresting to note that the b oth of the latter algorithms pro duce iden tical result as the m uc h discussed sw eep algorithm prop osed b y Bortfeld et al . [ 14 ]. Ah uja and Hamac her [ 3 ] form ulated the BOT problem as a net w ork-ro w problem. Their solution metho d sho ws that there are n umerous w a ys in whic h the shap e matrices can b e obtained for an y giv en ruence matrix whic h yields the same BOT. Th us, elucidating the fact that the BOT problem solution space is highly degenerate. Ma et al . [ 42 ] sho w ed ho w this nature of the problem can b e exploited to decomp ose a dose prole with limited n um b er of setups without compromising on the BOT used.
PAGE 82
71 2. With In terdigitation Constrain t: In some systems, the o v erlapping of the left and righ t lea v es of adjacen t b eamlet ro ws, often called inter digitation (or in terleaf motion or leaf collision constrain ts), is not allo w ed. Kamath et al . [ 35 ], and Boland et al . [ 11 ] prop osed algorithms that can solv e the BOT problem ecien tly in the presence of in terleaf motion constrain ts. Kamath et al . [ 35 ] also sho w ed that, in this case, the optimal unidir e ctional plan (i.e. lea v es mo v e only from left to righ t and nev er bac ktrac k) is no less MU ecien t than an y bidir e ctional plan (i.e. where lea v es can mo v e b oth left to righ t and righ t to left). 3. With In terdigitation and T ongue and Gro o v e Constrain ts: In most of the commercial systems, there is a T ongue and Gr o ove (T&G) arrangemen t at the in terface b et w een the lea v es. The area under this section can get seriously underdosed b ecause of the mec hanical arrangemen t. Kamath et al . [ 36 ] sho w ed that the BOT problem can b e solv ed ecien tly with the T&G constrain ts if the leaf motion is constrained to b e unidirectional. With the in terdigitation and T&G constrain ts, the algorithm b y Kamath et al . [ 36 ] is equiv alen t to the Ro d Push-and-Pull algorithm dev elop ed b y Sio c hi [ 72 ]. Th us, it also pro vides the pro of of optimalit y in BOT for the Ro d Push-andPull algorithm under in terdigitation and T&G constrain ts for unidirectional plans. A unidirectional plan is preferred o v er a bidirectional plan as simplicit y of the structure of leaf mo v emen t often results in less setup time; ho w ev er, a bidirectional plan allo ws greater rexibilit y and could pro vide smaller BOT and a smaller n um b er of segmen ts. Though the optimal unidirectional plan is not inferior to an y bidirectional plan with in terdigitation constrain ts imp osed, this migh t not hold true when the T&G constrain ts are imp osed. Consider the follo wing example,
PAGE 83
72 0 1 2 2 1 0 F or the ruence map in the matrix ab o v e, the b eamlets in second column should b e op ened together to a v oid T&G violation (NoteIn this example, w e do not imp ose in terdigitation constrain t). F or a unidirectional plan, the optimal decomp osition w ould b e r r = r r r + r r r r +r r r r 0r r 1r r 2r r 2r r 1r r 0r r 0r r 1r r 1r r 1r r 1r r 0r r 0r r 0r r 0r r 1r r 0r r 0r r 0r r 0r r 1r r 0r r 0r r 0r r Ho w ev er, for a bidirectional plan, the dose can b e deliv ered more ecien tly as, r = r r r + r r r r 0r r 1r r 2r r 2r r 1r r 0r r 0r r 1r r 1r r 1r r 1r r 0r r 0r r 0r r 1r r 1r r 0r r 0r r Curren tly , there do not exist an y ecien t metho ds to solv e the BOT problem with the T&G constrain ts for the bidirectional case. In this c hapter, w e prop ose a metho d to solv e this problem. 2) Minimizing T otal T reatmen t Time W e ha v e only discussed the MU eciency of the dose deliv ery so far. A problem whic h is of similar imp ortance in dening a \go o d" deliv ery plan is the n um b er of segmen ts used. In other w ords, it is desirable to reduce the n um b er of setups also. High n um b er of setups results in longer treatmen t time as it tak es nite time to mo v e from one setup to another. This includes the time required to mo v e the lea v es and also the v erication and recording (V&R) time. T o a v oid in v olun tary mo v emen t of patien t during long treatmen ts and high setup errors, few er n um b er of setups are desirable. A high n um b ers of setups also result in longer qualit y assurance (QA) time. It should b e noted here that one of the h urdles in applying IMR T in clinics is that a high QA time results in a higher time requiremen t of the p ersonnel in v olv ed.
PAGE 84
73 The algorithms discussed ab o v e do not consider the n um b er of setups needed to deliv er the ruence. It has b een sho wn b y Burk ard [ 18 ] and W o eginger [ 6 ] that the problem of minimizing n um b er of setups for MLC sequencing problem is NP-hard (a class of problems for whic h ecien t algorithms are not a v ailable and it is b eliev ed b y Op erations Researc h comm unit y that suc h algorithms do not exist at all); hence, it is unlik ely to b e solv able in p olynomial time. Langer et al . [ 39 ] prop osed a mixed-in teger programming form ulation of this problem. Unfortunately , this approac h cannot b e used to solv e problems of the magnitude whic h o ccur in real practice. Th us, heuristics ha v e b een prop osed to solv e the MLC problem whic h minimize the n um b er of setups required while also k eeping a tab on MU requiremen t of the plan. Xia and V erhey [ 85 ] and Sio c hi [ 72 ] prop osed suc h heuristics. The latter heuristic has b een implemen ted as the IMF AST pro cedure as part of the commercial radiation system COR VUS. These heuristics are, ho w ev er, based on iterativ ely extracting ap ertures from a giv en ruence matrix in an ad-ho c manner. As a result, the ap ertures extracted migh t not satisfy all the constrain ts (suc h as in terleaf motion and T&G constrain ts). In that case, the ap ertures require to b e c hanged iterativ ely to nd a solution whic h do es not violate an y constrain ts. Since the constrain t violation is dealt after the ap ertures are extracted and these t w o steps are not in tegrated, this could result in extracting v ery inecien t shap es. Moreo v er, in these metho ds, one needs to implicitly en umerate a large n um b er of solutions to nd a go o d solution. Th us, an ecien t implemen tation of the algorithms is desirable in whic h solutions can b e en umerated ecien tly . W e presen t a net w ork-based form ulation to in tegrate the extraction pro cess with constrain t enforcemen t. W e ha v e also dev elop ed v ery ecien t metho ds to p erform the extraction pro cess on these net w orks. This form ulation can em ulate most of the existing metho ds to solv e the MLC sequencing problem and can also
PAGE 85
74 p oten tially impro v e them. Our form ulation can serv e as a unied structure to compare the dieren t form ulations in the literature. W e also prop ose a no v el metho d to extract the ap ertures. Our metho d outp erforms the b est kno wn metho d in literature (Xia and V erhey [ 85 ]) b oth in terms of b eam-on-time and n um b er of setups used. The run time of our algorithm is on the order of fraction of a second. 4.2 Net w ork-Based Mo del In this c hapter, w e presen t a net w ork-based mo del to solv e the MLC sequencing problem. In the remainder of this c hapter, w e will use the terms network and gr aph in terc hangeably . 4.2.1 Shap e Matrix Graph Boland et al . [ 11 ] prop osed a net w ork mo del to form ulate the BOT problem. W e ha v e adopted their net w ork in our w ork. F or the sak e of completeness of this w ork, w e will discuss their mo del briery here. The MLC sequencing problems is translated on a graph whic h w e will refer to as the shap e matrix graph. The shap e matrix graph G s is dened in the follo wing w a y . The graph consists of m la y ers corresp onding to the m ro ws in the shap e matrix. Eac h ro w i con tains no des indexed as ( i , l , r ) where l and r represen t the p ossible p ositions of the left and the righ t lea v es resp ectiv ely . A no de ( i , l , r ) represen ts that in ro w i , the radiation will pass through the bixels n um b ered l ; l + 1 ; : : : ; r 1 and other bixels in ro w i will not b e exp osed. Since a left leaf can nev er mo v e past the righ t leaf in the corresp onding ro w, l r should hold. Th us, giv en that the left leaf can stop at an y of the n + 1 p ositions and righ t leaf can stop at an y p osition l ; : : : ; n + 1, eac h ro w will ha v e 1 = 2 ( n +1)( n +2) no des. The set of no des N s is N s := f ( i; l ; r ) : i = 1 ; : : : ; m; l = 1 ; : : : ; n; r = l ; : : : ; n; n + 1 g [ f S our ce; S ink g Eac h no de in the net w ork (excluding the sink no de) is connected to the no des in the next la y er. So, for i = 1 ; : : : ; m 1, the set of arcs b et w een la y ers is
PAGE 86
75 A ( i ) := f ( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 )) : ( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 ) 2 N s g : Sourcer 111r 113r 112r 133r 123r 122r 211r 213r 212r 233r 223r 222r Sinkr 411r 413r 412r 433r 423r 422r 311r 313r 312r 333r 323r 322r r 1 0r0 1r1 1r1 0r Figure 4{4: An example of shap e matrix graph. The inter digitation constrain ts dictate that a no de ( i; l ; r ) cannot b e paired up with a no de ( i + 1 ; l 0 ; r 0 ) in the next la y er if ( l 0 > r ) or ( r 0 < l ). (see Fig. 4{5 ) r lr r rr r rr r lr r r lr r rr r r'r r l'r r Figure 4{5: Demonstrating in terdigitation constrain ts. Th us, to eliminate in terdigitation, w e simply do not include the arcs violating these conditions. The arc set A ( ) is re-dened as, A ( i ) := f (( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 )) : ( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 ) 2 N s ; l 0 r ; r 0 l g All the no des in the rst la y er are connected from the source no de and all the no des in last la y er are connected to the sink no de. So, A (0) := f ( sour ce; (1 ; l ; r )) : (1 ; l ; r ) 2 N s g A ( m ) := f (( m; l ; r ) ; sink ) : ( m; l ; r ) 2 N s g It can b e sho wn (Boland et al . [ 11 ]) that an y path on this net w ork corresp onds to a deliv erable ap erture and vice-v ersa. This graph has O( mn 2 ) no des and O( mn 4 ) arcs.
PAGE 87
76 W e will no w discuss the form ulation of BOT minimization using this graph. Sa y w e extract shap e matrices S 1 ; S 2 ; : : : ; S k from the graph and their in tensities x 1 ; x 2 ; : : : ; x k , then the minim um BOT problem can b e form ulated as minimize K X k =1 x k sub ject to (MLC 2 ) K X k =1 S k x k = I (4.3) S k is a valid shap e matrix and x k 0 for al l k = 1 ; 2 ; : : : ; K : (4.4) The shap e matrices S 1 ; S 2 ; : : : ; S k are dened b y paths from the source to the sink in the shap e matrix graph. The constrain ts ( 4.3 ) ensure that the desired ruence is deliv ered. This problem can b e form ulated (see Boland et al . [ 11 ]) as a Linear Program (LP). Since LPs are p olynomially solv able, the form ulation ab o v e establishes that minim um BOT problem is p olynomially solv able to o. It should b e noted here that LPs do not guaran tee in tegralit y and one migh t get non-in teger solution to o. It should b e noted here that the decomp osition problem has the in tegral asp ect b ecause of prepro cessing (discretizing) the ruence maps b efore decomp osition. The non-in teger solutions are (usually) implemen table and th us do not p ose a problem. The shap e matrix graph is v ery v ersatile and can p oten tially include all of the constrain ts related to sev eral commercially a v ailable MLC systems. Boland et al . [ 11 ] sho w ed ho w ro w-con v exit y and in terleaf motion (or in terdigitation) constrain ts can b e included in the graph. W e will no w sho w ho w some other signican t constrain ts can b e included in the net w ork mo del ecien tly .
PAGE 88
77 4.2.2 Inclusion of Connectedness Constrain ts By connectedness constrained ap ertures, w e are referring to ap ertures where it is not allo w ed to ha v e a closed leaf pair inside the rectangular b ounding b o x of the ap erture (commonly applied on systems with a \bac k-up" ja w conguration). This, in particular, means that the b eamlet ro ws that do not co v er all b eamlets should b e consecutiv e (see Fig. 4{6 ). r r a) Connectedness constrain t violated. b)Connectedness constrain t not violated Figure 4{6: Example of connectedness constrain t. This additional constrain t can b e incorp orated in to the shap e matrix graph b y sligh t mo dications. These mo dications are based on the in tuition dev elop ed b y the transforming the connectedness constrain t on the shap e graph as sho wn in Fig. 4{7 . The rst mo dication is that w e only include no des whic h corresp ond to ro ws that exp ose at least one b eamlet in the shap e matrix graph. Put dieren tly , w e remo v e all the no des of the form ( i; l ; r ) where r = l from the graph. Then, to represen t the fact that the bac kup ja ws will b e used to co v er the unexp osed ro ws, w e create arcs from the source no de to all no des in set N s , as w ell as from eac h no de in set N s to the sink. It is easy to see that eac h path in this net w ork corresp onds to a connected ap erture. The mo died net w ork con tains O ( mn 2 ) no des and O ( mn 4 ) arcs.
PAGE 89
78 r Node in row kr r Sourcer r Sink r r Node in row k+1r r Node in row k+pr r Closed rshape rrows not rallowedr r r Figure 4{7: Connectedness constrain t. 4.2.3 Inclusion of T ongue and Gro o v e Constrain ts In all the commercial MLC systems, to a v oid leak age b et w een the lea v es, there is a T ongue-and-Gro o v e (T&G) arrangemen t b et w een the lea v es. Tonguer -r andr -r Groove regionr Figure 4{8: T ongue and gro o v e arrangemen t. If the lea v es are not sync hronized, then the T&G region can get seriously underdosed. It can easily b e observ ed that this region is exp osed only when bixels at b oth ends of the region are not co v ered. Th us, it is desirable to exp ose the adjacen t bixels together to reduce the underdosing of the T&G region to the maxim um p ossible exten t. In other w ords, the T&G arrangemen t requires that for an y set of p ositions ( i , j ) and ( i +1, j ), the corresp onding bixels are exp osed for min( I i;j I i +1 ;j ) time together.
PAGE 90
79 In an y plan, the ruence deliv ery at an y orthogonal pair of p ositions in the bixel map is constrained in the follo wing manner. F or an y pair ( i; j ) and ( i +1 ; j ), Case 1 : If f 0 < I i;j I i +1 ;j g Then, T&G constrain ts dictate that all the ap ertures whic h include the b eamlet ( i; j ) should also include b eamlet ( i +1, j ). This can b e ac hiev ed b y dropping all the arcs in the shap e matrix graph whic h corresp onds to shap e matrices whic h include the b eamlet ( i , j ) but not b eamlet ( i +1, j ). In other w ords, w e need to remo v e all the arcs f (( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 )) : i = 1 ; : : : ; m 1 ; ( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 ) 2 N s ; l j < r ; l 0 > j g and f (( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 )) : i = 1 ; : : : ; m 1 ; ( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 ) 2 N s ; l j < r ; j r 0 g This mak es sure that the b eamlet ( i; j ) can only b e exp osed when b eamlet ( i + 1 ; j ) is exp osed. Moreo v er, the constrain ts ( 4.1 ) ensure that b eamlet ( i; j ) is exp osed for I i;j units. These result in exp osing the T&G region b et w een the b eamlets for I i;j units. Case 2 : If f I i;j > I i +1 ;j > 0 g Then, T&G constrain ts dictate that all the ap ertures whic h include the b eamlet ( i +1, j ) should also include b eamlet ( i; j ). This can b e ac hiev ed b y dropping all the arcs in the shap e matrix graph whic h corresp onds to shap e matrices whic h include the b eamlet ( i; j ) but not b eamlet( i +1, j ). In other w ords, w e need to remo v e all the arcs f (( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 )) : i = 1 ; : : : ; m 1 ; ( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 ) 2 N s ; l > j; l 0 j < r 0 g and f (( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 )) : i = 1 ; : : : ; m 1 ; ( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 ) 2 N s ; j r ; l 0 j < r 0 g This mak es sure that the b eamlet ( i +1, j ) can only b e exp osed when b eamlet ( i , j ) is exp osed. Moreo v er, the constrain ts ( 4.1 ) ensure that b eamlet ( i + 1 ; j ) is
PAGE 91
80 exp osed for I i +1 ;j units. These result in to exp osing the T&G region b et w een the b eamlets for I i +1 ;j units. Th us, there is a one to one corresp ondence b et w een an y feasible set of ap ertures not violating T&G constrain ts and corresp onding paths from source to sink in the graph. W e ha v e sho wn that the connectedness constrain ts and T&G constrain ts can b e implemen ted in the shap e matrix graph b y simply mo difying the shap e graph arc set. These mo dications main tain the linearit y of the problem. Therefore, the BOT problem can b e form ulated and solv ed as a LP . W e can conclude our result in the follo wing theorem. Theorem 1. The minimum BOT pr oblem is p olynomial ly solvable with or without any of the c onne cte dness, inter digitation, and T&G c onstr aints (even when the le af motion is not r estricte d to b e unidir e ctional). 4.3 Heuristics for Ecien t Deliv ery As discussed in sections ab o v e, only MU eciency do es not determine the go o dness of the plan deliv ery . W e dev elop a net w ork ro w based heuristic for the MLC sequencing problem whic h implicitly accoun ts for b oth b eam-on time and n um b er of setups used. 4.3.1 Extracting Shap e Matrices As noted ab o v e, most of the heuristics to solv e the MLC sequencing problem are based on extracting ap ertures from a giv en ruence map iterativ ely un til the map is exhausted. W e will form ulate these metho ds on our graph. Before w e pro ceed further, let us dene some notations here. F or an y arc e , let U B ( e ) = upp er b ound on the ro w through arc e , T ail ( e ) = the no de from whic h the arcs e emanates, H ead ( e ) = the no de in whic h the arc e terminates,
PAGE 92
81 P ositionR Leaf ( e ) = the p osition of righ t leaf in the shap e ro w corresp onding to no de tail node ( e ), P ositionLLeaf ( e ) = the p osition of left leaf in the shap e ro w corresp onding to arc no de tail node ( e ). Let the L ength of an y arc in the graph denote the width of the shap e ro w corresp onding to the no de from whic h the arc emanates. F or all the arcs emanating from source no de, w e will dene the length as 0. In, other w ords, Leng th ( e ) = 0 where f e : e 2 A [ sour ce ] g Leng th ( e ) = P ositionR Leaf ( e ) P ositionLLeaf ( e ) where f e : e 2 A [ i ] ; i = 1 ; : : : ; m ) g Sourcer 0r 111r 113r 112r 133r 123r 122r 211r 213r 212r 233r 223r 222r Sinkr 1r 411r 413r 412r 433r 423r 422r 311r 313r 312r 333r 323r 322r 1r 1r 2rFigure 4{9: An example demonstrating length of the arcs. As p oin ted in Section 4.2.1 , a path from source to sink corresp onds to a v alid shap e matrix. As the length of the arcs represen t the width of op ening of corresp onding shap e ro ws, the longest path from source to sink on this graph will represen t the largest deliv erable ap erture. Our net w ork has a la y ered structure
PAGE 93
82 and it is acyclic in nature (for pro of, refer Boland et al . [ 11 ]). On an acyclic graph, the longest path problem can b e solv ed using r e aching (otherwise kno wn as a breadth-rst searc h also) algorithm. W e will dene the reac hing algorithm no w. Sa y , a ruence matrix I is giv en, and w e w an t to extract the largest shap e matrix of in tensit y from it. The reac hing algorithm examines the no des la y er b y la y er ev aluating the b est candidate paths. F or an y no de p in la y er i , w e dene the predecessor of the no de to b e the no de in la y er i -1 whic h yields in the longest path to the no de p from the source. The length of the path to no de p is the sum of length of arcs on the path from source to the predecessor no de and the length of the arc connecting the predecessor no de to it. In this pro cess, w e consider only those arcs whose U B is no less than . As the rst step, w e initialize the information for the no des in la y er 1 b y setting P r edecessor (1 ; l ; r ) = sour ce ; P athLeng th (1 ; l ; r ) = 0 : F or ( i = 2 ; : : : ; m ), w e dene P r edecessor ( p ) = ar g max f Leng th ( e ) + P athLeng th ( tail ( e )) : e 2 A [ i ] and head ( e ) = p g ; U B ( e ) g P athLeng th ( p ) = Leng th f ( P r edecessor ( p ) ; p ) g + P athLeng th ( P r edecessor ( p )) And nally , for the sink no de P r edecessor ( sink ) = ar g max f Leng th ( e ) + P athLeng th ( m; l ; r ) : e 2 A [ m ] ; head ( e ) = ( sink ) g ; U B ( e ) g P athLeng th ( sink ) = Leng th f ( P r edecessor ( sink ) ; sink ) g + P athLeng th ( P r edecessor ( sink )) : It can b e noted here that P athLeng th ( sink ) represen ts the length of the longest path in the graph. The predecessors of the no des can b e used to trace this path and nd the corresp onding shap e matrix.
PAGE 94
83 The reac hing op eration in v olv es examining eac h arc in the net w ork once. Since w e ha v e O ( mn 4 ) arcs in the net w ork, the algorithm will tak e O ( mn 4 ) time. W e will no w discuss ho w shap e matrices can b e extracted on this graph while making sure that the correct amoun t of ruence is deliv ered and also ho w v arious constrain ts will b e incorp orated. 4.3.2 Enforcing V arious Constrain ts in the Net w ork Reac hing is a v ery ecien t metho d. Ho w ev er, it is imp ortan t to main tain the structure of the net w ork to apply this metho d. An y additional side constrain ts whic h are not a part of the graphical represen tation of the net w ork can mak e the problem substan tially dicult to solv e. In the w a y the shap e graph matrix is form ulated, the problem has pure net w ork structure for the r ow-c onvexity , inter digitation , and c onne cte dness constrain ts. One only needs to re-dene the arc sets to incorp orate these constrain ts. Ho w ev er, w e need to imp ose side constrain ts to ensure that correct ruence is deliv ered b y the shap e matrices extracted. W e will also sho w in app endix A that the T&G constrain ts cannot b e enforced simply b y re-dening the arc set here and th us w e need to in tro duce side constrain ts for these as w ell. These constrain ts destro y the net w ork structure, restricting us from solving these problems using the reac hing metho d. Before w e pro ceed further, let us dene some more notations here. As denoted ab o v e, I represen ts the total ruence matrix to b e deliv ered. In our sc heme, as w e extract ap ertures one b y one, the ruence left to b e deliv ered decreases in eac h iteration. Let x k denote the in tensit y of the k th shap e matrix, and S k denote the k th shap e matrix. Then, the ruence left to b e deliv ered b efore the k th iteration can b e denoted b y I k where I k = I k 1 -x k 1 S k 1 and I 1 = I. F or the case where w e extract only one ap erture at a time, the net w ork structure can b e restored as discussed b elo w.
PAGE 95
84 Ensuring Correct Amoun t of Fluence is Deliv ered . Our approac h to ensure that correct ruence is deliv ered using the ap ertures can b e view ed as a metho d whic h has t w o goals; 1. Ensure that no more than the sp ecied ruence is deliv ered through a b eamlet. 2. Keep deliv ering dose un til all the ruence is deliv ered. An y no de ( i; l ; r ) represen ts the shap e ro w i , in whic h b eamlets ( i; j ), j = l ; : : : ; r 1 are exp osed. A t an y iteration k , for the arcs emanating from this no de, w e will imp ose an upp er b ound of max f I k i;j , l j < r g . This will mak e sure that goal (1) is ac hiev ed. It can b e easily observ ed that these constrain ts are equiv alen t to stating that P Kk =1 S k x k I . Goal (2) can b e ac hiev ed b y terminating the extraction pro cess only when all the ruence to b e deliv ered is accoun ted for. It can b e easily observ ed that these constrain ts are equiv alen t to stating that P Kk =1 S k x k I . Th us, goals (1) and (2) imply P Kk =1 S k x k = I T ongue and Gro o v e Constrain ts . In Section 4.2.3 , w e sho w ed ho w the arc set in shap e matrix graph can b e mo died to incorp orate the T&G constrain ts. One basic dierence b et w een the BOT form ulation and the form ulation discussed in this section is that the BOT problem is solv ed to nd all the shap e matrices together. But, in our heuristic, the extraction of a shap e matrix is not completely coupled with shap e matrices extracted in other iterations. More explicitly , w e can note here that the constrain ts 4.3 are not enforced here in the w a y as in form ulation MLC 2 . It can b e sho wn that w e migh t lead in to an infeasible case here (see app endix A). So, w e need to do more than just mo dify the arc set to imp ose the T&G constrain ts. W e will use follo wing lemma to deriv e suc h a constrain t.
PAGE 96
85 Lemma 1 { If f 0 < I i;j I i +1 ;j g , i.e. for case (1) of T&G violation, the b eamlet ( i + 1 ; j ) can deliv er at most dif f ij = I i +1 ;j I i;j amoun t of dose using the ap ertures whic h do not include b eamlet ( i; j ). Pro of: W e will pro v e this b y con tradiction. Let d T N G ij = min f I i;j ; I i +1 ;j = I i;j g . Sa y , b eamlet ( i + 1 ; j ) deliv ers dif f ij + amoun t of dose using ap ertures whic h do not include b eamlet ( i; j ). Th us, the b eamlets ( i; j ) and ( i + 1 ; j ) will together b e exp osed to deliv er I i +1 ;j ( dif f ij + ) = I i +1 ;j I i +1 ;j + I i;j = I i;j = d T N G ij amoun t of dose. The T&G constrain ts require that the T&G region receiv es d T N G ij amoun t of dose. Th us, the T&G constrain t are violated. Pro v ed b y con tradiction. Using Lemma 1, w e can a v oid T&G violation b y imp osing the additional constrain t that dose deliv ered b y shap e matrices whic h include b eamlet ( i; j ) but not b eamlet ( i + 1 ; j ) should b e less than dif f ij . P e x e d ij where f e : (( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 )) where ( l j < r ) and ( l 0 > j or r 0 j ) g : Ho w ev er, it can b e noted here that in an algorithm where only one shap e matrix is extracted at eac h iteration, ro w on at most one of the arcs in eac h T&G constrain t will b e p ositiv e. Th us, these constrain ts can b e alternativ ely written as, x e d ij where f e : (( i; l ; r ) ; ( i + 1 ; l 0 ; r 0 )) where ( l j < r ) and ( l 0 > j or r 0 j ) g : In other w ords, w e will imp ose an upp er b ound of d ij on the arcs emanating from no des whic h con tain bixel ( i; j ) but not terminating in no des con taining bixel ( i + 1 ; j ). An equiv alen t result can b e sho wn for case (2) violation of T&G. Th us, these are the necessary conditions for T&G constrain ts. If these conditions are satised at eac h step, it can b e easily observ ed that w e will get a T&G violation free solution. Th us, the conditions are sucien t to o. W e can conclude this observ ation in the follo wing theorem:
PAGE 97
86 Theorem 2. The c onditions in L emma 1 ar e ne c essary and sucient for T ongue and Gr o ove elimination. The T&G constrain ts will b e imp osed b y mo difying the net w ork in the follo wing w a y: Case 1 : If f 0 < I i;j I i +1 ;j g Remo v e the arcs violating T&G as sho wn in Section 4.2.3 . Imp ose the upp er b ound of dif f ij = I i +1 ;j I i;j on the arcs corresp onding to shap e matrices whic h include b eamlet ( i; j ) but not b eamlet ( i +1 ; j ). i.e., UB ( e ) min f UB ( e ), di ij g where f e : (( i; l ; r ), ( i +1 ; l 0 ; r 0 )) where ( l j < r ) and ( l 0 > j or r 0 j ) g . Case 2 : If f I i;j > I i +1 ;j > 0 g Remo v e the arcs violating T&G as sho wn in Section 4.2.3 . Imp ose the upp er b ound of dif f ij = I i;j I i +1 ;j on the arcs corresp onding to shap e matrices whic h include b eamlet ( i +1, j ) but not b eamlet ( i; j ). i.e., UB ( e ) min f UB ( e ), di ij g where f e : (( i; l ; r ), ( i +1 ; l 0 ; r 0 )) where ( l > j or r j ) and ( l 0 j < r 0 ) g . It should also b e noted here that if multiple upp er b ounds are imp osed on an arc, the lowest of these n um b er will represen t the actual upp er b ound on the arc. 4.3.3 Sc hemes for Extraction One of the adv an tages of using the net w ork approac h is that the extraction is done in a w a y suc h that the constrain t enforcemen t is in tegrated with the pro cess. Th us, w e nd the \b est" deliv erable ap erture p ossible at an y stage. Using the metho ds describ ed ab o v e, w e nd the upp er b ound on the arcs in the net w ork to ensure that the correct ruence is deliv ered and T&G constrain ts are not violated. The longest path problem is then solv ed on the mo died graph and the shap e matrices are extracted. After eac h iteration, the ruence left to b e deliv ered c hanges. So, w e need to up date the corresp onding upp er b ounds on the arcs. The up dates can b e p erformed v ery ecien tly b y up dating only the part of the graph
PAGE 98
87 whic h no w violates the upp er b ound denitions. W e can summarize the extraction pro cess as, Step 0. Dene shap e graph (no de and arc sets). Find length ( e ) for eac h e 2 E s . Step 1. Find ruence b ounds and T&G b ounds for the arcs. Find the upp er b ounds of the arcs. Step 2. Extract an ap erture on the mo died graph using the reac hing metho d. Step 3. Up date the ruence data for the graph. If no ruence left to b e deliv ered, terminate. Else, go to step 1. In the literature, v arious metho ds ha v e b een prop osed to dene what a \go o d" ap erture is. In the follo wing sections, w e will discuss some of the p opular metho ds of extraction. W e also prop ose a no v el metho d of extraction here. Sw eep T ec hnique (or Ro d Push-and-Pull Algorithm). This m uc h discussed algorithm prop osed b y Bortfeld et al . [ 14 ] pro duces sequences whic h are optimal in BOT for the unc onstr aine d case. Sio c hi [ 72 ] extended the sw eep metho d to enforce in terdigitation and T&G constrain ts. This extended routine is called a r o d push-and-pul l metho d. It can b e sho wn that the ro d push-and-pull routine generates solution optimal in BOT ev en with in terdigitation constrain ts (Kamath et al . [ 35 ]). This metho d, ho w ev er, do es not accoun t for n um b er of setups used and could generate plans with a v ery high n um b er of setups. The sw eep tec hnique can b e implemen ted using our shap e matrix graph also, but w e will a v oid this discussion here since another ecien t implemen tation already exists (Kamath et al . [ 35 , 36 ]). W e will briery discuss the time complexit y of this algorithm in terms of the n um b er of basic op erations required in the algorithm. Ah uja and Hamac her [ 3 ] sho w ed that sw eep algorithm tak es O( mn ) time to nd the shap e ro ws for the matrices. These ro ws can b e com bined to generate a shap e matrix in O( mn *max f I i;j g ) time. Th us, in unconstrained case, sw eep algorithm on
PAGE 99
88 the whole will tak e O( mn *max f I i;j g ) time. F or the time complexit y discussion of the ro d push and pull algorithm, w e refer the readers to Kamath et al . [ 35 , 36 ]. [ The IMF AST r outine uses the r o d push-and pul l algorithm only as a subr outine in their pr o c e dur e. Sinc e sucient details for the implementation of the IMF AST r outine ar e not available in the liter atur e, we c ould not implement this metho d to c omp ar e our r esults. However, even with the limite d details we have ab out the IMF AST metho d, we c an state that this metho d c an b e implemente d using our shap e gr aph. In the IMF AST metho d, extr action is not inte gr ate d with c onstr aint enfor c ement. Our appr o ach c an impr ove the IMF AST metho d by pr oviding a to ol to do so. ] Areal Reduction . Xia and V erhey [ 85 ] sho w ed that their ar e al r e duction tec hnique, on a v erage, w orks b est in practice in terms of the n um b er of setups used. An indep enden t study (Que et al . [ 55 ]) whic h compared the areal reduction tec hnique to other metho ds also came to the same conclusion. In this metho d, the net w ork is extracted using the follo wing sc heme. F or an y extract S k , an in tensit y lev el k is decided. Then, the new extract is dened as S k ( i; j ) = 8><>: e ; I k i;j k 0 ; other w ise (4.5) where k = 2 m 1 ; m = r ound l og 2 ( max f I k i;j g ) : (4.6) Xia and V erhey [ 85 ] prop osed a metho dology to remo v e in terdigitation from their solutions. Ho w ev er, T&G constrain ts w ere not included in this metho d. W e ha v e extended the areal reduction tec hnique using our graph to include the T&G constrain ts also. Our implemen tation completes the areal reduction metho d b y pro viding a metho d to ecien tly com bine the shap e ro ws. The additional
PAGE 100
89 adv an tage that our metho d pro vides is that there is no need to c hange shap es after extraction as these already satisfy the constrain ts. The extraction pro ceeds in the follo wing w a y . F or an y iteration k, the in tensit y lev el k is calculated using expression ( 4.5 ). Then, using reac hing metho d, the largest shap e matrix of in tensit y k is extracted. The graph is then up dated and the pro cess is rep eated un til all the ruence is deliv ered. The areal reduction metho d w as implemen ted on the shap e matrix graph. Ho w ev er, unlik e the unconstrained case (and the case with in terdigitation), it migh t not b e p ossible to extract a shap e matrix of in tensit y k dened b y expression ( 4.6 ) alw a ys. Consider a 3 1 shap e matrix, where I 11 = 16, I 12 =10, and I 13 = 6. Since max f I ij g =16, using expression 4.6 , w e get k = 8. Ho w ev er, a deliv erable ap erture whic h can b e extracted from this column matrix can b e of in tensit y no more than 6 due to the T&G restrictions. So, for the case with T&G constrain ts imp osed, w e mo died the metho d to calculate the in tensit y v alue k . W e calculate k using expression ( 4.5 ) and c hec k if there exists a feasible shap e matrix of in tensit y k . If not, w e mo dify the v alue of k b y one of the t w o sc hemes 1. k k = 2 2. k k 1 W e k eep reducing the v alue of k using one of these sc hemes un til a deliv erable ap erture is found. With our empirical testing, w e found that sc heme (i) consisten tly w orks b etter than sc heme (ii). So, w e only rep ort the results of sc heme (ii) in the results section. Eac h step of areal reduction requires O( mn 4 ) time. Since eac h iteration reduces the in tensit y of at least one bixel b y unit amoun t, w e need to p erform at most max P i P j I ij iterations. Th us, areal reduction metho d tak es O( mn * P i P j I ij ). It can b e noted that in most of the extractions, the shap e matrices co v ers large n um b er of bixels and th us, this estimate is a lo ose upp er b ound. F or the areal
PAGE 101
90 reduction metho d, the time complexities remain the same as the unconstrained case since the constrain ts w ere enforced simply b y applying upp er b ounds on arcs whic h do not ha v e an y eect on the net w ork size. W e ha v e not discussed the time in v olv ed in prepro cessing b efore eac h iteration of extraction pro cess. But, w e w ould lik e to p oin t out here that one need not build the shap e graph at eac h iteration from the scratc h. By only incremen tally up dating the graphs at eac h iteration to rerect the c hanges, substan tial sa vings can b e gained in the run-time. Highest Fluence Reducing Shap es (HFRS) Metho d . Que et al . [ 55 ] sho w ed that areal reduction metho d p erforms v ery w ell in terms of the n um b er of segmen ts used. The reason b eing that the extraction pro cess do es not simply try to nd the biggest ap erture in size; it has another aim whic h is to nd an ap erture whic h can deliv er a large amoun t of ruence. W e prop ose a no v el metho d to extract matrices based on the in tuition that the more ruence the ap erture can deliv er, the more desirable it is. In other w ords, at eac h step, w e extract the ap erture that causes the maxim um reduction in the total ruence left to b e deliv ered. Let S k denote the largest shap e matrix of in tensit y that can b e extracted at iteration k . Then, in HFRS sc heme, k = ar g max =1 ;::;max f I k ij g ( P i P j S k ( i; j ) ) W e form ulate this problem as a sequence of longest path problems whic h can b e v ery ecien tly solv ed. The implemen tation of the HFRS algorithm is quite similar to the implemen tation of areal reduction tec hnique. Ho w ev er, unlik e the Xia-V erhey metho d, the in tensit y lev el of the ap erture is also a v ariable of the extraction pro cess. In brief, in an y iteration, w e extract ap ertures for all p ossible v alues of , and compare these to nd the one whic h one reduces the ruence left to b e deliv ered b y the most. In other w ords, the \v alue" of the ap erture will b e the length of the longest path m ultiplied b y the in tensit y lev el of the ap erture.
PAGE 102
91 F or HFRS, in eac h iteration, w e need to nd the longest path for max f I k i;j g shap e matrices. Th us, the algorithm will tak e O( mn * P i P j I ij *max f I i;j g ) time o v erall. With in terdigitation and T&G constrain ts, the time complexities remain the same as the unconstrained case since the constrain ts w ere enforced simply b y applying upp er b ounds on arcs whic h do not ha v e an y eect on the net w ork size. The results of this approac h are discussed in more detail in the results section. 4.4 Maxim um Leaf Spread Constrain t In a V arian system, the lea v es rest on a carriage. The width of the carriage op ening is restricted b y the length of the lea v es (see Fig. 4{10 ). Carriager Least Extended Leaf Most Extended Leafr Figure 4{10: An example of the leaf setup As can b e observ ed in the diagram ab o v e, the range of motion of lea v es dep ends on the carriage p osition and the width of the op ening. The restriction on the eld size renders it compulsory to use m ultiple settings of carriage to deliv er dose to a large eld. Th us, the curren t practice in the clinic is to divide this eld in to segmen ts and sequence these segmen ts indep enden tly to deliv er the dose. Let us denote the columns in the matrices b y c 1 ; c 2 ; : : : ; c n . The planner nds a w a y of splitting the matrix using past exp erience suc h that these matrices need no further splitting and also so that the splits c hosen will yield in to an o v erall ecien t
PAGE 103
92 deliv ery . F or example, let us consider the shap e matrix b elo w and assume that the lea v es can extend up to at most four units horizon tally . 2r r 1r r 1r r 0r r 1r r 2r r 3r r 0r r 1r r 2r r 3r r 1r r 0r r 3r r 0r r 2r r 0r r 1r r 1r r 2r r 0r r 1r r 2r r 1r r r Figure 4{11: Example of a shap e matrix. There are n umerous w a ys in whic h the sequences can b e split. F or example, the matrix ab o v e can b e split as f ( c 1 ; c 2 ; c 3 ), ( c 4 ; c 5 ; c 6 ) g , f ( c 1 ; c 2 ; c 3 ; c 4 ), ( c 5 ; c 6 ) g , or f ( c 1 ; c 2 ), ( c 3 , c 4 ; c 5 ; c 6 ) g . r 2r r 1r r 1r r 0r r 1r r 2r r 3r r 0r r 1r r 2r r 3r r 1r r 0r r 3r r 0r r 2r r 0r r 1r r 1r r 2r r 0r r 1r r 2r r 1r r 2r r 1r r 1r r 0r r 1r r 2r r 3r r 0r r 1r r 2r r 3r r 1r r 0r r 3r r 0r r 2r r 0r r 1r r 1r r 2r r 0r r 1r r 2r r 1r r r 2r r 1r r 1r r 0r r 1r r 2r r 3r r 0r r 1r r 2r r 3r r 1r r 0r r 3r r 0r r 2r r 0r r 1r r 1r r 2r r 0r r 1r r 2r r 1r r r Setting 1(split at center)r r Setting 2 (split leftr -r biased)r r Setting 3 (split rightr -r biased)r r r Figure 4{12: Demonstrating p ossible w a ys in whic h the matrix can b e split. In this c hapter, w e explore the idea that w e can k eep some common region b et w een the split regions and let the algorithm determine ho w and in whic h split the bixels will b e exp osed. F or example, for the example ab o v e, our algorithm will consider the carriage settings of ( c 1 ; c 2 ; c 3 ; c 4 ) and ( c 3 ; c 4 ; c 5 ; c 6 ). In general, let us assume that there are n columns and eac h carriage setting allo ws a width of p units. If n 2 p , i.e., t w o settings will b e sucien t to co v er the whole plan at the gan try angle under consideration. Then, the sets will denoted b y ( c 1 ; c 2 ; : : : ; c p ) and ( c n p +1 ; c n p ; : : : ; c n ). It can b e noted here that an y bixel in columns c n p +1 ; : : : ; c p can b e exp osed in either of the settings. The algorithm for the extraction can b e mo died easily to accommo date for these c hanges. W e create t w o shap e graphs,
PAGE 104
93 corresp onding to sets ( c 1 ; c 2 ; : : : ; c p ), and ( c n p +1 ; c n p ; : : : ; c n ). A t an y iteration k , w e dene matrices I k ' and I k ' ' suc h that I k 0 i;j = 8><>: I i;j if j p 0, otherwise and I k 00 i;j = 8><>: I i;j if j > n p 0, otherwise Then, w e dene the shap e graph for b oth matrices and nd the b est ap erture in b oth sets. The bigger ap erture of these t w o will b e c hosen to b e the next ap erture to b e included in the MLC leaf sequencing les. Both graphs will b e up dated after the extraction and the pro cess con tin ues un til all the ruence is deliv ered. In the general case, if the eld to b e irradiated is v ery wide, w e migh t need more that t w o settings. In suc h a case, the denition of the rst and last setting remains same. W e prop ose to lo cate the in termediate setting of width p units eac h in the follo wing w a y . The columns of the shap e matrix can b e divided in to k equal sets and then the midp oin t of eac h set will dene the p osition at whic h the in termediate settings will b e cen tered. r r Split 1 Split 3r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r Split 2r r Figure 4{13: An example of 3 split. Eac h in termediate setting t will b e cen tered at closest in teger of ( t 1) n p + n 2 p . See an example, for a case with three splits ( n =10, p =4), in Fig. 4{13 . The results for this approac h are discussed in Section 4.5 .
PAGE 105
94 4.5 Results and Discussion T o test our mo dels, w e ha v e studied cases of head-and-nec k cancer patien ts. 3D treatmen t planning data for 9 head-and-nec k cancer patien ts w ere exp orted from a commercial patien t imaging and anatom y segmen tation system (V o xelQ, Philips Medical Systems) at the Univ ersit y of Florida, Departmen t of Radiation Oncology . This data w as then imp orted in to the Univ ersit y of Florida Optimized Radiation Therap y (UF OR T) treatmen t planning system and using our in-house optimizer the ruence maps w ere generated. The ruence maps th us generated w ere discretized to generate the data for studying our leaf sequencing mo dels. 4.5.1 Comparison of Metho ds for Unconstrained Case T able 4{1 sho ws the results for the case when no constrain ts are imp osed. T able 4{1: Summary of results of the sequencing algorithms in unconstrained case. Case Sw eep Areal Reduction HFRS BOT #Segs BOT #Segs BOT #Segs 1 377.8 166 465.2 94 475.8 87 2 371.3 140 465.2 94 426.0 78 3 398.2 194 513.3 118 473.5 120 4 414.8 192 525.6 121 462.9 102 5 483.2 157 595.1 98 510.6 94 6 458.9 159 550.2 98 497.9 85 7 371.5 168 480.1 99 448.4 93 8 372.5 185 463.4 111 439.6 96 9 379.8 160 454.3 96 446.4 86 (Note: BOT = Beam-On-Time, #Segs = Num b er of segmen ts/setups used). In terms of b eam on time and MU eciency , the HFRS algorithm w as in most of the cases (8 out of 9) substan tially b etter than the metho d of Xia and V erhey . A t the same time, the n um b er of segmen ts w ere alw a ys equiv alen t or b etter than this metho d. Th us, the HFRS sc heme is the preferred approac h as it pro duces a small n um b er of ap ertures while k eeping a tab on b eam-on-time as w ell.
PAGE 106
95 4.5.2 Comparison of Metho ds With In terdigitation and T ongue and Gro o v e Constrain ts T able 4{2 sho ws results with in terdigitation constrain ts imp osed on all three approac hes. T able 4{2: Summary of results of the sequencing algorithms with in terdigitation. Case Sw eep Areal Reduction HFRS BOT #Segs BOT #Segs BOT #Segs 1 435.5 186 627.5 134 575.8 129 2 464.9 182 703.9 136 641.0 120 3 494.0 241 860.2 209 767.3 192 4 572.9 240 882.0 191 749.4 173 5 592.3 234 969.3 177 803.3 164 6 531.5 192 763.9 147 670.7 114 7 451.2 199 689.3 158 645.1 148 8 422.7 201 578.1 149 608.0 142 9 423.8 185 617.2 140 534.6 119 Adding in terdigitation to the problem magnied the dierences b et w een the Areal Reduction and HFRS approac hes. The run-times of all of the routines implemen ted as net w ork ro ws w ere on the order of a fraction of a second, ev en when deliv ery constrain ts w ere enforced. T able 4{3 sho ws results with in terdigitation and T&G constrain ts imp osed on all three approac hes. T able 4{3: Summary of results of the sequencing algorithms with in terdigitation and T&G correction. Case Sw eep Areal Reduction HFRS BOT #Segs BOT #Segs BOT #Segs 1 452.1 190 905.2 169 564.0 142 2 473.2 182 970.0 178 574.6 132 3 528.5 244 1248.4 261 729.1 203 4 605.2 254 1288.3 259 734.4 194 5 614.0 250 1384.5 241 756.8 175 6 539.0 191 1048.3 190 661.4 136 7 483.8 210 1039.1 204 596.0 152 8 441.9 212 925.4 206 537.9 156 9 439.2 195 893.4 176 550.6 136
PAGE 107
96 The HFRS sc heme clearly outp erforms the areal reduction metho d here. By observing T ables 4{2 and 4{3 , w e can conclude when T&G corrections are applied, the HFRS algorithm results in using higher n um b er of segmen ts compared to the case without T&G constrain ts. Ho w ev er, in terestingly enough, w e found that the BOT needed to deliv er the dose w as less with T&G constrain ts imp osed. Th us, w e see in terpla y of the role of these t w o ob jectiv es in designing an ecien t deliv ery plan here. As w e ha v e empirically noticed, this in terpla y is only so conspicuous when the t w o ob jectiv es are pushed to the limits. So, w e surmise that this observ ation comes from the fact that our results ha v e v ery tigh t b ounds to the implicit ob jectiv es aimed at. 4.5.3 Results for Maxim um Leaf Separation Constrain t W e studied sev en b eam head-and-nec k patien t case whic h needed t w o splits to deliv er the map. The results are sho wn in T able 4{4 . T able 4{4: Summary of results for maxim um leaf separation. Case HFRS+ Split in to 2 equal halv es HFRS+ Common region Sc heme BOT #Segs BOT #Segs 1 380.0 22 340.0 20 2 340.0 19 340.0 19 3 420.0 24 330.0 21 4 440.0 20 340.0 18 5 380.0 21 380.0 20 6 450.0 21 350.0 20 7 400.0 24 270.0 18 The results indicate that b y allo wing a common region b et w een the splits, the plan qualit y can b e substan tially impro v ed b y reducing b oth the MU requiremen t and n um b er of setups.
PAGE 108
97 4.6 Conclusions and F urther Researc h Directions In summary , w e found that the net w ork form ulation of the leaf sequencing problem could b e used to in tegrate constrain t enforcemen t during extraction of deliv erable ap ertures. W e ha v e th us dev elop ed an ecien t metho d to solv e the BOT problem when T&G constrain ts are also imp osed. Ho w ev er, this metho d migh t yield non-in teger solution. It remains an op en problem to nd if, for in tegral solution, the BOT minimization problem with T&G constrain ts imp osed is p olynomially solv able or not. In some MLC systems, setup time (i.e. the time required to mo v e from one set to another and the V&R time for eac h setup) could b e the dominan t factor in determining total treatmen t time. In this w ork, w e did not discuss ho w to accoun t for this factor. Our results demonstrate that our algorithm pro duces less n um b er of setups whic h implicitly reduces the setup time. Setup time also dep ends on the starting and ending conguration of the lea v es in an y setup. Th us, an algorithm whic h can nd an optimal order of setup can reduce this time signican tly . With giv en set of shap e matrices, w e can form ulate and solv e a T ra v eling Salesman Problem (see Ah uja et al . [ 2 ]) to nd the order for whic h the setup time will b e minim um. W e do not consider in ter-leaf leak age, leaf transmission or output v ersus eldsize eects in the net w ork algorithm discussed ab o v e. In tro duction of these eects will p ossibly lead in to an in tractable non-linear system. As a feasible alternativ e, after designing the plan, w e can re-calculate these factors and then p erform a quic k segmen t w eigh t optimization to restore the qualit y of the plan. In next c hapter, w e dev elop a column generation approac h to generate eectiv e ap ertures to impro v e an y giv en solution. The solution obtained b y the segmen t w eigh t optimization can
PAGE 109
98 b e used as an input to this routine and a few more ap ertures can b e added to ll in the qualit y gap as needed.
PAGE 110
CHAPTER 5 APER TURE MODULA TION 5.1 In tro duction Conformal radiation therap y seeks to conform the geometric shap e of the deliv ered radiation dose as closely as p ossible to that of the in tended targets. In c onventional conformal radiation therap y , this is tried to ac hiev e b y deliv ering a radiation dose b y using a relativ ely small n um b er of dieren t b eam shap es, called ap ertur es , eac h with uniform in tensit y lev el. Suc h ap ertures are formed b y partially blo c king the radiation source. In this w a y w e ma y , for example, create an ap erture whose shap e conforms to the b eam's-ey e-view of the targets in the patien t as seen from a certain b eam orien tation. Since the geometry of the targets and critical structure is dieren t in eac h patien t, customized ap ertures ha v e to b e man ufactured for eac h patien t. Despite this practical limitation, there are sev eral adv an tages to con v en tional conformal radiation therap y . Firstly , the fact that only a limited set of ap ertures is used leads to a short patien t treatmen t time. Shorter treatmen t times not only limit the discomfort of the patien t, but also ensure that the exp osure to radiation of the patien t as a whole is limited. Additional adv an tages lie in the abilit y to accurately deliv er a planned treatmen t, due to the relativ e ease with whic h the cum ulativ e dose receiv ed b y the patien t for a giv en set of ap ertures and in tensities can b e determined. Finally , treatmen t plans using only a few (relativ ely large) ap ertures are t ypically insensitiv e to patien t setup errors as w ell as patien t motion during treatmen t. 5.1.1 T raditional Tw o-Phase Approac h to T reatmen t Planning The ruence map decomp osition problem is a non-trivial optimization problem, and is often solv ed as a separate routine (see e.g. Bortfeld et al . [ 14 ]; Xia and 99
PAGE 111
100 V erhey [ 85 ]; Sio c hi [ 72 ]; Boland et al . [ 11 ]; Kamath et al . [ 35 ]; Ah uja and Hamac her [ 3 ]). In last c hapter, w e did an in-depth analysis of some of these mo dels. W e also prop osed v ery ecien t net w ork-based mo del for the decomp osition problem. T o allo w a decomp osition of the ruence map in to a manageable n um b er of ap ertures, the b eamlet in tensities in a giv en b eam are rst discretized to lev els that are t ypically c hosen to b e m ultiples of 5-20% of their maximal v alue. It is easy to see that this discretization ma y cause a signican t deterioration of the treatmen t plan qualit y , in particular when not in tegrated in to the FMO problem. Moreo v er, despite the discretization step, the n um b er of ap ertures required to deliv er the ruence map is often relativ ely large (on the order of 20-30 p er b eam, for a total of up to 140-210 ap ertures), whic h ma y cause undesirably long treatmen t times, as w ell as an unnecessarily high lev el of sensitivit y of the treatmen t plan to patien t setup errors or patien t motion. Finally , a t w o-stage approac h is usually rigid, in the sense that a sound trade-o b et w een treatmen t time and treatmen t plan qualit y cannot b e made due to the decoupling of the decomp osition phase from the treatmen t plan optimization phase. 5.1.2 Ap erture Mo dulation Approac h to T reatmen t Planning An approac h to the FMO problem that deals explicitly with ap ertures rather than individual b eamlets is called ap ertur e mo dulation . Recen tly , Sieb ers et al . [ 71 ] ha v e prop osed to incorp orate the ap erture c hoice in to a heuristic lo cal searc h metho d for FMO. Bednarz et al . [ 9 ] dev elop ed a mixed-in teger linear programming mo del that incorp orates a predetermined p o ol of ap ertures from whic h to c ho ose. Finally , Shepard et al . [ 68 ] ha v e prop osed a sim ulated annealing algorithm for nding a treatmen t plan using only a limited n um b er of ap ertures. These approac hes ha v e b een applied to head-and-nec k and prostate cancer cases, and empirical exp erimen ts suggests that one can obtain high-qualit y treatmen t plans with only a limited n um b er of ap ertures. Ho w ev er, all these approac hes are
PAGE 112
101 heuristic in nature, and cannot guaran tee that an optimal solution will b e found. The k ey to obtaining the b est ap ertures seems to b e to formally in tegrate the FMO and decomp osition problems. In this c hapter, w e prop ose a new optimization form ulation of the ap erture mo dulation approac h to FMO. This form ulation is based on a con v ex optimization form ulation of the traditional b eamlet-based FMO problem as dev elop ed in Chapter 2. Ho w ev er, in con trast with these earlier w orks, in this dissertation w e represen t a treatmen t plan explicitly in terms of ap ertures and their asso ciated in tensities. This mo del th us in tegrates the rst stage problem of obtaining a set of ruence maps (i.e., a set of individual b eamlet in tensities) with the second stage problem of decomp osing these ruence maps in to a set of deliv erable ap ertures and in tensities. Ho w ev er, the n um b er of allo w able ap ertures, and therefore the n um b er of decision v ariables and constrain ts in the new problem form ulation, is t ypically enormous. T o o v ercome this, w e prop ose a column generation approac h to solv e the con v ex programming problem. A column generation approac h requires the rep eated solution of a subproblem called the pricing problem. Solving this problem iden ties one or more suitable ap ertures to add to a p o ol of allo w able ap ertures. The v ery nature of our algorithm th us allo ws us to study the eect of adding ap ertures on the qualit y of the treatmen t plan, thereb y enabling a b etter trade-o b et w een the n um b er of ap ertures (i.e., treatmen t time) and treatmen t plan qualit y . One of the goals of our w ork is to determine whether high qualit y treatmen t plans ma y b e obtained using a limited n um b er of ap ertures, whic h w ould mean that the higher qualit y and eciency of IMR T treatmen ts can b e ac hiev ed without sacricing the adv an tages of con v en tional conformal radiation therap y men tioned ab o v e. The set of feasible ap ertures from whic h the pricing problem ma y c ho ose is dictated b y the particular c haracteristics of the MLC system that is used to
PAGE 113
102 deliv er the treatmen t. W e study three dieren t v arian ts of the pricing problem, eac h incorp orating a set of constrain ts on deliv erable ap ertures that is commonly encoun tered in one or more of the commercially a v ailable IMR T systems. W e pro vide algorithms that solv e these subproblems in p olynomial time in the n um b er of b eamlets. This c hapter is organized as follo ws. In Section 5.2 , w e dev elop our basic mo del form ulation. In Section 5.3 , w e then dev elop a column generation approac h to solving the optimization problem, including a detailed discussion of v arian ts of the pricing problems and solution metho ds to solv e these. In Section 5.4 , w e enhance the mo del, incorp orating additional constrain ts that can con trol the qualit y of the obtained treatmen t plans. In Section 5.5 , w e discuss the results of our approac h, and w e end the c hapter with some concluding remarks and directions for future researc h in Section 5.6 . 5.2 Mo del F orm ulation The metho ds describ ed in this c hapter ha v e b een published in Romeijn et al [ 58 ]. In our mo dels, w e will use the notations describ ed in Section 2.2 . W e need to dene some additional notation for our discussion of ap erture-based approac h. In our ap erture mo dulation form ulation of the FMO problem, w e assume that w e are giv en a set of ap ertures, sa y K , and c ho ose the decision v ariables to b e the in tensities of these ap ertures, sa y y k ( k 2 K ). Note that the total n um b er of p oten tial ap ertures that needs to b e included in the mo del is v ery large. F or example, consider an MLC that allo ws all com binations of left and righ t leaf settings. With the m n bixel grid for a single b eam, there are th us 1 2 n ( n 1) + 1 p oten tial com binations of righ t and left leaf settings for eac h b eamlet ro w, for a total of ( 1 2 n ( n 1) + 1) m 1 p oten tial ap ertures p er b eam angle. With n = m = 10 and 7 b eams, this w ould yield appro ximately 3 10 17 ap ertures to consider.
PAGE 114
103 Nev ertheless, w e will form ulate the problem taking all p oten tial ap ertures in to accoun t, and then prop ose a column generation approac h for solving the problem. Let K ` denote the set of allo w able ap ertures for b eam ` 2 B , let K [ ` 2 B K ` denote the set of all ap ertures, and let A k N denote the set of b eamlets that are exp osed in ap erture k 2 K . As a b eam of radiation passes through a patien t, it dep osits an amoun t of energy along its path that decreases with the distance tra v eled, but is linearly dep enden t on the in tensit y of the ap erture. Denoting the dose receiv ed b y v o xel j in structure s from ap erture k at unit in tensit y b y D k j s w e can express the dose receiv ed b y v o xel j in structure s as a function of the v ector of ap erture in tensities y 2 R j K j as P k 2 K D k j s y k . Next, w e will express the so-called dose dep osition co ecien ts D k j s in terms of the bixels comprising eac h ap erture. T o this end, w e rst denote the dose receiv ed b y v o xel j in structure s from b eamlet i (or, more accurately , an ap erture exp osing this b eamlet only) at unit in tensit y b y D ij s . All traditional b eamlet-based FMO form ulations then pro ceed b y assuming that the amoun t of dose dep osited to a giv en v o xel p er unit in tensit y of an y ap erture can b e determined b y adding the dose dep osited b y eac h of the individual b eamlets comprising the ap erture. Ho w ev er, this assumption is not en tirely accurate since the actual dose deliv ered to eac h v o xel b y a giv en b eamlet dep ends on the shap e of the ap erture in whic h it is con tained. This dep endence cannot b e incorp orated in to traditional b eamlet-based FMO mo dels since the decomp osition of the ruence maps in to a set of deliv erable ap ertures is done in a second phase, outside of the mo del. A ma jor additional adv an tage of using an ap erture-based approac h is that it is p ossible to impro v e on the accuracy of the dose calculation b y accoun ting for this dep endency . T o this end, w e dene nonnegativ e in tensit y correction factors r k for all ap ertures k 2 K , represen ting the dep endence of the deliv ered dose on the shap e of the ap erture. W e can then express the dose dep osition co ecien ts as D k j s = r k P i 2 A k D ij s , leading to the
PAGE 115
104 follo wing expression for the dose z j s receiv ed b y v o xel j in structure s : z j s = X k 2 K r k X i 2 A k D ij s ! y k j = 1 ; : : : ; v s ; s = 1 ; : : : ; S: Our basic mo del for FMO uses a v o xel-based ob jectiv e function, where a p enalt y is assigned to eac h v o xel based on the dose it receiv es under a giv en set of b eamlet in tensities. A con v ex p enalt y function for v o xel j in structure s is denoted b y F j s . T ypically , these functions are c hosen to b e piecewise p olynomial functions, allo wing for a dieren t p enalization of dose receiv ed ab o v e or b elo w some threshold. In addition, the same p enalt y function is often, though not necessarily , applied to all v o xels that are in the same structure. F or v o xels that are in critical structures and not in targets, the p enalt y function is commonly c hosen to b e one-sided, where the p enalt y is equal to zero for dose b elo w a threshold, and p ositiv e ab o v e the threshold. An ap erture mo dulation form ulation of the FMO problem then reads minimize S X s =1 v s X j =1 F j s ( z j s ) sub ject to (AM 1 ) X k 2 K r k X i 2 A k D ij s ! y k = z j s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S (5.1) y k 0 k 2 K : (5.2) Note that if all ap ertures that consist of a single b eamlet are deliv erable (whic h is t ypically the case), the ap erture mo dulation form ulation (AM 1 ) is equiv alen t to the b eamlet-based form ulation of the FMO problem as presen ted in Chapter 2. Apart from b eing able to solv e for ap ertures and asso ciated in tensities directly without a p ostpro cessing phase, this equiv alence will allo w us to test the h yp othesis that using relativ ely few ap ertures often yields near-optimal solutions to the FMO problem.
PAGE 116
105 5.3 A Column Generation Approac h 5.3.1 In tro duction As w as noted ab o v e, the n um b er of v ariables in this form ulation is m uc h to o large for practical purp oses. Ho w ev er, it can b e exp ected that in the optimal solution to this problem relativ ely few ap ertures will actually b e used (i.e., ha v e p ositiv e in tensit y). In fact, empirical evidence of this conjecture exists, see e.g. Sieb ers et al . [ 71 ], Bednarz et al . [ 9 ], Shepard et al . [ 68 ]. This means that if w e w ould b e able to iden tify a set of \go o d" ap ertures, w e could obtain a (near) optimal solution to the problem b y only incorp orating these ap ertures. Our approac h to solving (AM 1 ) is to use column generation, whic h is a formalization of this idea. In this metho d, w e start b y c ho osing a limited set of ap ertures, and solv e the restricted (AM 1 ) using only that set of ap ertures. (Implicitly , this means that w e restrict the in tensities of the remaining ap ertures to b e equal to zero.) Solving the restricted problem to optimalit y , w e use the Karush-Kuhn-T uc k er (KKT)-conditions for optimalit y to c hec k whether the corresp onding dual solution is feasible for the dual problem of the complete problem (AM 1 ). If so, w e ha v e obtained a pair of feasible primal and dual solutions with the same ob jectiv e function v alue, and therefore the optimal solution. If the dual solution is infeasible, an y dual constrain t that is violated pro vides a candidate ap erture that ma y b e added to the restricted v ersion of (AM 1 ). Using the common terminology from the simplex metho d and column generation for linear programming problems, w e sa y that these ap ertures price out, or are protable, in the sense that adding them to the restricted problem will (lik ely) yield an impro v ed ruence map. W e can reoptimize the new restricted problem, and rep eat the pro cedure. Belo w, w e formally deriv e an approac h that either v eries that a particular ruence map is optimal, or pro vides one or more protable ap ertures.
PAGE 117
106 5.3.2 The Pricing Problem The KKT-conditions (see, e.g., Bazaraa et al . [ 8 ], Hiriart-Urrut y and Lemarec hal [ 33 ]) for optimalit y of (AM 1 ) are: j s 2 @ F j s ( z j s ) j = 1 ; : : : ; v s ; s = 1 ; : : : ; S (5.3) k = S X s =1 v s X j =1 r k X i 2 A k D ij s ! j s k 2 K (5.4) y k 0 k 2 K (5.5) z j s = X k 2 K r k X i 2 A k D ij s ! y k j = 1 ; : : : ; v s ; s = 1 ; : : : ; S (5.6) k y k = 0 k 2 K (5.7) k 0 k 2 K (5.8) where j s is the Lagrange m ultiplier corresp onding to the dose denition constrain t ( 5.1 ) for v o xel j in structure s , and k is the Lagrange m ultiplier corresp onding to the nonnegativit y constrain t ( 5.2 ) of the in tensit y of ap erture k . No w supp ose that w e ha v e a solution to the ap erture mo dulation problem with only a limited n um b er of ap ertures included, sa y ^ K K . That is, w e ha v e a solution ( ^ z ; ^ y ) where the elemen ts of ^ y corresp onding to omitted ap ertures ha v e b een set to zero and the remaining v alues are optimal giv en these constrain ts. W e then c ho ose the corresp onding v alues of ^ j s according to KKT condition ( 5.3 ): ^ j s 2 @ F j s ( ^ z j s ) j = 1 ; : : : ; v s ; s = 1 ; : : : ; S: (5.9) i.e., ^ j s is a subgradien t of the p enalt y function for v o xel j in structure s at the curren t dose. Since the solution is optimal to the restricted problem, the pair (( ^ z ; ^ y ) ; ( ^ ; ^ )) (with ^ dened b y ( 5.4 )) satises the KKT conditions for the restricted form ulation, i.e., it satises conditions ( 5.3 ) and ( 5.6 ), and it satises conditions ( 5.4 ), ( 5.5 ), ( 5.7 ), and ( 5.8 ) for k 2 ^ K . The pair is optimal for the complete problem if and only the remainder of the KKT conditions is satised
PAGE 118
107 as w ell. Since for all k 62 ^ K w e ha v e that ^ y k = 0, conditions ( 5.5 ) and ( 5.7 ) are satised for all k 2 K . It th us remains to v erify whether conditions ( 5.4 ) and ( 5.8 ) are satised. These conditions can b e summarized as S X s =1 v s X j =1 r k X i 2 A k D ij s ! ^ j s 0 k 62 ^ K or equiv alen tly , min k 2 K S X s =1 v s X j =1 X i 2 A k D ij s ! ^ j s 0 : As noted ab o v e, the condition is satised for all k 2 ^ K , whic h allo ws us to extend the range o v er whic h is minimized to K . Also note that the in tensit y correction factors r k disapp ear from the problem. Represen ting an ap erture b y a binary v ector w , where w i = 1 if b eamlet i is in the ap erture and w i = 0 otherwise, w e wish to c hec k that min w 2 W S X s =1 v s X j =1 X i 2 N D ij s w i ! ^ j s 0 where w e optimize o v er the set W of all binary v ectors w represen ting an allo w able ap erture. Finally , regrouping terms w e ma y rewrite this as min w 2 W X i 2 N S X s =1 v s X j =1 D ij s ^ j s ! w i 0 : (5.10) If this condition is satised, the curren t solution is optimal. Otherwise, an y v ector w that violates this condition corresp onds to an ap erture that is lik ely to impro v e the ruence map when added to the restricted v ersion of (AM 1 ). Clearly , the condition can b e v eried b y solving the optimization problem on the left-hand side of equation ( 5.10 ) to optimalit y . When applied to large-scale linear programming problems, the latter optimization subproblem, whic h either v eries that the curren t solution is optimal or pro vides a previously omitted v ariable that w ould impro v e the solution when
PAGE 119
108 added to the problem, is called the pricing pr oblem , since it determines the reduced cost of eac h (omitted) primal v ariable (or, equiv alen tly , the shado w price of the corresp onding nonnegativit y constrain t). The in tuition b ehind the pricing problem for the ap erture mo dulation problem is that, for eac h bixel, w e compute the net unit eect of increasing its in tensit y b y some small v alue. If the aggregate of these eects is nonnegativ e for all p ossible ap ertures, the curren t solution cannot b e impro v ed and is optimal. Otherwise, an y ap erture with negativ e aggregate eect can b e added to the set of ap ertures and impro v e the solution. As a nal remark, note that equation ( 5.9 ) ma y allo w for a c hoice of v alues for ^ j s when the set of subgradien ts of the p enalt y function for v o xel j in structure s at the curren t dose con tains more than a singleton. No w recall that the pricing problem can conclude that the curren t solution is optimal if the optimal solution v alue to the pricing problem is nonnegativ e. Th us, due to the con v exit y of the functions F j s the b est c hoice for ^ j s will b e ^ j s = ( F j s ) 0+ ( ^ z j s ) i.e., ^ j s is the righ t deriv ativ e of the p enalt y function for v o xel j in structure s at the curren t dose. 5.3.3 Solving the Pricing Problem The dicult y of the pricing problem dep ends critically on the structure of the feasible region W . First, note that w e will only allo w ap ertures that op en b eamlets at a single b eam angle. This means that w e in fact obtain a smaller scale pricing problem for eac h b eam ` 2 B : minimize X i 2 N ` S X s =1 v s X j =1 D ij s ^ j s ! w i sub ject to w 2 W `
PAGE 120
109 where W ` represen ts the set of ap ertures that are allo w ed for b eam ` . Supp ose rst that an y conguration of b eamlets in a b eam forms an allo w able ap erture. In that case, the feasible region W ` reduces to W ` = f 0 ; 1 g j N ` j and the pricing problem reduces to minimize m X r =1 n X c =1 S X s =1 v s X j =1 D ( `;r ;c ) j s ^ j s ! w ( `;r ;c ) sub ject to w ( `;r ;c ) 2 f 0 ; 1 g r = 1 ; : : : ; m ; c = 1 ; : : : ; n: The pricing problem then decomp oses b y b eamlet, and can trivially b e solv ed b y , for eac h b eam, simply c ho osing all b eamlets for whic h the co ecien t in the pricing problem is negativ e, i.e., w ( `;r ;c ) = 8><>: 1 if P Ss =1 P v s j =1 D ( `;r ;c ) j s ^ j s < 0 0 otherwise : The running time for solving a pricing problem for a giv en b eam is then O ( mn ). Ho w ev er, ap ertures th us obtained ma y not b e deliv erable using a MLC, and a so-called sequencing phase is needed to decomp ose the ap ertures and in tensities in to truly deliv erable ap ertures and in tensities. In the next three sections, w e dev elop algorithms to solv e the pricing problems for three nested sets of constrain ts on deliv erable ap ertures. These constrain ts are common in one or more of the commercially a v ailable MLC systems. F or all systems, the ph ysical limitations of the MLC (in particular, the leaf structure that is used to form the ap ertures) imply that, in eac h ro w of the b eamlet grid corresp onding to a particular b eam, the b eamlets that are exp osed should b e c onse cutive . Figure 5{1 sho ws examples of three dieren t ap ertures that satisfy this constrain t. The rst system that w e will consider is an MLC system where this is the only requiremen t that ap ertures should satisfy . Suc h a system is th us able to deliv er all three ap ertures illustrated in Fig. 5{1 . In some systems, o v erlapping of
PAGE 121
110 the left and righ t lea v es of adjacen t b eamlet ro ws, often called inter digitation , is not allo w ed. The constrain ts disallo wing in terdigitation are sometimes also called in terleaf motion constrain ts or leaf collision constrain ts. The second system that w e consider disallo ws in terdigitation, but allo ws the lea v es to b e closed in the middle of the b eam to prev en t leaf collisions. Suc h a system can th us deliv er ap ertures suc h as the ones in Figures 5{1 (a) and (b), but not an ap erture as the one in Fig. 5{1 (c) (where in that Fig. the lea v es violating the in terdigitation constrain t are indicated b y the dashed o v al). Finally , the third system that w e consider do es not allo w the lea v es to b e closed in the middle of an ap erture, i.e., it requires the ap erture to b e formed b y a c onne cte d set of b eamlets. This system can only deliv er ap ertures of the t yp e illustrated in Fig. 5{1 (a). Allo wing in terdigitation . The most widely used commercial MLC system is a system that allo ws in terdigitation. In the simple approac h discussed ab o v e, the pricing problems decomp ose b y b eamlet. Ho w ev er, if w e tak e in to accoun t that an ap erture is formed b y pairs of left and righ t lea v es for eac h b eamlet ro w, this decomp osition is not v alid an ymore. W e can reform ulate the pricing problem for b eam ` b y letting c 1 ( r ) and c 2 ( r ) denote the index of the last b eamlet that is blo c k ed b y the left leaf and the rst b eamlet that is blo c k ed b y the righ t leaf in ro w r of b eam ` , resp ectiv ely . When in terdigitation is allo w ed, the pricing problem b ecomes minimize m X r =1 c 2 ( r ) 1 X c = c 1 ( r )+1 S X s =1 v s X j =1 D ( `;r ;c ) j s ^ j s ! sub ject to c 1 ( r ) < c 2 ( r ) r = 1 ; : : : ; m c 1 ( r ) 2 f 0 ; : : : ; n g r = 1 ; : : : ; m c 2 ( r ) 2 f 1 ; : : : ; n + 1 g r = 1 ; : : : ; m:
PAGE 122
111 (a) (b) (c) Figure 5{1: Three t yp es of ap ertures. Ap erture (a) can b e deliv ered b y all three considered systems. Ap erture (b) can b e deliv ered b y systems that allo w disconnected ap ertures. Ap erture (c) can only b e deliv ered b y systems allo wing in terdigitation.
PAGE 123
112 It immediately follo ws that the pricing problem for b eam ` 2 B in this case decomp oses b y b eamlet ro w. In particular, for b eamlet ro w r in b eam ` w e need to solv e the problem minimize c 2 1 X c = c 1 +1 S X s =1 v s X j =1 D ( `;r ;c ) j s ^ j s ! sub ject to (P) c 1 < c 2 c 1 2 f 0 ; : : : ; n g c 2 2 f 1 ; : : : ; n + 1 g : If the optimal v alue of this problem is zero, all bixels in the b eamlet ro w are blo c k ed. Otherwise, denoting the optimal solution to the problem b y ( c 1 ; c 2 ), b eamlets 1 ; : : : ; c 1 in b eamlet ro w r in b eam ` are blo c k ed b y the left leaf, and b eamlets c 2 ; : : : ; n in the same b eamlet ro w are blo c k ed b y the righ t leaf. F or b eam ` , the collection of b eamlets c hosen for eac h ro w in that b eam forms an ap erture for that b eam. Solving the pricing problem in a straigh tforw ard w a y for eac h b eamlet ro w using simple en umeration tak es O ( n 2 ) time. W e can th us nd a new ap erture that prices out (or determine that none exist at this stage) in O ( mn 2 ) time p er b eam. Ho w ev er, w e next presen t an alternativ e approac h that reduces the computational complexit y of the pricing problem. Note that w e are lo oking for a consecutiv e set of b eamlets in a giv en ro w for whic h the sum of their v alues is minimal. W e can nd suc h a set of b eamlets b y passing through the n b eamlets in a giv en ro w and b eam from left to righ t only once, k eeping trac k of (i) the cum ulativ e v alue o v er all b eamlets considered so far ( v ), and (ii) the maxim um cum ulativ e v alue found so far ( v ). No w note that, at an y p oin t, the dierence v v b et w een
PAGE 124
113 these t w o is a candidate for the b est solution v alue found so far. The algorithm is more formally giv en b y Step 0. Set v = v = v = 0, c 1 = c 1 = 0, c 2 = c 2 = 1. Step 1. Set v = v + P Ss =1 P v s j =1 D ( `;r ;c 2 ) j s ^ j s . Step 2. If v v , set v = v and set c 1 = c 2 + 1. If v v < v , set v = v , c 1 = c 1 , c 2 = c 2 . Step 3. If c 2 < n , incremen t c 2 and return to Step 1. Otherwise, stop; c 1 and c 2 giv e the optimal ap erture for ro w r in b eam ` , and v the corresp onding reduced cost. Bates and Constable [ 7 ] and Ben tley [ 10 ] ha v e sho wn that this algorithm is correct, and that its running time is linear in the n um b er of b eamlets in a ro w. W e conclude that w e can nd a new ap erture that prices out (or determine that none exist at this stage) in O ( mn ) time p er b eam when in terdigitation is allo w ed. Disallo wing in terdigitation . Some MLC systems do not allo w in terdigitation, that is, the left leaf in a ro w cannot o v erlap the righ t leaf of an adjacen t ro w. This additional constrain t causes a dep endency b et w een the ro ws in a b eam that prev en ts the pricing problem from decomp osing b y b eamlet ro w. In this Section, w e will form ulate a net w ork ro w problem that ecien tly solv es the pricing problem for a giv en b eam in case in terdigitation is not allo w ed. Our net w ork ro w form ulation is related to a net w ork ro w form ulation of the problem of decomp osing a giv en ruence map in to ap ertures that w as in tro duced b y Boland et al . [ 11 ]. In this net w ork, w e create a no de for eac h p oten tial pair of left and righ t leaf settings in eac h b eamlet ro w. F or a xed b eam ` 2 B , this means that w e ha v e no des of the form ( r ; c 1 ; c 2 ), where, as in discussion for the case whic h allo ws in terdigitation, r ( r = 1 ; : : : ; m ) indicates the b eamlet ro w, c 1 ( c 1 = 0 ; : : : ; n ) indicates the last b eamlet co v ered b y the left leaf in ro w r , c 2 ( c 2 = 1 ; : : : ; n + 1) indicates the rst b eamlet co v ered b y the righ t leaf in ro w r , and c 1 < c 2 . In
PAGE 125
114 addition, w e dene a source no de, sa y 0, and a sink no de, sa y m + 1. The total n um b er of no des in the net w ork is O ( mn 2 ). W e next create arcs b et w een no des in consecutiv e b eamlet ro ws if the t w o ro ws ma y app ear in a single ap erture. Note that this will allo w us to prev en t in terdigitation, b y simply not including arcs that w ould cause in terdigitation. F or r = 1 ; : : : ; m 1, this means that w e create a directed arc from no de ( r ; c 1 ; c 2 ) to no de ( r + 1 ; c 01 ; c 02 ) only when c 01 < c 2 and c 1 < c 02 , whic h ensures that no leaf collisions tak e place. In addition, w e create directed arcs from the source no de to all no des for ro w r = 1, and also from all no des for ro w r = m to the sink no de. The total n um b er of arcs in the net w ork is O ( mn 4 ). Figure 5{2 illustrates the structure of the net w ork for a small case with n = 2 ro ws and m = 2 columns. Figure 5{2: Net w ork for the pricing problem in case in terdigitation is not allo w ed.
PAGE 126
115 This net w ork pro vides a one-to-one corresp ondence b et w een paths from the source 0 to the sink m + 1 and ap ertures that satisfy the in terdigitation constrain ts. If w e no w let all arcs leading to a no de of the form ( r ; c 1 ; c 2 ) ha v e cost c 2 1 X c = c 1 +1 S X s =1 v s X j =1 D ( `;r ;c ) j s ^ j s and all arcs to the sink ha v e cost 0, the cost of eac h path is equal to the cost of the corresp onding ap erture in the pricing problem. The candidate ap erture for b eam ` can no w b e found b y solving for the minim um cost path from the source to the sink in the net w ork. Since the net w ork is acyclic, the shortest path problem can b e solv ed in time prop ortional to the n um b er of arcs in the net w ork (see, e.g., Ah uja et al . [ 1 ]). F or eac h b eam ` , w e can th us nd a new ap erture that prices out (or determine that none exist at this stage) in O ( mn 4 ) time. Requiring connected ap ertures . The previous paragraph solv es the pricing problem when ap ertures need to satisfy in terdigitation constrain ts, but are allo w ed to consist of a disconnected set of b eamlets. Ho w ev er, there exist systems that cannot deliv er suc h ap ertures. F or suc h systems, it is not allo w ed that lea v es close in the middle of the ap erture. This, in particular, means that the b eamlet ro ws that do not co v er all b eamlets should b e consecutiv e. This additional constrain t can b e incorp orated in to the net w ork giv en as discussed in the case whic h disallo ws in terdigitation with sligh t mo dications. The rst mo dication is that w e only represen t b eamlet ro ws that exp ose at least one b eamlet b y no des in the net w ork. Put dieren tly , w e remo v e all no des of the form ( r ; c 1 ; c 1 + 1) from the net w ork describ ed in the case whic h disallo ws in terdigitation. Then, to represen t the fact that the rst and last ro w(s) in an ap erture ma y b e closed, w e create arcs from the source 0 to all no des, as w ell as from eac h no de to the sink m + 1. It is easy to see that eac h path in this net w ork corresp onds to a connected ap erture. Moreo v er, if w e again let all arcs leading to
PAGE 127
116 no de ( r ; c 1 ; c 2 ) ha v e cost c 2 1 X c = c 1 +1 S X s =1 v s X j =1 D ( `;r ;c ) j s ^ j s and all remaining arcs ha v e cost 0, the cost of eac h path is equal to the cost of the corresp onding ap erture in the pricing problem. As the net w ork in the case whic h disallo ws in terdigitation, the mo died net w ork con tains O ( mn 2 ) no des and O ( mn 4 ) arcs. F or eac h b eam ` , w e can th us nd a new ap erture that prices out (or determine that none exist at this stage) in O ( mn 4 ) time b y solving for the minim um cost path from the source to the sink in this net w ork. 5.4 Extending the Mo del with P artial V olume Constrain ts A to ol that is commonly used b y ph ysicians to judge the qualit y of a treatmen t plan is the so-called (cum ulativ e) Dose-V olume Histo gr am (D VH) . This histogram sp ecies, for a giv en target or critical structure, the fraction of its v olume that receiv es at least a certain amoun t of dose. Man y successful optimization mo dels for FMO include some t yp e of constrain ts (often called D VH constrain ts or partial v olume constrain ts) on the shap e and lo cation of the D VH of a structure. In this section, w e will extend our basic mo del with a class of partial v olume constrain ts. 5.4.1 Mo deling P artial V olume Constrain ts W e prop ose to emplo y the follo wing t yp e of partial v olume constrain ts: (P 1 ) The aver age dose receiv ed b y the subset of a target of relativ e v olume 1 receiving the lowest amoun t of dose m ust b e at least equal to L . (P 2 ) The aver age dose receiv ed b y the subset of a structure of relativ e v olume 1 receiving the highest amoun t of dose ma y b e no more than U . If, for a giv en structure s , v s is in tegral, these partial v olume constrain ts can b e reform ulated as b ounds on the sum of the (1 ) v s largest or smallest v o xel doses in a structure. This concept has b een applied to facilit y lo cation mo dels b y Andreata and Mason [ 5 ], who prop osed c ho osing a facilit y lo cation
PAGE 128
117 for whic h the sum of the k largest distances b et w een customers and facilit y , called the k -eccen tricit y , is minimized. They sho w ed that the k -eccen tricit y is a con v ex function, and prop osed algorithms for solving the corresp onding lo cation problems on tree graphs (see also Slater [ 73 ] who studied similar problems in a more abstract graph theoretical setting). More recen tly , T amir [ 75 ] and Ogryczak and Za w adzki [ 48 ] ha v e dev elop ed ecien t optimization approac hes for v arious lo cation mo dels emplo ying suc h an ob jectiv e. F urthermore, Ogryczak and T amir [ 47 ] ha v e dev elop ed a more general approac h for ecien tly minimizing the sum of the k largest of n function v alues. When form ulating partial v olume constrain ts for the FMO problem, care m ust b e tak en when v s is non-in tegral. In this case, w e could round up or do wn the n um b er of v o xels that w e are in terested in and apply the metho ds referred to ab o v e. Ho w ev er, this could ha v e undesirable round-o eects, in particular for smaller structures. This is due to the fact that although w e mo del the FMO problem on a discretization of the patien t using a nite n um b er of v o xels, the underlying structures are con tin uous ob jects. Therefore, w e ha v e instead c hosen to use the related concept of tail-a v erages, whic h applies to b oth con tin uous and discrete p opulations. This concept has in recen t y ears receiv ed m uc h atten tion in the risk managemen t literature (see Ro c k afellar and Ury asev [ 56 , 57 ] and Ogryczak and Ruszczy nski [ 46 ]), and are often referred to as Conditional V alue-atRisk (CV aR) or T ail V alue-at-Risk (TV aR). Applying these concepts to discrete distributions allo ws for a correction for non-in tegral v alues of v s b y including a fractional v o xel in the tail a v erage, and can as suc h b e exp ected to m uc h b etter appro ximate the con tin uous tail-a v erage than a rounding approac h. In Chapter 2, w e ha v e demonstrated the p o w er of suc h constrain ts in the con text of FMO. In addition, in con trast with other t yp es of constrain ts on the D VH prop osed in the literature (see, e.g., Langer et al . [ 37 ]; Bortfeld et al . [ 15 ]; Carol et al . [ 20 ]; W u and
PAGE 129
118 Mohan [ 83 ]; Lee et al . [ 40 , 41 ]), these constrain ts can b e incorp orated in to a con v ex programming form ulation of the problem. Constrain ts of the t yp e (P 1 ) can b e form ulated as lo w er b ound constrain ts on the so-called lo w er CV aR at lev el (for short, the lo w er -CV aR). Ro c k afellar and Ury asev [ 56 ] and Ogryczak and Ruszczy nski [ 46 ] sho w that the lo w er -CV aR can b e expressed as s = max 2 R ( 1 (1 ) v s v s X j =1 ( z j s ) + ) Similarly , constrain ts of the t yp e (P 2 ) can b e form ulated as upp er b ounds on the upp er -CV aR, whic h is equal to s = min 2 R ( + 1 (1 ) v s v s X j =1 ( z j s ) + ) : F or eac h structure, w e ma y no w dene a set of b ounds on lo w and high tail a v erages. F or targets, these will b e in the form of lo w er b ounds L s for 2 A s , where A s is a (nite) subset of (0 ; 1), and upp er b ounds U s for 2 A s , where A s is a (nite) subset of (0 ; 1) ( s = 1 ; : : : ; T ). F or the critical structures, these will b e in the form of only upp er b ounds U s for 2 A s , where A s is a (nite) subset of (0 ; 1) ( s = 1 ; : : : ; S ). 5.4.2 An Enhanced Problem F orm ulation T o ensure that a feasible solution to our enhanced mo del exists, w e will not incorp orate hard D VH constrain ts but instead allo w them to b e violated and incorp orate a p enalization of suc h violations in to the ob jectiv e function. If S A S (0 ; 1) denotes the set of lo w er CV aR constrain ts, w e let G s denote a con v ex and nonincreasing p enalt y function asso ciated with the lo w er -CV aR constrain t for structure s . Similarly , if S A + S (0 ; 1) denotes the set of upp er CV aR constrain ts, w e let G s denote a con v ex and nondecreasing p enalt y function asso ciated with the upp er -CV aR constrain t for structure s . This leads to the
PAGE 130
119 follo wing ap erture mo dulation form ulation of the FMO problem: minimize S X s =1 v s X j =1 F j s ( z j s ) + X ( s; ) 2 S A + G s ( s ) + X ( s; ) 2 S A G s ( s ) sub ject to (AM 2 ) X ` 2 B X k 2 K ` X i 2 A k D ij s ! y k = z j s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S s + 1 (1 ) v s v s X j =1 z j s s + s ( s; ) 2 S A + (5.11) s 1 (1 ) v s v s X j =1 s z j s + s ( s; ) 2 S A (5.12) y k 0 k 2 K : Note that it is easy to see that there will exist an optimal solution to (AM 2 ) for whic h constrain ts ( 5.11 ) and ( 5.12 ) are binding. 5.4.3 The Pricing Problem The KKT-conditions for optimalit y of (AM 2 ) are j s 2 @ F j s ( z j s ) + X ( s; ) 2 S A + : j =1 ;::: ;v s ; z j s s s + X ( s; ) 2 S A : j =1 ;::: ;v s ; z j s s s j = 1 ; : : : ; v s ; s = 1 ; : : : ; S (5.13) s 2 @ G s ( s ) ( s; ) 2 S A + (5.14) s 2 @ G s ( s ) ( s; ) 2 S A (5.15) k = S X s =1 v s X j =1 X i 2 A k D ij s ! j s k 2 K y k 0 k 2 K z j s = X k 2 K X i 2 A k D ij s ! y k j = 1 ; : : : ; v s ; s = 1 ; : : : ; S
PAGE 131
120 s s + 1 (1 ) v s v s X j =1 z j s s + ( s; ) 2 S A + s s 1 (1 ) v s v s X j =1 s z j s + ( s; ) 2 S A k y k = 0 k 2 K k 0 k 2 K s 0 ( s; ) 2 S A + (5.16) s 0 ( s; ) 2 S A (5.17) where j s and k are dened as in the basic mo del, and s and s are the Lagrange m ultipliers of the CV aR constrain ts ( 5.11 ) and ( 5.12 ). Note that the sign conditions ( 5.16 ) and ( 5.17 ) are automatically satised b y the monotonicit y assumptions on the p enalt y functions G s and G s , resp ectiv ely . Then pro ceeding as in the basic mo del, giv en a solution ( ^ z ; ^ y ) to the problem using a limited set of ap ertures ^ K K , w e use ( 5.13 )-( 5.15 ) to determine corresp onding v alues of ^ . As in the basic mo del, if an y set of subgradien ts con tains more than one elemen t, the righ t deriv ativ e (i.e., the largest v alue in the set of subgradien ts) is c hosen. With these c hoices, w e obtain a pricing problem of exactly the same form as in the basic mo del, whic h can then b e solv ed using the same approac hes. 5.5 Computational Results T o study the eect of the n um b er of ap ertures used on the qualit y of the treatmen t plans obtained, w e ha v e studied cases of head-and-nec k cancer patien ts where the primary reason for using IMR T is to preserv e saliv ary function while obtaining adequate target co v erage. 3D treatmen t planning data for ten headand-nec k cancer patien ts w ere exp orted from a commercial patien t imaging and anatom y segmen tation system (V o xelQ, Philips Medical Systems) at the Univ ersit y of Florida Departmen t of Radiation Oncology . This data w as then imp orted in to the Univ ersit y of Florida Optimized Radiation Therap y (UF OR T) treatmen t
PAGE 132
121 planning system and used to generate the data required b y the mo dels describ ed in Sections 5.2 and 5.4 . Our treatmen t planning system generated a v o xel grid with v o xels of size 4 4 4 mm 3 for the targets and critical structures except for unsp ecied tissue, for whic h w e used v o xels of size 8 8 8 mm 3 . T able 5{1 sho ws the problem dimensions for the ten cases. As men tioned in Section 2.3.1 , there is no fundamen tal basis for quan tifying the qualit y of a treatmen t plan. Therefore, to assess the abilit y of our mo del to generate high-qualit y treatmen t plans, w e ha v e ev aluated our treatmen t plans using a set of external planning criteria that w ere not included in our mo del. Note that in order to b e able to mak e an ob jectiv e assessmen t of mo del qualit y , it is imp ortan t that these external planning criteria w ere not included in our mo del, to a v oid bias of the conclusions. In particular, w e ha v e used the criteria published in t w o studies of the Radiation Therap y Oncology Group [ 65 , 66 ] with minor mo dications (where the mo dications made the criteria more demanding). The main goal w as to treat t w o targets: Planning T arget V olume 1 (PTV1) and Planning T arget V olume 2 (PTV2), whic h con tains PTV1. PTV1 consists of actual iden tied disease, while PTV2 is an expansion of PTV1 that also con tains regions where the disease ma y ha v e spread, as w ell as a margin to accoun t for the p ossibilit y of patien t motion and setup uncertain t y . PTV1 should b e treated to a prescription dose lev el of D PTV1 R x = 73 : 8 Gra y (Gy), and PTV2 should b e treated to a prescription dose lev el of D PTV2 R x = 54 Gy . More precisely , this means that at least 95% of b oth targets should receiv e the prescrib ed dose, no more than 1% of b oth targets should b e underdosed b y more than 7%, and no more than 20% of PTV1 should b e o v erdosed b y more than 10%. Eac h of the up to four saliv ary glands (LPG: left parotid gland, RPG: righ t parotid gland, LSG: left submandibular gland, and RSG: righ t submandibular gland) are considered spared if either no more than 50% of the v o xels in the gland receiv e a
PAGE 133
122 dose exceeding 30 Gy , or the mean dose to the gland do es not exceed 26 Gy . No v o xels in the spinal cord (SC) should receiv e a dose exceeding 45 Gy , and no v o xels in the brainstem (BS) should receiv e a dose exceeding 54 Gy . Less than 1% of the v o xels in unsp ecied tissue should receiv e more than 65 Gy to a v oid side eects of treatmen t due to an unacceptable lev el of o v erdosing of that structure. F or all ten cases, w e designed plans using sev en equispaced b eams at nominal settings, i.e., their orien tation w as not optimized for the individual patien ts. T able 5{1: Mo del sizes. case # b eamlets # v o xels 1 1,017 17,471 2 1,195 19,543 3 1,745 36,136 4 1,804 38,724 5 1,182 17,466 6 1,015 14,147 7 1,178 22,300 8 1,090 19,417 9 921 15,416 10 2,073 30,664 The structure-dep enden t p enalt y functions w ere c hosen to b e piecewise-linear appro ximations to piecewise-p olynomial functions of the follo wing form: F j s ( z ) = 1 v s U s max(0 ; T s z ) p Us + O s max (0 ; z T s ) p Os where the sup erscript U refers to underdosing and the sup erscript O refers to o v erdosing. The co ecien ts U s ; O s 0 and the p o w ers p Us ; p Os 1 to ensure con v exit y of the function. Finally , the parameters T s represen t thresholds, to allo w the shap e of the p enalt y functions to b e dieren t for underdosing and o v erdosing. W e emplo y ed piecewise-linear p enalt y functions in order to allo w us to solv e our mo dels using a linear programming solv er. The piecewise-linear functions used t w o segmen ts for underdosing and four segmen ts for o v erdosing for targets, and t w o segmen ts for o v erdosing for critical structures. In addition, w e dened lo w er b ounds
PAGE 134
123 L s and upp er b ounds U s on the dose for eac h structure, with v ery steep slop es (in the range 10 8 -10 12 ) in the piecewise linear functions for dose v alues b elo w the lo w er b ound and ab o v e the upp er b ound. Based on three of the ten patien t cases and the b eamlet-based mo del for FMO, w e tuned the problem parameters b y man ual adjustmen t. This resulted in the parameter v alues sho wn in T ables 5{2 and 5{3 . Similar to Tsien et al . [ 76 ], w e found that high p o w ers of dose dierence lead to excellen t results. W e next solv ed the b eamlet form ulation of the FMO problem with partial v olume constrain ts for all ten patien t cases using this set of parameters to obtain b enc hmarks. T able 5{2: V alues of the co ecien ts of the v o xel-based p enalt y functions. s T s L s U s p Us U s O s p Os PTV1 D PTV1 R x + 2 : 5 D PTV1 R x 0 : 5 7,500 12 D PTV1 R x + 5 : 5 7,500 6 PTV2 D PTV2 R x + 2 0 : 93 D PTV2 R x 75,000 12 D PTV1 R x + 5 : 5 75,000 6 Tissue 35 0 0 { D PTV1 R x + 5 : 5 8,000,000 5 LPG 0 0 0 { D PTV1 R x + 5 : 5 350,000 4 RPG 0 0 0 { D PTV1 R x + 5 : 5 350,000 4 LSG 0 0 0 { D PTV1 R x + 5 : 5 150,000 4 RSG 0 0 0 { D PTV1 R x + 5 : 5 150,000 4 SC 28 0 0 { 40 150 3 BS 33.6 0 0 { 48 200 3 T able 5{3: V alues of the co ecien ts for the CV aR constrain ts. L ower CV aR-c onstr aints Upp er CV aR-c onstr aints s L s U s PTV1 0.90 0 : 97 D PTV1 R x 0.99 1 : 07 D PTV1 R x PTV2 0.90 0 : 96 D PTV2 R x 0.95 1 : 1 D PTV2 R x LPG { { 0.60 26 RPG { { 0.60 26 LSG { { 0.60 26 RSG { { 0.60 26 Next, w e solv ed all ten problem instances using the column generation approac h to the ap erture mo dulation form ulation of the FMO problem as dev elop ed
PAGE 135
124 in this c hapter. Since the most commonly used MLC system allo ws for in terdigitation, w e ha v e used the pricing algorithm dev elop ed for the case whic h allo ws in terdigitation. Recall that, when the column generation algorithm is ran un til con v ergence, the optimal b eamlet-based FMO solution is obtained. T o determine the n um b er of ap ertures required b y our approac h to obtain a high-qualit y treatmen t plan, w e examined the solution obtained b y the column generation algorithm in eac h iteration. W e then ev aluated whic h of the external planning criteria w ere satised for eac h in termediate solution. The results for eac h of the ten cases are summarized in T able 5{4 . In particular, this table pro vides the n um b er of ap ertures in our solution that ensured that eac h of the external planning criteria w as satised. Note that w e ha v e added sev eral criteria that are stronger than the ones generally used in clinical practice, suc h as stronger o v erdosing criteria for PTV1 (no more than 5% or 10% (instead of 20%) of its v olume ma y b e o v erdosed b y more than 10%) and unsp ecied tissue (no more than 5% of its v olume ma y receiv e more than 50 Gy). In case 1, 95% co v erage of PTV2 with the prescription dose of 54 Gy w as not obtained (indicated b y ? ), ev en using the b eamlet-based FMO form ulation. In this case, 94% co v erage w as obtained using 108 ap ertures. In cases 4 and 7, the mean dose constrain t could not b e satised (indicated b y ). Ho w ev er, in these cases the partial v olume constrain t w as satised using a limited n um b er of ap ertures, ensuring sparing of these structures. Cases where a saliv a gland w as either not spareable (since it w as largely con tained in the target) or not in the vicinit y of a target (and therefore not included in the problem) are indicated b y {. Note that the abilit y of the mo dels to obtain high-qualit y treatmen t plans (as measured b y the external planning criteria) with a single parameter setting suggests that parameter tuning (or, equiv alen tly , a m ulti-criteria solution approac h) is not needed for eac h individual patien t.
PAGE 136
125 The last line in the table corresp onds to the n um b er of ap ertures that w as found using a traditional t w o-stage approac h, where an in tensit y prole found b y our b eamlet-based FMO mo del w as discretized and decomp osed in to ap ertures. T able 5{4: Num b er of ap ertures needed to satisfy eac h criterion. Case s Criterion 1 2 3 4 5 6 7 8 9 10 PTV1 99% 68.6 Gy 9 7 10 10 10 6 10 12 12 12 80% 81.2 Gy 16 19 17 17 18 17 19 17 18 23 90% 81.2 Gy 47 55 55 48 47 46 51 50 49 59 95% 81.2 Gy 111 126 130 102 108 111 108 116 119 118 PTV2 95% 54 Gy ? 48 56 43 38 84 209 85 175 43 99% 50.2 Gy 117 45 18 30 88 65 55 19 77 21 Tissue 95% 50 Gy 51 27 76 31 35 92 26 43 14 105 99% 65 Gy 12 19 18 10 17 17 11 15 12 12 LPG 50% 30 Gy 28 28 10 10 20 19 10 24 12 12 mean 26 Gy 36 28 10 10 45 34 10 43 12 12 RPG 50% 30 Gy 28 { 10 10 10 17 67 12 28 39 mean 26 Gy 34 { 10 10 10 46 15 52 55 LSG 50% 30 Gy 36 29 { 33 37 { 33 { 30 33 mean 26 Gy 80 30 { 54 53 { 39 { 52 36 RSG 50% 30 Gy { { { 78 19 { { { { { mean 26 Gy { { { 22 { { { { { SC 100% 45 Gy 73 30 37 44 26 80 34 40 83 39 BS 100% 54 Gy 9 15 { 10 10 19 10 12 12 42 FMO # ap ertures 166 140 194 192 170 157 168 185 160 157 T o compare the n um b er of ap ertures that w ere generated to the n um b er of ap ertures that w ere actually used, i.e., that receiv ed a p ositiv e in tensit y in the optimal solution giv en a particular n um b er of ap ertures, w e p erformed the follo wing linear regression c haracterizing this relationship for all 10 cases: Num b er of ap ertures used = 0 + 1 Num b er of ap ertures generated : T able 5{5 sho ws the estimates for the regression co ecien ts, as w ell as the go o dness-of-t. The p ositiv e in tercept 0 represen ts that a particular n um b er of ap ertures (apparen tly on the order of 2-4 p er b eam) should alw a ys b e used for treatmen t. Out of all additional ap ertures generated, a little o v er 1 out of 2, on
PAGE 137
126 a v erage, turns out to b e useful in the long run, while the other ones are ev en tually replaced b y one or more b etter ap ertures. T able 5{5: Regression results of the relationship b et w een generated and used ap ertures. Case 1 2 3 4 5 6 7 8 9 10 0 22. 0 27. 2 31. 4 18. 8 17. 5 21. 7 19. 3 25. 9 24. 1 18. 9 1 0. 55 0.60 0. 61 0. 54 0. 58 0.53 0. 59 0. 57 0.56 0. 61 R 2 (%) 99. 0 98. 0 97. 9 99. 2 99. 3 98. 3 99. 7 98. 7 98. 5 99. 1 F or one of the cases (case 5), Figures 5{3 5{6 illustrate the p erformance of our metho d in more detail. Figure 5{3 (a) sho ws the b eha vior of the ob jectiv e function v alue as a function of the n um b er of ap ertures generated, as w ell as the n um b er of ap ertures used (i.e., the n um b er of ap ertures with nonzero in tensit y in the corresp onding solution). Figure 5{3 (b) sho ws the n um b er of ap ertures used v ersus the n um b er of ap ertures generated. Figures 5{4 5{6 sho w ho w the n um b er of ap ertures used aects the lev el of satisfaction of the external criteria. F rom our results, w e conclude that all external planning criteria except for the criteria for PTV2 can b e satised with, on a v erage, 65 ap ertures (whic h corresp onds to 9-10 ap ertures p er b eam). In clinical practice, ph ysicians ma y decide to treat using these more ecien t plans at the cost of lo w er co v erage of PTV2. This trade-o b et w een treatmen t plan eciency and PTV2-co v erage needs to b e made on a case-b y-case basis based on the exten t of the disease and the individual lev el of clinical risk asso ciated with underdosing PTV2. T o obtain a go o d co v erage and a smaller lev el of underdosing of PTV2, an a v erage of 101 ap ertures is needed. In fact, in 7 out of the 10 cases few er than 100 ap ertures (appro ximately 14 p er b eam) suce. Ev en when the external planning criteria that are m uc h stricter than the published R TOG criteria are tak en in to accoun t, the a v erage required n um b er of ap ertures is only 131. Note that although it ma y tak e more ap ertures than this to ac hiev e the same ob jectiv e function v alue as the
PAGE 138
127 (a) (b) Figure 5{3: (a) Ob jectiv e function v alue as a function of n um b er of ap ertures generated and used, and comparison with the optimal v alue to the b eamlet FMO problem. (b) Num b er of ap ertures used as a function of n um b er of ap ertures generated. (a) (b) Figure 5{4: Co v erage of (a) PTV1 and (b) PTV2 as a function of n um b er of ap ertures used. (a) (b) Figure 5{5: Sparing of saliv a glands according to (a) D VH criterion (v olume > 26 Gy) and (b) mean dose as a function of n um b er of ap ertures used. (a) (b) Figure 5{6: Sparing of (a) spinal cord (v olume > 45 Gy) and brainstem (v olume > 54 Gy) and (b) unsp ecied tissue as a function of n um b er of ap ertures used.
PAGE 139
128 optimal v alue obtained using the b eamlet mo del for FMO, the sub optimal solutions obtained using this n um b er of ap ertures are of clinically acceptable qualit y . Most curren t commercial treatmen t planning systems t ypically require 20-30 ap ertures p er b eam, or a total of 140-210 ap ertures in order to ac hiev e a clinically acceptable treatmen t plan, whic h corresp onds to ab out 50-100% more ap ertures as our solutions. This is illustrated b y the last line in T able 5{4 , whic h sho ws that in all cases but one a traditional t w o-stage approac h yields signican tly more ap ertures than the ap erture mo dulation approac h. Our metho dology th us suggests that patien ts can b e treated in m uc h less time than in curren t practice. In addition, and p erhaps ev en more imp ortan tly , our approac h will pro vide the ph ysicians with a to ol to mak e a sound trade-o b et w een the n um b er of ap ertures (i.e., treatmen t time) and the lev el of satisfaction of the treatmen t plan criteria. These ndings ha v e b een published in Romeijn et al [ 58 ]. 5.6 Concluding Remarks and F uture Researc h In this c hapter, w e ha v e dev elop ed a new approac h to the design of radiation therap y treatmen t plans. Our approac h incorp orates deliv ery issues in to the design phase that ha v e to date usually b e handled in a second optimization phase. Using 10 problem instances obtained from clinical data, w e ha v e sho wn that patien t treatmen t times can b e reduced while main taining treatmen t plan qualit y b y in tegrating these t w o optimization phases. In our approac h, w e are curren tly using a simple tec hnique for nding an initial solution for the column generation pro cedure. F uture researc h will fo cus on dev eloping new heuristic tec hniques for nding b etter initial solutions, in order to reduce the n um b er of steps required b y the column generation pro cedure. Our mo del curren tly assumes that the b eam orien tations are xed in adv ance. Another direction of future researc h can b e incorp orating the b eam orien tation selection problem to o in to this optimization problem.
PAGE 140
APPENDIX A TONGUE AND GR OO VE CONSTRAINT VIOLA TION Consider the shap e matrix, 0 2 2 1 Here, I 1 ; 2 = 2, and I 2 ; 2 = 2. The T&G constrain ts dictate that b eamlet (2,2) should b e exp osed only when b eamlet (1,2) is exp osed. The shap e ro ws whic h corresp ond to the ap ertures whic h exp ose b eamlet (1,2) corresp ond to no des (1,0,3) and (1,1,3) in the graph, and the corresp onding ro ws for b eamlet (2,2) are no des (2,0,3) and (2,1,3). The shap e graph after applying the T&G correction and re-dening the arc sets is giv en in Fig. A{1 . F or clarit y , w e sho w only arcs corresp onding to elemen ts in column 2. In Fig. A{1 , dashed lines denote the path only whic h can b e used to select no des (2,0,3) or (2,1,3) to exp ose b eamlet (2,2) (Note that all the other arcs connecting these no des to the no des in la y er 1 will b e rejected b y T&G mo dication prop osed in section 4.2.3 ). After prepro cessing, ruence constrain ts will imp ose an upp er b ound of max(0,2) = 0 on arcs emanating from no de (1,0,3) and max(2) = 2 on arcs emanating from no de (1,1,3). A p ossible ap erture w ould no w b e the one corresp onding to path f source(1,1,3)-(2,0,2)-sink g with in tensit y 2. After this extraction, the matrix left w ould b e 0 0 0 1 whic h, after up dating upp er b ounds after extraction do es not ha v e a path connecting from source to sink. W e cannot terminate at this p oin t as not all 129
PAGE 141
130 1r ,r 0r ,r 1r Sourcer 2r ,r 0r ,r 1r 2r ,r 0r ,r 2r 2r ,r 0r ,r 3r 2r ,r 1r ,r 2r 2r ,r 1r ,r 3r 2r ,r 2r ,r 3r 1r ,r 2r ,r 3r 1r ,r 1r ,r 3r 1r ,r 1r ,r 2r 1r ,r 0r ,r 3r 1r ,r 0r ,r 2r Sinkr Figure A{1: An example of T&G constrain t. the ruence is deliv ered but w e cannot ev en extract a shap e as there is no path from source to sink. Th us, this instance b ecomes infeasible. So, to enforce T&G corrections in the example ab o v e, w e need to imp ose a b ound of 1 (2-1) on the arcs emanating from no des (1,0,2) and (1,0,3) but not terminating on no des (2,0,2) or (2,0,3).
PAGE 142
REFERENCES [1] Ah uja R K, Magnan ti T, and Orlin J. 1993. Network Flows: The ory, A lgorithms, and Applic ations . Pren tice-Hall, Englew o o d Clis, New Jersey [2] Ah uja R K. 2003. V ery large-scale neigh b orho o d (VLSN) searc h algorithms for reducing dose deliv ery times. Pr esentation at The First International Workshop on Optimization in R adiation Ther apy F ort Lauderdale [3] Ah uja R K, and Hamac her H W. 2005. Linear time net w ork ro w algorithm to minimize b eam-on-time for unconstrained m ultileaf collimator problems in cancer radiation therap y . Networks 45 (1):36-41 [4] Alb er M, Meedt G, Nusslin F and Reem tsen R. 2002. On the degeneracy of the IMR T optimization problem. Me dic al Physics 29 (11):2584{2589 [5] Andreatta G, and Mason F. 1985. k-Eccen tricit y and absolute k-cen trum of a tree. Eur op e an Journal of Op er ational R ese ar ch 19 :114{117 [6] Baatar D, Hamac her H W, Ehrgott M, and W o eginger G. 2004. Decomp osition of in teger matrices and m ultileaf collimator sequencing. R ep ort in Wirtschaftsmathematik h ttp://kluedo.ub.uni-kl.de/v olltexte/2004/1786 [7] Bates J, and Constable R. 1985. Pro ofs as programs. A CM T r ansactions on Pr o gr amming L anguages and Systems 7 :113{136 [8] Bazaraa M, Sherali H, and Shett y C. 1993. Nonline ar Pr o gr amming: The ory and A lgorithms John Wiley & Sons, 2nd edition, New Y ork [9] Bednarz G, Mic halski D, Houser C, Huq M S, Xiao Y, Anne P R et al. 2002. The use of mixed-in teger programming for in v erse treatmen t planning with pre-dened eld segmen ts. Physics in Me dicine and Biolo gy 47 (13):2235{2245 [10] Ben tley J. 1986. Pr o gr amming Pe arls Addison-W esley , Reading, MA [11] Boland N, Hamac her H W, and Lenzen F. 2004. Minimizing b eam-on time in cancer radiation treatmen t using m ultileaf collimators. Networks 43 (4):226240 [12] Bortfeld T, Burk elbac h J, Bo esec k e R, and Sc hlegel W. 1990. Metho ds of image reconstruction from pro jections applied to conformal radiotherap y . Physics in Me dicine and Biolo gy 25 (4):435{443 131
PAGE 143
132 [13] Bortfeld T, and Sc hlegel W. 1993. Optimization of b eam orien tations in radiation-therap y some theoretical considerations. Physics in Me dicine and Biolo gy 38 (2):291{304 [14] Bortfeld T R, Kahler D L, W aldron T J, and Bo y er A L. 1994. X-ra y eld comp ensation with m ultileaf collimators. International Journal of R adiation Onc olo gy Biolo gy Physics 28 :723{730 [15] Bortfeld T, Stein J, and Preiser K. 1997. Clinically relev an t in tensit y mo dulation optimization using ph ysical criteria. In Pr o c e e dings of the XIIth ICCR Salt Lak e Cit y , Utah, Ma y 27-30:1{4 [16] Bortfeld T, Kufer K-H, Monz M, Sc herrer A, Thiek e C, and T rinkhaus H. 2003. In tensit y mo dulated radiotherap y: a large scale m ulti-criteria programming problem. Berichte des F r aunhofer ITWM 43 F raunhofer Insitut T ec hnound Wirtsc haftsmathematik, Kaiserslautern, German y [17] Bro c kstein B and Masters G. 2003. He ad and Ne ck Canc er . Klu w er Academic Publishers, Boston, Massac h usetts [18] Burk ard R E. 2002. Op en researc h question section. Ob erwolfach c onfer enc e on Combinatorial optimization Ob erw olfac h, German y , No v em b er 24-29 [19] Cancer F acts & Figures Rep ort. 2002. American Cancer So ciet y , A tlan ta, GA [20] Carol M P , Nash R V, Campb ell R C, Hub er R and Sternic k E. 1997. The dev elopmen t of a clinically in tuitiv e approac h to in v erse treatmen t planning: partial v olume prescription and area cost function, Pr o c e e dings of the XII th ICCR ,Salt Lak e Cit y , Utah:317{319 [21] Cattaneo G M, Fiorino C, Lom bardi P , and Calandrino R. 1997. Optimizing the mo v emen t of a single absorb er for 1d non-uniform dose deliv ery b y (fast) sim ulated annealing. Physics in Me dicine and Biolo gy 42 :107{121 [22] Chao K S C, Blanco A, and Dempsey J F. 2003. A conceptual mo del in tegrating spatial information to assess target v olume co v erage for IMR T treatmen t planning. International Journal of R adiation Onc olo gy Biolo gy Physics 56 (5):1438{1449 [23] Cho P S, Lee S, Marks I I I R J, Redstone J A, and Oh S. 1997. Comparison of algorithms for in tensit y mo dulated b eam optimization: pro jections on to con v ex sets and sim ulated annealing. In Pr o c e e dings of the XIIth ICCR Salt Lak e Cit y , Utah, Ma y 27-30:310{312 [24] Cho P S, Lee S, Marks I I I R J, Oh S, Sutlief S G, and Phillips M H. 1998. Optimization of in tensit y mo dulated b eams with v olume constrain ts using t w o metho ds: cost function minimization and pro jections on to con v ex sets. Me dic al Physics 25 (4):435{443
PAGE 144
133 [25] Cho P S, and Marks R J. 2000. Hardw are-sensitiv e optimization for in tensit y mo dulated radiotherap y . Physics in Me dicine and Biolo gy 45 (2):429{440 [26] Cro oks S M, Pugac hev A, King C, and Xing L. 2002. Examination of the eect of increasing the n um b er of radiation b eams on a radiation treatmen t plan. Physics in Me dicine and Biolo gy 47 (19):3485{3501 [27] Das I J, and Cheng C W. 1999. T reatmen t plan ev aluation using dose-v olume histogram (D VH) and spatial dose-v olume histogram (zD VH). International Journal of R adiation Onc olo gy Biolo gy Physics 43 (5):1143{1150 [28] Das S K, and Marks L B. 1997. Selection of coplanar or noncoplanar b eams using three-dimensional optimization based on maxim um b eam separation and minimized non target irradiation. International Journal of R adiation Onc olo gy Biolo gy Physics 38 (3):643{655 [29] Das S, Cullip T, T racton G, Chang S, Marks L, Ansc her M et al. 2003. Beam orien tation selection for in tensit y-mo dulated radiation therap y based on target equiv alen t uniform dose maximization. International Journal of R adiation Onc olo gy Biolo gy Physics 55 (1):215{224 [30] Deasy J O. 1997. Multiple lo cal minima in radiotherap y optimization problems with dose-v olume constrain ts. Me dic al Physics 24 (7):1157{1161 [31] Dempsey J F, Deasy J O, Lomax A, Wiesmey er M, Bosc h W, and Lo w D A. 2001. T reatmen t plan review to ols incorp orating spatial dose-v olume information. International Journal of R adiation Onc olo gy Biolo gy Physics 51 Supplemen t 1:125. Abstract of oral presen tation at the 43rd Ann ual ASTR O Meeting [32] Hamac her H and Kufer K-H. 2002. In v erse radiation therap y planning a m ultiple ob jectiv e optimization approac h. Discr ete Applie d Mathematics 118 :145{161 [33] Hiriart-Urrut y J-B and Lemarec hal C. 1996. Convex A nalysis and Minimization A lgorithms I: F undamentals . Springer-V erlag, Berlin, German y [34] Holmes T and Mac kie T R. 1994. A ltered bac kpro jection dose calculation in v erse treatmen t planning. Me dic al Physics 21 :321{333 [35] Kamath S, Sahni S, Li J, P alta J, and Rank a S. 2003. Leaf sequencing algorithms for segmen ted m ultileaf collimation. Physics in Me dicine and Biolo gy 48 :307{324 [36] Kamath S, Sahni S, P alta J, Rank a S, and Li J. 2004. Optimal Leaf Sequencing with Elimination of T ongue-and-Gro o v e. Physics in Me dicine and Biolo gy 49 :N7{N19
PAGE 145
134 [37] Langer M, Bro wn R, Urie M, Leong J, Strac her M and Shapiro J. 1990. Large-scale optimization of b eam w eigh ts under dose-v olume restrictions. International Journal of R adiation Onc olo gy Biolo gy Physics 18 (4):887{893 [38] Langer M, Morrill S, Bro wn R, Lee O and Lane R. 1996. A comparison of mixed in teger programming and fast sim ulated annealing for optimizing b eam w eigh ts in radiation therap y . Me dic al Physics 23 (6):957{964 [39] Langer M, Thai V, and P apiez L. 2001. Impro v ed leaf sequencing reduces segmen ts or monitor units needed to deliv er IMR T using m ultileaf collimators. Me dic al Physics 28 :2450{2458 [40] Lee E K, F o x T, and Cro c k er I. 2000. Optimization of radiosurgery treatmen t planning via mixed in teger programming Me dic al Physics 27 :995{1004 [41] Lee E K, F o x T and Cro c k er I. 2003. In teger programming in in tensit ymo dulated radiation therap y treatmen t planning. A nnals of Op er ations R ese ar ch 119 :165{181 [42] Ma L, Bo y er A, Xing L, and Ma C-M. 1998. An optimized leaf-setting algorithm for b eam in tensit y mo dulation using dynamic m ultileaf collimators. Physics in Me dicine and Biolo gy 43 :1629-1643 [43] Mageras G S and Mohan R. 1993. Application of fast sim ulated annealing to optimization of conformal radiation treatmen ts. Me dic al Physics 20 :639{647 [44] Murt y K. 1976. Line ar Pr o gr amming and Combinatorial Optimization . John Wiley & Sons, New Y ork [45] Oelfk e U and Bortfeld T. 2001. In v erse planning for photon and proton b eams. Me dic al Dosimetry 26 (2):113{124 [46] Ogryczak W and Ruszczy nski A. 2002. Dual sto c hastic dominance and related mean-risk mo dels. SIAM Journal on Optimization 13 :60{78 [47] Ogryczak W and T amir A. 2003. Minimizing the sum of the k largest functions in linear time. Information Pr o c essing L etters 85 :117{122 [48] Ogryczak W and Za w adzki M. 2002. Conditional mean: a parametric solution concept for lo cation problems. A nnals of Op er ations R ese ar ch 110 :167{181 [49] Otto K, and Clark B G. 2002. Enhancemen t of IMR T deliv ery through MLC rotation. Physics in Me dicine and Biolo gy 47 (22):3997{4017 [50] Pugac hev A B, Bo y er A L, and Xing L. 2000. Beam orien tation optimization in in tensit y-mo dulated radiation treatmen t planning. Me dic al Physics 27 (6):1238{1245
PAGE 146
135 [51] Pugac hev A, and Lei X. 2001. Pseudo b eam's-ey e-view as applied to b eam orien tation selection in in tensit y-mo dulated radiation therap y . International Journal of R adiation Onc olo gy Biolo gy Physics 51 (5):1361{1370 [52] Pugac hev A, Li J G, Bo y er A L, Hanco c k S L , Le Q T, Donaldson S S et al. 2001. Role of b eam orien tation optimization in in tensit y-mo dulated radiation therap y . International Journal of R adiation Onc olo gy Biolo gy Physics 50 (2):551{560 [53] Pugac hev A, and Xing L. 2001. Computer-assisted selection of coplanar b eam orien tations in in tensit y-mo dulated radiation therap y . Physics in Me dicine and Biolo gy 46 (9):2467{2476 [54] Pugac hev A, and Xing L. 2002. Incorp orating prior kno wledge in to b eam orien tation optimization in IMR T. International Journal of R adiation Onc olo gy Biolo gy Physics 54 (5):1565{1574 [55] Que W. 1999. Comparison of algorithms for m ultileaf collimator eld segmentation. Me dic al Physics 26 :2390{2396 [56] Ro c k afellar R T and Ury asev S. 2000. Optimization of conditional v alue-atrisk. The Journal of R isk 2 (3):21{41 [57] Ro c k afellar R T and Ury asev S. 2002. Conditional v alue-at-risk for general loss distributions. Journal of Banking & Financ e 26 (7):1443{1471 [58] Romeijn H, Ah uja R, Dempsey J, and Kumar A. 2005. A column generation approac h to radiation therap y treatmen t planning using ap erture mo dulation. SIAM Journal of Optimization 15 (3):838{862 [59] Romeijn H, Ah uja R, Dempsey J, and Kumar A. 2005. A new linear programming approac h to radiation therap y treatmen t planning problems. T o app ear in Op er ations R ese ar ch [60] Romeijn, H E, Ah uja R K, Dempsey J F, Kumar A, and Li J G. 2003. A no v el linear programming approac h to ruence map optimization for in tensit y mo dulated radiation therap y treatmen t planning. Physics in Me dicine and Biolo gy 48 (21):3521{3542 [61] Romeijn H, Dempsey J, and Li J. 2004. A unifying framew ork for m ulticriteria ruence map optimization mo dels. Physics in Me dicine and Biolo gy 49 (10):1991{2013 [62] Rosen I, Lam K S, Lane R G, Langer M, and Morrill S M. 1995. Comparison of sim ulated annealing algorithms for conformal therap y treatmen t planning. International Journal of R adiation Onc olo gy Biolo gy Physics 33 (5):1091{1099 [63] Rosen I, Lane R G, Morrill S M, and Belli J A. 1990. T reatmen t plan optimization using linear programming. Me dic al Physics 18 (5):141{152
PAGE 147
136 [64] Ro wb ottom C G, Nutting C M, and W ebb S. 2001. Beam-orien tation optimization of in tensit y-mo dulated radiotherap y: clinical application to parotid gland tumors. R adiother apy and Onc olo gy 59 (2):169{177 [65] R TOG. 2000. H-0022: A phase I/I I study of conformal and in tensit y mo dulated irradiation for oropharyngeal cancer h ttp://www.rtog.org/mem b ers/proto cols/h0022/index.h tml [66] R TOG. 2002. H-0225: A phase I I study of in tensit y mo dulated radiation therap y (IMR T) +/c hemotherap y for nasopharyngeal cancer h ttp://www.rtog.org/mem b ers/proto cols/0225/p df le.h tml [67] Seco J, Ev ans P M, and W ebb S. 2002. An optimization algorithm that incorp orates IMR T deliv ery constrain ts. Physics in Me dicine and Biolo gy 47 (6):899{915 [68] Shepard D M, Earl M A, Li X A, Naqvi S, and Y u C. 2002. Direct ap erture optimization: a turnk ey solution for step-and-sho ot IMR T. Me dic al Physics 29 (6):1007{1018 [69] Shepard D M, F erris M C, Oliv era G H and Mac kie T R. 1999. Optimizing the deliv ery of radiation therap y to cancer patien ts. SIAM R eview 41 (4):721{ 744 [70] Shepard D M, Oliv era G H, Rec kw erdt P J and Mac kie T R. 2000. Iterativ e approac hes to dose optimization in tomotherap y . Physics in Me dicine and Biolo gy 45 (1):69{90 [71] Sieb ers J V, Lauterbac h M, Keall P J, and Mohan R 2002 Incorp orating m ulti-leaf collimator leaf sequencing in to iterativ e IMR T optimization Me dic al Physics 29 (6):952{959 [72] Sio c hi R A C. 1999. Minimizing static in tensit y mo dulation deliv ery time using an in tensit y solid paradigm. International Journal of R adiation Onc olo gy Biolo gy Physics 42 :671{680 [73] Slater P . 1978. Cen ters to cen troids in a graph. Journal of Gr aph The ory 2 :209{222 [74] Spirou S V and Ch ui C-S. 1998. A gradien t in v erse planning algorithm with dose-v olume constrain ts. Me dic al Physics 25 (3):321{333 [75] T amir A. 2001. The k-cen trum m ulti-facilit y lo cation problem. Discr ete Applie d Mathematics 109 :293{307 [76] Tsien C, Eisbruc h A, McShan D, Kessler M, Marsh R and F raass B. 2003. In tensit y-mo dulated radiation therap y (IMR T) for lo cally adv anced paranasal sin us tumors: Incorp orating clinical decisions in the optimization pro cess. International Journal of R adiation Onc olo gy Biolo gy Physics 55 (3):776{784
PAGE 148
137 [77] Ury asev S. 2000. Conditional V alue-at-Risk: Optimization algorithms and applications. Financial Engine ering News 14 :1{5 [78] W ebb S. 1989. Optimization of conformal radiotherap y dose distributions b y sim ulated annealing. Physics in Me dicine and Biolo gy 34 :1349{1370 [79] W ebb S. 1992. Optimization b y sim ulated annealing of three-dimensional conformal treatmen t planning for radiation elds dened b y m ultiple collimator. I I. Inclusion of t w o dimensional mo dulation of X-ra y in tensit y . Physics in Me dicine and Biolo gy 37 :1689{1704 [80] W ebb S. 1999. Conformal in tensit y-mo dulated radiotherap y (IMR T) deliv ered b y rob otic linac-testing IMR T to the limit? Physics in Me dicine and Biolo gy 44 (7):1639{1654 [81] W ebb S. 2001. Intensity-Mo dulate d R adiation Ther apy Institute of Ph ysics Publishing, Philadelphia [82] W ebb S and Lomax T. 2001. There is no IMR T? Physics in Me dicine and Biolo gy 46 (12):L7{L8 [83] W u Q and Mohan R. 2000. Algorithms and functionalit y of an in tensit y mo dulated radiotherap y optimization problem. Me dic al Physics 27 (4):701{ 711 [84] W u X, and Zh u Y. 2001. An optimization metho d for imp ortance factors and b eam w eigh ts based on genetic algorithms for radiotherap y treatmen t planning. Physics in Me dicine and Biolo gy 46 (4):1085{1099 [85] Xia P , and V erhey L J. Multileaf collimator leaf sequencing algorithm for in tensit y mo dulated b eams with m ultiple static segmen ts. Me dic al Physics 25 :1424{1434 [86] Xing L. 2003. Radiation Therap y|101. Pr esentation at The First International Workshop on Optimization in R adiation Ther apy F ort Lauderdale [87] Xing L and Chen G T Y. 1996. Iterativ e metho ds for in v erse treatmen t planning. Physics in Me dicine and Biolo gy 41 :2107{2123
PAGE 149
BIOGRAPHICAL SKETCH Arvind Kumar is a nativ e of India and got his early education at Gorakhpur, India. He obtained his bac helor's degree in c hemical engineering from Indian Institute of T ec hnology , Madras, India. Since 2001, he has b een pursuing his PhD degree in the Departmen t of Industrial and Systems Engineering at Univ ersit y of Florida. 138