Updating 90 Y Voxel S-Values including internal Bremsstrahlung: Monte Carlo study and development of an analytical model

Purpose: Internal Bremsstrahlung (IB) is a process accompanying β -decay but neglected in Voxel S-Values (VSVs) calculation. Aims of this work were to calculate, through Monte Carlo (MC) simulation, updated 90 Y-VSVs including IB, and to develop an analytical model to evaluate 90 Y-VSVs for any voxel size of practical interest. Methods: GATE (Geant4 Application for Tomographic Emission) was employed for simulating voxelized geometries of soft tissue, with voxels sides l ranging from 2 to 6 mm, in steps of 0.5 mm. The central voxel was set as a homogeneous source of 90 Y when IB photons are not modelled. For each l , the VSVs were computed for 90 Y decays alone and for 90 Y + IB. The analytical model was then built through fitting procedures of the VSVs including IB contribution. Results: Comparing GATE-VSVs with and without IB, differences between + 25% and + 30% were found for distances from the central voxel larger than the maximum β -range. The analytical model showed an agreement with MC simulations within ± 5% in the central voxel and in the Bremsstrahlung tails, for any l value examined, and relative differences lower than ± 40%, for other distances from the source. Conclusions: The presented 90 Y-VSVs include for the first time the contribution due to IB, thus providing a more accurate set of dosimetric factors for three-dimensional internal dosimetry of 90 Y-labelled radiopharmaceuticals and medical devices. Furthermore, the analytical model constitutes an easy and fast alternative approach for 90 Y-VSVs estimation for non-standard voxel dimensions.


Introduction
Internal dosimetry for nuclear medicine therapies has remarkable importance in the assessment of the efficacy of treatments and for the safety of organs at risk.Voxel-level three-dimensional (3D) dosimetry, in particular, enables not only the estimation of the average absorbed dose to volumes of interest (VOIs), but also the calculation of metrics, such as Dose Volume Histograms (DVHs) and isodose contours, giving information on the inhomogeneous spatial distribution of absorbed dose at sub-organ level [1,2].
The main calculation approaches currently employed for 3D internal dosimetry are: the Local Energy Deposition (LED) scheme [3]; the convolution of an activity map (or time-integrated activity map) with VSVs or with Dose Point-Kernels (DPKs) [4,5]; the direct MC simulation [6].Direct MC simulation is in principle the most accurate approach, but requires demanding computing resources [7,8]; VSVs convolution is consequently the most widely used approach in clinical practice [9], because of its calculation rapidity while guaranteeing adequate accuracy, especially with the increasingly refined methods proposed to account for the density inhomogeneities of the patient's body [10,11].
VSVs are defined as a 3D matrix representing the average absorbed dose per decay event in each voxel of a regular voxelized geometry of homogeneous material, immersed in an ideally infinite volume of the same material, caused by a radiation source homogeneously emitting in the central voxel, i.e. the voxel set in the center of the semi-infinite volume.VSVs are usually evaluated through a direct MC simulation of the entire emission spectrum of the radionuclide of interest, or through a calculation based on a set of simulations carried out with monoenergetic electrons and photons, at multiple energies [12]. 90Y is one of the most employed isotopes for RadioNuclide Therapies (RNT) [13]; it is an almost-pure β -emitter having a physical half-life of 64.2 h, an average emission energy of 0.93 MeV and an end-point energy of the beta spectrum of 2.28 MeV, corresponding to a maximum beta minus range of ~ 11 mm in water, and similarly in soft tissue [14,15].In addition, thanks to the emission of annihilation photons due to internal β + /β -pair production (with a branching ratio of 31.86 ± 0.47 × 10 − 6 [16]) and to the bremsstrahlung produced by its β particles, 90 Y enables tomographic imaging with PET/CT and SPECT/CT, respectively [17].Consequently, 90 Y is proficiently used in Selective Internal Radiation Therapy (SIRT) of hepatic carcinomas, also technically defined Trans-Arterial Radio-Embolization (TARE), in which resin or glass microspheres labelled with 90 Y are administered via arterial route to permanently reach the capillaries of the lesions [18,19].Other RNT treatments also exploit 90 Y, e.g.Peptide Radio-Receptorial Therapy (PRRT) of NeuroEndocrine Tumours (NETs) with 90 Y-DOTATOC [20] and radioimmunotherapy for hematologic malignancies with 90 Y-MoAbs (Zeva-lin®) [21].The internal dosimetry for such treatments has a great impact on their efficacy, since it affects the estimation of optimised administered activities from pre-therapy dosimetry [22,23], the assessment of pre-therapy estimations through post-therapy dosimetry [24,25], and the successive study of dose-effect correlations [26].
It has been recently highlighted that in dosimetric and radioprotection calculations for beta emitters, such as 90 Y, the phenomenon of IB is systematically neglected [27].To the knowledge of the authors, no Monte Carlo simulation software currently takes into account this process.IB consists in the emission of photons with a continuous energy spectrum simultaneously to the β decay, as a consequence of the interaction of the β particle with the electromagnetic field of the parent nucleus.Recent works showed that, for some irradiation scenarios, the contribution of IB to the deposited energy surrounding a radioactive source can be relevant for certain high-energy β-emitters, namely 32 P and 90 Y. IB photons can contribute up to 20% of the absorbed dose to the extremities of an operator handling an 90 Y source [27]; activity measurements on vial sources show a good agreement with MC simulations only if IB is included, whereas neglecting IB leads to an underestimation of the signal up to − 14% for 90 Y and − 17% for 32 P [28,29].
In the light of the above, the aims of this work were: i) to evaluate, via MC simulations, updated VSVs for 90 Y, including the IB contribution hitherto neglected; ii) to develop an analytical model to quickly calculate 90 Y VSVs, accounting for IB, for any cubic voxel dimension of practical interest.
Regarding point ii), it should be noted that the convolution of a VSVs matrix with an activity matrix requires that both have the same voxel size.Since SPECT and PET scanners reconstruct activity matrices in several voxel sizes, VSVs of that same dimension are needed, but generally only few standard values are published in literature [30][31][32].Consequently, VSVs for specific voxel dimensions should be calculated with direct MC.Alternatively, the activity matrix can be resampled to the VSVs matrix taken from literature, following an approach employed by some commercial software such as MIM [33].Another way is to adopt an analytical method deriving the VSVs for whatever voxel dimension.In Fernández et al. [34] some analytical methodologies were examined, based on the down-sampling of high-resolution VSVs, on VSVs interpolation and fits, affirming their promising accuracy and feasibility, but unfortunately without providing tabular data nor formulas to easily reproduce all their findings.Amato et al. [35] developed a fully documented method based on fitting MC-derived VSVs for monoenergetic electrons and photons as a function of the normalised distance from the source, interpolating the obtained parameters as a function of voxel dimension and energy, and merging the results to deduce VSVs for radionuclides.However, the accuracy of this method may be considered suboptimal when compared to direct MC results, depending on the considered isotope.Recently we proposed an analytical model for 177 Lu VSVs based on fitting the MC-derived VSVs for this isotope and successively fitting the obtained parameters as a function of voxel dimensions, obtaining very good agreement with direct MC results [32].Therefore, we developed a model for 90 Y, with the particular motivation of using as starting data for the model the updated MC-derived VSVs accounting for IB.

Settings of the Monte Carlo simulations
The first purpose of this study was the evaluation, through MC simulation, of 90 Y VSVs including IB contribution, for multiple voxel sizes.To this aim, GATE [36,37] was employed.GATE is a simulation toolkit based on the well validated GEANT4 (Geometry And Tracking 4) software [38][39][40] and widely exploited in the field of radiation applications in medical physics [41], including internal dosimetry [7,23,42,43].Specifically, GATE v9.1, relying on GEANT4 v10.07, was used.These codes do not take into account IB photons emitted together with the beta particles, as confirmed by the absence of any mention in the user documentation.
The VSVs for different voxel dimensions were calculated, using a cubic World having 1.0 m of side and made of soft tissue (material "TISSUE, SOFT ICRP", elemental composition from https://physics.nist.gov/, density 1.03 g⋅cm − 3 ).
Nine cubic voxelized geometries, defined as.mhd files, were generated by homemade Python 1 scripts, all made of 15✕15✕15 voxels, having voxel sides l ranging between 2 mm and 6 mm in steps of 0.5 mm.These files were converted in GATE into voxelized phantom volumes made of the same soft tissue material, through the Image-NestedParametrisedVolume algorithm [44].For each l value, two independent simulations, "Y" and "IB", were implemented, setting in both of them the central voxel as a homogeneous and isotropic particle source: • "Y" was aimed at simulating the "standard" 90 Y decay, not accounting for IB contribution.It was done by setting "ion" type primary radiation and activating the 90 Y decay with the G4RadioactiveDecay module, which implements the radionuclide emission spectra according to the ENSDF database [44,45].• "IB" was aimed at simulating only the IB contribution to the 90 Y decay, by setting photons as primary radiation, with the 90 Y IB energy spectrum defined in the model proposed by Italiano et al. [27] and validated by Auditore et al. [28].In these works, the IB spectrum was built from the experimental data by Venkataramaiah et al. [46], extrapolated at low energies with the theoretical model by Ford and Martin [47], and fitted according to Walrand et al. [48]: where E is the IB photon energy, E 0 is the end-point energy of the 90 Y beta spectrum (2.28 MeV) and the fit parameters are: a = 25.9, b = 10.0, c = 49.4,β = 0.141, γ = 2.84 [27,28].In the present work, the fit function was sampled with constant binning of 0.01 MeV to build an histogram, reported in Fig. 1, as it was done in [28].This histogram was used as input in GATE to set the spectral distribution of the photons and the source sampling probability, through the UserSpectrum module.
For both Y and IB simulations, the interactions of the emitted radiation with matter were simulated with the G4EmStandardPhysics_option4 Physics List [49], setting range cuts of 0.1 mm for the production of secondary particles from all the active processes (corresponding, in soft tissue, to energy thresholds of about 1 keV for photons and 87 keV for electrons), a value at least one order of magnitude lower than all the considered voxel dimensions of the VSVs geometrical setups, ensuring an accurate sampling of the energy spatial distribution.For each simulation, the absorbed dose (AD) and its statistical uncertainty, in terms of standard deviation of the mean [50], were scored in each voxel with GATE DoseActor.10 9 events for each simulation of both types were run, leading to AD uncertainties below 6% in each voxel for Y simulations and below 2% for IB simulations, for all the voxel sizes considered.A local workstation provided with Intel(R) Core(TM) i7-10700 K @ 3.80 GHz CPUs and 32 GB RAM was employed to run the simulations, spending from 50 h to 75 h for Y simulations and between 10 h and 15 h for IB simulations.

VSVs calculation from simulations outputs
The output AD matrices produced by MC simulations were saved in .mhd format and processed by home-made Python scripts for the following purposes: i) labelling each voxel with indices (i,j,k) ranging from (-7,-7,-7) to (7,7,7).ii) for each voxel size l, appropriately merge the AD of Y and IB simulations in each voxel, to account for IB contribution, by adding to the AD of Y the AD of IB times the integrated probability of IB emission for 90 Y (8.48⋅10 -3 photons/β-decay [27]).
The statistical uncertainties of the obtained total ADs were also evaluated, following the guidelines of Chetty et al. [50].iii) to convert the AD values (Gy) into S-Values (mGy/MBq⋅s), both for Y simulations and the ones including IB contribution as described in step ii), that will be referred to as "Y + IB".iv) to average the S-Values of the symmetrical voxels with respect to the centre of the central voxel, and re-calculate the statistical uncertainty (in terms of standard deviation of the mean) in each voxel accordingly.v) to produce text outputs with the S-Values reported in column as a function of the voxel indices and of the so called "normalised radius" R n , a dimensionless variable defined as: with R indicating the absolute distance from the centre of the (i,j,k) voxel to the centre of the central voxel (0,0,0).

Comparison between VSVs
In order to verify the reliability of the performed simulations and of the followed workflow, the S-Values calculated in this work with GATE MC -S GATE (i,j,k) -were compared with the S-Values published by Lanconelli et al. [31] -S Lan (i,j,k) -, freely available at https://www.medphys.it/down_svoxel.htm.S Lan (i,j,k) values had been calculated with the DOSXYZnrc code and had been in turn successfully compared with the S-Values calculated with MCNP4c and PENELOPE [31].The comparison with the results of the present work was carried out in terms of relative percent differences δ(i,j,k) (%) with respect to Lanconelli's S-Values in each voxel, for the l values available in both the databases, i.e. 3, 4, 5 and 6 mm: Since we performed the symmetric averaging of VSVs with respect to the centre of the central voxel (as described in step iv of Sec.2.1.2),we applied the symmetric averaging also to Lanconelli's VSVs before calculating the δ values.This ensured having a single averaged S-Value for each R n value in both datasets, reducing the possible δ fluctuationsdue to MC statistics -which can occur between different combinations of (i,j,k) coordinates equidistant from the central voxel.
This comparison with Lanconelli's S-Values was conducted both for GATE S-Values accounting for IB -S GATE Y+IB (i,j,k) -, i.e. obtained applying also the step ii) previously described, and for the ones not accounting for IB -S GATE Y (i,j,k) -, for which step ii) was not applied; note that IB was not considered in Lanconelli's calculations, nor in other previous works calculating VSVs.
In order to highlight the impact of including IB among the processes simulated for the calculation of 90 Y VSVs, for each voxel size the GATE S-Values obtained with and without IB contribution were compared in terms of relative percent difference ε(i,j,k) (%) defined as: 1 https://www.python.org/.

Description of the model
The second purpose of this study was the development of an analytical model to easily and rapidly evaluate 90 Y VSVs, including IB contribution, for whatever cubic voxel dimension of practical interest, specifically between 2 mm and 6 mm.The model relies on three successive steps: Fig. 2. Relative percent differences δ (Eq.( 4)) between Lanconelli VSVs [31] and GATE MC VSVs -including and not including IB contribution -, plotted as a function of R n for the l values available for comparison: 3, 4, 5 and 6 mm.
i) for each voxel size l, fitting the S-values obtained via GATE MC simulation including IB contribution (S GATE+IB ) as a function of R n , using Eq. ( 5); ii) fitting the fit parameters obtained from point i) as a function of l, each with a function between the ones that will be described in the following; iii) use the fit parameters obtained from point ii) in the function adopted at point i), so that for any l in the examined range the corresponding S-Values as a function of R n are calculated.
Step i) was carried out, adopting the function: where A 1 , A 2 , x 1 , x 2 , W 1 , W 2 , f and g are the fit parameters for each examined l.The first two terms are Gaussian functions, chosen to reproduce the contribution to the S-Values due to the β-particles emitted by 90 Y. Consequently, A 1 and A 2 represent the amplitude of the Gaussians, x 1 and x 2 their average value, W 1 and W 2 their widths.The third term of Eq. ( 5) is analogous to the one used in Eq. (3) of Pistone et al. [32], and allows to reproduce the contribution of the photons following 90 Y decays, coming both from "external" Bremsstrahlung, i.e. from the interactions of β-particles in the medium surrounding the radioactive source, and by internal Bremsstrahlung, as described in the previous sections.The parameter f determines the amplitude of this contribution, while g determines its shape and slope.
The S GATE Y+IB values were fitted using Eq. ( 5) for all the l values, and then the eight parameters were fitted as a function of l, performing thus step ii).Specifically, A 1 , x 1 , x 2 , W 1 and W 2 were fitted with 5th order polynomial functions; g was fitted with a 4th order polynomial function; A 2 was fitted with a mono-exponential function of the form y(l) = y 0 + a 0 ⋅ exp(-l/t 0 ); f was fitted with a bi-exponential function of the form y (l) = y 0 + a 0 ⋅ exp(-l/t 0 ) + a 1 ⋅ exp(-l/t 1 ).All the fits, both in step i) and ii), were performed with the software QtiPlot2 ver.0.9.8.9, using the Nelder-Mead Simplex algorithm.
Step iii) was implemented on a spreadsheet, by defining, for each voxel in the positive octant (i.e. for all the combinations of i,j,k from 0 to 7), the analytical S-Values -S An -as a function of the values of A 1 , A 2 , x 1 , x 2 , W 1 , W 2 , f and g, defined in turn as a function of the fit parameters deduced from step ii).
An optional density correction, according to the scheme described in Sec.2.4 of [51], was implemented in the spreadsheet.In this way, density-corrected S-Values, S An ρ , for the user-selected density ρ (g⋅cm − 3 ) can be calculated as follows: where 1.03 (g⋅cm − 3 ) is the density of the soft tissue used in the GATE MC simulations, described in section 2.1.1.

Comparison of the model with MC
In order to assess the accuracy of the proposed analytical model with respect to the MC results from which it was built, the S An values were compared with the S GATE+IB values for all the l values used in the simulations.Four additional l values were randomly selected: 2.39 mm, 3.84 mm, 4.41 mm and 5.19 mm; the MC simulations and the workflow described in Sec.2.1.2were carried out also for these l values, to test the accuracy of the model for voxel sizes different from the ones used to build the model itself, for either small, medium and large voxel sizes of practical interest.For all the l sizes considered, the comparison between analytical model results and MC results were made in terms of relative percent differences Δ(i,j,k) (%):

90 Y VSVs from GATE Monte Carlo simulations
The 90 Y Voxel S-Values obtained with the GATE MC simulations performed in this work are reported in the Supplementary materials in tabular form, for all the considered voxel dimensions l.As described in Sec.2.1.3,for the available l dimensions of 3, 4, 5 and 6 mm, both the GATE VSVs with and without the IB contribution were compared with the VSVs obtained by Lanconelli et al. [31].Fig. 2 reports the relative percent differences δ (Eq.( 3)) as a function of R n .For the S GATE Y values, without IB, δ lies within ± 5% for almost all the range of R n values, except for some isolated distances, depending on the l dimension.Considering that, as anticipated in Sec.2.1.1,the statistical uncertainties on S GATE Y values are within 6% for any R n value, the δ values obtained A very similar trend is observed for the relative percent difference ε (Eq.( 4)) between GATE MC VSVs including and not including IB contribution, reported in Fig. 3 as a function of R n : from R n > 2.5 in the case of l = 6.0 mm, up to R n > 6 in the case of l = 2.0 mm, ε reaches values between + 25% and + 30%, increasing up to + 27% and + 32% in correspondence to the farthest R n values considered, respectively for l = 2.0 mm and l = 6.0 mm.

90 Y VSVs from analytical model and comparison with Monte Carlo
The results of step i) described in Sec.2.2.1 are represented in the top panel of Fig. 4, where the GATE MC VSVs including IB contribution are plotted as a function of R n , together with the respective fits performed using Eq. ( 5).Error bars on the GATE MC VSVs were omitted for the sake of clarity, also because the final statistical uncertainties on S GATE Y+IB below 1% for all the data points, as already stated in Sec.3.1.
The bottom panel of Fig. 4 shows the comparison between the S-Values obtained from the analytical model and the ones from GATE MC, for the l values used to build the model; for each l the relative percent differences Δ (Eq.( 7) of Sec.2.2.2) are plotted as a function of R n .Δ in the central voxel lies within ± 5% for all the l values.For R n lower than the beginning of the photon tails, Δ values fall within ± 25% in general, with some exceptions for l ≥ 4.0 mm, for which Δ reaches up to about ± 40% for few R n values.Δ lies within ± 5% in correspondence to the photon tails of the VSVs, which, as already observed in Sec 3.1, start from different R n values depending on the voxel size l, from R n ≃ 2.5 for l = 6.0 mm, to R n ≃ 6 in the case of l = 2.0 mm.
All the fits in Fig. 4 converged with R 2 > 0.99, and the values of the fit parameters obtained are shown in Fig. 5 as a function of l, and are reported in Table A1 of the Annex.In the same figure the results of step ii) of Sec.2.2.1 are also reported, i.e., the fits as a function of l of the  5).(b) relative percent differences Δ (Eq.( 7)) between GATE MC VSVs including IB and VSVs evaluated with the analytical model proposed in this paper.Fig. 6 shows the comparison between analytical model and GATE MC for the additional l values introduced in Sec.2.2.2, plotting the respective VSVs vs. R n and their relative percent differences Δ.The behaviour of Δ for these l values is in line with the ones already described: Δ < ±5% in the central voxel and in the Bremsstrahlung tails, Δ < ±40% in the remaining R n values.Fig. 7 reports the Δ values of the comparison between the 90 Y VSVs from the pre-existing analytical model by Amato et al. [35] and the VSVs from GATE MC including IB.In this case, Δ lies between + 35% and − 70% for R n lower than the beginning of the tail, between − 5% and − 40% in the photon tail, and Δ(0,0,0) is within ± 5%.

Discussion
The comparison of S GATE Y and S Lan , both not accounting for IB, is fairly satisfactory, since δ < ±5% for almost all the R n considered in this study.Some exceptions are found for sporadic values of R n , for which δ is higher; this mostly happens in correspondence to R ≃ 11 mm, i.e., the value of 90 Y β-particles maximum range.These discrepancies can be attributed to differences in the code and simulation settings between GATE 9.1, used in this work, and DOSXYZnrc, used by Lanconelli et al. [31], namely: the different cross section data used by the two codes for physics processes; the size of the cubic world surrounding the voxelized scoring volume; the range cut for electrons and photons; the fact that in Lanconelli et al. monochromatic electron and photon sources were simulated and their results were then merged to reproduce radionuclide emission spectra.The computation of the energy deposition can be more sensitive to such differences at the mentioned distance from the source since the electronic contribution -dominating at low R and R n -fades, while the photon contribution dominates.Nevertheless, the overall reliability of GATE MC VSVs can be considered verified, in view of the good agreement for almost all R n values.
From the comparison between S GATE Y+IB vs. S Lan , the IB contribution to the 90 Y VSVs appears negligible for R n that, depending on l, correspond to R < 11 mm, the maximum range of 90 Y β particles; indeed, for these R n distances, δ values obtained comparing S GATE Y+IB vs. S Lan remain unchanged.δ remains unchanged also when comparing S GATE  7)).

D. Pistone et al.
Y+IB vs. S GATE Y , while the ε values of this latter comparison are nearly 0%.These results show that, within the range of β particles, the energy deposition is dominated by the βs' themselves, and adding the photon IB contribution, considered its integrated probability per β-decay, does not modify that evidence.On the contrary, the IB contribution becomes relevant for R n values corresponding to R ≥ 11 mm, for which the δ values of the comparison S GATE Y+IB vs. S Lan systematically rise up to + 30% within the R n values considered in this study; ε also increases up to between + 25% and + 30% for these R n values, slightly depending on l.Therefore, we can infer that IB significantly contributes to the Bremsstrahlung tail of VSVs, in addition to the external Bremsstrahlung contribution, already considered in the previous VSVs calculations published in literature.In view of these results, the updated 90 Y VSVs accounting for IB presented in this work improve the accuracy of the dosimetric estimations at distances where the photon contribution is dominant, i.e. for R ≥ 11 mm.Moreover, these results are in agreement with the work of Auditore et al. [52], investigating the impact of IB in the calculation of 90 Y DPKs in water, which showed that for radial distances from the point source larger than the maximum range of βs', DPKs Fig. 7. Relative percent differences Δ (Eq.( 7)) between GATE MC VSVs including IB and the VSVs evaluated with the analytical method by Amato et al. [35], for all the available voxel dimensions l considered in this study (the VSVs by Amato et al. were implemented for l > 3 mm).

Table A1
Values of the fit parameters (Sec.2.2.1), χ 2 /dof and R 2 for the fits of the S GATE Y+IB (R n ) with Eq. ( 5) for all the l values considered in this work.

Table A2
Values of the fit parameters, χ 2 /dof and R 2 for the fits, as a function of l, of A 1 , W 1 , W 2 , x 1 and x 2 with 5th order polynomial (y = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + a 4 x 4 + a 5 x 5 ) and of g with 4th order polynomial, as described in Sec.

Table A3
Values of the fit parameters, χ 2 /dof and R 2 for the fit of A 2 as a function of l with mono-exponential function of the form y(l) = y 0 + a 0 ⋅ exp(-l / t 0 ), as described in Sec. are underestimated between 20% and 34% when neglecting IB contribution.The results of the present work thus confirm the need of including IB among the processes related to β-decay available in MC codes.In the present study we limited to scoring geometries of 15✕15✕15 voxels of sides between 2.0 mm and 6.0 mm, thus having maximum distances from the central voxel between 16 mm and 48 mm.Within these distances, the photon tail of Bremsstrahlung photons (both internal and external) in soft tissue is only partially accounted for.
Estimating the discrepancies between VSVs calculated including or neglecting IB for larger geometries is beyond the scope of this work.
Concerning the analytical model developed to derive the 90 Y VSVs for any voxel dimension between 2 mm and 6 mm, the behaviour of the GATE MC VSVs as a function of R n , deserves a special attention.Since the maximum range of 90 Y βs' is of about 11 mm, thus larger than all the voxel dimensions l considered, the R n at which the electron contribution fades and the photon tail becomes dominant varies with l.Differently, other β-emitting radionuclides having β-ranges smaller than any of the considered voxel dimensions, such as 177 Lu (maximum range of 1.7 mm), exhibit the transition to the photonic tail at the same R n for all the l values [32].In addition, the shape of the β electronic part of 90 Y VSVs as a function of R n changes when varying l, with a steeper decrease for larger l, and smoother decrease with more variations in slope for smaller l.Nonetheless, a smooth variation with l of all the fit parameters of Eq. ( 5) was obtained, as visible in Fig. 5, testifying the reliability of the overall fitting algorithm.
The comparison of 90 Y VSVs from the analytical model (S An ) with the ones from the GATE MC simulations (S GATE Y+IB ) shows a good agreement for the central voxel (Δ(0,0,0) < ±5%), contributing the most to the self-dose deposition, and for the nearest neighbouring voxels (Δ(1,0,0), Δ(1,1,0) and all their symmetric combinations of i,j,k are within 10%), for all the l values examined, both among the ones used to build the model (Fig. 4) and the ones selected for the further verification (Fig. 6).Good agreement is also found for the photon tails (Δ < ±5%); variable and on average inferior agreement is found for the electron part of the VSVs, but in any case Δ is never higher than ± 40% and maximum values are found only for a few R n values.As can also be seen in Fig. 4, these R n values never correspond to the central voxel or the nearest neighbouring voxels, and the corresponding S-Values are one or two orders of magnitude lower than the S Value of the central voxel; consequently, the impact on the absorbed dose estimation on uptaking regions of a patient would be small.
Comparing the VSVs of the pre-existing general analytical model by Amato et al. [35] with the GATE MC VSVs of the present work in terms of Δ, as reported in Fig. 7, Δ is between + 35% and − 70% for R n within βs' maximum range, between − 5% and − 40% in the photonic tail (also given the neglection of IB in Amato et al.), within ± 5% in the central voxel and within ± 16% in the nearest neighbouring voxels.Therefore, the new 90 Y-specific analytical model introduced in the present work produces a general improvement in the agreement with MC with respect to the previous model by Amato et al.
The proposed analytical model shows to be a helpful tool for the calculation of 90 Y VSVs for user-defined non-standard voxel dimensions, whose accuracy compares favourably with the level of uncertainties practically achievable in clinical internal dosimetry [53].It constitutes an alternative to the resampling of activity maps and to the implementation of ad-hoc direct MC simulations for specific voxel dimensions, impractical in centres not having experts of MC simulations.In addition, for the first time it includes the contribution of Internal Bremsstrahlung, a process never considered before in existing analytical approaches to VSVs calculations, to the knowledge of the authors.The model, implemented in a simple spreadsheet, allows a user-friendly calculation of 90 Y VSVs, requiring only to choose the desired voxel dimensions and, at will, a desired density.

Conclusion
In this work, updated 90 Y VSVs were calculated with GATE MC simulation, including the process of Internal Bremsstrahlung, which always accompanies β-decay but is usually neglected in VSVs calculations.Including IB in MC simulations, the estimated VSVs show an increment between + 25% and + 30% with respect to the results obtained neglecting IB; the increment appears at distances from the source>11 mm, i.e., beyond the maximum range of 90 Y β particles, where the photon component is dominant.
An analytical model was also developed as a user-friendly and fast tool for calculating 90 Y VSVs for any voxel dimension l between 2 mm and 6 mm.For any l value examined, the model calculations agree with GATE MC estimations within ± 5% in the central voxel and beyond the maximum β range.This 90 Y-specific model, including for the first time the IB contribution, improved the accuracy with respect to the preexisting general model published by Amato et al. [33].

Table A4
Values of the fit parameters, χ 2 /dof and R 2 for the fit off as a function of l with bi-exponential function of the form y(l) = y 0 + a 0 ⋅ exp(-l / t 0 ) + a 1 ⋅ exp(-l / t 1 ) as described in Sec.

Fig. 3 .
Fig. 3. Relative percent differences ε (Eq.(2)) between GATE MC VSVs including and not including IB contribution, plotted as a function of R n for the different voxel sizes l values examined in this work (from 2 mm to 6 mm).

Fig. 4 .
Fig. 4. (a) Markers represent the S-Values calculated with GATE MC simulations including IB contribution (S GATE Y+IB ), for all the l values examined in the present work; the lines are the respective fits with the function of Eq. (5).(b) relative percent differences Δ (Eq.(7)) between GATE MC VSVs including IB and VSVs evaluated with the analytical model proposed in this paper.

D
.Pistone et al.

D
.Pistone et al.   parameters previously obtained.The parameters of this second set of fits are listed in Tables A2, A3 and A4 of the Annex.The step iii) described in Sec.2.2.1 was implemented in the spreadsheet provided as Supplementary material, which enables to exploit the presented analytical model for the calculation of 90 Y VSVs, including IB contribution, by simply inputting the desired values of l and ρ.