Comparison of electromagnetic and hadronic models generated using Geant 4 with antiproton dose measured in CERN.

After proposing the idea of antiproton cancer treatment in 1984 many experiments were launched to investigate different aspects of physical and radiobiological properties of antiproton, which came from its annihilation reactions. One of these experiments has been done at the European Organization for Nuclear Research known as CERN using the antiproton decelerator. The ultimate goal of this experiment was to assess the dosimetric and radiobiological properties of beams of antiprotons in order to estimate the suitability of antiprotons for radiotherapy. One difficulty on this way was the unavailability of antiproton beam in CERN for a long time, so the verification of Monte Carlo codes to simulate antiproton depth dose could be useful. Among available simulation codes, Geant4 provides acceptable flexibility and extensibility, which progressively lead to the development of novel Geant4 applications in research domains, especially modeling the biological effects of ionizing radiation at the sub-cellular scale. In this study, the depth dose corresponding to CERN antiproton beam energy by Geant4 recruiting all the standard physics lists currently available and benchmarked for other use cases were calculated. Overall, none of the standard physics lists was able to draw the antiproton percentage depth dose. Although, with some models our results were promising, the Bragg peak level remained as the point of concern for our study. It is concluded that the Bertini model with high precision neutron tracking (QGSP_BERT_HP) is the best to match the experimental data though it is also the slowest model to simulate events among the physics lists.

Introduction

Radiotherapy during its history has been looking for two ultimate goals; maximizing the tumor dose while minimizing the dose received by normal tissues. To attain these goals new beam delivery techniques such as intensity modulated radiation therapy (IMRT)[1] or new beam modalities such as proton and heavier ions called hadron therapy were introduced.[2345] The strength of hadron therapy lies in their unique physical and radiobiological properties of these particles. Hadrons can penetrate the tissues with less lateral straggling and deposit the maximum energy just before stopping. This is the so called Bragg peak, which is the characteristic of hadron therapy. With the use of hadrons the tumor can be irradiated while the damage to healthy tissues is less than using megavoltage X-rays.

In 1984, Gray and Kalogeropoulos[6] proposed the idea of Antiproton cancer treatment to take the physical and radiobiological advantages, which come from its annihilation reactions. The stopping power of high-energetic (>>2 MeV) antiprotons in tissue is similar to that of protons. When a beam of antiproton transverses a medium most of the particles’ kinetic energy will be lost before they come to rest. But when the antiproton stops, the particles will annihilate with nearby nuclei and result in releasing 1.9 GeV of energy according to the well known mass and energy conservation relation. Most of the annihilation energy is carried away by pions, gamma rays and neutrons, but part of the annihilation energy (up to 30 MeV) will be still deposited locally by recoiling nuclear fragments with limited range. This effect doubles the local energy deposition compared to protons.[7] It is expected that recoiling nuclei will contribute with an enhanced relative biological effectiveness in the vicinity of the annihilation vertex. So a very favorable radiobiology is expected: In the plateau region, antiprotons behave mostly like protons, with well-known radiobiology and low Relative Biological Effectiveness (RBE) but in the peak region because of recoiling nuclei with higher Linear Energy Transfer (LET) the biological effects will be increased and thereby increasing the peak to plateau ratio of the biological effect even larger compared to proton beam of same energy.

Attempts to find out reliable measurements of radiobiological behavior of antiproton therapy are limited because of the lack of reliable dosimetry. Antiproton dose measurement is complicated because of the mixed particle spectrum especially in the peak area. The pulsed form of the antiproton beam also makes the measurement difficult. For this reason Holzscheiter used the biologically effective dose ratio (BEDR) to interpret a substantial increase in effective dose within the spread out Bragg peak (SOBP) of a CERN antiprotons beam.[8] For such situations, validated simulation toolboxes could be helpful to answer remaining questions about using antiproton in cancer therapy. Furthermore simulating the adverse effects of ionizing radiation at the cellular and sub-cellular scale could be very helpful.

Among all available simulation codes, Geant4 (GEometry ANd Tracking) is providing acceptable flexibility and extensibility, which progressively leads to the development of novel Geant4 applications in research especially modeling the biological effects of ionizing radiation at the sub-cellular scale. One of the major advantages of Geant4 is the possibility of applying different physics based on the experiment characteristics. Besides unlimited combination of physics list, there are different standard packages, under development. Those can cover most routine cases of radiation physics.[9] Meanwhile, Geant4-DNA is a project initiated in 2001 by Dr P. Nieminen at the European Space Agency for the development of a computing platform allowing an estimation of the biological effects of ionizing radiation using the Geant4 toolkit, which is based on the object-oriented programming (C++).[10]

The current study tries to calculate depth dose of a beam of antiproton corresponding to the CERN antiproton beam energy by Geant4 employing the entire standard physics list available. Hereby after verifying the code, a connection can be made between our works to those from Geant4-DNA project to take advantage of radiobiological simulation capabilities of the Geant4.

Materials and Methods

In this study measurement data collected at the Antiproton Decelerator in CERN by Bassler et al., as the experimental benchmark.[11] Geant4 toolkit version 9.4.6.p01 was used to reproduce the dose profile of a beam of a 502 MeV/c (Δp/p = 0.251 MeV/c) antiprotons which is the energy range used in the measurements.[11] For all calculations, the beam profile characteristics were defined similar to the CERN beam and those indicated in the measurements study.[11] The geometry applied in our simulation consists of a 20 × 20 × 20 cm3 water tank to simulate the experimental beam at CERN. Because we did not find any significant improvements on statistical precision above 100000 primary particles (events), this number of events was used to get better CPU performance for the rest of simulations. The scoring mesh was a box placed at the entrance of the water phantom and in the centre of the beam. The length of the scoring region through the beam direction was 20 cm to cover maximum range of antiproton (which is around 13 cm). The effective diameter of the ionization chamber in the measurements study was 4 cm to include all the antiproton particles and low energy fragments generated in the annihilation events in the active region, especially in the Bragg peak.[11] For this reason the lateral dimensions of the scoring region (perpendicular to the beam direction) was defined to be equal to 4 cm to fulfill the measurement requirements. This scoring mesh consists of 800 bins in the beam direction. Each one had the dimension of 4 × 4 × 0.025 cm. It makes the simulation resolution along the beam axis equal to 0.025 cm.

For the simulation different standard physics lists were applied. These Physics lists apply a string model for the modeling the interactions of high energy hadrons. The shower physics at the lower energies are handled by one of the intranuclear cascade models or the precompound model. Nuclear capture process was handled by either the Chiral Invariant Phase Space (CHIPS) model or one of the intranuclear cascade models, (Bertini or Binary). For each physics list, total dose deposition calculated by Geant4 was compared to the measured data to find out which model can accurately predict antiproton interaction. Most of the secondary particles fluxes were considered to investigate the effect of different physics models on the acquired results. In this way it would be possible to understand the origin of the discrepancies between the simulation results and measurements.

Results

All calculated depth dose curves were normalized to the dose at 60 mm beneath the entrance.[11] The standard list is mainly based on one of two main string models which are Quark Gluon String (QGS) or FRITIOF (FTF) like string model. The simulated depth dose curves using different FTF based model lists are shown in Figure 1. Two different shower models Binary (BIC) and Bertini (BERT) were used, which are responsible to track low energy region. Compared to the measured dose, the FTF models can predict the dose before the Bragg peak (plateau dose) and also after it accurately. These models were able to locate the Bragg peak in the correct position. But in the high dose region (Bragg peak) they were not adequately accurate to predict the dose level. It was found that for a basic string models (QGS or FTF) applying different cascade models (Binary and Bertini) could not affect the simulation results significantly [Figure 2]. Despite of the study by Ribon et al.,[12] FTFP_BERT_TRV list could not give the reasonable results to predict the at-rest event using hFritiofCaptureAtRest model [Figure 1].[12] Also FTFP_BERT was used to investigate the effect of different electromagnetic options. For electromagnetic interactions there are different electromagnetic physics (EM) option available other than standard one such as EMV and EMX as the fastest and most accurate options respectively. As it was expected there was no change on the results while using different EM options for both FTF and QGS based physic lists [Figure 2].

Figure 1

Calculated depth dose profile of 502 MeV/C antiproton beam applying Geant4 physics lists (solid line) compared to measurements (dots). Vertical and horizontal axis represents relative dose (arbitrary unite) and equivalent depth of water (millimeter) respectively

Figure 2

Calculated depth dose profile of 502 MeV/C antiproton beam applying different Quark Gluon String Models (solid line) compared to the measurements (dots). Vertical and horizontal axis represents relative dose (arbitrary unite) and equivalent depth of water (millimeter) respectively

As it is shown in figure 2, no significant difference on depth dose curves was found when implementing different physics lists based on QGS models (QGSP, QGSP_BERT, QGSP_BIC, QGS_BIC and QGSC_CHIPS). Furthermore QGSC_CHIPS could not calculate the depth dose in the low dose region before and after the Bragg peak. It may be due to the improper use of QGSC model, which is appropriate for high energy region. But it is not appropriate for antiproton inelastic interactions in low energy region. Despite of CHIPS, QGSC_BERT and LHEP_EMV, all other lists underestimate the dose in the peak region. The CHIPS model using G4QCaptureAtRest model is not accurate to calculate the falling part of the peak in spite of its reasonable predication of plateau dose and rising arm of the Bragg peak. This result did not agree with those by Kossove[13] in which simulation of at rest antiproton using CHIPS model had been given reasonable results. This may be due to different energy range of antiproton used in this study.

Figure 3 illustrates particle flux for secondary particles such as gamma, pions (π +, π−, and π0), kaons (K +, K−, K0), neutron, electron, positron, deuteron, alpha and triton. It was shown that different models calculated different particle flux, which resulted in greater discrepancy on the calculated dose consequently. Figure 4 illustrates the relative antiproton flux to the entrance primary particle flux calculated by different physics lists. As it was shown, there is the significant difference on antiproton flux using CHIPS and QGSC_CHIPS models which agree with the calculated dose showed in Figure 1. Also there was not any significant difference on particle flux calculated by QGSP_BIC, QGSP_BERT, QGSP_BERT_EMV and QGSP_BERT_EMX which agree with the Percentage Depth Dose (PDD) curves in Figure 2.

Figure 3

Particle fluxes calculated by CHIPS (upper left), QGSC_CHIPS (lower left), FTF_BERT_TRV (upper right) and QGSP_BIC (lower right) models normalized to the primary antiproton flux at the entrance level. Vertical and horizontal axis represents relative flux to (arbitrary unite) and equivalent depth of water (millimeter) respectively

Figure 4

Antiproton flux calculated by different physics list normalized to the primary antiproton flux at the entrance level. Vertical and horizontal axis represents relative flux (arbitrary unite) and equivalent depth of water (millimeter) respectively

Discussion

Overall, none of the standard physics lists was able to predict the antiproton PDD. As it was shown the level of the plateau dose and the position of Bragg peak reproduced in a correct manner by almost all the applied physics lists. The results of this study was different from those attained in Keyes’ study using older version of Geant4 (version 9.3.p01) in terms of the dose level before and after the Bragg peak.[14] Despite of better results in our study some cases but the Bragg peak level remained as the point of concern for both studies.

From the particle flux study, we can conclude that different physics related to track the annihilation products could be one of the major causes to explain the different dose calculation by different lists. Regarding the maximum allowable energy transfer to annihilation products, choosing the right physics model responsible for these particles is very important. Since most of the annihilation dose will take away by pions, careful consideration of the physic lists must be chosen to track these particles.

Since this study dealt with low energy physics, there was no significant difference related to the main high energy physics models. In QGSP the high energy interaction creates an exited nucleus, which is passed to the Low Energy Parametrised (LEP) to model the nuclear de-excitation. QGSC is the same as QGSP except applying CHIPS, modeling the nuclear de-excitation to improve simulation of the nuclear de-excitation part of the interaction. QGSP_BERT and QGSP_BIC are like QGSP, LEP Bertini cascade and Binary cascade below ~10GeV thus replace the use of the LEP. This study showed that there was no clear and significant difference in different lists between QGS based physics with those based on FTF physics.

Another finding of this study was that there is not any significant difference in calculated dose when applying different cascade models (Binary and Bertini). Cao et al.,[15] suggested that that inconsistency of antiproton simulation by Geant4 can come from the wrong annihilation cross section which does not agree with our findings. The difference in peak dose level can be conferred mostly from different physics than wrong cross section data. Also there is no same behavior among different lists using the same cross section data e.g. LHEP and LHEP_EMV. The LHEP models are based on parameterised modeling of hadronic interactions consists of the High Energy Parametrised (HEP) and Low Energy Parametrised (LEP) models. They cover all the long-lived particles at all incident energies. The LHEP is the fastest physics list, but its models are usually not very accurate. Better and slower models exist, but they do not apply to all particles at all energies. Simulation results based on LHEP lists and LEP model showed that this model is not even reliable to calculate the dose in low dose region. Also LHEP_EMV is overestimating the Bragg peak level while it does not include any dedicated at rest model for antiproton and Kaon minus.

Conclusion

This study presents the comparisons between measurements and many hadronic models in Geant4 for an antiproton beam with 502 MeV/c kinetic energy. This work concludes that standard physics lists to predict the depth dose of 502 MeV/c antiproton beam are reliable just before and after the Bragg Peak i.e., in low dose region. The Bragg peak is placed to the correct position by near all the applied physic lists. Furthermore the rising side of the Bragg peak in simulated curves is interestingly matching with the measured data.

This study concluded that the Bertini model with high precision neutron tracking (QGSP_BERT_HP) could the best to match the experimental data though it is one of the slowest physic model. Nonetheless there are some discrepancies when applying QGSP_BERT_HP which could be from the wrong physics model applied to pions which are responsible for transferring most of the annihilation energy resulting from nuclear fragmentation. Finally it was concluded that the antiproton annihilation model needs to be improved from the current version of Geant4.