- Email: [email protected]
Monte Carlo dose calculations for high-dose–rate brachytherapy using GPU-accelerated processing
Brachytherapy 15 (2016) 387e398
Physics
Monte Carlo dose calculations for high-doseerate brachytherapy using GPU-accelerated processing Z. Tian1,*, M. Zhang2, B. Hrycushko1, K. Albuquerque1, S.B. Jiang1, X. Jia1,* 1
Department of Radiation Oncology, University of Texas Southwestern Medical Center, Dallas, TX Department of Radiation Oncology, Rutgers Cancer Institute of New Jersey, New Brunswick, NJ
2
ABSTRACT
PURPOSE: Current clinical brachytherapy dose calculations are typically based on the Association of American Physicists in Medicine Task Group report 43 (TG-43) guidelines, which approximate patient geometry as an infinitely large water phantom. This ignores patient and applicator geometries and heterogeneities, causing dosimetric errors. Although Monte Carlo (MC) dose calculation is commonly recognized as the most accurate method, its associated long computational time is a major bottleneck for routine clinical applications. This article presents our recent developments of a fast MC dose calculation package for high-doseerate (HDR) brachytherapy, gBMC, built on a graphics processing unit (GPU) platform. METHODS AND MATERIALS: gBMC-simulated photon transport in voxelized geometry with physics in 192Ir HDR brachytherapy energy range considered. A phase-space file was used as a source model. GPU-based parallel computation was used to simultaneously transport multiple photons, one on a GPU thread. We validated gBMC by comparing the dose calculation results in water with that computed TG-43. We also studied heterogeneous phantom cases and a patient case and compared gBMC results with Acuros BV results. RESULTS: Radial dose function in water calculated by gBMC showed !0.6% relative difference from that of the TG-43 data. Difference in anisotropy function was !1%. In two heterogeneous slab phantoms and one shielded cylinder applicator case, average dose discrepancy between gBMC and Acuros BV was !0.87%. For a tandem and ovoid patient case, good agreement between gBMC and Acruos BV results was observed in both isodose lines and dose-volume histograms. In terms of the efficiency, it took ~47.5 seconds for gBMC to reach 0.15% statistical uncertainty within the 5% isodose line for the patient case. CONCLUSIONS: The accuracy and efficiency of a new GPU-based MC dose calculation package, gBMC, for HDR brachytherapy make it attractive for clinical applications. 2016 American Brachytherapy Society. Published by Elsevier Inc. All rights reserved.
Keywords:
High-doseerate brachytherapy; Monte Carlo; Dose calculation; GPU
Introduction High-doseerate (HDR) brachytherapy has proved to be a highly successful treatment for cancers in a variety of disease sites (1), such as cervical cancer, breast cancer, prostate cancer, etc. By placing an applicator in the targeted area and
Received 9 November 2015; received in revised form 26 January 2016; accepted 27 January 2016. * Corresponding authors. Department of Radiation Oncology, University of Texas Southwestern Medical Center, Dallas, TX 75390. Tel.: 214-648-5031; fax: 214-645-2885. E-mail addresses: [email protected] (Z. Tian), [email protected] (X. Jia).
having a motor-controlled radioactive source dwell in a series of planned positions with preset time, HDR brachytherapy delivers a highly conformed dose to the target while sparing nearby critical organs. Accurate dose calculation in patients is a critical step in treatment planning for HDR brachytherapy to ensure good target coverage and sharp dose gradients for critical structure sparing. The standard dosimetry protocol currently used in clinical brachytherapy is based on the Association of American Physicists in Medicine Task Group report 43 (TG-43) (2, 3). This method approximates a three-dimensional (3D) dose distribution by superposition of interpolated and precalculated discrete dose matrices in an infinitely large water medium. Although this approach is straightforward and fast, the influence of material
1538-4721/$ - see front matter 2016 American Brachytherapy Society. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.brachy.2016.01.006
388
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
heterogeneity is ignored (4), causing dosimetric errors. For instance, the skin dose has been shown to be overestimated by as large as 15% (5e8). In breast brachytherapy, the TG43 approach tends to underestimate dose to lung by ~10% (9). In addition, neglecting dosimetric impacts of nontissue heterogeneities, such as applicator material and the shielding material, is another concern. It has been shown that ignoring the shielding of an HDR colorectal applicator leads to an overestimation of the target volume within 100% isodose lines by O5% (10), which implies potential tumor cold spots. Therefore, it is proposed in Task Group report 186 (4) to move toward more accurate model-based dose calculation algorithms for brachytherapy, such as collapsed-cone superposition/convolution methods (11, 12), deterministic solutions to the linear Boltzmann transport equation (13, 14), and Monte Carlo (MC) simulations (15). MC has been considered as one of the most accurate dose calculation method for determining radiation dose deposition in heterogeneous materials. Particle transport and interactions with the medium are simulated based on fundamental physics principles (4,15e17). However, MC simulations have been associated with extremely long calculation time, impeding applications in clinical practice. An accelerated CPU-based MC dose calculator (HDRMC) was recently developed (18) which used a fast rayetracing technique to transport photons through a 3D mesh of voxels representing the patient anatomy. It took ~5 min to calculate a typical HDR brachytherapy treatment plan consisting of 20e30 dwell positions. A track length estimator was used to improve the efficiency of the MC simulation (19). This track length estimator requires calculation of mass energy absorption coefficient before the simulation run, and the absorption energies used in the simulation must be identical to that used to determine the mass energy absorption (20). Recently, there have been tremendous research efforts in accelerating MC simulations using graphics processing unit (GPU) platforms (21e30). For HDR brachytherapy dose calculations, bGPUMCD (31) was developed by extending the photon transport part of the GPUMCD system (23) into the brachytherapy energy range. By using GPU’s rapid parallel processing power, GPU-friendly parallelization schemes, and a track length estimator, the dose calculation can be completed in seconds. This MC package was also extended to low-doseerate brachytherapy calculations (32,33). Grid-based Boltzmann equation solvers are another category of model-based dose calculation methods for brachytherapy. By directly solving the linear Boltzmann transport equation, the dose distribution is derived from fundamental physics principles and, hence, achieves guaranteed accuracy (14). This algorithm has been introduced for brachytherapy dose calculations in a clinic setting as Acuros BV (Varian Medical Systems, Palo Alto, CA). However, the calculation time for a treatment plan typically ranges in several minutes on a multicore CPU platform. This speed may be acceptable for plan dose calculation but becomes
burdensome when repeated calculations are needed, such as in inverse plan optimization. In this report, we will present our recent developments of a GPU-based MC dose engine for HDR brachytherapy, named gBMC (GPU-based brachytherapy Monte Carlo). We will show high computational efficiency without using variance reduction techniques except the Woodcock tracking method, a method that has been widely used in GPU-based photon transport simulation tools (34). To validate our gBMC package, we have benchmarked against the TG-43 results in a water phantom and against Acuros BV results in a heterogeneous slab phantom and shielded cylinder applicators. The feasibility of applying our gBMC package for clinical practice will also be shown in a typical HDR brachytherapy patient case.
Methods and materials Source modeling and physics The gBMC package was developed on top of our established GPU-based MC packages, gDPM for dose calculation for external radiotherapy (22, 24) and gCTD for imaging dose calculation of CT and cone-beam CT (35), both of which have been validated against the conventional CPU-based MC codes. The gBMC package currently focuses on 192Ir, the most commonly used radionuclide for HDR brachytherapy. The source model in gBMC was a phase-space file generated from the MC simulation of radioactive source decay and particle transport for a VariSource VS2000 192Ir HDR source (Varian Medical Systems) using GEANT4 (36). The source was modeled as a radioactive core inside a connecting cable with geometry and material composition following published data (37,38). For the 2 107 decay events sampled uniformly throughout the radioactive core, about 4.2 107 photons exited the surface of the cable, which were stored in the phase-space file. The energy spectrum of these recorded photons is shown in Fig. 1. In our gBMC package, we adopted the tissue compositions from the International Commission on Radiation Protection Publication 23 (39), which were used in Acuros BV. The cross-section data in our package for particle transport simulation are from NIST XCOM database (40). We would like to mention that our package supports a certain number of material types but allows continuous density variations for each voxel. The mass attenuation coefficients for these materials were stored. During the simulation, the actual density of a voxel is multiplied with the mass attenuation coefficient to yield the actual attenuation coefficient. Because the particle transport physics in the relevant energy range for HDR brachytherapy, indicated by Fig. 1, is well known, we only present it here briefly. In this lowenergy range, three main photon interactions, Compton scattering, Rayleigh scattering, and photoelectric absorption, were considered. Compton scattering was modeled using a
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
389
atomic relaxation. The photons were tracked until their energies fell below the photon absorption energy threshold, that is, 10 keV in our package, after which the energies were deposited locally. The Woodcock tracking method (34) was used to increase the simulation efficiency of photon transport by avoiding voxel boundary crossing calculations in heterogeneous volumes. Because of the low-energy range of HDR brachytherapy, the secondary electrons generated by Compton scattering had a sub-millimeter electron range (18, 31), usually smaller than a typical voxel size for MC dose calculation. Hence, the electron transport was completed ignored with the energies being locally deposited. In addition, the electron equilibrium can further reduce the impact of ignoring electron transport on the resulting dose. We would like to mention that cautions should be made in the materials with much smaller densities, for example lung, when using gBMC, where electron transport may be needed. Fig. 1. The energy spectrum of the source photons that were recorded in the phase-space file of 192Ir.
Simulation process and GPU implementation free electron approximation with the Klein-Nishina crosssection formula (41). Binding effects and atomic relaxation were not modeled. For Rayleigh scattering, its differential cross-section per unit solid angle was approximated by the classical Thomson DCS, which was for scattering by a free electron at rest, along with the atomic form factor calculated by Hubbel et al.(42,43). The photoelectric absorption was modeled by assuming that all the energy was locally deposited without generating secondary electrons, ignoring the
Preparation The flowchart of the simulation process is shown in Fig. 2. Fig. 2a illustrates the overall code structure. The physics and phantom/patient data, such as the crosssection data per unit mass for each material at different energies, voxelized phantom/patient density, and material volume, are loaded and transferred to the GPU texture memory at the beginning. They will be available to all GPU threads during simulation with fast access because of cached
Fig. 2. Flowchart of the gBMC dose package. (a) Main structure; (b) detailed steps of the simulation within each energy bin (Step 5 in the main structure).
390
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
Fig. 3. Dose of a single source in the lung slab phantom. (a1ea3) show the slab geometry in transverse, coronal, and sagittal views, respectively. The source position is pointed out by the red arrow in (a1) and (a3). (b1eb3) show the dose calculated by TG-43 for comparison. (c1ec3) show the dose calculated by Acuros BV and (d1ed3) show the dose calculated by our gBMC.
memory access. Texture memory also supports hardwarebased linear interpolation on the data. Specifications for a treatment plan, including source dwell positions and source direction are loaded from the plan file and transferred to the GPU constant memory because of the relatively small data sizes. Particles contained in the phase-space file are also loaded and divided into multiple energy bins stored in the GPU global memory. This is to allow particle transport simulation to be launched bin by bin and ensure particles simultaneously transported on the GPU are close in energy for a similar amount of simulation effort. This strategy avoids the existence of a long-life GPU thread that reduces the overall computational efficiency in parallel processing. Ten dose counters are allocated on GPU global memory. During dose calculations, a random counter is assigned for a history particle to score its dose deposition events. These 10 dose counters are used for estimation of the resulting uncertainty and may also help relieve GPU writing conflicts caused by simultaneously writing to the same memory address from different GPU threads (25). The next step is to preprocess the loaded data and generate necessary information for the actual simulation process. Given a user-specified total number of the particle histories to be simulated for a treatment plan, N, the number of histories for each energy bin is first determined as N NEi 5 N P pENi , where NpEi denotes the number of the i
pEi
source particles in the energy bin Ei in the phase-space file. We would like to point out that here NEi does not necessarily need to be an integer multiple of NpEi . Next, a lookup table is created to facilitate the assignment of dwell position of a source particle during dose calculations. A treatment plan is characterized by a series of Ns source dwell positions and corresponding dwell time Tj (j 5 1,.,N P s).P A table of cumulative probability PðSi Þ 5 Tj = Tj for a particle in the ith dwell position j!i
j!Ns
is created. In principle, each source particle can be assigned to a dwell position by sampling a random number, r, uniformly distributed between [0,1] and determining the dwell position index, i, that satisfies PðSi Þ#r!PðSiþ1 Þ. However, this strategy is not preferred in GPU computing because of the required searching operations involving frequent memory access and conditional statements, both of which leads to GPU thread divergence issues (30). To solve this problem, a lookup table is precalculated listing the dwell position index corresponding to an accumulated probability array p for p 5 0 : Dp : 1. With this lookup table, for a random number r˛½0; 1, the corresponding dwell position index value from the ½r=Dpth element of the table is determined. This lookup table was created for the current treatment plan and loaded into GPU constant memory at the beginning of MC simulation.
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
Dose calculation The actual dose calculation is performed in steps 5 and 6, which sequentially loops through all energy bins. Detailed steps of the simulation within each energy bin are shown in Fig. 2b. In GPU-based parallel processing, each time when the GPU kernel of particle transport was launched, M source particles were simulated simultaneously until the preset number of histories, NEi , for the current energy bin was reached. To avoid biasing the simulation, at the beginning of the MC simulation, a random start location, slEi ˛½0; NpEi 1, was assigned for each energy bin. Starting from this random location in the current energy bin, we began to load M source particles and transport them each time and looped to the start of the energy bin when its end was reached. Each GPU thread was responsible for the transport of one source particle. Before transporting a particle, a dwell position for this particle was assigned using the previously created lookup table. The particle’s position and direction are adjusted according to the dwell position information. Azimuthal particle redistribution is used to use the rotational symmetry of the source for latent variance reduction of the phase-space file (44). Particle transport simulation is then conducted until the particle either exits the patient/phantom or gets absorbed. Dose deposition is recorded during the transport process. Statistical uncertainty After all simulations are performed, statistical analysis is conducted over the dose recorded in the 10 counters to
391
obtain the average dose and statistical uncertainty for each voxel. There are two major uncertainty estimation methods, the batch method (45) and the history-by-history method (46). Although the history-by-history method has been reported to result in a smaller uncertainty on the estimated uncertainty itself than using the batch method with a low number of batches (46), this method may present serious challenges to parallel programming frameworks (47). Hence, in our gBMC package, we used the batch method to estimate the statistical uncertainty, which can be formulated as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 PNb i 5 1 Di D sD 5 ð1Þ Nb ðNb 1Þ Here Nb was the number of the dose counters, Di was the dose value for a voxel in the ith dose counter, and D was the mean dose value averaged over all the counters. With the statistical uncertainty of each voxel, the statistical uncertainty of the calculation was obtained by taking an average over voxels within 5% isodose line. Validation Several tests were performed to validate the gBMC package. The gBMC was first benchmarked against the TG-43 dose calculation method in an ‘‘infinitely’’ large water medium. Here we used an extremely large water phantom with a size of 100 100 100 cm3. A single source was placed at the center of the phantom. As we have
Fig. 4. Dose of a single source in the bone slab phantom. (a1ea3) show the slab geometry in transverse, coronal, and sagittal views, respectively. The source position is pointed out in (a1) and (a3). (b1eb3) show the dose calculated by TG-43 for comparison. (c1ec3) show the dose calculated by Acuros BV and (d1ed3) show the dose calculated by our gBMC.
392
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
introduced in ‘‘Preparation’’ section, the voxelized phantom/patient density and material volume are usually loaded and transferred to the GPU texture memory at the very beginning for MC simulation. However, because of the limited memory on a GPU card, we could not store the whole 3D voxelized phantom density volume and the corresponding material volume for this huge water phantom on GPU memory. Hence, gBMC was modified here for this case such that only the phantom dimension, size, its homogenous density value, and material type were stored on GPU. In addition, dose deposition was only recorded in the central 40 40 40 cm3 region, to be compared with that calculated by the TG-43 method. Here we consider Acuros BV as our reference dose calculation method and benchmarked our gBMC package against it in heterogeneous medium because Acuros BV
has been validated in many studies and is being increasingly used in clinical practice (14,48). Two slab phantom geometries were used. One was a water phantom with a 2cm-thick slab of lung equivalent material (Gammex Tissue Equivalent Materials, Gammex INC, Middleton, WI) inserted (lung slab phantom). The second one was a water phantom with a 1-cm-thick bone slab inserted (bone slab phantom). In both cases, one source position with 1-h dwell time was used for the dose calculation. The geometries of these two slabs used for the dose calculation were extracted from the CT scans, shown in Figs. 3 and 4 (a1ea3) in transverse, coronal, and sagittal views, with the source position pointed out. The CT numberto-material and density conversion curve, used both in Acuros BV and in our gBMC package, were defined using the RMI Gammex Tissue characterization phantom (Gammex Inc., Middleton, WI).
Fig. 5. Dose in water for a cylinder of 3.5 cm diameter with 270 shielding. (a1ea3) Show the voxelized digital phantom used for dose calculation in transverse, coronal, and sagittal views, respectively. The geometry in the dashed rectangular in (a2) is zoomed in for better display. The red points denote the source dwell positions included in the treatment plan. (b1eb3) Show the dose calculated by Acuros BV. (c1ec3) show the dose calculated by our gBMC package. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
393
0.1 0.1 0.1 cm . Each of them consisted of a uniform water medium with one of the 15 shielded cylinder applicators considered. These applicators included five different diameters (2.0, 2.3, 2.6, 3.0, and 3.5 cm) and three different partial shielding geometries (90 , 180 , and 270 ), as illustrated in Fig. 5 (a1ea3). In each case, dose of a treatment plan with eight dwell positions and equal dwell time (i.e., 100 s for each dwell position) was calculated. The results were compared with the dose calculated by Acuros BV. For final benchmarking of the gBMC package, dose calculations were performed for a patient case treated with a tandem and ovoid (T&O) CT-compatible applicator to demonstrate feasibility for clinical application. The patient was treated at our institution for cervical cancer. Patient data were collected for this study under the approval of our internal review board. A phantom used for dose calculation was created from the patient CT image having a size of 256 256 226 voxels and a resolution of 2.34 2.34 2 mm3. The material and density of the phantom were converted from the CT numbers using the CT number-to-material and density conversion curve mentioned earlier, except for the applicator, whose material type and density were overwritten to the actual values. There were 16 source dwell positions included in the treatment plan with eight positions in the tandem channel and 3
Fig. 6. Absolute dose differences in water for a single source between gBMC and TG-43 relative to the local dose.
Another group of cases to benchmark the gBMC package against Acuros BV was dose calculations with shielded cylinder applicators. Here, the Varian-shielded applicator set GM11004380 (Varian Medical Systems) was modeled. The geometry and material properties of the applicators were from published data (49); 15 voxelized digital phantoms were created with a resolution of
Fig. 7. Single-source dosimetry in water for varisource VS2000. (a) Radial dose function (bed) Anisotropic functions at r 5 1, 5, 20 cm, respectively, where r denotes the radial distance from the source. The relative differences between the gBMC results and the TG-43 results are also shown.
394
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
Fig. 8. Dose profiles along L1 in Fig. 3 and L2 in Fig. 4, respectively. The regions where the TG-43ecalculated dose deviated from the Acuros BV dose and the gBMC dose methods are emphasized by the dashed circles.
four in each ovoid channel. The prescription dose was 5.8 Gy for this patient. The efficiency of our gBMC package was tested on an NVidia GeForce GTX TITAN Black GPU card (Nvidia Corporation, Santa Clara, CA). The Acuros BV calculation used for comparison was performed on a desktop having a quad-core Intel Xeon E5620 CPU processor with 2.4-GHz frequency. We would like to mention that this efficiency comparison between the gBMC and the Acuros BV may not be very fair, considering that the Xeon E5620 CPU processor is relatively outdated compared with the TITAN Black GPU card (i.e., these two cards were released in the year 2010 and 2014, respectively).
Results Validation results Benchmark against TG-43 Dose calculation results comparing the TG-43 approach and gBMC in a water medium are presented in Fig. 6. Here the absolute differences between the two calculated doses relative to the local dose ðDgBMC DTG43 =DTG43 100%Þ is displayed in color wash on a plane on which the central axis of the radioactive source is. The TG-43 dose was calculated according to the two-dimensional dose-rate equation from the 1995 TG-43 protocol (2), using the dosimetry parameters calculated by Taylor and Rogers (16,50,51). It can be observed that the dose differences within 20 cm from the source are ! 2% except at the bottom left corner. Here, the relatively large discrepancies are ascribed to different handling of the connecting cable at one end of the radiative source between the two dose calculation methods. The radial dose function g(r) and the anisotropy functions F(r,q) at several different radial distances from the source (r 5 1, 5, and 20 cm) were derived from the gBMC-calculated dose distribution and compared with the corresponding data that we used in TG-43 calculation (Fig. 7). The differences in the radial dose function g(r) were !0.6%. The differences in the anisotropy functions F(r,q) were within 1% at most polar
angles except at the polar angles O173 , which corresponded to the connecting cable location. Benchmark against Acuros BV Slab phantoms. Dose distributions from a single source in the lung slab phantom calculated by Acuros BV and gBMC, respectively, are shown in Fig. 3 to give the readers an intuitive impression on the dosimetric effect of the tissue heterogeneity. The corresponding dose calculated using the TG-43 method by BrachyVision treatment planning system (Varian Medical Systems) is also shown for comparison purposes. Asymmetric dose distributions were observed from both the Acuros BV and gBMC calculations because of the inserted slab of the lung equivalent material. In contrast, TG-43 neglected this tissue inhomogeneity and resulted in a symmetric dose distribution. Similar results for the bone slab phantom can be found in Fig. 4. To better compare these three methods, the dose profiles along two lines which are perpendicular to the source and pass through its center (i.e. line L1 shown in Fig. 3 and line L2 shown in Fig. 4) are plotted in Fig. 8 in the log scale of the dose. It can be observed that the Acuros BV and the gBMC-calculated doses match well with each other. Although the TG-43ecalculated dose distribution shows an underestimation in the water which is after the lung slab along the particle penetration direction, an overestimation in the water which is after the bone slab along the particle penetration direction, and an overestimation around the boundary of water and air. These regions are emphasized with dashed circles in the figures. Quantitative comparisons were performed by calculating the absolute differences between the Acuros BV dose and the gBMC dose relative to Table 1 Dose differences between Acuros BV and gBMC relative to the local dose for heterogeneous slab phantoms Phantom
98% max difference (%)
Average difference (%)
Lung slab phantom Bone slab phantom
1.88 2.78
0.56 0.87
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
395
Fig. 9. Dose profiles along the lines L1, L2, and L3 in Fig. 5.
local dose. The results are listed in Table 1. Dose points very close to the source (within 2.5 mm) and the dose in air were excluded from evaluation. The 98% maximum difference was 1.88% for the lung slab and 2.78% for the bone slab, whereas the average difference was 0.56% and 0.87% for the two cases, respectively. These numbers demonstrated a good agreement between gBMC and the Acuros BV, especially in terms of capturing the attenuation through the inserted lung/bone slab and the lack of backscatter at water-air interfaces. Shielded cylinder applicator cases. Because page limitations, only the dose distributions for a shielded applicator of 3.5 cm diameter with 270 shielding angle are shown in Fig. 5. Only the dose outside the applicator was calculated by Acuros BV and shown here for comparison. The dose profiles along the three lines in Fig. 5b1 are plotted in Fig. 9, indicating a good agreement between the doses calculated by Acuros BV and gBMC. Table 2 lists the absolute dose differences between these two dose distributions relative to the reference dose, which is defined to be the dose 5 mm away from the applicator surface (outward direction and at unshielded region) and on the plane which is in the middle of the fourth and fifth dwell position. The average dose differences were !0.54%, and 98-percentile dose differences was !2.04% for all the 15 cases. In Acuros BV, the dose calculation grid has a resolution of 0.1 0.1 0.1 cm3 and that of the output dose grid was 0.25 0.25 0.1 cm3. Our gBMC dose was resampled to the same dose grid resolution for comparison. Some voxelization effects were seen around the applicator surface in Fig. 5. This can be ascribed to an option in Acuros that cleared dose inside the applicator after dose calculation. The display of a low-resolution dose array at a high-gradient region (from a high-dose level outside the applicator surface immediately to 0 set by Acuros) made the voxelization effect obvious. Dose difference close to the applicator tip between Acuros and gBMC were also observed. It may come from different source modeling along the source direction, as implied by the
relatively large discrepancies in this region shown in Figs. 6 and 7. Patient case The color wash dose calculated by gBMC for the T&O patient case is shown in the first row of Fig. 10. Isodose lines corresponding to 30%, 50%, 100%, and 120% of the prescription dose are shown in red in the second row. For comparison purposes, the isodose lines of the dose calculated by Acuros BV are also shown in blue. These two sets of isodose lines almost overlap with each other, indicating good agreements between the two dose distributions. The corresponding dose-volume histograms for this patient are shown in Fig. 11, which also show a satisfactory agreement. Dosimetric indices for evaluation of this treatment plan were calculated and presented in Table 3 for comparison. The absolute difference of clinical target volume D90 was 0.78%. The differences of D2cc were 2.17%, 0.96%, and 1.22% for sigmoid, rectum, and bladder, respectively. For this case, ~1.3 billion source histories were used in our MC simulation to compare with the Acuros calculated result. It took us ~47.5 s to simulate these 1.3 billion histories for this treatment plan with the gBMC package, achieving a small statistical uncertainty, that is, on average ~0.15% of the prescription dose within the 5% isodose line. In contrast, the computational time was ~176 s using Acuros BV on the quad-core CPU platform. We would like Table 2 Absolute dose differences between Acuros BV and gBMC for applicator shielding Applicator diameter (cm)
Difference (%) Shielding degree ( ) 2.0
2.3
2.6
3.0
3.5
Average 98-percentile Average 98-percentile Average 98-percentile
0.48 1.83 0.45 1.77 0.49 1.95
0.49 1.80 0.46 1.79 0.49 1.80
0.45 1.72 0.49 1.87 0.51 1.91
0.43 1.72 0.51 1.94 0.54 1.96
90 180 270
0.46 1.92 0.45 1.81 0.49 2.04
396
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
Fig. 10. Dose results for a tandem and ovoid brachytherapy patient case. The first row corresponds to the color wash dose calculated by our gBMC package, shown in transverse, coronal, and sagittal views. The second row shows the isodose lines corresponding to 120%, 100%, 50%, and 30% of the prescription dose. The red lines are for the dose calculated by gBMC and the blue lines are for the dose calculated by Acuros BV. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
to mention again that this efficiency comparison is not very fair because the Xeon E5620 CPU processor on which Acuros BV is installed is relatively outdated compared with the TITAN Black GPU card used for gBMC. In addition, the computational efficiency with Acuros is expected to scale with the number of the cores on the CPU processor available for the calculation (52). For some applications which may have a relatively loose criterion on the statistical uncertainty of the MC-calculated results, such as a secondary dose calculation for QA purpose, the dose calculation for inverse treatment planning, the MC simulation time can be shortened. For example, for this patient case, it only took gBMC about 1 s to achieve ~1% average statistical uncertainty.
and efficiency, compared with the Acuros BV method, demonstrating the feasibility for clinical applications. We would like to mention that because the gBMC package was developed on top of our established GPU-based MC packages, gDPM for dose calculation for external radiotherapy (22, 24), and gCTD for imaging dose calculation for CT and cone-beam CT (35), both of which have been validated against the conventional CPU-based MC codes, we did not benchmark our gBMC package against the CPU-based MC codes in this article. Instead, here we
Discussion and conclusions This article presents our recently developed GPU-based fast Monte Carlo dose calculation package for HDR brachytherapy, gBMC. This package has been benchmarked against the TG-43 dose calculation method in an idea infinitely large water phantom and the Acuros BV dose calculation method in several heterogeneous geometries. Satisfactory agreement was observed for all test cases, validating the accuracy of gBMC. This package was also tested for a T&O patient case and achieved both high accuracy
Fig. 11. Dose-volume histograms for the tandem and ovoid brachytherapy patient case in Fig. 10. Solid lines are for dose calculated by our gBMC package and dash lines are for Acuros BV.
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
397 192
Table 3 Dosimetric endpoints for the Tandem and ovoid patient case
Dose engine
D90 (Gy)
D2cc (Gy)
D2cc (Gy)
D2cc (Gy)
gBMC platform only includes the Varian VS2000 Ir source. Future versions will include more source simulations to generate phase-space files for other common HDR brachytherapy sources.
Acuros BV gBMC Difference (%)
10.26 10.18 0.78
2.77 2.83 2.17
4.15 4.19 0.96
5.74 5.81 1.22
References
CTV
Sigmoid
Rectum
Bladder
CTV = clinical target volume.
consider the Acruos BV as our reference dose calculation method for the complex heterogeneous geometries. That is because the Acuros BV has been validated in many studies and is being increasingly used in clinical practice (48,52). There may be a potential inaccuracy because of the voxelization of the applicator. Generally speaking, the smaller voxel size is, the more accurately we can model the applicator geometry. However, the computational efficiency was also significantly reduced with finer a resolution, if we would like to maintain a constant statistical uncertainty. This applicator voxelization problem for MC simulation has been studied extensively by Fonseca, Landry (53). In addition, the GPU memory is usually a limiting factor for the storage of a very large phantom with a fine resolution. It is for these considerations that we selected a modest voxel size that was found to be suitable for computations while maintaining sufficient accuracy. Although we have only presented applications of gBMC for forward dose calculations, there are potential applications in speeding up inverse treatment planning. In such a problem, repeated dose calculations for each dwell position are needed in which the dwell time is adjusted to yield a clinically acceptable total dose distribution. This is more computationally intensive than a typical forward dose calculation, and the current clinically used Acuros BV system may be too slow for this situation. Our ongoing effort is to integrate gBMC into an inverse optimization system, which is particularly important in cases where consideration of heterogeneities is clinically preferred, for example, in accelerated partial breast irradiation (54e56). After the success achieved in this work, we are in the process of improving the package to accommodate other clinical scenarios. The first direction is to incorporate geometry-modeling modules to enable particle transport simulation in a continuous geometry parameterized by analytical functions. This is an important step to support dose calculations in the presence of radioactive seeds and applicators that are defined by a continuous geometry. Also, we would like to extend gBMC into dose calculations for lower energies to support low-doseerate brachytherapy using implanted radioactive seeds. Although the physics data and transport model in the current gBMC platform are expected to be valid, comprehensive benchmarking in clinical cases are needed. Third, the current
[1] Thomadsen BR, Williamson JF, Rivard MJ, Meigooni AS. Anniversary paper: past and current issues, and trends in brachytherapy physics. Med Phys 2008;35:4708e4723. [2] Nath R, Anderson LL, Luxton G, et al. Dosimetry of interstitial brachytherapy sources: Recommendations of the AAPM Radiation Therapy Committee Task group No. 43. Energy 1995;1:3. [3] Rivard MJ, Coursey BM, DeWerd LA, et al. Update of AAPM Task Group No. 43 Report: A revised AAPM protocol for brachytherapy dose calculations. Med Phys 2004;31:633e674. [4] Beaulieu L, Tedgren AC, Carrier J-F, et al. Report of the Task Group 186 on model-based dose calculation methods in brachytherapy beyond the TG-43 formalism: Current status and recommendations for clinical implementation. Med Phys 2012;39:6208e6236. [5] Serago CF, Houdek PV, Pisciotta V, et al. Scattering effects on the dosimetry of iridium-192. Med Phys 1991;18:1266e1270. [6] Mangold CA, Rijnders A, Georg D, et al. 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 2001;58:43e51. [7] Kassas B, Mourtada F, Horton JL, et al. Dose modification factors for 192Ir high-dose-rate irradiation using Monte Carlo simulation. J Appl Clin Med Phys 2006;7:28e34. [8] Raffi JA, Davis SD, Hammer CG, et al. Determination of exit skin dose for I192r intracavitary accelerated partial breast irradiation with thermoluminescent dosimeters. Med Phys 2010;37:2693e 2702. [9] Lymperopoulou G, Papagiannis P, Angelopoulos A, et al. A dosimetric comparison of 169Yb and 192Ir for HDR brachytherapy of the breast, accounting for the effect of finite patient dimensions and tissue inhomogeneities. Med Phys 2006;33:4583e4589. [10] Yan X, Poon E, Reniers B, et al. Comparison of dose calculation algorithms for colorectal cancer brachytherapy treatment with a shielded applicator. Med Phys 2008;35:4824e4830. [11] Ahnesj€o A. Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media. Med Phys 1989;16:577e 592. [12] Mackie T, Scrimger J, Battista J. A convolution method of calculating dose for 15-MV x rays. Med Phys 1985;12:188e196. [13] Gifford KA, Horton JL Jr, Wareing TA, et al. Comparison of a finite-element multigroup discrete-ordinates code with Monte Carlo for radiotherapy calculations. Phys Med Biol 2006;51:2253. [14] Vassiliev ON, Wareing TA, McGhee J, et al. Validation of a new grid-based Boltzmann equation solver for dose calculation in radiotherapy with photon beams. Phys Med Biol 2010;55:581. [15] Rogers D. Fifty years of Monte Carlo simulations for medical physics. Phys Med Biol 2006;51:R287. [16] Taylor R, Yegin G, Rogers D. Benchmarking brachydose: Voxel based EGSnrc Monte Carlo calculations of TG-43 dosimetry parameters. Med Phys 2007;34:445e457. [17] Perez-Calatayud J, Ballester F, Das RK, et al. Dose calculation for photon-emitting brachytherapy sources with average energy higher than 50 keV: Report of the AAPM and ESTRO. Med Phys 2012;39:2904e2929. [18] Chibani O, Ma CC. HDRMC, an accelerated Monte Carlo dose calculator for high dose rate brachytherapy with CT-compatible applicators. Med Phys 2014;41:051712. [19] Williamson JF. Monte Carlo evaluation of kerma at a point for photon transport problems. Med Phys 1987;14:567e576.
398
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
[20] Salvat F, Fernndez-Varena J, Sampau J. Penelope-2011, a code system for Monte Carlo simulation of electron and photon transport, OECD. Issy les Moulineaux, France: NEA Data Bank; 2011. [21] Badal A, Badano A. Accelerating Monte Carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit. Med Phys 2009;36:4878e4880. [22] Jia X, Gu X, Sempau J, et al. Development of a GPU-based Monte Carlo dose calculation code for coupled electronephoton transport. Phys Med Biol 2010;55:3077. [23] Hissoiny S, Ozell B, Bouchard H, et al. GPUMCD: a new GPUoriented Monte Carlo dose calculation platform. Med Phys 2011; 38:754e764. [24] Jia X, Gu X, Graves YJ, et al. GPU-based fast Monte Carlo simulation for radiotherapy dose calculation. Phys Med Biol 2011;56:7017. [25] Jia X, Sch€ umann J, Paganetti H, et al. GPU-based fast Monte Carlo dose calculation for proton therapy. Phys Med Biol 2012;57:7783. [26] Jia X, Yan H, Cervi~no L, et al. A GPU tool for efficient, accurate, and realistic simulation of cone beam CT projections. Med Phys 2012;39:7368e7378. [27] Jahnke L, Fleckenstein J, Wenz F, et al. GMC: a GPU implementation of a Monte Carlo dose calculation based on Geant4. Phys Med Biol 2012;57:1217. [28] Townson RW, Jia X, Tian Z, et al. GPU-based Monte Carlo radiotherapy dose calculation using phase-space sources. Phys Med Biol 2013;58:4341. [29] Tian Z, Shi F, Folkerts M, et al. A GPU OpenCL based crossplatform Monte Carlo dose calculation engine (goMC). Phys Med Biol 2015;60:7419. [30] Jia X, Ziegenhein P, Jiang SB. GPU-based high-performance computing for radiation therapy. Phys Med Biol 2014;59:R151. [31] Hissoiny S, D’Amours M, Ozell B, et al. Sub-second high dose rate brachytherapy Monte Carlo dose calculations with bGPUMCD. Med Phys 2012;39:4559e4567. [32] Hissoiny S, Ozell B, Despres P, et al. Validation of GPUMCD for lowenergy brachytherapy seed dosimetry. Med Phys 2011;38:4101e4107. Magnoux V, Hissoiny S, et al. Fast GPU-based Monte [33] Bonenfant E, Carlo simulations for LDR prostate brachytherapy. Phys Med Biol 2015;60:4973. [34] Woodcock E, Murphy T, Hemmings P, et al. Techniques used in the GEM code for Monte Carlo neutronics calculations in reactors and other systems of complex geometry. In Proceedings of the Conference on the Applications of Computing Methods to Reactor Problems. 1965. [35] Jia X, Yan H, Gu X, et al. Fast Monte Carlo simulation for patientspecific CT/CBCT imaging dose calculation. Phys Med Biol 2012; 57:577e590. [36] Agostinelli S, Allison J, Amako K, et al. In: W Barletta, GEANT4-a simulation toolkit. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors, and Associated Equipment, Elsevier, 2003. 506: p. 250e303. [37] Rasmussen BE, Davis SD, Schmidt CR, et al. Comparison of airkerma strength determinations for HDR 192Ir sources. Med Phys 2011;38:6721e6729. [38] Zhang M, Zou W, Chen T, et al. Parameterization of brachytherapy source phase space file for Monte Carlo-based clinical brachytherapy dose calculation. Phys Med Biol 2014;59:455.
[39] Richmond C. ICRP publication 23. Int J Radiat Biol Relat studies Phys Chem Med 1985;48:285. [40] Berger M, Hubbell J, Seltzer S, et al. XCOM: Photon cross sections database. NIST Stand Reference Database 8; Gaithersburg, MD: National Institute of Standards and Technology; 2013. € [41] Klein O, Nishina Y. Uber die Streuung von Strahlung durch freie Elektronen nach der neuen relativistischen Quantendynamik von Dirac. Z f€ur Physik 1929;52:853e868. [42] Hubbell J, Veigele WJ, Briggs E, et al. Atomic form factors, incoherent scattering functions, and photon scattering cross sections. J Phys Chem Ref Data 1975;4:471e538. [43] Bohr, N., Atomic physics and human knowledge. Hoboken, NJ: John Wiley & Sons, 1958. [44] Sempau J, Sanchez-Reyes A, Salvat F, et al. Monte Carlo simulation of electron beams from an accelerator head using PENELOPE. Phys Med Biol 2001;46:1163e1186. [45] Rogers D, Kawrakow I, Seuntjens J, et al., NRC user codes for EGSnrc. NRCC Report PIRS-702 (Rev. B), Missoula, MT: Northen Rockies Coordination Center, 2003. [46] Walters B, Kawrakow I, Rogers D. History by history statistical estimators in the BEAM code system. Med Phys 2002;29: 2745e2752. [47] Renaud M-A, Roberge D, Seuntjens J. Latent uncertainties of the precalculated track Monte Carlo method. Med Phys 2015;42: 479e490. [48] Bush K, Gagne I, Zavgorodni S, et al. Dosimetric validation of Acuros XB with Monte Carlo methods for photon dose calculations. Med Phys 2011;38:2208e2221. [49] Petrokokkinos L, Zourari K, Pantelis E, et al. Dosimetric accuracy of a deterministic radiation transport based I192r brachytherapy treatment planning system. Part II: Monte Carlo and experimental verification of a multiple source dwell position plan employing a shielded applicator. Med Phys 2011;38: 1981e1992. [50] Taylor R, Rogers D. EGSnrc Monte Carlo calculated dosimetry parameters for Ir192 and Yb169 brachytherapy sources. Med Phys 2008;35:4933e4944. [51] Taylor R, Rogers D. More accurate fitting of I125 and Pd103 radial dose functions. Med Phys 2008;35:4242e4250. [52] Failla GA, Wareing T, Archambault Y, et al. Acuros XB advanced dose calculation for the Eclipse Treatment Planning System. Palo Alto, CA: Varian Medical Systems; 2010. [53] Fonseca GP, Landry G, White S, et al. The use of tetrahedral mesh geometries in Monte Carlo simulation of applicator based brachytherapy dose distributions. Phys Med Biol 2014;59:5921. [54] D’Amours M, Carrier J, Lessard E, et al. Th-c-aud A-03: A new approach for afterloading brachytherapy inverse planned dose optimization based on the accurate Monte Carlo method. Med Phys 2008;35:2969. [55] Hedtj€arn H, Carlsson GA, Williamson JF. Accelerated Monte Carlo based dose calculations for brachytherapy planning using correlated sampling. Phys Med Biol 2002;47:351. [56] Rivard MJ, Melhus CS, Granero D, et al. An approach to using conventional brachytherapy software for clinical treatment planning of complex, Monte Carlo-based brachytherapy dose distributions. Med Phys 2009;36:1968e1975.
Physics
Monte Carlo dose calculations for high-doseerate brachytherapy using GPU-accelerated processing Z. Tian1,*, M. Zhang2, B. Hrycushko1, K. Albuquerque1, S.B. Jiang1, X. Jia1,* 1
Department of Radiation Oncology, University of Texas Southwestern Medical Center, Dallas, TX Department of Radiation Oncology, Rutgers Cancer Institute of New Jersey, New Brunswick, NJ
2
ABSTRACT
PURPOSE: Current clinical brachytherapy dose calculations are typically based on the Association of American Physicists in Medicine Task Group report 43 (TG-43) guidelines, which approximate patient geometry as an infinitely large water phantom. This ignores patient and applicator geometries and heterogeneities, causing dosimetric errors. Although Monte Carlo (MC) dose calculation is commonly recognized as the most accurate method, its associated long computational time is a major bottleneck for routine clinical applications. This article presents our recent developments of a fast MC dose calculation package for high-doseerate (HDR) brachytherapy, gBMC, built on a graphics processing unit (GPU) platform. METHODS AND MATERIALS: gBMC-simulated photon transport in voxelized geometry with physics in 192Ir HDR brachytherapy energy range considered. A phase-space file was used as a source model. GPU-based parallel computation was used to simultaneously transport multiple photons, one on a GPU thread. We validated gBMC by comparing the dose calculation results in water with that computed TG-43. We also studied heterogeneous phantom cases and a patient case and compared gBMC results with Acuros BV results. RESULTS: Radial dose function in water calculated by gBMC showed !0.6% relative difference from that of the TG-43 data. Difference in anisotropy function was !1%. In two heterogeneous slab phantoms and one shielded cylinder applicator case, average dose discrepancy between gBMC and Acuros BV was !0.87%. For a tandem and ovoid patient case, good agreement between gBMC and Acruos BV results was observed in both isodose lines and dose-volume histograms. In terms of the efficiency, it took ~47.5 seconds for gBMC to reach 0.15% statistical uncertainty within the 5% isodose line for the patient case. CONCLUSIONS: The accuracy and efficiency of a new GPU-based MC dose calculation package, gBMC, for HDR brachytherapy make it attractive for clinical applications. 2016 American Brachytherapy Society. Published by Elsevier Inc. All rights reserved.
Keywords:
High-doseerate brachytherapy; Monte Carlo; Dose calculation; GPU
Introduction High-doseerate (HDR) brachytherapy has proved to be a highly successful treatment for cancers in a variety of disease sites (1), such as cervical cancer, breast cancer, prostate cancer, etc. By placing an applicator in the targeted area and
Received 9 November 2015; received in revised form 26 January 2016; accepted 27 January 2016. * Corresponding authors. Department of Radiation Oncology, University of Texas Southwestern Medical Center, Dallas, TX 75390. Tel.: 214-648-5031; fax: 214-645-2885. E-mail addresses: [email protected] (Z. Tian), [email protected] (X. Jia).
having a motor-controlled radioactive source dwell in a series of planned positions with preset time, HDR brachytherapy delivers a highly conformed dose to the target while sparing nearby critical organs. Accurate dose calculation in patients is a critical step in treatment planning for HDR brachytherapy to ensure good target coverage and sharp dose gradients for critical structure sparing. The standard dosimetry protocol currently used in clinical brachytherapy is based on the Association of American Physicists in Medicine Task Group report 43 (TG-43) (2, 3). This method approximates a three-dimensional (3D) dose distribution by superposition of interpolated and precalculated discrete dose matrices in an infinitely large water medium. Although this approach is straightforward and fast, the influence of material
1538-4721/$ - see front matter 2016 American Brachytherapy Society. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.brachy.2016.01.006
388
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
heterogeneity is ignored (4), causing dosimetric errors. For instance, the skin dose has been shown to be overestimated by as large as 15% (5e8). In breast brachytherapy, the TG43 approach tends to underestimate dose to lung by ~10% (9). In addition, neglecting dosimetric impacts of nontissue heterogeneities, such as applicator material and the shielding material, is another concern. It has been shown that ignoring the shielding of an HDR colorectal applicator leads to an overestimation of the target volume within 100% isodose lines by O5% (10), which implies potential tumor cold spots. Therefore, it is proposed in Task Group report 186 (4) to move toward more accurate model-based dose calculation algorithms for brachytherapy, such as collapsed-cone superposition/convolution methods (11, 12), deterministic solutions to the linear Boltzmann transport equation (13, 14), and Monte Carlo (MC) simulations (15). MC has been considered as one of the most accurate dose calculation method for determining radiation dose deposition in heterogeneous materials. Particle transport and interactions with the medium are simulated based on fundamental physics principles (4,15e17). However, MC simulations have been associated with extremely long calculation time, impeding applications in clinical practice. An accelerated CPU-based MC dose calculator (HDRMC) was recently developed (18) which used a fast rayetracing technique to transport photons through a 3D mesh of voxels representing the patient anatomy. It took ~5 min to calculate a typical HDR brachytherapy treatment plan consisting of 20e30 dwell positions. A track length estimator was used to improve the efficiency of the MC simulation (19). This track length estimator requires calculation of mass energy absorption coefficient before the simulation run, and the absorption energies used in the simulation must be identical to that used to determine the mass energy absorption (20). Recently, there have been tremendous research efforts in accelerating MC simulations using graphics processing unit (GPU) platforms (21e30). For HDR brachytherapy dose calculations, bGPUMCD (31) was developed by extending the photon transport part of the GPUMCD system (23) into the brachytherapy energy range. By using GPU’s rapid parallel processing power, GPU-friendly parallelization schemes, and a track length estimator, the dose calculation can be completed in seconds. This MC package was also extended to low-doseerate brachytherapy calculations (32,33). Grid-based Boltzmann equation solvers are another category of model-based dose calculation methods for brachytherapy. By directly solving the linear Boltzmann transport equation, the dose distribution is derived from fundamental physics principles and, hence, achieves guaranteed accuracy (14). This algorithm has been introduced for brachytherapy dose calculations in a clinic setting as Acuros BV (Varian Medical Systems, Palo Alto, CA). However, the calculation time for a treatment plan typically ranges in several minutes on a multicore CPU platform. This speed may be acceptable for plan dose calculation but becomes
burdensome when repeated calculations are needed, such as in inverse plan optimization. In this report, we will present our recent developments of a GPU-based MC dose engine for HDR brachytherapy, named gBMC (GPU-based brachytherapy Monte Carlo). We will show high computational efficiency without using variance reduction techniques except the Woodcock tracking method, a method that has been widely used in GPU-based photon transport simulation tools (34). To validate our gBMC package, we have benchmarked against the TG-43 results in a water phantom and against Acuros BV results in a heterogeneous slab phantom and shielded cylinder applicators. The feasibility of applying our gBMC package for clinical practice will also be shown in a typical HDR brachytherapy patient case.
Methods and materials Source modeling and physics The gBMC package was developed on top of our established GPU-based MC packages, gDPM for dose calculation for external radiotherapy (22, 24) and gCTD for imaging dose calculation of CT and cone-beam CT (35), both of which have been validated against the conventional CPU-based MC codes. The gBMC package currently focuses on 192Ir, the most commonly used radionuclide for HDR brachytherapy. The source model in gBMC was a phase-space file generated from the MC simulation of radioactive source decay and particle transport for a VariSource VS2000 192Ir HDR source (Varian Medical Systems) using GEANT4 (36). The source was modeled as a radioactive core inside a connecting cable with geometry and material composition following published data (37,38). For the 2 107 decay events sampled uniformly throughout the radioactive core, about 4.2 107 photons exited the surface of the cable, which were stored in the phase-space file. The energy spectrum of these recorded photons is shown in Fig. 1. In our gBMC package, we adopted the tissue compositions from the International Commission on Radiation Protection Publication 23 (39), which were used in Acuros BV. The cross-section data in our package for particle transport simulation are from NIST XCOM database (40). We would like to mention that our package supports a certain number of material types but allows continuous density variations for each voxel. The mass attenuation coefficients for these materials were stored. During the simulation, the actual density of a voxel is multiplied with the mass attenuation coefficient to yield the actual attenuation coefficient. Because the particle transport physics in the relevant energy range for HDR brachytherapy, indicated by Fig. 1, is well known, we only present it here briefly. In this lowenergy range, three main photon interactions, Compton scattering, Rayleigh scattering, and photoelectric absorption, were considered. Compton scattering was modeled using a
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
389
atomic relaxation. The photons were tracked until their energies fell below the photon absorption energy threshold, that is, 10 keV in our package, after which the energies were deposited locally. The Woodcock tracking method (34) was used to increase the simulation efficiency of photon transport by avoiding voxel boundary crossing calculations in heterogeneous volumes. Because of the low-energy range of HDR brachytherapy, the secondary electrons generated by Compton scattering had a sub-millimeter electron range (18, 31), usually smaller than a typical voxel size for MC dose calculation. Hence, the electron transport was completed ignored with the energies being locally deposited. In addition, the electron equilibrium can further reduce the impact of ignoring electron transport on the resulting dose. We would like to mention that cautions should be made in the materials with much smaller densities, for example lung, when using gBMC, where electron transport may be needed. Fig. 1. The energy spectrum of the source photons that were recorded in the phase-space file of 192Ir.
Simulation process and GPU implementation free electron approximation with the Klein-Nishina crosssection formula (41). Binding effects and atomic relaxation were not modeled. For Rayleigh scattering, its differential cross-section per unit solid angle was approximated by the classical Thomson DCS, which was for scattering by a free electron at rest, along with the atomic form factor calculated by Hubbel et al.(42,43). The photoelectric absorption was modeled by assuming that all the energy was locally deposited without generating secondary electrons, ignoring the
Preparation The flowchart of the simulation process is shown in Fig. 2. Fig. 2a illustrates the overall code structure. The physics and phantom/patient data, such as the crosssection data per unit mass for each material at different energies, voxelized phantom/patient density, and material volume, are loaded and transferred to the GPU texture memory at the beginning. They will be available to all GPU threads during simulation with fast access because of cached
Fig. 2. Flowchart of the gBMC dose package. (a) Main structure; (b) detailed steps of the simulation within each energy bin (Step 5 in the main structure).
390
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
Fig. 3. Dose of a single source in the lung slab phantom. (a1ea3) show the slab geometry in transverse, coronal, and sagittal views, respectively. The source position is pointed out by the red arrow in (a1) and (a3). (b1eb3) show the dose calculated by TG-43 for comparison. (c1ec3) show the dose calculated by Acuros BV and (d1ed3) show the dose calculated by our gBMC.
memory access. Texture memory also supports hardwarebased linear interpolation on the data. Specifications for a treatment plan, including source dwell positions and source direction are loaded from the plan file and transferred to the GPU constant memory because of the relatively small data sizes. Particles contained in the phase-space file are also loaded and divided into multiple energy bins stored in the GPU global memory. This is to allow particle transport simulation to be launched bin by bin and ensure particles simultaneously transported on the GPU are close in energy for a similar amount of simulation effort. This strategy avoids the existence of a long-life GPU thread that reduces the overall computational efficiency in parallel processing. Ten dose counters are allocated on GPU global memory. During dose calculations, a random counter is assigned for a history particle to score its dose deposition events. These 10 dose counters are used for estimation of the resulting uncertainty and may also help relieve GPU writing conflicts caused by simultaneously writing to the same memory address from different GPU threads (25). The next step is to preprocess the loaded data and generate necessary information for the actual simulation process. Given a user-specified total number of the particle histories to be simulated for a treatment plan, N, the number of histories for each energy bin is first determined as N NEi 5 N P pENi , where NpEi denotes the number of the i
pEi
source particles in the energy bin Ei in the phase-space file. We would like to point out that here NEi does not necessarily need to be an integer multiple of NpEi . Next, a lookup table is created to facilitate the assignment of dwell position of a source particle during dose calculations. A treatment plan is characterized by a series of Ns source dwell positions and corresponding dwell time Tj (j 5 1,.,N P s).P A table of cumulative probability PðSi Þ 5 Tj = Tj for a particle in the ith dwell position j!i
j!Ns
is created. In principle, each source particle can be assigned to a dwell position by sampling a random number, r, uniformly distributed between [0,1] and determining the dwell position index, i, that satisfies PðSi Þ#r!PðSiþ1 Þ. However, this strategy is not preferred in GPU computing because of the required searching operations involving frequent memory access and conditional statements, both of which leads to GPU thread divergence issues (30). To solve this problem, a lookup table is precalculated listing the dwell position index corresponding to an accumulated probability array p for p 5 0 : Dp : 1. With this lookup table, for a random number r˛½0; 1, the corresponding dwell position index value from the ½r=Dpth element of the table is determined. This lookup table was created for the current treatment plan and loaded into GPU constant memory at the beginning of MC simulation.
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
Dose calculation The actual dose calculation is performed in steps 5 and 6, which sequentially loops through all energy bins. Detailed steps of the simulation within each energy bin are shown in Fig. 2b. In GPU-based parallel processing, each time when the GPU kernel of particle transport was launched, M source particles were simulated simultaneously until the preset number of histories, NEi , for the current energy bin was reached. To avoid biasing the simulation, at the beginning of the MC simulation, a random start location, slEi ˛½0; NpEi 1, was assigned for each energy bin. Starting from this random location in the current energy bin, we began to load M source particles and transport them each time and looped to the start of the energy bin when its end was reached. Each GPU thread was responsible for the transport of one source particle. Before transporting a particle, a dwell position for this particle was assigned using the previously created lookup table. The particle’s position and direction are adjusted according to the dwell position information. Azimuthal particle redistribution is used to use the rotational symmetry of the source for latent variance reduction of the phase-space file (44). Particle transport simulation is then conducted until the particle either exits the patient/phantom or gets absorbed. Dose deposition is recorded during the transport process. Statistical uncertainty After all simulations are performed, statistical analysis is conducted over the dose recorded in the 10 counters to
391
obtain the average dose and statistical uncertainty for each voxel. There are two major uncertainty estimation methods, the batch method (45) and the history-by-history method (46). Although the history-by-history method has been reported to result in a smaller uncertainty on the estimated uncertainty itself than using the batch method with a low number of batches (46), this method may present serious challenges to parallel programming frameworks (47). Hence, in our gBMC package, we used the batch method to estimate the statistical uncertainty, which can be formulated as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 PNb i 5 1 Di D sD 5 ð1Þ Nb ðNb 1Þ Here Nb was the number of the dose counters, Di was the dose value for a voxel in the ith dose counter, and D was the mean dose value averaged over all the counters. With the statistical uncertainty of each voxel, the statistical uncertainty of the calculation was obtained by taking an average over voxels within 5% isodose line. Validation Several tests were performed to validate the gBMC package. The gBMC was first benchmarked against the TG-43 dose calculation method in an ‘‘infinitely’’ large water medium. Here we used an extremely large water phantom with a size of 100 100 100 cm3. A single source was placed at the center of the phantom. As we have
Fig. 4. Dose of a single source in the bone slab phantom. (a1ea3) show the slab geometry in transverse, coronal, and sagittal views, respectively. The source position is pointed out in (a1) and (a3). (b1eb3) show the dose calculated by TG-43 for comparison. (c1ec3) show the dose calculated by Acuros BV and (d1ed3) show the dose calculated by our gBMC.
392
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
introduced in ‘‘Preparation’’ section, the voxelized phantom/patient density and material volume are usually loaded and transferred to the GPU texture memory at the very beginning for MC simulation. However, because of the limited memory on a GPU card, we could not store the whole 3D voxelized phantom density volume and the corresponding material volume for this huge water phantom on GPU memory. Hence, gBMC was modified here for this case such that only the phantom dimension, size, its homogenous density value, and material type were stored on GPU. In addition, dose deposition was only recorded in the central 40 40 40 cm3 region, to be compared with that calculated by the TG-43 method. Here we consider Acuros BV as our reference dose calculation method and benchmarked our gBMC package against it in heterogeneous medium because Acuros BV
has been validated in many studies and is being increasingly used in clinical practice (14,48). Two slab phantom geometries were used. One was a water phantom with a 2cm-thick slab of lung equivalent material (Gammex Tissue Equivalent Materials, Gammex INC, Middleton, WI) inserted (lung slab phantom). The second one was a water phantom with a 1-cm-thick bone slab inserted (bone slab phantom). In both cases, one source position with 1-h dwell time was used for the dose calculation. The geometries of these two slabs used for the dose calculation were extracted from the CT scans, shown in Figs. 3 and 4 (a1ea3) in transverse, coronal, and sagittal views, with the source position pointed out. The CT numberto-material and density conversion curve, used both in Acuros BV and in our gBMC package, were defined using the RMI Gammex Tissue characterization phantom (Gammex Inc., Middleton, WI).
Fig. 5. Dose in water for a cylinder of 3.5 cm diameter with 270 shielding. (a1ea3) Show the voxelized digital phantom used for dose calculation in transverse, coronal, and sagittal views, respectively. The geometry in the dashed rectangular in (a2) is zoomed in for better display. The red points denote the source dwell positions included in the treatment plan. (b1eb3) Show the dose calculated by Acuros BV. (c1ec3) show the dose calculated by our gBMC package. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
393
0.1 0.1 0.1 cm . Each of them consisted of a uniform water medium with one of the 15 shielded cylinder applicators considered. These applicators included five different diameters (2.0, 2.3, 2.6, 3.0, and 3.5 cm) and three different partial shielding geometries (90 , 180 , and 270 ), as illustrated in Fig. 5 (a1ea3). In each case, dose of a treatment plan with eight dwell positions and equal dwell time (i.e., 100 s for each dwell position) was calculated. The results were compared with the dose calculated by Acuros BV. For final benchmarking of the gBMC package, dose calculations were performed for a patient case treated with a tandem and ovoid (T&O) CT-compatible applicator to demonstrate feasibility for clinical application. The patient was treated at our institution for cervical cancer. Patient data were collected for this study under the approval of our internal review board. A phantom used for dose calculation was created from the patient CT image having a size of 256 256 226 voxels and a resolution of 2.34 2.34 2 mm3. The material and density of the phantom were converted from the CT numbers using the CT number-to-material and density conversion curve mentioned earlier, except for the applicator, whose material type and density were overwritten to the actual values. There were 16 source dwell positions included in the treatment plan with eight positions in the tandem channel and 3
Fig. 6. Absolute dose differences in water for a single source between gBMC and TG-43 relative to the local dose.
Another group of cases to benchmark the gBMC package against Acuros BV was dose calculations with shielded cylinder applicators. Here, the Varian-shielded applicator set GM11004380 (Varian Medical Systems) was modeled. The geometry and material properties of the applicators were from published data (49); 15 voxelized digital phantoms were created with a resolution of
Fig. 7. Single-source dosimetry in water for varisource VS2000. (a) Radial dose function (bed) Anisotropic functions at r 5 1, 5, 20 cm, respectively, where r denotes the radial distance from the source. The relative differences between the gBMC results and the TG-43 results are also shown.
394
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
Fig. 8. Dose profiles along L1 in Fig. 3 and L2 in Fig. 4, respectively. The regions where the TG-43ecalculated dose deviated from the Acuros BV dose and the gBMC dose methods are emphasized by the dashed circles.
four in each ovoid channel. The prescription dose was 5.8 Gy for this patient. The efficiency of our gBMC package was tested on an NVidia GeForce GTX TITAN Black GPU card (Nvidia Corporation, Santa Clara, CA). The Acuros BV calculation used for comparison was performed on a desktop having a quad-core Intel Xeon E5620 CPU processor with 2.4-GHz frequency. We would like to mention that this efficiency comparison between the gBMC and the Acuros BV may not be very fair, considering that the Xeon E5620 CPU processor is relatively outdated compared with the TITAN Black GPU card (i.e., these two cards were released in the year 2010 and 2014, respectively).
Results Validation results Benchmark against TG-43 Dose calculation results comparing the TG-43 approach and gBMC in a water medium are presented in Fig. 6. Here the absolute differences between the two calculated doses relative to the local dose ðDgBMC DTG43 =DTG43 100%Þ is displayed in color wash on a plane on which the central axis of the radioactive source is. The TG-43 dose was calculated according to the two-dimensional dose-rate equation from the 1995 TG-43 protocol (2), using the dosimetry parameters calculated by Taylor and Rogers (16,50,51). It can be observed that the dose differences within 20 cm from the source are ! 2% except at the bottom left corner. Here, the relatively large discrepancies are ascribed to different handling of the connecting cable at one end of the radiative source between the two dose calculation methods. The radial dose function g(r) and the anisotropy functions F(r,q) at several different radial distances from the source (r 5 1, 5, and 20 cm) were derived from the gBMC-calculated dose distribution and compared with the corresponding data that we used in TG-43 calculation (Fig. 7). The differences in the radial dose function g(r) were !0.6%. The differences in the anisotropy functions F(r,q) were within 1% at most polar
angles except at the polar angles O173 , which corresponded to the connecting cable location. Benchmark against Acuros BV Slab phantoms. Dose distributions from a single source in the lung slab phantom calculated by Acuros BV and gBMC, respectively, are shown in Fig. 3 to give the readers an intuitive impression on the dosimetric effect of the tissue heterogeneity. The corresponding dose calculated using the TG-43 method by BrachyVision treatment planning system (Varian Medical Systems) is also shown for comparison purposes. Asymmetric dose distributions were observed from both the Acuros BV and gBMC calculations because of the inserted slab of the lung equivalent material. In contrast, TG-43 neglected this tissue inhomogeneity and resulted in a symmetric dose distribution. Similar results for the bone slab phantom can be found in Fig. 4. To better compare these three methods, the dose profiles along two lines which are perpendicular to the source and pass through its center (i.e. line L1 shown in Fig. 3 and line L2 shown in Fig. 4) are plotted in Fig. 8 in the log scale of the dose. It can be observed that the Acuros BV and the gBMC-calculated doses match well with each other. Although the TG-43ecalculated dose distribution shows an underestimation in the water which is after the lung slab along the particle penetration direction, an overestimation in the water which is after the bone slab along the particle penetration direction, and an overestimation around the boundary of water and air. These regions are emphasized with dashed circles in the figures. Quantitative comparisons were performed by calculating the absolute differences between the Acuros BV dose and the gBMC dose relative to Table 1 Dose differences between Acuros BV and gBMC relative to the local dose for heterogeneous slab phantoms Phantom
98% max difference (%)
Average difference (%)
Lung slab phantom Bone slab phantom
1.88 2.78
0.56 0.87
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
395
Fig. 9. Dose profiles along the lines L1, L2, and L3 in Fig. 5.
local dose. The results are listed in Table 1. Dose points very close to the source (within 2.5 mm) and the dose in air were excluded from evaluation. The 98% maximum difference was 1.88% for the lung slab and 2.78% for the bone slab, whereas the average difference was 0.56% and 0.87% for the two cases, respectively. These numbers demonstrated a good agreement between gBMC and the Acuros BV, especially in terms of capturing the attenuation through the inserted lung/bone slab and the lack of backscatter at water-air interfaces. Shielded cylinder applicator cases. Because page limitations, only the dose distributions for a shielded applicator of 3.5 cm diameter with 270 shielding angle are shown in Fig. 5. Only the dose outside the applicator was calculated by Acuros BV and shown here for comparison. The dose profiles along the three lines in Fig. 5b1 are plotted in Fig. 9, indicating a good agreement between the doses calculated by Acuros BV and gBMC. Table 2 lists the absolute dose differences between these two dose distributions relative to the reference dose, which is defined to be the dose 5 mm away from the applicator surface (outward direction and at unshielded region) and on the plane which is in the middle of the fourth and fifth dwell position. The average dose differences were !0.54%, and 98-percentile dose differences was !2.04% for all the 15 cases. In Acuros BV, the dose calculation grid has a resolution of 0.1 0.1 0.1 cm3 and that of the output dose grid was 0.25 0.25 0.1 cm3. Our gBMC dose was resampled to the same dose grid resolution for comparison. Some voxelization effects were seen around the applicator surface in Fig. 5. This can be ascribed to an option in Acuros that cleared dose inside the applicator after dose calculation. The display of a low-resolution dose array at a high-gradient region (from a high-dose level outside the applicator surface immediately to 0 set by Acuros) made the voxelization effect obvious. Dose difference close to the applicator tip between Acuros and gBMC were also observed. It may come from different source modeling along the source direction, as implied by the
relatively large discrepancies in this region shown in Figs. 6 and 7. Patient case The color wash dose calculated by gBMC for the T&O patient case is shown in the first row of Fig. 10. Isodose lines corresponding to 30%, 50%, 100%, and 120% of the prescription dose are shown in red in the second row. For comparison purposes, the isodose lines of the dose calculated by Acuros BV are also shown in blue. These two sets of isodose lines almost overlap with each other, indicating good agreements between the two dose distributions. The corresponding dose-volume histograms for this patient are shown in Fig. 11, which also show a satisfactory agreement. Dosimetric indices for evaluation of this treatment plan were calculated and presented in Table 3 for comparison. The absolute difference of clinical target volume D90 was 0.78%. The differences of D2cc were 2.17%, 0.96%, and 1.22% for sigmoid, rectum, and bladder, respectively. For this case, ~1.3 billion source histories were used in our MC simulation to compare with the Acuros calculated result. It took us ~47.5 s to simulate these 1.3 billion histories for this treatment plan with the gBMC package, achieving a small statistical uncertainty, that is, on average ~0.15% of the prescription dose within the 5% isodose line. In contrast, the computational time was ~176 s using Acuros BV on the quad-core CPU platform. We would like Table 2 Absolute dose differences between Acuros BV and gBMC for applicator shielding Applicator diameter (cm)
Difference (%) Shielding degree ( ) 2.0
2.3
2.6
3.0
3.5
Average 98-percentile Average 98-percentile Average 98-percentile
0.48 1.83 0.45 1.77 0.49 1.95
0.49 1.80 0.46 1.79 0.49 1.80
0.45 1.72 0.49 1.87 0.51 1.91
0.43 1.72 0.51 1.94 0.54 1.96
90 180 270
0.46 1.92 0.45 1.81 0.49 2.04
396
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
Fig. 10. Dose results for a tandem and ovoid brachytherapy patient case. The first row corresponds to the color wash dose calculated by our gBMC package, shown in transverse, coronal, and sagittal views. The second row shows the isodose lines corresponding to 120%, 100%, 50%, and 30% of the prescription dose. The red lines are for the dose calculated by gBMC and the blue lines are for the dose calculated by Acuros BV. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
to mention again that this efficiency comparison is not very fair because the Xeon E5620 CPU processor on which Acuros BV is installed is relatively outdated compared with the TITAN Black GPU card used for gBMC. In addition, the computational efficiency with Acuros is expected to scale with the number of the cores on the CPU processor available for the calculation (52). For some applications which may have a relatively loose criterion on the statistical uncertainty of the MC-calculated results, such as a secondary dose calculation for QA purpose, the dose calculation for inverse treatment planning, the MC simulation time can be shortened. For example, for this patient case, it only took gBMC about 1 s to achieve ~1% average statistical uncertainty.
and efficiency, compared with the Acuros BV method, demonstrating the feasibility for clinical applications. We would like to mention that because the gBMC package was developed on top of our established GPU-based MC packages, gDPM for dose calculation for external radiotherapy (22, 24), and gCTD for imaging dose calculation for CT and cone-beam CT (35), both of which have been validated against the conventional CPU-based MC codes, we did not benchmark our gBMC package against the CPU-based MC codes in this article. Instead, here we
Discussion and conclusions This article presents our recently developed GPU-based fast Monte Carlo dose calculation package for HDR brachytherapy, gBMC. This package has been benchmarked against the TG-43 dose calculation method in an idea infinitely large water phantom and the Acuros BV dose calculation method in several heterogeneous geometries. Satisfactory agreement was observed for all test cases, validating the accuracy of gBMC. This package was also tested for a T&O patient case and achieved both high accuracy
Fig. 11. Dose-volume histograms for the tandem and ovoid brachytherapy patient case in Fig. 10. Solid lines are for dose calculated by our gBMC package and dash lines are for Acuros BV.
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
397 192
Table 3 Dosimetric endpoints for the Tandem and ovoid patient case
Dose engine
D90 (Gy)
D2cc (Gy)
D2cc (Gy)
D2cc (Gy)
gBMC platform only includes the Varian VS2000 Ir source. Future versions will include more source simulations to generate phase-space files for other common HDR brachytherapy sources.
Acuros BV gBMC Difference (%)
10.26 10.18 0.78
2.77 2.83 2.17
4.15 4.19 0.96
5.74 5.81 1.22
References
CTV
Sigmoid
Rectum
Bladder
CTV = clinical target volume.
consider the Acruos BV as our reference dose calculation method for the complex heterogeneous geometries. That is because the Acuros BV has been validated in many studies and is being increasingly used in clinical practice (48,52). There may be a potential inaccuracy because of the voxelization of the applicator. Generally speaking, the smaller voxel size is, the more accurately we can model the applicator geometry. However, the computational efficiency was also significantly reduced with finer a resolution, if we would like to maintain a constant statistical uncertainty. This applicator voxelization problem for MC simulation has been studied extensively by Fonseca, Landry (53). In addition, the GPU memory is usually a limiting factor for the storage of a very large phantom with a fine resolution. It is for these considerations that we selected a modest voxel size that was found to be suitable for computations while maintaining sufficient accuracy. Although we have only presented applications of gBMC for forward dose calculations, there are potential applications in speeding up inverse treatment planning. In such a problem, repeated dose calculations for each dwell position are needed in which the dwell time is adjusted to yield a clinically acceptable total dose distribution. This is more computationally intensive than a typical forward dose calculation, and the current clinically used Acuros BV system may be too slow for this situation. Our ongoing effort is to integrate gBMC into an inverse optimization system, which is particularly important in cases where consideration of heterogeneities is clinically preferred, for example, in accelerated partial breast irradiation (54e56). After the success achieved in this work, we are in the process of improving the package to accommodate other clinical scenarios. The first direction is to incorporate geometry-modeling modules to enable particle transport simulation in a continuous geometry parameterized by analytical functions. This is an important step to support dose calculations in the presence of radioactive seeds and applicators that are defined by a continuous geometry. Also, we would like to extend gBMC into dose calculations for lower energies to support low-doseerate brachytherapy using implanted radioactive seeds. Although the physics data and transport model in the current gBMC platform are expected to be valid, comprehensive benchmarking in clinical cases are needed. Third, the current
[1] Thomadsen BR, Williamson JF, Rivard MJ, Meigooni AS. Anniversary paper: past and current issues, and trends in brachytherapy physics. Med Phys 2008;35:4708e4723. [2] Nath R, Anderson LL, Luxton G, et al. Dosimetry of interstitial brachytherapy sources: Recommendations of the AAPM Radiation Therapy Committee Task group No. 43. Energy 1995;1:3. [3] Rivard MJ, Coursey BM, DeWerd LA, et al. Update of AAPM Task Group No. 43 Report: A revised AAPM protocol for brachytherapy dose calculations. Med Phys 2004;31:633e674. [4] Beaulieu L, Tedgren AC, Carrier J-F, et al. Report of the Task Group 186 on model-based dose calculation methods in brachytherapy beyond the TG-43 formalism: Current status and recommendations for clinical implementation. Med Phys 2012;39:6208e6236. [5] Serago CF, Houdek PV, Pisciotta V, et al. Scattering effects on the dosimetry of iridium-192. Med Phys 1991;18:1266e1270. [6] Mangold CA, Rijnders A, Georg D, et al. 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 2001;58:43e51. [7] Kassas B, Mourtada F, Horton JL, et al. Dose modification factors for 192Ir high-dose-rate irradiation using Monte Carlo simulation. J Appl Clin Med Phys 2006;7:28e34. [8] Raffi JA, Davis SD, Hammer CG, et al. Determination of exit skin dose for I192r intracavitary accelerated partial breast irradiation with thermoluminescent dosimeters. Med Phys 2010;37:2693e 2702. [9] Lymperopoulou G, Papagiannis P, Angelopoulos A, et al. A dosimetric comparison of 169Yb and 192Ir for HDR brachytherapy of the breast, accounting for the effect of finite patient dimensions and tissue inhomogeneities. Med Phys 2006;33:4583e4589. [10] Yan X, Poon E, Reniers B, et al. Comparison of dose calculation algorithms for colorectal cancer brachytherapy treatment with a shielded applicator. Med Phys 2008;35:4824e4830. [11] Ahnesj€o A. Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media. Med Phys 1989;16:577e 592. [12] Mackie T, Scrimger J, Battista J. A convolution method of calculating dose for 15-MV x rays. Med Phys 1985;12:188e196. [13] Gifford KA, Horton JL Jr, Wareing TA, et al. Comparison of a finite-element multigroup discrete-ordinates code with Monte Carlo for radiotherapy calculations. Phys Med Biol 2006;51:2253. [14] Vassiliev ON, Wareing TA, McGhee J, et al. Validation of a new grid-based Boltzmann equation solver for dose calculation in radiotherapy with photon beams. Phys Med Biol 2010;55:581. [15] Rogers D. Fifty years of Monte Carlo simulations for medical physics. Phys Med Biol 2006;51:R287. [16] Taylor R, Yegin G, Rogers D. Benchmarking brachydose: Voxel based EGSnrc Monte Carlo calculations of TG-43 dosimetry parameters. Med Phys 2007;34:445e457. [17] Perez-Calatayud J, Ballester F, Das RK, et al. Dose calculation for photon-emitting brachytherapy sources with average energy higher than 50 keV: Report of the AAPM and ESTRO. Med Phys 2012;39:2904e2929. [18] Chibani O, Ma CC. HDRMC, an accelerated Monte Carlo dose calculator for high dose rate brachytherapy with CT-compatible applicators. Med Phys 2014;41:051712. [19] Williamson JF. Monte Carlo evaluation of kerma at a point for photon transport problems. Med Phys 1987;14:567e576.
398
Z. Tian et al. / Brachytherapy 15 (2016) 387e398
[20] Salvat F, Fernndez-Varena J, Sampau J. Penelope-2011, a code system for Monte Carlo simulation of electron and photon transport, OECD. Issy les Moulineaux, France: NEA Data Bank; 2011. [21] Badal A, Badano A. Accelerating Monte Carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit. Med Phys 2009;36:4878e4880. [22] Jia X, Gu X, Sempau J, et al. Development of a GPU-based Monte Carlo dose calculation code for coupled electronephoton transport. Phys Med Biol 2010;55:3077. [23] Hissoiny S, Ozell B, Bouchard H, et al. GPUMCD: a new GPUoriented Monte Carlo dose calculation platform. Med Phys 2011; 38:754e764. [24] Jia X, Gu X, Graves YJ, et al. GPU-based fast Monte Carlo simulation for radiotherapy dose calculation. Phys Med Biol 2011;56:7017. [25] Jia X, Sch€ umann J, Paganetti H, et al. GPU-based fast Monte Carlo dose calculation for proton therapy. Phys Med Biol 2012;57:7783. [26] Jia X, Yan H, Cervi~no L, et al. A GPU tool for efficient, accurate, and realistic simulation of cone beam CT projections. Med Phys 2012;39:7368e7378. [27] Jahnke L, Fleckenstein J, Wenz F, et al. GMC: a GPU implementation of a Monte Carlo dose calculation based on Geant4. Phys Med Biol 2012;57:1217. [28] Townson RW, Jia X, Tian Z, et al. GPU-based Monte Carlo radiotherapy dose calculation using phase-space sources. Phys Med Biol 2013;58:4341. [29] Tian Z, Shi F, Folkerts M, et al. A GPU OpenCL based crossplatform Monte Carlo dose calculation engine (goMC). Phys Med Biol 2015;60:7419. [30] Jia X, Ziegenhein P, Jiang SB. GPU-based high-performance computing for radiation therapy. Phys Med Biol 2014;59:R151. [31] Hissoiny S, D’Amours M, Ozell B, et al. Sub-second high dose rate brachytherapy Monte Carlo dose calculations with bGPUMCD. Med Phys 2012;39:4559e4567. [32] Hissoiny S, Ozell B, Despres P, et al. Validation of GPUMCD for lowenergy brachytherapy seed dosimetry. Med Phys 2011;38:4101e4107. Magnoux V, Hissoiny S, et al. Fast GPU-based Monte [33] Bonenfant E, Carlo simulations for LDR prostate brachytherapy. Phys Med Biol 2015;60:4973. [34] Woodcock E, Murphy T, Hemmings P, et al. Techniques used in the GEM code for Monte Carlo neutronics calculations in reactors and other systems of complex geometry. In Proceedings of the Conference on the Applications of Computing Methods to Reactor Problems. 1965. [35] Jia X, Yan H, Gu X, et al. Fast Monte Carlo simulation for patientspecific CT/CBCT imaging dose calculation. Phys Med Biol 2012; 57:577e590. [36] Agostinelli S, Allison J, Amako K, et al. In: W Barletta, GEANT4-a simulation toolkit. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors, and Associated Equipment, Elsevier, 2003. 506: p. 250e303. [37] Rasmussen BE, Davis SD, Schmidt CR, et al. Comparison of airkerma strength determinations for HDR 192Ir sources. Med Phys 2011;38:6721e6729. [38] Zhang M, Zou W, Chen T, et al. Parameterization of brachytherapy source phase space file for Monte Carlo-based clinical brachytherapy dose calculation. Phys Med Biol 2014;59:455.
[39] Richmond C. ICRP publication 23. Int J Radiat Biol Relat studies Phys Chem Med 1985;48:285. [40] Berger M, Hubbell J, Seltzer S, et al. XCOM: Photon cross sections database. NIST Stand Reference Database 8; Gaithersburg, MD: National Institute of Standards and Technology; 2013. € [41] Klein O, Nishina Y. Uber die Streuung von Strahlung durch freie Elektronen nach der neuen relativistischen Quantendynamik von Dirac. Z f€ur Physik 1929;52:853e868. [42] Hubbell J, Veigele WJ, Briggs E, et al. Atomic form factors, incoherent scattering functions, and photon scattering cross sections. J Phys Chem Ref Data 1975;4:471e538. [43] Bohr, N., Atomic physics and human knowledge. Hoboken, NJ: John Wiley & Sons, 1958. [44] Sempau J, Sanchez-Reyes A, Salvat F, et al. Monte Carlo simulation of electron beams from an accelerator head using PENELOPE. Phys Med Biol 2001;46:1163e1186. [45] Rogers D, Kawrakow I, Seuntjens J, et al., NRC user codes for EGSnrc. NRCC Report PIRS-702 (Rev. B), Missoula, MT: Northen Rockies Coordination Center, 2003. [46] Walters B, Kawrakow I, Rogers D. History by history statistical estimators in the BEAM code system. Med Phys 2002;29: 2745e2752. [47] Renaud M-A, Roberge D, Seuntjens J. Latent uncertainties of the precalculated track Monte Carlo method. Med Phys 2015;42: 479e490. [48] Bush K, Gagne I, Zavgorodni S, et al. Dosimetric validation of Acuros XB with Monte Carlo methods for photon dose calculations. Med Phys 2011;38:2208e2221. [49] Petrokokkinos L, Zourari K, Pantelis E, et al. Dosimetric accuracy of a deterministic radiation transport based I192r brachytherapy treatment planning system. Part II: Monte Carlo and experimental verification of a multiple source dwell position plan employing a shielded applicator. Med Phys 2011;38: 1981e1992. [50] Taylor R, Rogers D. EGSnrc Monte Carlo calculated dosimetry parameters for Ir192 and Yb169 brachytherapy sources. Med Phys 2008;35:4933e4944. [51] Taylor R, Rogers D. More accurate fitting of I125 and Pd103 radial dose functions. Med Phys 2008;35:4242e4250. [52] Failla GA, Wareing T, Archambault Y, et al. Acuros XB advanced dose calculation for the Eclipse Treatment Planning System. Palo Alto, CA: Varian Medical Systems; 2010. [53] Fonseca GP, Landry G, White S, et al. The use of tetrahedral mesh geometries in Monte Carlo simulation of applicator based brachytherapy dose distributions. Phys Med Biol 2014;59:5921. [54] D’Amours M, Carrier J, Lessard E, et al. Th-c-aud A-03: A new approach for afterloading brachytherapy inverse planned dose optimization based on the accurate Monte Carlo method. Med Phys 2008;35:2969. [55] Hedtj€arn H, Carlsson GA, Williamson JF. Accelerated Monte Carlo based dose calculations for brachytherapy planning using correlated sampling. Phys Med Biol 2002;47:351. [56] Rivard MJ, Melhus CS, Granero D, et al. An approach to using conventional brachytherapy software for clinical treatment planning of complex, Monte Carlo-based brachytherapy dose distributions. Med Phys 2009;36:1968e1975.