Comparison of TG‐43 dosimetric parameters of brachytherapy sources obtained by three different versions of MCNP codes

Monte Carlo simulations are widely used for calculation of the dosimetric parameters of brachytherapy sources. MCNP4C2, MCNP5, MCNPX, EGS4, EGSnrc, PTRAN, and GEANT4 are among the most commonly used codes in this field. Each of these codes utilizes a cross‐sectional library for the purpose of simulating different elements and materials with complex chemical compositions. The accuracies of the final outcomes of these simulations are very sensitive to the accuracies of the cross‐sectional libraries. Several investigators have shown that inaccuracies of some of the cross section files have led to errors in  125I and  103Pd parameters. The purpose of this study is to compare the dosimetric parameters of sample brachytherapy sources, calculated with three different versions of the MCNP code — MCNP4C, MCNP5, and MCNPX. In these simulations for each source type, the source and phantom geometries, as well as the number of the photons, were kept identical, thus eliminating the possible uncertainties. The results of these investigations indicate that for low‐energy sources such as  125I and  103Pd there are discrepancies in gL(r) values. Discrepancies up to 21.7% and 28% are observed between MCNP4C and other codes at a distance of 6 cm for  103Pd and 10 cm for  125I from the source, respectively. However, for higher energy sources, the discrepancies in gL(r) values are less than 1.1% for  192Ir and less than 1.2% for  137Cs between the three codes. PACS number(s): 87.56.bg


I. INTRODUCTION
Brachytherapy has been a cornerstone treatment technique for the management of various malignancies such as gynecological, prostate, and breast tumors. In low-dose-rate (LDR) brachytherapy treatments, 125 I, 103 Pd, 192 Ir, and 137 Cs sources have been utilized for many years. Task Group 43 (TG-43) of the American Associations of Physicists in Medicine (AAPM) (1)(2) has recommended the dosimetric evaluation of low energy brachytherapy sources using experimental technique and/or Monte Carlo simulations. (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) A subsequent Working Group from AAPM has recommended the dosimetric evaluation of high-energy sources (16) using a similar method. Several different Monte Carlo codes have been utilized for the dosimetric evaluation of brachytherapy sources including MCNP4C2, (17)(18) MCNP5, (19) MCNPX, (20) EGS4, (21) A.1 125 I The IsoAid Advantage 125 I (model IAI-125A) source was introduced in the North American market in 2002. (9) Dimensions of this source for the Monte Carlo simulations are taken from the study by Meigooni et al. (9) The IsoAid Advantage seed contains a 3 mm long silver rod X-ray marker with a diameter of 0.50 mm. The silver rod is coated with a 1 μm thick layer of AgI containing 125 I. The active silver rod is encapsulated in 0.05 mm thick titanium (Ti) capsule with an outside diameter of 0.80 mm. The laser-welded end-caps have a maximum thickness of 0.10 mm. Two hemispheres are used for modeling the end-caps, one with 0.40 mm radius Ti hemisphere overlapped with a 0.35 mm radius air. The center of the air hemisphere was shifted by 0.05 mm relative to the Ti hemisphere. The overall source length is 4.50 mm and the active length is 3.0 mm. The cylindrical source element is free to move approximately 0.150 mm along the longitudinal direction of the source and 0.100 mm in the radial direction. The schematic diagram of the simulated source is shown in Fig. 1(a).

A.2 103 Pd
The 103 Pd source (Model Best 2335) was obtained from Best Industries. The dimensions of the source were taken from previous investigations. (11) The Best 2335 source consists of a cylindrical tungsten X-ray marker with a diameter of 0.5 mm and a length of 1.2 mm. There are three polymer spheres with 0.5 mm diameter located on either side of the X-ray marker. The centers of these pellets are located at distances of 0.9 mm, 1.5 mm, and 2.1 mm relative to the center of the seed. The relative chemical compositions of these spheres by weight are: 89.73% C, 7.85% H, 1.68% O, and 0.740% N. In addition, the mass density of the polymer is assumed to be 1.00 g/cm 3 . The polymer spheres are uniformly coated in 103 Pd. The spheres and the X-ray marker are double encapsulated in titanium shell with total thickness of 0.080 mm and outer diameter of 0.800 mm. The overall source length is 5.00 mm and the active length is 4.76 mm. The maximum possible displacement of a source sphere is 0.570 mm along the seed axis and up to 0.070 mm in the radial direction. The schematic diagram of the simulated source is shown in Fig. 1(b).

A.3 192 Ir
The source geometry and dimensions of the VariSource VS2000 were taken from the investigation of Angelopoulos et al., (13) Papagiannis et al., (14) and Llisp et al. (15) The VS2000 consists of two 2.50 mm long sources pellets. Each source pellet is made up of a 2.16 mm long cylindrical section with a 0.34 mm diameter and hemispherical ends with the same diameter. The two sources are placed laser drilled hole made at the end of a solid Ni/Ti wire (55.6% / 44.4%) with a diameter of 0.59 mm. The end of the wire was sealed with a 1 mm thick plaque. This source end was modeled as a hemisphere with its center shifted 3.205 mm relative to the center of the source. During the simulation, the length of the wire was extended by about 5.0 cm from the center of the source to include the impact of the source wire on the dose distribution. The active length of this source is 5.0 mm. The schematic diagram of the simulated source is shown in Fig. 1(c).

A.4 137 Cs
Spherical 137 Cs pellet sources of the LDR Selectron remote afterloading system, with 2.50 mm in diameter encapsulated in 0.5 mm stainless steel material, was simulated for this source. (8) The active cores of the sources contain spherical borosilicate glass with 1.50 mm diameter. For treatment of gynecological diseases, different combinations of active pellet sources and dummy pellets are moved through a plastic catheter and stopped at the end of the applicator by a stainless stopping screw.
In this study, a source train containing five active Amersham pellets was considered at the end of the cylindrical applicator to simulate a linear 137 Cs source, with 1 cm active length. The schematic diagram of the simulated source is shown in Fig. 1(d).

B. TG-43 dose calculations
The AAPM Task Group 43 (TG-43 )(1-2) introduced a formalism for calculation of dose distribution around a low-energy sealed brachytherapy source. Also, there is a report from the AAPM and ESTRO for dose calculation around photon-emitting brachytherapy sources with average energy higher than 50 keV. (16) We have used this report for dosimetric evaluation of 137 Cs and 192 Ir. This new report explains the differences of low-energy sources and high-energy sources for distances very close to the sources. At larger distances, there are not significant differences between the dosimetry of the sources using the two recommendations. Since the goal of this project is toward the dosimetry of the source for distances greater than 1 cm, the updated TG-43 protocol (2) (i.e., TG-43U1) will be used for all the sources.
According to the TG-43U1 formalism, (2) the dose rates around brachytherapy sources are obtained from: where S K is the air-kerma strength of the source, Λ is the dose rate constant, G L (r,θ) is the geometry function, g L (r) is the radial dose function, and F(r,θ) is the 2D anisotropy function. The (r 0 ,θ 0 ) is the reference point relative to the source with r 0 = 1 cm and θ 0 = π/2. The subscript "L" has been added in TG-43U1 (2) to denote the linear source approximation. The air-kerma strength, S K is the air-kerma rate K δ (d) multiplied by the square of the distance for energy levels greater than δ: L is the ratio of the dose rate at the reference point, P(r 0 ,θ 0 ) and S K . The unit of L is cGyh -1 U -1 : The geometry function accounts for the activity distribution within the source and the distance between the point of interest and the source. TG-43U1 defines the geometry function as where β is the angle in radians, subtended by the point of interest, P(r,θ), to the tips of the active length of a hypothetical line; L is the active length of the source. The radial dose function, g L (r), takes into account the effects of photon scattering and attenuation and is defined as follows: The 2D anisotropy function describes the variations in dose as a function of polar angle relative to the transverse plane and is defined as:

C. MCNP code
MCNP is a general purpose Monte Carlo radiation transport code that can be used for neutron, photon, electron, or coupled neutron/photon/electron transport. (30) Three versions of the MCNP code have been used in this study: MCNP4C2, MCNP5, and MCNPX. When a version of MCNP has been released, newly introduced features may be based on those from previous versions or these features may be novel to MCNP. However, the sources of these new features are not always explicit.

C.1 MCNP4C2
MCNP4C2 is the first major release of MCNP since version MCNP4B (February 4, 1997). (31) One of the major new features of MCNP4C2 relative to version MCNP4B is the use of ENDF/B-VI improvements. Monte Carlo simulations are, in general, only as accurate as the underlying cross-sectional libraries. (31) Several investigators (7,32) have discussed the problems with the photon cross-sectional data that were included in MCNP (31) simulations. The effect of these problems is especially visible in dose calculations involving low-energy photon emitting seeds. In 2002, the photon cross-sectional data included with the standard distribution of MCNP was completely revised (X-5 Monte Carlo Team, 2003). (31) For MCNP4C2 there are two photon interaction data libraries Labeled as ZZZ000.nnP with nn = 01 and nn = 02. (29) For the ZAID = ZZZ000.01P library, the photon interaction tables were created for elements from z = 1 through z = 94 except z = 84, 85, 87, 88, 89, 91, and 93. ZAID is Nuclide Identification number. This number is used to identify the element or nuclide desired. The form of the number is ZZZAAA.nnX, where ZZZ is the atomic number of the element or nuclide, AAA is the mass. These tables were based on evaluated data from ENDF (33) from 1 keV to 100 MeV. However, in the ZAID = ZZZ000.02P library, the values were a subset of the ZAID = ZZZ000.01P library with pair production thresholds added for the Storm-Israel (34) data. Data above 15 MeV for the Storm-Israel data and above 100 MeV for the ENDF data come from adaptation of the Livermore Evaluated Photon data library (EPDL (35) ) and are valid for energies up to 100 GeV. (29) For this project we used data library for nn = 02. (34)

C.2 MCNP5
The MCNP5 code incorporates an arbitrary three-dimensional configuration for materials in geometric cells bounded by first-and second-degree surfaces and fourth-degree elliptical tori. For photons, the code accounts for incoherent and coherent scattering, the possibility of fluorescent emission after photoelectric absorption, and the absorption in electron-positron pair production. (28) There are four photon transport libraries maintained by X-5 and distributed with MCNP: MCPLIB, MCPLIB02, MCPLIB03, and MCPLIB04. (35) MCPLIB04 was officially released in 2002. (29) The cross section, form factor, and fluorescence data are all derived from the ENDF/B-VI.8 data library that are derived from EPDL97. (34) Cross-sectional data are given for incident photon energies from 1 keV to 100 GeV. Fluorescence data are derived from the atomic relaxation data available in ENDF/B-VI.8, but use the storage and sampling scheme defined by Everett and Cashwell. (33) Electron interaction data tables are required both for problems in which electrons are actually transported, and for photon problems in which the thick-target bremsstrahlung model is used. Electron data tables are identified by ZAIDs of the form ZZZ000.nnE, and are selected by default when the problem mode requires them. There are two electron interaction data libraries: el01 (ZAID endings of .01e) and el03 (ZAID endings of .03e). (30) The electron libraries contain data on an element-by-element basis for atomic numbers from Z equal 1 to 94. The library data contain energies for tabulation, radiative stopping power parameters, bremsstrahlung production cross sections bremsstrahlung energy distributions, K-edge energies, and Auger electron production energies. It also contains parameters for the evaluation of the Goudsmit-Saunderson theory (29) for angular deflection based on the Riley cross-sectional calculation. This file includes Mott correction factors to the Rutherford cross sections used in the Goudsmit-Saunderson theory. The el03 library includes the atomic data of Carlson that was used in the density effect calculation. (30)

C.3 MCNPX
MCNPX is a general purpose Monte Carlo radiation transport code that tracks nearly all particles, as well as photons and electrons, at almost all energies. It is a superset of MCNP4C3 (36) and has many capabilities beyond MCNP4C3 such as heavy ion transport and long file names. It was first released to the public in 1999 as version 2.1.5. Many new tally sources and variance reduction options have been developed and new tallies have been created specific to the intermediate and high-energy physics ranges. (36) The MCNPX photon transport libraries are identical to the MCNP5 libraries that were defined in the previous section. The major difference of MCNP5 and MCNPX is their energy ranges. The energy ranges for MCNP5 are 1 keV to 1 GeV for electrons, and 1 keV to 100 GeV for photons.

C.4 Revised MCNP4C
Another code that has been generated in this project (referred here after as MCNP4C-revised) is the original MCNP4C2 code with the cross-sectional library taken from the MCNP5 and MCNPX library. Therefore, the cross-sectional library is changed from ENDF/B-VI version 5 to ENDF/B-VI.8. This adjustment is made in order to identify if the differences between MCNP4C2 and other codes are simply the differences in their cross-sectional library or if changes in the code are responsible for these differences.

D. Phantom geometry
To obtain the dose rate constant, Λ, radial dose function, g L (r), and the anisotropy function, F(r, θ), the simulations were performed in a spherical water phantom with 30 cm radius. Each source type was placed at the center of the phantom, and spherical tally cells were defined at different distances from the source center. The radius of the tally cells was 0.02 cm for distances less than 2 cm and 0.05 cm for distances greater than 3 cm. For each simulation, 10 9 particle histories were considered to ensure that the relative errors ( that represents the estimated relative error (%) at the 1 σ level) of MCNP simulations at all distances are less than 1%. To obtain air-kerma strength, each source type was placed at the center of an air phantom, the kerma rates were obtained in spherical tally cells, using F6-tally. The air-kerma rates, obtained for distances 1 to 10 cm, were multiplied by the square of the distance to obtain the air-kerma strength. The average values of the air-kerma strength were considered as the S K value for each source. In all simulations *F4 and F6 tallies were used. The results of the *F4 tally were multiplied by the values of mass absorption coefficients, to convert them to dose values. The TG-43 parameters (i.e., dose rate constant, g L (r), and F(r, θ)) were obtained according to Eq. (2) to (6). The uncertainty analysis of each parameter was obtained using the error propagation method. F(r, θ) is obtained at distances from 1 cm to 10 cm and the angles of 0° to 180° by 10° spacing. Radial dose functions have been obtained at distances of 1 cm to 6 cm of the source. The photon energy spectrum for 125 I and 103 Pd were taken from the TG43U1 report, (2) and the spectrum for 192 Ir was taken from national nuclear data center. (37) For 137 Cs, 662 keV monoenergetic gamma rays were considered. Table 1 shows the dose rate constant of 103 Pd, 125 I, 192 Ir, 137 Cs sources obtained by the four different MCNP codes and also the published values by different investigators. (9,11,13,38,39) For a better visibility of the differences between different Monte Carlo codes, these results were compared with the values obtained from the MCNP5 code (the most recently used code). Differences between the data from all the codes are excellent (0.1%) except MCNP4C2 that is 4.4% and 1.1% for 103 Pd and 125 I, respectively. It should be noted that the results of the simulations by the MCNP4C-revised code was similar to those of MCNP5 and MCNPX codes. Therefore, the main source of error for the dose rate constant of the 103 Pd and 125 I from MCNP4C2 code was the use of the cross-sectional file. The uncertainty of dose rate constant, originates from two parameters: S K and D · (r 0 , θ 0 ). The error of S K is calculated by the quadratic summation of the uncertainties from a point dose and the values calculated by deviations from the fitted line source. The uncertainty of D · (r 0 , θ 0 ) is calculated statistical fluctuations by the MCNP codes. Then, the uncertainty of the dose rate constant is calculated following the error propagation formula. (32) (11) 0.69 -0.6% TG-43U1S1 (38) 0.685 -1.3% Taylor & Rogers (39) 0.650 -6.3% MCNP4C2 0.935±0.6% 1.1% MCNPX 0.924±0.6% -0.1% MCNP5 0.925±0.6% ------125 I MCNP4C-Revised 0.928±0.6% 0.3% Meigooni et al (9) 0.99 7% TG-43U1S1 (38) 0.981 6% Taylor & Rogers (39) 0  Table 2 shows the simulated radial dose functions of the 103 Pd, 125 I, 192 Ir, 137 Cs sources as a function of radial distance, obtained by the four different Monte Carlo codes, using *F4 and F6 tally options. The percentage difference between the results of these codes relative to the MCNP5 data sets for 103 Pd and 192 Ir are shown in Fig. 2. These results indicate that for the F6 tally option, the percentage differences between the values obtained by the MCNP5 and MCNPX Monte Carlo codes are less than 6.7% for all distances. However, the percentage differences between the data obtained by the MCNP4C2 and MCNP5 codes increases with increasing distance from the source. The percentage differences between MCNP4C2 and MCNP5 codes at 6 cm distance from the 103 Pd source reach approximately 27% and 28% for F6 and *F4 tallies, respectively. However, the percentage differences between the MCNP4C2 and MCNP5 for 125 I, 192 Ir, and 137 Cs sources at a distance of 10 cm are less than 10.5%, 0.5%, and 1%, respectively. Interestingly, the results of the simulations by the MCNP4C-revised code are similar to those of MCNP5 and MCNPX codes. Therefore, the main source of error for the radial dose function of the low-energy sources from MCNP4C2 code was its cross-sectional   (11) for 103 Pd, Meigooni et al. (9) for 125 I , and Angelopoulos et al. (13) for 192 Ir.

C. Anisotropy function
The values of 2D anisotropy function, F(r,θ), for the 103 Pd, 125 I, 192 Ir, 137 Cs sources were obtained by MCNP4C2, MCNP5, MCNPX, and MCNP4C-revised Monte Carlo codes at distances of 1 cm to 10 cm. Figure 4 shows the comparison of the 2D anisotropy functions of these sources by the four codes using *F4 tally and also published data by other investigators. (9,11,13,38) Similar results were obtained by F6 tally. No significant differences between different codes had been observed.

IV. DISCUSSION & CONCLUSION
Four versions of Monte Carlo codes (MCNP4C2, MCNPX, MCNP5, and MCNP4C-revised) were used to obtain the TG-43 dosimetry parameters for low-and high-energy brachytherapy sources. To demonstrate the variation of the results between the codes, all the simulated parameters were compared with the data from the MCNP5 code. The results of these investigations for the high energy brachytherapy sources (i.e., 192 Ir and 137 Cs) indicated that the differences among the data from the four codes are less than 5.6%. However, in the case of low-energy brachytherapy sources such as 103 Pd and 125 I, the difference between the radial dose functions obtained by MCNP4C2 and those obtained by two other versions of the code (MCNPX and MCNP5) changes as a function of distance. The differences were found to be more than 27% for 103 Pd at r = 6 cm, and more than 10% for 125 I at r = 10 cm. The values of the dose rate constants for MCNP4C2 were approximately 4.4% and 1.1% larger than MCNP5 value for 103 Pd and 125 I, respectively. No significant differences were found on 2D anisotropy functions. Such differences show that the cross-sectional library used by the MCNP4C2 code for low-energy gamma rays were inaccurate. Therefore, the simulations were performed using the MCNP4C-revised code, which means that the MCNP4C2 cross section was changed to that used in the MCNPX and MCNP5 (ENDF/B-VI.8) codes. A comparison of results obtained with  (11) for 103 Pd, Meigooni et al. (9) for 125 I, and Angelopoulos et al. (13) for 192 Ir . the MCNP4C-revised and MCNP5 codes show that the average differences in the dose rate constant are 1.2%, 0.1%, 0.2%, and 0.1% for 137 Cs,192 Ir, 125 I, and 103 Pd respectively.
The maximum difference for *F4 and F6 tallies obtained by MCNPX and MCNP4C-revised in all points around the source were not more than 0.01%, 0.21%, 0.67%, and 0.41% for 137 Cs,192 Ir, 125 I, and 103 Pd respectively.
The results indicate that the only difference between MCNP4C-revised and MCNPX is the difference in their cross-sectional libraries and, as such, the results obtained with the two codes are similar. However, the difference between the results obtained with the MCNP5 and MCNP4C-revised codes is greater than the difference between the MCNPX and MCNP4Crevised codes.