Home
Search
Collections
Journals
About
Contact us
My IOPscience
Parameterization of brachytherapy source phase space file for Monte Carlo-based clinical brachytherapy dose calculation
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 Phys. Med. Biol. 59 455 (http://iopscience.iop.org/0031-9155/59/2/455) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 128.97.90.221 This content was downloaded on 02/05/2014 at 21:38
Please note that terms and conditions apply.
Institute of Physics and Engineering in Medicine Phys. Med. Biol. 59 (2014) 455–464
Physics in Medicine and Biology
doi:10.1088/0031-9155/59/2/455
Parameterization of brachytherapy source phase space file for Monte Carlo-based clinical brachytherapy dose calculation M Zhang, W Zou, T Chen, L Kim, A Khan, B Haffty and N J Yue Department of Radiation Oncology, Rutgers Cancer Institute of New Jersey, Rutgers University, New Brunswick, NJ 08901, USA E-mail: [email protected] Received 23 August 2013, revised 28 October 2013 Accepted for publication 4 November 2013 Published 30 December 2013 Abstract
A common approach to implementing the Monte Carlo method for the calculation of brachytherapy radiation dose deposition is to use a phase space file containing information on particles emitted from a brachytherapy source. However, the loading of the phase space file during the dose calculation consumes a large amount of computer random access memory, imposing a higher requirement for computer hardware. In this study, we propose a method to parameterize the information (e.g., particle location, direction and energy) stored in the phase space file by using several probability distributions. This method was implemented for dose calculations of a commercial Ir-192 high dose rate source. Dose calculation accuracy of the parameterized source was compared to the results observed using the full phase space file in a simple water phantom and in a clinical breast cancer case. The results showed the parameterized source at a size of 200 kB was as accurate as the phase space file represented source of 1.1 GB. By using the parameterized source representation, a compact Monte Carlo job can be designed, which allows an easy setup for parallel computing in brachytherapy planning. (Some figures may appear in colour only in the online journal) 1. Introduction Accurately calculating the radiation dose for patients is critical in brachytherapy treatment planning. Currently, the standard clinical brachytherapy dosimetry protocol is based on Task Group report 43 (TG-43), published by the American Association of Physicists in Medicine (AAPM) (Nath et al 1995, Rivard et al 2004), which is designed to determine dose distribution in an infinitely sized uniform water phantom. A patient anatomy with brachytherapy applicators can be approximated by an infinitely sized water phantom if the following two conditions are 0031-9155/14/020455+10$33.00
© 2014 Institute of Physics and Engineering in Medicine Printed in the UK & the USA
455
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
met: (1) there are no large heterogeneities inside the region of interest, and (2) there is adequate tissue surrounding the region of interest to provide the same scattering condition found in an infinitely sized water phantom. However, clinically, these conditions are rarely met. For example, there is no tissue beyond the skin in breast cancer brachytherapy to provide the scattering condition. The skin dose calculation following the TG-43 protocol would assume a full scattering condition and would very likely overestimate the true skin dose. Investigations (Serago et al 1991, Mangold et al 2001, Wallner et al 2004, Kassas et al 2006, Raffi et al 2010) have shown that the overestimation could be as large as 15%. Evidently, more accurate dose calculation methods are required in brachytherapy; especially for cases in which heterogeneous tissues and materials are involved. By simulating individual interactions between the impinging particles and the material, the Monte Carlo method can accurately determine the radiation dose in heterogeneous materials. The Monte Carlo method has been used to generate the TG-43 type dosimetric parameters, such as the radial dose function, the anisotropy function, and the dose rate constant (Williamson et al 2005, DeWerd et al 2011, Dezarn et al 2011, Chiu-Tsao et al 2012). However, several challenges exist in the implementation of the Monte Carlo method in clinical treatment planning. One of the challenges is the simulation of the radiation source. In general, there are two approaches: (1) construction of the source geometry, including its encapsulation, materials and structure, and (2) utilization of a phase space file. In the first approach, it is often difficult to accurately and precisely define the complicated source geometry, e.g., radioisotope core and encapsulation, along with the voxelized patient geometry used for dose scoring. Furthermore, with this approach, the dose calculation process can be time consuming. Recently, Enger et al (2012) proposed a way to construct the radiation source in a parallel space to the space containing the actual patient geometry using the Geant4 Monte Carlo toolkit (Agostinelli et al 2003). Although it is a viable solution, it requires a customized source code that is not publicly available. The second approach requires the user to first prepare a phase space file containing information on all the particles emitted from a radiation source. Then, in the subsequent patient dose calculation, the particles stored in the phase space file would be randomly sampled and transported. On one hand, it is a desirable solution since particles from the source are generated once and can be used in future calculations to significantly reduce computing time. On the other hand, the phase space file needs to contain information on a large amount of particles in order to yield reliable results, leading to a file size of a several gigabytes (GB). To expedite calculation, the phase space file needs to be loaded into computer RAM prior to the dose calculation. The large size of the phase space file often becomes a limiting factor for the practical implementation of this approach. In this paper, we proposed a method to reduce the size of the phase space file by parameterizing the location, direction and energy distribution of particles stored in the phase space file. The parameterization was evaluated and implemented for the phase space file of a high dose rate (HDR) brachytherapy source. The accuracy and efficiency of this method was tested by comparing the dosimetric parameters calculated in both a water phantom and a clinical brachytherapy case between a phase space file represented radiation source and its parameterized counterpart. 2. Methods and materials 2.1. Phase space file preparation
The Geant4 Monte Carlo toolkit (Agostinelli et al 2003) was used throughout this study to generate the phase space file and related dose calculations for the Varian Ir-192 HDR source 456
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 1. The sketch (left) and the schematic diagram (right) of the Varian VariSource brachytherapy radiation source.
(VariSource, Varian Medical Systems, Palo Alto, CA). The detailed source geometry was constructed following the published literature (Rasmussen et al 2011) as shown in figure 1. The Ir-192 isotope was assumed to be uniformly distributed inside the core with a radius of 0.17 mm and length of 5 mm. The G4RadioactiveDecay model in Geant4 was used to simulate the decay of Ir-192. The phase space file scored all the photons emitted from the source surface, covering 8 mm from the tip of the source, as shown in figure 1. Based on our simulation, the number of photons emitted beyond the scoring region at the tail of the source was less than 0.2% of all the photons emitted from the active Ir-192 core, and therefore ignored. Any photon passing the scoring region was stored in the phase space file and was represented with seven coordinates: three for its location, three for its direction and one for its energy. In total, 2 × 107 Ir-192 decay events were simulated. It resulted in a phase space file of 1.1 GB containing 4.2 × 107 photons. 2.2. Phase space file parameterization
The particle’s location was described using the Cartesian coordinate system with the z-axis along the source central axis (figure 1). The positive z direction pointed from the cable to the tip of the source. The x- and y-axis were on the transverse plane of the source. The particle’s direction was described by a spherical coordinate system. The polar angle (ϕ) was the angle between the particle’s direction and the positive z-axis. The azimuthal angle (θ ) was the angle between the particle direction’s x–y projection and the positive x-axis. To parameterize the phase space file, we used several probability distributions generated from the phase space file to characterize the particles stored in the file. First of all, we assumed that the energy spectrum of the particles emerging from the source was universal, regardless of its location and direction. Proof of this will be provided in the Results section. The energy spectrum (E-distribution) with a bin size of 5 keV has been generated from the phase space file. It was known that there would be a different number of particles emerging from the surface of the source, with more at the center of the source and less toward both ends. Therefore, a Z-distribution was generated, which characterized the number of particles emerging from the surface of the source at different z, with a bin size of 0.1 mm along the z direction covering the entire 8 mm length of the source. Because the source geometry has rotational symmetry along the z-axis, the outgoing particles also follow this rotational symmetry. Therefore, uniform 457
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
distributions can be used to characterize the particles’ location and direction on the x-y plane. Finally, the only coordinate left uncharacterized was the polar angle of the particles’ direction. Due to the different amount of attenuation from the source encapsulation along the z-axis, the polar angle distribution would vary with z. For the same reason, it would also be different for different energy groups. Two Phi-distributions respectively describing the polar angular distribution of the high energy group (E > 250 keV) and the low energy group (E < = 250 keV) were established for each z location with a bin size of 0.1 mm. The process of generating E-distribution, Z-distribution and Phi-distribution from the phase space file was the parameterization. These three distributions were used instead of the phase space file to generate source particles in Monte Carlo simulations. The following procedures were carried out to generate a source particle from these distributions during a simulation: (1) The energy of a particle was randomly sampled from the E-distribution. (2) The location of a particle along the source direction (z) was randomly sampled from the Z-distribution. (3) The azimuthal angle (θ ) of a particle’s direction was independent of the other six coordinates which could be sampled uniformly from 0 to 2π . (4) The location of a particle in the transverse plane of the source (x–y) depended on the azimuthal angle (θ ) of the particle’s direction. To ensure the particle location being on the surface of the source, its location in the transverse plane was assumed to be uniformly distributed on the arc centered at (0 mm, 0 mm) with the radius of the source (0.17 mm) ranging from θ –π /2 to θ + π /2. (5) The polar angle (ϕ) of a particle’s direction was randomly sampled from the Phidistribution after z and E were known. 2.3. Phantom evaluations
Two test simulation phantoms were used to evaluate the accuracy of the parameterized source: a simple simulated water phantom and a clinical breast cancer patient’s CT scan with placed MammoSite applicator (Proxima Therapeutics Inc., Alpharetta, GA). We compared the Monte Carlo calculation results between the phase space file represented source (S_phase) and the parameterized source (S_para). The simulated water phantom was 10 × 10 × 10 cm3 in size. The radial dose function and anisotropy function as defined in TG-43 were calculated with the brachytherapy source simulated in the center of the phantom. The line-source approximation was used to calculate the geometry function. In total, 4 × 108 particles were sampled for each calculation. This demonstrated that the statistical uncertainty of the calculated dose was less than 1% within 5 cm of the source on the central transverse plan, and progressively increased toward the tip of the source. The CT images of a clinical breast cancer patient treated with MammoSite applicatorbased HDR brachytherapy was also used for the comparison of dose calculation results (figure 2). An in-house developed software code (Zhang et al 2011) was used to convert patient CT scans into a data format suitable for Geant4. The target (PTV_EVAL) was contoured following the Radiation Therapy Oncology Group protocol (RTOG 2005). Skin and rib were contoured as the critical organs. Skin volume was defined as 3 mm expansion inside the body surface. The clinical plan parameters, including source dwell locations and dwell times, were imported from the treatment planning system. For this particular plan, there were 22 dwell locations, and 5 × 107 particles were generated for each dwell location calculation. The total dose was calculated as the sum of doses contributed from each dwell location with statistical 458
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 2. A breast cancer patient with MammoSite applicator. PTV_EVAL (red), skin
(green), and rib (yellow) were contoured for dose evaluation.
uncertainty of less than 1% within 5 cm of the source, and increasing toward the tip of the source. 3. Results and discussion Figure 3 shows the spectrum of photons emitted from three locations along the source: the tip, the middle and the tail. There were minimal differences between the energy spectrums, which supports our assumption that the energy spectrum of a photon can be treated as one universal distribution, regardless of its location and direction. Figure 4 shows the E-distribution and Z-distribution for the tested HDR source. Figure 5 shows the Phi-distribution of particles at three different locations along the z-axis: the tip, the middle and the tail of the source. These three distributions utilized 220 kB of memory space during each calculation. This was a huge memory reduction compared to the 1.1 GB used by the S_phase represented source model. Although there were additional computations required to generate a source particle from the parameterized source representation (S_para), e.g., sampling several distributions, we did not find any computation speed difference between those two source representations. Figure 6(a) shows the comparison of the computed radial dose function between S_phase and S_para. The difference was smaller than 2% within a 50 mm range. Figures 6(b)–(e) presents the comparison of the anisotropy function. The discrepancy between two source representations increased with theta close to 0◦ and 180◦ , which corresponded to the areas close to the tip and the tail of the source. Figure 7 shows the local percentage difference of the anisotropy function between S_phase and S_para. As illustrated in the figure, the difference was less than 5%, except in areas within a 10◦ angle of the source z-axis. The possible cause 459
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 3. Energy spectrum of photons emitted from the tip, the middle and the tail of
the source. (a)
(b)
Figure 4. E-distribution (a) and Z-distribution (b) of the tested HDR source.
of the larger discrepancy between the two source representations at the tip and the tail of the source may be due to the fact that fewer particles were scored in the phase space file at those locations. With fewer particles, the parameterized Phi-distribution would have a greater amount of uncertainties, which may misrepresent the real scenario. The variance reduction technique should be used in future studies to generate a better representation of particles emitted from the tip and the tail of the source. Figure 8 shows the isodose distributions in the target region computed with S_phase and S_para, respectively, in the MammoSite treatment. The relative isodose values were normalized to the corresponding maximum dose. In most areas, dose distributions computed from S_para correlated very well to those generated from S_phase. The only area that showed a relatively large dose difference between S_para and S_phase was at the tip of the applicator, as indicated by the arrow in figure 8. This result was expected, since the largest difference between S_phase and S_para in the water phantom was at the tip of the source within a 10◦ angle of the source axis (figure 7). Dose volume histograms (DVHs) of PTV_EVAL, skin and rib are presented in figure 9. Following clinical conventions, the Monte Carlo calculated raw dose value was normalized so that 90% of the PTV_EVAL was covered by 90% of the prescription dose. Table 1 further summarizes the differences in the clinically relevant DVH 460
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 5. Phi-distribution of high energy group (E > 250 keV) and low energy group
(E < 250 keV) particles emitted from the tip, the middle and the tail of the tested HDR source.
(a)
(c)
(b)
(d)
(e)
Figure 6. Comparison of the radial dose function (a) and anisotropy function at r =
5 mm (b), 10 mm (c), 20 mm (d), and 50 mm (e) between S_phase and S_para.
points of these structures. It was observed that the differences between the dose distributions calculated from S_phase and S_para were less than 1%, except with the maximum dose on the skin, which showed a difference of 1.59%. In this work, we proposed a method to parameterize the phase space file to reduce the memory consumption of the Monte Carlo calculation. It has been successfully implemented for the VariSource Ir-192 brachytherapy source. Testing the efficacy of this method on different brachytherapy sources remains an interesting research topic for the future. The major bottleneck to implementing the Monte Carlo method in clinical applications is the unacceptably 461
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 7. The local percentage difference of the anisotropy function between S_phase
and S_para.
Figure 8. The isodose distribution in the target region from the S_phase (red) and S_para
(green). Table 1. Selected dose-volume parameters for S_phase and S_para calculated dose distribution. The dose was normalized to have D90%(PTV_EVAL) = 90% of the prescription in S_phase calculated dose distribution.
Dmax (rib) Dmax (skin) D90% (PTV_EVAL) V150% (PTV_EVAL) V200% (PTV_EVAL)
S_phase
S_para
Abs. Diff.
156.48% 149.86% 90.00% 24.22% 7.13%
157.26% 148.27% 90.02% 24.82% 7.41%
−0.78% 1.59% −0.02% −0.6% −0.28%
long computing time. Reducing memory consumption by using the parameterized source representation would not directly accelerate the computation itself. However, it allows easy setup of parallel computation on a computer network for the Monte Carlo dose calculations of 462
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 9. DVH curves of PTV_EVAL, skin and rib from the MammoSite breast cancer
treatment calculated from S_phase and S_para.
brachytherapy plans with multiple source dwell positions. The Monte Carlo calculations can be naturally parallelized based on individual source dwell position, and their computational times could be significantly reduced. Using the parameterized source presentation, each individual computational job would be compact and easily loadable onto a calculation node. It minimizes the possibility of calculation node memory overflow by an otherwise large phase space file. Recently, the Monte Carlo method has been developed on graphic processing unit (GPU) architecture (Jia et al 2012). The limited onboard memory on the GPU card compared to CPU random access memory would limit the file size that can be loaded for a given job. By using parameterized source configuration instead of the phase space file represented source, multiple jobs can be run in parallel on GPU. In that sense, the parameterized source approach overcomes the computer memory limitation and allows the implementation of parallel Monte Carlo brachytherapy dose calculation over a computer network or a GPU to speed up the computation. 4. Conclusion In this work, we proposed a method to use several probability distributions to parameterize information stored in a phase space file of a brachytherapy source. The methodology was evaluated on an HDR source. Compared to the phase space file represented source, the parameterized source representation significantly reduces the memory consumption without sacrificing speed and accuracy. It may potentially simplify the parallel computing setup process on a computer network for brachytherapy plan dose calculations. Acknowledgment None of the authors have any potential conflict of interest to the work presented in this paper. References Agostinelli S et al 2003 GEANT4-a simulation toolkit Nucl. Instrum. Methods Phys. Res. A 506 250–303 Chiu-Tsao S et al 2012 Dosimetry of 125I and 103Pd COMS eye plaques for intraocular tumors: Report of Task Group 129 by the AAPM and ABS Med. Phys. 39 6161–84 463
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
DeWerd L et al 2011 A dosimetric uncertainty analysis for photon-emitting brachytherapy sources: Report of AAPM Task Group No. 138 and GEC-ESTRO Med. Phys. 38 782–93 Dezarn W A et al 2011 Recommendations of the American Association of Physicists in Medicine on dosimetry, imaging, and quality assurance procedures for 90Y microsphere brachytherapy in the treatment of hepatic malignancies Med. Phys. 38 4824–36 Enger S A et al 2012 Layered mass geometry: a novel technique to overlay seeds and applicators onto patient geometry in Geant4 brachytheapy simulations Phys. Med. Biol. 57 6269–77 Jia X et al 2012 A GPU tool for efficient, accurate, and realistic simulation of cone beam CT projections Med. Phys. 39 7368–78 Kassas B et al 2006 Dose modification factors for 192Ir high-dose-rate irradiation using Monte Carlo simulations J. Appl. Clin. Med. Phys. 7 28–34 Mangold C A et al 2001 Quality control in interstitial brachytherapy of the breast using pulsed dose rate: treatment planning and dose delivery with an Ir-192 afterloading system Radiother. Oncol. 58 43–51 Nath R et al 1995 Dosimetry of interstitial brachytherapy sources: recommendations of the AAPM Radiation Therapy Committee Task Group No. 43 Med. Phys. 22 209–34 Radiation Therapy Oncology Group (RTOG) 2005 A randomized phase III study of conventional whole breast irradiation (WBI) versus partial breast irradiation (PBI) for women with Stage 0, I, or II breast cancer NSABP Protocol B-39, RTOG Protocol 0413 (Philadelphia, PA: RTOG) Raffi J A et al 2010 Determination of exit skin dose for 192Ir intracavitary accelerated partial breast irradiation with thermoluminescent dosimeters Med. Phys. 37 2693–702 Rasmussen et al 2011 Comparison of air-kerma strength determinations for HDR 192Ir sources Med. Phys. 38 6721–9 Rivard M J et al 2004 Update of AAPM Task Group No. 43 Report: a revised AAPM protocol for brachytherapy dose calculations Med. Phys. 31 633–74 Serago C F et al 1991 Scattering effects on the dosimetry of iridium-192 Med. Phys. 18 1266–70 Wallner P et al 2004 Workshop on partial breast irradiation: state of the art and the science Bethesda, MD, December 8–10, 2002 J. Natl Cancer Inst. 96 175–84 Williamson J et al 2005 Recommendations of the American Association of Physicists in Medicine regarding the impact of implementing the 2004 task group 43 report on dose specification for 103Pd and 125I interstitial brachytherapy Med. Phys. 32 1424–40 Zhang M et al 2011 Introducing an on-line adaptive procedure for prostate image guided intensity modulate proton therapy Phys. Med. Biol. 56 4947–65
464
Search
Collections
Journals
About
Contact us
My IOPscience
Parameterization of brachytherapy source phase space file for Monte Carlo-based clinical brachytherapy dose calculation
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 Phys. Med. Biol. 59 455 (http://iopscience.iop.org/0031-9155/59/2/455) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 128.97.90.221 This content was downloaded on 02/05/2014 at 21:38
Please note that terms and conditions apply.
Institute of Physics and Engineering in Medicine Phys. Med. Biol. 59 (2014) 455–464
Physics in Medicine and Biology
doi:10.1088/0031-9155/59/2/455
Parameterization of brachytherapy source phase space file for Monte Carlo-based clinical brachytherapy dose calculation M Zhang, W Zou, T Chen, L Kim, A Khan, B Haffty and N J Yue Department of Radiation Oncology, Rutgers Cancer Institute of New Jersey, Rutgers University, New Brunswick, NJ 08901, USA E-mail: [email protected] Received 23 August 2013, revised 28 October 2013 Accepted for publication 4 November 2013 Published 30 December 2013 Abstract
A common approach to implementing the Monte Carlo method for the calculation of brachytherapy radiation dose deposition is to use a phase space file containing information on particles emitted from a brachytherapy source. However, the loading of the phase space file during the dose calculation consumes a large amount of computer random access memory, imposing a higher requirement for computer hardware. In this study, we propose a method to parameterize the information (e.g., particle location, direction and energy) stored in the phase space file by using several probability distributions. This method was implemented for dose calculations of a commercial Ir-192 high dose rate source. Dose calculation accuracy of the parameterized source was compared to the results observed using the full phase space file in a simple water phantom and in a clinical breast cancer case. The results showed the parameterized source at a size of 200 kB was as accurate as the phase space file represented source of 1.1 GB. By using the parameterized source representation, a compact Monte Carlo job can be designed, which allows an easy setup for parallel computing in brachytherapy planning. (Some figures may appear in colour only in the online journal) 1. Introduction Accurately calculating the radiation dose for patients is critical in brachytherapy treatment planning. Currently, the standard clinical brachytherapy dosimetry protocol is based on Task Group report 43 (TG-43), published by the American Association of Physicists in Medicine (AAPM) (Nath et al 1995, Rivard et al 2004), which is designed to determine dose distribution in an infinitely sized uniform water phantom. A patient anatomy with brachytherapy applicators can be approximated by an infinitely sized water phantom if the following two conditions are 0031-9155/14/020455+10$33.00
© 2014 Institute of Physics and Engineering in Medicine Printed in the UK & the USA
455
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
met: (1) there are no large heterogeneities inside the region of interest, and (2) there is adequate tissue surrounding the region of interest to provide the same scattering condition found in an infinitely sized water phantom. However, clinically, these conditions are rarely met. For example, there is no tissue beyond the skin in breast cancer brachytherapy to provide the scattering condition. The skin dose calculation following the TG-43 protocol would assume a full scattering condition and would very likely overestimate the true skin dose. Investigations (Serago et al 1991, Mangold et al 2001, Wallner et al 2004, Kassas et al 2006, Raffi et al 2010) have shown that the overestimation could be as large as 15%. Evidently, more accurate dose calculation methods are required in brachytherapy; especially for cases in which heterogeneous tissues and materials are involved. By simulating individual interactions between the impinging particles and the material, the Monte Carlo method can accurately determine the radiation dose in heterogeneous materials. The Monte Carlo method has been used to generate the TG-43 type dosimetric parameters, such as the radial dose function, the anisotropy function, and the dose rate constant (Williamson et al 2005, DeWerd et al 2011, Dezarn et al 2011, Chiu-Tsao et al 2012). However, several challenges exist in the implementation of the Monte Carlo method in clinical treatment planning. One of the challenges is the simulation of the radiation source. In general, there are two approaches: (1) construction of the source geometry, including its encapsulation, materials and structure, and (2) utilization of a phase space file. In the first approach, it is often difficult to accurately and precisely define the complicated source geometry, e.g., radioisotope core and encapsulation, along with the voxelized patient geometry used for dose scoring. Furthermore, with this approach, the dose calculation process can be time consuming. Recently, Enger et al (2012) proposed a way to construct the radiation source in a parallel space to the space containing the actual patient geometry using the Geant4 Monte Carlo toolkit (Agostinelli et al 2003). Although it is a viable solution, it requires a customized source code that is not publicly available. The second approach requires the user to first prepare a phase space file containing information on all the particles emitted from a radiation source. Then, in the subsequent patient dose calculation, the particles stored in the phase space file would be randomly sampled and transported. On one hand, it is a desirable solution since particles from the source are generated once and can be used in future calculations to significantly reduce computing time. On the other hand, the phase space file needs to contain information on a large amount of particles in order to yield reliable results, leading to a file size of a several gigabytes (GB). To expedite calculation, the phase space file needs to be loaded into computer RAM prior to the dose calculation. The large size of the phase space file often becomes a limiting factor for the practical implementation of this approach. In this paper, we proposed a method to reduce the size of the phase space file by parameterizing the location, direction and energy distribution of particles stored in the phase space file. The parameterization was evaluated and implemented for the phase space file of a high dose rate (HDR) brachytherapy source. The accuracy and efficiency of this method was tested by comparing the dosimetric parameters calculated in both a water phantom and a clinical brachytherapy case between a phase space file represented radiation source and its parameterized counterpart. 2. Methods and materials 2.1. Phase space file preparation
The Geant4 Monte Carlo toolkit (Agostinelli et al 2003) was used throughout this study to generate the phase space file and related dose calculations for the Varian Ir-192 HDR source 456
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 1. The sketch (left) and the schematic diagram (right) of the Varian VariSource brachytherapy radiation source.
(VariSource, Varian Medical Systems, Palo Alto, CA). The detailed source geometry was constructed following the published literature (Rasmussen et al 2011) as shown in figure 1. The Ir-192 isotope was assumed to be uniformly distributed inside the core with a radius of 0.17 mm and length of 5 mm. The G4RadioactiveDecay model in Geant4 was used to simulate the decay of Ir-192. The phase space file scored all the photons emitted from the source surface, covering 8 mm from the tip of the source, as shown in figure 1. Based on our simulation, the number of photons emitted beyond the scoring region at the tail of the source was less than 0.2% of all the photons emitted from the active Ir-192 core, and therefore ignored. Any photon passing the scoring region was stored in the phase space file and was represented with seven coordinates: three for its location, three for its direction and one for its energy. In total, 2 × 107 Ir-192 decay events were simulated. It resulted in a phase space file of 1.1 GB containing 4.2 × 107 photons. 2.2. Phase space file parameterization
The particle’s location was described using the Cartesian coordinate system with the z-axis along the source central axis (figure 1). The positive z direction pointed from the cable to the tip of the source. The x- and y-axis were on the transverse plane of the source. The particle’s direction was described by a spherical coordinate system. The polar angle (ϕ) was the angle between the particle’s direction and the positive z-axis. The azimuthal angle (θ ) was the angle between the particle direction’s x–y projection and the positive x-axis. To parameterize the phase space file, we used several probability distributions generated from the phase space file to characterize the particles stored in the file. First of all, we assumed that the energy spectrum of the particles emerging from the source was universal, regardless of its location and direction. Proof of this will be provided in the Results section. The energy spectrum (E-distribution) with a bin size of 5 keV has been generated from the phase space file. It was known that there would be a different number of particles emerging from the surface of the source, with more at the center of the source and less toward both ends. Therefore, a Z-distribution was generated, which characterized the number of particles emerging from the surface of the source at different z, with a bin size of 0.1 mm along the z direction covering the entire 8 mm length of the source. Because the source geometry has rotational symmetry along the z-axis, the outgoing particles also follow this rotational symmetry. Therefore, uniform 457
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
distributions can be used to characterize the particles’ location and direction on the x-y plane. Finally, the only coordinate left uncharacterized was the polar angle of the particles’ direction. Due to the different amount of attenuation from the source encapsulation along the z-axis, the polar angle distribution would vary with z. For the same reason, it would also be different for different energy groups. Two Phi-distributions respectively describing the polar angular distribution of the high energy group (E > 250 keV) and the low energy group (E < = 250 keV) were established for each z location with a bin size of 0.1 mm. The process of generating E-distribution, Z-distribution and Phi-distribution from the phase space file was the parameterization. These three distributions were used instead of the phase space file to generate source particles in Monte Carlo simulations. The following procedures were carried out to generate a source particle from these distributions during a simulation: (1) The energy of a particle was randomly sampled from the E-distribution. (2) The location of a particle along the source direction (z) was randomly sampled from the Z-distribution. (3) The azimuthal angle (θ ) of a particle’s direction was independent of the other six coordinates which could be sampled uniformly from 0 to 2π . (4) The location of a particle in the transverse plane of the source (x–y) depended on the azimuthal angle (θ ) of the particle’s direction. To ensure the particle location being on the surface of the source, its location in the transverse plane was assumed to be uniformly distributed on the arc centered at (0 mm, 0 mm) with the radius of the source (0.17 mm) ranging from θ –π /2 to θ + π /2. (5) The polar angle (ϕ) of a particle’s direction was randomly sampled from the Phidistribution after z and E were known. 2.3. Phantom evaluations
Two test simulation phantoms were used to evaluate the accuracy of the parameterized source: a simple simulated water phantom and a clinical breast cancer patient’s CT scan with placed MammoSite applicator (Proxima Therapeutics Inc., Alpharetta, GA). We compared the Monte Carlo calculation results between the phase space file represented source (S_phase) and the parameterized source (S_para). The simulated water phantom was 10 × 10 × 10 cm3 in size. The radial dose function and anisotropy function as defined in TG-43 were calculated with the brachytherapy source simulated in the center of the phantom. The line-source approximation was used to calculate the geometry function. In total, 4 × 108 particles were sampled for each calculation. This demonstrated that the statistical uncertainty of the calculated dose was less than 1% within 5 cm of the source on the central transverse plan, and progressively increased toward the tip of the source. The CT images of a clinical breast cancer patient treated with MammoSite applicatorbased HDR brachytherapy was also used for the comparison of dose calculation results (figure 2). An in-house developed software code (Zhang et al 2011) was used to convert patient CT scans into a data format suitable for Geant4. The target (PTV_EVAL) was contoured following the Radiation Therapy Oncology Group protocol (RTOG 2005). Skin and rib were contoured as the critical organs. Skin volume was defined as 3 mm expansion inside the body surface. The clinical plan parameters, including source dwell locations and dwell times, were imported from the treatment planning system. For this particular plan, there were 22 dwell locations, and 5 × 107 particles were generated for each dwell location calculation. The total dose was calculated as the sum of doses contributed from each dwell location with statistical 458
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 2. A breast cancer patient with MammoSite applicator. PTV_EVAL (red), skin
(green), and rib (yellow) were contoured for dose evaluation.
uncertainty of less than 1% within 5 cm of the source, and increasing toward the tip of the source. 3. Results and discussion Figure 3 shows the spectrum of photons emitted from three locations along the source: the tip, the middle and the tail. There were minimal differences between the energy spectrums, which supports our assumption that the energy spectrum of a photon can be treated as one universal distribution, regardless of its location and direction. Figure 4 shows the E-distribution and Z-distribution for the tested HDR source. Figure 5 shows the Phi-distribution of particles at three different locations along the z-axis: the tip, the middle and the tail of the source. These three distributions utilized 220 kB of memory space during each calculation. This was a huge memory reduction compared to the 1.1 GB used by the S_phase represented source model. Although there were additional computations required to generate a source particle from the parameterized source representation (S_para), e.g., sampling several distributions, we did not find any computation speed difference between those two source representations. Figure 6(a) shows the comparison of the computed radial dose function between S_phase and S_para. The difference was smaller than 2% within a 50 mm range. Figures 6(b)–(e) presents the comparison of the anisotropy function. The discrepancy between two source representations increased with theta close to 0◦ and 180◦ , which corresponded to the areas close to the tip and the tail of the source. Figure 7 shows the local percentage difference of the anisotropy function between S_phase and S_para. As illustrated in the figure, the difference was less than 5%, except in areas within a 10◦ angle of the source z-axis. The possible cause 459
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 3. Energy spectrum of photons emitted from the tip, the middle and the tail of
the source. (a)
(b)
Figure 4. E-distribution (a) and Z-distribution (b) of the tested HDR source.
of the larger discrepancy between the two source representations at the tip and the tail of the source may be due to the fact that fewer particles were scored in the phase space file at those locations. With fewer particles, the parameterized Phi-distribution would have a greater amount of uncertainties, which may misrepresent the real scenario. The variance reduction technique should be used in future studies to generate a better representation of particles emitted from the tip and the tail of the source. Figure 8 shows the isodose distributions in the target region computed with S_phase and S_para, respectively, in the MammoSite treatment. The relative isodose values were normalized to the corresponding maximum dose. In most areas, dose distributions computed from S_para correlated very well to those generated from S_phase. The only area that showed a relatively large dose difference between S_para and S_phase was at the tip of the applicator, as indicated by the arrow in figure 8. This result was expected, since the largest difference between S_phase and S_para in the water phantom was at the tip of the source within a 10◦ angle of the source axis (figure 7). Dose volume histograms (DVHs) of PTV_EVAL, skin and rib are presented in figure 9. Following clinical conventions, the Monte Carlo calculated raw dose value was normalized so that 90% of the PTV_EVAL was covered by 90% of the prescription dose. Table 1 further summarizes the differences in the clinically relevant DVH 460
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 5. Phi-distribution of high energy group (E > 250 keV) and low energy group
(E < 250 keV) particles emitted from the tip, the middle and the tail of the tested HDR source.
(a)
(c)
(b)
(d)
(e)
Figure 6. Comparison of the radial dose function (a) and anisotropy function at r =
5 mm (b), 10 mm (c), 20 mm (d), and 50 mm (e) between S_phase and S_para.
points of these structures. It was observed that the differences between the dose distributions calculated from S_phase and S_para were less than 1%, except with the maximum dose on the skin, which showed a difference of 1.59%. In this work, we proposed a method to parameterize the phase space file to reduce the memory consumption of the Monte Carlo calculation. It has been successfully implemented for the VariSource Ir-192 brachytherapy source. Testing the efficacy of this method on different brachytherapy sources remains an interesting research topic for the future. The major bottleneck to implementing the Monte Carlo method in clinical applications is the unacceptably 461
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 7. The local percentage difference of the anisotropy function between S_phase
and S_para.
Figure 8. The isodose distribution in the target region from the S_phase (red) and S_para
(green). Table 1. Selected dose-volume parameters for S_phase and S_para calculated dose distribution. The dose was normalized to have D90%(PTV_EVAL) = 90% of the prescription in S_phase calculated dose distribution.
Dmax (rib) Dmax (skin) D90% (PTV_EVAL) V150% (PTV_EVAL) V200% (PTV_EVAL)
S_phase
S_para
Abs. Diff.
156.48% 149.86% 90.00% 24.22% 7.13%
157.26% 148.27% 90.02% 24.82% 7.41%
−0.78% 1.59% −0.02% −0.6% −0.28%
long computing time. Reducing memory consumption by using the parameterized source representation would not directly accelerate the computation itself. However, it allows easy setup of parallel computation on a computer network for the Monte Carlo dose calculations of 462
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
Figure 9. DVH curves of PTV_EVAL, skin and rib from the MammoSite breast cancer
treatment calculated from S_phase and S_para.
brachytherapy plans with multiple source dwell positions. The Monte Carlo calculations can be naturally parallelized based on individual source dwell position, and their computational times could be significantly reduced. Using the parameterized source presentation, each individual computational job would be compact and easily loadable onto a calculation node. It minimizes the possibility of calculation node memory overflow by an otherwise large phase space file. Recently, the Monte Carlo method has been developed on graphic processing unit (GPU) architecture (Jia et al 2012). The limited onboard memory on the GPU card compared to CPU random access memory would limit the file size that can be loaded for a given job. By using parameterized source configuration instead of the phase space file represented source, multiple jobs can be run in parallel on GPU. In that sense, the parameterized source approach overcomes the computer memory limitation and allows the implementation of parallel Monte Carlo brachytherapy dose calculation over a computer network or a GPU to speed up the computation. 4. Conclusion In this work, we proposed a method to use several probability distributions to parameterize information stored in a phase space file of a brachytherapy source. The methodology was evaluated on an HDR source. Compared to the phase space file represented source, the parameterized source representation significantly reduces the memory consumption without sacrificing speed and accuracy. It may potentially simplify the parallel computing setup process on a computer network for brachytherapy plan dose calculations. Acknowledgment None of the authors have any potential conflict of interest to the work presented in this paper. References Agostinelli S et al 2003 GEANT4-a simulation toolkit Nucl. Instrum. Methods Phys. Res. A 506 250–303 Chiu-Tsao S et al 2012 Dosimetry of 125I and 103Pd COMS eye plaques for intraocular tumors: Report of Task Group 129 by the AAPM and ABS Med. Phys. 39 6161–84 463
M Zhang et al
Phys. Med. Biol. 59 (2014) 455
DeWerd L et al 2011 A dosimetric uncertainty analysis for photon-emitting brachytherapy sources: Report of AAPM Task Group No. 138 and GEC-ESTRO Med. Phys. 38 782–93 Dezarn W A et al 2011 Recommendations of the American Association of Physicists in Medicine on dosimetry, imaging, and quality assurance procedures for 90Y microsphere brachytherapy in the treatment of hepatic malignancies Med. Phys. 38 4824–36 Enger S A et al 2012 Layered mass geometry: a novel technique to overlay seeds and applicators onto patient geometry in Geant4 brachytheapy simulations Phys. Med. Biol. 57 6269–77 Jia X et al 2012 A GPU tool for efficient, accurate, and realistic simulation of cone beam CT projections Med. Phys. 39 7368–78 Kassas B et al 2006 Dose modification factors for 192Ir high-dose-rate irradiation using Monte Carlo simulations J. Appl. Clin. Med. Phys. 7 28–34 Mangold C A et al 2001 Quality control in interstitial brachytherapy of the breast using pulsed dose rate: treatment planning and dose delivery with an Ir-192 afterloading system Radiother. Oncol. 58 43–51 Nath R et al 1995 Dosimetry of interstitial brachytherapy sources: recommendations of the AAPM Radiation Therapy Committee Task Group No. 43 Med. Phys. 22 209–34 Radiation Therapy Oncology Group (RTOG) 2005 A randomized phase III study of conventional whole breast irradiation (WBI) versus partial breast irradiation (PBI) for women with Stage 0, I, or II breast cancer NSABP Protocol B-39, RTOG Protocol 0413 (Philadelphia, PA: RTOG) Raffi J A et al 2010 Determination of exit skin dose for 192Ir intracavitary accelerated partial breast irradiation with thermoluminescent dosimeters Med. Phys. 37 2693–702 Rasmussen et al 2011 Comparison of air-kerma strength determinations for HDR 192Ir sources Med. Phys. 38 6721–9 Rivard M J et al 2004 Update of AAPM Task Group No. 43 Report: a revised AAPM protocol for brachytherapy dose calculations Med. Phys. 31 633–74 Serago C F et al 1991 Scattering effects on the dosimetry of iridium-192 Med. Phys. 18 1266–70 Wallner P et al 2004 Workshop on partial breast irradiation: state of the art and the science Bethesda, MD, December 8–10, 2002 J. Natl Cancer Inst. 96 175–84 Williamson J et al 2005 Recommendations of the American Association of Physicists in Medicine regarding the impact of implementing the 2004 task group 43 report on dose specification for 103Pd and 125I interstitial brachytherapy Med. Phys. 32 1424–40 Zhang M et al 2011 Introducing an on-line adaptive procedure for prostate image guided intensity modulate proton therapy Phys. Med. Biol. 56 4947–65
464
