Introduction
Radiation therapy has a key role in treating different curable tumours, as well as tumours that can be inhibited or treated temporarily.Reference Haberer 1 The radiation therapy techniques have been evolved in recent years and new methods have been developed. However, all of these techniques aim at concentrating on the absorbed dose in a specific tissue and preventing the healthy, and specifically the sensitive, organs from exposure to radiation.Reference Owadally and Staffurth 2 Despite modern and advanced radiation therapy techniques, the control and treatment of benign tumours continue to be a major challenge in cancer treatment, necessitating the improvement of therapeutic techniques. In this regard, the introduction of new types of ionising beams is a promising achievement in the field of radiation therapy.Reference Walton 3
Cancer treatment with charged particle and high-energy beams of heavy ions, such as proton and carbon, is called ‘Hadron therapy’. Robert Wilson investigated the depth dose profile of proton beams and observed that these particles deposited the greatest amount of energy right before a peak called the Bragg peak, due to their specific physics.Reference Wilson 4 This high ionisation density at the end of the particle range results in greater dose deposition in deep tumours relative to the conventional radiation therapy techniques.Reference Kraft 5
As a result, the greatest advantage of ion beams over photon beams (i.e., inverse depth dose profile) is that they lead to better dose deposition in the intended region.Reference Amaldi and Kraft 6 On the other hand, due to the sharp fall-off at the end of proton beam range, treatment planning with these beams is very sensitive to range changes and range uncertainties. This is because this sharp fall-off is a great advantage of the proton beam because treatment planning is highly precise; hence there is no exit dose that can result in delivering more treatment dose to the tumour and at the same time less dose to the adjacent normal tissues and organs. This makes proton therapy more efficient comparing traditional radiation therapy. Treatment planning with proton beams needs very precise dose calculation.Reference Walton 3 Accuracy of clinical proton beam delivery depends highly on accurate modeling of the proton range in body tissues. Currently, the water phantom is the benchmark phantom most similar to the human body that is used for the quality control of the absorbed dose management. The basic data of the dose distribution are usually measured in a water phantom, which closely approximates the radiation absorption and scattering properties of soft tissues.
Since the nature of all interactions in the Monte Carlo simulation is based on theoretical models and empirical cross-section data for electromagnetic (EM) and nuclear interactions, this method is regarded as the most accurate simulator of particle interactions and calculator of dose in radiation therapy.Reference Paganetti, Jiang and Parodi 7 Accuracy of the Monte Carlo simulation depends on different parameters, such as the cutoff energy, step size and meshing. Moreover, the simulation accuracy depends on the applied physics. Geant4 (CERN, Geneva, Switzerland) toolkit allows users to employ different physics models.Reference Paganetti 8
Although the dose models in treatment planning systems (TPSs) completely depend on experiment, the Monte Carlo simulation has a significant role in this regard. The Monte Carlo particle transport codes, such as Geant4,Reference Agostinelli, Allison and Amako 9 MCNPXReference Waters 10 and FLUKA,Reference Ferrari, Sala and Fasso 11 are often used as the golden standard, based on which TPS is confirmed. As a result, the validation of nuclear interaction models used by these codes is very important.Reference Hall, Makarova and Paganetti 12
The Geant4 Application for Tomographic Emission (GATE) simulation code is a Monte Carlo simulation platform, which has been under development by the OpenGATE collaboration since 2001 and published in 2004. This code offers advantages of the valid physics models of GEANT4 simulation toolkit and can describe complex geometries.Reference Assie, Breton and Buvat 13 The GATE code was developed first for Positron Emission Tomography (PET) and Single Photon Emission Computed Tomography (SPECT) applications.Reference Jan, Santin and Strul 14 , Reference Santin, Strul and Lazaro 15 Since the development of version 6·0, some new GATE-specific tools with radiotherapy applications were added to it. However, there are scant studies into the GATE code applications in radiotherapy, specifically proton therapy. As a result, after careful validation studies into the GATE code, it can be used as a suitable tool for both imaging and radiotherapy studies.
This study aimed at validating the GATE Monte Carlo simulation code by obtaining the proton beam range with different energies in a water phantom. It also intended to investigate energy deposition and calculate proton beam range by measuring the Bragg curve of proton beams using the Monte Carlo calculations in the energy range used in the clinic.
Material and Methods
GATE Monte Carlo simulation
The physics instructions in the GATE code have been written based on the Geant4 library, which includes different physics models for all interactions in different particle and photon energies. Therefore, the GATE code uses the scripting mechanism instead of programming. To simulate the proton pencil beam profile in the current study, we used the GATE8 simulation code, which is based on the Geant4-version 10·3·3.Reference Agostinelli, Allison and Amako 9 , Reference Jan, Santin and Strul 14 , Reference Allison, Amako and Apostolakis 16
Simulation settings
Simulations were performed by utilising a PC-based including 12 independent Intel Xeon central processing units (CPUs) and 24 GB RAM.
Simulation geometry
The assigned geometry to the simulation code was a cubic water phantom with dimensions of 40×40×40 cm3. According to Figure 1, the axis of incident beam was selected to be z-axis. The source of proton emitting incident particles was located at 1 cm from one side of the phantom. The monoenergetic proton pencil beams collided with phantom in the energy range of 5–30 MeV with 5 MeV energy steps and for the energy range of 30–250 MeV with 10 MeV energy steps. The reason for choosing these energies was the coverage of the energy range used in proton therapy systems.
Figure 1 Geometry of the phantom, the cubic water phantom extending along the z-axis.
Pencil beam model
In addition, the pencil beam model (PBS) was used and the Gaussian pencil beam shape was considered. The proton beams within the proton energies, between 5 and 250 MeV, were used.
Since, in this study, the aim is to evaluate the validity of the GATE through the comparing the proton ranges from the simulation with the National Institute of Standards and Technology (NIST) data, the proton beams should have the characteristics in accordance with those NIST data are calculated with, that is, the monoenergetic protons without spread in energy.
The direction of the incident beam was along the z-axis, and the lateral directions were along the x- and y-axis. The ‘spot size’ is determined with σ x and σ y (standard deviation in x- and y-direction in isocenter) and was set at 3 mm. In addition, the realistic beam divergence was considered negligible (σ = 3 mrad). The monoenergetic proton pencil beam was used perpendicular to the xy-plane towards the end of the phantom. Proton histories for all energies were selected of order 106. The particle transport simulation continued until the statistical uncertainties in all energy ranges reached lower than 0·3%. Table 1 represents the properties of beam line.
Table 1 Properties of beam line
Physics-list selection
In Geant4, all physics processes are described with corresponding cross-section tables and sometimes there are different models for a certain physics process. All of these physics models and cross-section tables under 10 GeV are accessible in the GATE code. Specifically, new models have been developed to describe the Hadron interactions.Reference Jan, Benoit and Becheva 17 Currently, Geant4-version 10.3.3 provides different ‘packages’ as the physics list, allowing the user to select the most suitable model. To apply the suitable physics model for simulation, we investigated different models. To this end, the physics models in the reference physics list were used, including QGSP_BERT_HP_EMY (Bertini Cascade), QGSP_BIC_EMY (Binary Cascade) and QBBCReference Ivantchenko, Ivanchenko and Molina 18 recommended for medical applications and ion beam therapy. 19
The quark-gluon string precompound (QGSP) model was implemented to handle collision of high-energy hadrons. Two different shower models Binary Light Ion Cascade (BIC) and Bertini cascade (BERT) were used, which are responsible to track low energy region. QGSP_BERT_HP is identical to QGSP_BERT except that neutrons of 20 MeV and lower use the High-Precision neutron models and cross sections to describe elastic and inelastic scattering, capture and fission. For medical physics applications, the EM standard package with the Option 3 (EMY) parameters list is recommended by the GEANT4 EM Standard working group. Opt3 refers to processes recommends reference parameters to reach a high level of accuracy, it uses a set of EM processes with accurate simulation of gamma and charged particle transport.
QBBC (includes combinations of QGSP, BIC, BIC-Ion, BERT, … models) is recommended for applications where accurate simulation for low-energy transport of protons and neutrons is needed. It usually produces the best agreement in the energy range below 1 GeV for thin target experiments.
To obtain more precise results, all projected particles should be traced. In our simulation, proton, neutron, photon, electron, deuteron, triton, helium and alpha (secondary) particles were traced.Reference Jia, Hadizadeh and Mowlavi 20
Cut values (often referred to as ‘SetCut’) is a criteria specifying minimum secondary production ranges, internally converted into energy, to avoid the dependence on material. Below this threshold, the post step (discrete) processes do not occur and particles only undergo continuous processes, such as ionisations and excitations.
Results from the Monte Carlo simulation depend on the selected step size. It is worth noting that the selected step size should be small. As a result, there is a slight difference in the cross section at the beginning and end of the step. On the other hand, a large step size reduces the calculation time.Reference Paganetti 8 We set the Step max at 0·005 mm based on the Jia et al.’s study.Reference Jia, Romano and Cirrone 21
Actor selection
The required actors in this study were dose and fluence. To estimate the absorbed energy in the phantom and the incident proton fluence, the phantom meshing of 0·1 mm was conducted along the transmitted beam (Z). As a result, voxels were defined 40×40 cm2, perpendicular to the beam, and 0·01 cm along the beam direction.
Measuring range
In this study, R x was used to identify the specifications of the Bragg peak shape. R x Indicates the depth, at which the dose reaches x% of its maximum value and R peak is depth of the Bragg peak.
Due to the range straggling, not all protons with equal energy have the same range. As a result, the range needs to be defined for a beam of protons that results in a broadened Bragg peak or a spread-out Bragg peak (SOBP). Ideally, the range is defined at the position where the dose has decreased to 80% of the maximum dose, that is, in the distal dose fall-off.Reference Jia, Romano and Cirrone 21 The reason for this choice is the fact that for a monoenergetic proton beam, 80% fall-off position coincides with the mean projected proton range, that is, the range at which 50% of the protons have stopped. Furthermore, the 80% fall-off position is independent of the broadened energy beam.Reference Paganetti 22
GATE comparison with NIST (validation measurement)
To validate, the simulation results were compared to the NIST data. NIST has published data pertinent to the range and stopping power, assuming Continuous-Slowing-Down Approximation (CSDA).Reference Berger, Coursey and Zucker 23
Our validation included a comparison of depth dose profiles in water and fluence profiles in the overall energy range of beams used at the clinic. Accuracy of the model was determined by comparing 80% range of the end part of the dose curve for the dose profiles. Then, the output of the GATE code was compared to NIST data and simulation parameters were examined again. The relative percentage difference (RPD) and the least squares were used to validate and measure the agreement and overlapping of proton range results and CSDA range data in the NIST library.
In this study, we first obtained RPD of proton range with NIST data. Then, the range difference of all three physics with NIST data was compared to absolute uncertainty of the range. In addition, we compared the range difference of each physics with the absolute uncertainty of the range. Then, the least squares method was used to perform a more accurate investigation.
SetCut
To have a good agreement between NIST data and simulation results, it is crucial to select appropriate parameter values, particularly for the secondary production threshold (‘set cut’), and for the maximum allowed stepping (‘step max’).
In the next stage, the optimal SetCut in a selected energy (100 MeV) was addressed. To this end, we used QGSP_BIC_EMY physics because it had the least difference with NIST data. The SetCut values of 1 mm, 0·1 mm and 0·01 mm were selected and accuracy of the proton range calculation and simulation efficiency was addressed. To investigate the range accuracy, results were compared with NIST data and RPD was calculated. Then, the simulation range was obtained, using the following equation:
$$\eta {\equals}{1 \over {\sigma ^{2} {\times}T}}$$
where T is the simulation time and σ is the statistical uncertainty of simulation which is defined as mean uncertainty of all voxels between the entrance point and proton range.Reference Grevillot, Frisson and Zahra 24
Results
Figure 2 presents discrepancies between different R x from Monte Carlo calculation and R CSDA from NIST data. As it is seen, the least difference with NIST data is at 80% depth of maximum dose. As a result, we selected the depth of 80% of the maximum dose as the comparison reference.
Figure 2 Deviations of the characteristic beam range parameters from the Continuous-Slowing-Down Approximation (CSDA) range.
Evaluation of results obtaining from all three physics
Table 2 shows the results of the proton range for all three examined physics:
Table 2 Proton range results for all three physics
Results for proton range with different physics are presented in Figure 3a. For providing a better comparison, Figure 3b shows the difference between proton range with NIST data based on the incident energy in different physics.
Figure 3 (a) Comparing the results of all three physics obtained from the Geant4 Application for Tomographic Emission (GATE) code and National Institute of Standards and Technology (NIST) database. (b) The difference in the results from all three physics with NIST data.
The results of the least squares method for all three physics are shown in Table 3. Although no significant difference was observed between QGSP_BIC_EMY and QGSP_BERT_HP_EMY in terms of statistical uncertainty and difference between range of each physics in all energy values, Table 3 shows that applying the least squares method resulted in the least difference between QGSP_BIC_EMY results and NIST data. As a result, we recommended QGSP_BIC_EMY in the current study.
Table 3 The results of the least squares method
The results of QGSP_BIC_EMY physics
Results of this study were compared to the NIST data using the GATE simulation. The QGSP_BIC_EMY results showed a good agreement within 1% for energies higher than 20 MeV and within 0·1% for energies higher than 70 MeV.
Table 4 presents the results from the measurement of proton range in a water phantom at different energies in the QGSP_BIC_EMY physics list. As it can be seen, the proton range in water directly increases with the incident energy. In this study, the depth of 80% of the maximum dose was selected for comparison with NIST data (Figure 2). Results pertinent to the RPD with NIST data are presented in Table 4:
Table 4 Characteristic beam range parameters Rpeak, R90, R80, R10 in water and comparison with NIST data and along with relative percentage difference of 80% of the maximum dose
The results of the calculated depth–dose profile and proton beam fluence in the water phantom for QGSP_BIC_EMY physics list are shown in Figure 4.
Figure 4 (a) Proton depth–dose profile. (b) Normalised fluence profile of proton beam in the energy range of 5–250 MeV with QGSP_BIC_EMY physics.
The range–energy relation can be presented as follows by fitting the proton range results with QGSP_BIC_EMY physics in various energies (according to Figure 3a):
$$R\left( {{\rm g\,/\,cm}^{2} } \right){\equals}0.00257\left[ {E\left( {{\rm MeV}} \right)} \right]^{{1.73887}} $$
The results of dose distribution and proton fluence in water phantom are evaluated for each of examined energy as follow (Figure 5).
Figure 5 Dose profile and fluence of proton beam in different selected energies; (a) 50 MeV; (b) 100 MeV; (c) 200 MeV; (d) 250 MeV.
The results of SetCut
The results of SetCut are as follow (Table 5):
Table 5 Different range obtained for a 100 MeV proton beam at different SetCut values
Based on Table 5, the RPD is minimum at SetCut = 0·1 mm, at which the proton range results have the best agreement with NIST data. We also investigated the simulation efficiency and performed calculations for these three SetCut values, according to Equation (1):
According to Table 6, the simulation efficiency reduced with reducing SetCut. With the reduction of SetCut from 0·1 to 0·01 mm, the simulation time significantly increased and resulted in efficiency reduction; whereas, no significant increase in simulation time was observed with reducing SetCut from 1 to 0·1 mm and thus a slight reduction was observed in the simulation efficiency.
Table 6 Simulation efficiency for various SetCut
Discussion
This study aimed at validating the GATE simulation code and investigating the range of charged proton particles. These simulations were carried out in a water phantom and clinical energies between 5 and 250 MeV. We also showed the role of different physics in improving the Bragg peak depth.
It has been proven that the depth at which the proton fluence distribution halved to the baseline value R 50 (practical range), coincides with the depth at which the dose reduced to 80% of its maximum value d 80 (mean range).Reference Gottschalk 25 As a result, we compared the NIST data with 80% depth of the maximum dose.
It can be seen that in all energies, the difference in the range value obtained from all three physics and NIST data is higher than the statistical uncertainty.
According to Figure 3b, the difference between these physics is greater at higher energies. This is because, a change in physics results in a change only in the nuclear part, and nuclear interactions have a greater cross section at higher energies.
In all investigated energy ranges, the difference between range results obtained from QGSP_BIC_EMY and QBBC and QGSP_BERT_HP_EMY and QBBC is higher than the statistical uncertainty. As a result, QBBC physics has a greater deviation than the NIST data. On the other hand, comparison of QGSP_BIC_EMY and QGSP_BERT_HP_EMY physics showed that their difference in terms of obtained range is lower than the statistical range uncertainty. As a result, none of them outperformed the other. Therefore, QGSP_BIC_EMY and QGSP_BERT_HP_EMY proved valid in terms of range calculation accuracy.
To perform more accurate investigations, the least squares method was employed. By applying this method, the QGSP_BIC_EMY physics showed the lowest difference with NIST data. As a result, QGSP_BIC_EMY was recommended. It is emphasised that both QGSP_BIC_EMY and QGSP_BERT_HP_EMY physics were within the statistically uncertainty domain.
Moreover, the difference between the obtained range from QGSP_BIC_EMY physics and NIST data in all energies (except below 15 MeV) is lower than 1%. This difference was 1–3% in energies lower than 15 MeV.
In this study, the range–energy relation (Equation (1)) was obtained by fitting the QGSP_BIC_EMY physics data. A similar study was conducted by Bozkurt, using the MCNPX code.Reference Bozkurt 26 However, we investigated the overall energy range used in radiation therapy and obtained the fit model using a greater bunch of data.
When a high-energy proton radiation interacts with the patient body, secondary particles are produced during the nuclear interactions. Tracing all of these secondary particles until they reach zero-point energy significantly prolongs the simulation time. It is assumed in such codes as MCNPX that when the particle energy is minimised (cut-off energy), the overall kinetic energy of the particle is locally deposited, and the tracing and simulation processes are terminated. However, the Geant4 considers a range cut-off of the secondary particle production and the energy cut-off of secondary particle production. In other words, the particle is traced until it reaches the zero-point energy; however, the secondary particles are no longer produced below this energy cut-off. From this point onward, the initial particle energy reaches zero point based on the CSDA, thereby reducing the simulation time.Reference Jia, Romano and Cirrone 21 Therefore, the number of secondary particles increases with reducing the energy cutoff, thereby improving the simulation accuracy; however, it significantly affects the calculation time.Reference Zahra, Frisson and Grevillot 27
In the Monte Carlo simulation, since the secondary particle production cutoff can affect the energy loss and simulation results, tracing all particles should be avoided to improve the calculation efficiency. Energy of untraced particles should be locally deposited to ensure energy preservation.Reference Paganetti 8
The optimal SetCut was obtained at 0·1 mm, as not only the RPD was minimum, but also simulation efficiency significantly reduced with SetCut reduction from 0·1 to 0·01 mm. Although the efficiency is slightly better at SetCut = 1 mm, uncertainty is lower. Therefore, an increase in simulation accuracy at the cost of a slight increase in simulation time seems reasonable. As a result, the optimal SetCut was obtained 0·1 mm.
In general, this study investigated various physics lists in Greant4 toolkit. In this study, the QGSP_BIC_EMY physics showed the lowest difference with NIST data. As a result, to use the GATE code for radiation therapy, the QGSP_BIC_EMY physics and SetCut of 0·1 mm are recommended. Nevertheless, both QGSP_BIC_EMY and QGSP_BERT_HP_EMY physics are within the statistically uncertainty domain.
Conclusion
In this study, we validated the Gate Monte Carlo simulation code in proton therapy applications. This study first examined different physics lists and showed that the results obtained using QGSP_BIC_EMY physics are in the best agreement with NIST database data. Then, using this physics, we examined the SetCut values, and we found that the value of 0·1 mm was the optimal value for proton therapy.
Acknowledgements
This work was supported by Semnan University of Medical Sciences.
