Monte Carlo characterization of 169Yb as a high‐dose‐rate source for brachytherapy application by FLUKA code

Higher initial dose rate and simplifying HDR room treatment of 169Yb element among other brachytherapy sources has led to investigating its feasibility as high‐dose‐rate seed. In this work, Monte Carlo calculation was performed to obtain dosimetric parameters of 169Yb, Model M42 source at different radial distances according to AAPM TG‐43U1 and HEBD Report about HDR sources in both air vacuum and spherical homogeneous water phantom. The deposited energy resulted by FLUKA as Monte Carlo code using binning estimators around 169Yb source was converted into radial dose rate distribution in polar coordinates surrounding the brachytherapy source. The results indicate a dose rate constant of 1.14±0.04cGy.h−1.U−1 with approximate uncertainty of 0.04%, air kerma strength, 1.082±2.6E−06U.mCi−1 and anisotropy function ranging from 0.386 to 1.00 for radial distances of 0.5–10 cm and polar angles of 0°–180°. Overall, FLUKA dosimetric outputs were benchmarked with those published by Cazeca et al. via MCNP5 as one of validate dosimetry datasets related to 169Yb HDR source. As a result, it seems that FLUKA code can be applicable as a valuable tool to Monte Carlo evaluation of novel HDR brachytherapy sources. PACS number: 87.15.ak


I. IntrodUCtIon
Radionuclides of lanthanide group are gaining importance as emerging therapeutic agents in nuclear medicine and radiation therapy. (1,2) 169 Yb is a lanthanide element that can be produced by neutron activation of 168 Yb in a nuclear reactor. (3) 169 Yb with a half-life of 32.02 ± 0.009 days decays by electron capture to 169 Tm and emits a low-energy photon spectrum with average energy of 93 keV. (4,5) It has an equivalent dose distribution similar to 192 Ir in water, but a much smaller half-value layer of 0.2 mm in lead and average tenth value thickness of 1 mm. (4) Therefore, these valuable features result in improved radiation protection consisting of simplified HDR room treatment with low level of installation costs, customized shielding of dose-limiting anatomic structures, and fabrication of shielded applicator to modulate the desirable dose distribution to achieve the specific requirements during adaptive brachytherapy treatment planning. (4) Moreover, it has a slightly higher initial dose rate of 12.5 cGy.h -1 , thus self-attenuation of medium, particularly soft tissue, would be negligible in comparison with other high-dose-rate (HDR) sources, (6) and penetration power would help to candidate 169 Yb as a valuable HDR sources in various treatment techniques of tumors such as intravascular brachytherapy (IVBT), prostate, breast, tongue, and gynecological applications. (4,7,8) According to the American Association of Physicists in Medicine (AAPM) TG-43U1 recommendations, before using each new source in treatment planning calculations, the dosimetric characteristics of the source must be determined to provide reliable data. (9)(10)(11)(12) The dosimetry dataset of 169 Yb as a novel HDR source has not been implemented to new updated reports and it has not been commercially fabricated yet. Therefore various models of 169 Yb have been investigated by experts, using different dosimetry characterization methods, ranging from Monte Carlo calculations to experimental measurements, with various level of perception of their results. (11) For instance, Cazeca et al. (5) evaluated dosimetry parameter based on Monte Carlo method like MCNP5, and Chiu-Tsao et al. (3) in 2000 investigated 169 Yb application on IVBT and presented results on the intravascular brachytherapy. In 2013, a paper was published by Saidi et al. (13) concerning experimental measurement and Monte Carlo methods of a type of eye plaques loaded with IRA1-103 Pd seeds for the radiation therapy application of choroidal melanoma. In some studies, authors compared different sources with each other to measure dose distribution and select appropriate source to special application radiotherapy. Mainegra and Capotez (14) in 1997 compare dose rate constants for 125 I, 103 Pd, 192 Ir, and 169 Yb brachytherapy sources via EGS4 as a Monte Carlo code. In 2005, a dosimetric comparison of 169 Yb versus 192 Ir for HDR prostate brachytherapy was done by Lymperopoulou et al. (15) and according to their results, the shielding requirements for the 169 Yb energies are minimal relative to those for 192 Ir. Furthermore, in some studies various techniques have been reviewed in order to reach dosimetric fundamentals for routine clinical use of HDR photon emitting brachytherapy sources, such as a publication by Li et al. (16) in 2007. In the current study 169 Yb source, Model M42, has been simulated via FLUKA as a Monte Carlo method and TG-43U1 brachytherapy dosimetry parameters such as anisotropy function, radial dose function, air kerma strength, and dose-rate constant have been obtained. Ultimately, to validate the data resulted by FLUKA, they have been compared with the results of similar source model as was used by Cazeca et al. based on MCNP5 scores.

A. 169 Yb source description
The source model employed in this research is similar to that applied in the investigation performed by Cazeca et al. (5) for the Monte Carlo analysis. The 169 Yb HDR source, Model M42, consists of a cylinder loaded by Ytterbium oxide of 7.1 mg.mm -3 and encapsulated by two long cylinders. The inner cylinder is made of Titanium of 4.51 mg.mm -3 and outer one is made of 304 stainless steel of 7.80 mg.mm -3 composition of 16% Cr, 10% Ni, 0.03% S, 0.75% Si, 0.08% C, 2% Mn, 0.045% P, 69.095% Fe, and 2% Mo. This source placed at the center of a 50 cm radius spherical phantom filled with homogenous water. The schematic view of 169 Yb, Model M42, is presented in Fig. 1 along with the coordinate system used. Due to the addition of an attached cable, the source assembly is asymmetric about the transverse plane principally, but symmetric about the coaxial axis.

B. Monte Carlo calculation method
In order to simulate photon transport of hypothetical 169 Yb HDR source situated once in center of a water phantom and in the air vacuum another time. To extract proportional scores, version 2011.2.14 (updated in 17th October 2012) of FLUKA code was utilized in the present research. This code is a multipurpose Monte Carlo code which has been developed for high accurate simulation the interaction and propagation in matter of about 60 different particles, including photons and electrons from 1 keV to thousands of TeV, neutrinos, muons of any energy, hadrons of energies up to 20 TeV and all the corresponding antiparticles, neutrons down to thermal energies and heavy ions. Moreover, time evolution and tracking of emitted radiation from unstable residual nuclei can be performed online by this code. Decay scoring, one of the valuable capabilities of this code, has been carried out in this study. (17)(18)(19)(20)

B.1 Essential setting
The lanthanide source of 169 Yb emits a photon spectrum that could be defined in the source routine and then it should be attached to FLUKA input. In the recent releases, new particle name has been implemented to beam card called ISOTOPE. By choosing this option and determining the type of isotope via the HI-PROBE card, the code applies relevant spectrum to input file, conveniently.
Since, in the current research, radioisotope decay is considered as a radiation source; therefore, appropriate cards of decay such as RADDECAY should be applied in the transport section of Flair as a FLUKA GUI. Furthermore, the card of DCYSCORE should be embedded in the scoring setup and then linked to the dose estimator of USRBIN.
In addition, the type of interaction, as well as the adequate thresholds, can be set in the physic and transport sections and to enable pure electromagnetic intersection input file was carried out in electromagnetic FLUKA (EMF) cascade mode. (17) Finally, this code was set at 3 × 10 8 of initial photon histories to minimize relative statistical error of simulation which varied from 3 × 10 -4 % to 8.31 × 10 -1 % with decreasing distance from source.

B.2 Geometry definition
According to common layout of detectors in the literature, (5) the group of detectors is usually composed of ring-shaped array formed by a chain of concentric spherical shells intersected by a chain of concentric cones. Figure 2 shows the detector arrangement defined by FLUKA geometry, in which 0.5 cm thickness of shells in 1 cm increments from 1 cm to 10 cm and angular aperture of cones in steps of 3° from 2° to 89°.
In spite of MCNP code, the number of bodies and regions in FLUKA does not have low limitation for geometry definition, which of course was not utilized in current work because of increasing run-time and decreasing simulating accuracy. (17) Thus, use of region binning method would be suggested. In fact, this novel approach was performed by definition of useful virtual cylindrical binning on the geometry in USRBIN cards in structure of the mesh R-Phi-Z (Fig. 3) as a deposited energy estimator in FLUKA code in unit of GeV per source particle and air-kerma strength were calculated in separate runs, in different structure.
Ultimately, dataset delivered by FLUKA estimators has been selected by replacement of the calculated position of detectors in phantom of water. To obtain the dose in various detectors, output results were divided by volume of each virtual bin, and then selected data were converted to cGy.h -1 .mCi -1 .
It should be noted that various data between two or four common virtual surfaces should be averaged. The input file has been run by a personal computer with Dual-Core CPU and Linux as an operating system; as a result, radiation dose in given detectors would be accessible with average error of well under 2.5 × 10 -2 %.

B.3 Obtaining dosimetry dataset
In the most relevant literatures, updated AAPM Task group 43 (TG-43U1) recommendations have been referred to evaluate high-dose-rate sources. In August 2012, AAPM and ESTRO's working group of High Energy Brachytherapy source Dosimetry (HEBD) issued a separate full report focused on dose calculation for photon-emitting brachytherapy sources with average energy higher than 50 keV. (21) As claimed by this report, the TG-43U1 guidelines and recommendations are applicable to HDR sources, unless configurations of sources are specified in new report based on some characteristics such as active length of source.
According to both standard methodologies with respect to cylindrical sources, the proposed formula for two-dimensional dose rate is: (11,21) (1) In this equation, dose rate is accessible via a given FLUKA estimator to evaluate absorbed dose distribution which utilized in other parameters such as dose rate constant, Λ, was calculated by dividing the dose rate, D • (r, ) θ , at 1 cm in water by the calculated air-kerma rate, K • (d) at 1 m in air. (11,21) Due to photon collisions, the produced electrons are low energy, and they are absorbed locally without any scattering or absorption in dosimetry measurement point, (22) and it could be assumed that dose is equal to kerma at the demand point. The air-kerma rate has been calculated by separate input file in which its geometry configured in a 1 m radius vacuum phantom surrounded by an air-filled ring detector made by an intersection of two spherical shells with inner and outer radius of 97.5 and 102.5, respectively, and two cones with 88° and 92° angles of aperture. (5) Consequently, by 1 × 108 histories, air-kerma rate was estimated 1.981 × 10 -4 cGy.h -1 .mCi -1 with statistical uncertainty of about 2.63 × 10 -6 %. Air-kerma strength (S k ) would be obtained by replacement of air-kerma rate in various positioned detectors based on distance, (Eq. (2)): (11,21) Obviously, air kerma rate follows inverse square law, owing to attenuation and scattering of photons in air. Moreover other dosimetric characteristics, including geometry function, G L , radial dose function (g L (r)), and anisotropy function (F(r,Θ)), were calculated in similar approach to that described in TG-43U1 protocol and HEBD report.

III. rESULtS & dISCUSSIon
The procedure applied to obtain dosimetry dataset for 169 Yb source, Model M42, in current issue has been evaluated by simulation of similar source characterized by Cazeca et al. (5) with MCNP5 and could be applied as an appropriate benchmark of obtained results.
As mentioned in previous section, scoring data from FLUKA estimators in water phantom were converted to dose rate, D • (r, ) θ , at the radial distance, r, and the polar angle, Θ, and have been presented at Table 1. The results of dose rate of polar angle 90° in various distances up to 10 cm are illustrated in Fig. 4, in comparison to reference data. Clearly, it could be seen that a complete conformity between two datasets of current research and the Cazeca study. The radial dose function, g L (r), was achieved by using dose rate parameter based on standard formalism (11,21) for 169 Yb HDR source with active length, L, of 0.6 cm.
The difference between the two studies is tabulated in Table 2. This difference could be related to variety of cross-sectional data utilized by FLUKA and MCNP5. Furthermore, by using Eq. (3), this radial dose function was fitted to a fifth order polynomial for commercial treatment planning purpose (Fig. 5): g L (r) = a 0 + a 1 r + a 2 r 2 + a 3 r 3 + a 4 r 4 + a 5 r 5 Where a 0 = 9.19 × 10 -1 , a 1 = 9.52 × 10 -2 , a 2 = -1.39 × 10 -2 , a 3 = 2.77 × 10 -3 , a 4 = -4.12 × 10 -4 , and a 5 = 1.96 × 10-5 , and define R² = 0.999. The 2D anisotropy function, F(r,Θ), of 169 Yb source, as another dosimetry parameter can be calculated by using of dose rate by means of FLUKA simulation, is reported in Table 3 in various radial distance, r, from 0.5 to 10 cm relative to source center, and polar angle, Θ, from 0° to 180° with respect to the source long axis. In fact, this factor indicates various doses in the longitudinal plane of source, because of different distribution of radioactivity along with cylindrical source, self-absorption, and self-filtration of seed layers. (22) Besides, the elliptical isodose curves in Fig. 6 demonstrate dose distribution around source. The significant dose attenuation can be seen in both end caps which applied in treatment planning purposes.  The aforementioned air-kerma strength, S k , was calculated as 1.082 ± 2.6E-06 U. mCi -1 (1 U = 1 μGy m 2 h -1 ), and is in good agreement with the study by Cazeca and colleagues of 1.08 ± 0.03 U.mCi -1 . The dose rate constant value, Λ, as accuracy criteria of present Monte Carlo simulation in comparison with reference work was obtained by 1.14 ± 0.04 cGy.h -1 .U -1 , which is close to 1.12 ± 0.04 cGy.h -1 .U -1 found by Cazeca et al. (5)

IV. ConCLUSIonS
The dosimetry of 169 Yb HDR brachytherapy source, Model M42, was investigated in this study in association with FLUKA estimators, according to AAPM standard methodology. In summary, simulation results of FLUKA code are consistent with those published by Cazeca et al. (5) via MCNP5 as one of validate dosimetry datasets related to 169 Yb HDR source. Actually, it seems that FLUKA code with high quality of GUI can be facilitated by implementing standard protocols' procedures similar to MCNP. Therefore, it can be applicable as a valuable alternative tool to Monte Carlo evaluation of novel HDR brachytherapy sources.