Tilburg University
Improving the quality, efficiency and robustness of radiation therapy planning and delivery through mathematical optimization
Balvert, Marleen
Publication date:
2017
Document Version
Publisher's PDF, also known as Version of record Link to publication in Tilburg University Research Portal
Citation for published version (APA):
Balvert, M. (2017). Improving the quality, efficiency and robustness of radiation therapy planning and delivery through mathematical optimization. CentER, Center for Economic Research.
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal
Take down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
Improving the quality, efficiency and robustness
of radiation therapy planning and delivery
through mathematical optimization
Proefschrift
ter verkrijging van de graad van doctor aan Tilburg University op gezag van de rector magnificus, prof.dr. E.H.L. Aarts, in het openbaar te verdedigen ten overstaan van een door het college voor promoties aangewezen commissie in de aula van de Universiteit op
vrijdag 31 maart 2017 om 14.00 uur door
Marleen Balvert
Promotor: prof. dr. ir. Dick den Hertog
Copromotores: dr. ir. Aswin Hoffmann
dr. David Craft
Overige leden: prof. dr. Ben Heijmen
prof. dr. ir. Edwin van Dam dr. Juan Vera Lizcano dr. ir. Marleen Keijzer
Improving the quality, efficiency and robustness of radiation therapy planning and delivery through mathematical optimization
Copyright c 2017 Marleen Balvert
Acknowledgements
After three and a half years, my dissertation is finished, and my time as a PhD student is over. I have loved being a PhD student: over the past years, I have learned an incredible amount and met great people, some of whom I would like to thank.
I cannot but start with my supervisors. Dick, you are a great example to me. We share the philosophy that one should contribute to the benefit of society, which brought us both to treatment plan optimization. You put your PhD students high on your priority list. Your positive way of giving feedback is remarkable: even when a text I had written was completely red with your remarks, you said “Great work!” I said: “But...?” and you said “No but. Let’s go through the comments.” This positive attitude makes you very approachable and creates a safe environment for us to make mistakes that we can learn from. You are a kind, compassionate and reliable person, someone I feel safe to confide in and like to share my successes with. Without the slightest exaggeration I can say that I could not have wished for a better, kinder and more devoted promotor.
Aswin, I admire your devotion to research and that you speak two languages: that of a clinical physicist and that of a mathematician. You introduced me to the world of a clinical physicist, which was completely alien to me. You taught me how to properly structure a text and revise my own work critically, from the paper as a whole to the smallest detail. But your most valuable lesson to me lies in your comments on my introductory chapter, which pushed me to think beyond one-paper projects and look at the bigger picture. Thank you for your valuable guidance throughout the past six years.
David, when we met at a conference in Portugal, we clicked right away. At the conference you talked about “The VMAT Problem”, and shortly after that, I joined you in your quest. I very much enjoyed our discussions, scribbling away on the whiteboard. You also showed me some nice places and events in Boston and Cambridge. I admire how you balance enjoying life and doing what is best for the environment and for others. Thank you for the great time in Boston and the wise lessons you taught me.
mem-bers of my dissertation committee: Ben Heijmen, Edwin van Dam, Juan Vera Lizcano and Marleen Keijzer. Thank you for reading this thesis and providing comments, and for taking the effort to come to Tilburg for the defense. I would also like to thank my co-authors, who have made all the work over the past years possible. Stefan van Hooff, Frank Verhaegen, Daniela Trani and Patrick Granton from MAASTRO Clinic, Steven Petit and Sebastiaan Breedveld from Erasmus MC and Cees Roos from TU Delft, thank you for the pleasant cooperation. Sincere acknowledgements to Bram Gorissen for his supervision during my bachelor’s and master’s thesis: I have learned a lot from your critical view on even the smallest details. Many thanks to the people from the radiation oncology group at MGH, who made my stay not only educative, but very pleasant as well. Special thanks to Tom Madden for being the drive behind Friday afternoon drinks and other social events, and for proofreading these acknowledgements.
Over the past few years, I not only entered the world of scientific research, but also the world of education. This opportunity was given to me by Henk, Herbert, Susan and especially Marieke (Quant). Marieke, you gave me my first opportunity to give a lecture in front of a big group of students, which was pretty scary at the time. In the years that followed you were my mentor in lecturing, making exams and coordinating a course. Thank you for everything you have taught me about teaching. Many of my colleagues have made the past years very enjoyable. I would like to thank Anja, Lenie and Heidi for their administrative support and pleasant chats. Ahmadreza, Alaa, Amparo, Bas, Bas, Frans, Hettie, Jan, Liz, Maria, Mario, Olga, Shobeir, Sybren and Viktoriya, thank you for the lunches, tea breaks, Sinterklaas evenings, GSS events, game evenings, drinks at Kandinsky and, most importantly, cake afternoons. I would like to thank four people in particular. Nick, your cheerful-ness is infectious! It’s always nice talking with you, and not a single topic is awkward. I’ll never forget “Who is the EOR mole?” which we organized together with Anja. Krzysztof, although we are opposites in many ways, you’re a great friend. You’re always up for a “toppertje”, and you don’t hesitate to ask critical questions which others may avoid. Trevor, you’re an amazing office mate: we talked about research, life or just plain nonsense. When you were in Singapore, the office felt empty, and I’ll miss having you as an office mate. Marieke, we did so much together that people kept mixing us up, though I still don’t understand why. Yes, we are both Dutch, blond, the same height, enjoy running and baking, started our PhDs in 2013, love teaching, had boyfriends with the same first name, got our BKO, had tea breaks together... but that’s all! I really enjoyed having you as a colleague, and missed our chats from the moment you left Tilburg.
iii is more to life than work. Therefore, I’d like to thank my friends from econometrics for being such a great distraction from research through our many gatherings. The same goes for my friends from secondary school: Beau, Maike, Manouk, Sanne and Sanne. We’ve done a lot together, from going crazy during carnaval up to helping each other getting through breakups with boyfriends. Ladies, thank you for being such sweet friends for the past fourteen years.
There are four people I’d like to pay special attention to, as they have been great supporters in many ways, and have closely followed my progress as a PhD student. Beau, I am incredibly lucky to have you as a friend. We share our ups and downs: I am always anxious to tell you good news as soon as possible, and you’re always there to listen and support me when I need it. Thank you for all these years of friendship. Now we are both leaving Tilburg, a city that became a home to us both, but what I’ll miss most is living close to you.
I would like to take the opportunity to express my deepest gratitude to my sister. Eline, you’ve been more than my dearest friend all your life, and our connection grew even stronger when we outgrew childhood. We understand each other like no-one else does (and that includes our sense of humor). You are my home, and you’re the one person whom I can tell everything. I know I can always rely on you. Thank you for being the lovely sister you are.
Finally, and most importantly, I would like to thank my mom and dad. You have always and in every possible way supported Eline and me. You’ve encouraged our ambition to get the best out of ourselves, both at school and in our personal development, but never pushed us. You gave us a home where every issue, no matter how difficult, could be discussed, and where we learned to support each other and share in each other’s happiness and sadness. Without a safe home like that, I could never have become who I am now, and I will always be very grateful to you.
Marleen Balvert
Contents
1 Introduction 1
1.1 Radiation treatment modalities 2
1.1.1 High-dose rate brachytherapy 3
1.1.2 Intensity-modulated radiation therapy 3
1.1.3 Volumetric modulated arc therapy 6
1.1.4 Small animal radiotherapy 6
1.2 Inverse treatment plan optimization problem 6
1.2.1 Multi-criteria optimization 8
1.2.2 Objective functions 9
1.2.3 Optimization via fluence maps 12
1.3 Uncertainties in radiotherapy 13
1.3.1 Sources of uncertainty 13
1.3.2 Uncertainties in treatment plan optimization 14
1.4 Summary of contributions 17
1.5 Discussion and conclusions 21
1.6 Future perspectives 23
1.7 Overview of research papers 25
1.8 Disclosure 26
2 A framework for inverse planning of beam-on times for 3D small
animal radiotherapy using interactive multi-objective optimization 27
2.1 Introduction 28
2.1.1 Motivation 28
2.1.2 Inverse treatment planning in radiotherapy 29
2.1.3 Pre-clinical treatment planning 30
2.2 Methods and materials 32
2.2.1 Inverse planning optimization model 32
2.2.2 Treatment planning 34
2.2.3 Treatment planning and delivery times 34
2.3 Results 36
2.3.1 Dose distributions 36
2.3.2 Trade-off 40
2.3.3 Treatment planning and delivery times 40
2.4 Discussion 41
2.5 Conclusion 44
2.A Example automated optimization process 45
3 Dwell time modulation restrictions do not necessarily improve
treatment plan quality for prostate HDR brachytherapy 51
3.1 Introduction 52
3.1.1 Clinical procedure and workflow 53
3.1.2 High-dose subvolumes 53
3.1.3 Uncertainty in dwell locations 54
3.1.4 Dwell time modulation restrictions 55
3.1.5 Aim of the paper 55
3.2 Methods and materials 56
3.2.1 Dwell time optimization models 56
3.2.2 Modulation restrictions 56
3.2.3 Patient data and simulations 57
3.2.4 Plan quality indicators 58
3.3 Results 59
3.3.1 Interpretation of the figures 59
3.3.2 Relative dwell time difference restricted in (LD) 60
3.3.3 General results 62
3.4 Discussion 62
3.5 Conclusion 64
4 Robust optimization of dose-volume metrics for prostate HDR
brachytherapy incorporating target volume delineation
uncertain-ties 67
4.1 Introduction 68
4.1.1 Review of methods accounting for delineation uncertainties 69
4.1.2 Aim and contribution of the paper 71
4.2 Treatment plan optimization model 72
4.2.1 Dose prescription and plan evaluation 72
4.2.2 Nominal treatment planning model 72
4.2.3 Accounting for delineation uncertainties of the PTV 74
vii
4.3 Computational experiments 79
4.3.1 Patient and source data 79
4.3.2 Experiment setup 79
4.3.3 Numerical results 82
4.4 Discussion 87
4.5 Conclusion 89
4.A Full nominal treatment plan optimization model 90
5 Fast approximate delivery of fluence maps for IMRT and VMAT 93
5.1 Introduction 93
5.2 Materials and Methods 95
5.2.1 Treatment planning model 95
5.2.2 Four step solution approach 97
5.2.3 Step 1: treatment plan optimization for a single arc segment 98
5.2.4 Step 2: delivery time allocation 100
5.2.5 Step 3: fluence map optimization for the full arc 101
5.2.6 Step 4: alternate improvement of dose rate schedule and leaf
trajectories 102
5.2.7 Test cases 102
5.3 Results 103
5.3.1 The single map case 103
5.3.2 The full arc case 107
5.4 Discussion 111
5.5 Conclusion 116
5.A The need for non-unidirectional leaf trajectories 116
5.B Gradient and Hessian 119
6 Robust dose-painting accounting for uncertainties in dose-response
parameters and patient positioning 121
6.1 Introduction 122
6.2 Methods 124
6.2.1 Dose-painting optimization based on TCP 124
6.2.2 Accounting for uncertainties in dose painting 125
7 A universal and structured way to derive dual optimization
prob-lem formulations 141
7.1 Introduction 141
7.2 Preliminaries 143
7.3 Fenchel’s dual problem of a convex optimization problem 145
7.4 Recipe for setting up Fenchel’s dual problem 146
7.5 Some applications 148
7.5.1 Linear optimization 148
7.5.2 Conic optimization 149
7.5.3 Quadratically constrained quadratic optimization 149
7.5.4 Second order cone optimization 150
7.5.5 Geometric optimization 152
7.5.6 `p-norm optimization 153
7.5.7 Radiotherapy treatment planning 155
7.6 Concluding remarks 161
7.A Some examples of conjugate functions 162
7.A.1 Linear function 162
7.A.2 Log-sum-exp function 162
7.A.3 Arbitrary norm 163
7.A.4 Quadratic function 164
7.A.5 Conjugate of 1pxp, x ≥ 0, p > 1 164
7.A.6 Conjugate of px1p, x > 0, p > 0 165
7.B Some identities for conjugate functions 166
7.B.1 Conjugate of f (x) = maxifi(x) 166
7.B.2 Deriving conjugates via the adjoint I 166
7.B.3 Deriving conjugates via the adjoint II 167
7.B.4 Linear substitution rule 168
7.B.5 Composite function 168
7.B.6 Conjugates with ‘branches’ 169
7.B.7 Support function of a cone 170
7.B.8 Discontinuity of the closure of a perspective function 170
7.C Lagrange’s duality theorem for convex optimization 171
7.D Fenchel duality 171
7.E Derivation of Fenchel’s dual problem 172
7.F Tables of conjugates and support functions 174
A Supplementary data to Chapter 3 179
A.1 Relative dwell time difference restricted 180
ix
A.1.2 (LDV) model 187
A.2 Absolute dwell time difference restricted 194
A.2.1 (LD) model 194
A.2.2 (LDV) model 201
A.3 Quadratic dwell time difference restricted 208
List of Figures
1.1 A schematic view of the treatment for prostate HDR brachytherapy (a) and the template that is used for this treatment (b). From Gorissen
(2014), with permission. 4
1.2 Illustration of a linear accelerator suitable for both IMRT and volu-metric modulated arc therapy (VMAT, see Section 1.1.3), treating a patient with a brain tumour. The gantry, the collimator and the treat-ment couch can be rotated. From K¨ufer et al. (2005), with permission. 5 1.3 While the leaves block out some of the bixels as in the aperture in (a),
the beam is switched on to deliver the fluence map in (b). Consecu-tively applying these apertures, either in a step-and-shoot fashion or through dynamic delivery, gives the total fluence map shown in (c). 5 1.4 An example of a cumulative dose-volume histogram for the target
vol-ume (solid line) and an OAR (dashed line), where the prescription dose to the target volume is 60 Gy. The minimum dose to the hottest 50% of the OAR, or the median OAR dose, D50%(OAR), and the volume
receiving at least 100% of the prescribed dose, V100%(TV), are indicated. 11
1.5 Illustration of an MLC (top), the total resulting fluence map (middle) and its contribution to the total dose distribution in a transversal slice of the patient’s scan (bottom). The fluence map shapes the dose distribution. (Image courtesy of RaySearch AB, Stockholm, Sweden.) 12 2.1 Illustration of Pareto optimality with two objectives. In figure (a),
solutions 3 and 4 are not Pareto optimal, solutions 1 and 2 are Pareto optimal, and therefore lie on the Pareto frontier. Figure (b) shows the complete Pareto frontier. Note that the region below the Pareto
frontier represents infeasible solutions. 30
2.2 Schematic overview of the different phases in precision small animal
2.3 Beam configurations of the four treatment cases. Beam directions are indicated with arrows. Dynamic arc beams are illustrated by drawing the initial and final beam positions, connected with curved arrows. Beam index numbers correspond to the beam-on time data in table 2.2. 36 2.4 Overview of the dose distributions resulting from automated and
man-ual treatment plan optimization and their differences. Window and level settings for visualization of both the CT and overlaid dose dis-tributions are equal for all images. PTVs are delineated in red, the right kidney in yellow, the left kidney in blue, the GI in white, and
the spine in green. 38
2.5 DVHs obtained using manual (gray) and automated (black) beam-on time optimization. The vertical dashed line indicates the prescription
dose. 39
2.6 A trade-off between target coverage (PTV V95%) as a percentage of
the total volume and sparing of (a) the GI (D1%) for case 1, (b) the
GI (D5%) for case 2, (c) the left kidney (D5%) for case 3 and (d) the
spine (D5%) for case 4, measured in Gy. 40
3.1 Graphs for patient 1 showing: (a) penalty and (b) D90%(PTV)
gen-erated by (LD) extended with DTMR-R. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard
deviation from the mean. 61
3.2 Graphs for patient 1 showing: high-dose volume performance indica-tors (a) DHI and (b) DVHc generated by (LD) extended with DTMR-R. The solid lines in graph (a) represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 62 4.1 Two delineations of a prostatic target volume based on a transversal
ultrasound imaging scan. Both delineations yield a different set of calculation points residing in the structure: the blue points reside in both structures, the yellow and red points only in the yellow and red
delineation, respectively. 69
4.2 Example of sets of calculation points that are likely to receive a high dose ( ˜In
C) or a low dose ( ˜I f
C). 77
4.3 Vectorized figure of minimal (light gray), delineated (medium gray)
xii
4.4 Transrectal ultrasound scan with the PTV delineated in green. The red dwell positions are included in the original dataset, the green ones
are added to obtain the extended set of dwell positions. 81
4.5 Histogram of (a) V100%(CTV) as a percentage of the total CTV
vol-ume, (b) D90%(CTV) as a percentage of the prescribed dose and (c)
V200%(CTV) as a percentage of the total CTV volume for patient 1,
measured over all scenarios for the nominal (black), margin model, (white), robust model (gray), margin model with extra dwell posi-tions (white striped) and robust model with extra dwell posiposi-tions (gray striped). The vertical dashed line indicates the desired minimum level for V100%(CTV) and D90%(CTV) and the desired maximum level for
V200%(CTV). 83
4.6 Histogram of (a) V100%(CTV) as a percentage of the total CTV
vol-ume, (b) D90%(CTV) as a percentage of the prescribed dose and (c)
V200%(CTV) as a percentage of the total CTV volume for patient 2,
measured over all scenarios for the nominal (black), margin model, (white), robust model (gray), margin model with extra dwell posi-tions (white striped) and robust model with extra dwell posiposi-tions (gray striped). The vertical dashed line indicates the desired minimum level for V100%(CTV) and D90%(CTV) and the desired maximum level for
V200%(CTV). 84
4.7 Histogram of (a) V100%(CTV) as a percentage of the total CTV
vol-ume, (b) D90%(CTV) as a percentage of the prescribed dose and (c)
V200%(CTV) as a percentage of the total CTV volume for patient 3,
measured over all scenarios for the nominal (black), margin model, (white), robust model (gray), margin model with extra dwell posi-tions (white striped) and robust model with extra dwell posiposi-tions (gray striped). The vertical dashed line indicates the desired minimum level for V100%(CTV) and D90%(CTV) and the desired maximum level for
V200%(CTV). 85
5.1 Relationship between bixel indices and leaf positions. This illustration is for B = 5 bixels. The leaves are drawn as gray rectangles and they are both shorter than they would actually be: the left leaf would extend leftward and the right leaf would extend rightward, as shown by the dots, for the opening shown. We assume that all leaves are physically long enough to cover the entire row of bixels if needed. 96 5.2 An example of a fluence map where optimization is more difficult for
5.3 Example of a random leaf trajectory (solid lines) for three arc segments with 12 time steps allocated to each of them. The arc segments are indicated by the dashed lines. The arc segments are subdivided into three, one and three subsegments, respectively, (dotted lines). The leaf trajectories of these seven subsegments are a left-right sweep, an open-out trajectory, a close-in trajectory, a left-right sweep, a random trajectory, a left-right sweep and a right-left sweep, respectively. 102 5.4 Prostate case: Solution to the fluence map matching problem for
var-ious delivery times. The y-axis of the graph shows the objective func-tion that is minimized,P
ij(fij−gij)2, which is labelled as ssdif (sum of
squared differences). Fluence map units are MU. The dose rate plots show the dose rate versus time for each delivery. The solid line is for optimized dose rates and the dashed line is for fixed dose rate of 10 MU/s (the maximum allowed value for the optimized case). Note that for very short delivery times the maximum dose rate is used always
due to the need to get in sufficient fluence. 104
5.5 Head and neck case example (for description, see Figure 5.4 legend). 105 5.6 Two sample leaf trajectories for a difficult (high SPG) head and neck
fluence map row. The left plot shows the fluence row that the optimizer is trying to match along with the resulting fluence from the optimal 5 second and 11 second solutions. The right figures show these solutions and also indicate with the heat map the dose rate chosen at each time step (compare with dose rates as plotted in Figure 5.5). Here, white corresponds to the maximum dose rate of 10 MU/s. The left leaf is
depicted in light gray, the right leaf in dark gray. 106
5.7 Trade-off between delivery time and ssdif for the prostate case. Results obtained at step 3 with optimized dose rates (solid) and dose rates equal to 10 MU/s (dashed) are presented in gray. Improvements from these results obtained in step 4 are indicated in black. The lower bound based on individually optimized arc segments is shown by the
dotted line. 108
5.8 Fluence maps for arc segment 11 of the prostate case, optimized using a total arc delivery time (Ttotal) of 60s, 90s, 120s, 150s and 180s. Fluence
maps obtained with individual arc segment optimization are optimized with delivery times equal to the time allotted to this arc segment in the full arc treatment plans (6s, 9s, 15s, 18s and 18s, respectively).
xiv
5.9 DVH plots for the PTV (black) (PTV 56 from the CORT data set) and rectum (gray) for the plans obtained with various total times (Ttotal,
in seconds) and the sliding window approach. 110
5.10 The leaf trajectories for row 8 of the prostate case using variable (a) and constant maximum dose rates (b) for treatment durations from 1 (left-most plot) up to 3 minutes (right-most plot) with 30 second increments. The horizontal lines demarcate the time allotted to each
arc segment. 110
5.11 Trade-off between delivery time and ssdif for the head & neck case. Results obtained at step 3 with optimized dose rates (solid) and dose rates equal to 10 MU/s (dashed) are presented in gray. Improvements from these results obtained in step 4 are indicated in black. The lower bound based on individually optimized arc segments is shown by the
dotted line. 111
5.12 Fluence maps for arc segment 15 of the head and neck case, optimized using a total arc delivery time (Ttotal) of 150s, 330s, 510s and 690s.
Fluence maps obtained with individual arc segment optimization are optimized with delivery times equal to the time allotted to this arc seg-ment in the full arc treatseg-ment plans (4s, 6s, 16s and 20s, respectively).
The corresponding ssdif is indicated below each map. 112
5.13 The leaf trajectories for row 22 of the head and neck case using variable (a) and constant maximum dose rates (b) for treatment durations from 75 (left-most plot) up to 345 seconds (right-most plot) with 90 second
increments. 113
5.14 An example of bidirectional (a) trajectories for the left and right leaf that can be transformed into unidirectional (b) trajectories without
changing the fluence map. 117
5.15 The total fluence delivered to a bixel is equal to the area enclosed by
the leaf trajectories. 118
5.16 The only optimal VMAT treatment plan for a leaf pair with desired fluence map [5 5.5 1 5.5 10 5.5 1 5.5 5] when time is limited to T = 8. 118 5.17 Figure (a) shows a non-unidirectional VMAT plan that gives fluence
map [10 17 2 17 10 0 0 0 0]. With the given dose rates, it is not possible to find unidirectional leaf trajectories that yield the fluence
6.1 Sampled scenarios of radiosensitivity parameters α1 and α2 for the
phantom case showing (a) their individual probability distributions with the sampled points and (b) the corresponding subvolume control
probability clouds curves representing ± 1 SD for α. 129
6.2 (a) Sampled scenarios of radiosensitivity parameters for the CTV ex-cluding GTV, for the resistant and sensitive part of the GTV and (b) the corresponding control probability clouds representing ± 1 SD for α .130 6.3 Dose distributions for the phantom of the phantom for different
radio-sensitivty distributions. Positional uncertainty is ignored here. (a) No uncertainty in radio-sensitivity is considered; (b) Gaussian uncertainty distributions for α are incorporated; (c) uniform uncertainty distribu-tions for α are incorporated. The outer white contour represents the resistant subvolume 1 and the inner contour represents the sensitive
subvolume 2. 132
6.4 Dose profiles for the phantom case along the (a) right-left, (b) anterior-posterior and (c) caudal-cranial direction through the center of the target, without considering positioning uncertainty. The light and dark shaded areas show the resistant subvolumes 1 and the sensitive subvolume 2 respectively. Treatment plans are generated while ig-noring uncertainty in radio-sensitivity (green), and while accounting for Gaussian (blue) and uniform (red) uncertainty in radio-sensitivity. The ideal dose distributions obtained with the stylistic model are in-dicated with the dashed lines. DVHs for the resistant (solid) and
sensitive (dashed) subvolume are shown in (d). 133
6.5 Dose profiles for the phantom case along the (a) right-left, (b) anterior-posterior and (c) caudal-cranial direction through the center of the tar-get, while accounting for positioning uncertainty. The light and dark shaded areas denote subvolumes 1 and 2 respectively. Treatment plans are generated while ignoring radiosensitivity uncertainty (green) and accounting for uncertainty in radiosensitivity using a Gaussian distri-bution on α (blue). All plans account for positioning uncertainties. For comparison the corresponding dose distributions that ignored
xvi
6.6 Dose distributions for the NSCLC patient obtained when (a) ignoring uncertainty; (b) accounting for uncertainty in radiosensitivity only; (c) accounting for positional uncertainty only; (d) accounting for both radio-sensitivity and positioning uncertainties. The contours shown represent the CTVprim (red), the GTVprim (black), the resistant part
of the GTVprim (white) and PTVnodes (pink). 136
A.1 Objective function value for all patients as function of the maximum relative dwell time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the
mean. 180
A.2 DHI for all patients as function of the maximum relative dwell time difference γ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 181
A.3 PTV DVHc for all patients. 182
A.4 D90%(PTV) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 183 A.5 V100%(PTV) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 183 A.6 V150%(PTV) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 184 A.7 V200%(PTV) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 184 A.8 D10%(rectum) for all patients as function of the maximum relative
A.9 D2cc(rectum) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 185 A.10 D10%(urethra) for all patients as function of the maximum relative
dwell time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean.186 A.11 D0.1cc(urethra) for all patients as function of the maximum relative
dwell time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean.186 A.12 V100%(PTV) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 187 A.13 DHI for all patients as function of the maximum relative dwell time
difference γ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 188
A.14 PTV DVHc for all patients. 189
A.15 D90%(PTV) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 190 A.16 V150%(PTV) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 190 A.17 V200%(PTV) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 191 A.18 D10%(rectum) for all patients as function of the maximum relative
xviii
A.19 D2cc(rectum) for all patients as function of the maximum relative dwell
time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 192 A.20 D10%(urethra) for all patients as function of the maximum relative
dwell time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean.192 A.21 D0.1cc(urethra) for all patients as function of the maximum relative
dwell time difference γ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean.193 A.22 Objective function value for all patients as function of the absolute
dwell time difference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean.194 A.23 DHI for all patients as function of the absolute dwell time difference
θ. The solid lines represent the minimum, mean and maximum values,
and the dotted line is the pre-plan value. The grey area denotes values
at most one standard deviation from the mean. 195
A.24 PTV DVHc for all patients. 196
A.25 D90%(PTV) for all patients as function of the absolute dwell time
dif-ference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
values at most one standard deviation from the mean. 197
A.26 V100%(PTV) for all patients as function of the absolute dwell time
dif-ference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
values at most one standard deviation from the mean. 197
A.27 V150%(PTV) for all patients as function of the absolute dwell time
dif-ference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
values at most one standard deviation from the mean. 198
A.28 V200%(PTV) for all patients as function of the absolute dwell time
dif-ference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
A.29 D10%(rectum) for all patients as function of the absolute dwell time difference θ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 199 A.30 D2cc(rectum) for all patients as function of the absolute dwell time
difference θ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 199 A.31 D10%(urethra) for all patients as function of the absolute dwell time
difference θ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 200 A.32 D0.1cc(urethra) for all patients as function of the absolute dwell time
difference θ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 200 A.33 V100%(PTV) for all patients as function of the absolute dwell time
dif-ference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
values at most one standard deviation from the mean. 201
A.34 DHI for all patients as function of the absolute dwell time difference
θ. The solid lines represent the minimum, mean and maximum values,
and the dotted line is the pre-plan value. The grey area denotes values
at most one standard deviation from the mean. 202
A.35 PTV DVHc for all patients. 203
A.36 D90%(PTV) for all patients as function of the absolute dwell time
dif-ference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
values at most one standard deviation from the mean. 204
A.37 V150%(PTV) for all patients as function of the absolute dwell time
dif-ference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
values at most one standard deviation from the mean. 204
A.38 V200%(PTV) for all patients as function of the absolute dwell time
dif-ference θ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
xx
A.39 D10%(rectum) for all patients as function of the absolute dwell time difference θ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 205 A.40 D2cc(rectum) for all patients as function of the absolute dwell time
difference θ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 206 A.41 D10%(urethra) for all patients as function of the absolute dwell time
difference θ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 206 A.42 D0.1cc(urethra) for all patients as function of the absolute dwell time
difference θ. The solid lines represent the minimum, mean and maxi-mum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 207 A.43 Objective function value for all patients as function of the maximum
sum of squared differences ρ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the
mean. 208
A.44 DHI for all patients as function of the maximum sum of squared differ-ences ρ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes
values at most one standard deviation from the mean. 209
A.45 PTV DVHc for all patients. 210
A.46 D90%(PTV) for all patients as function of the maximum sum of squared
differences ρ. The solid lines represent the minimum, mean and max-imum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 211 A.47 V100%(PTV) for all patients as function of the maximum sum of squared
differences ρ. The solid lines represent the minimum, mean and max-imum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 211 A.48 V150%(PTV) for all patients as function of the maximum sum of squared
A.49 V200%(PTV) for all patients as function of the maximum sum of squared differences ρ. The solid lines represent the minimum, mean and max-imum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean. 212 A.50 D10%(rectum) for all patients as function of the maximum sum of
squared differences ρ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean.213 A.51 D2cc(rectum) for all patients as function of the maximum sum of
squared differences ρ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean.213 A.52 D10%(urethra) for all patients as function of the maximum sum of
squared differences ρ. The solid lines represent the minimum, mean and maximum values, and the dotted line is the pre-plan value. The grey area denotes values at most one standard deviation from the mean.214 A.53 D0.1cc(urethra) for all patients as function of the maximum sum of
List of Tables
2.1 Dose-volume metrics using manual (Man) and automatic (Aut) beam-on time optimizatibeam-on. Dose values are given in Gy, and V95% is given
as a percentage of the total PTV volume. σD/Dmean is the standard
deviation of the dose to the voxels in the PTV divided by the mean
PTV dose. L=left, R=right. 37
2.2 Beam-on times, time spent on plan optimization and time required to deliver treatment plans for the manual (Man) and automated (Aut) optimization approaches. Beam-on times are provided per beam, and the total beam-on time per treatment plan. The time required to deliver treatment plans is split into the total beam-on time for the plan, and time to reposition the gantry and switch the X-ray tube on and off, referred to as non-beam-on time. All times are denoted as
min:s. 42
2.3 Dose-volume metrics using automatic beam-on time optimization where the optimization weights for the OARs are all equal to 0.25, and the PTV optimization weights vary from 0.1 to 0.9 with a step size of 0.1. Dose values are given in Gy, and V95% is given as a percentage of the
total PTV volume. 45
2.4 Dose-volume metrics using automatic beam-on time optimization where the optimization weight for the left kidney is 0.5, the weights for the three remaining OARs are all equal to 1/6, and the PTV optimization weights vary from 0.15 to 0.25 with a step size of 0.01. Dose values are given in Gy, and V95% is given as a percentage of the total PTV
volume. 47
2.6 Dose-volume metrics using automatic beam-on time optimization where the optimization weight for the left kidney is 0.5, the wieghts for the three remaining OARs are all equal to 1/6, and the PTV optimization weights vary from 0.05 to 0.14 with a step size of 0.01. Dose values are given in Gy, and V95% is given as a percentage of the total PTV
volume. 49
3.1 Dose-volume criteria, based on the protocol by Hoskin et al. (2007). 59 4.1 Dose-volume criteria based on the protocol by Hoskin et al. (2007). 72
4.2 Problem sizes for (LDV ) and (RC). 76
4.3 Tissue structure volumes (cc). 80
4.4 Number of dose calculation points. 80
4.5 Summary of DVH statistics for the OARs for each scenario, obtained with five different planning models. All values are given in Gy. dps =
dwell positions. 86
4.6 Comparison of treatment plans generated with (LDV ) and plans
gen-erated with a relaxation of (LDV ) (Rel.). 86
4.7 Solution times (s) for all solution approaches and three patients, dps
= dwell positions. 87
4.8 Model parameters used in the numerical experiments. 91
6.1 Evaluation of treatment plans obtained while ignoring positional un-certainty for different distributions of unun-certainty in radio-sensitivity. The rows denote the different dose distributions, the column denote the different uncertainty scenarios used to estimate the corresponding
TCP values. 132
6.2 DVH statistics and expected TCP obtained for the different dose dis-tributions. The expected TCP is calculated using Gaussian uncer-tainty for α in the nominal and worst case positioning scenarios. 135 7.1 Examples of convex conjugate functions available from the literature.
In lines 6 and 7 k.k denotes any norm, and k.k∗ the dual of this norm. In lines 8 - 10, the numbers p and q are related according to 1
p+
1
q = 1.
In lines 13 and 14 P† denotes the generalized inverse of matrix P . In lines 15 and 16, Sn
+ and Sn++ denote the set of positive semidefinite
and positive definite n × n matrices, respectively. BV = Boyd and Vandenberghe (2004), R = Rockafellar (1970), GH = Gorissen and
xxiv
7.2 Rules for calculating the conjugate function. The functions f (x), fi(x)
and h(x) in this table are real valued functions, f (x) and fi(x) convex
and h(x) concave. In line 5 we assume ∩m
i=1ri.dom fi 6= ∅; if fi(x)
is linear for some i then the corresponding ri.dom fi can be replaced
by dom fi in this condition. We call this the sum-rule for conjugate functions. In line 6, x1 : xmis a partition of x. We will refer to line 7 as the linear substitution rule. In lines 9 and 10, f♦ denotes the adjoint
of f : R++ → R, which is defined as f♦(x) := xf (1/x). In line 11,
S denotes the unit simplex. BV = Boyd and Vandenberghe (2004),
R = Rockafellar (1970), HU = Hiriart-Urruty (2006), G = Gushchin
(2008), BT = Ben-Tal et al. (1991). 176
7.3 Support functions. In line 3, K is a pointed cone and if K is nonlin-ear, then it is assumed that S contains a strictly feasible point. The functions fi and the sets Si in lines 5-7 are closed convex. In line 5
we assume ∩m
i=1ri.dom fi 6= ∅; if fi(x) is linear for some i then the
corresponding ri.dom fi can be replaced by dom fi in this condition.
Similarly, in line 6 we assume ∩mi=1ri Si 6= ∅, and if Si is polyhedral
then the corresponding ri Si can be replaced by Si. BV = Boyd and
Vandenberghe (2004), BT = Ben-Tal et al. (2015), R = Rockafellar
Chapter 1
OARs Organs at risk
CT Computed Tomography
MRI Magnetic Resonance Imaging
US Ultrasound
TV Target volume
EBRT External beam radiotherapy
BT brachytherapy
HDR High-dose rate
IMRT Intensity-modulated radiotherapy VMAT Volumetric modulated arc therapy MLC Multi-leaf collimator
MCO Multi-criteria optimization
DVH Dose-volume histogram
TCP Tumor control probability
NTCP Normal tissue complication probability EUD Equivalent uniform dose
DTMR Dwell time modulation restriction
GTV Gross target colume
CTV Clinical target volume
Chapter 2
PTV Planning target volume
CBCT Cone-beam Computed Tomography
MC Monte Carlo
HU Hounsfield Units
ICRU International Commission on Radiation Units & Measurements GI Gastrointestinal tract
xxvi
Chapter 3
LD Linear dose (model)
LDV Linear dose-volume (model)
DHI Dose homogeneity index
IPSA Inverse Planning by Simulated Annealing HIPO Hybrid Inverse Plann Optimization TRUS Transrectal Ultrasound
Chapter 5
Linac Linear accelerator
MU Monitor Units
SPG Sum of positive gradients DAO Direct aperture optimization
Chapter 6
NSCLC Non-small cell lung cancer
SD Standard deviation
SUV Standardized uptake volume
MLD Mean lung dose
VCP Volume control probability
Chapter 7
CHAPTER 1
Introduction
Cancer is one of the most prevalent diseases in the world, with 14.1 million patients newly diagnosed in 2012, causing 8.2 million deaths in that same year (Ferlay et al., 2013). These numbers are expected to increase to 22.2 million cancer patients and 13 million deaths by 2030 (Bray et al., 2012). There is thus an urgent need for the (further) development of cancer treatments.
a starting point and the task is to find the design parameters such that the dose objectives are met as well as possible.
This dissertation covers various topics in treatment plan optimization. The first chapter gives an introduction to radiation therapy and some of the basics in treatment plan optimization. This chapter starts with a brief description of a typical workflow followed in radiation therapy and a discussion of the four radiation treatment modal-ities considered in this work (Section 1.1). Section 1.2 provides an overview of the modeling approaches commonly used in the treatment plan optimization literature and practice. The quality of a treatment plan is prone to various sources of geometric and dosimetric uncertainties, which can have a large impact on treatment outcome as radiotherapy is a high-precision treatment. This is the topic of Section 1.3, fol-lowed by a short introduction to three common methods to incorporate uncertainties in the planning process: margins, stochastic programming and robust optimization. The chapter is concluded by an overview of the contributions in this dissertation in Section 1.4, concluding remarks in Section 1.5 and future perspectives in Section 1.6. Section 1.7 comprises the research articles on which this work is based and Section 1.8 contains a disclosure.
1.1
Radiation treatment modalities
Prior to radiation delivery, a treatment plan needs to be designed and optimized. First a computed tomography (CT), magnetic resonance imaging (MRI) or ultra-sound (US) scan is made to depict the patient’s anatomy. Subsequently a radiation oncologist delineates the anatomical structures of interest (i.e. the target volume (TV) and the OARs) on the scan. These delineations form the geometrical basis for treatment planning. Treatment plan optimization is an interdisciplinary pro-cess where the radiation technologists, radiation oncologist and medical physicist are involved to design a treatment plan for an individual patient. This is often a trial-and-error process, which may take a substantial amount of time. Once the treatment plan has been approved, actual radiation delivery can start. Dependent on the de-livery method, the type of cancer and the patient, the dede-livery can be repeated in a number of fractions at intervals ranging from a couple of hours up to a few days.
Radiation treatment modalities 3 modality depends on characteristics of the tumour, such as its volume (BT is only suitable for relatively small tumours) and accessibility. In some cases, BT is used to give a boost dose in addition to the EBRT. There are various types of EBRT (e.g. stereotactic radiotherapy, intensity modulated radiotherapy) and brachytherapy (e.g. low-, high- and pulsed dose rate). For an extensive overview, see Hevezi (2003). The treatment modalities that are considered in this dissertation are further discussed in the remainder of this section.
1.1.1
High-dose rate brachytherapy
High-dose rate (HDR) BT obtained its name from the high dose rate radioactive source that is temporarily placed inside or close to the tumour. This treatment is amongst others used for prostate cancer, for which it has been shown to result in a highly conformal dose distribution and high efficacy (Yamada et al., 2012).
The patient is placed in dorsal lithotomy position and a template containing a grid of around 40 holes is positioned in front of the patient’s perineum (Figure 1.1). Out of these 40 holes, 15 up to 20 (depending on the tumour volume in mainly the lateral directions) are selected for inserting catheters into the prostate, usually perpendicular to the template. A small radioactive source can now be moved through each catheter consecutively by a remote afterloader device. The source stops at predetermined dwell locations for a predetermined amount of time (dwell time) to deliver a certain amount of dose to the tissue. Delivery is generally repeated several times at time intervals of at least six hours, while the catheters remain in situ. Note that the treatment procedure varies among hospitals.
As longer dwell times result in a higher dose deposited to the nearby tissue, the number of catheters, their locations and the dwell times together shape the dose distribution and thus constitute the degrees of freedom in the treatment planning process. Treatment planning for HDR BT for prostate cancer is the focus of Chapters 3 and 4.
1.1.2
Intensity-modulated radiation therapy
Prostate
Catheters Template
US probe
(a) (b)
Figure 1.1 – A schematic view of the treatment for prostate HDR brachyther-apy (a) and the template that is used for this treatment (b). From Gorissen (2014), with permission.
IMRT is the topic of Chapters 5 and 6.
A radiation beam can be virtually divided into a 2D grid of beam elements per-pendicular to the beam direction. “Intensity-modulated” refers to the characteristic feature of IMRT to modulate the intensity (or rather fluence) over the bixels. This is achieved by using a multi-leaf collimator (MLC, see Figure 1.3a for a schematic figure), which is a device that is positioned in the linear accelerator head. The MLC is made up of thin tungsten “leaves” lined up on either side of the beam. Each leaf pair can (partially) block a row of bixels (Figure 1.3a) so that radiation can only pass through the aperture formed by the leaves (Figure 1.3b).
Consecutively applying different apertures results in a modulation of the fluence over the bixels. This can be done either in a step-and-shoot fashion, where the beam is switched off while the leaves move to their positions for the next aperture, or using dynamic delivery, where the leaves move continuously while leaving the beam switched on. The total fluence delivered at a bixel is determined by the total amount of time a bixel is exposed multiplied with the fluence rate. An overview of the fluence at each bixel is presented in a fluence map, which is often denoted as an intensity map (Figure 1.3c) or an intensity matrix.
Radiation treatment modalities 5
Figure 1.2 – Illustration of a linear accelerator suitable for both IMRT and volumetric modulated arc therapy (VMAT, see Section 1.1.3), treating a pa-tient with a brain tumour. The gantry, the collimator and the treatment couch can be rotated. From K¨ufer et al. (2005), with permission.
Additionally, one can optimize the couch position and the collimator angle, though these degrees of freedom have not yet been included in treatment plan optimization methods. For a more elaborate description of IMRT and the underlying planning problem, see Webb (2003) and Bortfeld (2006).
(a) (b) (c) 0 MU 5 MU 10 MU 15 MU 20 MU
1.1.3
Volumetric modulated arc therapy
Volumetric modulated arc therapy (VMAT, first proposed by Yu 1995) is very similar to IMRT, with the difference that fluence is delivered while the gantry makes a continuous full-arc rotation around the patient instead of delivering fluence from a discrete set of preselected beam angles. The MLC leaves move continuously across the beam opening as well. This introduces the gantry rotation speed as an additional degree of freedom in the plan optimization procedure.
An advantage of VMAT over IMRT may be that the same plan quality can be achieved with a shorter delivery time, as the beam does not need to be switched off to change the device settings (Unkelbach et al., 2015). An even lower delivery time can be achieved if a (minor) reduction in plan quality is accepted. This trade-off between delivery time and plan quality is inherent to both IMRT and VMAT, which is the topic of Chapter 5.
1.1.4
Small animal radiotherapy
In order to get a better understanding of the biological effects of radiotherapy, pre-clinical experiments are conducted. These experiments are often performed on small animals, e.g. rats and mice. The linear accelerators that are used in the clinic are unsuitable for high-precision pre-clinical radiotherapy: because the animals are much smaller than human patients, one cannot use the same (high) energy level. This requires a different type of energy, for which micro-irradiators embedded in dedicated systems for pre-clinical irradiation have recently been developed (Verhaegen et al., 2011; Granton et al., 2014). These are low-voltage X-ray devices where kilovoltage X-ray beams are delivered from a discrete set of preselected directions. A range of different circular and rectangular beam sizes can be used by manual insertion of fixed aperture collimators into a collimator holding assembly. Treatment plan optimization for such systems is the topic of Chapter 2.
1.2
Inverse treatment plan optimization problem
Inverse treatment plan optimization problem 7 form the basis for treatment planning by the medical physicist, who will consider a plan optimal whenever the plan meets the protocol requirements. A mathematician will only consider a treatment plan optimal when it gives the global optimum to a pre-defined objective function, preferably with a guarantee on its optimality. Inverse treatment planning methods aim at finding treatment plans that are mathematically optimal with respect to some objective function, though due to lengthy calculations one often needs to accept a near optimal solution. The key issue in the development of inverse treatment planning models is to define an objective function which best represents the clinical goals, for which the dose prescription protocol can provide guidelines.
To obtain the dose distribution that would result from a given treatment plan, the volumes of the anatomical structures delineated on the images of the patient’s anatomy are discretized into small volume elements, called voxels. Alternatively, one could define so-called dose calculation points over the structure volumes. As these two methods are analogous, for notational convenience only voxels are mentioned in the following. For given bixels (for EBRT) or dwell positions (for HDR BT), the dose per unit fluence (for EBRT) or per second (for EBRT or HDR BT) to a voxel is pre-calculated. This dose rate is stored in a matrix ˙D, where ˙Dij is the dose rate
from bixel (for EBRT) or dwell position (for HDR BT) j to voxel i. The dose rate can either be given in Gy per unit fluence if we wish to optimize the total fluence delivered (for EBRT), or in Gy per second if we wish to optimize the beam-on times (for small animal radiotherapy) or dwell times (for HDR BT). When optimizing the fluence, the decision variable xj denotes the fluence at bixel j, and the dose to a voxel i is calculated by ˙dTi x, where ˙diis the ith column of ˙D containing the dose rates from
all bixels to voxel i. When optimizing beam-on times or dwell times, the decision variable tj denotes the beam-on time of beam j or the dwell time at dwell position j, and the total dose delivered to voxel i can be calculated as ˙dTi t. The doses in all voxels combined give the dose distribution, which is used for both treatment plan optimization and evaluation. In the remainder of this section, we will use the decision variable t to denote radiation times (for HDR BT and small animal radiotherapy), which can be replaced by the variable x if bixel intensities are to be optimized (for EBRT).
results in a multi-criteria optimization problem. Below, a discussion of multi-criteria optimization in radiation treatment planning is given. Furthermore, three types of objectives in treatment plan optimization models are discussed. Finally, optimization via fluence maps is explained.
1.2.1
Multi-criteria optimization
As noted before, treatment plan optimization is inherently a multi-criteria optimiza-tion (MCO) problem. We are interested in Pareto-optimal soluoptimiza-tions, i.e., soluoptimiza-tions where none of the objectives can be improved without deteriorating at least one of the other objectives. The two most common approaches to MCO problems that identify a Pareto-optimal solution are the weighted-sum method and the epsilon-constraint method (Miettinen, 1999). The weighted sum method optimizes a weighted sum of the conflicting objectives, where different weights result in different trade-offs between the objectives. In the epsilon-constraint method, one of the objectives is optimized while placing a constraint on the values of the remaining objectives.
Many Pareto-optimal solutions are not usable, as they do not yield a desirable balance between the various objectives (for example, not treating the patient is Pareto optimal, as it fully spares the OARs). One needs to find the Pareto optimal plan that yields the most preferred balance between the conflicting objectives. However, when applying the weighted sum method to radiotherapy treatment planning, it is difficult to a priori determine the weights that yield the most preferred trade-off between the objectives as there is no clinical interpretation for their values. A similar issue arises for the epsilon-constraint method: since the trade-off between the objectives is different for each individual patient, it is not known a priori to what extent the dose to the various OARs can be minimized, and hence, one cannot set the constraint bounds a priori. In order to find the most preferred trade-off, it is thus necessary to identify multiple Pareto optimal solutions and pick the most preferred one a posteriori. This has led to the development of Pareto front generation for treatment planning (K¨ufer et al., 2005; Craft et al., 2006; Hoffmann et al., 2006; Bokrantz, 2013). Optimizing a weighted sum of the conflicting objectives for various (relative) weights allows the user to investigate different trade-offs and select the most preferable solution a posteriori. Also with the epsilon constraint method one can construct a set of Pareto-optimal solutions by varying the bounds on the constrained objectives.
Inverse treatment plan optimization problem 9 methods generate a single Pareto optimal solution. The advantage is that the planner does not need to navigate the Pareto front in search for the most preferred trade-off, but is presented a solution that meets the predetermined dose requirements in order of priority. This results in more uniformity amongst plans optimized by different planners. It may however be a challenging task to define a prioritized list of objectives without knowing their interdependence in advance.
1.2.2
Objective functions
There are several ways to translate the dose prescriptions into an optimization model and planning objectives. The simplest approach is to prescribe a minimum (mean) dose to (voxels of) the tumour, and a maximum (mean) dose to (voxels of) each of the structures. Alternatively, one can use dose-volume histogram (DVH) metrics, providing insight in the fraction of a structure that receives a particular dose. These metrics are typically used in practice for dose prescriptions and treatment plan eval-uation, as they are strongly correlated to treatment outcome (mostly side effects). However, the true goals are of course to maximize the probability of tumour control (TCP) and to limit the probability of normal tissue complications (NTCP), for which we could use radiobiology-based (i.e., dose-response based) objectives. As there is still quite some ambiguity regarding the true values of the parameters that shape these functions (Nahum and Kutcher, 2007; Marks et al., 2010), their use remains limited.
Penalty-based dose objectives One can prescribe a minimum dose to the tumour voxels combined with a desired maximum on the dose to both tumour and OAR voxels. Penalty-based objective functions aim to minimize a penalty based on the difference between these prescriptions and the planned dose. If the planned dose to a single voxel lies between the lower and the upper bound (where the lower bound is zero for OARs), the penalty to this voxel equals zero, while it is a non-decreasing function of the over- or underdose to that voxel. The penalty pi to a voxel i can thus
be written as: pi(t) = 0 if Li ≤ ˙d T i t ≤ Ui fi(Li− ˙d T i t) if ˙d T i t < Li gi( ˙dTi t − Ui) if ˙d T i t > Ui,
where fi and gi are increasing functions, and Li and Ui are the lower- and upper
bounds on the dose to voxel i, respectively.
various structures. However, there is no clinical interpretation for the values of the weights.
The main advantage of penalty-based objectives is their simplicity. In HDR BT, the functions fi(t) and gi(t) are often chosen to be linear, resulting in a piece-wise
linear function for pi(t) (Lessard and Pouliot, 2001) and thus a linear programming
model. However, the penalty values often show limited correlation to DVH metrics (Holm, 2011; Gorissen et al., 2013), which complicates the task of determining the penalty weight values. For IMRT it is more common to minimize the sum of squared differences between the desired and the planned dose (Bortfeld, 2006), which has the advantage that the problem size does not increase with the number of voxels (Lahanas et al., 2003). Due to the simplicity of the penalty-based models, even large instances of several thousands of bixel weights can be solved fast.
DVH metrics In current clinical practice, treatment plans are evaluated based on DVH metrics since they show a strong correlation with treatment outcome (Marks et al., 2010). A dose-volume histogram presents the fraction of a structure’s volume as a function of the administered dose, and DVH metrics can be derived from this.
Dx%(S) and Dycc(S) denote the minimum dose received by the hottest x% and y cc
of the volume of structure S, respectively. Alternatively, Vw%(S) and VzGy(S) denote
the fraction of the volume of structure S receiving at least w% of the prescription dose and z Gy, respectively. An example of a cumulative dose-volume histogram is given in Figure 1.4. The solid and dashed lines indicate how much dose is delivered to the volume fraction of the TV and an OAR, respectively, and the prescribed tumour dose is 60 Gy. In this example, we observe that D50%(OAR) = 20 Gy and
V100%(TV) = 90%.
Dose prescriptions are often formulated as requirements on DVH metrics. For example, one can require V100%(TV) ≥ 90%, or D2cc(Rectum) ≤ 7 Gy. Such
require-ments can be included in the treatment planning model as constraints. Alternatively, one can optimize DVH-metrics, e.g. maximize V100%(TV) under the constraint that
D2cc(Rectum) ≤ 7 Gy.
DVH-based models directly optimize the metrics that are used as quality measures for a treatment plan. As a result, they are more comprehensible than penalty-based models and yield treatment plans that better meet the DVH criteria (Langer et al., 1990; Spirou and Chui, 1998; Alterovitz et al., 2006; Holm, 2011; Gorissen et al., 2013). However, translating these objectives into a mathematical optimization program inevitably results in a mixed integer program (Siauw et al., 2011; Gorissen et al., 2013), which requires longer solution times than the penalty-based optimization model.
Inverse treatment plan optimization problem 11 Dose (Gy) 0 20 40 60 80 Relativ e v olume (%) 25 50 75 100 D50%(OAR) V100%(TV)
Figure 1.4 – An example of a cumulative dose-volume histogram for the target volume (solid line) and an OAR (dashed line), where the prescription dose to the target volume is 60 Gy. The minimum dose to the hottest 50% of the OAR, or the median OAR dose, D50%(OAR), and the volume receiving at least 100% of the prescribed dose, V100%(TV), are indicated.
DVH prescriptions are surrogates for the true goal of cancer treatment, which is to control the tumour while minimizing the risk of side-effects. We could directly include TCP and NTCP into the optimization model. Recall that these can both be described as functions of the dose distributions. The functions are often assumed to be S-shaped, with parameters that vary for different types of tumours, clinical endpoints for toxicity (measurable treatment outcomes, e.g. whether or not a patient suffers from a particular side-effect) and for individual cases. The TCP parameters depend on various aspects, such as the tumour size, the tumour cell density, the proliferation rate, the type of radiation and the level of hypoxia (Nahum and Kutcher, 2007). NTCP parameters depend on the clinical endpoint considered and on the volume-effect. The latter is organ-specific: for some organs, delivering a low dose to the whole organ is preferable over delivering a high dose to a small part of its volume, while for other organs this is the other way around.
Also the equivalent uniform dose (EUD), developed by Niemierko (1997), has been considered a useful objective function for inverse treatment plan optimization (Wu et al., 2002). The EUD is the dose that, if homogeneously delivered to the tumour, would yield the same fraction of surviving tumour cells as the planned non-homogeneous dose distribution. Romeijn et al. (2004) note that as some formulations of TCP are an increasing function of EUD, including the EUD for the respective structures yields the exact same optimal solutions as TCP.
