|Home | About | Journals | Submit | Contact Us | Français|
The aim of this work was to develop a framework for modeling organ effects within TOPAS (TOol for PArticle Simulation), a wrapper of the Geant4 Monte Carlo toolkit that facilitates particle therapy simulation. The DICOM interface for TOPAS was extended to permit contour input, used to assign voxels to organs. The following dose response models were implemented: The Lyman-Kutcher-Burman model, the critical element model, the population based critical volume model, the parallel-serial model, a sigmoid-based model of Niemierko for Normal Tissue Complication Probability (NTCP) and Tumor Control Probability (TCP), and a Poisson-based model for TCP. The framework allows easy manipulation of the parameters of these models and the implementation of other models.
As part of the verification, results for the parallel-serial and Poisson model for x-ray irradiation of a water phantom were compared to data from the AAPM Task Group 166. When using the task group dose-volume histograms (DVHs), results were found to be sensitive to the number of points in the DVH, with differences up to 2.4%, some of which are attributable to differences between the implemented models. New results are given with the point spacing specified. When using Monte Carlo calculations with TOPAS, despite the relatively good match to the published DVH’s, differences up to 9% were found for the parallel-serial model (for a maximum DVH difference of 2%) and up to 0.5% for the Poisson model (for a maximum DVH difference of 0.5%). However, differences of 74.5% (in Rectangle1), 34.8% (in PTV) and 52.1% (in Triangle) for the critical element, critical volume and the sigmoid-based models were found respectively.
We propose a new benchmark for verification of organ effect models in proton therapy. The benchmark consists of customized structures in the spread out Bragg peak (SOBP) plateau, normal tissue, tumor, penumbra and in the distal region. The DVH’s, DVH point spacing, and results of the organ effect models are provided. The models were used to calculate dose response for a Head and Neck patient to demonstrate functionality of the new framework and indicate the degree of variability between the models in proton therapy.
The superior accuracy of Monte Carlo simulations over analytical methods for dose calculation in radiotherapy is advantageous in proton therapy [Paganetti2012], perhaps more important than in x-ray therapy [TG105]. In either case, dose differences of 5% in radiotherapy can result in differences up to 10–20% in tumor control probability (TCP) or up to 20–30% in normal tissue complication probability (NTCP) [TG105]. The determination of the best treatment plan is aided through dose-effect modeling by evaluating the dose-volume histograms (DVH) for the target and nearby critical structures or organs. Commercial planning systems generally provide DVHs. If Monte Carlo simulation is not available in the planning system, the simulation may be performed separately. The subsequent step in organ effect modeling is to evaluate the therapeutic ratio using TCP and NTCP calculated from the DVHs or full dose distributions to determine the best plan [Armstrong2005]. Most often, the analysis of DVHs or the calculation of NTCP and TCP is performed using separate software for organ effect modeling such as DRESS [ElNaqa2006], CERR [Deasy2003], BIOPLAN [SanchezNieto2000], BIOSUITE [Uzan2009], HART [Pyakuryal2010] and others. None of these incorporate calculation of the dose distribution.
The TOPAS software was introduced in 2012 as an accurate, easy to use Monte Carlo tool for proton therapy simulation ([Perl2012] [Shin2012] [Schuemann2012]). The accuracy of TOPAS has been demonstrated using various experimental data [Testa2013]. Along with the accuracy, the availability of TOPAS to simulate time dependent calculations [Shin2012] makes it a powerful and easy to use Monte Carlo tool. The implementation of biological modeling in TOPAS is a step towards a full toolkit for patient specific dose and effect calculation in proton therapy.
The AAPM Task Group 166 [TG166] provides a benchmark for verification of the implementation of TCP/NTCP methods. A benchmark phantom configuration for TCP/NTCP calculation for photon therapy is described and reference data is supplied. However, due to the differences in dose distributions between photons and protons, this phantom configuration is not suitable for proton therapy and a new benchmark configuration is needed.
This work describes the addition of TCP/NTCP models to the TOPAS toolkit. A framework for adding dose response models to Monte Carlo simulation systems is described. A number of organ effect models were added to TOPAS using this framework, including many commonly used models. Models were validated against the TG-166 benchmark data for photon therapy and new information to improve the accuracy of the verification along with new data for models that were not benchmarked in TG-166 is supplied. A new benchmark based on the TG-166 photon benchmark is devised with data to support verification of organ models implemented by others. In order to show the new features of TOPAS for TCP/NTCP calculation, results of the various TCP/NTCP models implemented in TOPAS are compared for a head and neck patient with a realistic treatment plan using accurately commissioned proton fields.
TOPAS has the ability to read 3D information including CT data from DICOM files by means of the Grassroots DICOM Library (GDCM)[GDCM]. The DICOM reading functionality has been extended to read structure information. The dose distribution for each structure can be delineated and the corresponding DVH calculated. The dose at each organ is classified as follows. At initialization of TOPAS, for each organ, the coordinates that compose their corresponding contours are read from the DICOM-RT file, a winding number algorithm is used to determine whether a voxel (of the score grid) is enclosed by the polygon created with the coordinates of that contour. If it is enclosed, its location is storage in a Boolean grid with value set to True. At finalization, the corresponding DVH is calculated from the full dose distribution by assigning each voxel to a dose bin, to which it contributes with a count of the DVH [Sempau2000] and by using the Boolean grid as a filter for each organ considered.
The TOPAS modeling extension provides the capability to add dose response models to calculate the TCP and NTCP values using the full dose distribution calculated with TOPAS or DVHs (differential or cumulative) input separately. The values for the parameters of the TCP/NTCP models are defined by the user within the framework of the TOPAS parameter system [Perl2012]. Default values are provided. The user can override these values for any of the parameters using their own parameter files. The user defines the model and associate parameters in the input parameter file and TOPAS can calculate the TCP or NTCP from two scenarios. Scenario 1: If the dose calculation in an organ of interest is requested, the dose is calculated and TOPAS classifies those voxels contained in the organ. If the user does not request DVHs, TOPAS calculates the TCP or NTCP from the full dose distribution in the selected organ. In this case there is no re-binning in the dose axis of the DVH, with the number of bins in the DVH equal to the number of voxels in the corresponding organ with dose values larger than zero, while the voxel volume is set to unity. Otherwise, TOPAS uses the differential DVHs. Scenario 2: TOPAS reads external DVHs (in comma separated values format) or dose distributions for each organ of interest and calculates the NTCP or TCP values. If cumulative DVHs are imported they are converted to differential DVHs prior the dose response calculation. Summarizing, Scenario 1 uses the data produced by TOPAS at the end of a simulation (particle tracking plus TCP or NTCP calculation) whereas Scenario 2 uses input data to calculate TCP or NTCP (no particle tracking).
In addition to models that have been implemented in TOPAS, users may implement other models and may also include their own databases for the prebuilt models. The Figure 1 shows a sample of the code to implement a new model. TOPAS takes care of most of the underlying work so that the user who writes a new outcome modeling class just has to implement that which is unique about their own model.
In this section, we provide the equations of the main models we tested in order to avoid conflict with variations that even mathematically equivalent, computationally present certain differences. Several analytical methods are implemented in TOPAS to calculate NTCP or TCP, i.e. for NTCP, the Lyman-Kutcher-Burman [Burman1991], the parallel-serial or s-model [Kallman1992], the Critical Volume (population based) [Stavrev2001], the critical element [Niemierko1992] and the Niemierko model [Gay2007], and for TCP, a Poisson-based model [Warkentin2004] and the Niemierko model [Gay2007].
TD50 is the position of the 50% probability dose point, m is a parameter to control the slope of the dose response and gEUDa is the generalized uniform equivalent dose with a value a. vi is the fractional organ volume receiving a dose Di. To account for differences in dose per fraction, the dose Di can be substituted with the equivalent dose delivered in N fractions of 2 Gy [TG166] given by:
where α and β are the coefficients of the linear quadratic model.
The Critical Element model (CE) [Niemierko1991], the Critical Volume model [Niemierko1992] (a population-based variation is included in TOPAS and fully described in [Stavrev2001]), and the parallel-serial or S-model [Kallman1992] are included. These models are based on the idea that organs are composed of functional subunits. The CE model considers that these subunits are arranged in a serial way, while the CV model considers a parallel arrangement. The S-model is a generalization of the idea of functional subunits. Basically this model merges the CV and CE models into one single model. A value s determines the ‘degree’ of relatively seriality of the subunits. For s equal to 0, this model is equivalent to the CV model, while for s equal to unity it is equivalent to the CE model. Nevertheless, values outside of the range 0 and 1 are often found [Savreva2002]. The NTCP [Kallman1992] is then given by:
where P(Di) is calculate by means of the Poisson model written as [Kallman1992]
where γ is the normalized dose-response gradient at the dose, where the absolute dose-response gradient is the steepest.
Another NTCP model included is the proposal of Niemierko for parameterization of dose response using the logistic function for NTCP or TCP [Gay2007]. The model depends on the gEUDa value described in equation 3. The value for the tolerance of dose determines the kind of calculation that is performed [Gay2007]. The NTCP and TCP are given by:
where TCD50 states for tumor control dose at 50%.
This TCP model is based on the Poisson probability distribution as described in [Kallman1992]. In the TOPAS implementation, the probability per voxel is calculated with equation 10 [Warkentin2004], and the TCP is calculated as the product over all voxels as:
For LKB, CV and Poisson models, predefined databases with parameters corresponding to certain organs or tumor sites are implemented in TOPAS. The data were obtained from references [Burman1991] [Stavrev2001] [Okunieff1995] for LKB, CV and Poisson models, respectively. In addition, users may revise these model parameters by including customized TOPAS parameter values in their own parameter files, which will override the predefined values of the databases.
To verify the implementation of the models in TOPAS, the NTCP and TCP values calculated with TOPAS were compared to the corresponding values calculated with Pinnacle system version 9.2 (Philips Medical Systems, Andover, MA) for the benchmark water phantom and structures configuration for photon beams in the AAPM TG-166 task group report [TG166]. For this work a water phantom was created consisting of 100 × 100 × 100 voxels of 5 × 5 × 2.5 mm3 dimension with 4 basic structures (Table VI in TG-166). The phantom was created using Pinnacle and exported to a DICOM-RT file read by TOPAS for dose calculation and subsequent TCP/NTCP calculation. The geometry with structures is shown in Figure 2. The beam source consisted of accurate phase space files from a Siemens 6 MV x-ray beam modeled with the BEAMnrc software [Rogers1995] using accurately determined source and geometry details [Sawkey2009]. This is the same beam commissioned at UCSF in Pinnacle on the 4 UCSF Siemens Primus linacs, matched to within 1% out to the penumbra. The field was set to 20 × 20 cm2 at the 100 cm source to surface distance. The dose calculation was performed with TOPAS up to a statistical uncertainty (calculated using the method recommended in [Rogers2000]) lower than 0.7% for voxels with dose exceeding 50% of the maximum. The production cut of secondary particles was set to 0.05 mm with a max step length of 0.5 mm. The source phase space was recycled 10 times. The final dose distribution was normalized to a dose of 72 Gy in the 20 × 20 × 0.5 cm3 region located at 6 cm depth along the central axis of the water phantom. The dose normalization region encompassed several voxels to minimize statistical fluctuation in this region.
The TCP or NTCP values for the 4 structures (for organs listed in the Table VI in [TG-166]) were calculated as follows:
Summarizing, item 2 was used to estimate the variations in NTCP values when different number of bins in the DVH is used. Item 3 was used to compare the NTCP TOPAS calculations by using the same DVH as Pinnacle.
The previous section described a benchmark configuration for TCP/NTCP suggested by the TG166. This setup was defined for photon radiotherapy. Because of differences between dose conformality between photon therapy and proton therapy this benchmarking (mainly the particle source) of dose response models using the phantom from report [TG166] is not suitable for proton therapy. To our knowledge, a setup to benchmark dose response models for proton therapy is currently not available. In this section we propose a benchmark phantom for proton therapy.
As in [TG166] the phantom implemented in TOPAS consist of a water phantom but with different structures configuration. Four structures were considered to account for the dose at proximal, target, distal, and penumbra regions. The parameters of the structures to be used in the models were chosen to avoid 0 % and 100% values in NTCP or TCP calculations. The detailed configuration is given in the Figure 3 and Table 1. The water phantom has 20 × 20 × 20 cm3 dimension with voxel resolution of 0.5 × 0.5 × 2 mm3. A squared proton field of 10 cm in x and y, 10 cm of range and 5 cm of modulation impinges on the surface of the phantom from the negative z-axis. The source model was calculated based on the Francis H Burr Proton Therapy treatment head [Paganetti2004]. To reduce the computation time at the treatment head, a geometrical particle split was used [Ramos2012]. Phase space files were scored at the end of the treatment head and subsequently used for the dose calculation. In total, 72 million source protons were simulated to achieve a statistical uncertainty below 1% at the flat region of the spread-out Bragg peak (SOBP).
The final dose distribution was re-scaled to obtain a dose of 78 Gy in 39 fractions at the center of the SOBP (7.5 cm from the entrance beam in the water phantom on the central z-axis).
Figure 4 shows the cumulative DVHs calculated with TOPAS and Pinnacle, the latter calculated with the adaptive convolution option. The DVHs of the four structures considered were included. Differences of about 0.5% between both systems were found for the structures PTV, Rectangle1 and Rectangle2 and up to 2% for the structure Triangle. These differences were mainly due to the algorithms used for dose calculation. For Monte Carlo calculations the statistical uncertainties in dose distributions were below the 0.7% for voxels with values larger than the 50% of the maximum dose. Further, a maximum difference of 1% was found between the DVHs from [TG166] and those calculated in this work with Pinnacle for the structures PTV, Rectangle1, Rectangle2 and Triangle.
The TCP and NTCP values are shown in the Table 2. The TCP values were calculated for the structures PTV and Rectangle1 only, as in TG-166. Further, the gEUDa value for a=1 is included in the table. ‘LKB with Eq. 4’ denotes the dose fractionation correction prior to evaluating the dose response; ‘LKB without Eq. 4’ denotes that dose fractionation was not considered in the dose response evaluation.
The NTCP of the dose distributions calculated with the adaptive convolution algorithm in Pinnacle was within 0.5% of the TG-166 result, in agreement with the TG-166 comparison to values calculated with their spreadsheet. The LKB without Eq. 4 matched better with the Pinnacle LKB model, suggesting the model in Pinnacle is incomplete, lacking the correction. This underscores the importance of verification testing as discussed in TG-166.
There is a wide range of values from the different dose response models. For the LKB model, a difference of up to 50% between the NTCP values calculated with S-model and LKB has been observed [TG166], consistent with the results for Rectangle2.
The TCP/NTCP results for the PTV and rectangular structures calculated with the models implemented in TOPAS and Pinnacle agree within 2 standard deviations of the statistical precision of the TOPAS calculation or 1% of the TCP/NTCP result. However, for the triangle, there are significant differences between TOPAS and Pinnacle for the LKB without Eq. 4 (a 9% difference) and S-Model (13%) outside of the statistical precision of the TOPAS dose calculation. These differences are attributed to the (unknown) point spacing in the DVH used in Pinnacle and to the difference in the dose distributions calculated with TOPAS and Pinnacle. For example, differences up to 4% in NTCP with respect to the TG166 data were found for the three algorithms (adaptive convolution, fast convolution and collapse cone convolution) available in the Pinnacle
Results from TOPAS in Table 2 used dose distributions calculated with TOPAS. To eliminate the dependence of the verification testing on the Monte Carlo calculation, further TCP/NTCP calculations were done with the DVHs from Pinnacle imported directly into TOPAS.
We used an arbitrary dose point spacing of 29.7 cGy (corresponding to 200 points) for comparison purposes, which was not necessarily the spacing used internally by Pinnacle. Table 3 shows the differences in percent between Pinnacle and TOPAS in this case. For this particular study, the difference could be mainly due to a difference in the number of points in the DVH. Cozzi et al [Cozzi2000] have shown that a refinement as small as 2% in the DVH point spacing leads to a difference up to 10% in TCP/NTCP values. For the TCP model, our implementation of the Poisson model agreed within 0.5% with the reference given in [TG166]. For the S-Model in Table 3, the maximum difference of 3% between the Pinnacle and TOPAS results (5.4% corresponding to 100 points) was obtained for the structure Rectangle2. This difference was also in part due to the Poisson distribution P(Di) in Equation 6, which differs mathematically from the distribution implemented in Pinnacle [TG166]. On the other hand, the Niemierko, CE and CV models were compared against the S-Model from Pinnacle, and differences up to 74.5% (Critical Element model in Rectangle1) were found. The differences are largely attributable to the sensitivity of the model results to DVH point spacing. This shows that the benchmark results in [TG166] should specify the DVH point spacing, or preferably not use rebinning (see section 2.1).
The cumulative DVHs (bin width of 0.01 Gy) for the proposed benchmark phantom are shown in Figure 5. The corresponding statistical uncertainties were below 0.8% for those voxels with dose values larger than the 10% of the maximum dose. In Table 4, the NTCP, TCP and gEUDa=1 values are shown. The statistical uncertainties are only due to Monte Carlo. The values for the LKB without Eq 4 and gEUDa=1 without Eq. 4 are included for completeness. The S-model was taken as reference to validate the LKB model with Eq. 4, the CV model, the CE model, and the Niemierko model. This is a convenient reference as the S-model is used in TG-166 to validate biological modeling implemented in commercial planning systems. Interestingly, the differences largely fall within the statistical precision of the TOPAS calculation. The largest differences were 2.0%±4.6% for LKB with Eq. 4, 1.6%±4.9% for Niemierko model, 3.7%±1.6% for CV model and 12.05%±0.38% for CE Model. Although the models evaluated in this work have been generally restricted to photon therapy, with limited development of proton specific models, they should be useful for protons. The validation for photons in this study applies equally to protons as there is no distinction between protons and photons for these models, other than the model parameters, listed in the Table 1 and in references [Burman1991] and [Stavrev2001]. Thus the results for the proposed benchmark geometry serve as validation data.
This section shows a short example of how to use the TCP/NTCP extension implemented in TOPAS. In this example, the LKB model is calculated for the organ Brain stem. The user inputs the DICOM and dose scoring information along with the parameters of the model as shown in Figure 6. The figure also shows the cumulative DVHs for selected organs of a Head and Neck patient treated with protons as calculated with TOPAS. The corresponding parameters for the LKB and CV models can be found in [Burman1991] and [Stavrev2001], respectively. Table 5 shows the parameters for the S-Model, CE, and Niemierko models. The prescribed dose was 60 Gy at the target delivered in 30 fractions. Table 6 shows the NTCP and TCP values (corrected for fractionation assuming 2 Gy per fraction and α/β of 3 Gy) for the corresponding organs. The TCP values were calculated only for GTV and CTV regions. As shown, the critical element and the S-Model agreed for the brain stem NTCP value, because the seriality of the organ is equal to 1. The missing values in the cells are because the parameter values for the Critical Volume model for optical nerve are not defined in the literature.
In this work, commonly used NTCP and TCP models were implemented in TOPAS. The framework developed for this purpose is suitable for implementation of other models.
In [TG166] the dose distribution was calculated with analytical algorithms available in several systems (Pinnacle system included). We compared TOPAS results to the Pinnacle system. NTCP results for the LKB model, S-model and a Poisson-based model agreed with the benchmark data provided in [TG166] for x-ray therapy. The largest differences are due to differences in the dose calculation algorithm and the number of points in the DVHs. The Pinnacle result reported in TG-166 was likely calculated with the adaptive convolution algorithm, although this was not stated in the publication, and the LKB model implemented in Pinnacle version 9.2 appears to lack a correction for fractionation. When comparing the Critical Element, Critical Volume and Niemierko models, maximum differences of 74.5% (in Rectangle1), 34.8% (in PTV) and 52.1% (in Triangle) were found respectively. This is attributed to differences in the DVH point spacing, which was not given in [TG166]. The point spacing is explicit in our results.
We provided a benchmark configuration suitable for reporting validation data for proton beams for NTCP and TCP calculation. The structure was configured to account for 4 regions of interest in proton therapy: at plateau, at maximum, at penumbra and at range straggling region of a spread out Bragg peak. We provided new reference NTCP/TCP values for all 6 models implemented in TOPAS. Interestingly, the differences between the S-Model and either the LKB model, Critical Volume and Niemierko models, largely fall within the statistical precision of the TOPAS calculation. The largest differences were 1.97%±4.6% for LKB with Eq. 4, 1.57%±4.9% for Niemierko model, 3.74%±1.6% for CV model and 12.05%±0.38% for CE Model.
Finally, for the head and neck patient, the framework implemented in TOPAS allows the simultaneous comparison of all implemented models in a practical scenario, with different criteria for evaluating the TCP/NTCP. In this particular example, the Critical Element, Niemierko and S-Model showed similar results for the Brain stem organ (7–12%).
We would like to thank Dr. Josephine Chen at the University of California San Francisco for help with the Pinnacle treatment planning system. This work was supported by National Cancer Institute Grant R01CA140735