Physics elements of an algorithm for brachytherapy dose calculation in homogeneous media for 192Ir source

Author’s Accepted Manuscript Physics Elements of an Algorithm for Brachytherapy Dose Calculation in Homogeneous Media for 192Ir Source Eshraq Ababneh, Saed Dababneh, Shada WadiRamahi, Jamal Sharaf www.elsevier.com/locate/radphyschem

PII: DOI: Reference:

S0969-806X(18)30183-X https://doi.org/10.1016/j.radphyschem.2018.04.004 RPC7811

To appear in: Radiation Physics and Chemistry Received date: 3 March 2018 Accepted date: 2 April 2018 Cite this article as: Eshraq Ababneh, Saed Dababneh, Shada Wadi-Ramahi and Jamal Sharaf, Physics Elements of an Algorithm for Brachytherapy Dose Calculation in Homogeneous Media for 192Ir Source, Radiation Physics and Chemistry, https://doi.org/10.1016/j.radphyschem.2018.04.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Physics Elements of an Algorithm for Brachytherapy Dose Calculation in Homogeneous Media for 192Ir Source Eshraq Ababneha , Saed Dababneha,∗, Shada Wadi-Ramahib , Jamal Sharafc a

Department of Physics, Faculty of Science, Al-Balqa Applied University, P.O. Box 2587, Amman 11941, Jordan. b Biomedical Physics Department, King Faisal Specialist Hospital and Research Center, MBC03, P.O. Box 3354, Riyadh 11211, Saudi Arabia. c Department of Physics, The University of Jordan, Amman, Jordan.

Abstract Unlike for external beam treatment, Model Based Dose Calculation Algorithms for Brachytherapy are still under development, and have been recently encouraged by the American Association of Physicists in Medicine. This delay is at least partially caused by the extreme complications associated with the multidimensional considerations underlying Brachytherapy dose calculations. These algorithms require detailed mathematical description, based on extensive radiation physics considerations, of the spatial and angular distribution of the points where the decay photon first interacts in matter, called TERMA points, along with the functions describing the radial and angular dose spread around each TERMA point, called Kernel. This work contributes to this global effort by providing these functions for the individual photon lines in the decay of the exactly modeled Nucletron microSelectron v3 192 Ir source. The underlying physics of all photon and electron interaction mechanisms has been extensively considered. Furthermore, the effect of the absorbing media has been also investigated. Ultimately, three-dimensional integration of dose Kernel contribution at a specified position, running all over TERMA points in the irradiated volume, has been conducted to generate the dose at that specific location. This methodology is much less time, data storage, and CPU consuming compared to repeating ∗

Corresponding author: Tel. +962-79-5606613 Email address: [email protected] (Saed Dababneh)

Preprint submitted to Radiation Physics and Chemistry

April 4, 2018

the extensive Monte Carlo simulations. The work is merely intended to provide a launching point and basic key physics elements for anticipated future quantitative clinical use. Keywords: Brachytherapy, HDR, Monte Carlo Simulation, Geant4, Dose PACS: 87.55.D-, 87.55.Gh, 87.55.K-, 87.55.ne, 87.53.Jw, 87.53.Bn 1. Introduction Monte Carlo based Brachytherapy dose calculation for 192 Ir sources has been repeatedly reported in recent litaerature. Dose-rate constant, radial dose functions, geometry factors, and anisotropy functions, utilized in the AAPM Task Group 43 dose calculation formalism, were calculated using Monte Carlo simulation (Angelopoulos (2000)). In another work, Melhus and Revard examined underlying characteristics in Brachytherapy dosimetry parameters for medical radionuclides using Monte Carlo methods, in which sources were modeled as unencapsulated point or line sources in liquid water to negate variations due to materials and construction (Melhus (2006)). Mikell and Mourtada tackled the issue of dosimetric impact of an 192 Ir brachytherapy source cable length modeled using a grid-based Boltzmann transport equation solver (Mikell (2010)). Van der Laarse et al. proposed a method to model 192 Ir wires of any length and shape (Van der Laarse (2008)). Taylor and Rogers presented results of EGSnrc Monte Carlo calculations of the dose distribution surrounding a high dose rate 169 Y b Brachytherapy source and 14 high dose rate and pulsed dose rate 192 Ir Brachytherapy sources (Taylor (2008)). The literature described above calculated dose in Brachytherapy based on the American Association of Physicists in Medicine (AAPM) Task Group 43 (TG-43 and TG-43U1) (Rivard (2007) and Rivard (2004)), which tabulate data points that require interpolations and extrapolations, while a model based on basic principles provides the possibility that dose can be calculated at any point. In addition, unlike the TG-43, model based calculations, after being established, can be expanded to include the effect of scatter and absorption from heterogeneities. The TG-43 also assumes an infinite medium, thus overestimates dose at interfaces. The current work thus aims at providing basic physics elements that can be further explored for developing model-based algorithms. The accurate and fast calculation of a three-dimensional dose distribution within the patient is one of the most central procedures in modern radiation 2

oncology. The conflict between high speed and high accuracy is one of the crucial challenges for the development of modern dose calculation algorithms. In the beginning, dose algorithms for high energy external photon beams were developed for the ultimate homogeneous patient; a patient completely consisting of water (Brady, 2006). However, to account directly for the underlying physical processes responsible for the energy deposition within the patient, one has to introduce models which describe explicitly some aspects of the related energy transport (Mayles, 2007). This is usually accomplished by the introduction of dose Kernels, which describe in different levels the energy transport and deposition caused by a defined set of primary photon tissue interactions. Model based algorithms for external photon beams have started since the 1980’s, and accuracy was developed slowly. For Brachytherapy, however, development of Model Based Dose Calculation Algorithms (MBDCA) has been recently encouraged by the American Association of Physicists in Medicine (AAPM) Task Group 186 (TG-186) in their report (Beaulieu, 2012). Several literatures discuss the available current dose algorithm models, such as Ahnesj¨o (1989), Carlsson Tedgren and Ahnesj¨o (2008), Carlsson Tedgren and Ahnesj¨o (2003) and Beaulieu (2012). The efficiency of the Collapsed Cone (CC) superposition algorithm has been optimized for low-energy sources by approximating the point Kernels as isotropic (Mayles, 2007). However, this is not a valid approximation, particularly for 192 Ir, since the photon angular distribution is more forwardly directed (Poon, 2009). The method is thus not yet feasible for 192 Ir treatment planning. Because of the previous limitations on using CC approximation, and the fact that it is a very slow methodology for dose calculation, the AAPM report underlines specific research that physicists should tackle, in order to provide “required data” to move forward with MBDCAs. This work has been thus motivated by these considerations. We aim at providing basic arguments, physics based calculations, and mathematical description for the ”TERMA point” (points of first decay photon interaction) and Kernel distribution functions for the individual photon energies resulting from the decay of 192 Ir. Accurate calculation of scatter dose requires methods capable of integrating dose contributions over 3D and 2D geometries. Monte Carlo MC simulation constitutes the only practical method for the calculation of energy deposition Kernels that can consider all interactions of importance for the transport of secondary particles (Ababneh, 2009). The methodology and 3

calculations adopted in this work heavily depend on data generated from validated extensive MC simulations. 3D integration of dose Kernel contribution at a specified position, resulting from TERMA points distributed all over the irradiated volume, has been conducted in this work. In summary, due to the lack of applied algorithmic models for Brachytherapy, this work employs extensive, validated and detailed Monte Carlo simulations to generate the dose Kernels. By fitting the energy deposition Kernels and TERMA points, a mathematical description is provided. This helps later on in the convolution of the TERMA point and Kernel over the irradiated volume, ultimately aiming at obtaining analytical description for the absorbed dose. 3D integration of dose Kernel contribution at a specified position, running all over TERMA points in the irradiated volume, has been conducted in this work to generate the dose at that specific location. This methodology is much less time, data storage, and CPU consuming compared to repeating the extensive MC simulations, and provides a major step towards a more direct algorithm. Methods for dose calculation should not neglect the effects of heterogeneities (Beaulieu, 2012). For this purpose, the methodology developed in this work can be used to calculate TERMA points and Kernels separately, taking into account different standard media. This study thus provides an important point of departure to study the effect of heterogeneity in the future. The exact decay scheme of the radioisotope is also included in this work, rather than approximations currently employed assuming monoenergetic sources. 2. Materials and Methods 2.1. Geometry Modeling and Simulation The GEANT4 code built in this work has been implemented on Linux operating system platform (Ubuntu 16.04 LTS). A detailed description of the source geometry modeling, physics and essential validation, is presented in an earlier work published by the authors (Ababneh et al., 2014). The geometry in this work, however, is modeled as a non-voxelated phantom sphere with diameter 240 cm instead of a cube, for symmetry considerations (see Fig. 1(a)). The Nucletron Classic High Dose Rate HDR 192 Ir source is placed at the center of the phantom sphere and modeled as a cylinder with a pure 192 Ir core of density 22.42 g/cm3 , encased in stainless steel. The internal construction and dimensions of the source are depicted in Fig. 1(b). 4

Although the full decay of the source is simulated, the analysis considers the most intense 192 Ir primary decay photons; namely the 295 keV, 308 keV, 316 keV, 468 keV and 604 keV photons. In addition to water, bone as another therapeutically relevant material has been also used as a “phantom” material.

(a)

(b)

Figure 1: (a) A spherical non-voxelized phantom of diameter = 240 cm. (b) Schematic diagram of the 192 Ir source.

Each particle in an event acquires an Identification Number ID. The simulation code stores the event, position, parent ID, track ID and the energy deposition in output ASCII files. Independent C++ (compiler for C++ is the GNU g++ version 4.8.2) and ROOT (version 5.34) codes have been especially written to analyze the stored data. The analysis codes also generate plots in order to investigate dose, Kernel and TERMA points by using ROOT’s different draw and fitting options (ROOT, 2013). 2.2. TERMA Point Calculation In addition to the standard concept of the total energy released in material (TERMA), this work introduces the concept of the TERMA “point” (XT , YT , ZT ) as the point in the phantom where the decay photon first interacts by e.g. Compton scattering or photoelectric effect. By using track information, we extract TERMA point information from the primary decay gamma (parent ID = 1). After storing TERMA point information, analysis proceeds to study the behavior of the TERMA point as a function of the radial distance (r) from the origin point of the phantom where the source 5

is centered, together with the relevant angular distribution, θt and φt . The next step is to construct a mathematical function that has the best fit to the series of TERMA points. 2.3. Kernel Calculation Once energy is transferred to the medium starting at the TERMA point, the liberated electrons continue to travel and interact within the medium. All generated electrons need to be accounted for. The Kernel function describes the energy deposited by these electrons (or further transferred by any associated bremsstrahlung) around the relevant TERMA point. A comprehensive analysis model has been implemented for this purpose, which is outlined in what follows. 2.3.1. Kernel Electrons The TERMA point vector is defined as the vector pointing from the phantom origin to the TERMA point. The transported electron position relative to the parent TERMA point position (XK , YK , ZK ) is determined, so that the coordinate origin for these electrons becomes the unified TERMA point as follows: XK = XElectron − XT , YK = YElectron − YT , ZK = ZElectron − ZT ,

(1)

where XElectron , YElectron , ZElectron are in the original global coordinate system. Then, in order to obtain the distribution of Kernel electrons relative to their TERMA point vector, the polar angle θ and the azimuthal angle φ are sampled. Through multiple Coulomb interactions, the electron energy is deposited locally with a fraction being carried away by photons created in bremsstrahlung events. The produced Kernel electrons thus include; electrons resulting from photoelectric effect, liberated electrons by (multiply) Compton scattered photons, and “remote” electrons resulting from a bremsstrahlung photon interaction in the medium. 2.3.2. Non-Kernel Electrons In the analysis procedure, a special file (Depos.dat) contains all electrons generated from each run, regardless of the creator process (including decay 6

electrons). Once Kernel electrons are determined, the remaining entries in the Depos.dat file are categorized as Non-Kernel electrons. These electrons include: • Decay electrons, with parent ID = 1, originating inside the iridium core and occasionally ending in the phantom volume close to the source. • Electrons in the phantom related to decay electrons. The decay electrons produce bremsstrahlung inside the iridium source, thus producing photons that reach the phantom volume and interact via photoelectric absorption or encounter multiple Compton events until losing all energy with the photoelectric process. Obviously, these phantom electrons are classified as Non-Kernel because they are not connected to a TERMA point. 2.3.3. Kernel Fitting The Kernel spatial distribution exhibits spikes at short distances near the TERMA points. A noticeable single spike within r = 1 mm from the TERMA point was recorded for lower energy decay photons, namely the 295 keV , 308 keV , 316 keV and 468 keV . At higher energies, like 604 keV , the spike was within r = 2 mm. This variation is due to the electron rangeenergy dependence, as will be elaborated in more detail as part of future work. During spatial Kernel fitting, these Kernel spikes within r = 1 mm and r = 2 mm, for lower and higher energy decay photons, respectively, have been excluded from the fitting and then added in a separate process for the purpose of comprehensive dose calculation. 2.4. Dose Calculation By using a special script in combination with TERMA point data file supplied as input, dose due to secondary electrons is calculated at a certain point i in a 1 mm3 volume. The dose calculation follows the following main steps: • Calculate the radial and angular position of i with respect to each TERMA point in the space (Xi , Yi , Zi ), and then supply them as input to the Kernel function.

7

Figure 2: Multiple integral over the surface of a sphere volume element of 1 mm3 .

• Three-dimensional integration is carried out in spherical coordinates for the space around i in a 1 mm3 volume (Fig. 2). Generally, a triple integral is converted from rectangular to spherical coordinates, and carried out using the appropriate limits of integration for the variables r, θ and φ. The Kernel dose contribution is calculated for all the five main 192 Ir decay energies from Kernel electrons in a volume equal to 1 mm3 (Eq. 2). 1 R r+0.5 R θ+ sin−12 ( r1 ) R φ+ 2rsin(θ)

r−0.5

DoseKernelf unction =

θ−

1) sin−1 ( r 2

1 φ− 2rsin(θ)

K(r, θ)drdθdφ . (2)

2∗π

• The excluded spikes of first Compton electrons at close distances from their related TERMA points are then added to the dose as follows: – Calculate the number of TERMA points per mm3 for each incident photon energy individually using the TERMA point function: −1 ( 1 ) rt 2 sin−1 ( r1 ) t θt − 2

R rt +0.5 R θt + sin rt −0.5

T ERM A P oints =

1 R φt + 2rt sin(θ t) 1 t sin(θt )

φt − 2r

T (rt , θt )drt dθt dφt .

2 ∗ π ∗ DecayEnergy (3)

– Determine the spike contribution for each TERMA point within the selected interaction site depending on the decay photon energy. 8

– Multiply the previous two quantities to obtain the total energy deposition in a 1 mm3 volume. • Finally, calculate the final Kernel deposited energy at this part of the phantom as the sum of all Kernel function contributions added to the spike contribution. The obtained results (in MeV/mm3 ) are converted to dose (in Gy) using a conversion factor equal to 1.6 × 10−6 (Ababneh et al., 2014), taking into account the volume size and the material characteristics (water in this case). • Calculate the dose at user selected sites. A numerical comparison is then provided in order to examine the consistency between the calculated dose using this approach and the simulated data that represent all decay line intensities. 3. Results and Discussion 3.1. TERMA Points In this section, “TERMA point” distributions are considered and examined for the most intense 192 Ir primary decay photons; namely the 295 keV, 308 keV, 316 keV, 468 keV and 604 keV photons. The TERMA point radial distribution is an exponentially decaying curve, taking into account that the primary intensity is being attenuated downstream because of the increasing photon interactions. The results are depicted in Figure 3(a). A “TERMA point ray or vector” connects the phantom origin to the TERMA point. It sweeps out an orbit with varying circumference with theta (being measured relative to the Y axis). Consequently, the number of TERMA points sharply increases until it peaks around 90◦ (Figure 4(a)). The additional attenuation at 90◦ due to the capsule head of the source placed along the Z axis causes a slight bulge in the peak. It is worth mentioning here that this angular behaviour, as discussed above, depends on the circumference of the circle at a given theta. Since the final dose will be calculated in a small volume at a given location in the phantom, this implies the need to modify the theta dependence in the final calculation. This issue is ultimately taken into consideration for dose calculation by dividing the relevant fitting equation by sin(θ). On the other hand, the TERMA point distribution is in principle independent of the azimuth angle. However, two clear features around the azimuth 9

295 keV 308 keV 316 keV 468 keV 604 keV

Normalized TERMA

Normalized TERMA

1

10 −1

1

Fitting Curve 316 keV

0.8

0.6

10 −2 0.4

0.2

0

100

200

300

400

500

600

700

800

900 0 0

Radius (mm)

200

400

(a)

600

800

1000

1200 Radius (mm)

(b) Normalized TERMA

Water-Bone TERMA 1 316keV-Water 316keV-Bone

10−1

10−2

10−3 0

200

400

600

800

1000 1200 Radius (mm)

(c) Figure 3: (a) The normalized TERMA point radial distributions for the most intense 192 Ir primary decay photons. (b) The corresponding fitted radial distribution for 316 keV photons. (c) The normalized TERMA point radial distribution for 316 keV photons in water compared to bone.

angles 90◦ and 270◦ are due to absorption at the source capsule heads in the positive and negative Z axis, respectively (Figure 5). This bulge is relatively stronger than in the θ distribution since only small fraction of the θ = 90◦ cases are opposite to the capsule head on the Z axis. 3.1.1. Radial Fitting The radial behavior of TERMA point distribution in water phantom exhibits exponential decay, as represented by: T (r) = p0 exp(−p1 r),

(4)

where r is the TERMA point radial distance (mm) from the phantom origin, and the two coefficients (p0 , p1 ) are the fitting parameters. An example is depiced in Figure 3(b). 10

Normalized TERMA

295 keV

1

308 keV 316 keV 468 keV 604 keV

0.8

0.6

0.4

0.2

0 0

20

40

60

80

100

120

140

160

180

Theta

TERMA points x 0.316 MeV

(a) ×103 Fitting Curve 800

TERMA Points

700 600 500 400 300 200 100 0

0

20

40

60

80

100

120

140

160

180 Theta

(b) Figure 4: (a) The normalized TERMA point angular distributions for the most intense 192 Ir primary decay photons. (b) The fitted TERMA point angular distribution for 316 keV photons (number of primary decays = 1 × 108 ).

Based on the fact that Compton interaction depends on the electron density, bone exhibits stronger attenuation than water. This is reflected on the TERMA point distribution behavior, where energy released in bone is stronger than in water. An example is illustrated in Fig. 3(c) and Table 1 for 316 keV photons. The fitting parameters for the other 192 Ir photon energies in water phantom are listed in Table 2. It is worth mentioning here that throughout the fitting processes, non of the parameters was fixed, and consequently all parameters were subject to 11

Normalized TERMA

1

0.8

0.6

0.4

295 keV 308 keV 316 keV 468 keV

0.2

604 keV

0 0

50

100

150

200

250

300

350 Angle (Phi)

Figure 5: The normalized TERMA point azimuthal distributions for the most intense 192 Ir primary decay photons.

variation in the fitting process iterations. Table 1: TERMA point radial fitting parameters for 316 keV photons. Name p0 p1

Water Parameters Value Error 2.36689 × 105 7.40208 × 101 1.15686 × 10−2 2.55892 × 10−6

Name p0 p1

Bone Parameters Value Error 4.33432 × 105 1.40132 × 102 2.07398 × 10−2 4.65519 × 10−6

Table 2: TERMA point radial fitting parameters in water for the other most intense 192 Ir primary decay photons. Name p0 p1 Name p0 p1

295 keV Parameters Value Error 7.67541 × 104 4.26752 × 101 1.18577 × 10−2 4.66323 × 10−6 468 keV Parameters Value Error 1.93001 × 105 6.19005 × 101 9.92395 × 10−3 2.25110 × 10−6

Name p0 p1 Name p0 p1

12

308 keV Parameters Value Error 8.30445 × 104 4.40445 × 101 1.16788 × 10−2 4.38028 × 10−6 604 keV Parameters Value Error 3.97388 × 104 2.66006 × 101 8.91551 × 10−3 4.21814 × 10−6

3.1.2. Angular Fitting After a comprehensive thorough investigation, a best match (Figure 4(b)) has been achieved using a fitting function of the form: p0 θ − p1 2 )) ∗ exp(−0.5 ( p2 2 ∗ π p2 θ − p4 2 p3 +√ ∗ exp(−0.5 ( )) p5 2 ∗ π p5 p6 θ − p7 2 +√ ∗ exp(−0.5 ( )) p8 2 ∗ π p8 + p9 .

T (θ) = √

(5)

The parameters of TERMA point angular fitting for 192 Ir photons are listed in Table 3 for 1 × 108 primary decays. 3.1.3. Two-Dimensional Fitting The behavior of the TERMA point azimuthal distribution is considered constant for ease of calculation, then the TERMA point fitting as a function of r and θ is constructed as in Equation 6. T (r, θ) = exp(−p0 r) ∗ p1 ∗ θ2  θ − p3 2 p2 )) ∗ exp(−0.5 ( × √ p4 2 ∗ π p4 p5 θ − p6 2 +√ ∗ exp(−0.5 ( )) p7 2 ∗ π p7  θ − p10 2 p8 ∗ exp(−0.5 ( +√ )) . p9 2 ∗ π p9

(6)

The fitting for 316 keV photons is depicted in Figure 6. The fitting parameters for all selected 192 Ir photons are also presented in Table 4. After thorough investigation, Eq. 6 appears to be the most appropriate distribution function for the generation of TERMA point distribution. 3.2. Kernel Results The Kernel defines the spread of energy as a function of radial distance, polar and azimuthal angle relative to the TERMA point “vector” direction. Keeping in mind the spikes described earlier in this paper which were excluded from the Kernel fitting function, the basic radial behavior of the fitted 13

Table 3: TERMA point angular fitting parameters for 316 keV photons and for the other most intense 192 Ir

primary decay photons (1 × 108

Name p0 p1 p2 p3 p4 p5 p6 p7 p8 p9 Name p0 p1 p2 p3 p4 p5 p6 p7 p8 p9

192 Ir

primary decays). 316 keV Parameters Name Value Error p0 -1.50387 × 108 2.95761 × 106 p1 6.42523 × 101 4.03967 × 10−1 p2 -4.93884 × 101 5.28137 × 10−1 p3 -6.86650 × 107 13.37649 × 105 p4 1.93004 × 102 1.38441 p5 -9.20639 × 101 7.7849 × 10−1 p6 -5.84234 × 107 2.47070 × 106 2 p7 1.39457 × 10 2.46660 × 10−1 p8 -3.54727 × 101 3.51303 × 10−1 5 p9 -6.25451 × 10 1.06714 × 104 295 keV Parameters 308 keV Parameters Value Error Name Value Error -3.68112 × 107 5.31781 × 105 p0 -4.50622 × 107 7.10946 × 105 5.88578 × 101 3.48511 × 10−1 p1 6.11931 × 101 3.75180 × 10−1 -4.33557 × 101 2.79089 × 10−1 p2 -4.56145 × 101 3.33298e × 10−1 9.12385 × 107 1.23225 × 105 p3 9.74313 × 107 1.22161 × 105 5.76922 × 102 7.94326 p4 5.13161 × 102 4.92006 3.57834 × 102 8.71111 × 10−2 p5 2.61902 × 102 6.75071 -2.83181 × 107 3.65045 × 105 p6 -2.80762 × 107 4.70725 × 105 1.35326 × 102 4.03549 × 10−1 p7 1.37414 × 102 3.84377 × 10−1 -3.92322 × 101 1.94415 × 10−1 p8 -3.81653 × 101 2.19278 × 10−1 -1.84547 × 105 1.59055 × 103 p9 -2.06361 × 105 2.28347 × 103 468 keV Parameters 604 keV Parameters Value Error Name Value Error -1.40319 × 108 7.79955 × 105 p0 -4.47350 × 107 5.67422 × 105 6.33794 × 101 1.61114 × 10−1 p1 4.31461 × 101 9.99812 × 10−2 -4.96189 × 101 1.63996 × 10−1 p2 -6.23659 × 101 3.82677 × 10−1 -6.01300 × 107 1.45883 × 105 p3 1.07243 × 106 6.42241 × 105 2 2 2.12416 × 10 2.43656 p4 4.68795 × 10 7.24457 × 10−1 -1.15498 × 102 3.78207 × 10−1 p5 8.09276 × 102 7.12038 × 10−1 7 5 7 -7.03110 × 10 7.10057 × 10 p6 -5.44795 × 10 4.44981 × 105 1.39124 × 102 1.55014 × 10−1 p7 1.37946 × 102 2.52264 × 10−1 -3.95755 × 101 7.855401 × 10−2 p8 -6.51859 × 101 3.37769 × 10−1 -6.05263 × 105 5.41365 × 103 p9 -2.80305 × 105 2.64848 × 103

14

Table 4: TERMA point two-dimensional fitting parameters for 316 keV photons and the other most intense

192 Ir

Name p0 p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 Name p0 p1 p2 p3 p4 p5 p6 p7 p8 p9 p10

primary decays = 1 × 108 . 316 keV Parameters Name Value Error p0 1.15840 × 10−2 2.45895 × 10−6 p1 1.43170 × 10−4 1.40142 × 10−6 6 p2 2.78197 × 10 3.46904 × 104 p3 1.58511 × 101 1.45155 p4 8.59644 × 101 6.66771 × 10−1 p5 -3.39882 × 107 4.03155 × 105 p6 1.75577 × 103 2.41762 × 101 p7 6.97789 × 103 8.79270e × 101 p8 4.38572 × 107 4.16344 × 105 p9 -9.97752 × 101 8.30422 × 10−1 p10 5.88334 × 101 1.29449 × 10−1 295 keV Parameters 308 keV Parameters Value Error Name Value Error 1.18890 × 10−2 4.74771 × 10−6 p0 1.17068 × 10−2 4.48742 × 10−6 1.02872 × 10−3 4.63350 × 10−6 p1 1.09865 × 10−3 6.49775 × 10−6 -1.78833 × 108 7.08733 × 105 p2 7.09450 × 108 4.22949 × 106 -1.01632 × 105 5.11153 × 102 p3 -6.26870 × 102 1.24777 3.85264 × 104 1.77056 × 102 p4 1.74907 × 102 2.19588 × 10−1 -4.16925 × 108 2.02485 × 106 p5 -4.31841 × 108 2.69151 × 106 -2.48980 × 106 1.14933 × 104 p6 -3.94555 × 104 2.30350 × 102 1.24546 × 105 5.14438 × 102 p7 4.48719 × 104 2.83138 × 102 1.50003 × 108 6.68071 × 105 p8 7.77970 × 108 4.78488 × 105 -5.69554 × 102 1.08535 p9 -5.87067 × 103 3.57845 × 101 2 −1 3 1.82396 × 10 2.14204 × 10 p10 4.13486 × 10 2.49244 × 101 468 keV Parameters 604 keV Parameters Value Error Name Value Error 9.94786 × 10−3 2.25531 × 10−6 p0 8.99125 × 10−3 4.04425 × 10−6 4.76861 × 10−5 3.33428 × 10−8 p1 3.24146 × 10−3 2.44690 × 10−5 4.47741 × 105 1.57779 × 103 p2 -6.96950 × 107 1.41423 1.15997 × 102 7.37775 × 10−2 p3 9.99383 × 10−2 1.41419 2.90851 × 101 3.74090 × 10−2 p4 4.09397 × 106 1.42751 2.58685 × 107 2.91406 × 104 p5 -1.61345 × 105 1.04718 × 101 -4.58600 × 101 7.57674 × 10−2 p6 4.25559 × 103 6.50374 × 101 6.38381 × 101 5.59302 × 10−2 p7 2.39663 × 104 3.76049 × 102 1.12916 × 107 1.30334 × 104 p8 1.25822 × 106 1.94708 × 101 -2.58569 × 101 2.12036 × 10−1 p9 -3.01250 × 102 8.88619 × 10−1 2.12724 × 101 1.15399 × 10−1 p10 1.42033 × 102 2.17120 × 10−1

photons. Number of

192 Ir

15

TERMA points x 0.316 MeV

104 3

10

102 10 1 0

200

400 Rad 600 ius ( mm800 ) 1000

1200

0

20

40

60

100 80 Theta

120

140

160

180

Figure 6: 316 keV TERMA point two-dimensional fitting for water phantom. The fitting curve (for 1 × 108 192 Ir primary decays) is displayed in red.

Kernel (presented in Figure 7(a)) spreads out and peaks around r = 100 mm from the TERMA point. Then, the Kernel declines to zero at deeper depths away from the TERMA point. Obviously, for 604 keV photons, the higher electron energies allow for the Kernel to penetrate to deeper distances, on the expense of the number of Kernel electrons at shorter distances. Figure 8(a) shows that 295 keV, 308 keV and 316 keV photons have a noticeable peak around 40◦ -45◦ relative to the TERMA point ray direction. For higher energies, namely 468 keV and 604 keV photons, the peaking is shifted toward smaller angles, 30◦ -40◦ , consistent with the fact that Compton scattered photons in this case tend to favour forward angles. This peaking at forward angles is obviously on the expense of Kernel dose contribution from these high energy photons in opposite direction to the TERMA point ray direction (backward angles). The energy dependence on the azimuthal direction (Figure 9) shows nearly constant behavior with a slight attenuation at angles 90◦ and 270◦ because of source capsule head effects at each end. This effect is spread out due to the fact that the azimuth angle is measured around the TERMA point ray instead of the phantom Y-axis. The 2D behavior in Figure 10 for 316 keV photons shows asymmetrical two dimensional surface. 16

Normalized Kernel

1 316 keV-Kernel 295 keV-Kernel 308 keV-Kernel 468 keV-Kernel 604 keV-Kernel

10− 1

10− 2

10− 3 0

200

400

600

800

1000 1200 Radius (mm)

Kernal (MeV)

(a)

50000

Fitting Curve Simulated Kernel

40000

30000

20000

10000

0

0

200

400

600

800

1000

1200 Radius (mm)

(b) Figure 7: (a) Normalized Kernel radial distribution for the most intense 192 Ir primary decay photons. (b) Fitted Kernel radial distribution for 316 keV photons obtained from 1 × 108 primary 192 Ir decays.

The Kernel is distributed intensively toward forward angles (20◦ − 60◦ ) at distances near the TERMA point (less than 100 mm). The Kernel surface is decreasingly distributed at higher angles and distances. 3.2.1. Radial Fitting Once the radial distribution is defined for the most intense five 192 Ir energies, their curves are fitted with a Poisson-like function as given in Equation 7. r

K(r) =

( pp01 ) p1 ∗ exp(− pp01 ) Γ( pr1 + 1)

r

+

( pp23 ) p3 ∗ exp(− pp23 ) Γ( pr3 + 1)

.

(7)

In this equation, Γ is the well-known gamma function. Then, the set of parameters are derived and listed in Table 5. An example for 316 keV photons 17

Normalized Kernel

1

316 keV Kernel 295 keV Kernel 308 keV Kernel 468 keV Kernel 604 keV Kernel

0.8

0.6

0.4

0.2

0 0

20

40

60

80

100

120

140

160

180 Theta

Kernel (MeV)

(a) ×103

Fitting Curve

1200

Simulated Kernel 1000

800

600

400

200

0

0

20

40

60

80

100

120

140

160

180 Theta

(b) Figure 8: (a) Normalized Kernel polar distributions for the most intense 192 Ir primary decay photons. (b) Fitted Kernel polar distribution for 316 keV photons obtained from 1 × 108 primary 192 Ir decays.

is depicted in Fig. 7(b). 3.2.2. Angular Fitting The Landau distribution (Landau (1944) and Schorr (1974)) combined with other functions are used to fit the angular behavior of Kernel. An excellent fitting is produced using Eq. 8, and the case of 316 keV photons is considered as an example depicted in Fig. 8(b). K(θ) = p0 ∗ Landau(θ, p1 , p2 ) + p3 ∗ exp(−p4 /θ) ∗ sin(p4 /θ).

(8)

Extracted parameters for 316 keV photons as well as for the other energies 18

Normalized Kernel

1

0.8

316 keV Kernel 295 keV Kernel 308 keV Kernel

0.6

468 keV Kernel 604 keV Kernel

0.4

0.2

0 0

50

100

150

200

250

Figure 9: Normalized Kernel azimuthal distributions for the most intense decay photons.

300

192

Ir primary

Table 5: Kernel radial fitting parameters for the 316 keV photons and the other most intense primary decay photons.

316 Name p0 p1 p2 p3

keV Parameters Value Error 128.464 0.542722 107.295 0.372625 151.300 0.699386 182.676 0.804350

295 Name p0 p1 p2 p3 468 p0 p1 p2 p3

keV Parameters Value Error 152.002 1.34922 177.085 1.52396 128.004 1.04121 103.387 0.627818 keV Parameters 127.183 0.920719 133.605 0.543127 152.085 1.28098 213.037 1.57354

are presented in Table 6.

19

308 Name p0 p1 p2 p3 604 p0 p1 p2 p3

350 Phi

keV Parameters Value Error 151.967 1.258900 180.051 1.434680 127.951 0.969579 105.694 0.614991 keV Parameters 110.615 1.080470 135.512 1.136050 184.845 1.795870 196.763 2.139290

192 Ir

Table 6: Kernel angular fitting parameters for the 316 keV photons and the other most intense

192 Ir

photons.

Name p0 p1 p2 p3 p4 Name p0 p1 p2 p3 p4

295 keV Value 2.41638 × 5.57980 × 2.43013 × -3.23738 × 6.58241 × 468 keV Value 7.46362 × 4.19238 × 1.48421 × -1.83723 × 4.59258 ×

316 keV Parameters Name Value Error p0 6.78358 × 106 2.24418 × 103 p1 4.15749 × 101 1.06817 × 10−2 p2 1.76555 × 101 6.64088 × 10−3 p3 -3.43253 × 106 1.07526 × 104 2 p4 4.68572 × 10 2.05382 × 10−1 Parameters 308 keV Parameters Error Name Value Error 106 2.65076 × 103 p0 2.33307 × 106 1.29281 × 103 101 5.12430 × 10−2 p1 4.21741 × 101 1.85460 × 10−2 101 2.87768 × 10−2 p2 1.79036 × 101 1.15200 × 10−2 105 1.14058 × 103 p3 -1.20251 × 106 6.39196 × 103 101 2.16859 × 10−1 p4 4.67767 × 102 3.51634 × 10−1 Parameters 604 keV Parameters Error Name Value Error 106 2.48709 × 103 p0 1.78165 × 106 1.23792 × 103 101 8.3445 × 10−3 p1 3.38022 × 101 1.66668 × 10−2 101 5.09907 × 10−3 p2 1.42941 × 101 1.05376 × 10−2 106 9.65974 × 103 p3 -1.81677 × 105 2.11333 × 103 102 3.39325 × 10−1 p4 3.76169 × 102 1.26195

3.2.3. Two-Dimensional Fitting The objective of curve fitting is to theoretically describe experimental data with a model function and to find the parameters associated with this model. Our study case in this section is to find the Kernel equation as a function of r and θ for the most intense 192 Ir decay energies, namely the 295 keV, 308 keV, 316 keV, 468 keV and 604 keV photons. The Kernel fitting finds the coefficients (parameters), which makes the function match the data as closely as possible. After thorough investigation, the best possible function describing the Kernel data (Fig. 10) was found as in Equation 9, and the resulting parameters are listed in Table 7. r − p1 2 p0 ∗ exp(−0.5 ∗ ( ) ) ∗ θ ∗ Landau(θ, p3 , p4 ). K(r, θ) = √ p2 2 ∗ π p2

(9)

3.3. Dose Calculation Making use of the fitting procedures described above, the constructed system for dose calculation at a given location in the phantom is based on the Kernel contribution from all TERMA points in the phantom. This new 20

Table 7: Kernel two-dimensional fitting parameters for the 316 keV photons and the other most intense 192 Ir

photons.

Name p0 p1 p2 p3 p4 Name p0 p1 p2 p3 p4

316 keV Parameters Name Value Error p0 2.59070 × 105 3.66997 × 102 p1 -1.16740 × 101 3.49612 × 10−1 p2 2.29520 × 102 1.57661 × 10−1 p3 3.03090 × 101 1.43628 × 10−2 1 p4 1.15922 × 10 4.03359 × 10−3 295 keV Parameters 308 keV Parameters Value Error Name Value Error 7.76312 × 104 1.84506 × 102 p0 8.73913 × 104 2.00526 × 102 -3.56152 5.75398 × 10−1 p1 -6.40488 5.57842 × 10−1 2.23933 × 102 2.65693 × 10−1 p2 2.26187 × 102 2.52912 × 10−1 3.08885 × 101 2.61821 × 10−2 p3 3.05374 × 101 2.48722 × 10−2 1.19158 × 101 7.33585 × 10−3 p4 1.17209 × 101 6.90878 × 10−3 468 keV Parameters 604 keV Parameters Value Error Name Value Error 3.40406 × 105 6.28965 × 102 p0 9.27299 × 104 3.78385 × 102 -6.01190 × 101 5.06699 × 10−1 p1 -8.22493 × 101 1.16574 2.61656 × 102 2.05449 × 10−1 p2 2.77480 × 102 4.46773 × 10−1 2.71003 × 101 1.31147 × 10−2 p3 2.52227 × 101 2.54074 × 10−2 9.75533 3.52200 × 10−3 p4 8.65717 6.58610 × 10−3

approach is tested at different Y and Z coordinates in planes located perpendicular to two points on the X-axis. A comparison is then carried out with the simulated dose scored by Ababneh et al. (2014). The simulated dose was calculated in 1 mm thick water slices positioned and centered at X = 5 mm and X = 20 mm, Figures 11 and 12, respectively. Generally, the dose calculated using the methodology described in this paper (Table 8) is in fair agreement with the directly simulated ”total“ dose (Figures 11 and 12). The fitting procedures is one main factor affecting the calculated dose. In addition, the dose calculation is foreseen to be further improved by tuning the dose contribution from the spikes excluded from the Kernel fitting, which are then added separately. This contribution is the focus of our next paper. Finally, it is worth mentioning here that the physics and calculations described in this paper do not take into consideration limited organ size and heterogeneities. Nevertheless, we provide the necessary elements to be utilized in any further detailed model based dose calculation.

21

3

Kernel (MeV)

10

102 10

1 0

200

400 Rad 600 ius ( mm)800

1000

1200

0

20

40

60

80 100 Theta

120

140

160

180

Figure 10: 316 keV Kernel two-dimensional fitting for water phantom. The fitting curve (for 1 × 108 192 Ir primary decays) is displayed in red.

Table 8: Dose at different positions in the phantom. The calculated dose is compared with the directly simulated ”total” dose as depicted in Figures 11 and 12. X (mm) Y (mm) Z (mm) Calculated Dose (Gy) 5 0.5 0.5 2.5 × 10−4 5 6.5 6.5 1.1 × 10−4 5 10.5 10.5 5.5 × 10−5 20 0.5 0.5 1.4 × 10−5 20 20.5 20.5 1.1 × 10−5 20 50.5 50.5 2.5 × 10−6

22

Simulated Dose (Gy) 2.4 × 10−4 8.0 × 10−5 4.0 × 10−5 1.6 × 10−5 6.0 × 10−6 2.0 × 10−6

Percentage Error 4% 27% 27% 14% 45% 20%

Total

-3

× 10 0.25

30 0.2

20 10

0.15

0 0.1

-10 -20

0.05 -30 -40 -40

-30

-20

-10

0

10

20

30

Primary

× 10

-3

0.22 20

0.2 0.18

10

0.16 0.14

0

0.12 0.1

-10

0.08 0.06

-20

0.04 0.02

-30 -30

-20

-10

0

10

20

30

Scatter

× 10

150

-6

16 100

14 12

50

10 0

8 6

-50

4 -100 2 -150

-100

-50

0

50

100

0

Figure 11: Y − Z dose spatial distribution (the dose scale is in Gy) for a 1 mm water slice around X = 5 mm. The number of primary decays was 105 × 106 (Ababneh et al., 2014).

23

Total 150

×10 18

100

16

-6

14 50

12 10

0

8 -50

6 4

-100

2 -150

-100

-50

0

50

0

100

Primary

× 10

-6

14

100

12 50 10 8

0

6 -50 4 -100

2

-100

-50

0

50

0

100

Scatter

× 10

150

-6

4 100

3.5 3

50

2.5 0 2 1.5

-50

1 -100 0.5 -150 -150

-100

-50

0

50

100

150

0

Figure 12: Y − Z dose spatial distribution (the dose scale is in Gy) for a 1 mm water slice around X = 20 mm. The number of primary decays was 105 × 106 (Ababneh et al., 2014).

24

4. Conclusions In this work, we addressed main basic physics components of developing model based dose calculation algorithm from MC data for HDR 192 Ir Brachytherapy. The investigated approach is developed based on GEANT4 MC as well as a number of dedicated ROOT and C++ analysis codes developed in this work. An analysis tool has been built and intensively used to calculate the TERMA point distribution and the related Kernel. The spatial, angular, and azimuthal one- and two-dimensional functions describing the distributions, as well as the corresponding parameters have been obtained. Three-dimensional integration of Kernel contribution, at a specified volume and location, from each TERMA point has been performed. Then, 3D integration of TERMA point distribution at that location has been also conducted to account for electrons generated at close distances from the TERMA point, which were originally excluded from the Kernel fitting. Ultimately, the described approach generates the dose at that specific location, which was then benchmarked against directly measured MC total doses. This methodology is much less time, data storage, and CPU consuming compared to repeating the extensive MC simulations. The work is merely intended to provide a launching point and basic key physics elements for anticipated future quantitative clinical use.

25

Ababneh, E., 2009. Evaluation of Scatter Dose Contribution of 192 Ir in Brachytherapy by Monte Carlo Simulation (Master Dissertation). Al-Balqa Applied University, Al-Salt, Jordan. Ababneh, E., Dababneh, S., Qatarneh, Sh., Wadi-Ramahi, Sh., 2014. Enhancement and validation of Geant4 Brachytherapy application on clinical HDR 192 Ir source. Radiat. Phys. Chem. 103, 57-66. Ahnesj¨o, A., 1989. Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media. Med. Phys. 16, 577-592. Angelopoulos, A., Baras, P., Sakelliou, L., Karaiskos, P., Sandilos, P., 2000. Monte Carlo dosimetry of a new 192 Ir high dose rate brachytherapy source. Med. Phys. 27, 2521-2527. Beaulieu, L., Carlsson Tedgren, A., Carrier, J., Davis, S., Mourtada, F., Rivard, M., Thomson, R., Verhaegen, F., Wareing, T., 2012. Report of the task group 186 on model-based dose calculation methods in brachytherapy beyond the TG-43 formalism: current status and recommendations for the clinical implementation. Med. Phys. 39, 6208-6236. Brady, L.W., Heilmann, H.-P., Molls, M., 2006. New Technologies in Radiation Oncology. Schlegel, W., Bortgeld, T., Grosu, A.-L., (Editors). Springer. . Carlsson Tedgren, ˚ A.K., Ahnesj¨o, A., 2003. Accounting for high Z shields in brachytherapy using collapsed cone superposition for scatter dose calculation. Med. Phys. 30, 2206-2217. Carlsson Tedgren, ˚ A., Ahnesj¨o, A., 2008. Optimization of the Computational Efficiency of a 3D, Collapsed Cone Dose Calculation Algorithm for Brachytherapy. Med. Phys. 35, 1611-1618. Landau, L., 1944. On the Energy Loss of Fast Particles by Ionisation. J. Phys. USSR. 8, 201-205. Mayles, P., Nahum, A., Rosenwald, J.C., (Editors), 2007. Handbook of Radiotherapy Physics, Theory and Practice, 1st ed. Taylor and Francis Group: CRC Press.

26

Melhus, C.S., Rivard, M.J., 2006. Approaches to calculating AAPM TG43 brachytherapy dosimetry parameters for 137 Cs, 125 I, 192 Ir, 103 P d, and 169 Y b sources. Med. Phys. 33, 1729-1737. Mikell, J.K., Mourtada, F., 2010. Dosimetric impact of an 192 Ir brachytherapy source cable length modeled using a grid-based Boltzmann transport equation solver. Med. Phys. 37, 4733-4743. Poon, E., 2009. Patient-specific dose calculation methods for high-dose-rate iridium-192 brachytherapy (Ph.D. Dissertation). McGill University, Montreal, Canada. Rivard, M.J., Butler, W.M., DeWerd, L.A., Huq, M.S., Ibbott, G.S., Meigooni, A.S., Melhus, C.S., Mitch, M.G., Nath, R., Williamson, J.F., 2007. Supplement to the 2004 update of the AAPM Task Group No. 43 Report. Med. Phys. 34, 2187-2205. Rivard, M.J., Coursey, B.M., DeWerd, L.A., Hanson, W.F., Huq, M.S., Ibbott, G.S., Mitch, M.G., Nath, R., Williamson, G.F., 2004. Update of AAPM Task Group No. 43 Report: A revised AAPM protocol for brachytherapy dose calculations. Med. Phys. 31, 633-674. ROOT. . Schorr, B., 1974. Programs for the Landau and the Vavilov distributions and the corresponding random numbers. Comput. Phys. Commun. 7, 215-224. Taylor, R.E., Rogers, D.W., 2008. EGSnrc Monte Carlo calculated dosimetry parameters for 192 Ir and 169 Y b brachytherapy sources. Med. Phys. 35, 4933-4944. Van der Laarse, R., Granero, D., P´erez-Calatayud, J., Meigooni, A.S., Ballester, F., 2008. Dosimetric characterization of 192 Ir LDR elongated sources. Med. Phys. 35, 1154-1161.

27