The accuracy of absorbed dose estimates in tumours determined by Quantitative SPECT: A Monte Carlo study

Abstract

Background. Dosimetry in radionuclide therapy estimates delivered absorbed doses to tumours and ensures that absorbed dose levels to normal organs are below tolerance levels. One procedure is to determine time-activity curves in volumes-of-interests from which the absorbed dose is estimated using SPECT with appropriate corrections for attenuation, scatter and collimator response. From corrected SPECT images the absorbed energy can be calculated by (a) assuming kinetic energy deposited in the same voxel where particles were emitted, (b) convolve with point-dose kernels or (c) use full Monte Carlo (MC) methods. A question arises which dosimetry method is optimal given the limitations in reconstruction- and quantification procedures. Methods. Dosimetry methods (a) and (c) were evaluated by comparing dose-rate volume histograms (DrVHs) from simulated SPECT of ¹¹¹In, ¹⁷⁷Lu, ¹³¹I and bremsstrahlung from ⁹⁰Y to match true dose rate images. The study used a voxel-based phantom with different tumours in the liver. SPECT reconstruction was made using an iterative OSEM method and MC dosimetry was performed using a charged-particle EGS4 program that also was used to determined true absorbed dose rate distributions for the same phantom geometry but without camera limitations. Results. The DrVHs obtained from SPECT differed from true DrVH mainly due to limited spatial resolution. MC dosimetry had a marginal effect because the SPECT spatial resolution is in the same order as the energy distribution caused by the electron track ranges. For ¹³¹I, full MC dosimetry made a difference due to the additional contribution from high-energy photons. SPECT-based DrVHs differ significantly from true DrVHs unless the tumours are considerable larger than the spatial resolution. Conclusion. It is important to understand limitations in quantitative SPECT images and the reasons for apparent heterogeneities since these have an impact on dose-volume histograms. A MC-based dosimetry calculation from SPECT images is not always warranted.

In radionuclide therapy (RNT) the goal is to deliver the highest possible absorbed dose to tumours while ensuring that absorbed dose levels to normal- or risk-organs are kept below tolerance levels so that side effects are minimised [1]. For high-dose treatments, such as myeloablative radioimmunotherapy including stem cell support used for treatment of B-cell lymphoma [2], or radiopeptide therapy used for treatment of neuroendocrine tumours [3], the absorbed doses to normal organs may indeed reach tolerance levels. At our institution we are currently performing a dose escalation study using ⁹⁰Y Zevalin^TM where the amount of activity administered for the therapy is calculated based on the predicted absorbed dose to the first dose limiting organ, based on a pre-therapy ¹¹¹In Zevalin^TM administration [4]. Within this study, we are also performing dosimetry during the therapy based on ⁹⁰Y bremsstrahlung imaging [4,5], to confirm the predicted absorbed dose distribution. Concurrently, we are also performing dosimetry for radiopeptide therapy using ¹⁷⁷Lu-Dotatate [6]. For both ⁹⁰Y Zevalin^TM and ¹⁷⁷Lu-Dotatate, the data concerning tumours reported in the literature are few. This is related to the difficulty in determining the absorbed doses to small structures, given the limited spatial resolution of the imaging system.

Absorbed dose determination in RNT relies on scintillation camera imaging at multiple times after administration to obtain a time-activity curve (TAC) within the volume-of-interest (VOI). The scintillation camera imaging can either be in two-dimensional (2D) planar imaging mode or by three-dimensional (3D) SPECT tomographic imaging mode. The 3D SPECT methodology is a potentially better procedure for dosimetry since it involves more consistent corrections for non-homogeneous attenuation, scatter and collimator response [7]. Moreover, using contemporary hybrid SPECT/computed tomography (CT) systems the mass of the voxel or tissue of interest can be determined from the CT study. From quantitatively accurate SPECT images of the activity distribution, the absorbed energy can be calculated by principally three different ways. Firstly, the kinetic energy released by charged particles can be assumed to be locally absorbed with in the same voxel as the decay. This approach is essentially a rescaling of the voxel activity value and is thus fast. Secondly, a convolution procedure can be applied with an appropriate 3D point-dose kernel that describes the spatial distribution of released energy relative to a point-source [8,9]. In practice, this method is applied as a filtering procedure in the frequency domain with an invariant kernel regarding the spatial shape or interaction properties of surrounding tissues. Thirdly, an explicit 3D Monte Carlo simulation of the energy transport of all charged particles and photons can be performed using a set of patient-specific SPECT images that describe the activity, together with a registered set of CT images describing the density distribution [10,11]. Theoretically, the last method provides the most accurate absorbed dose estimate since it also takes into account differences in interaction probabilities between different tissues. From the absorbed energy distribution the absorbed dose is calculated by determining the voxel-masses obtained from the CT images of the patient. Within a VOI, a dose-volume histogram (DVH) can finally be determined that describes either the frequency or the integral frequency of the voxel-based absorbed dose values.

It must, however, be remembered that a reconstructed activity image is degraded by the limited spatial resolution of the SPECT imaging system and hence the image will not show the real activity distribution even if state-of-the-art correction methods are applied [12]. A natural question therefore arises on how well the absorbed dose distribution estimated from the SPECT images really reflects the true absorbed dose distribution within a VOI and to what extent the distribution is affected by the image reconstruction and other related quantification procedures. The aim of this work is to investigate the dosimetric inaccuracies obtained for different tumour sizes located in the liver. The investigation is performed using Monte Carlo simulated SPECT images in a realistic patient-like geometry. The first and third of the above-mentioned methods of determining the absorbed dose rate are examined in this work, and are applied to images of ¹¹¹In, ¹⁷⁷Lu, ¹³¹I and ⁹⁰Y, being of interest within RNT.

Material and methods

Two Monte Carlo programs were used in this study. The first program was the SIMIND Monte Carlo code [13] that simulates realistic scintillation camera imaging of an activity distribution in a patient-like computer voxel phantom. Realistic SPECT projections were simulated to be reconstructed into tomographic 3D images. The second code was a program [14] for the calculation of the absorbed dose distribution, using a set of 3D activity images and registered density images as input. This program is based on the EGS4 transport code [15] for charged-particle Monte Carlo simulations. The EGS4 program was used both to calculate the absorbed dose rate from the images obtained from the simulated SPECT projection from SIMIND, but also to calculate the ‘true’ absorbed dose rate distribution using the same computer voxel-phantom and activity distribution as used in the SIMIND SPECT simulation. This simulation was, however, performed from the true activity distributions defined as high resolution images without passing the imaging system.

To mimic a clinical patient-like geometry the XCAT anthropomorphic computer software was used [16] to create realistic patient density and activity distributions. The underlying information in XCAT to obtain voxel-images of the phantom is based on NURB (Non-Uniform Rational B-splines) surfaces which makes the software flexible with the ability to create different types of patient-like phantoms. In this work, three spherical tumours of different diameters were defined in the liver region. The volumes of the tumours 1–3 were 2.2 cm³, 17.9 cm³ and 60.7 cm³, corresponding to diameters of 16.2 mm, 32.4 mm and 48.8 mm, respectively. Simulations were made for two cases; (a) with activity in the tumours only, and (b) with activity also defined in the liver. For case (b) the activity concentration ratio between the tumours and the liver was set to 5:2. The total activity in the phantom was defined so that the activity in each tumour was always 1 MBq. For (b), this gave a total phantom activity of 237 MBq, 29.8 MBq and 9.2 MBq for tumours 1–3, respectively. The image matrix size of the XCAT voxel-based phantom was set to 256³ cubic voxels with a side of 1.62 mm. These images were then read into the SIMIND code with the purpose of simulating realistic SPECT projections. The simulation was carried out so that SPECT projections were generated with a very low noise level. However, it must be stated that since the noise obtained from the Monte Carlo simulation do not have the same characteristics as in a real measurement, the noise levels used in this work do not correspond to the noise obtained in practice for a source of 1 MBq.

SPECT projections were simulated from a photon fluence of ¹¹¹In, ¹⁷⁷Lu, ¹³¹I and from bremsstrahlung photons from ⁹⁰Y. A SPECT scintillation camera with a 2.54 cm thick NaI(Tl) crystal was simulated. For ¹⁷⁷Lu and ¹¹¹In a medium-energy general-purpose collimator was simulated whereas for ⁹⁰Y and ¹³¹I, a high-energy general-purpose collimator was used. The energy resolution of the camera was set to 9.5% (FWHM) at 140 keV. Full decay schemes were used thus including the principal energy used for imaging, and also higher photon energies important for the component of septal penetration [17] in the projection images. The SPECT parameters were 60 projections in a 360° full rotation mode. For ¹¹¹In, ¹⁷⁷Lu, and ¹³¹I the energy window was set to 20% centred on the principal photo-peaks. For the ⁹⁰Y bremsstrahlung spectrum that does not have a principal photon energy, a 60% energy window was used centred around 150 keV. The matrix size was 64 × 64 with a pixel size of 6.5 mm for all radionuclides simulated.

SPECT reconstruction was performed using an OSEM iterative method [18] that includes correction for non-homogeneous attenuation, scatter and collimator response. Image reconstruction was made using six iterations and six subsets for all radionuclides and tumour cases. The attenuation correction method used density images obtained from the original XCAT phantom but interpolated to the same voxel and matrix size as used for the SPECT simulation. No miss-registration between SPECT and CT images was thus introduced. Scatter correction was performed using the ESSE method [18] that utilises pre-calculated scatter kernels for the appropriate radionuclide to model the scatter in the SPECT projections. Included in the reconstruction was also a compensation for the degradation in spatial resolution, mainly caused by the geometrical properties of the collimator. For ¹³¹I and ⁹⁰Y, septal penetration and contribution from photons back-scattered in a volume behind the crystal were also included. The back-scatter volume aims to simulate an equivalent to scatter in light-guide and the PM-tubes and was defined as a slab of Lucite with a thickness of 8 cm.

The SIMIND program only simulates photon interactions and, as a consequence, bremsstrahlung photon emission from electron interactions is not possible to simulate in an explicit way. However, simulation of bremsstrahlung SPECT imaging can be performed by using a pre-calculated energy spectrum obtained from a separate electron Monte Carlo program. It is important that this pre-calculated bremsstrahlung spectrum is as accurate as possible with proper photon abundances. In this work, an ⁹⁰Y source in a small tissue-equivalent sphere was simulated in MCNPX to obtain an bremsstrahlung spectrum [19]. This spectrum was then applied into the SIMIND code for the ⁹⁰Y simulations.

The total activity in each spheres were calculated using VOIs in the images was made from the initially known volumes and locations of the tumours. Thus, no manual segmentation or subjective outlining was made. To partly compensate for the spatial resolution the calculation of the activity was based on VOIs that were increased by a shell of a thickness equal to one voxel (i.e. 6.5 mm).

The absorbed dose rate in each voxel was calculated from the reconstructed quantitative SPECT images by two models; either by assuming that all the kinetic energy released from the emitted electrons was locally absorbed within the voxel (denoted below as “EA-SPECT”) or by performing a full Monte Carlo simulation (denoted “MC-SPECT”) using the EGS4 program. In the first case, the absorbed energy from the β-particles was calculated from the average β-particle energy. The difference between these two methods is then that the MC-SPECT calculation takes into account the distribution of energy due to interactions along the particle tracks and also includes bremsstrahlung photon emissions. As the ¹¹¹In radionuclide is a commonly used tracer isotope for ⁹⁰Y-imaging, this radionuclide was also evaluated but with a difference to the other radionuclides in that the absorbed dose rate was calculated using the decay scheme for ⁹⁰Y, using the ¹¹¹In SPECT activity images as input, mimicking the procedure currently applied in monoclonal antibody imaging.

The true absorbed dose rate distribution, denoted as “MC-TRUE” in the text below, was calculated with the same resolution as the original 256³ matrix XCAT phantom with cubic voxels with voxel sides of 1.62 mm. The Monte Carlo calculation of the radiation transport used a complete decay scheme for all radionuclides including emission of gamma photons and bremsstrahlung photons. For ¹⁷⁷Lu, ¹³¹I and ⁹⁰Y, the β-particles have an average energy of 149 keV, 192 keV and 935 keV and an abundance of 78.6%, 89.4% and 100%, respectively. These energies correspond approximately to a particle range in water of 1.9 mm, 3.6 mm and 11.8 mm, respectively. The cut-off values, below which the particle transport is terminated and remaining kinetic energy is absorbed, were set to 30 keV and 10 keV for electrons and photons, respectively. Absorbed dose rate maps were calculated using the individual voxel-masses. Integral dose-rate-volume histograms (DrVH) were determined from the absorbed dose rate images by calculating the frequency of the voxels having a particular absorbed dose rate value.

Results

Table I shows the true and the obtained activity values from the SPECT images for the tumours and total-body, for all radionuclides simulated. Results are shown for cases (a) and (b) described above, i.e. without activity in the liver, and with activity in the liver with a tumour-to-liver ratio of 5:2. Since the activity in the tumours was set to 1 MBq for all tumours, the total-body activity is different for the three tumour simulations and the results are given in Table I as Total Body 1–3. By comparison of the true and obtained activities it is seen that for total body, the activity recovery is within 15% for all radionuclides for both the case with and without activity in the liver. For the tumours in the case without activity in the liver, the recovery is within 5% for the smallest tumour 1 while for the two larger tumours 2 and 3, there is an underestimation of between 2% and 13% which is more pronounced for ⁹⁰Y.

Table I. Calculated activity normalised to the true activity obtained from quantitative SPECT images for total body and within the VOI for each of tumour 1–3. “True” represent the correct activity.

Radionuclide	Total Body 1	Total Body 2	Total Body 3	Tumour 1	Tumour 2	Tumour 3
No activity in liver
True (MBq)	1.00	1.00	1.00	1.00	1.00	1.00
^90Y	0.97	0.94	0.93	0.95	0.91	0.87
^177Lu	0.98	0.98	0.97	0.97	0.95	0.91
^131I	1.14	1.10	1.10	1.03	0.98	0.94
^111In	0.98	0.97	0.97	0.96	0.94	0.91
Activity in liver
True (MBq)	236.8	29.8	9.2	1.0	1.0	1.0
^90Y	0.89	0.89	0.89	1.84	1.22	0.93
^177Lu	0.98	0.99	0.98	2.16	1.32	0.99
^131I	1.09	1.09	1.09	2.26	1.38	1.01
^111In	0.97	0.98	0.98	2.15	1.31	0.98

For the case with activity in the liver there is an over-estimation for the tumour 1. These results are is most likely related to the spill-over of counts from the liver region into the tumour region for the smaller tumour and a reversed spill-out effect for the larger tumour, due to the partial-volume effect. For tumour 3 with liver activity, the activity recovery is very good with little contribution from spill-over.

Figure 1 shows tomographic images of the absorbed dose rate in a transversal slice through the liver with tumour 2. It is clear that the absorbed dose rate in the MC-TRUE image is very well-defined even for ⁹⁰Y which has the highest β-particle energy. For most of the radionuclides the absorbed energy and related absorbed dose rate is close to the source but for ¹³¹I it can be seen that a relatively large fraction of energy is deposited outside the source region. This is due to the energy deposition from high-energy photons (364 keV, 637 keV and 723 keV) that have a large mean-free path length. The pattern of absorbed dose rate values outside the liver differs somewhat between EA-SPECT and MC-SPECT as seen for instance for ¹⁷⁷Lu and ¹¹¹In where for MC-SPECT there is a distribution of values also outside the liver.

Figure 1. Absorbed dose rate images calculated from a) the assumption of local average energy deposition from charged particles within the voxel where the decay occur (EA-SPECT), b) a full Monte Carlo simulation of the energy deposition using activity distribution from quantitative SPECT images (MC-SPECT) and c) a full Monte Carlo from the true high-resolution activity distribution (MC-TRUE).

Figure 2 show dose-rate volume histograms (DrVH) for the case of no activity present in the liver. For ¹⁷⁷Lu and ¹³¹I, the DrVH for MC-TRUE exhibits a steep decent in absorbed dose rate values in the higher range indicating that the absorbed dose rate per voxel is uniform over most of the volume. For ⁹⁰Y a broader distribution of absorbed dose rate values is seen. This can be explained by the fact that the β-particles emitted by ⁹⁰Y have the highest kinetic energy (2.23 MeV) and can thus escape the VOI with a significant amount of energy. For ¹¹¹In, where the absorbed dose rate calculation in performed for ⁹⁰Y, the distribution is almost identical to that obtained from the ⁹⁰Y quantitative SPECT images. This supports the idea of using ⁹⁰Y bremsstrahlung images as a valid alternative to tracer ¹¹¹In images, even for dosimetry for small structures such as tumours.

Figure 2. Dose-rate Volume Histograms (DrVH) for the case with activity located only in the tumours. Solid lines represent EA-SPECT, short dashed lines represent MC-SPECT and long dashed lines represent MC-TRUE. Note that the range on the x-scale is not equal for the different graphs but remains the same between this figure and Figure 3.

Thus, generally it can be seen that the DrVHs obtained from SPECT are significantly less sharp than the true DrVHs, a result mainly an effect of limited spatial resolution. As the tumour size increase (tumour 1 → 3) the DrVHs determined from SPECT come closer to the true DrVHs but the shapes are still different. It should be remembered that the dimensions of tumour 1 is in the same order as the spatial resolution. Moreover, the maximum absorbed dose rate values as determined from SPECT images are higher than the true values, likely because of noise propagation in the reconstruction method. Generally, explicit MC calculation of the absorbed dose rate from the SPECT images has a marginal effect except for ¹³¹I where the photons contribute significantly.

Figure 3 show similar DrVHs for the case of activity present in the liver. The difference between MC-SPECT and EA-SPECT is generally non-discernible, except for ¹³¹I. Thus explicit Monte Carlo calculation has little effect compared to the assumption of local voxel energy deposition except for ¹³¹I where there is a significant energy deposition from high-energy photons.

Figure 3. Dose-rate Volume Histograms (DrVH) for the case with activity both located in the tumour and in the surrounding liver and with a tumour-to-liver activity concentration of 5:2. Solid lines represent EA-SPECT, dotted lines represent MC-SPECT and dashed lines represent MC-TRUE.

Comparing the results of the true DrVHs of Figures 2 and 3, it is seen that the difference is small except for ¹³¹I. For pure β-emitters, such as ⁹⁰Y, or radionuclides with a small photon fraction in the decay, such as ¹⁷⁷Lu, the gain of using Monte Carlo based dosimetry calculation is thus questionable for structures located in soft-tissue, even for relatively long-ranged β-particles such as those emitted by ⁹⁰Y. For ¹³¹I, the absorbed dose rate contribution from photons in non-neglectable, and a Monte Carlo based calculation is warranted.

Figure 4 display profiles through the centre of the tumour for case (b) with activity in the liver and for the four different radionuclides. From the figure it is evident that there is a large difference between the curve shapes for the small tumour 1, but that there is a better correspondence between the curve shapes for the larger tumours.

Figure 4. Profiles through the tumours for case b) with activity present in the liver are shown for the four radionuclides. Solid lines represent EA-SPECT, dotted lines represent MC-SPECT and dashed lines represent MC-TRUE.

Discussion

An accurate dosimetry is essential both for planning radionuclide therapies but also for follow-up of studies and research concerning dose-response. Today, quantitative SPECT combined with registered anatomical CT images provide data that are useful for calculations of the absorbed energy and the absorbed dose. Absorbed dose calculation from SPECT images can be performed by (a) assuming that the average energy from the charged particles are deposited within the voxel where the decay occurred, (b) a convolution procedure based on point-dose kernels or (c) by performing a full Monte Carlo dosimetry calculation. The calculation time and complexity differs significantly between these methods, (a) being the simplest and fastest method. All methods rely on a high degree of accuracy in the quantitative SPECT images. The present study focuses on the comparison of methods (a) and (c) for spherical tumours of different diameters located in the liver, for radionuclides of interest in RNT. Results show that dosimetry using explicit Monte Carlo calculations is not always warranted, a result that may appear counter-intuitive. The main reason is that the blurring effect caused by the limited spatial resolution of the SPECT imaging system (mainly caused by the collimator properties and the distance to the object) is in the same order of magnitude, or even higher, as compared to the particle ranges of the charged electrons. For ⁹⁰Y, ¹⁷⁷Lu and ¹¹¹In→⁹⁰Y, this is seen in the DrVHs of Figures 2 and 3, where there is almost no difference between the Monte Carlo based calculations and the results from the assumption of local energy deposition. For ¹³¹I, there is a significant contribution of photon energy deposition, which makes a Monte Carlo based calculation warranted. The effect is more pronounced for smaller objects such as the smallest tumour.

The results further show that the determined activity within the three tumours also depends on the surrounding activity distribution (in this study the presence or absence of liver activity). This difference is caused by spill-over and spill-out of counts into and out from the tumour VOI. For the results in Table I we used a simple approach of adding a shell of voxels to the known voxel volume for each tumour. A partial-volume correction (PVC) for this effect can also possibly be made by experimentally determining recovery coefficients for tumour-mimicking spheres of different dimensions within a phantom, but the main problems is that spill-out and spill-over also depends on the tumour shape and on the distribution of the surrounding activity [20]. A more generalised method would be to apply image-based PVC [21,22].

It can be seen from Figure 4 and especially for the graphs showing tumour 3 that a small dip in the profile matching the centre of the tumour. The reasons for this underestimation are related to the collimator-response correction. For high-contrast objects with a sharp boundary-edge a phenomenon called Gibbs effect [23] can introduce ringing artefacts in the collimator-response correction. Despite the fact that the total counts are well preserved in the collimator-response compensation, the count distribution is pushed towards the edges of the tumour VOI and therefore potentially suffers more from the partial volume effects than the smallest tumour. One way of omitting this effect is, of course, to disregard the collimator-response correction but for ¹³¹I and ⁹⁰Y imaging such a correction can also include septal penetration and is therefore very important because emitted high-energy photons significantly contribute to the total counts in the image due to penetration and therefore needs to be compensated for.

The absorbed dose rate distribution heavily depends on the accuracy of the quantitative SPECT images regarding both the spatial resolution and signal to noise ratio. These parameters depend on the choice of reconstruction method and the parameters used for the corrections of attenuation and scatter. It is known that the convergence in an iterative MLEM/OSEM algorithm is faster for large objects than for small objects and that noise may be amplified when using a large number of iterations. This means that the counts distribution may be larger within a VOI due to noise propagation and amplification but the average absorbed dose value this can converge towards a correct mean value. This is, of course, good if the objective is to get the average absorbed dose but if the variation of absorbed doses within a VOI is to be taken into consideration then noise amplification should be considered. For example, the concept of EUD (Equivalent Uniform Dose) has recently been introduced for biological interpretation of tumour absorbed dose distributions, for correlation with treatment response [1]. It is then especially important to keep the limitations in the DVHs in mind, so that heterogeneities in the SPECT images that are really caused by a limited system spatial resolution or noise are not interpreted as variability in radiopharmaceutical uptake related to a biological variation in tumour cell viability. An example of a regularisation method is the MAP iterative reconstruction method [24] where a known constraint, such as a restriction in the local voxel-value variations, are built into the reconstruction as a penalty-factor. These methods, although not commercially available at the moment, may reduce the noise amplification when using a large number of iterations [25]. In this work, all reconstructions have been made with the same number of iterations and subsets and could potentially have been optimised for each tumour/liver combination. Methods for regularisation of the collimator response compensation to avoid effects on the DrVH such as the Gibbs ringing artefacts are also important to further investigate. However, such an investigation is beyond the scope of this paper.

A full Monte Carlo simulation of the absorbed dose seems justified when the particle range is much longer than the spatial resolution of the image from which the absorbed dose calculations are conducted. For most β-particles and electrons located in soft tissue this is, however, not the case when considering SPECT images obtained from clinical SPECT systems. The maximum range of ⁹⁰Y is about 11.8 mm and the mean range of ∼5 mm. This is, however, the total path-length but since the direction of the particle change during the slowing-down process the mean radial distance from the decay (the projected path-length) will be even less. In soft tissue, it is thus only the presence of high-energy photons that make a complete 3D Monte Carlo based absorbed dose calculation warranted for clinical SPECT images. The lung-region may be a region where Monte Carlo calculations can be justified due to the low density. For high-resolution clinical PET and pre-clinical μSPECT/μPET systems where, in the case of μSPECT, the resolution is perhaps less than 1 mm then the particle tracks and related cross-dose will be important.

Acknowledgements

The Swedish Cancer Foundation, Berta Kamprad's Foundation, The Swedish Research Council and Gunnar Nilsson's Foundation.

Declaration of interest: The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.