Monte Carlo study of TG‐43 dosimetry parameters of GammaMed Plus high dose rate 192Ir brachytherapy source using TOPAS

Abstract Purpose To develop a simulation model for GammaMed Plus high dose rate 192Ir brachytherapy source in TOPAS Monte Carlo software and validate it by calculating the TG‐43 dosimetry parameters and comparing them with published data. Methods We built a model for GammaMed Plus high dose rate brachytherapy source in TOPAS. The TG‐43 dosimetry parameters including air‐kerma strength S K, dose‐rate constant Λ, radial dose function g L(r), and 2D anisotropy function F(r,θ) were calculated using Monte Carlo simulation with Geant4 physics models and NNDC 192Ir spectrum. Calculations using an old 192Ir spectrum were also carried out to evaluate the impact of incident spectrum and cross sections. The results were compared with published data. Results For calculations using the NNDC spectrum, the air‐kerma strength per unit source activity S K/A and Λ were 1.0139 × 10‐7 U/Bq and 1.1101 cGy.h−1.U−1, which were 3.56% higher and 0.62% lower than the reference values, respectively. The g L(r) agreed with reference values within 1% for radial distances from 2 mm to 20 cm. For radial distances of 1, 3, 5, and 10 cm, the agreements between F(r,θ) from this work and the reference data were within 1.5% for 15° < θ < 165°, and within 4% for all θ values. The discrepancies were attributed to the updated source spectrum and cross sections. They caused deviations of the S K/A of 2.90% and 0.64%, respectively. As for g L(r), they caused average deviations of −0.22% and 0.48%, respectively. Their impact on F(r,θ) was not quantified for the relatively high statistical uncertainties, but basically they did not result in significant discrepancies. Conclusion A model for GammaMed Plus high dose rate 192Ir brachytherapy source was developed in TOPAS and validated following TG‐43 protocols, which can be used for future studies. The impact of updated incident spectrum and cross sections on the dosimetry parameters was quantified.

heterogeneous medium, but they are currently regarded as only supplements to water-based dose calculation formalism. 3,4 Monte Carlo (MC) method is considered the "gold standard" for dose calculation in radiation therapy. Precise dose predictions can be achieved using MC method, especially in highly complex and heterogeneous environments such as human tissue.
Except for the general purpose MC codes used in brachytherapy dose calculations such as Geant4 5 and EGSnrc, 6,7 several MC dose calculation engines for brachytherapy applications have been developed including gBMC 8 and RapidBrachyMCTPS 9 based on Geant4, as well as Brachy-Dose 10 and egs_brachy 11 based on EGSnrc. MC dose calculations have been used as the ground truth in validating novel applications for brachytherapy. Intensity-modulated brachytherapy (IMBT) methods are commonly validated by MC dose calculations. 12 Skinner et al. 13 192 Ir source using Geant3 MC code, 18 where r is the distance from the center of the active source, and θ is the polar angle relative to the source longitudinal axis. r 0 and θ 0 denote the reference distance and angle, and are specified to be 1 cm and 90°, respectively.
The air-kerma strength S K , dose-rate constant Λ, radial dose function g L (r), 2D anisotropy function F(r,θ), and geometry function G L (r,θ) as well as their calculation methodologies were defined in TG-43U1 protocol. A fifth-order polynomial fit to the g L (r) data is commonly used.

2.A.2 | Air-kerma strength
Air-kerma strength S K is defined as the air-kerma rate _ K δ d ð Þ in vacuo at distance d located on the transverse plane of the source due to photons of energy greater than δ, multiplied by d 2 , and has units of cGy.cm 2 .h −1 (these unit combinations are also denoted by U), In this work, _ K δ d ð Þ per initial photon, _ k δ d ð Þ, is calculated using the following equation 24 where E i (MeV) is the midpoint of each energy bin, ϕ E i ð Þ is the photon fluence per initial photon at energy E i , μ en ρ E i ð Þ is the mass energy absorption coefficient at energy E i , and E is the bin size. The airkerma strength per unit source activity S K /A is then calculated from where 2.363 is the average number of photons emitted from one 192 Ir decay. 24 In this study, dry air was used in the calculations of air-kerma strength. 17 The photon fluence was calculated in a  Table 1.

2.C | Monte Carlo code and the simulation configuration
The TOPAS MC code 23 is an advanced and user-friendly extension to Geant4, [26][27][28] which can be used in simulation studies of various forms of radiotherapy. TOPAS version 3.4 was used in this work with Geant4 version 10.05.p02. The physics modules used in this work were "g4em-standard_opt4," "g4h-phy_QGSP_BIC_HP," "g4decay," "g4ion-binarycascade," "g4h-elastic_HP," and "g4stopping." The production threshold for all particles was taken as 10 keV following a previous recommendation. 20 The range cutoff for all particles was taken as 0.05 mm. The maximum step size varied from 0.05 to 1 mm for different voxel sizes.
We used two photon only 192 Ir spectra from National Nuclear Data Center (NNDC) 21 respectively. The statistical uncertainties of these parameters were smaller than 0.1%.

3.B | Radial dose function
The radial dose functions g L (r) calculated at radial distances from 2 mm to 20 cm are presented in Fig. 3. Two sets of data from previous works 7,18 are also shown for comparison reason. The g L (r) calculated using NNDC spectrum was on average 0.22% lower than that  Table 2 presents the obtained 2D anisotropy functions F(r,θ) using NNDC spectrum. The 2D anisotropy functions F(r,θ) calculated using two incident spectra at radial distances of 1, 3, 5, and 10 cm are shown in Fig. 4, as well as data from a previous published work. 20 Calculations inside the source or within 2 mm from the surface of the source are not accurate 7 ; thus, these results are not presented.

3.C | 2D anisotropy function
For these radial distances, the agreement between both our calculations and the reference data was within 1.5% for 15°<θ < 165°. The relative differences were larger (<4%) for θ values closer to 0°and 180°. The average deviations for these radial distances were −0.25%, −0.23%, −0.57%, and −0.38% for calculation using NNDC spectrum, and −0.33%, −0.31%, −0.46%, and −0.34% for calculation using Duchemen and Coursol's spectrum, respectively. The statistical uncertainties of F(r,θ) were smaller than 1% for 15°<θ < 165°, and within 5% for θ values closer to 0°and 180°. The deviation between calculations could be attributed to different source spectra and cross sections. A previous research 31 indicated that different 192 Ir source spectra can cause S k differences up to 2%. The relative difference was reported to be as high as 1% 32

| DISCUSSION
for Compton scattering attenuation coefficient in water between Geant4 "g4em-standard_opt4" physics model used in our study and XCOM photon cross sections used in Taylor and Rogers' work. In this work, the S k /A value calculated using new incident spectrum showed 3.56% relative difference from the reference, and that calculated using the old spectrum also used by Taylor and Rogers showed 0.64% difference. This indicated that the majority of the S k difference (2.90%) was attributed to different source spectra, while a small portion (0.64%) was caused by different cross sections.
The functions g L (r) and F(r,θ) near the source (r ≤ 2 mm) were not accurate thus not displayed because electronic equilibrium may not exist and the dose contribution from the beta spectrum of 192 Ir average energy of 181 keV is ignored using photon spectrum. 33,34 The use of new source spectrum caused average decrease of 0.22% for g L (r), which was observed from comparison between our calculations using different source spectra. But by comparing our work and Taylor and Rogers' both using the old spectrum, we observed that new cross sections seemed to have caused average increase of 0.48% for g L (r).
Less incident photons were simulated in F(r,θ) calculations because they were very time consuming, which caused the relatively high statistical uncertainties. Larger statistical uncertainties of F(r,θ) for θ values close to 0°and 180°were due to the small voxel size near the axis of the source. This can be improved in further studies by using more simulation histories. In this case, we conclude that the updated 192 Ir spectrum and cross sections did not result in significant discrepancies of F(r,θ); however, their impact cannot be quantified in this study.

DATA AVAILABILITY STATE MENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.