ORIGINAL ARTICLE Year : 2012  Volume : 8  Issue : 4  Page : 610618 Dosimetry parameters calculation of two commercial iodine brachytherapy sources using SMARTEPANTS with EPDL97 library Navid Ayoobian^{1}, Kamal Hadada^{2}, Barry D Ganapol^{3}, ^{1} Department of Nuclear Engineering, College of Mechanical Engineering, Shiraz University, Shiraz, Iran ^{2} Department of Nuclear Engineering, College of Mechanical Engineering, Shiraz University, Shiraz, Iran; Department of Aerospace Mechanical Engineering, University of Arizona, Tucson AZ ^{3} Department of Aerospace Mechanical Engineering, University of Arizona, Tucson AZ Correspondence Address: Aim: Simulating Many Accumulative Rutherford Trajectories Electron Photon and Neutral Transport Solver (SMARTEPANTS) is a discrete ordinates S _{N} Boltzmann/SpencerLewis solver that was developed during 19881993 by William Filippone and his students. The code calculates particle fluxes, leakage currents as well as energy and charge deposition for coupled electron/photon in xyz geometries both in forward and in adjoin modes. Originally, SMARTEPANTS was designed to utilize CEPXS crosssection library for shielding calculation in satellite electronics. The aim of this study was to adapt SMARTEPANTS to use a new photon crosssection library from Evaluated Photon Data Library, 1997 version (EPDL97) for intravascular brachytherapy ^{125} Isimulations. Materials and Methods: A MATLAB (MathworkNatick, Massachusetts) program was written to generate an updated multigroupLegendre crosssection from EPDL97. The new library was confirmed by simulating intravascular brachytherapy Best^{®} Model 2301 and Intersource ^{125} I dosimetry parameters using SMARTEPANTS with different energy groups (g), Legendre moments (L) and discrete ordinate orders (S). Results: The dosimetry parameters for these sources were tabulated and compared with the data given by AAPM TG43 and other reports. The computation time for producing TG43 parameters was about 29.4 min in case of g = 20, L = 7 and S = 16. Conclusion: The good agreement between the results of this study and previous reports and high computational speed suggest that SMARTEPANTS could be extended to a realtime treatment planning system for ^{125} I brachytherapy treatments.
Introduction The general introduction to Discrete Ordinates Method (DOM) During the past decade, numerous Monte Carlo (MC) codes have been available for medical dosimetry; however Discrete Ordinates Method codes (DOM) were just a few. The adaption of Acuros XB and Acuros BV, [1],[2] Boltzmann Transport Solvers, by Varian (Varian System, Palo Alta, CA, USA) after gaining FDA approval is a motivation for DOM codes to compete with MC in medical dosimetry. A major stumbling block in DOM dosimetry applications has been the complicated geometry of the human body, in which the phase space needed to be discretized very finely, requiring extensive memory and CPU. In today's computing, where powerful CPU and almost unlimited memory resources are available, DOM codes are able to use fine discretization of phase space and are much faster than MC and as accurate as MC codes. The DOM, SN , is a deterministic solution of the general Boltzmann/SpencerLewis equations governing particle transport. The method is well known and routinely used for reactor physics applications. However, DOM is rarely applied to medical radiation physics problems, where MC solutions of the Boltzmann equation are almost universally used. [2],[3],[4],[5],[6],[7],[8] The term DOM refers to the manner in which the angular domain is differenced. The energy variable is discretized by means of the multigroup energy approach using precalculated multigroup crosssection libraries. For discretization of the spatial variable, a variety of potential differencing schemes (e.g., diamond difference) are employed. After discretization, the original equation is replaced by a series of simultaneous linear difference equations that are solved iteratively. A general derivation of the multigroup energy approach and the spatial and angular discretizations can be found in the textbook by Lewis and Miller. [9] The solution of discrete ordinates equations approaches the exact solution of the Boltzmann equation as the space, energy and angle discrete in sizes approach differential size. The deterministic method has been applied to a variety of medical radiation physics problems, including brachytherapy calculations. [2],[3],[4],[5],[6],[7],[8] The aforementioned studies calculated dose distributions around a single source or applicator. Zhou [6],[7] parallelized their implementation of the deterministic method and calculated dose resulting from an 81 seed 125 I prostate implant. ATTILA [10] (Transpire Inc., Gig Harbor, WA, USA), an efficient, generalized geometry transport code developed at the Los Alamos National Laboratory, solves the threedimensional linear Boltzmann transport equation (BTE). ATTILA solves the BTE by discretizing all variables, energy (multigroup method), space (finite element) and angle (discrete ordinates or SN method) and then iteratively solving the differential form of the BTE for the phase space solution everywhere in the computational domain. Gifford [5] calculated the dose distribution around a Fletcher Suit Delclos ovoid loaded with Selectron 137 Cs pellets (Nucletron Trading BV, Veenendaal, The Netherlands) in a water phantom with a nine energy group crosssection set and S18 angular order. ATTILA is also benched against MCNPX for heterogeneities (ICBT tandem and ovoid applicators and their constituents) with a threegroup energy crosssection set and reduced angular order. [8] Simulating Many Accumulative Rutherford Trajectories Electron Photon and Neutral Transport Solver (SMARTEPANTS) [11],[12] is a diamond difference SN Boltzmann/SpencerLewis solver in xyz geometry for coupled electron/photon transport of up to four types of chargedneutral particles. The code has been mainly developed by William L. Filippone and his students at the University of Arizona during 1988/1993. The SMARTEPANTS approach is based on the GoudsmitSaunderson solution to the infinite medium SpencerLewis equation. The SMARTEPANTS code uses several variations of the SN diamond differencing algorithm to solve coupled electron/photon charged and neutral particle transport problems. Energy and charge deposition, particle fluxes and leakage currents are determined in xyz geometries. SMARTEPANTS were benched against CEPXS/ONEDANT [13] (a discrete ordinate code), ITS [14] (MC code) and analytical benchmarks of point and finite line beam sources in 3D geometry. [15] The code starts by reading the geometry, energy and angular mesh (SN ) provided in the input deck. At each energy, the code reads corresponding in group and group to group scattering crosssections from the crosssection libraries. CEPXS [13],[16] provides the scattering crosssections output and is rearranged into the proper crosssection library format usable by SMARTEPANTS. CEPXS is a multigroupLegendre crosssection generating code that was developed by Sandia National Laboratories. The multigroupLegendre crosssections produced by CEPXS enable coupled electronphoton transport calculations to be performed with the deterministic codes, over the photons energy range from 100 MeV to 1.0 keV. The SMARTEPANTS capabilities include its speed, accuracy and benchmark results. [17] Despite all appealing features, the code utilized an old crosssection library from CEPXS. In this paper, we present the improved photon crosssection library by implementing the early crosssection library, Evaluated Photon Data Library, 1997 version(EPDL97) [18] into the SMARTEPANTS code. EPDL97 is designed for use in photon transport calculations at the Lawrence Livermore National Laboratory. This library includes photon interaction data for all elements with atomic number between Z=1 (hydrogen) and 100 (fermium), over the energy range 1 eV, or threshold, to 100 GeV. Units are barns and MeV. The specific data in EPDL97 are the following: Coherent scattering, Incoherent, Total photoelectric reaction, Photoelectric reaction (by subshell), Pair production reaction and Triplet production reaction. EPDL97 does not include any photonuclear data. The updated code is used to calculate dosimetric characteristics of two commercially available iodine brachytherapy sources. The code's results are validated using the computational and experimental data for 125 I brachytherapy sources, available in the publication of AAPM Task Group No. 43 [19] and other investigator reports. The linearized BTE [20] Our study is focused on photon transport, and the DOM for photon transport is presented here in brief. For a homogeneous medium, the linear BTE can be written as: [INLINE:1] The first term on the lefthand side is interpreted as the rate of production of particles in phase space due to streaming while the second term represents the loss of particles due to all types of collisions. The first term of the righthand side is the rate of gain of particles in phase space resulting from scattering from all phase space elements at r, while the second term is the gain due to external source. Two physical assumptions are inherent in the above equation as follows: Scattering is assumed to be dependent only on the scattering angle [INSIDE:1] which is valid in isotropic media;To preserve linearity, macroscopic crosssections are assumed to be independent of the flux. The multigrouplegendre approximation [16] The multigroup method involves a discretization of the particle energy domain into energy intervals or groups: [INLINE:2] where E1> E2> E3 >…>E G+1 and EG+1 is the cutoff energy. By convention, the higher group numbers are associated with lower particle energies. The multigroup approximation is realistic only if the crosssections do not vary greatly in energy within a group. Hence, the structure of the energy grid can impact the accuracy of a prediction. The multigroup angular flux,[INLINE:3] is defined as: [INLINE:4] where w(E) is the multigroup weighting function. Because a weighting function must be arbitrarily chosen, a unique set of multigroup crosssections does not exist. In this work, w(E)=cg for Eg≥E≥Eg+1 , where cg is a groupdependent constant. Hence, the multigroupLegendre expansion coefficients become: [INLINE:5] Form factor and scattering function [21] In the case of Rayleigh scattering (coherent scattering), the photons are scattered by the bound electrons in a process where the atom is eventually neither excited nor ionized. The atom is therefore left in its ground state itself after the scattering process. The photon is scattered in some direction, making an angle ϑ with the incident direction with a primal energy. In a simplified treatment, called the form factor approximation, the differential scattering crosssection, σR , is given by the formula: [INLINE:6] where f (q, Z) is called the atomic form factor, r0 is the classical electron radius, Z is the atomic number, ϑ is the scattering angle and q is the corresponding momentum transfer. Sometimes, another quantity represented by x (momentum transfer) is used in place of q. This new quantity is defined by the equation: [INLINE:7] where λ is the incident photon wavelength. x is given in (Ε) 1 when the wavelength is given in (Å)1 Angstrom units. In the process of Compton scattering (incoherent scattering), the incident photon collides with an atomic electron. In this case, the photon loses a part of its energy to the electron. As a result, the electron is ejected from the atom. The differential crosssection for the Compton scattering of unpolarized photon from a free electron has been derived by Klien and Nishina: [INLINE:8] where α is the gamma energy expressed in unit of the rest energy of the electron. When the effect of electron binding is taken into account, the differential crosssection for the Compton scattering gets modified by a multiplicative function as follows: [INLINE:9] The function S(x, Z) is called the incoherent scattering function. Dose Calculation Formalism [19],[22] Based on the TG43 protocol, [19] the absorbed dose rate distribution around a sealed brachytherapy source for line source approximation can be determined using the following formalism: [INLINE:10] where Sk is the air kerma strength of the source, Λ is the dose rate constant, GL(r, ϑ) is the geometry function, gL(r) is the radial dose function and F(r, ϑ) is the anisotropy function. The above quantities are defined and discussed in detail in TG43. [19] The subscript ''L'' has been added to denote the line source approximation. r0 =1 cm and ϑ0 =0 o denote the reference distance and angle, respectively [Figure 1].{Figure 1} The dose rate constant was obtained from: [INLINE:11] GL(r, ϑ) is the geometry function (/cm 2 ) that describes the dose distribution in the absence of scattering and attenuation within the source itself and any surrounding materials. It is defined as: [INLINE:12] where β is the angle, in radians, subtended by the tips of the hypothetical line source with respect to the calculation point, P(r, ϑ), and L is the active source length [Figure 1]. The radial dose function, gL(r), describes the attenuation in tissue of the photons emitted from the brachytherapy source. The radial dose function is defined as: [INLINE:13] Materials and Methods In applying what has been discussed in the above sections, we wrote a MATLAB (MathworkNatick, Massachusetts) program to generate the multigroupLegendre photon crosssection from EPDL97. The program starts by reading the atomic numbers, weight fractions, energy group number, energy group type (linear or logarithm), maximum energy, cutoff energy, Legendre moment number (L) and density for each material. [Figure 2] shows the photon crosssection in the Graphic User Interface GUI control panel.{Figure 2} In the "MATERIAL INPUT" panel, the number of elements (one to five), element's name and weight fraction (varying between 0 and 1) for the corresponding element are inputted. The other input data such as energy group number (g), energy group type (li for linear and lo for logarithmic), maximum and cutoff energy, Legendre moment number (L) and density are interned in the "INPUT DATA" panel. The form factor and scattering function con cepts are used to generate differential multigroup scattering crosssection. As the "RUN" push button key is clicked, the following quantities are produced: coherent/ incoherent/ photoelectric/ pair production/ triplet production/ absorption/ total multigroup cross sections; multigroupLegendre selfscattering/ downscattering cross sections; and energy deposition. By clicking each push button key in the "Cross Section Plot" panel, the corresponding Figure will be plotted in the axes area, e.g. in [Figure 2], by clicking the "Incoherent" key the multigroup incoherent scattering crosssection for water in the energy interval of 0.001 to 2 MeV is plotted. The necessary crosssection files (called cpx.MG and cpx.CSD, original name by W. Fillopne [11],[12] ) are generated by the MATLAB program for SMARTEPANTS. For the iodine study, we used 5, 10, 15 and 20 energy group approximations (g) when the Legendre coefficients (L) were 5 and 7. Encapsulated 125 I is used in radiation brachytherapy for prostate cancer and brain tumors. The 125 I physical halflife is 60.159 days, and it emits gamma rays of 35.5 keV. The Best® Model 2301 125 I source is comprised of a doubleencapsulated titanium source surrounding a tungsten Xray marker coated with an organic carbon layer impregnated with 125 I. The source has a physical length of 5mm and an outer diameter of 0.8mm. The internal cavity of the source capsule has a diameter of 0.64 mm. The cylindrical tungsten Xray marker has a physical length of 3.7mm and a diameter of 0.25mm, and is coated with a 0.1mmthick organic matrix. The thickness of the coated organic matrix at each of the two ends of the marker was 0.15mm. The titanium wall thickness was 0.08mm. The thickness of the two ends of the capsule was 0.08 mm [Figure 3]. [23]{Figure 3} The InterSource 125 I has a physical length of 4.5 mm and an outer dimension of 0.81 mm. This source is manufactured by placing a 0.045mmthick Pt/Ir alloy (90% Pt and 10% Ir) Xray marker and three cylindrical bands of an insoluble organic matrix containing iodine between two hollow titanium cylindrical tubes. The inner and outer walls of the source are 0.04mm thick [Figure 4]. [24]{Figure 4} The cubic mesh structure is used to introduce two iodine sources (Best® Model 2301 and Intersource) and water phantom in SMARTEPANTS. In this meshing method, the volumes of the sources are exact. For example, [Figure 5] shows the top view of the Best® Model 2301 source meshing that was used in the code.{Figure 5} The source was simulated at the center of a spherical water phantom 15 cm in diameter and the dose rate distribution in water was calculated for different values of g (g=5, 10, 15, 20 and 25), L (L=5 and 7) and S (S=4, 8, 12 and 16). Sk was determined by calculating air kerma rate at 5 cm (to closely simulate the NIST measured value) distance and correcting for the inverse square of the distance to obtain the value at 1 cm. [23],[24] Results MATLAB software used to extract crosssections from the EPDL and prepare the input crosssection library (cpx.MG and cpx.CSD) to be used by the SMARTEPANTS code. As an example, the fivegroup scattering moment crosssections with L=5 for titanium tubes are listed in [Table 1].{Table 1} The doserate constant, Λ, for the Best® Model 2301 source and InterSource have been found to be 1.025 and 1.019, respectively, in case of g=10, S=16 and L=7. [Table 2] illustrates the calculated dose rate constants (in case of g=10 and different values of S and L) in water medium. The revisions proposed in various references for these sources as well as other commercially available sources were also considered in this table. Same values as [Table 2] but for g=5, 10, 15, 20, S=16 and L=7 are tabulated in [Table 3]. [Table 3] and [Table 4] illustrates the values of radial dose function for g=10, L=7 and S=12, 16. {Table 2}{Table 3}{Table 4} A comparison between the radial dose functions obtained in this work (for g=10, L=7 and S=16) and those given in the literature as recommended and published values (Nath et al., [19] Sowards et al., [23] Meigooni et al.[24] and Meigooni et al.[25] ) is presented in [Figure 6].{Figure 6} The anisotropy function, F(r, ϑ), of the sources in water phantom is calculated for g=10, L=7 and S=16 at 100 angle increments relative to the source longitudinal axis at distances of 1, 2, 3, 4, 5, 6 and 7 cm from the center of the seed using the SMARTEPANTS code. The calculated and recommended anisotropy functions are shown in [Table 5] and [Table 6] for the Best® Model 2301 and InterSource 125I, respectively. [Figure 7] illustrates the variation of F(r, ϑ) with respect to that for different r between 2 and 7 cm. {Figure 7}{Table 5}{Table 6} Finally, [Table 7] displays the computation time in minutes on an Intel Core i5 2.8 GHzprocessor in the case of InterSource, g=10, L=7 and S=12/16. The computation time for producing TG43 parameters by the MCNP5 MC code is reported to be 49 h for 9,500,000 source particles using a single Pentium 4 PC (2.4 GHz). [31]{Table 7} Conclusion SMARTEPANTS has been shown to have several capabilities that makes it an applicable SN code in photon radiation transport calculations. Despite all appealing features, due to its old crosssection library, the code has not been used for medical dosimetry. By implementing the use of an updated EPDL97 photon data library, which is currently used in many dosimetry codes, SMARTEPANTS could be validated and benched with current dosimetric parameters. Dose rate constants, radial dose and anisotropy functions of the two commonly used brachytherapy sources, Best® Model 2301 and Intersource 125 I, have been determined by different investigators using theoretical and experimental methods according to the TG43 original protocol. In this study, the simulations were performed in water phantom using the SMARTEPANTS code for different values of g, L and S. The dosimetry parameters of this source have been presented in this work and the results were compared with the previous investigations. Excellent agreement was observed between this study and the reference results. The calculated values of the dose rate constant, radial dose and anisotropy function were compared with previous investigations for different values of g, L and S. These calculated results for g=10, L=7 and S=16 are shown to be in good agreement with the previous studies [Table 2],[Table 3],[Table 4],[Table 5] and [Table 6] and [Figure 7]. The maximum percentage of deviation between the anisotropy function computed in the presentstudy and the reference data (Sowards et al., [23] Meigooniet al.[24] ) for the Best® Model 2301 and Intersource 125 I were observed to be 2% and 1.9% at r = 2 cm and ϑ = 6070°, respectively. The values of the dose rate constant, radial dose and anisotropy function of the sources could be used to determine the dose distribution in any point around the short line 125 I source according to the formalism of TG43. [Figure 6] compares the radial dose functions for the Best® Model 2301and Intersource 125 I with available data. For small values of r (r<1 cm), there are some differences between SMARTEPANTS and the recommended data. The differences could be explained by the mesh types used in different approaches. Because of source geometry, a cylindrical mesh would reduce the differences near the source. The good agreement between the results of this study and previous reports suggests that SMARTEPANTS can accurately reproduce TG43 dosimetry parameters for intravascular brachytherapy 125 I sources. According to the accuracy and computational speed, the code could be extended to a realtime treatment planning system for 125 I brachytherapy treatments. Finally, as [Table 7] implies, the computation time of SMARTEPANTS is much less than the current MC dosimetry codes for achieving the same accuracy. However, the code could be further optimized and accelerated for the speed considered necessary in a treatment planning software.[32] References


