A comparison of methods for adapting 177 Lu dose-voxel-kernels to tissue inhomogeneities

In 177 Lu radionuclide therapies, dosimetry is used for determining patient-individual dose burden. Standard approaches provide whole organ doses only. For assessing dose heterogeneity inside organs, voxel-wise dosimetry based on 3D SPECT/CT imaging could be applied. Often, this is achieved by convolving voxel-wise time-activity-curves with appropriate dose-voxel-kernels (DVK). The DVKs are meant to model dose deposition, and can be more accurate if modelled for the specific tissue type under consideration. In literature, DVKs are often not adapted to these inhomogeneities, or simple approximation schemes are applied. maps, mass-density


Personalized dosimetry
The increasing number of 177 Lu -labelled radionuclide therapies, such as e.g. 177 Lu-labelled somatostatin antagonists and 177 Lu-labelled prostate-specific membrane antigen (PSMA) inhibitors has led to an increased interest in patient and radionuclide-specific dosimetry. This follows a general trend towards personalized medicine. Thus considering patient-specific radionuclide therapy, plain spatially constant soft-tissue kernels have to be replaced by DVKs, which include information about individual differences in tissue heterogeneity , Traino et al 2013, Gustafsson et al 2015.
In clinical practice, positron emission tomography (PET) or single photon emission computed tomography (SPECT) images are employed to assess the 3D distribution of a radiopharmaceutical in a patient, while individual maps of tissue compositions can be derived from x-ray computed tomography (CT). With the advent of hybrid SPECT/CT and PET/CT imaging, it is possible to obtain tissue density information in conjunction with the distribution of radioactivity inside the human body. The spatial resolution of PET and SPECT imaging is in the range of 5-20 mm. Consequently, non-uniformities in the activity distribution on a microscopic scale can be detected. In Sgouros et al (2008), Grimes (2013) an overview of imaging-based patient-specific dosimetry methods is given. For estimating the absorbed radiation dose, the spatial and temporal distribution of radioactivity in the patient's body is determined (Sgouros et al 2004, Jackson et al 2013 and organ-, volume-of-interest-(VOI), or voxel-wise time-activity-curves (TAC) are generated. By temporal integration, a spatial distribution of the number of radioactive decays is obtained. Finally, based on the number of nuclear disintegrations, the absorbed dose can be determined with multiple techniques of varying complexity such as S-values, dose-point-kernels (DPK), dose-voxel-kernels (DVK), and full-scale Monte Carlo (MC) calculations.

Organ-wise dosimetry
The most common method in dosimetry is the application of S-values when it comes to determine radiation doses for large, i.e. 'macroscopic' target volumes like whole organs. In radiation dosimetry in general, hence also when using the Medical Internal Radiation Dose Committee (MIRD) formalism (Bolch et al 1999), sourceand target volumes are defined, in which distributions of radioactivity and absorbed dose are assumed to be isotropic and homogeneous. Consequently, an S-value summarises the complex radiation-matter interactions in a large target volume, e.g. a whole organ, and describes the resulting absorbed energy dose in the target volume per nuclear disintegration in the source volume. Often, tabulated, pre-calculated S-values are used, which were determined from a Monte-Carlo simulation of radiation-matter interactions within a standard male or female body-phantom (Petoussi-Henss et al 2007). Disadvantages of this approach are that spatial variations of the distribution of radioactivity, of tissue-composition, and thus of absorbed dose inside the macroscopic source and target volumes are not taken into account. The concept of dose-kernels extends the MIRD formalism, which estimates dose in entire organs, to smaller, i.e. 'microscopic' source and target volumes. Thus in the limit of ever smaller spatial scales, DVKs ultimately could be regarded as approximating a continuous spatial distribution of point sources (Strifari et al 2006, Dewaraja et al 2010, Lanconelli et al 2012. Nevertheless, organ-wise S-values and voxel-based dose kernels follow the same principle, since in both cases the dose in a specified target volume is computed as the superposition of radiation contributions from each source volume.

Voxel-wise dosimetry
If the source distribution is modelled as a collection of radiation point sources, the energy dose per decay D E = dE 2 /dmdN (Gy/deacy), where E (J) denotes absorbed radiation energy, m (kg) denotes mass and N denotes the number of nuclear disintegrations, caused by an individual isotropic point source in a homogeneous medium, is called a dose-point-kernel (Bochkarev et al 1972, Pérez et al 2011. DPKs represent the radial distribution of absorbed dose around an isotropic point source in an infinite homogeneous propagation medium (Bardiès et al 2003). Thus we are considering a point source and an infinitely extended target region. DPKs can be computed using MC techniques and transform time-integrated-activity (TIA) maps into spatial absorbed-dose distributions. However, convolutions of TIA maps with DPKs are only a good approximation in regions of homogeneous tissue density. A number of publications compare DPK databases generated with different MC simulation software such as GATE, FLUKA, MCNP4C, CGSnrc, and GEANT4/GAMOS , Papadimitroulas et al 2012. In the actual MC simulation, radionuclide-specific (e.g. I-131, Y-90, 177 Lu ) DPKs are determined for specific tissues like e.g. bone, lung, soft-tissue, and water.
Instead of calculating a continuous DPK, a discrete dose-voxel-kernel (DVK) can be determined (Pacilio et al 2009(Pacilio et al , 2015. The voxel size can be chosen arbitrarily, because the kernel can be scaled, as shown in Reiner et al (2009), Dieudonné et al (2010, Fernández et al (2013), Pacilio et al (2015). To obtain the distribution of absorbed dose based on the patient's anatomy, the TIA map, i.e. the spatial distribution of radioactive decays, has to be convolved with various tissue-specific DVKs. Note, however, that DVK-based dose estimations are only appropriate for radionuclides which predominantly emit massive particles of short range in condensed matter. The radioisotope 177 Lu , which is considered here, is a β-emitter, whose electrons have an average range of less than 3 mm in tissue. However, as with every nuclear disintegration, γ-radiation is also emitted, which itself contributes to the absorbed energy dose in the target volume. Unfortunately, the range of γ-radiation is infinite, hence any finite-sized voxel-kernel will inevitably underestimate the amount of absorbed radiation energy in any tissue. Thus DVKs have to balance to diverging requiements: the smaller the voxel scale the higher the spatial resolution and the better tissue heterogeneity can be represented, but the more difficult it is to account for dose contributions stemming from γ-radiation. The latter is frequently scattered and may enter volumes, which are relatively far from the loci of nuclear disintegrations, hence enter the VoI from outside.

Other studies on DVK-dosimetry
To the best of our knowledge, only very few studies, which compare results of dosimetry based on either convolution with DVKs or full Monte-Carlo simulations, are available so far. Many studies aim at introducing new algorithms or technologies, thus have not been evaluated on a larger number of patients. For example, Dieudonné et al presented an algorithm for re-scaling DVKs to arbitrary voxel sizes (Dieudonné et al 2010). They did not adapt for inhomogeneities in tissue-type and -density. They compared the results of their DVK method to a full Monte-Carlo simulation for two patients treated with Y-90 and I-131, respectively. In subsequent work of Dieudonné et al, the authors proposed an algorithm for adapting DVKs to tissue inhomogeneities (Dieudonné et al 2013). Again, results of the DVK-based dosimetry were compared to those of a full MC simulations. Most likely due to the pioneering nature of their study, evaluation was carried out for three patients only, one 131 I, one 177 Lu, and one 90 Y. In another study, which took into account different tissue types, DVKs for several materials were calculated and used for dose calculation (Scarinci et al 2013). Unfortunately, no comparison to a MC gold standard was given, therefore any possible improvement could not be quantified. In the study of Moghadam et al (2016), the DVKs were scaled according to tissue properties in order to account for inhomogeneities and the results of DVK-based dosimetry were compared to those of full MC simulations. However, the method was evaluated on a computational Zubal phantom for three different dose distributions, mimicking 90 Y, 177 Lu, and 32 P TIA activity distributions.

Aims of the study
In summary, so far different DVK methods for 177 Lu radiotherapy have only been evaluated on a very small number of patients or for simple phantoms. Furthermore, a direct comparison of the proposed methods has not been carried out yet. Thus, the error, which could realistically be expected from a dosimetry based on convolution with DVKs, most likely could not be estimated reliably based on existing studies. Additionally, conclusive answers to the questions, which DVK-algorithm achieves best results, and if the additional effort of adapting DVKs to tissue inhomogeneities is justified, could not be provided. Consequently, the aim of the current study is to compare the results of different DVK-based methods for determining patient-and tissue-specific absorbed dose due to 177 Lu radiotherapy in a large number of patients.

Data and methods
The estimation of energy dose D(r T ) = dE(r T )/dm(r T ) J/kg = A(r S ) ⊗ d(r T ← r S ) transmitted by radiation from nuclear disintegrations at source location r S and absorbed at target location r T needs an estimate of two factors: a time-integrated radioactivity A(r S ) and a kernel d(r T ← r S ). The latter represents the energy transfer from the source to the target region, thus summarizing the plethora of physical radiation-matter interactions along the path of the emitted particle or quai-particle. The dose in any volume-of-interest (VOI) is then obtained via a convolution (denoted as ⊗ here) of both factors. This study is about the second factor and about DVK methods for absorbed dose estimation in radiotherapies performed with the radionuclide 177 Lu . The latter represents an electron emitting nuclide (β-emitter) with concomitant photon radiation (γ-radiation). Hence methods applied and results obtained may be of relevance only for radionuclids with similar properties.

Patient cohort
In order to obtain a representative and realistic estimate of the deviations caused by the different DVK methods relative to the reference MC simulation, we decided to carry-out the comparison on real patient data sets instead of artificial phantom data. For this, we randomly selected 26 177 Lu -therapy patients from our clinical routine program which met the following criteria.
• At least one therapeutic application (≈7 GBq) of 177 Lu , either labelled with DOTATOC or PSMA • Availability of SPECT/CT data 24 h p.i.
For the 26 (23 m, 3 f) selected patients, the average age was 64 ± 10 years, ranging from 45 a to 80 a.The average injected radioactivity of 177 Lu was 6.6 ± 0.4 GBq, ranging from 5.8 GBq to 7.3 GBq. Twelve patients have been treated with PSMA-labelled 177 Lu and 14 patients with DOTATOC. The retrospective study has been approved by the Institutional Review Board of the University Hospital Erlangen.

Image acquisition
For the Monte-Carlo simulations, as well as for DVK convolutions, several 3D imaging datasets were used. Generally, they were taken from a hybrid SPECT/CT acquisition at approximately 24 h p.i.

SPECT acquisition
Data were acquired on a Siemens Symbia T2 SPECT/CT system. Acquisition and reconstruction were performed based on our standard quantitative 177 Lu protocol, which is described in detail in Sanders et al (2014). Image calibration was performed using the package xSPECT Quant from Siemens Healthineers. Further details can be found in Tran-Gia et al (2018). Therefore, it is herein only briefly outlined.
• SPECT acquisition using medium energy collimators, 3 • angular sampling, 15 min total acquisition time • Iterative ordered subset expectation-maximisation (OSEM) reconstruction of the 208 keV photopeak data with 16 iterations, 8 subsets, matrix 128 × 128 • Point-spread-function modelling in reconstruction • Triple energy window based scatter correction • CT-based attenuation correction • No post-reconstruction smoothing The SPECT image was used for deriving the number of radioactive decays per voxel. We will refer to this information as Decay Map in the following.

CT acquisition
Additionally, a low-dose CT was carried out for SPECT attenuation correction, calculation of the voxel-wise mass-density, as well as determination of the tissue class. The CT covered the same field-of-view as the SPECT and was acquired and reconstructed using the following parameters.
• Slice collimation of 2 × 5 mm, pitch of 1.8, time per rotation of 0.8 s, tube voltage of 130 kVp, Siemens CareDose 4D tube-current modulation with 30 mAs reference • Filtered-Back Projection (FBP) (Bruyant 2002) reconstruction with B08s and B41s Kernels, 512 × 512 matrix, 2.5 mm slice thickness • The B08s image was used for attenuation correction of the SPECT data • The B41s image was used for deriving Density and Tissue Map The respective methods for calculating Density Map and Tissue Map were identical for both full MC simulations and DVK calculations, and are described below. As CT-images have a higher spatial resolution than SPECT images, the former were re-sampled to a common voxel size of 4.8 × 4.8 × 4.8 mm 3 , which was also used for all MC simulations and DVKs. Hence both sets of images were co-registered to the CT image set and the latter was re-scaled to match the lower spatial resolution of the SPECT images.

Image post-processing
SPECT images were manually co-registered to the corresponding CT images. Images were post-processed to generate spatial distributions (at voxel resolution) of time-integrated activities, simply named decay maps, as deduced from SPECT images, corresponding spatial distributions of mass density, so-called density maps, as deduced from re-scaled CT images and a local characterization of respective tissues as deduced from a classification of voxel mass densities (see table 5 in the appendix).

Decay maps
The number of 177 Lu decays in any voxel was derived for every patient and one therapy cycle from corresponding absolute quantitative SPECT images. Activity densities were obtained in Bq/ml from these SPECT images (e.g. as described in Sanders et al (2014)) at multiple time points, e.g. 4 h, 24 h, 48 h, 72 h, 144 h after administering the radiopharmaceutical. By voxel-wise integration of the respective time-activity curve (TAC), the total number of decays, i.e. the time-integrated activity (TIA), could be obtained for each voxel. The challenges related to the voxel-wise integration of the rather noisy TACs are not the focus of this work (see however (Götz et al 2019a)), rather we assume the decay map as known. Furthermore, we choose to state results as relative differences between the compared methods only and not absolute dose values per injected amount of activity. This is mainly because of strong fluctuations in absolute dose values if estimated voxelwise. The difficulties of voxelwise integration of time activity curves (TACs), deduced from SPECT images, is discussed in Götz et al (2019a) in greater detail, hence is omitted here.
In figure 2, a SPECT image of a representative patient is shown.

Density maps
For calculating patient-specific maps of mass density, the CT part of the SPECT/CT acquisition was used. In every CT, the attenuation coefficient µ of tissue for polychromatic x-rays is measured and usually given in Hounsfield units, which expresses the degree of attenuation relative to water and air (see equation (1)).
Since attenuation is mainly caused by interactions between photons and electrons, µ is strongly correlated to electron density, which itself is correlated to mass density ρ. For this, a CT image can be used to estimate a mass density map, usually by using a piecewise linear relation, with different slopes for HU 0 and HU 0. It can be derived for every CT-scanner by special calibration phantoms. A system-specific CT calibration was not carried out, however. Instead, mass density was calculated as defined in Annkah et al (2014).
where ρ(x) is the mass density in units of g cm −3 belonging to a voxel with HU number x. In figure 1, a CT image slice, reflecting different tissue density values, is illustrated.

Tissue maps
In reality, especially for absorption of beta-particles, not only the mass-density determines absorptive properties but also the varying atomic composition of tissues (e.g. effective mass-and atomic numbers) in every voxel. Unfortunately, an exact determination of these parameters with a single CT acquisition is not possible, but several methods for a rough estimate of tissue type exist. We used a classification based on each voxel's mass density, as proposed in Schneider et al (2000). The assigned tissue-type, composition of tissue-type, and respective range of mass densities can be found in the table 4 and in the appendix in table A1. A representative example can be seen in figure 1.

Reference data
For generating reference data, full Monte-Carlo simulations of radiation transport were carried out for every patient, employing the GEANT4-based GAMOS 5.0.0 toolkit (Arce et al 2014). The procedure is briefly described in the following.
• A voxelized geometry was generated based on the tissue classes and mass-density of each voxel, determined as described before. The size of voxels was (4.8 mm) 3 , identical to that of decay, density, and tissue-type maps. • The decay modes of 177 Lu and emitted gamma and beta energy spectra were defined according to Arce et al (2014). • In total, 1.2 · 10 10 decays of 177 Lu were performed using GAMOS on a high-performance cluster with 176 nodes. • The relative spatial frequency of decays was taken from the decay map, as described before.
• Deposited energy was scored per voxel in units of Gy/decay.

Statistical uncertainty
A general discussion of error computation during dose estimation can be found in a companion paper (Götz et al 2019b) and will not be repeated here. However, because the MC simulations of the 26 patients were considered the gold standard for the respective patient, it was important to know the precision of this gold standard. The main parameter which influences the statistical uncertainty in the absorbed dose, as calculated by the MC simulation, was the number of simulated decays of the respective radionuclide. Due to limited memory available during the simulations, only M ≈ 10 7 decays could be considered in a single run. In order to increase the total number of decays, and thus decrease the statistical uncertainty of the dose values, we simulated a total number of N runs, which were computed in parallel on the cluster. Random number seeds were strictly different for individual runs, in order to ensure statistically independent results. The standard error SE can be calculated based on equation (3).
where N is the total number of runs, M the total number of decays in a single run, d m,n the dose deposited in a given voxel by decay m during run n. A derivation of equation (3) is provided in the appendix A.2.1. Since tracking sum and sum of squares of d m,n has comparably high computational demands, we carried out the error calculations for one patient only, which we considered representative for our patient population. The value of the SE and the relative SE for the full MC-simulation of this patient is provided in the Results section.

Dose-voxel-kernel methods
For most common therapeutic radionuclides used in Nuclear Medicine, the dominant effect of energy deposition by electrons is due to collisional energy-loss, as can be derived from e.g. the data provided in the EStar database (NIST 0000). The energy transfer along the path of the particle can in principle be computed by the Berger-Seltzer formula (Seltzer and Berger 1982). But if the energy deposited in each voxel of an image should be determined, one would have to sum up all electrons which move inside an individual voxel and integrate over its respective electron paths. Since the velocity distribution of the electrons varies spatially, not only depending on the radionuclide, but also on the traversed path, analytical calculations of deposited dose are impractical in case of inhomogeneous media and off-center voxels. Thus, Monte-Carlo simulations provide the only practical means for estimating deposited energy.
As an alternative to full Monte-Carlo simulations of the particle transport inside a human body, energy dose can be estimated voxel-wise by convolution of the decay map A(r S ) with a dose-voxel-kernel (DVK) d VOI (r T ← r S ). The latter effectively contains the average dose per voxel and decay, which is deposited in a cube of voxels (=VOI) for any nuclear disintegration occurring at a random position in the central-voxel of the cube. Given a certain radionuclide ( 177 Lu in our case), such a DVK is specific for the size and tissue composition of the individual voxels of the cube and can be obtained by Monte-Carlo simulations. Since tissue-type and density vary continuously throughout a patient's body, one would have to simulate an infinite number of different DPKs at the expense of computation time. Alternatively, multiple techniques exist for generating, for heterogeneous tissues, discrete approximations of DVKs with finite voxel size from a limited number of basis DVKs. The latter have been computed in this study from MC simulations for 200 different tissue types spanning the range of mass densities occurring inside the human body. These basis DVKs have then been used to estimate tissuespecific DVKs by proper re-scaling. The kernel scaling methods used in this comparative study are listed below and explained in more detail below.
• Constant Kernel: A homogeneous soft-tissue kernel is used, which is constant for the entire image slice of the patient. This procedure neglects tissue inhomogeneities in the region covered by the DVK. • Center-Voxel Kernel: The DVK is scaled according to the mass density of the voxel of the density map, which corresponds to the center voxel of the DVK. This takes into account variations of mass density of the center voxel, but neglects inhomogeneities in all 9 3 − 1 off-center voxels covered by the DVK. • Density Scaled Kernel: Each voxel of the DVK is scaled according to the individual mass density of the respective voxel in the density map. This takes into account tissue inhomogeneities to some extent, but neglects the influence of the traversed path of the energy-transmitting radiation.
• Percentage Scaled Kernel: Each voxel of the DVK is scaled according to properties of the respective voxel itself and the path from the center voxel in the density map. This takes into account tissue inhomogeneities and the traversed path.
Considering density and percentage scaling, see (Götz et al 2019b) for further explanations. An advantage of these methods is that only a very limited number (here 200) of DVKs needs to be known beforehand. We calculated these DVKs from MC simulations, as explained below. Alternatively, one could have taken appropriate DVKs listed in literature (Grassi et al 2015).

Calculation of basis DVKs
For calculating the DVKs, a Monte-Carlo simulation of 2 · 10 8 177 Lu decays of a Lu nuclide placed at random positions inside the center voxel of a 9 × 9 × 9 cube of voxels (each (4.8 mm) 3 large) of certain tissue-type and density was simulated and the deposited energy scored. This was carried out using the GAMOS Monte-Carlo toolkit on a high-performance compute cluster.

Statistical uncertainty of DVK-dosimetry
The accuracy and precision of DVK dosimetry is not only influenced by the method for adapting the DVK to heterogeneous tissue, but also by the statistical uncertainty and size of the DVK. Note that with decreasing voxel size the poor counting statistics puts any serious dose estimation in jeopardy. The standard error (SE) for individual voxels of the DVK were calculated as described before in section Statistical Uncertainty, using equation appendix A.2.1 with M = 10 7 and N = 20. The resulting relative standard errors of the estimated DVKs are provided in the Results section.

DVK-size
In order to determine the effect of the finite DVK size, an additional simulation of a much larger kernel was carried out, but because of the substantial computational load for one patient only. First, a full MC simulation was performed, as described in 2.4, except that the density and tissue-class of voxels inside the patient were set to correspond to soft tissue, i.e. ρ = 1.04 g cm −3 and t = ST, respectively. This resulting soft tissue dosedistribution was used as reference. Second, one large (81 × 81 × 81 voxels) soft tissue DVK with ρ = 1.04 g cm −3 and t = ST was calculated, as described in 2.5.1. Third, multiple smaller-sized DVKs were subsequently generated by truncating the large DVK. Finally, dose distributions were obtained by convolving the individual activity map with the variously sized DVKs. The deviations of these dose distributions from the reference dose distributions, obtained with full MC simulations, were determined and are provided in the Results section 3.2. However note that for the subsequent analysis of adapting DVKs to different densities and tissue types, the kernel size was kept constant at 9 × 9 × 9 voxels, each with a voxel size of 4.8 × 4.8 × 4.8 mm 3 , which represents a common size in DVK dosimetry (Dieudonné et al 2013, Fernández et al 2013.

Constant scaling
In this method, the DVK used for convolution with the decay map was taken from MC simulations with soft tissue and kept constant for the entire SPECT field-of-view, as also described in Dieudonné et al (2010). All voxels of the DVK thus were assumed to consist of soft tissue (t = ST) as defined in table 4 at a constant density of ρ = 1.04 g cm 3 , as seen in equation (4).
We will refer to this method by the term 'constant'-scaling in the following.

Center-voxel scaling
In this method, all voxels of the DVK were assumed to have the same tissue type as the center voxel of the cube, and the corresponding basis DVK was used to perform the convolution operation during dose estimation. More specifically, proper DVKs were estimated with MC simulations in the following way.
• The mass-density of each voxel of the patient was taken from its density map and discretized to the closest value of ρ which fulfilled ρ = r · 0.01 g · cm −3 with r ∈ {1, 2, . . . , 200}. • For each of these R = 200 discrete densities ρ r , a dedicated homogeneous basis DVK d r,MC (r T ← r S ) was computed with a MC simulation employing 2 · 10 8 decays and assuming that all voxels of the kernel possess the same discrete density. Thus 200 basis DVKs resulted, which could be employed in dose estimations of heterogeneous tissues in the human body. • The center-voxel scaling method now determined from the density map for every VOI and related kernel size the corresponding mass density and tissue type of the center voxel of the VOI. The related homogeneous basis DVK has then been employed to perform the convolution operations needed for the dose estimation of the VOI.
The advantage of this method is that basis kernels need to be determined only once on the outset, while during dose estimation no kernel re-scaling, and its concomitant computation time, is necessary anymore. We will refer to this scaling method simply by using the term 'center-voxel'-scaling in the following.

Density scaling
The energy dose D(r) in a voxel centered at location r is defined as energy imparted to matter per unit mass by ionizing radiation: where dE(r) is the energy deposited in a voxel of volume dV centered at spatial target location r, (r) represents the corresponding deposited energy density, m the mass of the voxel, ρ m (r) the spatially varying mass density of the voxel and V the voxel volume. Both, deposited energy and density depend heavily on tissue specific parameters. In consequence, also the individual voxels of a DVK vary, since the voxels represent the average deposited energydose transported by radiation following any nuclear disintegration in the central voxel. A simplification for calculating DVKs for an inhomogeneous medium was proposed in Dieudonné et al (2013). It was assumed that the deposited energy density (r) was invariant to tissue type t and only dependent on the position r. For example, it can be set equal to the corresponding soft-tissue values, i.e.
which leads to When applied to DVKs, this leads to equation (8) for scaling individual voxels of the DVK according to tissue type t and related voxel density ρ(i, j, k).
Here d ST (i, j, k) represents a DVK which is computed in a MC simulation with all voxels containing soft tissue with spatially invariant mass density ρ ST (i, j, k) = 1.04 g cm −3 . Note that this DVK scaling method only needs computation time (in the ms range) to re-scale the soft tissue kernel according to the local density. We will refer to this scaling method by the term 'density'-scaling in the following. As previously noted, this scaling method assumes that the deposited energy density remains constant, i.e. t ≈ ST , and thus is largely independent from the tissue composition of the voxel. It effectively neglects that the deposited energy not only depends on properties of the respective target voxel, but on the properties of all voxels which have been traversed by the radiation from source to target voxel.

Percentage scaling
Rather than taking into account different tissue densities for DVK calculations via a scaling method as discussed above, in this work a new method is proposed to account for tissue heterogeneity following ideas first expressed in Moghadam et al (2016). The method works as follows.
• Consider a cube with 9 3 voxels, whereby the voxel size plays only a minor role, and assume a radioactive isotope is placed in the center voxel located at r c = a · (ijk) T ∈ VOI in a lattice coordinate system with lattice constant a = 4.8 mm. Assume that the cube consists of a homogeneous tissue r corresponding to one of the R = 200 discrete tissue types, for which basis DVKs have been computed with MC simulations. • Next consider a voxel at location r = a · (cde) T c, d, e 9 belonging to the homogeneous DVK d r (VOI) and compute the voxel-wise ratio δ r := d r (r c )/d r (r). • Next determine the tissue type at any location r T inside a target VOI, e.g. an entire human body or a specific organ, and choose the related basis DVK from the set of R = 200 basis DVKs. • Finally, the homogeneous basis DVKs d r (i,j ,k) will be used to computed new DVKs for heterogeneous tissue densities. To do so, the relative change δ t of dose per decay in the homogeneous basis DVK is used to re-scale the dose per decay in any voxel of the considered heterogeneous DVK with tissue type t according to whereby the following shifts along the grid coordinates are considered: Note that the dose per decay in every voxel of the VOI now depends on the dose deposited in any neighboring voxel closer to the center voxel along a path that directly connects both voxels. Thereby the relative decrease of dose per decay δ r is considered an invariant property of any DVK and is used to scale the dose in all voxels of the heterogeneous kernel by the same relative amount. Note also that the computation of new kernels from the set of basis kernels is computationally expensive, however. We will refer to this new scaling method by the term 'percentage'-scaling in the following.

Evaluated metrics
As mentioned in Data and Methods, since we did not integrate the real TAC for each voxel but used one single SPECT image of each patient, which we took as representative for the relative spatial distribution of the number of decays, we could only calculate differences in absorbed dose between the individual methods, but not absolute absorbed dose values. We motivate this largely be the fact that mean absorbed dose D(r T i , T D ) (Gy) to target voxel r T i from radiation emitted from source locations r S j , j ∈ S over a period T D is given by which under mild restrictions simplifies to where ⊗ indicates a convolution operation in real space. The number of nuclear disintegrations per unit time in the source tissue a(r S j , t) can be obtained experimentally from quantitative SPECT imaging. If integrated over time, the total number of nuclear disintegrations at voxel loaction r S j results, the so-called time-integrated activity (TIA). The quantity d(r T i ← r S j , t) describes a kernel summarizing the energy transfer from the source location to the target location. In most practical applications it can be considered a time-independent factor. Since integration is a linear operation, the spatially varying TIA, hence the estimated dose, is proportional to the local activity at voxel level spatial resolution.
To characterize the dissimilarities between two dose distributions, the following metrics were calculated.
• Relative organ-wise deviation ∆ C (Organ) for DVK-method C with the dose value D C,n and related dose kernel d C,n of voxel n, calculated using DVK-method C, and D MC,n the dose value of the same voxel computed by the full Monte-Carlo simulation. N Organ is the total number of voxels in the organ and C ∈ {constant, center − voxel, density, percentage}. Note that if all dose values are normalized by the total dose of the gold standard, the percentage organ-wise deviation becomes where D X,n = D X,n / NOrgan n=1 D X,n . Note that positive deviations indicate an overestimation of absorbed dose by the DVK method, while negative deviations, accordingly, indicate an underestimation. The latter is to be expected because of the limited size of the kernels used.
• Relative per tissue-type deviation ∆ C (Tissue) for DVK-method C. It is defined in analogy to the organ-wise deviation ∆ C (Organ).
Note that dose values are non-negative by definition, hence value compensations cannot occur. Results were tested for significant differences using a two-sided Wilcoxon signed-rank test. Levels of significance were set to p < 0.05 (indicated by one asterisk, * ) and p < 0.001 (indicated by two asterisks, * * ).

Organ delineation
In order to calculate organ-wise deviations, the respective organs had to be segmented. A medical expert manually delineated the organs liver, both kidneys, spleen, and tumors using the software ITK-snap (Yushkevich et al 2006). This was carried out on the original CT data sets. The resulting segmentation was subsequently downsampled to match the voxel-size of the dose maps. In figure 3, the segmentation for the same patient, already used previously for illustrative purposes, is shown.

MC reference data
For one representative patient, the standard error SE of the MC simulations was determined as described in 2.4. Given a voxel inside the body, the estimated values of the absorbed dose at Cartesian grid points inside this voxel result in a histogram representing a uniform distribution. The mean absorbed energy dose of such voxels is D = (0.72 ± 0.41) × 10 −5 (Gy/decay). If we take into account that N = 1200 samples were drawn to obtain this mean value, a relative standard error of rSE = 1.67% results. Altogether, the statistical error of the MC simulations was found to be rSE 2 % for voxels belonging to the evaluated organs. The largest deviation from the mean absorbed dose is obtained in air surrounding the patient with a maximal relative standard error of rSE = 4.50 %. Figure 4 exhibits the deposited dose per-voxel and the related rSE estimate at one representative slice position.

MC-based basis DVK
The rSE of the voxel values of a typical kernel (t = ST, ρ = 1.04), when computed with MC simulations employing a total number of 2 · 10 8 simulated decays, amounts to ∼0.48%. The error is smaller for voxels near the center of the DVK, compared to more distant ones, as can be seen in figure 5 for a Kernel with t = ST, ρ = 1.04.  Phys. Med. Biol. 64 (2019) 245011 (20pp) Having obtained proper DVKs from MC simulations, the error of dose estimates resulting from the DVK convolution with the decay map needs to be computed in a next step. Let d t (i,j ,k) denote the DVKs for a tissue cube V of size 9 × 9 × 9 = 729 voxels. The absolute error for each voxel is designated as e(i, j, k), and the timeintegrated radioactivity distribution A is assumed to be spatially uniform. The relative error for the calculated dose value for one voxel D(r), centered at r, can be estimated through the following equation: In figure 3, the occurrence of relative standard errors per voxel, resulting from DVK convolutions and estimated for all R = 200 computed basic kernels, is illustrated in form of a histogram. The mean error is 0.0113 %, which is very small because of the high number of decays which were simulated. The rSE is of the same size as the rSE resulting from full MC simulations. Consequently, the dose distribution can be estimated with an accuracy that is only limited by the accuracy of the TIA (see (Götz et al 2019a) for further details).

Impact of DVK size
With increasing size of the DVK used for convolving the time-integrated activity distribution, the deviation relative to values from the full MC simulations decrease, as can be seen in table 1.

Average deviation per organ
In table 2, the mean organ-wise deviations of DVK-dosimetry relative to gold standard MC dosimetry are presented in detail. Also relative deviations of every DVK-based dose estimate from the corresponding MC gold standard are illustrated in figure 6. Altogether, the DVK methods underestimate mean organ-specific energy dose by about ∼7%-10% on average. Note that from the range of observed values it becomes obvious that occasionally absorbed dose even becomes overestimated by the tested methods. The entries to table 2 in fact    represent a grand average over the cohort of 26 patients. When analyzing the differences for all organs together, using a Wilcoxon rank-sum test, it was found that there are no significant differences between the DVK scaling methods when estimating mean dose organ-wise.

Average deviation per tissue-type
In table 3, grand average deviations of DVK dosimetry relative to full MC dosimetry are presented, resolved with respect to tissue type and DVK method. Two patients were excluded from the analysis for class 9 (cortical bone), since large metallic implants were present in the CT images, which led to extended metal artifacts and prevented a meaningful determination of density information. Obviously mean percentage deviations vary strongly with tissue type and less so with DVK method. Again all results consistently reveal an underestimation of mean dose, but the range of dose estimates is considerable and frequent overestimation is observed as well. Despite the large patient cohort also sizable fluctuations around the mean are observed, whereby SDs are smallest for low density tissue and largest for high density tissue. However, the number of voxels, belonging to each tissue class, differs. Therefore the percentage ratios of voxels from one tissue class relative to the total number of voxels within the body are collected in table 4. Also the percentage doses of the voxels from one tissue class are shown.
When analysing these results for the tissue classes lung, adipose, muscle, and cortical bone for significant differences, employing a Wilcoxon rank sum test, the following ordering results (remember that ** indicates a significance level of p = 0.001 and n.s. indicates 'not significant'): These results are also illustrated as box plots in figure 7. The graphs clearly reveal relatively narrow distributions of deviations for low and high density tissues, while tissue with medium densities show rather broad distributions. It is also remarkable that the DVK method 'density' consistently overestimates dose, indicating that a simple linear density scaling deems inappropriate for high density tissues.

Discussion
The presented investigation studied the use of several methods for adapting a set of 200 basis DVKs to tissue inhomogeneities. The homogeneous basis DVKs were established with MC simulations for 200 typical tissue densities. Their adaptation to tissue heterogeneities was meant to improve dosimetry based on DVKs and thus increase dosimetry accuracy without the need to perform elaborate and time-consuming MC simulations for each voxel size, density, tissue-type, and radionuclide in clinical use. Our study encompasses 26 177 Lu therapy patients from clinical routine, for whom the results from DVK dosimetry and full MC dosimetry were compared for several organs and tissue types. The MC simulations were performed employing the GEANT4/GAMOS toolkit, which subsequently acted as gold standard with respect to the absorbed energy dose in each voxel. As full MC simulations are impractical in daily clinical use, the computational load of the different DVK methods and their use in tissues of different densities is of interest. We provide in table 5 a summary of typical computation times needed.

Statistical uncertainty of MC simulation results
We found that the statistical uncertainty, as measured using the rSE of dose values of individual voxels estimated with full MC simulations, was consistently lower than about 2%, when simulating N = 1.2 · 10 10 decays. Our estimate for the precision of dose calculations by DVK-convolution, with homogeneous basis DVKs estimated via MC simulations as well, is in a similar range. Altogether, the achieved precision is similar to comparable studies from literature, e.g. Moghadam et al (2016) state a precision below 2%, but did not disclose how this had been determined. In the study of Dieudonné et al (2013), the authors provide a precision of 4.7% for individual voxels, but also did not describe the method for determining this number. Overall, the accuracy of our simulations should be sufficient for enabling a meaningful comparison of different DVK-scaling methods to results of full MC simulations encompassing the entire body of individual patients.

Influence of size of dose-voxel-kernels
In our evaluation of the influence of DVK size, it was shown for dose determination in soft-tissue that the DVK-based methods systematically underestimate deposited dose, compared to values obtained from full MC simulations. Most likely, this is due to an insufficient sampling of all radiation-matter interactions due to the limited size of a DVK, whereas radiation transport in the MC simulations is only limited by the comparably larger size of the simulation volume, i.e. the entire patient. Consequently, dose underestimation should decrease with increasing size of the DVK. This observation is also corroborated by results reported in literature. For instance, Moghadam et al (2016) used a cube of 5 3 voxels with a voxel size of (4 mm) 3 für soft tissue und bone, while Botta et al (2011) used a DPK of (10 cm) 3 in order to further minimize this systematic error. Table 4. Prototypical materials ordered by their characteristic densities. The material Soft-Tissue, is defined by the atomic composition and aggregated density range of materials marked by * . In the fourth column the percentage ratios of voxels from one tissue class referred to the total number of voxels within the body are illustrated. In the fifth column the percentage dose in each tissue class is shown.

Deviation of DVK-dosimetry in organs
When analyzing the results of organ-wise dose deviations between the considered DVK methods and full MC simulations, it can be seen that, on average, the presented DVK scaling methods obtain similar results for the most critical organs kidneys and spleen. In a previous study (Götz et al 2019b), the total error of an S-value, estimated organ-wise based on the standard MIRD method, i.e. using the software OLINDA, and employing a standard body phantom, relative to a full MC simulation was determined for the same patient cohort. For entire kidneys, the grand average total error of the dose estimate based on S-values amounted to 16.7% ± 12.8%. However, considering the dose estimates obtained in this study, employing a soft-tissue kernel and its adaptation via density scaling, significantly smaller deviations to the gold standard could be achieved, yielding rSE of 7%-10%. It is noteworthy that a relatively high inter-patient variability was observed in these estimates of percentage deviations of absorbed dose between the gold standard and DVK-based methods. When grand averages were considered, this resulted in standard deviations of a similar size as the mean deviation itself (see table 2) and indicated that also frequent overestimation of dose occurs. Generally, an underestimation of dose is seen for all DVK methods when applied to soft tissues. On average, the dose is underestimated by about ∼5%-8%. This is similar to results from other groups, which compared DVK-and MC simulations-dosimetry. For example, Dieudonne et al report an average underestimation of between 2.2% and 5.6% for organs liver, spleen, and renal pelvis/ medulla. Unfortunately, their results were based on one 177 Lu patient only. In a simulation study based on the Zubal phantom, Moghadam et al (2016) report an average deviation of 5.2% between MC simulation and DVK dosimetry for tumorous and normal tissue. The trend of underestimation, seen in our study and by others, could most likely be explained by the finite size of DVKs, as has been shown in our experiment using different kernel sizes (section 3.2).

Figure 7.
Relative deviations between gold standard and four different methods to addapt DVKs to tissue density. The box plots exhibit results for the tissue classes: lung tissue, adipose tissue, muscle tissue and cortical bone.

Deviation of DVK-dosimetry per tissue-type
Since other tissue classes besides previously discussed soft tissue organs are of interest for 177 Lu dosimetry, e.g. in patients with lung, lymph node, or bone metastases, we also evaluated the deviations between dose estimates from DVK methods and full MC simulations for other tissue types, namely lung, adipose, muscle and cortical bone tissue. In contrast to whole organ dosimetry, where all DVK methods showed similar mean deviations with no statistically significant differences, here decent differences could be observed. At first sight it is obvious that the larger the differences in tissue densities are compared to soft tissue, the larger are the deviations of DVK-based estimates from the gold standard. For instance, for the lower-density tissue classes, such as lung and adipose tissue, deviations are pretty large but well-defined with small SDs and when comparing DVK methods, center-and percentage scaling achieved the smallest deviations. Also all estimates for this tissue type yield consistently negative deviations from the gold standard, thus underestimate absorbed dose. The reason for this observation might be twofold. On the one hand, low mass density offers only a small number of interaction partners to emitted radiation, hence only a small amount of energy can be dissipated to matter. On the other hand, dependence on density might not be linear anymore. For muscle tissue, which is slightly denser than soft tissue, again percentage scaling achieves the best results, but the deviations obtained by the other methods only differ by a small margin. With increasing tissue density, mean deviations again grow large. However, with high tissue densities our statistical analysis reveals that the distributions of observed deviations become large also indicating frequent overestimations of absorbed dose in such tissues. Again a non-linear dependence on density might explain the bad performance of the DVK methods investigated. Surprisingly also, the rather simple constant DVK method, although not adapting to higher densities, achieved the best results in cortical bone, whereas density scaling exelled in case of compact bone tissue. According to performed hypothesis testing, these results are highly significant statistically. In summary, it has to be mentioned that the researchers which introduced these DVK adaption methods (except the percentage scaling) in literature, usually evaluated their algorithms with a very small number of patients (e.g. Moghadam et al (2016) one artificial phantom, Dieudonné et al (2013) one patient) and mainly focused on the error in soft-tissue organs. Generally, errors increase as the tissue properties and densities differ from soft-tissue. In addition, inter-patient variability was again pronounced, as can be seen in table 3.
Altogether, when looking at the results for different types of tissues, no single best method exists for all situations and conditions. In order to reduce the error of DVK-based dosimetry, one has to carefully choose the best method according to the tissue and structure of primary dosimetric interest. The results presented in this study may provide guidelines for achieving optimal dose predictions.

Limitations
Although we aimed at thoroughly addressing all relevant issues, some limitations of our study still remain. First of all, the determination of uncertainty in the full Monte-Carlo simulation was performed for one patient only, due to limited computational resources. Although we considered the patient representative, we cannot be absolutely certain that this assessment is justified. However, the results from other groups indicate that our conclusion regarding this point is plausible. Similarly, the influence of the kernel size for DVK dosimetry was evaluated for one patient and one tissue type only (soft-tissue). Nevertheless, these results should be representative for dosimetry of most relevant organs in 177 Lu dosimetry. As DVK-truncation could have a different effect for other tissue types and densities, future studies specifically addressing this point are needed.

Conclusion
We compared the impact of DVK estimation methods to corresponding full MC simulations while estimating absorbed dose in 177 Lu -based radiotherapy. Dose estimation was based on convolution of decay maps with various DVKs, which were adapted to tissue heterogeneities and result were compared to full Monte-Carlo simulations. Results were obtained on a large patient cohort of 26 subjects to assure sufficient statistics. For dosimetry of soft tissues like the liver, kidneys, and spleen, the DVK-method 'constant', employing a soft tissue kernel, achieved results closest to the full Monte-Carlo simulation. Though all other DVK methods achieved similar results, as there were no significant statistical differences between methods, the fast and easy to apply soft tissue DVK method is recommended in such settings. Also all DVK methods consistently underestimated mean absorbed organ dose with respect to full-MC simulations by an amount of 5%-8 %. It was shown that for softtissues this can, at least partially, be attributed to the insufficient kernel-size typically used in DVK-dosimetry. For other tissue types, no single optimal method exists. Generally, for lower-than-soft-tissue densities, percentage and center scaling achieve better results. Surprisingly enough, for a highly dense material as for instance cortical bone, the simple soft-tissue kernel method exhibits the lowest deviation with respect to the gold standard. No obvious explanation for this astonishing observation is yet available to us. Remarkably also, a considerable amount of inter-patient variability of estimated dose has been observed. This biological variability effectively underlines that data of a sufficient number of patients need to be used when evaluating new algorithms in this field.
• A grand average deposited energy dose which, in abuse of notation, could be called a population mean