Evaluation of a deterministic grid-based Boltzmann solver (GBBS) for voxel-level absorbed dose calculations in nuclear medicine

To evaluate the 3D Grid-based Boltzmann Solver (GBBS) code ATTILA® for coupled electron and photon transport in the nuclear medicine energy regime for electron (beta, Auger and internal conversion electrons) and photon (gamma, x-ray) sources. Codes rewritten based on ATTILA are used clinically for both high-energy photon teletherapy and 192Ir sealed source brachytherapy; little information exists for using the GBBS to calculate voxel-level absorbed doses in nuclear medicine. We compared DOSXYZnrc Monte Carlo (MC) with published voxel-S-values to establish MC as truth. GBBS was investigated for mono-energetic 1.0, 0.1, and 0.01 MeV electron and photon sources as well as 131I and 90Y radionuclides. We investigated convergence of GBBS by analyzing different meshes (M0,M1,M2), energy group structures (E0,E1,E2) for each radionuclide component, angular quadrature orders (S4,S8,S16), and scattering order expansions (P0–P6); higher indices imply finer discretization. We compared GBBS to MC in (1) voxel-S-value geometry for soft tissue, lung, and bone, and (2) a source at the interface between combinations of lung, soft tissue, and bone. Excluding Auger and conversion electrons, MC agreed within  ≈5% of published source voxel absorbed doses. For the finest discretization, most GBBS absorbed doses in the source voxel changed by less than 1% compared to the next finest discretization along each phase space variable indicating sufficient convergence. For the finest discretization, agreement with MC in the source voxel ranged from  −3% to  −20% with larger differences at lower energies (−3% for 1 MeV electron in lung to  −20% for 0.01 MeV photon in bone); similar agreement was found for the interface geometries. Differences between GBBS and MC in the source voxel for 90Y and 131I were  −6%. The GBBS ATTILA was benchmarked against MC in the nuclear medicine regime. GBBS can be a viable alternative to MC for voxel-level absorbed doses in nuclear medicine. However, reconciliation of the differences between GBBS and MC at lower energies requires further investigation of energy deposition cross-sections.


investigated convergence of GBBS by analyzing different meshes (M M M
, ,

Introduction
Voxel-based dosimetry models (VBDM) are the norm in radiation oncology practice where external beam and sealed source brachytherapy procedures are performed. Several nuclear medicine researchers have implemented their own VBDM, but unfortunately such VBDM are not widely used within nuclear medicine departments (Sgouros et al 1990, Furhang et al 1996, Dieudonne et al 2010, Tagesson et al 1996, Descalle et al 2003, Dewaraja et al 2005, Chiavassa et al 2006. Nuclear medicine instead continues to rely on anthropomorphic phantoms that are not patient-specific, assume uniform activity distributions throughout organs, and do not include tumor models (Sgouros and Hobbs 2014). Some researchers have added simple tumor models these phantoms (Clairand et al 1999, Johnson et al 1999, Stabin et al 2005. Studies have shown that VBDMs using patientspecific anatomy including tumors improves treatment planning and refines predictions of response and toxicities (Dewaraja et al 2010, Strigari et al 2010, Senthamizhchelvan et al 2012, Dewaraja et al 2014. Implementation of VBDMs in nuclear medicine range from simplified complete local deposition to solutions of the fully coupled transport equations for electrons, positrons, and photons (Furhang et al 1996, Tagesson et al 1996, Ljungberg et al 2002, Pasciak and Erwin 2009, Pacilio et al 2015, Pasciak et al 2014, Mikell et al 2015. VBDMs more sophisticated than complete local deposition include dose voxel kernels (Bolch et al 1999) and collapsedcone convolution (Sanchez-Garcia et al 2014). Solutions to the linear Boltzmann transport equation (LBTE) have traditionally been solved stochastically through Monte Carlo (MC). However, deterministic methods also exist to solve the LBTE (Lewis 1993).
Deterministic solvers of the LBTE are relatively new to the radiotherapy community. Only recently have grid-based Boltzmann solvers (GBBS) been adapted for clinical use with external photon beams and 192 Ir sealed source brachytherapy (Gifford et al 2006, Vassiliev et al 2008, Gifford et al 2010, Mikell and Mourtada 2010, Vassiliev et al 2010, Han et al 2011, Han et al 2012. They have even recently been investigated for radiation transport in the presence of magnetic fields (Aubin et al 2015).
ATTILA ® (Varian Medical Systems, Palo Alto, CA) is a general purpose GBBS code, and it has been investigated for applications to megavoltage photon beams and high energy sealed source gamma emitters ( 137 Cs,192 Ir) (Gifford et al 2006, Vassiliev et al 2008. It has not been studied for nuclear medicine applications, where there is a primary electron (betas, Auger electrons, internal conversion electrons) and primary photon (x-rays, gammas) source. In the nuclear medicine therapy regime, the primary electron sources are ultimately responsible for the majority of energy deposition; the x-ray and gamma sources transfer energy to electrons that transport and deposit energy, but this energy deposition is typically orders of magnitude less than that caused by the primary electron sources. This regime is in contrast to previous studies with GBBS that have focused on the energy deposition of electrons liberated by high-energy megavoltage photon beams or scored KERMA for sealed-source brachytherapy. When compared to external beam and sealed source brachytherapy, unsealed sources produce absorbed doses with much larger gradients and increased heterogeneity. It is important to evaluate the GBBS at the lower (relative to external beam) energies found in nuclear medicine to ensure that its transport and cross sections are sufficient for unsealed source voxel-level dosimetry.
The goal of this investigation was to evaluate the GBBS ATTILA for voxel-level dosimetry in nuclear medicine. To accomplish this, MC voxel-S-values were calculated with DOSXYZnrc and benchmarked against published voxel-S-values. Then, taking MC as the truth, absorbed doses calculated using the GBBS were compared to those estimated from MC for monoenergetic photons and electrons, 90 Y, and 131 I in the following configurations: (1) voxel-S-values in uniform soft tissue, bone, and lung; and (2) a single source voxel at the interface of materials. This work represents the first step in showing the GBBS is suitable for absorbed dose calculations in nuclear medicine.

DOSXYZnrc Monte Carlo simulations
The MC method was used as the gold standard for comparisons with the GBBS code. MC simulations were performed using DOSXYZnrc (Walters and Rogers 2002) which is a user code of EGSnrc (version 4.2.4.0) (Kawrakow and Rogers 2000). EGSnrc is a general purpose radiation transport code with improved low energy support compared with EGS4 (Kawrakow 2000). DOSXYZnrc allows scoring in the familiar Cartesian voxel geometry ubiquitous in medical imaging.
The simulation geometry employed was the 3 mm voxel-S-value geometry found in MIRD Report 17 (Bolch et al 1999). Briefly, this phantom geometry consisted of a 3D infinite distribution of soft tissue segmented into 3 mm isotropic voxels, where the center voxel was uniformly filled with activity. The absorbed dose to each target voxel from the source voxel was then reported as a function of radial distance. The source voxel center to target voxel center was defined as a radial distance 2 , where ∆ = 3 mm and ′ ′ ′ i j k , , represent 0-based indices relative to the source voxel.
Radionuclide spectra were obtained from RADTABS software (Eckerman and Endo 2008). For each radionuclide, independent simulations were performed for the following emission components: (1) photons: gamma and x-ray; (2) continuous beta spectra ( ) β − ; and (3) mono-energetic electrons: Auger and internal conversion (Auger + IC). The simulated Auger component was the collapsed Auger electrons listed in the *.RAD file output by RADTABS. The continuous beta spectrum is often the dominant source of energy deposition in the source voxel, but Auger + IC may provide substantial contributions. Beyond the range of electron sources, the gamma and x-ray emissions are dominant. We did not simulate the gamma and x-ray emissions or Auger and conversion electrons for 90 Y because their contributions are negligible for dosimetry.
As an initial check, we compared our MC simulations with published data for voxel-S-values (Bolch et al 1999, Pacilio et al 2009, Dieudonne et al 2010, Lanconelli et al 2012, Amato et al 2013. For soft tissue we compared our MC simulation results with published data for 131 I, 90 Y, 153 Sm, 177 Lu, and 99m Tc (Bolch et al 1999, Lanconelli et al 2012, Amato et al 2013. For bone we compared our MC simulation results for 131 I, 90 Y, 153 Sm, and 177 Lu with Lanconelli et al (2012). The published data were generated using DOSXYZnrc (Lanconelli et al 2012), Geant4 (Amato et al 2013), and EGS4 (Bolch et al 1999). We tabulated the emission contributions (gamma + x-ray, beta, Auger + IC) from our MC simulations as a percentage of the total voxel-S-value at each voxel for the radionuclides investigated. When comparing with Lanconelli et al (2012), we performed comparisons using two different MC voxel-S-values. The first was our standard calculation that included contributions from all three components including gamma + x-ray, beta, and Auger + IC; the second excluded contributions from the Auger + IC.
For quantitative comparisons we followed the analysis of Pacilio et al (2009) and investigated differences in the source voxel (0 0 0), nearest neighbor (0 0 1), and along the diagonal (0 1 1). For qualitative comparisons we plotted the voxel-S-values by collapsing them to one dimension as a function of the radius, r i . The voxel-S-values were sorted by radial distance (r i ) such that r r r n . In addition to the absolute voxel-S-values, we plotted percent differences of MC relative to published data as a function of r i .
For comparing GBBS with MC, we investigated three different materials including soft tissue (Cristy and Eckerman 1987) (ρ = 1.04 g cc −1 ), cortical bone (Snyder et al 1975) (ρ = 1.85 g cc −1 ), and lung (Snyder et al 1975) (ρ = 0.26 g cc −1 ). The two most common radionuclides, 90 Y and 131 I, for internal radionuclide therapy were simulated. Mono-energetic 1 MeV, 100 keV, and 10 keV electrons and photons were simulated to span the energy range found in typical nuclear medicine procedures.
All electrons and photons were tracked down to kinetic energies of 1 keV. The simulation parameters were set to all the advanced options including bound Compton scattering, Rayleigh scattering, atomic relaxations, electron impact ionization, XCOM photon cross sections (Berger et al 2010), spin effects, exact boundary crossing, and PRESTA-II (Walters and Rogers 2002).
Our MC simulations were performed using 1E + 09 source particles for each emitted radiation: (1) gamma + x-ray, (2) beta, and (3) Auger + IC. The full voxel-S-values were constructed by adding the individual emitted radiation components together with weightings taken from RADTABS (Eckerman and Endo 2008). In addition, a bar plot was generated to visually compare the differences among the multiple MC codes in the source voxel.
as well as high energy gamma ( 137 Cs, 192 Ir) emitting sealed sources (Gifford et al 2006, Vassiliev et al 2008. ATTILA discretizes space, angle, and energy to solve the LBTE for photons and the linear Boltzmann-Fokker-Planck transport equation (LBFPTE) for charged particles. It solves for the energy and angular dependent particle flux r E , , ( ) Ψ Ω → throughout space. The GBBS ATTILA has the following properties: (1) support for electrons, positrons, photons, and neutrons; (2) fully-coupled transport where electrons and positrons generate photons, and photons generate electrons and positron; (3) unstructured tetrahedral mesh to represent geometry; (4) linear discontinuous finite element in space; (5) multi-group discretization in energy; (6) discrete ordinates differencing in angle; and (7) spherical harmonics expansion of scattering source (Gifford et al 2006, Vassiliev et al 2008, McGhee et al 2009.
The energy dependent scalar flux, r E , ( ) Φ → , is computed by ATTILA and used to calculate reaction rates corresponding to absorbed dose (energy reaction rate) for electron groups or kinetic energy released in material (KERMA) for photons. Energy dependent reaction rate cross sections were calculated by the cross section generator described later. Details of the LBTE, LBFPTE, and calculating reaction rates as post-processing steps are published elsewhere (Gifford et al 2006, Vassiliev et al 2008. ATTILA requires multi-group energy cross sections, which were generated using ZERKON version 1.0.0 (Varian Medical Systems, Palo Alto, CA). ZERKON is based on CEPXS (Lorence et al 1989). CEPXS has been benchmarked against several other codes suggesting it sufficiently accounts for electron and photon interactions (Lorence et al 1990). For cross section generation, full-coupling was used for both photon and electron sources with a Legendre expansion of order 7. Cross sections were generated for soft tissue, lung, and bone matching the atomic composition and densities used in the MC simulations. Cross section files were generated for each radionuclide emission component and each energy group discretization.

Voxel-S-value simulations
For the initial evaluation of ATTILA in the nuclear medicine regime, we chose to calculate voxel-S-values, as published by the Society of Nuclear Medicine and Molecular Imaging's Medical Internal Radionuclide Dosimetry (MIRD) committee (Bolch et al 1999). Voxel-Svalues are well known in the nuclear medicine community, and publications offer tabulated results for quantitative comparisons (Bolch et al 1999, Lanconelli et al 2012, Amato et al 2013. Most publications to date have reported voxel-S-values only in soft tissue. In this work, we performed calculations in soft tissue, bone, and lung. 2.3.1. ATTILA solver settings. We used square Chebychev Legendre quadrature sets, with Galerkin scattering treatment, and diagonal transport correction. For photon sources, we set the diffusion synthetic acceleration (DSA) to simplified_WLA. DSA was turned off for the electron sources. Further detail on DSA can be found in the literature (Wareing et al 2001). All electron and photon energy cuts were 1 keV.

2.3.2.
Sensitivity analysis for ATTILA. The GBBS requires sufficient discretization in space, angle, energy, and scatter source to converge to a solution. In this work, we wanted to ensure that the space, angle, energy, and scatter source (also known as phase space variables) converged sufficiently in and near the source; this was the rationale behind the sensitivity analysis. To evaluate the effect of the phase space parameters, we generated calculations for three tetrahedral meshes ( respectively for the Square-Chebychev quadrature set, and up to a 7th order spherical harmonics polynomial expansion of the scattering source (P 0 -P 6 ). Higher numbers on the subscripts indicate finer discretizations.
To quantify the convergence in terms of the four discretization parameters, we calculated the ratio in the source voxel absorbed dose along each variable of the finest discretization to the next finest for each variable: M M E P S ME P S  ) with ≈16 000 tetrahedrons, and a maximum tetrahedral edge length of 0.5 mm was selected to create the fine mesh (M 2 ) with ≈64 000 tetrahedrons (figure 1).

Energy groups.
We investigated three energy group discretizations for each radionuclide/source component. The coarse (E 0 ), intermediate (E 1 ), and fine (E 2 ) discretizations had approximately 30, 60, and 90 groups for each particle respectively, yielding approximately 60, 120, and 180 total groups for photons and electrons. Energy groups were distributed logarithmically using an approximately constant number of groups per decade down to 1 keV. The energy group widths are described in table 1. Photon energy group structure was adjusted for gamma and x-ray emission photo-peaks contributing >99% of the total gamma + x-ray emission energy for 131 I (80.2, 284, 326, 364, 503, 637, 643, and 723 keV). Energy group widths of 1 keV centered (e.g. 363.5-364.5 keV) on the photo-peaks were added to the photon group structure. Similar adjustments to the group structure were made for the 131 I Auger + IC emission component. Group structures for monoenergetic electron and photon sources were also similarly adjusted using 1 keV group widths for 1 MeV and 100 keV while a 0.1 keV group width was used for 10 keV sources.

Comparison metrics.
Similar to the validation of our MC, we performed quantitative comparisons of ATTILA (coarsest and finest discretizations) in and around the source (0 0 0, 0 0 1, 0 1 1) by calculating percent differences of the GBBS relative to our MC. The ground truth in the comparisons was our MC. Qualitative comparisons of percent differences were performed graphically for the finest discretization using radial plots of the voxel-S-values. All points were not plotted in the percent difference figures to maintain clarity. Only the average statistical uncertainty from MC and the maximum and minimum percent differences at each radial distance were plotted.

Single source voxel at interface simulations
To evaluate the GBBS with non-uniform materials, we performed simulations at an interface of two materials. This consisted of a single voxel of activity on one side of the interface. The interfaces studied were lung → soft tissue (L_S), soft tissue → lung (S_L), bone → soft tissue (B_S), and soft tissue → bone (S_B); the first material represents the material containing the source voxel. The simulation geometry is illustrated in figure 2. The scoring geometry matched MC. The same multi-group cross sections and source spectra were used as in the voxel-S-value simulations. However, a full voxel geometry (not octant) and tetrahedral mesh consisting of ≈120 000 tetrahedrons were generated (supplemental_figures.docx: supplemental figure 1 (stacks.iop.org/PMB/61/4564/mmedia)). The energy groups, scatter expansion, and angular quadrature set order were identical to the finest discretization in the uniform material simulations.
Similar to the procedure for uniform voxel-S-values, we performed quantitative compariso ns of ATTILA using tabulations near the source voxel. On the source side, we investigated the source voxel (0 0 0), its neighbor (0 0 1), and along the diagonal (0 1 1). We investigated voxels across the interface from the source (0 0 1*) and diagonally across the interface from the source (0 1 1*). Local percent differences were evaluated with the ground truth set to our MC. A notable difference in this analysis was that we only investigated voxels with a boundary on the interface. Evaluated voxels existed in the plane defined by z = 0 or z = +1. We used the r i in the two planes on either side of the interface as shown in figure 2. Qualitative comparisons were performed graphically using radial plots of voxel-S-values and percent differences for voxels in the z = 0 or z = +1 plane.

DOSXYZnrc MC comparison with published data
Our MC statistical uncertainty was ⩽0.05% in the source voxel (0 0 0) for all investigated sources in soft tissue. For all materials and all sources, the source voxel uncertainty  was ⩽0.09%; similarly, the neighboring voxel (0 0 1) was ⩽0.4%. Excluding 10 keV electrons, the statistical uncertainty in the diagonal voxel (0 1 1) was ⩽0.7%. Statistical uncertainty for 10 keV electrons in the diagonal voxel was 10%, 30%, and 4% for soft tissue, bone, and lung, respectively. Figure 3 shows 131 I plots comparing our DOSXYZnrc MC simulations with published data in soft tissue. Additional plots for 153 Sm, 177 Lu, 90 Y, and 99m Tc are available in the supplemental data (MC_compared_to_published_data\MC_validation_plots.docx). The absolute voxel-S-value graphs were in good qualitative agreement on the logarithmic scale, where the absorbed dose can change by 4 orders of magnitude from the source voxel to the bremsstrahlung tail. The local percent differences in figure 3 indicate excellent agreement with published data for 131 I. Specifically, the MC 131 I absorbed doses near the source were within 3% of Lanconelli et al (2012) and Amato et al (2013) while the differences at distance were within 5%. Compared to EGS4 used by Bolch et al (1999), the MC 131 I source voxel had good agreement (<3%); larger differences of 25% were observed for the 001 and 011 voxels, but these differences were also observed by both Lanconelli et al (2012) and Pacilio et al (2009). At distances further from the source, the differences were within ≈20%.
Absorbed dose differences in the source voxel between MC and published data were within 5.3% for 131 I, 90 Y, and 99m Tc (figure 4). However, MC differences with Lanconelli et al (2012) were ≈8% and ≈23% for 177 Lu and 153 Sm, respectively, whereas MC differences with Amato et al (2013) were only ≈3% and ≈5%, respectively. Excluding the Auger + IC component from our MC simulations improved agreement with Lanconelli et al (2012) to within 3%. Figure 2. Illustration of the interface simulation geometry with 3 mm voxels. The source material was defined for voxel centers with z ⩽ 0, whereas the other material was defined for voxel centers having z ⩾ +1. The source voxel is at the center of the z = 0 plane. We evaluated the voxel-S-values and percent differences for radii confined to the planes z = 0 and z = +1. The ordering of radii in each planar comparison is shown as ordered subscripts. Similar levels of agreement were found in bone. Supplemental data contains plots (MC_com-pared_to_published_data\MC_validation_plots.docx) and tabulated data (MC_compared_to_ published_data\MC_validation_*.txt files) for 90 Y, 131 I, 153 Sm, 177 Lu, and 99m Tc in soft tissue, bone, and lung.

ATTILA voxel-S-value simulations
3.2.1. Sensitivity results. Figure 5 illustrates convergence of the GBBS voxel-S-values for a 1 MeV electron source in soft tissue as the space (M M M , , 2 ), energy (E E E , , 0 1 2 ), angle (S S S , , 4 8 16 ), and scattering source expansion (P 0 -P 6 ) were refined. The effect of the angular discretization was apparent at larger radii with increased angular quadrature order leading to reduced variations in the absorbed dose at distance (figure 5(a)). Negative fluxes were occasionally returned by the coarsest discretization which resulted in negative voxel-S-values. Negative values are an indication that the phase-space variables have not been sufficiently refined. A few negative values still existed with the intermediate discretization, but the negative values disappeared for the finest discretization.
Effects of the mesh, energy and scatter source expansion discretizations for S 16 were observed as changes in the source voxel-S-value for 1 MeV electron (figures 5(c) and (d)). Increasing the energy group discretization did change the solution, but given P 3 or higher then the difference between E 1 and E 2 was ≈0.5%. The difference between E 0 and E 1 given P 3 or higher was ≈2%. The spatial discretization differences between M 1 and M 2 were less than 0.1%. However, differences between M 1 and M 0 ranged from 4.5% for P 0 to 1.5% for P 3 and higher (figures 5(c) and (d)).
For the finest discretizations, most GBBS absorbed doses in the source voxel changed by less than 1% compared with the next finest discretization along each phase space variable. This indicated that our finest discretization was sufficiently refined. In the source voxel for all sources and materials investigated M ∆ and S ∆ ranged from 0.999 to 1.005. E ∆ ranged from 0.964 to 1.001; most E ∆ were from 0.990 to 1.001 with worse values for photons in bone and for 100 keV photons in lung. P ∆ ranged from 1.000 to 1.005 with the exception of 1 MeV photon in lung which was 1.015.

Comparison with DOSXYZnrc MC.
Qualitatively, the finest discretization of GBBS and MC in figure 6 exhibited good agreement. The agreement worsened at larger radii (figure 5(a)), most likely due to ray effects and increased MC statistical uncertainty. However, it should be noted that ray effects can be reduced through the use of a first scattered distributed source calculation or further increasing the number of angles, , McGhee et al 2009) but this was not performed in the current study. Plots in bone and lung are available as supplemental data (supplemental_figures.docx: supplemental figure 2). Figure 7 illustrates quantitative differences between GBBS and our DOSXYZnrc MC by plotting the minimum and maximum local percent differences, assuming MC as the truth, at a given radius for the finest GBBS calculations. The overall trend is for the GBBS to produce absorbed doses less than MC. The magnitude of the differences was generally <20% within 2 cm of the source voxel. For electron sources, larger differences occurred near the end of the beta range for 90 Y and 1 MeV in soft tissue ( figure 7(a)). The 131 I Auger + IC and 10 keV sources showed the largest percent differences, but the MC statistical uncertainty was also largest for these components (figures 7(c)-(e)). Reporting GBBS KERMA instead of absorbed dose for 100 keV and 10 keV photons improved agreement with MC (figures 7(d) and (f)); the source voxel agreement changed from −9.9% to 0.3% for 100 keV and −15.3% to −1.4% for 10 keV. Figure 8 complements figure 7 by comparing both GBBS and published data relative to MC in soft tissue. The magnitude of differences in the source voxel are shown in figure 8(a) for both the finest and coarsest GBBS calculations; the finest GBBS was within 6% of MC. The magnitude of differences were, on average, slightly larger with GBBS, but the differences were comparable to the magnitude of differences seen between other Monte Carlo codes which were on the order of a few to 5% in the source voxel.
For 90 Y and 131 I, the differences between the coarsest and finest discretized GBBS calculations with MC were approximately 5% and 1%, respectively. The effect of discretization was larger at distance. For example, at 10.3 mm (not shown) for 131 I the coarsest and finest GBBS calculations' agreement with MC was −93% and −6.5%, respectively. Figure 8(b) illustrates that higher energy emissions in the source voxel require finer discretization to converge, and that differences between converged GBBS and MC become increasingly negative at lower energies. The converged GBBS consistently underestimated MC in the source voxel (figures 7 and 8(b)). The bias was approximately −7% for 1.0 MeV photons & electrons and worsened to approximately −15% for 10 keV photons & electrons across all materials studied.

ATTILA and MC interface simulations
The MC and finest GBBS absorbed doses to voxels in both planes along the interface were in good qualitative agreement. Figure 9 shows the soft tissue-lung (S_L) interface. The remaining Convergence as a function of scatter source expansion P N for the 1 MeV electron source voxel-S-value (square = M 2 , triangle = M 1 , circle = M 0 , red = E 2 , green =E 1 , orange = E 0 ) is shown in (c) with a zoomed view in (d). In (d) notice that moving from M 1 to M 2 makes little difference ( 1%), moving from 5 to 6 in P N makes little difference ( 1%) for E 1 and E 2 , and moving from E 1 to E 2 results in about a 0.5% difference; this tells us that our GBBS solution, at least in the source voxel, is converged. interface plots (L_S, B_S, S_B) are available as supplemental data (supplemental_figures. docx: supplemental figures 6 and 7). An artifact can be seen near the end of the 1 MeV electron range in figure 9. This dipping artifact was also seen with the 1 MeV electron source for the voxel-S-value in lung (supplemental_figures.docx: supplemental figure 2(e)). The artifact occurred because the solution (angular flux) was changing too quickly at that spatial location. The dip became wider, deeper, and included negative values for the next coarsest mesh with other phase-space variables held constant. Thus, we expect further mesh refinement to resolve the dipping artifact.
For the interface simulations, local percent differences between GBBS and MC were similar to the differences between GBBS and MC in the uniform material simulations. To conserve space, the interface percent difference plots are not in the manuscript. Instead, they are available as supplemental data (supplemental_figures.docx: supplemental figures 8 and 9).

Discussion
To benchmark our MC models, we first compared our DOSXYZnrc MC voxel-S-values with published data generated using DOSXYZnrc, Geant4, and EGS4. Pacilio et al (2009) characterized differences of a few percent (−3.5% to +4%) for mono-energetic electrons and differences up to 7% for mono-energetic photons in a 3 mm source voxel. In general, our comparisons agreed with the observations of Pacilio et al (2009), but we did note larger discrepancies between our MC estimates and Lanconelli et al (2012) for radionuclides with nonnegligible Auger + IC emissions ( 153 Sm, 177 Lu). This was surprising because we both used DOSXYZnrc. A possible explanation for the difference is that Lanconelli et al (2012)  For voxel-S-values, the range of absorbed doses can span several orders of magnitude over a few millimeters. Figure 7 shows large percent differences (>20%) further from the source,  but the magnitude of these differences is similar to that reported by other publications (Pacilio et al 2009, 2015, Amato et al 2013. Previous work in external beam and sealed source brachytherapy found minimal differences between ATTILA and other transport codes. Gifford et al (2006) scored photon KERMA in a plane around an ovoid containing a 137 Cs source and found differences between ATTILA and MCNPX to be on the order of 2-5%. The same study compared ATTILA with EGS4 for an 18 MV percent depth dose curve in a heterogeneous phantom and found the largest difference was ATTILA overestimating the EGS4 values by 2.2%. An additional study by Gifford et al (2008) compared ATTILA with MCNPX photon KERMA around an 192 Ir source and   (1 MeV), green lines and diamonds (100 keV), and purple lines and squares (10 keV) represent both mono-energetic electron and photon sources. 131 I components are represented by the color blue with X's (beta) and triangles (Auger + IC) for electrons and X's (x-ray + gamma) for photons. 90 Y is represented by orange pentagons and doesn't have a photon source. The legend is identical to figure 6. found 98% of voxels to be within 5% of MCNPX; although there were localized differences of −7% beyond the source tip and differences over +5% near the source shown graphically. Vassiliev et al (2008) compared ATTILA with EGSnrc for clinical patient cases treated with 6 MV photon beams and reported >98% of voxels had passing gamma indices (Low et al 1998) for 3%/3 mm.
Previous publications with the GBBS focused on external beam and high energy sealed source brachytherapy. In this work, we found that differences exist between MC and GBBS in the nuclear medicine energy regime. Gifford et al (2006) suggested that differences in cross sections used by the GBBS and MC over the energy range of interest in external beam and high energy brachytherapy were minimal. Our results suggest that the differences in cross sections become noticeable at lower electron energies including the nuclear medicine energy regime.
We hypothesize that these differences manifest in the energy reaction rate cross sections, which are used to convert the local electron spectra to absorbed dose. Comparing both KERMA and absorbed dose GBBS calculations with MC calculated absorbed doses for low energy photon sources (figures 7(d) and (f)), where KERMA is a good approximation to absorbed dose, showed that the KERMA reaction rate agreed within 3.3% of the MC calculated absorbed dose over all materials studied. Thus, for low energy photons, calculating KERMA instead of absorbed dose with the GBBS yielded better agreement with MC absorbed doses.
To determine whether the differences in absorbed dose were from transport or the energy reaction rate cross sections, we ran ATTILA simulations with the electron energy transport cut > maximum energy in the problem (i.e. no transport performed). We did this for monoenergetic electrons of 1 MeV, 100 keV, and 10 keV. For these simulations we expected the energy deposited in the source voxel to equal the mono-energetic electron energy. However, the observed differences without transport were −6.5%, −7.5%, and −13.5% for 1 MeV, 100 keV, and 10 keV, respectively. The differences persisted without transport, providing additional evidence that the energy reaction rate cross sections may be the reason for differences, but further investigation is required to reconcile the differences.
Similar to MC, convergence with the GBBS is not constant throughout space. GBBS absorbed doses converge 'faster' closer to the source as seen in figure 5(a). MC absorbed doses also converge 'faster' for voxels near the source as shown in figure 7. Ignoring variance reduction methods, reducing the statistical uncertainty of voxels far from the source requires simulating many additional source particles with MC. For the GBBS to reduce errors at distance, additional angular discretization (increased angular quadrature order) or a firstscattered distributed source with ray-tracing must be used. For patient-specific calculations, the discretization at distance and corresponding ray effects will be less important because the source will be distributed spatially, and the absorbed dose in a voxel will be dominated by the local activity. Yoriyaz et al (2009) investigated absorbed dose fractions for spheres of material and found maximum differences between MCNP and GEANT4 to be 5% for photons and 10% for electrons. Our GBBS results in the source voxel, for the full 90 Y and 131 I spectra shown in figure 8, are similar in magnitude to the differences they reported. Pacilio et al (2009) estimated differences on the order of 5% in the source voxel and much larger in surrounding voxels. It is important to note that the uncertainties in cross sections increase at lower energies, so different models can give quite different results(International Commission on Radiation Units and Measurement 1984, Berger et al 2005. The magnitude can seem alarming to those with backgrounds in traditional radiation oncology transport, but one has to take these differences in the proper context which we have tried to give by looking at differences between other modern radiation transport codes (figure 8).
Voxel-S-values are typically calculated once and then used in a convolution. As such, these are reference values and their computation time is irrelevant. Forcing at least 6 tetrahedrons per voxel to preserve voxel boundaries is not an efficient use of tetrahedrons. However, such a setup ensured the benchmarking elucidated differences underlying the transport and crosssections used in the GBBS and not geometric differences. Geometry matching was essential for evaluating voxel-S-values which have very large absorbed dose gradients; previous studies of the GBBS did not follow such a methodology.
In addition, we ensured the source voxel was sufficiently converged by comparing multiple phase-space variable discretizations. An important distinction of this work relevant to nuclear medicine is that previous studies exclusively used photon sources with either KERMA reaction rates for sealed sources or energy deposition reaction rates for megavoltage photon beams with partial-coupling. In this work we investigated both electron and photon sources with full-coupling.
We captured source spectral shapes and significant discrete radiations through the energy group discretizations. Previous studies of GBBS for radiotherapy used much coarser energy group structures than we did and reported excellent agreement. This suggests that optimization (minimizing number of energy groups) of group structure may lead to good agreement with MC, but it is unclear if such solutions will be converged in the energy discretization variable. It is likely possible to combine all the emission components (gamma + x-ray, beta, Auger + IC electrons) into a single simulation source for the GBBS with full-coupling. This would be accomplished by mapping the corresponding emission to the correct energy group and weighting the component to account for branching ratios.
One potential advantage of the GBBS that has been realized in radiation oncology has been its speed and accuracy. Providing absorbed dose maps for nuclear medicine scans in a timely matter will help the adoption of voxel-level dosimetry for nuclear medicine departments. Our GBBS calculation times on a 24-core (two 12-core AMD 6174) machine ranged from ≈2 min for the coarsest discretizations to ≈20 min for the intermediate discretizations. The finest calcul ations took up to a few days. The coarsest calculations required a couple gigabytes of RAM whereas the finest required around 60 GB of RAM.
A modern deterministic code such as ATTILA has several adaptable (space & energy) options for controlling levels of discretization. Nuclear medicine absorbed dose calculations depend on radionuclide spectra and spatial distributions of activity, material, and density. An advantage of ATTILA is its use of an unstructured tetrahedral mesh. However, in this work the advantage of accurately approximating general geometries through unstructured tetrahedral meshing was negated; we forced tetrahedrons to not cross voxel boundaries to avoid geometric source specification errors and absorbed dose scoring errors when comparing with MC.
The results of this study lead naturally to several additional areas of investigation including the use of GBBS for patient-specific absorbed doses from unsealed sources. Investigation of GBBS computation time is also warranted and should address the following: (1) relaxing the tetrahedral mesh to not match voxels exactly; (2) coarsening energy group structures for each radionuclide; (3) setting energy group dependent S N and P N ; and (4) combining the Auger + IC, beta, and gamma + x-ray into a single fully-coupled source.
The GBBS ATTILA was benchmarked against MC in the nuclear medicine regime. The magnitude of differences between the GBBS and MC was similar to that seen among modern MC codes; this suggests that the GBBS is a viable alternative to MC for voxel-level absorbed doses in nuclear medicine. Reconciliation of the differences between GBBS and MC at lower energies requires further investigation of energy deposition cross-sections.