Introduction

With the increased usage of Monte Carlo (MC) packages such as MCNP, Geant 4, EGSnrc, EGS5 and PENELOPE in a large array of conditions for dose calculation in radiation therapy and radioiodinetherapy, many studies have reported an excellent agreement between the derived results of these MC packages and experimental measurements [ 1 - 4 ]. In this regard, EGSnrc was shown to “pass the Fano cavity test at 0.1% level” [ 1 ]. Despite this accuracy, MC calculation is rarely used in radiation therapy clinics due to long computation time, which is needed to achieve sufficient statistical accuracy of the results. However, it has been stated that MC calculation can be considered as a gold standard of dose benchmarking especially in pencil-beam simulations of electron beams [ 5 , 6 ]. Moreover, it has been reported that MC can be accurately used in convolution-superposition techniques of photon beam calculations [ 7 ]. It should be noted that in MC calculation techniques, pre-computed data are employed to incorporate physics-rich elements in the dose calculation process. Nonetheless, in MC calculations a trade-off is needed between accuracy and computation time. Therefore, in complex heterogeneous geometries where enough computation time was not met, up to 10% discrepancy between the MC-derived results and experimental measurements have been observed [ 8 , 9 ].

In order to overcome the computation time limitation, a fast MC platform namely, Voxel-based Monte Carlo (VMC) [ 10 ] and Dose Planning Method (DPM) [ 11 ] were developed for the specific purpose of radiotherapy (RT) dose calculations. In such packages, simpler treatment of certain particle interactions is considered to speed up the calculations. This may lead to a decrease in the absolute accuracy of simulations; however, a speed of a factor of 50 can be achieved compared to general-purpose solutions.

Calculation time can be significantly reduced by applying advanced parallel process on the computer architecture of the MC package [ 5 ]. One straightforward way to ease this issue is to perform the computations in a parallel approach by using cluster architectures. Therefore, it is possible to increase the speedup factor by distributing the total computation load to the existing computing units [ 12 , 13 ]. Over the years, a significant number of reports have been published concerning the implementation of various MC and dose calculation packages on a variety of computing architectures, including Graphics Processing Unit (GPU), multi-core Central Processing Unit (CPU) and its clusters, and cloud computing [ 14 - 17 ].

The use of MC-based software seems to be essential in the modern RT practice especially in deriving patient dose, small field dosimetry, and quality assurance of treatment planning algorithms [ 18 - 20 ]. Therefore, GPU programming can be an essential tool in high-speed MC dose calculation of modern RT techniques namely, Intensity Modulated Radiation Therapy (IMRT), Stereotactic Body Radiation Therapy (SBRT), (Stereotactic Radiosurgery) SRS, etc.

The aim of this work was to create a novel GPU-based Fast Monte Carlo Photon Dose Code (MCPDC) for fast and accurate dose calculation in radiotherapy treatment planning. This package can be used for routine radiotherapy treatment planning systems with a high accuracy in complex heterogeneous geometries.

Material and Methods

Random Number Generators

In this analytical study, random number generators are critical parts of MC packages. Park and Miller random number series generator is more applicable due to simplicity and sufficient generation period [ 21 ]. The generator outputs are floating point numbers in the range of 0 to 1 [ 21 ].

In this work, the implemented random number generators of the created MCPDC have two independent random number subroutines for running on both host and device. The host routine has a unique seed number whereas the device routine has a vector (array) of seed numbers on any running threads for making differences between series in any threads.

Photon Simulation

Generally, analog simulation technique is widely used for the photon transport in MC packages. The photon is followed until it leaves desired medium in the simulated geometry or its energy becomes less than the pre-defined photon cutoff energy (Pcut).

Photon Interactions

Since the physic theory of photoelectric interactions is rather complex, MCPDC converts photon with certain energy to electrons with the same kinetic energy. Such an approximation is acceptable because the binding energy of electrons is significantly lower than that of photoelectrons.

In MCPDC, if pair production interactions occur, photons energy is randomly distributed between electron and positron and then, the electron is transported by its subroutine and positronis gone through Continuous Slowing Down Approximation (CSDA) model.

A condensed-history algorithm for electron transport requires complex physics of electron interactions with matter. As an alternative, MCPDC used a technique borrowed from fast electron Macro Monte Carlo (MMC), for electron transport algorithm [ 22 ]. This method avoids the use of complex electron transport algorithms.

On the first step, sphere geometry with fixed 0.05 cm in radius was simulated in EGSnrc code. Electrons with fixed energies were entered into the defined sphere. Output was defined based on any exited particle form sphere. Position, transport vector and energy of the exited particles were scored in a text file. Next, the text file was analyzed, and the needed parameters were derived for application in the electron subroutine.

GPU Simulation using CUDA

GTX 970 (GeForce GTX 970, NVIDIA Corporation GM204) was used as a GPU with acceptability for this work. The card was installed on a computer with Linux CentOS version 6.6 Operating System (OS). CUDA (Release 6.5, Version 6.5.16, CUDA compilation tools, NVIDIA Corporation) was installed on the OS as programing framework. The MCPDC was written in C programing language based on CUDA context.

Validation of MCPDC

MC-derived Validation

The MCPDC MC-derived validation was performed by EGSnrc (Version 3, National Research Council Canada, Ottawa, Canada) DOSXYZnrc. DOSXYZnrc is an EGSnrc-based MC simulation code for calculating dose distributions in a rectilinear voxel phantom [ 23 , 24 ].

For all simulations, the material data file named “700icru.pegsdata”, prepared in PEGS4 (Preprocessor for Electron Gamma Shower), was used. This data file has been previously validated in RT EGSnrc MC code simulations [ 1 , 23 , 24 ]. In the mentioned file, the threshold kinetic energy for secondary electron production is set on 189 keV. Therefore, using this file in EGSnrc code results in global ECUT and PCUT of 0.70 and 0.20 MeV, respectively. In addition, the default values of ESTEPE = 0.25 and S_max = 0.5 were used. The boundary- crossing algorithm in the EGSnrc/DOSXYZnrc was set on PRESTA-II.

On this step, poly-energetic pencil and rectangular collimated isotropic point source photon beams were simulated in MCPDC and DOSXYZnrc. For both beams, a tuned energy spectrum of the source (0.2 to 6.7 MeV) was obtained by analyzing the International Atomic Energy Agency (IAEA) published phase space file of the Siemens Primus linac (Siemens Healthcare Company, USA).

For MC-derived validation, virtual phantoms with different layers of water, bone, lung and air materials were simulated in both MCPDC and DOSXYZnrc. The anthropomorphic phantoms were also used for tuning the MCPDC, in terms of material heterogeneity compared to DOSXYZnrc code. The voxel size for this phantom was set as 0.001mm³.

DOSXYZnrc was run on a PC server, a part of a cluster server, with 2×16 AMD (Opteron 2.2 GHz) and 64 Giga Bytes of RAM that was installed at Isfahan University of Medical Sciences. It must be noted that, natively, “the DOSXYZnrc does not support multi-core architectures” and therefore has not been modified [ 24 , 25 ].

Experimental Measurement Validation

Percentage Depth Dose (PDD) and Lateral Profiles of the Siemens linac were measured using Photon Diode Detector (PFD3G Model, Sweden). The resolution of the measurements for the PDD and profile was 1 mm. The results of the measurements were used for the validation of the energy spectrum, which was derived from IAEA publications, and also the MCPDC and DOSXYZnrc.

A semi-automatic rotary asymmetric phantom was designed and produced using bone, lung and soft tissue equivalent materials including acrylic (ρ = 1.19 gr/cm³), polyurethane (ρ =0.2 gr/cm³) and Teflon (ρ =1.8 gr/cm³), respectively. The external and internal geometries of the produced phantom were designed based on the measured dimensions of the patients’ CT data at Department of Radiology, Isfahan University of Medical Sciences, Iran. A servomotor was used for phantom rotations at 90° intervals. For this step, Mapcheck (Sun Nuclear, Melbourne, Florida, United States) dosimeter was used for the experimental measurement of the dose distribution produced by water phantom.

Figure 1 shows the experimental setup of the work, including the designed anthropomorphic phantom, Mapcheck dosimeter and Primus linac head. The water phantom was also used for evaluating the MCPDC calculations in terms of media heterogeneity. For validation, a digital CT-based model of the produced phantom was generated using the CT-create program in DOSXYZnrc [ 24 ]. The generated digital phantom was used for MC-derived validation and also for MCPDC dose calculations in the experimental validation step.

Figure 1. Setup of the experimental measurement including the linac head, Mapcheck dosimeter and the designed and produced rotary phantom.

Since Map-check is a two-dimensional dosimetry tool, for the full validation of MCPDC, the produced phantom was exposed from 3 different directions. For each direction, the Map-check measurements in the form of dose matrix were compared with MCPDC calculations.

To avoid the variability inherent to the experimental dosimetry, ion chamber (Scanditronix / Wellhofer Farmer Type Chamber FC65-P, Sweden) measurements were performed for three independent experiments. Mean values and standard deviations were calculated, and used for the evaluation of the MCPDC calculated results. All experimental measurements were done in the Radiation Therapy section, Seyed al-shahada University Hospital, Isfahan University of Medical Sciences, Isfahan, Iran.

Statistical Data Analysis

Gamma index was used to compare MC-derived and experimental measurement data using a formula introduced by Low et al. [ 26 ]. It is a factor which accounts for both distances to an agreement (DTA) and dose difference [ 26 ], as follows:

${\Gamma }_r\left(D_c,r_c\right)=\sqrt{\frac{{(r_r-r_c)}^2}{\Delta d^2}+\frac{{(D_c-D_r)}^2}{\Delta D^2}}$

Where Δd and ΔD are DTA and dose difference, respectively. r_r is reference position and r_c is the position of the comparison point. D_r and D_c are dose at reference and comparing points, respectively. It should be noted that in case of Γ_r (D_c, r_c ) ≤ 1, the test is passed and otherwise it is failed [ 26 ].

Results

Using the innovative definition of photon and electron interactions, the average calculation time for the MCPDC was 5.4 sec for 5e7 source particle history, significantly less than the DOSXYZnrc which was 400 min.

Figure 2 shows the results of MCPDC simulation with fully tuned energy spectrum compared to DOSXYZnrc and experimental measurements in a water phantom. As can be seen from this figure considering the MCPDC and DOSXYZnrc, in the build-up region, a discrepancy of up to 2% of the maximum dose or 2 mm was found for most of the points (more than 80%) in water within 18 mm of the interface. Less than 1% point-to-point difference was observed for depths above 19 mm. Similar results were seen for MCPDC derived results compared to the experimental measurements.

Figure 2. Results of MCPDC simulation with fully tuned energy spectrum compared to DOSXYZnrc and experimental measurements in the water phantom.

Figures 3 and 4 show the MCPDC and DOSXYZnrc results including PDD and dose distribution, respectively, in the water virtual phantoms with a layer of bone and lung heterogeneity. Table 1 reflects the results of statistical analysis of the dose discrepancies between MCPDC and DOSXYZnrc at different depths of the virtual phantoms.

Figure 3. PDD of MCPDC and DOSXYZnrc in the water virtual phantoms with a layer of bone (a) and lung (b) heterogeneity.

Figure 4. Dose distribution of MCPDC and DOSXYZnrc in anthropomorphic phantom with a layer of bone (a) and lung (b) heterogeneity.

Figure 5 illustrates dose distribution of MCPDC and DOSXYZnrc on a CT slice of the anthropomorphic phantom.

Figure 5. Dose distribution of MCPDC and DOSXYZnrc on a CT slice of the anthropomorphic phantom.

Considering the phantom rotations at 90° intervals, experimental dose measurements using the Mapcheck and MCPDCMC-derived results are shown in Figure 6 (a, b and c). This figure demonstrates an overlay of 6 MV photon beam isodose contour maps the MCPDC and Mapcheck measured dose matrix. The reference for normalizations on these two sets of isodose lines was the central point of the Mapcheck dosimeter of the 10 × 10 cm² beam. Table 2 shows the results of statistical analysis of the dose discrepancies between MCPDC derived results and Mapcheck dose measurements on different steps of the rotation.

Figure 6. Results of MCPDC simulations and experimental measurements using Mapcheck dosimeter in different views, considering the phantom rotations at 0° (a), 90° (b) and 180° (c) of the produced rotary phantom.

Discussion

In this work, a novel GPU-based MC code (called MCPDC) is presented; it can be performed for RT plan dose calculation with a good accuracy (less than 1%). This newly developed GPU-based MC code implements the semi-empirical physical models, derived from EGSnrc, in the range of radiation therapy photon energies, which have been shown to demonstrate accurate results. In addition, MCPDC has the same capability of using CT-based phantom data as DOSXYZnrc and other relevant methods of dose measurements [ 27 ].

In order to validate the results of MCPDC, in addition to virtual inhomogeneous phantoms, a semi-automatic rotary asymmetric phantom was designed and fabricated. The results were compared with those of DOSXYZnrc code in the virtual phantom and for most of the voxels (>96%) less than 3% dose difference or 3 mm DTA was found (Table 1). For the produced phantom, compared to the Mapcheck dose measurements less than 3% dose difference or 3 mm DTA was observed (Table 2). It should be noted that in typical clinical use, “the fraction of points that exceed 3% dose difference and 3 mm DTA can be extensive” [ 26 ], so compared to other published reports, this criterion is acceptable for clinical evaluations. Low et al. have evaluated the gamma dose distribution comparison method for clinical applications [ 26 ]. They have shown that using 5% dose difference and 2–3 mm DTA is suitable for RT clinical plan evaluations [ 26 ]. In Task Group 120, it has been stated that to create the pass/fail acceptance criteria for the results of an array detector such as Mapcheck, used in this study, careful consideration should be given [ 28 , 29 ]. In this regard, pass rates of 90% of the evaluated points are suggested when using 3 mm DTA and 3% dose difference criterion [ 28 , 29 ]. This is in good agreement with previous results that easily passed the criteria since 97% of the points within this criterion [ 5 ]. Whereas, 70.3% of the points existed within a tighter criterion of 2 mm DTA and 2% dose difference, reflecting high accuracy of the MCPDC in dose calculation (Table 2).

Findings comparing MCPDC and DOSXYZnrc revealed that MCPDC is considerably faster by a factor of at least 4360 times for poly-energetic photon beams. It should be noted that such speed was obtained using GTX 970 graphic card. Using more efficient card may produce better performance in terms of running time. The results of running time here are comparable with results of History et al. that simulated a mono-energetic photon beam of 15 MeVon a GTX 480 card [ 25 ].

Despite the achieved accuracy, in order to reflect better results, MCPDC may be run on a high-performance GPU card of a modern workstation. However, this can be more sustainable in comparison with a large cluster server including a significant number of CPU cores (> 256 cores). In this regard, as an instance, the power consumption of a single modern workstation with a high-performance GPU card in the worst case is about 800 Watts, significantly less than any server in a cluster form.

Generally, GPU programing has more constraints than CPU programing such as total memory limitation and no access to other system components [ 30 ]. Total memory in GPU is limited to a few Gigabytes (4 Gigabytes in this project); however, for CPU memory, it is significantly higher than this value [ 30 , 31 ].

Very fast calculation speed and high accuracy in dose calculation might allow the MCPDC to be used in routine RT clinics as a central component of a treatment plan verification system, and as a dose calculation engine for MC-based planning. Furthermore, using novel semi-empirical physical models with detailed photon transportation modeling not only gives us higher confidence in our physical dose calculations but also will allow us to perform accurate GPU-based MC calculations.

Conclusion

In this work, a novel GPU-based fast Monte Carlo photon dose calculating method called MCPDC is presented for fast and accurate dose calculation in radiotherapy treatment planning practice. The MCPDC accuracy was evaluated with various methods, including Mapcheck measurement and DOSXYZnrc simulation in both homogeneous and heterogeneous phantoms.

MCPDC is currently being developed for electron dose calculation module and graphic user interface. In addition, optimizations toward using efficient variance reduction techniques are in progress. Moreover, future work on the applicability of the improved version of the MCPDC in transit dosimetry of Megavoltage CT is in process.

References

1. Rogers D W. Fifty years of Monte Carlo simulations for medical physics. Phys Med Biol. 2006; 51:R287-301. DOI | PubMed
2. Shahbazi-Gahrouei D, Ayat S. Comparison of three methods of calculation, experimental and monte carlo simulation in investigation of organ doses (thyroid, sternum, cervical vertebra) in radioiodine therapy. J Med Signals Sens. 2012; 2:149-52. Publisher Full Text | PubMed
3. Shahbazi-Gahrouei D, Saeb M, Monadi S, Jabbari I. Clinical implications of TiGRT algorithm for external audit in radiation oncology. Adv Biomed Res. 2017; 6:117. DOI
4. Shahbazi-Gahrouei D, Ayat S. Determination of Organ Doses in Radioiodine Therapy using Monte Carlo Simulation. World J Nucl Med. 2015; 14:16-8. Publisher Full Text | PubMed
5. Karbalaee M, Shahbazi-Gahrouei D, Tavakoli M B. An Approach in Radiation Therapy Treatment Planning: A Fast, GPU-Based Monte Carlo Method. J Med Signals Sens. 2017; 7:108-13. Publisher Full Text | PubMed
6. Jelen U, Sohn M, Alber M. A finite size pencil beam for IMRT dose optimization. Phys Med Biol. 2005; 50:1747-66. DOI | PubMed
7. Khosravi M, Shahbazi-Gahrouei D, Jabbari K, Nasri-Nasrabadi M, Baradaran-Ghahfarokhi M, Siavashpour Z. Photoneutron contamination from an 18 MV Saturne medical linear accelerator in the treatment room. Radiat Prot Dosimetry. 2013; 156:356-63. DOI | PubMed
8. Krieger T, Sauer O A. Monte Carlo- versus pencil-beam-/collapsed-cone-dose calculation in a heterogeneous multi-layer phantom. Phys Med Biol. 2005; 50:859-68. DOI | PubMed
9. Akmali Z, Shahbazi-Gahrouei D, Mosleh Shirazi M A, Baradaran-Ghahfarokhi M, Fallahian N, Sherkat S. Design and fabrication of a 4-dimensional of respiratory phantom for studying tumor movement in radiotherapy by MR imaging. J Isfahan Med Sch. 2015; 33:631-42.
10. Gardner J, Siebers J, Kawrakow I. Dose calculation validation of Vmc++ for photon beams. Med Phys. 2007; 34:1809-18. DOI | PubMed
11. Sempau J, Wilderman S J, Bielajew A F. DPM, a fast, accurate Monte Carlo code optimized for photon and electron radiotherapy treatment planning dose calculations. Phys Med Biol. 2000; 45:2263-91. DOI | PubMed
12. Wang H, Ma Y, Pratx G, Xing L. Toward real-time Monte Carlo simulation using a commercial cloud computing infrastructure. Phys Med Biol. 2011; 56:N175-81. Publisher Full Text | DOI | PubMed
13. Pratx G, Xing L. Monte Carlo simulation of photon migration in a cloud computing environment with MapReduce. J Biomed Opt. 2011; 16:125003. Publisher Full Text | DOI | PubMed
14. Zhou B, Yu C X, Chen D Z, Hu X S. GPU-accelerated Monte Carlo convolution/superposition implementation for dose calculation. Med Phys. 2010; 37:5593-603. Publisher Full Text | DOI | PubMed
15. Tian Z, Zhang M, Hrycushko B, Albuquerque K, Jiang S B, Jia X. Monte Carlo dose calculations for high-dose-rate brachytherapy using GPU-accelerated processing. Brachytherapy. 2016; 15:387-98. DOI | PubMed
16. Tian Z, Graves Y J, Jia X, Jiang S B. Automatic commissioning of a GPU-based Monte Carlo radiation dose calculation code for photon radiotherapy. Phys Med Biol. 2014; 59:6467-86. DOI | PubMed
17. Magnoux V, Ozell B, Bonenfant E, Despres P. A study of potential numerical pitfalls in GPU-based Monte Carlo dose calculation. Phys Med Biol. 2015; 60:5007-18. DOI | PubMed
18. Benmakhlouf H, Sempau J, Andreo P. Output correction factors for nine small field detectors in 6 MV radiation therapy photon beams: a PENELOPE Monte Carlo study. Med Phys. 2014; 41:041711. DOI | PubMed
19. Cranmer-Sargison G, Weston S, Evans J A, Sidhu N P, Thwaites D I. Implementing a newly proposed Monte Carlo based small field dosimetry formalism for a comprehensive set of diode detectors. Med Phys. 2011; 38:6592-602. DOI | PubMed
20. Mosleh-Shirazi M A, Karbasi S, Shahbazi-Gahrouei D, Monadi S. A Monte Carlo and experimental investigation of the dosimetric behavior of low- and medium-perturbation diodes used for entrance in vivo dosimetry in megavoltage photon beams. J Appl Clin Med Phys. 2012; 13:3917. DOI | PubMed
21. Press W, Flannery B, Teukolsky S, Vetterling W. Numerical Recipies in C: The Art of Scientific Computing. Cambridge Press: Cambridge ; 1988.
22. Neuenschwander H, Mackie T R, Reckwerdt P J. MMC--a high-performance Monte Carlo code for electron beam treatment planning. Phys Med Biol. 1995; 40:543-74. DOI | PubMed
23. Rogers D W, Faddegon B A, Ding G X, Ma C M, We J, Mackie T R. BEAM: a Monte Carlo code to simulate radiotherapy treatment units. Med Phys. 1995; 22:503-24. DOI | PubMed
24. Walters B, Kawrakow I, Rogers D W O. DOSXYZnrc Users Manual. Ionizing Radiation Standards. National Research Council of Canada: Ottawa K1A 0R6 ; 2013.
25. Hissoiny S, Ozell B, Bouchard H, Despres P. GPUMCD: A new GPU-oriented Monte Carlo dose calculation platform. Med Phys. 2011; 38:754-64. DOI | PubMed
26. Low D A, Dempsey J F. Evaluation of the gamma dose distribution comparison method. Med Phys. 2003; 30:2455-64. DOI | PubMed
27. Rezaee V, Shahbazi-Gahrouei D, Monadi S, Saeb M. Evaluation of error doses of treatment planning software using solid anthropomorphic phantom. J Isfahan Med Sch. 2016; 34:908-13.
28. Ezzell G A, Burmeister J W, Dogan N, LoSasso T J, Mechalakos J G, Mihailidis D. IMRT commissioning: multiple institution planning and dosimetry comparisons, a report from AAPM Task Group 119. Med Phys. 2009; 36:5359-73. DOI | PubMed
29. Low D A, Moran J M, Dempsey J F, Dong L, Oldham M. Dosimetry tools and techniques for IMRT. Med Phys. 2011; 38:1313-38. DOI | PubMed
30. Jason S, Edward K. CUDA by example: an introduction to general-purpose GPU programming. 2010.
31. Cook S. CUDA programming: a developer’s guide to parallel computing with GPUs. Elsevier: Amsterdam, Boston ; 2012.