Home
Search
Collections
Journals
About
Contact us
My IOPscience
A global quality assurance system for personalized radiation therapy treatment planning for the prostate (or other sites)
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 Phys. Med. Biol. 59 5575 (http://iopscience.iop.org/0031-9155/59/18/5575) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 137.189.170.231 This content was downloaded on 19/02/2015 at 11:36
Please note that terms and conditions apply.
Institute of Physics and Engineering in Medicine Phys. Med. Biol. 59 (2014) 5575–5591
Physics in Medicine & Biology doi:10.1088/0031-9155/59/18/5575
A global quality assurance system for personalized radiation therapy treatment planning for the prostate (or other sites) Obioma Nwankwo, Dwi Seno K Sihono, Frank Schneider and Frederik Wenz Department of Radiation Oncology, Universitätsmedizin Mannheim, University of Heidelberg, Theodor-Kutzer-Ufer 1–3, D-68167, Mannheim, Germany Received 6 November 2013, revised 14 May 2014 Accepted for publication 28 July 2014 Published 29 August 2014 Abstract
Introduction: the quality of radiotherapy treatment plans varies across institutions and depends on the experience of the planner. For the purpose of intra- and inter-institutional homogenization of treatment plan quality, we present an algorithm that learns the organs-at-risk (OARs) sparing patterns from a database of high quality plans. Thereafter, the algorithm predicts the dose that similar organs will receive in future radiotherapy plans prior to treatment planning on the basis of the anatomies of the organs. The predicted dose provides the basis for the individualized specification of planning objectives, and for the objective assessment of the quality of radiotherapy plans. Materials and method: one hundred and twenty eight (128) Volumetric Modulated Arc Therapy (VMAT) plans were selected from a database of prostate cancer plans. The plans were divided into two groups, namely a training set that is made up of 95 plans and a validation set that consists of 33 plans. A multivariate analysis technique was used to determine the relationships between the positions of voxels and their dose. This information was used to predict the likely sparing of the OARs of the plans of the validation set. The predicted doses were visually and quantitatively compared to the reference data using dose volume histograms, the 3D dose distribution, and a novel evaluation metric that is based on the dose different test. Results: a voxel of the bladder on the average receives a higher dose than a voxel of the rectum in optimized radiotherapy plans for the treatment of prostate cancer in our institution if both voxels are at the same distance to the PTV. Based on our evaluation metric, the predicted and reference dose to the bladder agree to within 5% of the prescribed dose to the PTV in 18 out of 33 cases, while the predicted and reference doses to the rectum agree to within 5% in 28 out of the 33 plans of the validation set. 0031-9155/14/185575+17$33.00 © 2014 Institute of Physics and Engineering in Medicine Printed in the UK & the USA
5575
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Conclusion: We have described a method to predict the likely dose that OARs will receive before treatment planning. This prospective knowledge could be used to implement a global quality assurance system for personalized radiation therapy treatment planning. Keywords: personalized radiotherapy treatment planning, dose prediction algorithm, IMRT/VMAT treatment plan optimization, normal tissue sparing, knowledge based radiation therapy treatment planning S Online supplementary data available from stacks.iop.org/PMB/59/5575/ mmedia (Some figures may appear in colour only in the online journal) 1. Introduction The treatment plan quality across institutions that utilize similar treatment planning and delivery systems should be comparable. However, a previous study (Chung et al 2008) concluded that an experienced institution could produce superior Intensity Modulated Radiotherapy (IMRT) plans compared to a less experienced centre. Additionally, an objective decision on the quality of a treatment plan requires knowledge of what is dosimetrically achievable for the given anatomy (Appenzoller et al 2012). A recent publication (Good et al 2013) proposed knowledge-based radiation therapy (KBRT) treatment planning as a means of homogenizing treatment plan quality across institutions, by transferring planning expertise from the experienced to the less experience institutions. The KBRT treatment planning techniques (Wu et al 2009, Chanyavanich et al 2011, Zhu et al 2011a, Petit et al 2012, Yuan et al 2012, Appenzoller et al 2012, Good et al 2013, Wang et al 2013) aims at incorporating the knowledge that is learned from prior plans in the making and acceptance of future ones. Some advantages of utilizing prior knowledge in the treatment-planning and acceptance process include consistency of treatment plan quality across institutions (Chanyavanich et al 2011), less dependency of plan quality on the experience of planners, reduction in the treatment planning time and improved efficiency in the treatment plan optimization and acceptance procedures (Petit et al 2012). The KBRT treatment planning techniques vary in approach. Some approaches rely on the premise that the initial plan optimization constraints, which resulted in the desired coverage of the planning target volume (PTV) and yielded sufficient sparing of the organs at risk (OARs) in prior plans, could be reused for the making of new plans, provided that the geometries of the treated region are similar in both cases (Wu et al 2009, Chanyavanich et al 2011, Petit et al 2012, Good et al 2013). One approach (Chanyavanich et al 2011, Good et al 2013) relies on the beam’s eye view of the treated region for the assessment of the similarity of geometries, while another approach uses the overlap volume histogram (OVH) (Wu et al 2009, Petit et al 2012, Wang et al 2013) for this evaluation. Both approaches are based on the search and retrieval of prior treatment plans from a database of reference plans. An alternative to the search and retrieve approach has been described (Zhu et al 2011a, Yuan et al 2012, Yang et al 2013). These methods predict the likely dose volume histograms (DVHs) of the OARs of interest from their distance-to-target histograms (DTHs) (Zhu et al., 2011a, Yuan et al 2012). The DTH is a similar anatomical descriptor as the OVH (Kazhdan and McNutt 2011, Zhu et al 2011b). A learning algorithm was first trained on a database of DTHs and DVHs to derive organ-specific models. Thereafter, the models could predict the 5576
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
DVH of similar organs from their DTHs. The predicted DVHs could be used to specify treatment plan objectives or to objectively assess the quality of a plan. The quality of treatment plans that were generated by the KBRT approaches were reported to be non-inferior to the plans that were made by expert planners (Chanyavanich et al 2011, Good et al 2013). A limitation of the reviewed KBRT approaches is that at best, they predict the likely DVH of the OARs, which does not provide any spatial information about the dose distribution within the organ of interest. This work describes an alternative approach to the problem. The organ-specific models that are derived through our method predict the dose to the individual voxels of the OARs. We also describe an evaluation metric that is based on the dose difference test for comparing the predicted and reference (calculated with the treatment planning system (TPS)) dose of the OARs. Prostate plans were used for this study because firstly, the prostate is a popular site for testing KBRT methods based on literature analysis (Chanyavanich et al 2011, Zhu et al 2011a, Good et al 2013, Wang et al 2013, Yang et al 2013). Additionally, the prostate is the most frequently planned site for Intensity Modulated Radiation Therapy (IMRT)/volumetric modulated Arc Therapy (VMAT) techniques in our institution (16.5% for the period reviewed). Thus, prostate patients will benefit the most from a clinical implementation of IMRT/VMAT KBRT treatment planning in our institution. 2. Materials and method 2.1. Introduction to key variables
The following terms and notations are used frequently in this work and require introduction. Figure 1 can also aid the understanding of the introduced terms. 1. Distance-to-PTV (ri ): is the shortest Euclidean distance between a voxel of an OAR and the surface of the PTV, and i(i = 1…n) is the index of the voxel in the organ. 2. Slice level (Z i ): specifies the distance of a voxel of an OAR to the path of the treatment beam. This variable is calculated as the absolute minimum difference between the z-coordinate of the voxel and the z-coordinates of all the voxels of the PTV, multiplied by the slice thickness of the computerized tomography (CT) scan. In other words, Z i = min Z i OAR− {Z PTV} × T (1)
ZiOAR is the z-coordinate of the OAR voxel whose slice level is being computed, while set {Z PTV} contains the z-coordinates of all the voxels of the PTV. The variable T denotes the slice thickness of the CT scan, which was 0.5 cm for the analyzed plans. The z-coordinate defines the cranio-caudal direction of the body in the reference coordinate system. We assume voxels with Zi = 0, henceforth referred to as in-field voxels to be located in the direct path of the treatment beam, while voxels with Zi > 0, henceforth referred to as out-of-field voxels are assumed to be outside the direct path of the treatment beam (Yuan et al 2012) (see figure 1). 3. Distance-to-PTV bin ( [r j − 1, r j ) ): is a left-closed and right-open interval, where j(j = 1…M) is the index of a bin whereas rj − 1 and rj are real numbers that represent the lower and upper values of the distance-to-PTV bin respectively. Voxels with approximately the same distance-to-PTV value (bin size = 0.5 mm) are grouped into the same distance-to-PTV bin. 4. Mean-dose-at-distance (Dj Z ): is the mean dose (D ) of all the voxels that (i) belong to the same (jth) distance-to-PTV bin and (ii) also have the same slice level (Z) value. The formula for calculating this variable is provided in figure 1. 5577
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Figure 1. Illustration of the training method. The out-of-field voxel (Z i>0) are outside
the direct path of the treatment beam. The mean dose-at-distance function describes the relationship between the mean dose of the in-field voxels of an OAR and the distance (ri) between the voxels and the surface of the PTV. The slice weight function relates the dose of a voxel to its position relative to the path of the treatment beam
5. Mean-dose-at-distance function (f ): is a function of all the mean-dose-at-distance values of the in-field voxels (i.e. (Dj Z = 0)) against the centre, Rj, of the corresponding distance-to-PTV bin, where Rj = (r j − 1 + r j ) / 2. The mean-dose-at-distance function has an associated mean dose standard deviation function (g), which as the name implies, specifies the standard deviation that is associated with the mean dose-at-distance values. Both functions are further described in later sections. 5578
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
6. Unweighted dose (Ui ): is calculated from functions f and g and a random number generator. If the voxel is an in-field voxel, then this value is equal to the predicted voxel dose (described later). Otherwise, a weighting factor (known as the dose weight factor) that depends on the slice level of the voxel is used to scale the unweighted dose. The computation of this variable is further described in later sections. 7. Slice weight factor (wZ ): is the mean weight of the dose of all the voxels on a given slice level (Z = 0, 0.5…α, where α is the maximum possible value of Z) relative to the mean weight of the dose of all the in-field voxels (Z = 0). This derivation of this variable is described later. 8. Slice weight function (p): is the function of the slice weight factors against the corresponding slice level (i.e. wZ against the corresponding Z). The function has an associated slice weight standard deviation function (q). Both functions are further discussed in later sections. 9. Dose weight factor (wi ): is the weight that should be applied to the unweighted dose, Ui of a voxel on the basis of its distance to the path of the treatment beam (Z i ). This variable is calculated for each voxel from functions p and q and a random number generator. 10. Predicted voxel dose (Di p): is dose that is predicted for the ith voxel of an OAR by our algorithm on the basis of the position of the voxel (relative to the PTV and treatment beam) and prior experience (derived from studying reference treatment plans). It is derived applying the dose weight factor that was computed for the voxel to its unweighted dose. In other words, Di p= Uf × wi 11. Reference dose (Di p): is the dose that was calculated for the ith voxel of the OAR in the original treatment plan by the TPS. 2.2. Database of reference plans
One hundred and twenty eight (128) patient plans were selected from a database of prostate cancer cases that were treated at the Department of Radiation Oncology, University Medical Centre Mannheim, Germany. The target volumes in these Volumetric Modulated Arc Therapy (VMAT) (Stieler et al 2009) plans were the prostate and seminal vesicles. The plans were made with Monaco® (CMS, Elekta) TPS (versions 2.04.00, 3.20.01 and 3.30.00). The prescribed dose to the PTV was variable for the plans. The planning objectives were specified with an institutional template. A template for the specification of the planning goals for prostate cancer treatment in our institution is presented in table 1. The constraints of the OARs are appropriately scaled when the prescribed dose to the target is changed. The maximum number of rotations of the gantry per plan and control points per arc in the template are 2 and 120 respectively. These objectives are modified if the resultant plan is deemed to be undesirable. All the 128 plans were approved by both a physician and physicist for real patient treatment and are therefore considered as the best achievable compromise between the delivery of the prescribed dose to the PTV and the adequate sparing of the adjacent OARs. 2.3. Extraction of information from treatment plans
The plans were imported into MATLAB® (Mathworks Inc., Natick, MA, USA) computing environment using the Computational Environment for Radiotherapy Research (CERR) software (Deasy et al 2003). The following variables were extracted from each plan 1. The coordinates (x, y, z ) of each voxel of an OAR. The rectum and bladder are the two OARs that were considered in this work. This information was only used for visualization. 5579
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Table 1. A MONACO® planning template for prostate VMAT treatment at the
University Medical Centre Mannheim
Structure
Cost function
Reference dose (cGy)
Isoconstraint
P + SB
Target EUD Quadratic overdose Quadratic underdose
6000.0 6000.0
6000.0 70.0 70.0
Bladder
Parallel Serial
3000.0
40.0 3500.0
Rectum
Parallel Quadratic overdose
3000.0 3000.0
30.0 50.0
Patient
Quadratic overdose Quadratic overdose Conformality
3500.0 2200.0
120.0 60.0 0.7
P + SB = prostate and seminal vesicles
2. The reference dose, Di r, (normalized to the prescribed dose to the PTV) of the voxel. The subscript, i ( i = 1... n ), specifies the index of the voxel in the organ. The normalization of the dose renders the analysis to be independent of the absolute value of the prescribed dose to the PTV. The reference dose is the dose that was calculated for the voxel by MONACO® TPS. 3. The distance-to-PTV of the voxel (see the preceding section). Following the convention that was used in previous works (Wu et al 2009, Zhu et al 2011a, Yuan et al 2012), this value is signed negative if the voxel is shared by both the OAR and PTV. 4. The slice level of the voxel (see the preceding section). Two matrices, one for each OAR (rectum and bladder) were therefore computed from each plan (2 × 128 files). A row of the matrix represents a voxel while the columns specifies its parameters (coordinates, dose, distance-to-PTV and slice level). 2.4. The training set and the validation set
The 128 reference plans were randomly divided into two groups. The first group is henceforth referred to as the training set and consists of 95 plans, whereas the second group, which is henceforth referred to as the validation set is made up of 33 plans. For each OAR, the 95 matrices of the training set were merged to derive a single larger matrix of the organ. An organspecific model that relates the distances of voxels of the OARs to the PTV and path of the treatment beam, to their average dose was thereafter computed from this matrix. The plans of the validation set were used to assess the predictive accuracy of the derived model, and to evaluate the usefulness of this algorithm as a guiding tool for treatment planning. Each organspecific model is defined by two main functions (already described in the preceding section), namely the mean-dose-at-distance function and the slice weight function. Both functions have their associated standard deviation functions. The derivation of these functions is described in the following sections. 2.5. The mean dose-at-distance function (slice level = 0)
Figure 1 explains how the mean dose-at-distance function and the slice weight function are computed for each OAR from the training set. Consider a voxel as the multivariate object, (r , Z , D )i, where r is its distance-to-PTV, Z is the slice level and D is the reference dose (predicted by the 5580
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
TPS) of the voxel. The subscript i ( i = 1... N ) is the index of the voxel in the training set (after the merging of the 95 matrices). The mean dose-at-distance function of an organ is computed from the in-field voxels (i.e. having slice level = 0) of the OAR. The voxels were first binned (0.5 mm bin size according to their distance-to-PTV) value. The mean dose-at-distance value, Dj Z = 0 (also see preceding section), is thereafter computed for each distance-to-PTV bin. The formula for calculating this variable is given in figure 1. The mean-dose-at-distance function, f, is the function of the mean dose-at-distance values of the in-field voxels (the set containing all Dj Z = 0 values) against the centre of the corresponding distance-to-PTV bin (the set containing all Rj values). This function relates the mean dose of all the in-field voxels within a distance-toPTV bin to the median distance of the voxels to the surface of the PTV. The mean dose standard deviation function, g, is the function of the standard deviation Z=0 σ values against the centre of the corresponding distance( j ) of the mean dose-at-distance to-PTV bin. The variable, σj Z = 0 is calculated according to equation (2) σj Z = 0=
P
1 ∑(Dj Z = 0 − Di )2 P − 1 i=1
if ri ∈ [r j − 1, r j ) and Zi = 0
(2)
the variables i and P specifies the index of the in-field voxel in the jth distance-to-PTV bin and the total number of the in-field voxels in the bin respectively. 2.6. The slice weight function
We consider the out-of field voxels to be outside the path of the treatment beam (see the illustration in figure 1) in a coplanar beam arrangement, which is the case for prostate plans in our institution. We theorize that these voxels receive lower dose in optimized treatment plans when compared to the dose of the in-field voxels that are at the same distance to the PTV. We also hypothesize that the dose of the out-of-field voxels in a treatment plan can be determined by weighting the dose of the in-field voxels (if this is known) and vice versa. The slice weight function, p, specifies the weights that are required to scale the dose of a voxel according to its slice level. In contrast to the computation of the mean dose-at-distance function, all voxels (out-of field and in-field) of the training set of the organ of interest are used to compute the slice weight function according to the following steps. Firstly, the voxels of the training set (of the organ) were binned (0.5 mm bin size) according to their distance-to-PTV value. Voxels within a given distance-to-PTV bin were further separated according to their slice level. The mean dose-at-distance value Dj Z , is thereafter computed for each slice level (Z = 0.5, 1...Z max ) of a distance-to-PTV bin (j). This operation is performed for the entire distance-to-PTV bin. The mean dose-at-distance value of each slice level of a distance-to-PTV bin is then normalized to mean dose-at-distance value of the in-field voxels of that same bin. In other words Dj Z Dj Z = Z = 0 (3) Dj
A distance bin is ignored if it does not contain an in-field voxel since the described normalization cannot be performed. The slice weight factor, wZ, is calculated for each slice level as the average of the mean dose-at-distance values of that slice level in all the distance-to-PTV bins. The mathematical description of this operation (calculation of wZ) is provided in figure 1. The slice weight function, p, is the function of the slice weight factors against their corresponding slice levels. The slice-weight standard deviation function, q, specifies the standard deviation σZ, of the slice weight factors. The variable σZ is calculated according to equation (4) 5581
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
σZ =
M
1 ∑(wZ − Dj Z )2 M − 1 j=1
(4)
j(j = 0, 0.5…M ) is the index of the distance-to-PTV bin while Z (Z = 0.5, 1...Z max ) is the slice level. The mean-dose-at-distance function and the slice weight function of an OAR, together with their standard deviation functions (g and q) represent the model of learned relationship between the distance of a voxel of the organ to the PTV and treatment beam respectively, and the average dose of the voxel in optimized radiotherapy treatment plans. 2.7. Predicting the dose to a voxel from its position
The distance-to-PTV and the slice level (r , Z )i of a voxel are the two variables that are required to predict the dose of the voxel in an optimized treatment plan. Both variables are computed from the planning CT data. The subscript i(i = 1...n ) is the index of the voxel in the OAR. Firstly, the mean dose-at-distance value for the distance-to-PTV bin of the voxel is calculated from the mean dose-at-distance function according to equation (5) D j = f (ri ) (5) Z=0
Thereafter, standard deviation, σj Z = 0, of this value is computed from the mean dose standard deviation function, g, as shown in equation (6) σj Z = 0 = g (ri ) (6)
The unweighted dose, Ui, of an in-field voxel with a distance-to-PTV value equal to ri in an optimized treatment plan is estimated from Dj Z = 0 and σj Z = 0 according to equation (7) Z=0 U , σj Z = 0) (7) l = RNG1 (Dj
Function RNG1 randomly generates a number with a mean centered at Dj Z = 0 and a standard deviation that is defined by σj Z = 0 (RNG1 = normrnd() function in MATLAB). The unweighted dose is scaled if the voxel is an out-of-field voxel. The weighting factor depends on the slice level of the voxel. Firstly, the slice weight factor is calculated from the slice weight function according to equation (8) w (8) Z = p (Z i )
wZ and Zi are the slice weight factors and the slice level of the voxel (the input to the function) respectively. The standard deviation, σZ, of the slice weight factor is calculated from the slice weight standard deviation function q, as specified by equation (9) σ(9) Z = q (Z i )
The voxel-specific dose weight factor, wi, is calculated from the slice weight factor and its standard deviation according to equation (10) w (10) i = RNG2 (wZ , σZ )
The function RNG2 is similar to the previously described RNG1. The predicted dose, Di P, of the voxel is derived by applying the dose weight factor of the voxel to its unweighted dose as expressed in equation (11) DiP = Ui × wi (11) 5582
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
The value of the dose weight factor (wi ) is always equal to 1 for in-field voxels and by implication P U if Z i = 0 i = Di (12)
2.8. Plan evaluation metrics: the agreement scores
The comparison of the predicted and reference dose distributions starts at the voxel level. Firstly, the dose difference, Qi, between the reference, Di r, and predicted, Di P, dose to the ith voxel is calculated as the signed difference between both values P r Q i = Di −Di (13)
The organ agreement score, Sl organ, is the mean of the dose difference values of all the (n) voxels of the organ as described by equation (14) n
1 organ S(14) = ∑ Qi l n i=1
The subscript index l (l = 1, 2 ...V ) specifies the index of the organ in the patient. The sign on the organ agreement score indicates the direction of any difference; a negative score implies that the predicted dose to the organ is higher than the reference dose and vice versa. The plan agreement score, Splan, is calculated as the sum of the agreement scores for all the (V) OARs that are considered in the analysis V
organ S(15) plan = ∑ Sl l=1
2.9. Refinement of the training data and optimization of the organ models
The initial organ models that were generated from the training data (95 training data for each organ) were used as an ‘outlier detector’ to identify and remove suboptimal data from the training set of the organ. Refined models were thereafter generated from this modified training data. The outlier detection and removal process is as follows. The initial model was used to predict the dose to the individual plan data that make up the training set. The organ agreement score was thereafter computed for each plan. Since the predicted dose is subtracted from the reference dose in the computation of the organ agreement score (equations (13) and (14)), a positive score implies that the reference dose to the organ is higher than the dose predicted by the model. It should be borne in mind that the prediction of the model represents the average sparing pattern of the plans of the training set. A positive score (depending on the magnitude) therefore implies that the sparing achieved for the organ in the plan is lower than the group (training set) average. Our approach to optimizing the training data and refining the model is to remove plans whose organ agreement scores exceed a certain positive threshold value. 2.10. Validation of the model
The refined models were used to predict the likely 3D dose distribution in the OARs of the validation set based on their anatomies. The DVHs, organ agreement scores and the plan 5583
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
i
120
ii
1 Bladder Rectum
Bladder Rectum 0.9
100
0.7 80
Slice weight factor
Mean dose−at−distance [%]
0.8
60
40
0.6
0.5
0.4
0.3
0.2 20 0.1
0 −5
0
5
10
Distance−to−PTV [cm]
15
0 −2
0
2
4
Slice level [cm]
6
8
Figure 2. Subplot i shows the mean dose-at-distance functions of the rectum and blad-
der for prostate plan. The functions relate the mean dose that is received by voxels of these organs in optimized treatment plans for the prostate to the distances between the voxels and the PTV. The standard deviations of these mean values are plotted as error bars. The slice weight functions (subplot ii) relates the dose of voxels to the distances between the voxels a and the treatment beam direction. The standard deviations of the slice weight factors are plotted as error bars. These functions represent the organ-specific models for predicting the likely 3D dose distribution in these organs in optimized radiotherapy plans for the treatment of prostate cancer.
agreement scores were computed for the plans for the qualitative and quantitative comparison of the two (reference and predicted) datasets. 3. Results 3.1. Refined training set
The training set was refined by excluding plans whose organ agreement scores are >10%. The refined training set of the bladder and rectum contained 90 and 71 training data respectively. In other words, the final models of the bladder and rectum were generated from 90 and 71 plans respectively. 3.2. Relationship between position of a voxel and its mean dose (training set)
Figure 2 presents the relationship between the distance of voxels to the surface of the PTV and path of the treatment beam, and their average dose in optimized plans for the treatment of prostrate cancer. Subplot i shows that the mean dose-at-distance function of the bladder and rectum and their associated standard deviations (as error bars). The mean dose mean doseat-distance function of the bladder is higher than that of the rectum. This implies that for the treatment of prostate in our institution, a voxel of the bladder will receive a higher dose than a voxel of the rectum in an optimized plan if both voxels are at the same distance to the PTV. 5584
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Table 2. Comparison of the predicted and reference dose distribution by means of the
agreement scores. The organ agreement score quantifies the agreement between the predicted and reference 3D dose to an OAR. The plan agreement score is the sum of all the organ agreement scores of all the organs that are considered.
Plan index
Bladder agreement score (%)
Rectum agreement score (%)
Plan agreement score (%)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
−7.16 5.13 −6.30 4.16 8.03 1.77 −2.12 3.12 1.03 −1.15 4.03 −9.06 3.12 3.10 7.48 12.19 −1.58 8.19 0.91 4.51 4.22 1.35 3.83 −9.89 10.04 5.87 −2.36 −6.23 −2.82 6.37 6.36 −0.54 −0.26
6.83 4.45 3.00 −0.58 4.55 1.81 7.33 −1.35 0.93 2.05 2.42 −3.40 −1.21 −0.23 −1.63 4.36 2.75 −4.82 −4.11 4.28 1.94 −2.83 8.22 −1.54 −6.57 −0.35 3.39 −4.28 7.00 4.03 4.33 4.12 −2.34
−0.33 9.58 −3.30 3.59 12.57 3.58 5.21 1.77 1.97 0.90 6.45 −12.46 1.90 2.87 5.85 16.55 1.18 3.37 −3.20 8.79 6.16 −1.47 12.05 −11.43 3.46 5.52 1.03 −10.51 4.18 10.39 10.69 3.58 −2.61
Subplots ii present the slice weight function of the OARs and their standard deviations (as error bars). The function relates the weight of the dose of a voxel to its distance to the path of treatment beam. This means that if two voxels of an OAR are at the same distance to the PTV but are located on different slice levels, then given the dose of either voxel, the unknown dose can be estimated from the known dose and the slice weight function. For example, the slice weight factor for the bladder is 1 when the slice level is 0 (in-field voxel) and 0.2 (ignoring the standard deviation) when the slice level is 1.5 cm. This implies that if two voxels (a and b) of the bladder are at the same distance to the target, but with voxel a in the path of the treatment beam (Z = 0) and voxel b at a distance (slice level) of 1.5 cm away from the path of the beam, then the dose of both voxels are related by equation (16) (16) Db = Da × 0.2 Da and Db represents the unweighted doses of voxels a and b respectively (the unweighted dose of an in-field voxel is however equal its predicted dose as previously described in equation (12)). Polynomial approximations of these functions (figure 2) define the organ-specific models that were derived from the refined training set.
5585
O Nwankwo et al
0.5
0
0
50
Dose [%]
100
1.03
1
9
0.5
0
0
50
Dose [%]
100
3.12
1
13
0.5
0
0
50
Dose [%]
100
0
50
Dose [%]
1.77
1
6
0.5
0
0
50
Dose [%]
100
−1.15
1
10
0.5
0
0
50
Dose [%]
100
3.10
1
14
0.5
0
0
50
Dose [%]
100
0.5
0
0
50
Dose [%]
100
−2.12
1
7
0.5
0
0
50
Dose [%]
100
4.03
1
11
0.5
0
0
50
Dose [%]
100
7.48
1
15
0.5
0
Volume [Fraction]
Volume [Fraction]
100
3
Volume [Fraction]
5
Volume [Fraction]
8.03
1
0
−6.30
Volume [Fraction]
100
1
0
50
Dose [%]
Volume [Fraction]
50
Dose [%]
0.5
Volume [Fraction]
0
Volume [Fraction]
0
Predicted Reference
2
Volume [Fraction]
0.5
Bladder
5.13
1
Volume [Fraction]
1
Volume [Fraction]
−7.16
1
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Phys. Med. Biol. 59 (2014) 5575
100
4.16
1
4
0.5
0
0
50
Dose [%]
100
3.12
1
8
0.5
0
0
50
Dose [%]
100
−9.06
1
12
0.5
0
0
50
Dose [%]
100
12.19
1
16
0.5
0
0
50
Dose [%]
100
Figure 3. Validation of the organ model of the bladder using the plans of the validation set. The DVHs were computed from the predicted and reference 3D dose distributions of the plans. Based on the bladder agreement score, plan 9 showed the best agreement between the compared datasets while plan 16 showed the worst agreement.
3.3. Comparison of the predicted and reference datasets (validation set)
Table 2 quantitatively compares the predicted and reference dose distributions for the plans of the validation set by means of the organ agreement and the plan agreement scores, while figures 3 and 4 compares 16 plans of both datasets using their DVHs. Judging from the bladder agreement score (second column of table 2), plans 16 and 25 have worst agreement between the predicted and reference dose to the bladder, while plans 32 and 33 show the best agreement. The rectum agreement scores suggests that plans 4 and 14 have the best agreement between the predicted and reference datasets dose to the bladder. The plan agreement score suggests that Plan 1 has the best agreement between the predicted and reference datasets when the total dose to both organs are considered. The DVHs (figures 3 and 4) were computed from the predicted and reference 3D dose distributions in of OARs. Only the DVH of the first 16 plans of the validation set are presented for the sake of space. The integer values in the subplots specifies the index of the plan while the title of the plots show the agreement score for the organ. Figure 3 compares the predicted and reference DVHs of the bladder for 16 plans of the validation set. Based on the bladder agreement scores, plans 9 and 10 show the best agreement between the predicted and reference dose (values closest to zero), while plans 12 and 16 show the worst agreement (values farthest from zero). 5586
O Nwankwo et al
5
0.5
0
0
50
Dose [%]
100
0.93
1
9
0.5
0
0
50
Dose [%]
100
−1.21
1
13
0.5
0
0
50
Dose [%]
100
0
50
Dose [%]
1.81
1
6
0.5
0
0
50
Dose [%]
100
2.05
1
10
0.5
0
0
50
Dose [%]
100
−0.23
1
14
0.5
0
0
50
Dose [%]
100
0.5
0
0
50
Dose [%]
100
7.33
1
7
0.5
0
0
50
Dose [%]
100
2.42
1
11
0.5
0
0
50
Dose [%]
100
−1.63
1
15
0.5
0
Volume [Fraction]
Volume [Fraction]
100
3
Volume [Fraction]
1
Volume [Fraction]
4.55
0
3.00
Volume [Fraction]
100
1
0
50
Dose [%]
100
Volume [Fraction]
50
Dose [%]
0.5
Volume [Fraction]
0
Volume [Fraction]
0
2
Volume [Fraction]
0.5
Rectum
4.45
1
Volume [Fraction]
Predicted 1 Reference
Volume [Fraction]
6.83
1
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Phys. Med. Biol. 59 (2014) 5575
−0.58
1
4
0.5
0
0
50
Dose [%]
100
−1.35
1
8
0.5
0
0
50
Dose [%]
100
−3.40
1
12
0.5
0
0
50
Dose [%]
100
4.36
1
16
0.5
0
0
50
Dose [%]
100
Figure 4. Comparison of the reference and predicted DVHs of the rectum of the valida-
tion set. The organ agreement score (title of subplots) is a quantitative measure of the accuracy between the predicted and reference doses.
The quantitative scores agree with the visual ranking of the DVHs of these plans. This concordance reflects the accuracy of this metric in quantifying the level of agreement between the compared datasets. Figure 4 presents the comparison of the predicted and reference DVHs of the rectum. The rectum agreement score suggests that plan 14 has the best agreement between the predicted and reference data, while plans 1 and 7 show the least agreement. These results closely agree with the visual ranking of the agreement between the predicted and reference DVHs, thus providing further proof that the organ agreement score is a useful metric for quantifying the agreement between the predicted and reference dose distributions. A uniqueness of our approach relative to the reviewed methods is the ability of our algorithm to predict the achievable 3D dose distribution in the OARs of interest. Figure 5 shows the reference and predicted dose distributions for the OARs of Plan 8. The bladder and rectum agreement score are 3.12 and −1.35 respectively. The bladder (upper subfigures) and rectum (lower subfigures) of this plan consistes of 3.830 7 × 104 and 1.218 1 × 104 voxels respectively. Visual analysis (organ agreement scores and plan agreement score provides quantitative assessment) shows good agreement between both datasets. The predicted and
5587
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Figure 5. Visualization of the predicted and reference 3D dose distributions of the OARs of Plan 8. The reference dose distribution was calculated with the TPS while the predicted dose was calculated with our algorithm on the basis of the anatomical relationships between the organs and the target volume.
reference 3D dose distributions of all 33 plans of the validation set are provided as a supplementary data (stacks.iop.org/PMB/59/5575/mmedia). 4. Discussion Our aim in this work is to learn the relationship between the anatomies of organs at risk (OARs) and their dose in optimized radiotherapy treatment plans. This information was learned from a database of high quality previously treated patient plans. Thereafter, this knowledge is used to predict the likely dose of similar OARs on the basis of their anatomies. This prospective knowledge could be used for the anatomy-specific specification of treatment planning objectives with respect to normal tissue sparing, as well as a quality assurance (QA) tool to aid the objective assessment of the quality of radiotherapy plans on the basis of anatomy. The use of templates to specify planning objectives does not consider the individual patient anatomical variations. In addition, an objective decision on the quality of a radiotherapy treatment plan requires the quantitative knowledge of what can be achieved for that particular anatomy (Appenzoller et al 2012). Prospective knowledge of the achievable OAR sparing for a given anatomy is therefore an essential requirement for the implementation of an objective and individualized radiation therapy treatment-planning and acceptance scheme. We 5588
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
have described a method to learn the OAR sparing patterns from a database of high quality plans. This knowledge is then used to predict the likely dose that the organ will receive in an optimized treatment plan within a specific framework (similar planning conditions as the reference plans). Our approach is unique because it predicts the 3D dose in the OARs unlike previously described methods. The method that was used to derive the organ-specific models described in this work is similar to the technique that was used to create a virtual source model of a kilo-voltage radiotherapy device from a phase space file (Nwankwo et al 2013). Both are multivariate analyses methods, which operate by binning multivariate observations according to one variable (distance-to-PTV value of voxels in this work and the x-coordinates of photons in the cited work). Thereafter, the mean value of a second variable (dose in this work and direction cosine in the cited work) is computed for all the observations that fall within a bin, thus establishing a relationship between both variables (i.e. mean dose versus distance-to-PTV). The visual and quantitative comparison of the predicted and reference dose distributions show that this approach is accurate in predicting the 3D dose in the OARs, and that the evaluation metrics are accurate in quantifying the agreement between the compared datasets. It is possible to achieve a better than average sparing for an OAR at the expense of a competing organ through the adjustment of the prescribed dose to the OARs (increasing the penalty dose to one OAR while reducing or keeping constant the assigned penalty values of the competing OAR). This factor should be taken into account when comparing the predicted and reference dose distribution. It is for this reason that the plan agreement score is necessary. Plan 1 illustrates the usefulness of the plan agreement score metric. The bladder agreement score of Plan 1 is −7.16% while the rectum agreement score is 6.83% (table 2). This implies that the predicted to the bladder is higher than the reference dose, and that the reference dose to the rectum is higher than the prediction by our algorithm. This is suggestive of a prioritized sparing of the bladder at the expense of the rectum in the reference plan. This possible prioritized sparing of an OAR at the expense of another is accounted for by the plan agreement score. There are several potential applications of our algorithm in radiation therapy. One such application is the comparison of the OAR sparing of parallel treatment techniques. A previous study (Wolff et al 2009) used the mean dose to the anterior and posterior portions of the rectum, in addition to other metrics to compare the OAR sparing of parallel treatment techniques (VMAT, step and shoot IMRT, serial tomotherapy and 3D conformal radiotherapy). However, a restriction on using the mean organ dose and dose-to-volume in such analysis is that the same patient data should be used for the different techniques. This could be challenging if a large number of plans are to be analyzed. The mean-dose-at-distance function could be used to compare different treatment techniques even if the patient data is variable. This is possible because a superior treatment technique will deliver lower dose to the OARs (for the same PTV coverage) relative to an inferior technique. The difference in the achievable OAR sparing of the compared techniques will manifest in their respective organ-models (meandos-at-distance and slice weight functions). Such comparison does not require similar patient data because anatomy (distance-to-PTV) is already included in the model, thereby making it possible to objectively compare plans from different techniques, irrespective of whether or not the patient data is the same. Additionally and on the basis of the foregoing, this approach could be used to compare the treatment plan quality and or techniques across institutions. The approach that we have described in this work is not limited to the prostate and can be applied to other treatment sites. An important step in knowledge based radiation therapy treatment planning is the introduction of the overlap volume histogram (Wu et al 2009), distanceto-target-histogram and distance-to-PTV (Zhu et al 2011a). These authors demonstrated that 5589
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
the complex 3D OAR-target relationship could be reduced to the OVH/DTH of the organ, and further to the distance-to-PTV values of the voxels that make up the organ. Our approach of using the mean value of prior observation (mean dose of voxels at a given distance to PTV) to estimate the likelihood of future outcome can be applied to predict the likely dose of OARs at other treatment sites, for example in the head and neck region. Our future works will focus on the clinical validation and applications of our algorithm. It is possible to generate organ-models of desired properties by the iterative removal or addition of plans of certain characteristics to training set. The prediction of the model can be made increasing strict by removing plans with high positive organ agreement scores from the training set while keeping plans with low or negative agreement scores in the training set. Figure 2 provides information regarding the likely limitations of the current model. We will describe this likely limitation with the mean-dose-at-distance function of the bladder. The maximum distance-to-PTV value that the model was trained on is approximately 13.0 cm. Thus, the predictive accuracy of the model for voxels with distance-to-PTV values >13.0 cm is uncertain and depends on the behavior of the polynomial fit function beyond the available training data. A robust and generalized model requires the use of as large a training data as possible. The overprediction by the models may be the result of superior planning expertise, compromised coverage of the target to achieve a better than average sparing of the OARs, or the influence of features that were not captured by the model. 5. Conclusion We have described an algorithm to predict the likely 3D dose distribution in OARs before treatment planning. This prospective knowledge could be used to implement an objective and individualized quality assurance system for personalized radiation therapy treatment planning. Supplementary data Predicted and reference 3D dose distributions for all the plans of the validation set. Acknowledgment The authors wish to acknowledge the contribution of Ms Kerstin Siebenlist, physicists and physicians of University Medical Centre Mannheim, who either planned or reviewed the analyzed plans. Obioma Nwankwo also acknowledge to have been introduced to knowledge based radiation therapy treatment planning technique in 2010 by Q J Wu, X Zhu, and T Li of Duke University Medical Centre, Durham. References Appenzoller L M, Michalski J M, Thorstad W L, Mutic S and Moore K L 2012 Predicting dose-volume histograms for organs-at-risk in IMRT planning Med. Phys. 39 7446 Chanyavanich V, Das S K, Lee W R and Lo J Y 2011 Knowledge-based IMRT treatment planning for prostate cancer Med. Phys. 38 2515 Chung H T, Lee B, Park E, Lu J J and Xia P 2008 Can all centers plan intensity-modulated radiotherapy (IMRT) effectively? An external audit of dosimetric comparisons between 3D conformal radiotherapy and IMRT for adjuvant chemoradiation for gastric cancer Int. J. Radiat. Oncol. Biol. Phys. 71 1167–74 5590
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Deasy J O, Blanco A I and Clark V H 2003 CERR: a computational environment for radiotherapy research Med. Phys. 30 979 Good D, Lo J, Lee W R, Wu Q J, Yin F-F and Das S K 2013 A knowledge-based approach to improving and homogenizing intensity modulated radiation therapy planning quality among treatment centers: an example application to prostate cancer planning Int. J. Radiat. Oncol. Biol. Phys. 87 176–81 Kazhdan M, McNutt T, Taylor R, Wu B and Simari P 2011 Comment on ‘A planning quality evaluation tool for prostate adaptive IMRT based on machine learning’ Med. Phys. 38 2820 Nwankwo O, Clausen S, Schneider F and Wenz F 2013 A virtual source model of a kilo-voltage radiotherapy device Phys. Med. Biol. 58 2363 Petit S F, Wu B, Kazhdan M, Dekker A, Simari P, Kumar R, Taylor R, Herman J M and McNutt T 2012 Increased organ sparing using shape-based treatment plan optimization for intensity modulated radiation therapy of pancreatic adenocarcinoma Radiother. Oncol. 102 38–44 Stieler F, Wolff D, Lohr F, Steil V, Abo-Madyan Y, Lorenz F, Wenz F and Mai S 2009 A fast radiotherapy paradigm for anal cancer with volumetric modulated arc therapy (VMAT) Radiat. Oncol. 4 25–48 Wang Y, Zolnay A, Incrocci L, Joosten H, McNutt T, Heijmen B and Petit S 2013 A quality control model that uses PTV-rectal distances to predict the lowest achievable rectum dose, improves IMRT planning for patients with prostate cancer Radiother. Oncol. 107 352–7 Wolff D, Stieler F, Welzel G, Lorenz F, Abo-Madyan Y, Mai S, Herskind C, Polednik M, Steil V and Wenz F 2009 Volumetric modulated arc therapy (VMAT) vs. serial tomotherapy, step-and-shoot IMRT and 3D-conformal RT for treatment of prostate cancer Radiother. Oncol. 93 226–33 Wu B, Ricchetti F, Sanguineti G, Kazhdan M, Simari P, Chuang M, Taylor R, Jacques R and McNutt T 2009 Patient geometry-driven information retrieval for IMRT treatment plan quality control Med. Phys. 36 5497 Yang Y, Ford E C, Wu B, Pinkawa M, van Triest B, Campbell P, Song D Y and McNutt T R 2013 An overlap-volume-histogram based method for rectal dose prediction and automated treatment planning in the external beam prostate radiotherapy following hydrogel injection Med. Phys. 40 011709 Yuan L, Ge Y, Lee W R, Yin F F, Kirkpatrick J P and Wu Q J 2012 Quantitative analysis of the factors which affect the interpatient organ-at-risk dose sparing variation in IMRT plans Med. Phys. 39 6868 Zhu X, Ge Y, Li T, Thongphiew D, Yin F-F and Wu Q J 2011a A planning quality evaluation tool for prostate adaptive IMRT based on machine learning Med. Phys. 38 719 Zhu X, Li T, Yin F-F, Wu Q J and Ge Y 2011b Response to ‘Comment on A planning quality evaluation tool for prostate adaptive IMRT based on machine learning’ Med. Phys. 38 2821
5591
Search
Collections
Journals
About
Contact us
My IOPscience
A global quality assurance system for personalized radiation therapy treatment planning for the prostate (or other sites)
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 Phys. Med. Biol. 59 5575 (http://iopscience.iop.org/0031-9155/59/18/5575) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 137.189.170.231 This content was downloaded on 19/02/2015 at 11:36
Please note that terms and conditions apply.
Institute of Physics and Engineering in Medicine Phys. Med. Biol. 59 (2014) 5575–5591
Physics in Medicine & Biology doi:10.1088/0031-9155/59/18/5575
A global quality assurance system for personalized radiation therapy treatment planning for the prostate (or other sites) Obioma Nwankwo, Dwi Seno K Sihono, Frank Schneider and Frederik Wenz Department of Radiation Oncology, Universitätsmedizin Mannheim, University of Heidelberg, Theodor-Kutzer-Ufer 1–3, D-68167, Mannheim, Germany Received 6 November 2013, revised 14 May 2014 Accepted for publication 28 July 2014 Published 29 August 2014 Abstract
Introduction: the quality of radiotherapy treatment plans varies across institutions and depends on the experience of the planner. For the purpose of intra- and inter-institutional homogenization of treatment plan quality, we present an algorithm that learns the organs-at-risk (OARs) sparing patterns from a database of high quality plans. Thereafter, the algorithm predicts the dose that similar organs will receive in future radiotherapy plans prior to treatment planning on the basis of the anatomies of the organs. The predicted dose provides the basis for the individualized specification of planning objectives, and for the objective assessment of the quality of radiotherapy plans. Materials and method: one hundred and twenty eight (128) Volumetric Modulated Arc Therapy (VMAT) plans were selected from a database of prostate cancer plans. The plans were divided into two groups, namely a training set that is made up of 95 plans and a validation set that consists of 33 plans. A multivariate analysis technique was used to determine the relationships between the positions of voxels and their dose. This information was used to predict the likely sparing of the OARs of the plans of the validation set. The predicted doses were visually and quantitatively compared to the reference data using dose volume histograms, the 3D dose distribution, and a novel evaluation metric that is based on the dose different test. Results: a voxel of the bladder on the average receives a higher dose than a voxel of the rectum in optimized radiotherapy plans for the treatment of prostate cancer in our institution if both voxels are at the same distance to the PTV. Based on our evaluation metric, the predicted and reference dose to the bladder agree to within 5% of the prescribed dose to the PTV in 18 out of 33 cases, while the predicted and reference doses to the rectum agree to within 5% in 28 out of the 33 plans of the validation set. 0031-9155/14/185575+17$33.00 © 2014 Institute of Physics and Engineering in Medicine Printed in the UK & the USA
5575
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Conclusion: We have described a method to predict the likely dose that OARs will receive before treatment planning. This prospective knowledge could be used to implement a global quality assurance system for personalized radiation therapy treatment planning. Keywords: personalized radiotherapy treatment planning, dose prediction algorithm, IMRT/VMAT treatment plan optimization, normal tissue sparing, knowledge based radiation therapy treatment planning S Online supplementary data available from stacks.iop.org/PMB/59/5575/ mmedia (Some figures may appear in colour only in the online journal) 1. Introduction The treatment plan quality across institutions that utilize similar treatment planning and delivery systems should be comparable. However, a previous study (Chung et al 2008) concluded that an experienced institution could produce superior Intensity Modulated Radiotherapy (IMRT) plans compared to a less experienced centre. Additionally, an objective decision on the quality of a treatment plan requires knowledge of what is dosimetrically achievable for the given anatomy (Appenzoller et al 2012). A recent publication (Good et al 2013) proposed knowledge-based radiation therapy (KBRT) treatment planning as a means of homogenizing treatment plan quality across institutions, by transferring planning expertise from the experienced to the less experience institutions. The KBRT treatment planning techniques (Wu et al 2009, Chanyavanich et al 2011, Zhu et al 2011a, Petit et al 2012, Yuan et al 2012, Appenzoller et al 2012, Good et al 2013, Wang et al 2013) aims at incorporating the knowledge that is learned from prior plans in the making and acceptance of future ones. Some advantages of utilizing prior knowledge in the treatment-planning and acceptance process include consistency of treatment plan quality across institutions (Chanyavanich et al 2011), less dependency of plan quality on the experience of planners, reduction in the treatment planning time and improved efficiency in the treatment plan optimization and acceptance procedures (Petit et al 2012). The KBRT treatment planning techniques vary in approach. Some approaches rely on the premise that the initial plan optimization constraints, which resulted in the desired coverage of the planning target volume (PTV) and yielded sufficient sparing of the organs at risk (OARs) in prior plans, could be reused for the making of new plans, provided that the geometries of the treated region are similar in both cases (Wu et al 2009, Chanyavanich et al 2011, Petit et al 2012, Good et al 2013). One approach (Chanyavanich et al 2011, Good et al 2013) relies on the beam’s eye view of the treated region for the assessment of the similarity of geometries, while another approach uses the overlap volume histogram (OVH) (Wu et al 2009, Petit et al 2012, Wang et al 2013) for this evaluation. Both approaches are based on the search and retrieval of prior treatment plans from a database of reference plans. An alternative to the search and retrieve approach has been described (Zhu et al 2011a, Yuan et al 2012, Yang et al 2013). These methods predict the likely dose volume histograms (DVHs) of the OARs of interest from their distance-to-target histograms (DTHs) (Zhu et al., 2011a, Yuan et al 2012). The DTH is a similar anatomical descriptor as the OVH (Kazhdan and McNutt 2011, Zhu et al 2011b). A learning algorithm was first trained on a database of DTHs and DVHs to derive organ-specific models. Thereafter, the models could predict the 5576
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
DVH of similar organs from their DTHs. The predicted DVHs could be used to specify treatment plan objectives or to objectively assess the quality of a plan. The quality of treatment plans that were generated by the KBRT approaches were reported to be non-inferior to the plans that were made by expert planners (Chanyavanich et al 2011, Good et al 2013). A limitation of the reviewed KBRT approaches is that at best, they predict the likely DVH of the OARs, which does not provide any spatial information about the dose distribution within the organ of interest. This work describes an alternative approach to the problem. The organ-specific models that are derived through our method predict the dose to the individual voxels of the OARs. We also describe an evaluation metric that is based on the dose difference test for comparing the predicted and reference (calculated with the treatment planning system (TPS)) dose of the OARs. Prostate plans were used for this study because firstly, the prostate is a popular site for testing KBRT methods based on literature analysis (Chanyavanich et al 2011, Zhu et al 2011a, Good et al 2013, Wang et al 2013, Yang et al 2013). Additionally, the prostate is the most frequently planned site for Intensity Modulated Radiation Therapy (IMRT)/volumetric modulated Arc Therapy (VMAT) techniques in our institution (16.5% for the period reviewed). Thus, prostate patients will benefit the most from a clinical implementation of IMRT/VMAT KBRT treatment planning in our institution. 2. Materials and method 2.1. Introduction to key variables
The following terms and notations are used frequently in this work and require introduction. Figure 1 can also aid the understanding of the introduced terms. 1. Distance-to-PTV (ri ): is the shortest Euclidean distance between a voxel of an OAR and the surface of the PTV, and i(i = 1…n) is the index of the voxel in the organ. 2. Slice level (Z i ): specifies the distance of a voxel of an OAR to the path of the treatment beam. This variable is calculated as the absolute minimum difference between the z-coordinate of the voxel and the z-coordinates of all the voxels of the PTV, multiplied by the slice thickness of the computerized tomography (CT) scan. In other words, Z i = min Z i OAR− {Z PTV} × T (1)
ZiOAR is the z-coordinate of the OAR voxel whose slice level is being computed, while set {Z PTV} contains the z-coordinates of all the voxels of the PTV. The variable T denotes the slice thickness of the CT scan, which was 0.5 cm for the analyzed plans. The z-coordinate defines the cranio-caudal direction of the body in the reference coordinate system. We assume voxels with Zi = 0, henceforth referred to as in-field voxels to be located in the direct path of the treatment beam, while voxels with Zi > 0, henceforth referred to as out-of-field voxels are assumed to be outside the direct path of the treatment beam (Yuan et al 2012) (see figure 1). 3. Distance-to-PTV bin ( [r j − 1, r j ) ): is a left-closed and right-open interval, where j(j = 1…M) is the index of a bin whereas rj − 1 and rj are real numbers that represent the lower and upper values of the distance-to-PTV bin respectively. Voxels with approximately the same distance-to-PTV value (bin size = 0.5 mm) are grouped into the same distance-to-PTV bin. 4. Mean-dose-at-distance (Dj Z ): is the mean dose (D ) of all the voxels that (i) belong to the same (jth) distance-to-PTV bin and (ii) also have the same slice level (Z) value. The formula for calculating this variable is provided in figure 1. 5577
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Figure 1. Illustration of the training method. The out-of-field voxel (Z i>0) are outside
the direct path of the treatment beam. The mean dose-at-distance function describes the relationship between the mean dose of the in-field voxels of an OAR and the distance (ri) between the voxels and the surface of the PTV. The slice weight function relates the dose of a voxel to its position relative to the path of the treatment beam
5. Mean-dose-at-distance function (f ): is a function of all the mean-dose-at-distance values of the in-field voxels (i.e. (Dj Z = 0)) against the centre, Rj, of the corresponding distance-to-PTV bin, where Rj = (r j − 1 + r j ) / 2. The mean-dose-at-distance function has an associated mean dose standard deviation function (g), which as the name implies, specifies the standard deviation that is associated with the mean dose-at-distance values. Both functions are further described in later sections. 5578
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
6. Unweighted dose (Ui ): is calculated from functions f and g and a random number generator. If the voxel is an in-field voxel, then this value is equal to the predicted voxel dose (described later). Otherwise, a weighting factor (known as the dose weight factor) that depends on the slice level of the voxel is used to scale the unweighted dose. The computation of this variable is further described in later sections. 7. Slice weight factor (wZ ): is the mean weight of the dose of all the voxels on a given slice level (Z = 0, 0.5…α, where α is the maximum possible value of Z) relative to the mean weight of the dose of all the in-field voxels (Z = 0). This derivation of this variable is described later. 8. Slice weight function (p): is the function of the slice weight factors against the corresponding slice level (i.e. wZ against the corresponding Z). The function has an associated slice weight standard deviation function (q). Both functions are further discussed in later sections. 9. Dose weight factor (wi ): is the weight that should be applied to the unweighted dose, Ui of a voxel on the basis of its distance to the path of the treatment beam (Z i ). This variable is calculated for each voxel from functions p and q and a random number generator. 10. Predicted voxel dose (Di p): is dose that is predicted for the ith voxel of an OAR by our algorithm on the basis of the position of the voxel (relative to the PTV and treatment beam) and prior experience (derived from studying reference treatment plans). It is derived applying the dose weight factor that was computed for the voxel to its unweighted dose. In other words, Di p= Uf × wi 11. Reference dose (Di p): is the dose that was calculated for the ith voxel of the OAR in the original treatment plan by the TPS. 2.2. Database of reference plans
One hundred and twenty eight (128) patient plans were selected from a database of prostate cancer cases that were treated at the Department of Radiation Oncology, University Medical Centre Mannheim, Germany. The target volumes in these Volumetric Modulated Arc Therapy (VMAT) (Stieler et al 2009) plans were the prostate and seminal vesicles. The plans were made with Monaco® (CMS, Elekta) TPS (versions 2.04.00, 3.20.01 and 3.30.00). The prescribed dose to the PTV was variable for the plans. The planning objectives were specified with an institutional template. A template for the specification of the planning goals for prostate cancer treatment in our institution is presented in table 1. The constraints of the OARs are appropriately scaled when the prescribed dose to the target is changed. The maximum number of rotations of the gantry per plan and control points per arc in the template are 2 and 120 respectively. These objectives are modified if the resultant plan is deemed to be undesirable. All the 128 plans were approved by both a physician and physicist for real patient treatment and are therefore considered as the best achievable compromise between the delivery of the prescribed dose to the PTV and the adequate sparing of the adjacent OARs. 2.3. Extraction of information from treatment plans
The plans were imported into MATLAB® (Mathworks Inc., Natick, MA, USA) computing environment using the Computational Environment for Radiotherapy Research (CERR) software (Deasy et al 2003). The following variables were extracted from each plan 1. The coordinates (x, y, z ) of each voxel of an OAR. The rectum and bladder are the two OARs that were considered in this work. This information was only used for visualization. 5579
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Table 1. A MONACO® planning template for prostate VMAT treatment at the
University Medical Centre Mannheim
Structure
Cost function
Reference dose (cGy)
Isoconstraint
P + SB
Target EUD Quadratic overdose Quadratic underdose
6000.0 6000.0
6000.0 70.0 70.0
Bladder
Parallel Serial
3000.0
40.0 3500.0
Rectum
Parallel Quadratic overdose
3000.0 3000.0
30.0 50.0
Patient
Quadratic overdose Quadratic overdose Conformality
3500.0 2200.0
120.0 60.0 0.7
P + SB = prostate and seminal vesicles
2. The reference dose, Di r, (normalized to the prescribed dose to the PTV) of the voxel. The subscript, i ( i = 1... n ), specifies the index of the voxel in the organ. The normalization of the dose renders the analysis to be independent of the absolute value of the prescribed dose to the PTV. The reference dose is the dose that was calculated for the voxel by MONACO® TPS. 3. The distance-to-PTV of the voxel (see the preceding section). Following the convention that was used in previous works (Wu et al 2009, Zhu et al 2011a, Yuan et al 2012), this value is signed negative if the voxel is shared by both the OAR and PTV. 4. The slice level of the voxel (see the preceding section). Two matrices, one for each OAR (rectum and bladder) were therefore computed from each plan (2 × 128 files). A row of the matrix represents a voxel while the columns specifies its parameters (coordinates, dose, distance-to-PTV and slice level). 2.4. The training set and the validation set
The 128 reference plans were randomly divided into two groups. The first group is henceforth referred to as the training set and consists of 95 plans, whereas the second group, which is henceforth referred to as the validation set is made up of 33 plans. For each OAR, the 95 matrices of the training set were merged to derive a single larger matrix of the organ. An organspecific model that relates the distances of voxels of the OARs to the PTV and path of the treatment beam, to their average dose was thereafter computed from this matrix. The plans of the validation set were used to assess the predictive accuracy of the derived model, and to evaluate the usefulness of this algorithm as a guiding tool for treatment planning. Each organspecific model is defined by two main functions (already described in the preceding section), namely the mean-dose-at-distance function and the slice weight function. Both functions have their associated standard deviation functions. The derivation of these functions is described in the following sections. 2.5. The mean dose-at-distance function (slice level = 0)
Figure 1 explains how the mean dose-at-distance function and the slice weight function are computed for each OAR from the training set. Consider a voxel as the multivariate object, (r , Z , D )i, where r is its distance-to-PTV, Z is the slice level and D is the reference dose (predicted by the 5580
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
TPS) of the voxel. The subscript i ( i = 1... N ) is the index of the voxel in the training set (after the merging of the 95 matrices). The mean dose-at-distance function of an organ is computed from the in-field voxels (i.e. having slice level = 0) of the OAR. The voxels were first binned (0.5 mm bin size according to their distance-to-PTV) value. The mean dose-at-distance value, Dj Z = 0 (also see preceding section), is thereafter computed for each distance-to-PTV bin. The formula for calculating this variable is given in figure 1. The mean-dose-at-distance function, f, is the function of the mean dose-at-distance values of the in-field voxels (the set containing all Dj Z = 0 values) against the centre of the corresponding distance-to-PTV bin (the set containing all Rj values). This function relates the mean dose of all the in-field voxels within a distance-toPTV bin to the median distance of the voxels to the surface of the PTV. The mean dose standard deviation function, g, is the function of the standard deviation Z=0 σ values against the centre of the corresponding distance( j ) of the mean dose-at-distance to-PTV bin. The variable, σj Z = 0 is calculated according to equation (2) σj Z = 0=
P
1 ∑(Dj Z = 0 − Di )2 P − 1 i=1
if ri ∈ [r j − 1, r j ) and Zi = 0
(2)
the variables i and P specifies the index of the in-field voxel in the jth distance-to-PTV bin and the total number of the in-field voxels in the bin respectively. 2.6. The slice weight function
We consider the out-of field voxels to be outside the path of the treatment beam (see the illustration in figure 1) in a coplanar beam arrangement, which is the case for prostate plans in our institution. We theorize that these voxels receive lower dose in optimized treatment plans when compared to the dose of the in-field voxels that are at the same distance to the PTV. We also hypothesize that the dose of the out-of-field voxels in a treatment plan can be determined by weighting the dose of the in-field voxels (if this is known) and vice versa. The slice weight function, p, specifies the weights that are required to scale the dose of a voxel according to its slice level. In contrast to the computation of the mean dose-at-distance function, all voxels (out-of field and in-field) of the training set of the organ of interest are used to compute the slice weight function according to the following steps. Firstly, the voxels of the training set (of the organ) were binned (0.5 mm bin size) according to their distance-to-PTV value. Voxels within a given distance-to-PTV bin were further separated according to their slice level. The mean dose-at-distance value Dj Z , is thereafter computed for each slice level (Z = 0.5, 1...Z max ) of a distance-to-PTV bin (j). This operation is performed for the entire distance-to-PTV bin. The mean dose-at-distance value of each slice level of a distance-to-PTV bin is then normalized to mean dose-at-distance value of the in-field voxels of that same bin. In other words Dj Z Dj Z = Z = 0 (3) Dj
A distance bin is ignored if it does not contain an in-field voxel since the described normalization cannot be performed. The slice weight factor, wZ, is calculated for each slice level as the average of the mean dose-at-distance values of that slice level in all the distance-to-PTV bins. The mathematical description of this operation (calculation of wZ) is provided in figure 1. The slice weight function, p, is the function of the slice weight factors against their corresponding slice levels. The slice-weight standard deviation function, q, specifies the standard deviation σZ, of the slice weight factors. The variable σZ is calculated according to equation (4) 5581
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
σZ =
M
1 ∑(wZ − Dj Z )2 M − 1 j=1
(4)
j(j = 0, 0.5…M ) is the index of the distance-to-PTV bin while Z (Z = 0.5, 1...Z max ) is the slice level. The mean-dose-at-distance function and the slice weight function of an OAR, together with their standard deviation functions (g and q) represent the model of learned relationship between the distance of a voxel of the organ to the PTV and treatment beam respectively, and the average dose of the voxel in optimized radiotherapy treatment plans. 2.7. Predicting the dose to a voxel from its position
The distance-to-PTV and the slice level (r , Z )i of a voxel are the two variables that are required to predict the dose of the voxel in an optimized treatment plan. Both variables are computed from the planning CT data. The subscript i(i = 1...n ) is the index of the voxel in the OAR. Firstly, the mean dose-at-distance value for the distance-to-PTV bin of the voxel is calculated from the mean dose-at-distance function according to equation (5) D j = f (ri ) (5) Z=0
Thereafter, standard deviation, σj Z = 0, of this value is computed from the mean dose standard deviation function, g, as shown in equation (6) σj Z = 0 = g (ri ) (6)
The unweighted dose, Ui, of an in-field voxel with a distance-to-PTV value equal to ri in an optimized treatment plan is estimated from Dj Z = 0 and σj Z = 0 according to equation (7) Z=0 U , σj Z = 0) (7) l = RNG1 (Dj
Function RNG1 randomly generates a number with a mean centered at Dj Z = 0 and a standard deviation that is defined by σj Z = 0 (RNG1 = normrnd() function in MATLAB). The unweighted dose is scaled if the voxel is an out-of-field voxel. The weighting factor depends on the slice level of the voxel. Firstly, the slice weight factor is calculated from the slice weight function according to equation (8) w (8) Z = p (Z i )
wZ and Zi are the slice weight factors and the slice level of the voxel (the input to the function) respectively. The standard deviation, σZ, of the slice weight factor is calculated from the slice weight standard deviation function q, as specified by equation (9) σ(9) Z = q (Z i )
The voxel-specific dose weight factor, wi, is calculated from the slice weight factor and its standard deviation according to equation (10) w (10) i = RNG2 (wZ , σZ )
The function RNG2 is similar to the previously described RNG1. The predicted dose, Di P, of the voxel is derived by applying the dose weight factor of the voxel to its unweighted dose as expressed in equation (11) DiP = Ui × wi (11) 5582
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
The value of the dose weight factor (wi ) is always equal to 1 for in-field voxels and by implication P U if Z i = 0 i = Di (12)
2.8. Plan evaluation metrics: the agreement scores
The comparison of the predicted and reference dose distributions starts at the voxel level. Firstly, the dose difference, Qi, between the reference, Di r, and predicted, Di P, dose to the ith voxel is calculated as the signed difference between both values P r Q i = Di −Di (13)
The organ agreement score, Sl organ, is the mean of the dose difference values of all the (n) voxels of the organ as described by equation (14) n
1 organ S(14) = ∑ Qi l n i=1
The subscript index l (l = 1, 2 ...V ) specifies the index of the organ in the patient. The sign on the organ agreement score indicates the direction of any difference; a negative score implies that the predicted dose to the organ is higher than the reference dose and vice versa. The plan agreement score, Splan, is calculated as the sum of the agreement scores for all the (V) OARs that are considered in the analysis V
organ S(15) plan = ∑ Sl l=1
2.9. Refinement of the training data and optimization of the organ models
The initial organ models that were generated from the training data (95 training data for each organ) were used as an ‘outlier detector’ to identify and remove suboptimal data from the training set of the organ. Refined models were thereafter generated from this modified training data. The outlier detection and removal process is as follows. The initial model was used to predict the dose to the individual plan data that make up the training set. The organ agreement score was thereafter computed for each plan. Since the predicted dose is subtracted from the reference dose in the computation of the organ agreement score (equations (13) and (14)), a positive score implies that the reference dose to the organ is higher than the dose predicted by the model. It should be borne in mind that the prediction of the model represents the average sparing pattern of the plans of the training set. A positive score (depending on the magnitude) therefore implies that the sparing achieved for the organ in the plan is lower than the group (training set) average. Our approach to optimizing the training data and refining the model is to remove plans whose organ agreement scores exceed a certain positive threshold value. 2.10. Validation of the model
The refined models were used to predict the likely 3D dose distribution in the OARs of the validation set based on their anatomies. The DVHs, organ agreement scores and the plan 5583
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
i
120
ii
1 Bladder Rectum
Bladder Rectum 0.9
100
0.7 80
Slice weight factor
Mean dose−at−distance [%]
0.8
60
40
0.6
0.5
0.4
0.3
0.2 20 0.1
0 −5
0
5
10
Distance−to−PTV [cm]
15
0 −2
0
2
4
Slice level [cm]
6
8
Figure 2. Subplot i shows the mean dose-at-distance functions of the rectum and blad-
der for prostate plan. The functions relate the mean dose that is received by voxels of these organs in optimized treatment plans for the prostate to the distances between the voxels and the PTV. The standard deviations of these mean values are plotted as error bars. The slice weight functions (subplot ii) relates the dose of voxels to the distances between the voxels a and the treatment beam direction. The standard deviations of the slice weight factors are plotted as error bars. These functions represent the organ-specific models for predicting the likely 3D dose distribution in these organs in optimized radiotherapy plans for the treatment of prostate cancer.
agreement scores were computed for the plans for the qualitative and quantitative comparison of the two (reference and predicted) datasets. 3. Results 3.1. Refined training set
The training set was refined by excluding plans whose organ agreement scores are >10%. The refined training set of the bladder and rectum contained 90 and 71 training data respectively. In other words, the final models of the bladder and rectum were generated from 90 and 71 plans respectively. 3.2. Relationship between position of a voxel and its mean dose (training set)
Figure 2 presents the relationship between the distance of voxels to the surface of the PTV and path of the treatment beam, and their average dose in optimized plans for the treatment of prostrate cancer. Subplot i shows that the mean dose-at-distance function of the bladder and rectum and their associated standard deviations (as error bars). The mean dose mean doseat-distance function of the bladder is higher than that of the rectum. This implies that for the treatment of prostate in our institution, a voxel of the bladder will receive a higher dose than a voxel of the rectum in an optimized plan if both voxels are at the same distance to the PTV. 5584
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Table 2. Comparison of the predicted and reference dose distribution by means of the
agreement scores. The organ agreement score quantifies the agreement between the predicted and reference 3D dose to an OAR. The plan agreement score is the sum of all the organ agreement scores of all the organs that are considered.
Plan index
Bladder agreement score (%)
Rectum agreement score (%)
Plan agreement score (%)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
−7.16 5.13 −6.30 4.16 8.03 1.77 −2.12 3.12 1.03 −1.15 4.03 −9.06 3.12 3.10 7.48 12.19 −1.58 8.19 0.91 4.51 4.22 1.35 3.83 −9.89 10.04 5.87 −2.36 −6.23 −2.82 6.37 6.36 −0.54 −0.26
6.83 4.45 3.00 −0.58 4.55 1.81 7.33 −1.35 0.93 2.05 2.42 −3.40 −1.21 −0.23 −1.63 4.36 2.75 −4.82 −4.11 4.28 1.94 −2.83 8.22 −1.54 −6.57 −0.35 3.39 −4.28 7.00 4.03 4.33 4.12 −2.34
−0.33 9.58 −3.30 3.59 12.57 3.58 5.21 1.77 1.97 0.90 6.45 −12.46 1.90 2.87 5.85 16.55 1.18 3.37 −3.20 8.79 6.16 −1.47 12.05 −11.43 3.46 5.52 1.03 −10.51 4.18 10.39 10.69 3.58 −2.61
Subplots ii present the slice weight function of the OARs and their standard deviations (as error bars). The function relates the weight of the dose of a voxel to its distance to the path of treatment beam. This means that if two voxels of an OAR are at the same distance to the PTV but are located on different slice levels, then given the dose of either voxel, the unknown dose can be estimated from the known dose and the slice weight function. For example, the slice weight factor for the bladder is 1 when the slice level is 0 (in-field voxel) and 0.2 (ignoring the standard deviation) when the slice level is 1.5 cm. This implies that if two voxels (a and b) of the bladder are at the same distance to the target, but with voxel a in the path of the treatment beam (Z = 0) and voxel b at a distance (slice level) of 1.5 cm away from the path of the beam, then the dose of both voxels are related by equation (16) (16) Db = Da × 0.2 Da and Db represents the unweighted doses of voxels a and b respectively (the unweighted dose of an in-field voxel is however equal its predicted dose as previously described in equation (12)). Polynomial approximations of these functions (figure 2) define the organ-specific models that were derived from the refined training set.
5585
O Nwankwo et al
0.5
0
0
50
Dose [%]
100
1.03
1
9
0.5
0
0
50
Dose [%]
100
3.12
1
13
0.5
0
0
50
Dose [%]
100
0
50
Dose [%]
1.77
1
6
0.5
0
0
50
Dose [%]
100
−1.15
1
10
0.5
0
0
50
Dose [%]
100
3.10
1
14
0.5
0
0
50
Dose [%]
100
0.5
0
0
50
Dose [%]
100
−2.12
1
7
0.5
0
0
50
Dose [%]
100
4.03
1
11
0.5
0
0
50
Dose [%]
100
7.48
1
15
0.5
0
Volume [Fraction]
Volume [Fraction]
100
3
Volume [Fraction]
5
Volume [Fraction]
8.03
1
0
−6.30
Volume [Fraction]
100
1
0
50
Dose [%]
Volume [Fraction]
50
Dose [%]
0.5
Volume [Fraction]
0
Volume [Fraction]
0
Predicted Reference
2
Volume [Fraction]
0.5
Bladder
5.13
1
Volume [Fraction]
1
Volume [Fraction]
−7.16
1
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Phys. Med. Biol. 59 (2014) 5575
100
4.16
1
4
0.5
0
0
50
Dose [%]
100
3.12
1
8
0.5
0
0
50
Dose [%]
100
−9.06
1
12
0.5
0
0
50
Dose [%]
100
12.19
1
16
0.5
0
0
50
Dose [%]
100
Figure 3. Validation of the organ model of the bladder using the plans of the validation set. The DVHs were computed from the predicted and reference 3D dose distributions of the plans. Based on the bladder agreement score, plan 9 showed the best agreement between the compared datasets while plan 16 showed the worst agreement.
3.3. Comparison of the predicted and reference datasets (validation set)
Table 2 quantitatively compares the predicted and reference dose distributions for the plans of the validation set by means of the organ agreement and the plan agreement scores, while figures 3 and 4 compares 16 plans of both datasets using their DVHs. Judging from the bladder agreement score (second column of table 2), plans 16 and 25 have worst agreement between the predicted and reference dose to the bladder, while plans 32 and 33 show the best agreement. The rectum agreement scores suggests that plans 4 and 14 have the best agreement between the predicted and reference datasets dose to the bladder. The plan agreement score suggests that Plan 1 has the best agreement between the predicted and reference datasets when the total dose to both organs are considered. The DVHs (figures 3 and 4) were computed from the predicted and reference 3D dose distributions in of OARs. Only the DVH of the first 16 plans of the validation set are presented for the sake of space. The integer values in the subplots specifies the index of the plan while the title of the plots show the agreement score for the organ. Figure 3 compares the predicted and reference DVHs of the bladder for 16 plans of the validation set. Based on the bladder agreement scores, plans 9 and 10 show the best agreement between the predicted and reference dose (values closest to zero), while plans 12 and 16 show the worst agreement (values farthest from zero). 5586
O Nwankwo et al
5
0.5
0
0
50
Dose [%]
100
0.93
1
9
0.5
0
0
50
Dose [%]
100
−1.21
1
13
0.5
0
0
50
Dose [%]
100
0
50
Dose [%]
1.81
1
6
0.5
0
0
50
Dose [%]
100
2.05
1
10
0.5
0
0
50
Dose [%]
100
−0.23
1
14
0.5
0
0
50
Dose [%]
100
0.5
0
0
50
Dose [%]
100
7.33
1
7
0.5
0
0
50
Dose [%]
100
2.42
1
11
0.5
0
0
50
Dose [%]
100
−1.63
1
15
0.5
0
Volume [Fraction]
Volume [Fraction]
100
3
Volume [Fraction]
1
Volume [Fraction]
4.55
0
3.00
Volume [Fraction]
100
1
0
50
Dose [%]
100
Volume [Fraction]
50
Dose [%]
0.5
Volume [Fraction]
0
Volume [Fraction]
0
2
Volume [Fraction]
0.5
Rectum
4.45
1
Volume [Fraction]
Predicted 1 Reference
Volume [Fraction]
6.83
1
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Volume [Fraction]
Phys. Med. Biol. 59 (2014) 5575
−0.58
1
4
0.5
0
0
50
Dose [%]
100
−1.35
1
8
0.5
0
0
50
Dose [%]
100
−3.40
1
12
0.5
0
0
50
Dose [%]
100
4.36
1
16
0.5
0
0
50
Dose [%]
100
Figure 4. Comparison of the reference and predicted DVHs of the rectum of the valida-
tion set. The organ agreement score (title of subplots) is a quantitative measure of the accuracy between the predicted and reference doses.
The quantitative scores agree with the visual ranking of the DVHs of these plans. This concordance reflects the accuracy of this metric in quantifying the level of agreement between the compared datasets. Figure 4 presents the comparison of the predicted and reference DVHs of the rectum. The rectum agreement score suggests that plan 14 has the best agreement between the predicted and reference data, while plans 1 and 7 show the least agreement. These results closely agree with the visual ranking of the agreement between the predicted and reference DVHs, thus providing further proof that the organ agreement score is a useful metric for quantifying the agreement between the predicted and reference dose distributions. A uniqueness of our approach relative to the reviewed methods is the ability of our algorithm to predict the achievable 3D dose distribution in the OARs of interest. Figure 5 shows the reference and predicted dose distributions for the OARs of Plan 8. The bladder and rectum agreement score are 3.12 and −1.35 respectively. The bladder (upper subfigures) and rectum (lower subfigures) of this plan consistes of 3.830 7 × 104 and 1.218 1 × 104 voxels respectively. Visual analysis (organ agreement scores and plan agreement score provides quantitative assessment) shows good agreement between both datasets. The predicted and
5587
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Figure 5. Visualization of the predicted and reference 3D dose distributions of the OARs of Plan 8. The reference dose distribution was calculated with the TPS while the predicted dose was calculated with our algorithm on the basis of the anatomical relationships between the organs and the target volume.
reference 3D dose distributions of all 33 plans of the validation set are provided as a supplementary data (stacks.iop.org/PMB/59/5575/mmedia). 4. Discussion Our aim in this work is to learn the relationship between the anatomies of organs at risk (OARs) and their dose in optimized radiotherapy treatment plans. This information was learned from a database of high quality previously treated patient plans. Thereafter, this knowledge is used to predict the likely dose of similar OARs on the basis of their anatomies. This prospective knowledge could be used for the anatomy-specific specification of treatment planning objectives with respect to normal tissue sparing, as well as a quality assurance (QA) tool to aid the objective assessment of the quality of radiotherapy plans on the basis of anatomy. The use of templates to specify planning objectives does not consider the individual patient anatomical variations. In addition, an objective decision on the quality of a radiotherapy treatment plan requires the quantitative knowledge of what can be achieved for that particular anatomy (Appenzoller et al 2012). Prospective knowledge of the achievable OAR sparing for a given anatomy is therefore an essential requirement for the implementation of an objective and individualized radiation therapy treatment-planning and acceptance scheme. We 5588
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
have described a method to learn the OAR sparing patterns from a database of high quality plans. This knowledge is then used to predict the likely dose that the organ will receive in an optimized treatment plan within a specific framework (similar planning conditions as the reference plans). Our approach is unique because it predicts the 3D dose in the OARs unlike previously described methods. The method that was used to derive the organ-specific models described in this work is similar to the technique that was used to create a virtual source model of a kilo-voltage radiotherapy device from a phase space file (Nwankwo et al 2013). Both are multivariate analyses methods, which operate by binning multivariate observations according to one variable (distance-to-PTV value of voxels in this work and the x-coordinates of photons in the cited work). Thereafter, the mean value of a second variable (dose in this work and direction cosine in the cited work) is computed for all the observations that fall within a bin, thus establishing a relationship between both variables (i.e. mean dose versus distance-to-PTV). The visual and quantitative comparison of the predicted and reference dose distributions show that this approach is accurate in predicting the 3D dose in the OARs, and that the evaluation metrics are accurate in quantifying the agreement between the compared datasets. It is possible to achieve a better than average sparing for an OAR at the expense of a competing organ through the adjustment of the prescribed dose to the OARs (increasing the penalty dose to one OAR while reducing or keeping constant the assigned penalty values of the competing OAR). This factor should be taken into account when comparing the predicted and reference dose distribution. It is for this reason that the plan agreement score is necessary. Plan 1 illustrates the usefulness of the plan agreement score metric. The bladder agreement score of Plan 1 is −7.16% while the rectum agreement score is 6.83% (table 2). This implies that the predicted to the bladder is higher than the reference dose, and that the reference dose to the rectum is higher than the prediction by our algorithm. This is suggestive of a prioritized sparing of the bladder at the expense of the rectum in the reference plan. This possible prioritized sparing of an OAR at the expense of another is accounted for by the plan agreement score. There are several potential applications of our algorithm in radiation therapy. One such application is the comparison of the OAR sparing of parallel treatment techniques. A previous study (Wolff et al 2009) used the mean dose to the anterior and posterior portions of the rectum, in addition to other metrics to compare the OAR sparing of parallel treatment techniques (VMAT, step and shoot IMRT, serial tomotherapy and 3D conformal radiotherapy). However, a restriction on using the mean organ dose and dose-to-volume in such analysis is that the same patient data should be used for the different techniques. This could be challenging if a large number of plans are to be analyzed. The mean-dose-at-distance function could be used to compare different treatment techniques even if the patient data is variable. This is possible because a superior treatment technique will deliver lower dose to the OARs (for the same PTV coverage) relative to an inferior technique. The difference in the achievable OAR sparing of the compared techniques will manifest in their respective organ-models (meandos-at-distance and slice weight functions). Such comparison does not require similar patient data because anatomy (distance-to-PTV) is already included in the model, thereby making it possible to objectively compare plans from different techniques, irrespective of whether or not the patient data is the same. Additionally and on the basis of the foregoing, this approach could be used to compare the treatment plan quality and or techniques across institutions. The approach that we have described in this work is not limited to the prostate and can be applied to other treatment sites. An important step in knowledge based radiation therapy treatment planning is the introduction of the overlap volume histogram (Wu et al 2009), distanceto-target-histogram and distance-to-PTV (Zhu et al 2011a). These authors demonstrated that 5589
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
the complex 3D OAR-target relationship could be reduced to the OVH/DTH of the organ, and further to the distance-to-PTV values of the voxels that make up the organ. Our approach of using the mean value of prior observation (mean dose of voxels at a given distance to PTV) to estimate the likelihood of future outcome can be applied to predict the likely dose of OARs at other treatment sites, for example in the head and neck region. Our future works will focus on the clinical validation and applications of our algorithm. It is possible to generate organ-models of desired properties by the iterative removal or addition of plans of certain characteristics to training set. The prediction of the model can be made increasing strict by removing plans with high positive organ agreement scores from the training set while keeping plans with low or negative agreement scores in the training set. Figure 2 provides information regarding the likely limitations of the current model. We will describe this likely limitation with the mean-dose-at-distance function of the bladder. The maximum distance-to-PTV value that the model was trained on is approximately 13.0 cm. Thus, the predictive accuracy of the model for voxels with distance-to-PTV values >13.0 cm is uncertain and depends on the behavior of the polynomial fit function beyond the available training data. A robust and generalized model requires the use of as large a training data as possible. The overprediction by the models may be the result of superior planning expertise, compromised coverage of the target to achieve a better than average sparing of the OARs, or the influence of features that were not captured by the model. 5. Conclusion We have described an algorithm to predict the likely 3D dose distribution in OARs before treatment planning. This prospective knowledge could be used to implement an objective and individualized quality assurance system for personalized radiation therapy treatment planning. Supplementary data Predicted and reference 3D dose distributions for all the plans of the validation set. Acknowledgment The authors wish to acknowledge the contribution of Ms Kerstin Siebenlist, physicists and physicians of University Medical Centre Mannheim, who either planned or reviewed the analyzed plans. Obioma Nwankwo also acknowledge to have been introduced to knowledge based radiation therapy treatment planning technique in 2010 by Q J Wu, X Zhu, and T Li of Duke University Medical Centre, Durham. References Appenzoller L M, Michalski J M, Thorstad W L, Mutic S and Moore K L 2012 Predicting dose-volume histograms for organs-at-risk in IMRT planning Med. Phys. 39 7446 Chanyavanich V, Das S K, Lee W R and Lo J Y 2011 Knowledge-based IMRT treatment planning for prostate cancer Med. Phys. 38 2515 Chung H T, Lee B, Park E, Lu J J and Xia P 2008 Can all centers plan intensity-modulated radiotherapy (IMRT) effectively? An external audit of dosimetric comparisons between 3D conformal radiotherapy and IMRT for adjuvant chemoradiation for gastric cancer Int. J. Radiat. Oncol. Biol. Phys. 71 1167–74 5590
O Nwankwo et al
Phys. Med. Biol. 59 (2014) 5575
Deasy J O, Blanco A I and Clark V H 2003 CERR: a computational environment for radiotherapy research Med. Phys. 30 979 Good D, Lo J, Lee W R, Wu Q J, Yin F-F and Das S K 2013 A knowledge-based approach to improving and homogenizing intensity modulated radiation therapy planning quality among treatment centers: an example application to prostate cancer planning Int. J. Radiat. Oncol. Biol. Phys. 87 176–81 Kazhdan M, McNutt T, Taylor R, Wu B and Simari P 2011 Comment on ‘A planning quality evaluation tool for prostate adaptive IMRT based on machine learning’ Med. Phys. 38 2820 Nwankwo O, Clausen S, Schneider F and Wenz F 2013 A virtual source model of a kilo-voltage radiotherapy device Phys. Med. Biol. 58 2363 Petit S F, Wu B, Kazhdan M, Dekker A, Simari P, Kumar R, Taylor R, Herman J M and McNutt T 2012 Increased organ sparing using shape-based treatment plan optimization for intensity modulated radiation therapy of pancreatic adenocarcinoma Radiother. Oncol. 102 38–44 Stieler F, Wolff D, Lohr F, Steil V, Abo-Madyan Y, Lorenz F, Wenz F and Mai S 2009 A fast radiotherapy paradigm for anal cancer with volumetric modulated arc therapy (VMAT) Radiat. Oncol. 4 25–48 Wang Y, Zolnay A, Incrocci L, Joosten H, McNutt T, Heijmen B and Petit S 2013 A quality control model that uses PTV-rectal distances to predict the lowest achievable rectum dose, improves IMRT planning for patients with prostate cancer Radiother. Oncol. 107 352–7 Wolff D, Stieler F, Welzel G, Lorenz F, Abo-Madyan Y, Mai S, Herskind C, Polednik M, Steil V and Wenz F 2009 Volumetric modulated arc therapy (VMAT) vs. serial tomotherapy, step-and-shoot IMRT and 3D-conformal RT for treatment of prostate cancer Radiother. Oncol. 93 226–33 Wu B, Ricchetti F, Sanguineti G, Kazhdan M, Simari P, Chuang M, Taylor R, Jacques R and McNutt T 2009 Patient geometry-driven information retrieval for IMRT treatment plan quality control Med. Phys. 36 5497 Yang Y, Ford E C, Wu B, Pinkawa M, van Triest B, Campbell P, Song D Y and McNutt T R 2013 An overlap-volume-histogram based method for rectal dose prediction and automated treatment planning in the external beam prostate radiotherapy following hydrogel injection Med. Phys. 40 011709 Yuan L, Ge Y, Lee W R, Yin F F, Kirkpatrick J P and Wu Q J 2012 Quantitative analysis of the factors which affect the interpatient organ-at-risk dose sparing variation in IMRT plans Med. Phys. 39 6868 Zhu X, Ge Y, Li T, Thongphiew D, Yin F-F and Wu Q J 2011a A planning quality evaluation tool for prostate adaptive IMRT based on machine learning Med. Phys. 38 719 Zhu X, Li T, Yin F-F, Wu Q J and Ge Y 2011b Response to ‘Comment on A planning quality evaluation tool for prostate adaptive IMRT based on machine learning’ Med. Phys. 38 2821
5591
