The Role of Lung Density in the Voxel-Based Dosimetry of 90 Y-TARE Evaluated with the Voxel S-Value (VSV) Method and Fast Monte Carlo Simulation

.


Introduction
TransArterial RadioEmbolization (TARE), also known as Selective Internal Radiation Therapy (SIRT), is a well-established technique used in order to treat unresectable primary and metastatic liver tumors.TARE involves the direct infusion of 90 Yttrium ( 90 Y)-labelled microspheres into the tumor via a hepatic artery catheter [1]. 90Y is preferentially taken up by hyper vascular tumors, such as hepatocellular carcinoma (HCC), which receive most of their blood supply from the hepatic artery, as opposed to normal liver parenchyma, which is supplied by the portal system [2].Radioembolization consists of a planning procedure (using a surrogate gamma-emitting radionuclide) followed by a treatment procedure (using the radioactive microspheres) [3,4]; Technetium-99m Macro-Aggregated Albumin ( 99m Tc-MAA) is used as an imaging surrogate of radioactive microspheres, and it is administered at the catheter intended for treatment with microspheres.Since arteriovenous shunts may provide a direct vascular pathway to the lungs for 90 Y microspheres, causing, in certain conditions in the lungs, irradiation [5], the distribution of 99m Tc-MAA is used to calculate the predicted Lung Shunt Fraction (LSF).Radiation-induced Pneumonitis (RP) is in fact one of the major risks of radioembolization procedures.The recommendations, expressed as LSF, do not provide direct information on doses absorbed in the lungs.However, the threshold doses are based on a maximum dose of 30 Gy to a lung mass of 1.0 kg for a single treatment and 50 Gy for multiple treatments [6].These values were set because patients who received a Lung Mean Dose (LMD) greater than 30 Gy for a single treatment or a cumulative lung dose greater than 50 Gy for multiple treatments developed RP [7,8].In order to calculate the absorbed dose, TARE uses the Medical Internal Radiation Dose (MIRD) formalism, whereby all activity is localized within the volume of the perfused liver and the administered activity is converted to the absorbed dose by assuming that 1 GBq of administered activity per kg of water-equivalent tissue (liver) mass gives an absorbed dose of 49.38 ± 0.5 Gy, usually rounded up to 50 Gy (Equation ( 1)).

D[Gy
In the presence of lung shunting, if the total activity is divided between the liver and lungs volumes, the dose to the lungs can be calculated directly from Equation (1), knowing the lungs' mass and the partitioned activity of the lungs (Equation ( 2)).

D[Gy
Traditionally, the MIRD scheme has been applied to organs and sub-organs with the assumption of uniform activity distributions [9,10]; within each organ or major component of an organ, the source was considered to be uniformly distributed.In general, both the tumor and the normal tissues acquire non-uniform activity distributions.More accurate and sophisticated methods to better assess the dose distribution in 90 Y-TARE have been based on this formalism.Several voxel-based dosimetry models have been developed using activity concentrations from treatment imaging including the Local Deposition Model (LDM) [11,12], Dose Point Kernel convolution (DPK) [13,14], Voxel S-Value (VSV) kernel convolution [15,16], and collapsed cone convolution [17,18].However, all these methods assume that 90 Y microspheres release the dose in a water-equivalent medium (i.e., the liver).Since lung density is quite different from liver density, the lung's absorbed dose could vary considerably, especially at the liver/lungs interface.To the best of our knowledge, a few previous studies [19,20] have reported the use of voxel-based dosimetry methods to determine the absorbed dose to the lungs, and have investigated the effects of applying these different methods in the liver/lungs interface region.
In this scenario, a Monte Carlo simulation would certainly be the most robust method for dose estimation; the simulation of radiation transport through any medium is possible using several MC codes, and provides an accurate dose estimation [21].Even if computer hardware is becoming more powerful, absorbed dose calculations based on MC simulations using patient-specific characteristics are still time-consuming, especially when attempting to assign absorbed dose at the voxel level for reporting dose-volume histograms.Despite continuous improvements in hardware and software, fast approaches (DPK and VSV) are still more attractive.Recently, software implementing a modified version of the MC code dose-planning method, optimized to run on Graphics Processing Units (GPUs), have become available [22][23][24].This improves the computational efficiency of MC simulations to facilitate large-scale research use and clinical applications [17,25,26].These software allow us to evaluate the dosimetric effect of using water-equivalent or patient-specific densities in 90 Y TARE dosimetry, especially for the lungs' region.The aim of this work is to compare the dosimetric results obtained by two dedicated software implementing a water-equivalent dose calculation and a fast MC simulation, with special emphasis on possible differences at the lungs' level and their impact on 90 Y-TARE dosimetry.In order to assess this purpose, a first comparison was performed on phantom simulation images, and then on a retrospective selection of 24 patients with primary HCC who are treated with 90 Y-TARE.

Phantom Preparation and Data Acquisition
For this study, an IEC (International Electrotechnical Commission) anthropomorphic torso phantom was used to simulate the thorax and internal anatomy of an adult patient.It consisted of a 1.2 L liver compartment, a 10.3 L background, and two lung compartments, left and right, of 0.9 L and 1.1 L, respectively.Distilled water was used to fill the phantom background and liver compartment.A Siemens Biograph mCT PET/CT was used for the 90 Y data acquisition, made of 4 mm lutetium oxyorthosilicate (LSO)-based crystals and arranged in blocks of 13 × 13 elements and 4 photomultiplier tubes per block.The system consists of 4 detector rings, each formed by 48 blocks with 624 detector elements.The coincidence window was set at 5.4 ns.The axial and transaxial Fields Of View (FOV) were 162 mm and 700 mm, respectively; the CT-extended FOV was 780 mm.The CT scan was a low-dose CT (130 kV and mAs greater than 100).The images were then reconstructed using a 256 × 256 matrix and 3D TOF-OSEM (1 iteration, 21 subsets), and were filtered with a Gaussian function filter of 6 mm Full Width at Half Maximum (FWHM).This reconstruction modality was proved to have the best quantitative accuracy [27].

Mathematical Phantom Simulation
The logical steps followed in the preparation of the phantom simulations are shown in the flow chart reported in Figure 1.Starting from the original phantom images acquired on the PET scanner, the following operations were performed to recreate a more realistic situation.The lung compartment was manually expanded to simulate a scenario in which liver and lungs are adjacent.Then, in a first simulation (Setup 1), only the activity concentration of the liver compartment was manually fixed at 1.38 MBq/mL, while in a second simulation (Setup 2), also the lung compartment was set to have a homogeneous concentration of 0.138 MBq/mL to simulate a 10% LSF clinical situation.Once the phantom images were obtained, the following Volumes Of Interest (VOIs) were drawn: Lungs, whole liver (Liver) and liver-lung edge (L-L Edge).The L-L Edge region was created by expanding the liver structure by 1 cm within the lung structure to reproduce a typical blurred region due to patient's respiratory motion.The expansion value of 1 cm was chosen considering the Continuous Slowing Down Approximation (CSDA) range for 938 keV electrons in water [28].
region due to patient's respiratory motion.The expansion value of 1 cm was chosen considering the Continuous Slowing Down Approximation (CSDA) range for 938 keV electrons in water [28].The first step was to obtain the CT image series, while the second step was to prepare the PET series by masking the original with fixed activity concentrations.

Patient Selection and Data Acquisition
Twenty-four consecutive patients with primary liver tumors (HCC) treated in our institution with 90 Y-TARE, using both glass and resin microspheres, were included in this retrospective study (notified to the ethic committee and approved by institutional review board with prot.N. 1534/24).The demographic characteristics of the cohort are summarized in Table 1.Inclusion criteria for the study were the following: 1.Primary liver tumors (HCC) eligible for 90 Y-TARE, according to an internal tumor board; 2. Pre-treatment high-resolution imaging (triple phase contrast CT or MRI); 3. Lung shunt fraction assessed by 99m Tc-MAA Single Photon Emission Computed Tomography (SPECT) <20%; 4. Absence of abdominal extrahepatic shunts assessed by 99m Tc-MAA SPECT; 5.At least nine months of follow-up; 6. Signed informed consent.Dose assessments were performed using 90 Y-PET/CT images acquired after treatment using the same acquisition protocol as the phantom.The first step was to obtain the CT image series, while the second step was to prepare the PET series by masking the original with fixed activity concentrations.

Patient Selection and Data Acquisition
Twenty-four consecutive patients with primary liver tumors (HCC) treated in our institution with 90 Y-TARE, using both glass and resin microspheres, were included in this retrospective study (notified to the ethic committee and approved by institutional review board with prot.N. 1534/24).The demographic characteristics of the cohort are summarized in Table 1.Inclusion criteria for the study were the following: 1.
Primary liver tumors (HCC) eligible for 90 Y-TARE, according to an internal tumor board; 2.
Pre-treatment high-resolution imaging (triple phase contrast CT or MRI); 3.
At least nine months of follow-up; 6.
Signed informed consent.Dose assessments were performed using 90 Y-PET/CT images acquired after treatment using the same acquisition protocol as the phantom.

Dosimetry Software
MIM 90 Y-SurePlan (version 7.2.8,MIM software Inc., Cleveland, OH, USA) was used to obtain Voxel-S-Value (VSV) approach dose values, in opposition to Torch software (version 1.3.0.34,Voximetry Inc, Madison, WI, USA), which was used to gain dose values in a MC approach.

MIM 90 Y-SurePlan
For β-emitters such as 90Y, MIM Software supports dosimetry using either LDM or VSV kernel convolution, and can be performed on either bremsstrahlung SPECT or PET images.In this work, VSV kernel convolution was used.The VSV convolution method in SurePlan is an approach based on the 90 Y scheme defined in MIRD Pamphlet 17 for non-uniform distribution of radioactivity [29].In MIRD Pamphlet 17, MC methods are used to calculate the dose deposited in 3 × 3 × 3 mm water target voxels surrounding a uniform 90 Y source in the central voxel.The 3 × 3 × 3 mm VSV kernel for 90 Y is convolved with the PET activity concentration image (in Bq/mL), which has been decay-corrected back to the time of injection, to calculate the absorbed dose in units of Gy [30].

Torch
Torch is configured to use the parallel processing capabilities of GPUs to handle the successive steps of image registration, contour propagation, kinetic modelling, and radiation transport.A key component of this workflow is the software's proprietary, GPU-accelerated MC algorithm.A modified version of the MC code dose-planning method, optimized for use on GPUs, was employed; electron transport was performed using the condensed history method, where large energy transfers are treated in an analogue manner and small energy transfers are considered by the continuous deceleration approximation [31].After each step, the angular distribution of the electrons is determined using step size independent multiple scattering theory.Photons are transported with a standard analogue approach, taking into account photoelectric absorption and Compton scattering.

Extracted Data and Statistical Analysis
Mean Absorbed Dose (MAD) of the structures: Lungs, Liver and L-L Edge were extracted for both phantom and patients.Indeed, MAD of the structures: Left Lung (Lung L), Right Lung (Lung R), Healthy Liver (H-L) and tumor target were extracted and compared only for patients.H-L structures were generated by Boolean subtraction of the Tumor contour from the Liver one.The Lungs were segmented using the region growing tool in MIM 90 Y-SurePlan, while the Liver was segmented using an atlas-based automatic contour workflow in MIM 90 Y-SurePlan.All patient contours were reviewed and manually adjusted by a physician.Percentage differences between MAD values obtained using the MC method and those obtained using the VSV method were calculated, taking the MC method as the reference.The agreement between the two methods was evaluated using the Passing-Bablok regression [32].R (analytical software, version 2023.09.0 Build 463) was used for statistical calculations [33].

Mathematical Phantom Simulation Results
Figure 2 shows the differences in absorbed dose distribution between MIM and Torch software for setup 1 and 2. The mean percentage dose differences between VSV and MC were found to be 1.2% and 0.5% (with absolute dose variations of 0.7 Gy and 0.3 Gy, respectively) for the Liver, −56% and −70% (with absolute variations of −0.3 Gy and −16.2 Gy, respectively) for the Lungs, and −48% and −60% (with absolute variations of −4.3 Gy and −16.5 Gy, respectively) for the L-L Edge region in the two simulations.The MAD values are summarized in Table 2 for the first and second setup simulations.

Mathematical Phantom Simulation Results
Figure 2 shows the differences in absorbed dose distribution between MIM and Torch software for setup 1 and 2. The mean percentage dose differences between VSV and MC were found to be 1.2% and 0.5% (with absolute dose variations of 0.7 Gy and 0.3 Gy, respectively) for the Liver, −56% and −70% (with absolute variations of −0.3 Gy and −16.2 Gy, respectively) for the Lungs, and −48% and −60% (with absolute variations of −4.3 Gy and −16.5 Gy, respectively) for the L-L Edge region in the two simulations.The MAD values are summarized in Table 2 for the first and second setup simulations.

Patient Analysis Results
The mean percentage dose differences between VSV and MC were found to be 7.0 ± 0.8%, 4.1 ± 2.9% and 6.7 ± 1.4% for water-equivalent structures: Liver, H-L and Tumor, respectively.For air-equivalent structures, the mean percentage dose differences were found to be −61.2± 7.7% for the Lung L and −61.1 ± 6.2% for both Lung R and Lungs.Considering the L-L Edge structure, the percentage mean dose differences between VSV and MC were −50.8 ± 5.3%.The box plots of the percentage differences between VSV and

Patient Analysis Results
The mean percentage dose differences between VSV and MC were found to be 7.0 ± 0.8%, 4.1 ± 2.9% and 6.7 ± 1.4% for water-equivalent structures: Liver, H-L and Tumor, respectively.For air-equivalent structures, the mean percentage dose differences were found to be −61.2± 7.7% for the Lung L and −61.1 ± 6.2% for both Lung R and Lungs.Considering the L-L Edge structure, the percentage mean dose differences between VSV and MC were −50.8 ± 5.3%.The box plots of the percentage differences between VSV and MC results are shown in Figure 3. Passing-Bablok regression plots for the agreement assessment between the VSV method and MC one for each water-equivalent structure extracted from the patient analysis are shown in Figure 4 and display a mean slope for these structures equal to 1.07 ± 0.01.
Figure 5, instead, shows Passing-Bablok regression plots for the agreement evaluation between the VSV method and MC one for each air-equivalent structure extracted from the patient analysis, displaying a mean slope equal to 0.39 ± 0.01.The L-L Edge structure is also included in Figure 5, with a regression slope of 0.49.More details about Passing-Bablok regression analysis results are listed in Table 3.
MC results are shown in Figure 3. Passing-Bablok regression plots for the agreement assessment between the VSV method and MC one for each water-equivalent structure extracted from the patient analysis are shown in Figure 4 and display a mean slope for these structures equal to 1.07 ± 0.01.
Figure 5, instead, shows Passing-Bablok regression plots for the agreement evaluation between the VSV method and MC one for each air-equivalent structure extracted from the patient analysis, displaying a mean slope equal to 0.39 ± 0.01.The L-L Edge structure is also included in Figure 5, with a regression slope of 0.49.More details about Passing-Bablok regression analysis results are listed in Table 3.

Discussion
This work fits into the 90 Y dosimetry studies by evaluating the differences in absorbed dose when different approaches are applied at the voxel level.The aim was to compare the dosimetric results obtained from two dedicated software packages implementing waterequivalent dose calculations and fast MC simulations, respectively.There are several methods to calculate voxel-based absorbed doses for 90 Y-TARE, but there are not so many works in the literature that compare different techniques and methodologies, especially at the lung level.To achieve this objective, a first phantom analysis was performed, followed by a verification on a selected cohort of patients with the same histological disease (HCC).The results of the phantom analysis show that when the dose is calculated within uniform water-equivalent density regions, the differences between the values obtained using the VSV method and those obtained using the MC method are very small (1.2% and 0.5% for setup 1 and setup 2, respectively).These values are very similar to those reported by Dieudonne et al. [20] and Mikell et al. [19], as well as for differences evaluated in airequivalent density regions.Mikell et al. reported an underestimation of the MAD in the right lung of about −60% and about −70% for the absorbed dose deep in the lung when a soft tissue kernel (i.e., VSV method) is used with respect to the MC.This lung dose underestimation behavior is confirmed for both phantom simulations and patient analysis, with mean percentage differences for lung mean dose of −56% and −67% in the phantom simulations and −61.1% in the patient study.
Statistical analysis through the Passing-Bablok regression method showed that no systematic differences (confidence intervals of the intercept reported in Table 3 contain the value zero) between MC and VSV for almost all structures except for the H-L appear to be present.Moreover, at least a proportional difference (confidence intervals of the slope reported in Table 3 do not contain the value one) between the two methods for all structures is evident.These proportional differences are small for water-equivalent structures (slopes very close to one), and are more significant for air-equivalent structures, where a greater discrepancy between the two methods is observed, in line with the dosimetric evaluations.
Lung dosimetry plays an important role in 90 Y-TARE, and, to our knowledge, very few works have investigated the importance of a correct dosimetric model and its implications at the lung level (including the liver-lungs interface) [34][35][36].There have been major advances in liver and tumor dosimetry models that allow clinicians to plan treatments prospectively with higher administered activities to increase the likelihood of local tumor control after radioembolization.Unfortunately, lung dosimetry and dose constraints are very rudimentary [37].The clinical data that justified lung dose limits of 30 and 50 Gy are based on studies dating back to 25 years ago [37] on planar imaging data for a standard lungs mass of 1 kg [38].A reliable state-of-the-art lungs-absorbed dose limit for radioembolization is not yet available.Using dedicated software that implements a MC code would be important, especially when there are significant density differences and motion artefact issues.In fact, when the correct medium density is applied, the activity at the liver-lungs interface has a significant impact on the lungs' MAD, and locally at the L-L Edge region.
Nevertheless, the phantom study showed that the liver MAD is not so much influenced by the possible presence of activity in the lung: by adding a representative percentage (10%) of activity, the liver MAD changed from 61.40 Gy to 62.50 Gy (an increase of about 1.8%), and from 62.11 Gy to 62.79 Gy (an increase of about 1.1%) when using the MC and VSV methods, respectively.Obviously, the patient's respiratory motion could potentially affect the definition of the liver-lungs boundary region and its effective activity distribution.However, the small difference between the MC and VSV methods in the liver region confirms that VSV can provide accurate dose estimation in cases where the lung shunt is limited and motion artefacts are not so evident.The use of a free-breathing CT scan within the PET/CT could be a potential limitation as well as the accuracy of 90 Y-PET [39].It is well known that the extremely low positron yield of 32 × 10 −6 per 90 Y disintegration [40] has been a unique challenge in 90 Y PET imaging.This limitation is overcome in the phantom study by manually setting the activity concentration in each voxel of the PET image series, but remains in the patient analysis.Furthermore, this work confirmed that at the lungs level the MAD values could potentially vary significantly by using a MC instead of a soft-tissue dosimetric model (i.e., VSV), but no clinical implications were explored and future works must be conducted to include evaluations on the lungs' dose-response constraint, both in the pre-treatment evaluation with 99m Tc-MAA and in the post-treatment evaluation with 90 Y microspheres.In fact, the results of this study highlight that the actual dose constraints of 30 Gy and 50 Gy need to be reviewed and updated with new retrospective studies based on MC calculations, given the differences in the lungs' dose when using an MC instead of a soft-tissue dosimetric model.MC methods are considered the gold standard in dose assessment and the increase in modified versions of MC codes operating on GPUs represents a significant improvement.

Conclusions
Both VSV and MC allowed for accurate radiation dose estimation with small differences (<7%) in regions of uniform water-equivalent density (i.e., within the liver).Larger differences (>50%) were observed for air-equivalent regions in both phantom simulations and patient analysis.Although further studies should be carried out to assess the adequacy of the actual mean dose limits (30 Gy for single treatment and 50 Gy for multiple treatments), these differences should be taken into account when planning treatments, as lungs' dosimetry is fundamental to the 90 Y-TARE procedure.

Figure 1 .
Figure 1.Flowchart describing the logical steps taken to prepare images for phantom simulations.The first step was to obtain the CT image series, while the second step was to prepare the PET series by masking the original with fixed activity concentrations.

Figure 1 .
Figure 1.Flowchart describing the logical steps taken to prepare images for phantom simulations.The first step was to obtain the CT image series, while the second step was to prepare the PET series by masking the original with fixed activity concentrations.

Figure 2 .
Figure 2.An axial view with phantom liver and lungs compartments showing the differences between Setup 1 and 2 when using MIM 90 Y-SurePlan (a,c) and Torch (b,d).

Figure 2 .
Figure 2.An axial view with phantom liver and lungs compartments showing the differences between Setup 1 and 2 when using MIM 90 Y-SurePlan (a,c) and Torch (b,d).

Figure 3 .
Figure 3. Box plots of percentage differences between MAD calculated using the VSV method and MC one for selected VOIs.

Figure 3 .
Figure 3. Box plots of percentage differences between MAD calculated using the VSV method and MC one for selected VOIs.

Figure 4 .
Figure 4. Passing-Bablok regression for comparison of MAD values between VSV method and MC one for water-equivalent structures (a) Liver, (b) H-L, and (c) Tumor.

Figure 4 .
Figure 4. Passing-Bablok regression for comparison of MAD values between VSV method and MC one for water-equivalent structures (a) Liver, (b) H-L, and (c) Tumor.

Figure 5 .
Figure 5. Passing-Bablok regression for comparison of MAD values between VSV method and MC one for air-equivalent structures (a) Lung R, (b) Lung L, (c) Lungs, and (d) L-L Edge.

Figure 5 .
Figure 5. Passing-Bablok regression for comparison of MAD values between VSV method and MC one for air-equivalent structures (a) Lung R, (b) Lung L, (c) Lungs, and (d) L-L Edge.

Table 1 .
Baseline characteristics of the patients selected for the analysis.

Table 1 .
Baseline characteristics of the patients selected for the analysis.

Table 2 .
MAD values extracted from phantom analysis for selected VOIs by using MIM and Torch, respectively.

Table 2 .
MAD values extracted from phantom analysis for selected VOIs by using MIM and Torch, respectively.

Table 3 .
Passing-Bablok regression analysis results reporting for each fitting parameter (intercept and slope) the relative lower and upper confidence interval.

Table 3 .
Passing-Bablok regression analysis results reporting for each fitting parameter (intercept and slope) the relative lower and upper confidence interval.