In Silico Validation of MCID Platform for Monte Carlo-Based Voxel Dosimetry Applied to ⁹⁰Y-Radioembolization of Liver Malignancies

Featured Application

Novel treatment planning system performing Monte Carlo-based voxel dosimetry for a tailored radioembolization of liver malignancies.

Abstract

The aim was the validation of a platform for internal dosimetry, named MCID, based on patient-specific images and direct Monte Carlo (MC) simulations, for radioembolization of liver tumors with ⁹⁰Y-labeled microspheres. CT of real patients were used to create voxelized phantoms with different density and activity maps. SPECT acquisitions were simulated by the SIMIND MC code. Input macros for the GATE/Geant4 code were generated by MCID, loading coregistered morphological and functional images and performing image segmentation. The dosimetric results obtained from the direct MC simulations and from conventional MIRD approach at both organ and voxel level, in condition of homogeneous tissues, were compared, obtaining differences of about 0.3% and within 3%, respectively, whereas differences increased (up to 14%) introducing tissue heterogeneities in phantoms. Mean absorbed dose for spherical regions of different sizes (10 mm ≤ r ≤ 30 mm) from MC code and from OLINDA/EXM were also compared obtaining differences varying in the range 7–69%, which decreased to 2–9% after correcting for partial volume effects (PVEs) from imaging, confirming that differences were mostly due to PVEs, even though a still high difference for the smallest sphere suggested possible source description mismatching. This study validated the MCID platform, which allows the fast implementation of a patient-specific GATE simulation, avoiding complex and time-consuming manual coding. It also points out the relevance of personalized dosimetry, accounting for inhomogeneities, in order to avoid absorbed dose misestimations.

1. Introduction

Radioembolization (RE) is a clinical therapy for the treatment of primary or secondary hepatic tumors. The procedure is based on the administration of ⁹⁰Y-loaded microspheres through the hepatic artery, which is the major blood supplier for liver malignancies [1]. The vascular targeting makes RE extremely selective: high dose can be delivered to neoplastic areas while preserving nearby healthy tissue. Several approaches, empirical and dosimetric, have been proposed to establish the activity to be administered [2]. Due to their simplicity, empirical models have been applied for several years, but at present they are not considered as adequate for patient-specific treatments. In fact, the clinical benefit evidence of dosimetry-based approaches [3,4,5,6,7,8,9,10] has led the scientific community to recognize the importance of accurate absorbed dose evaluation and to focus on the absorbed dose–biological effectiveness relationship.

Different dosimetric approaches for RE have been applied throughout the years: the MIRD scheme at the organ level [11], the partition model [12], the local energy deposition method (LDM) [13,14,15], convolution calculations by voxel S-values [16,17,18,19], and direct Monte Carlo (MC) simulations [20,21]. Methods based on the MIRD approach at the organ level are mostly widespread due to their ease of use, but they present two main limitations: the regions of interest are considered to have homogeneous density, with uniform activity distributions. These unrealistic hypotheses can be overcome introducing the image-based 3D voxel dosimetry. Predictive patient-specific dosimetry is derived simulating the therapeutical procedure with the injection of ^99mTc-macroaggregated albumin particles (^99mTc-MAA) to mimic the ⁹⁰Y-microspheres distribution inside the liver. The invasive procedure and the use of nonidentical particles may affect the perfect correspondence of pre vs. post-therapeutical activity distribution, depending on patient specific vascularity, type of disease, experience, and shrewdness. However, many authors [6,9,22,23] have shown the good representativity of ^99mTc-MAA activity distribution for treatment planning, which is verified afterwards by post ⁹⁰Y-PET or ⁹⁰Y-SPECT imaging. 3D approaches based on LDM or convolution of S-voxel values allow to account for the nonuniform activity distribution derived from the SPECT images, whereas they still assume a homogeneous density. MC simulation, unlike the previous methods, takes into account both tissue inhomogeneities and nonuniform activity distribution and is considered as the gold standard.

In the last two decades, various Monte Carlo-based internal dosimetry programs or routines have been applied [24,25,26,27,28,29]. Despite being considered the gold standard, they all differ in complexity and have not been integrated in the clinical practice yet, except in few centers for research purposes only. The main reason is that they usually require very high computer performances and calculation time, which are not compatible with the daily clinical routine. A second reason is often the lack of validation processes.

This work deals with the physical validation of a novel treatment planning system (TPS) named MCID (MC Internal Dosimetry tool) [30,31], performing Monte Carlo-based voxel dosimetry, applied to ⁹⁰Y-radioembolization of liver malignancies. Patient’s CT and SPECT can be imported in MCID software, which creates different macros, i.e., sequences of scripted commands, for the simulation with GATE/Geant4 [32,33]. Each macro contains various settings about a specific aspect of the simulation (e.g., geometry of the system, source type, physics of the simulation etc. [33]). MCID prepares the macros using some settings defined by the user through the platform (e.g., numbers of primaries, type of radionuclide etc.) and some default settings, such as the definition of the physics of the simulation. Most importantly, the activity source of the simulation is defined through the loaded SPECT of the patient, while the geometry of the simulation is defined after the segmentation of the loaded patient CT performed by the user on the platform (i.e., the morphological image is converted into a density map). Default and user-defined settings can be eventually modified, if necessary, manually changing the macros. These features allow the preparation of a personalized simulation, accounting for the specific patient morphology and activity distribution, in a very short time, avoiding manual coding related difficulties. An additional aim of this work was also to investigate the impact of tissue inhomogeneities on the dosimetric evaluation for the RE treatment and the potential improvement of a MC approach in this therapy.

2. Materials and Methods

For the validation workflow voxelized phantoms were defined, starting from real patient images. The choice of using phantoms was expressly made in order to set activity ground truth and to easily define different density/activity combinations by properly modifying the phantom. It is important to remark that the workflow is strictly related to the validation process with computational phantoms; in future clinical applications, real patient images will be directly loaded on MCID.

2.1. Workflow

The workflow developed for the validation and the evaluation of the TPS can be summarized as follows: (i) phantoms creation, (ii) simulation of the SPECT projections, (iii) attenuation map generation and tomographic reconstruction, (iv) creation of the input file for the MC simulation, (v) conversion of the GATE output file in absorbed dose images, (vi) calculation of absorbed dose images by convolution of voxel S-values (the latter calculated independently from the TPS tested here) for comparison and validation purposes. The software packages (and respective references) used for each step are reported in the corresponding subsection.

2.1.1. Phantoms Creation

The voxelized phantoms used in this study were constructed starting from a real patient CT (512 × 512 matrix with a pixel size of 1.367 × 1.367 mm² and a spacing between slices of 3.27 mm, high quality full diagnostic scan). CT contouring was performed with ITK-SNAP software [34], using a semiautomatic segmentation technique based on 3D active contours. The contoured volumes of interest (VOIs) were liver, soft tissue, lungs, adipose tissue, cortical bone, trabecular bone, air, breast prosthesis. Activity maps, with the same voxel and matrix size of the CT scan, were also constructed using ITK-SNAP, which allows the definition of the activity concentration ratio among different regions.

2.1.2. SPECT Simulation

^99mTc-SPECT projections were simulated with the SIMIND MC code [35] using the contoured CT and the activity map created with ITK-SNAP as input. The main parameters chosen for the simulation are reported in

Table 1. Main parameters set for the SPECT simulation.

Parameter	Value
Crystal material	NaI
Crystal length	39 cm
Crystal thickness	0.9530 cm
Crystal width	53 cm
Energy resolution (at 140 keV)	9.7%
Intrinsic spatial resolution (at 140 keV)	0.32 cm
Final matrix dimension	128 × 128 × 120
Pixel size in the final image	0.41 cm
Photons’ histories per projection angle	108
Collimator type	Hexagonal parallel holes

.

These parameters mimic a MEDISO Anyscan gamma camera, which MCID was initially optimized for. The simulated activity was 100 MBq (as suggested in [1]) and an acquisition time for projection angle of 20 s was simulated. Poissonian noise was added to the simulation in order to obtain a more realistic functional image. A resampled CT (128 × 128) was also obtained.

2.1.3. Map of Attenuation and Reconstruction

Attenuation maps and reconstructed projections were obtained from MCID, using the resampled CT and the projections obtained from SIMIND. The reconstruction was performed using the ordered-subset expectation maximization (OSEM) algorithm, with two iterations and 10 subsets. Corrections for the collimator-detector response (CDR), attenuation, and scatter were applied, the latter obtained through the ESSE algorithm [36]. The reconstructed images had cubic voxels of 4.1 mm size and matrix size of 128 × 128.

2.1.4. Input File Elaboration for GATE and Simulation

A new version of the MCID tool, designed for the GATE/Geant4 10.4 code, was used in this study. MCID is a Windows Executable that can also be run on Windows virtual machines. MCID can generate input files for GATE starting from a morphological image (CT) and a coregistered functional image (SPECT or PET). First, the CT image is converted into a density image, assigning to each voxel density and atomic composition. In our case, the CT image was the resampled one, obtained from SIMIND. Tissue segmentation was realized assigning a constant density value for each HU (Hounsfield Units) range we defined on MCID. The materials selected for the segmentation are reported in the description of the simulated scenarios (§ 2.2). The Spectrum source description was set: ⁹⁰Y β^- spectrum was discretized in steps of 10 keV [37]. 10⁸ primaries were also selected for the simulations. All the files generated by MCID, used as input for GATE, are automatically stored in two folders, named data and mac. Each .mac file contains a group of GATE commands which define the settings of the simulation, such as the geometry and the density of the system, the characteristics of the sources, the physics of the simulation, the features of the output, the level of verbosity, and the visualization settings. Moreover, some .mac files recall specific files stored in the data folder (such as the patient density and activity maps). A list of some .mac files, along with set parameters of most interest, is reported in

Table 2. Some .mac files generated by MCID (Monte Carlo Internal Dosimetry tool) and respective main set parameters.

File Name	Function/Set Parameters
_Y90_spectrum.mac	Recalls all the files, # of primaries
actor.mac	Types of output *
geometry.mac	Geometry and density (CT)
physic.mac	Physics list (emstandard_opt3 **)
sourceY90.mac	Source type (SPECT)

* (absorbed dose map, deposited energy map, squared deposited energy map, uncertainty map); ** (standard electromagnetic library, highly accurate for electrons).

.

The file _Y90_spectrum.mac recalls the other files, therefore the user can start the simulation just running the command gate _Y90_spectrum.mac on the terminal. The simulations were performed on a desktop computer with Intel Core i7 5960X and 3 GHz, 16 GB RAM, and Linux Ubuntu 18.04 LTS.

2.1.5. Conversion of the Output GATE File in Absorbed Dose Images with Requested Units

Gate provides as output absorbed dose maps expressed in Gy units [32,33]. However, to be usable, i.e., quantified coherently, these have to be rescaled considering the effective activity at the injection time, the number of primaries simulated, and the abundancy of emission channels. The user needs to insert a calibration factor to rescale for the original activity at injection time:

$$f=\frac{A_{{0,}^{90}Y}}{C_{tot}\cdot {\lambda }_{{0,}^{90}Y}},$$

(1)

where $A_{{0,}^{90}Y}$ is the ⁹⁰Y activity at the injection time, $C_{tot}$is the total SPECT counts in the whole liver and ${\lambda }_{{0,}^{90}Y}$ is the physical decay constant of ⁹⁰Y, while MCID will automatically rescale for the numbers of primaries and the abundancy of emission channels. The ratio between $A_{{0,}^{90}Y}$ and $C_{tot}$ consists in a relative calibration factor (converting voxel counts into activity) [1,12], whereas considering microspheres permanently trapped in capillaries (as commonly assumed), biokinetics estimation are not needed, and the physical decay constant of ⁹⁰Y is used for calculations of number of decays. The ⁹⁰Y activity at the injection time simulated in this work is 2 GBq.

2.2. Simulated Cases

Voxelized phantoms were designed, simulating three clinical cases with different activity and density maps. The materials selected for the segmentation of the different scenarios are generically reported in

Table 3. Tissues used for the segmentation of validation scenarios.

Tissue/Material	Density (g/cm3)
Air	1.205∙10−3
Adipose tissue	0.920
Compact bone	1.850
Lung	0.296
Soft tissue	1.040

, even though some variations were experienced (according to the scenario) and described in the following. The ICRU atomic composition was always assigned [38].

2.2.1. Uniform Liver (UL) Case

The Uniform Liver (UL) case has liver with homogeneous density, uniform activity inside the liver and no activity outside of it. Liver homogeneity and uniform activity distribution are the hypotheses of the standard MIRD approach at the organ level, so this scenario was considered in order to compare the average absorbed dose calculated by the MC-based TPS with the average absorbed dose calculated by the classical AAPM (American Association of Physicists in Medicine) formula [1]:

$$D\left(Gy\right)=49.38 (Gy kg/GBq)\cdot \frac{A_0 \left(GBq\right)}{m \left(kg\right)}$$

(2)

where $A_0$is the ⁹⁰Y activity at the injection time and $m$ is the liver mass. The multiplicative factor 49.38 $Gy kg/GBq$ accounts for the physical characteristics of ⁹⁰Y (half-life and average energy emitted per nuclear transition) and has a relative statistical uncertainty of 0.1% [1]. No production of bremsstrahlung and energy completely released inside the mass of interest are also assumed for the calculation of the multiplicative factor.

Examples of slices of the contoured CT, activity map and reconstructed SPECT image for this case are reported in Figure 1.

Liver and breast prosthesis regions were considered as soft tissue.

2.2.2. Spherical Regions (SR) Case

The Spherical Regions (SR) case has homogeneous liver and activity in three spherical regions only. These spheres, placed inside the liver, are named “SS” (Small Sphere), “MS” (Medium Sphere), and “BS” (Big Sphere) and their radii are 10, 20, and 30 mm, respectively. Slices of the contoured CT, activity map and reconstructed SPECT image for this case are reported in Figure 2.

We developed this scenario in order to compare the average absorbed dose for each sphere obtained from the GATE output with that obtained from OLINDA-EXM [39] S-factor for ⁹⁰Y and spheres of soft tissues with unit density. From the table of S-factors associated to sphere masses included into OLINDA-EXM, the following relationship between S-factor (expressed in $\frac{mGy}{MBq\cdot s}$) and sphere mass m (expressed in g) was obtained and used to derive the corresponding dosimetric factor for a sphere of arbitrary mass:

$$S\left[\frac{mGy}{MBq\cdot s}\right]=0.12456\cdot m{\left[g\right]}^{-0.97267}$$

(3)

As S-factors are referred to spheres with unit density, the soft tissue density for this scenario was set to 1.00 g/cm³.

2.2.3. Nonuniform Liver (NUL) Case

The Nonuniform Liver (NUL) case can be divided in two subcases:

• NUL-a, presenting a liver with homogeneous density and activity placed inside both spherical regions and liver with activity concentration ratio of 5:1, respectively;
• NUL-b, presenting a liver with nonhomogeneous density and activity placed inside the spherical regions (possible tumor lesions) and liver with activity concentration ratio of 5:1, respectively (Figure 3). For this scenario, four tissues were added to the segmentation list in Table 3: tumor (1.200 g/cm³), water (0.998 g/cm³), liver (1.050 g/cm³), and trabecular bone (1.140 g/cm³). The spheres were segmented in MCID as tumor, with different density as compared to the surrounding liver. The left prosthesis for this subcase was segmented as water.

NUL-a case presents a nonuniform activity distribution and was built to evaluate the absorbed dose differences between the S-voxel convolution method and MC (the TPS).

NUL-b is a realistic scenario and was developed in order to quantify the absorbed dose distribution in tumoral regions with different density as compared to the healthy parenchyma, using the S-voxel convolution method and MC (the TPS).

2.3. Convolution with Voxel S-Values Calculated Independently

For validation purposes, absorbed dose images were also obtained applying convolution technique on the reconstructed SPECT images. The voxel S-values used in the convolution were calculated with the EGSnrc code [17], independently from the tested TPS, considering a 4.1 mm voxel size and ICRU soft tissue [38] as medium (density = 1.04 g/cm³). The S-values calculation was set as described in [17]. The source description was the same as that used for the GATE simulation in the TPS, to increase comparability of dosimetric results. Absorbed dose images obtained by convolution were used to validate and test the TPS. Since the convolution technique is based on the application of MIRD approach at the voxel level, the results obtained with convolution will be labeled as “MIRD” hereinafter.

A schematic representation of the described workflow is reported in Figure 4.

2.4. Statistical Uncertainty on the GATE Simulation

A map of relative statistical uncertainties of the absorbed dose in each voxel is obtained as output from each GATE simulation. These values were used for the assessment of the relative statistical uncertainty on the mean absorbed dose in one region of interest (ROI) through the following expression, according to [40]:

$$U_{roi,rel}=\sqrt{\frac{1}{n}\cdot {\displaystyle \sum }_{k=1}^n{\left(u_{k,rel}\right)}^2}$$

(4)

where $U_{roi,rel}$ is the relative statistical uncertainty on the mean absorbed dose to the ROI, n is the number of the voxels inside the ROI, and $u_{k,rel}$ is the relative statistical uncertainty of the absorbed dose in the voxel k.

3. Results

The computational time for each MC simulation was around 5 h. The relative statistical uncertainty on the absorbed dose value in a single voxel was below 1% for the range 150 Gy–700 Gy, below 2% for the range 50 Gy–150 Gy and below 10% for the range 1 Gy–10 Gy.

3.1. UL Case

The results for the mean absorbed dose calculation with the MIRD approach at the organ level and direct MC simulation are reported in

Table 4. Comparison between average absorbed dose values to the liver, calculated from the Gate image (D_gate) and from the MIRD approach (D_mird) at the organ level. The relative difference (RD) is defined as $=\frac{{\overline{D}}_{gate}-{\overline{D}}_{mird}}{{\overline{D}}_{mird}}.$

Dgate (Gy)	Dmird (Gy)	RD (%)
36.79	36.69	0.27

.

The relative statistical uncertainty on D_gate, assessed using (4), was 2% while the relative statistical uncertainty on D_mird, calculated from the uncertainty on the factor 49.38 $Gy kg/GBq$ in (2) (§ 2.2.1), was 0.1%. Therefore, the mean absorbed dose values to the liver obtained from the two methods are compatible.

3.2. SR Case

The comparison between average absorbed dose values for each sphere is reported in

Table 5. Mean absorbed dose values for each sphere, calculated from the Gate image (D_gate) and from the application of the OLINDA/EXM factors (D_olinda). The relative difference (RD) is defined as $\mathrm{RD}=\frac{{\overline{D}}_{gate}-{\overline{D}}_{olinda}}{{\overline{D}}_{olinda}}$.

	Dgate (Gy)	Dolinda (Gy)	RD (%)
BS	666	625	6.57
MS	558	604	−7.62
SS	176	571	−69.20

.

The relative statistical uncertainties on the mean absorbed dose from GATE simulations for BS, MS and SS are 0.5%, 0.6% and 0.8%, respectively. Statistical uncertainties for OLINDA/EXM S-factors were not available.

In this case, differences are more evident, above all for SS, whose result is probably affected by partial volume effects (PVEs). In order to verify this assumption and reduce these effects, the initial activity map created with ITK-SNAP was used as input for MC simulation, skipping SPECT simulation and reconstruction. The obtained results are reported in

Table 6. Mean absorbed dose values for each sphere, calculated from the Gate image (D_gate) with a reduced workflow and from the application of the OLINDA/EXM factors (D_olinda). The relative difference (RD) is defined as $\mathrm{RD}=\frac{{\overline{D}}_{gate}-{\overline{D}}_{olinda}}{{\overline{D}}_{olinda}}$.

	Dgate (Gy)	Dolinda (Gy)	RD (%)
BS	613	625	−1.92
MS	588	604	−2.65
SS	517	571	−9.46

.

The relative statistical uncertainties on the mean absorbed dose in the three spheres obtained from the GATE simulations was 1%.

3.3. NUL Case

Several absorbed dose profiles for each subcase were extracted from different transversal slices. One example profile for the NUL-a case, presenting homogeneous liver and activity in three spherical regions and liver in concentration ratio 5:1, respectively, is reported in Figure 5. An image of voxel-by-voxel relative difference (RD) for the same slice from which the profile was selected is also reported in Figure 5.

All the profiles selected for the NUL-a case showed a relative difference within 3% between the absorbed dose images calculated by MC simulation and convolution of voxel S-values. Relative differences for the entire liver confirmed that this result is valid for all liver slices, except for some boundary voxels, actually external to the liver and characterized by low dose values (less than few Grays).

Absorbed dose profiles for the NUL-b case, presenting nonhomogeneous liver and nonuniform activity, are reported in Figure 6 and Figure 7.

The relative differences between the absorbed dose images were up to 14% in the spherical regions, having a different density (1.200 g/cm³) as compared to the surrounding liver (1.050 g/cm³).

4. Discussion

The validation of MCID platform is demonstrated at both organ and voxel level. In particular, for the UL scenario, the comparison between mean absorbed doses to liver assessed with the MIRD approach at the organ level and with the MC-based TPS showed a very good agreement (RD = 0.27%). This result highlights the efficiency of the developed dosimetric routine: under the hypothesis of homogeneous density and uniform activity distribution in each tissue, both methods are equivalent, as expected. It is interesting to point out that D_mird is a merely theoretical quantity, while D_gate depends on the image quality, e.g., partial volume effects (PVEs), which in this first case appear negligible due to the big size of the observed object (i.e., the whole liver). The effects due to image blurring become relevant instead when dealing with smaller objects, as in the SR case. This scenario allowed a comparison between the average absorbed doses to each sphere. While the BS and the MS present a RD < 8% (absolute value) between the two methods, the SS shows a dramatic RD of −69.2%, caused by the PVEs affecting the SPECT simulation with SIMIND. One contribution to PVEs rises from the matrix resampling (512 × 512 to 128 × 128), an additional contribution derives from the impulsive response of the imaging system, in the SIMIND simulation: spill-out effects affecting voxels cause a change in the activity quantified by imaging, only partially recovered with CDR corrections during the tomographic reconstruction. D_gate depends on the activity estimated from the image, while D_mird is only related to the theoretical initial activity value. The drastic impact of PVEs is evidenced by the results obtained without the SIMIND simulation: for the BS and the MS, RD reduced to 2% and 3% (absolute values), respectively, while for the SS, RD reduced to 9% (absolute value). Remaining discrepancies could be related to some differences between the MC codes used for calculating the OLINDA-EXM sphere S-factors and the updated code used in this work, and to possible mismatching in source description. Relative differences for the BS and the MS are in agreement with the results shown in [41], where the authors presented an analogous comparison for various sphere diameters: their smallest sphere had a 27 mm diameter, resulting in RD = −5%, so the RD here obtained for SS (20 mm diameter, see Table 6) also appears reasonable and comparable. Dosimetry in lesions with size equal (or less) to the maximum ⁹⁰Y β^- range (around 12 mm in soft tissue) should be treated carefully, also due to limited resolving power of SPECT imaging. An analysis of dosimetry in small lesions and a correction factor for the MIRD standard equation are proposed in [42].

A validation of the TPS at the voxel level was presented in the NUL-a case. The differences between the two methods (MC simulation and S-voxel convolution) for each absorbed dose profile were always within 3%.

Finally, the NUL-b scenario shows the importance of a personalized dosimetry in heterogeneous tissues, accounting for both activity distribution and density inhomogeneities, which are often overlooked in present internal dosimetry evaluations. The analysis of different absorbed dose profiles reveals that in the tumoral spherical regions absorbed dose values derived from the MC simulation are lower (up to 14% in absolute value) than absorbed dose values calculated with the convolution method, due to a density increase in lesions of about 20% with respect to the surrounding tissue. Therefore, it is necessary to include morphological patient-specific information in the treatment planning system, including careful calibration of the CT and possibly high quality CT systems (to allow HU-based density estimation for each voxel) to avoid inaccuracies based on the assumption of homogeneous tissues. The use of the highest possible quality CT could be of special value to improve the information especially in inhomogeneous tumors or e.g., in bone metastases.

The concept of voxel dosimetry is bright and suitable, but it is a highly complex association between image reconstruction, segmentation, density, PVE, activity recovery, and absorbed dose calculation. All these issues can concur to limitations which still need to be well understood and solved.

The scientific community is investing appreciable efforts to improve and assess the reliability and accuracy of dosimetry at the voxel level in volumes of interest of various scenarios.

Concerning image quantification, relevant studies are ongoing to highlight and compensate imaging and reconstruction deficiencies that may lead to unrealistic absorbed dose distributions for different organ substructures, lesions or voxel dimensions [43].

As regards the potential impact of density and inhomogeneities, the results of this study, although theoretical, based on a MC approach represent a proof of concept and challenge analysis of more complex situations with real patient data of different clinical situations. This is in fact the topic of a current study of some of the authors, and preliminary results will be soon presented.

5. Conclusions

The validated MCID platform allows the fast implementation of a personalized MC dosimetry, based on patient imaging data, avoiding complex and time-consuming manual coding to the user. The system could be easily integrated in the clinical practice, considering the total computational time of the MC simulation (~ 5 h), suitable with clinical routine. For liver with homogenous density, the comparison between absorbed dose assessed from MC simulations and from the MIRD approach at the organ and voxel level showed a very good agreement. For smaller regions, PVEs from imaging deeply influence the final absorbed dose evaluation, therefore a careful analysis and PVEs corrections are required. Finally, patient heterogeneities should be considered in the treatment planning in order to obtain accurate dosimetric estimations.