Quantification of internal dosimetry in PET patients: individualized Monte Carlo vs generic phantom-based calculations.

PURPOSE: The purpose of this work is to calculate individualized dose distributions in patients undergoing 18 F-FDG PET/CT studies through a methodology based on full Monte Carlo (MC) simulations and PET/CT patient images, and to compare such values with those obtained by employing nonindividualized phantom-based methods.
METHODS: We developed a MC-based methodology for individualized internal dose calculations, which relies on CT images (for organ segmentation and dose deposition), PET images (for organ segmentation and distributions of activities), and a biokinetic model (which works with information provided by PET and CT images) to obtain cumulated activities. The software vGATE version 8.1. was employed to carry out the Monte Carlo calculations. We also calculated deposited doses with nonindividualized phantom-based methods (Cristy-Eckerman, Stabin, and ICRP-133).
RESULTS: Median MC-calculated dose/activity values are within 0.01-0.03 mGy/MBq for most organs, with higher doses delivered especially to the bladder wall, major vessels, and brain (medians of 0.058, 0.060, 0.066 mGy/MBq, respectively). Comparison with values obtained with nonindividualized phantom-based methods has shown important differences in many cases (ranging from -80% to + 260%). These differences are significant (p < 0.05) for several organs/tissues, namely, remaining tissues, adrenals, bladder wall, bones, upper large intestine, heart, pancreas, skin, and stomach wall.
CONCLUSIONS: The methodology presented in this work is a viable and useful method to calculate internal dose distributions in patients undergoing medical procedures involving radiopharmaceuticals, individually, with higher accuracy than phantom-based methods, fulfilling the guidelines provided by the European Council directive 2013/59/Euratom.

Introduction

Radiopharmaceuticals are widely used in medicine for diagnosis and therapy. , The European Association of Nuclear Medicine Internal Dosimetry Task Force has recently reported a significant expansion in the use of radiopharmaceuticals. Positron emission tomography (PET) represents the gold standard of functional imaging in many applications, such as oncology, neurology, and cardiology. , , The use of PET is increasing; five studies per 1000 habitants are performed worldwide each year. The most used positron‐emitting radiopharmaceutical in PET is 18F‐fluoro‐2‐deoxy‐D‐glucose (FDG).

Radiation dosimetry of patients undergoing procedures involving radiopharmaceuticals is far from the accuracy achieved in other radiation applications. Absorbed dose is associated with tumor control and toxicity in therapy, and with cancer induction risk in diagnostic procedures. To quantify such risks, it is important to know accurately the absorbed doses. Therefore, there is a need to develop precise, individualized methods for internal dosimetry in this field. This need has been highlighted in the recent European Council directive 2013/59/Euratom (Art. 56). The directive establishes minimum safety requirements for ionizing radiation exposure in the member states of the EU.

Traditionally, internal dosimetry in this field is determined from dosimetric factors calculated in phantoms, the so‐called S‐factors, relating disintegrations in a given organ to doses deposited in another organ. Internal dose distributions can be calculated applying such S‐factors to a distribution of activities in each organ (either measured or calculated with a nonindividualized biokinetic model). The calculation of S‐factors can be nowadays performed through precise Monte Carlo (MC) dose calculations. , However, internal dose distributions obtained with these methods have limitations, as they are not individualized for each patient, phantoms may not be a good representation of patients, and the spatial heterogeneity of activities in organs and tumors is not considered. , , ,

More precise methods to address internal dosimetry are convolution/superposition methods, which rely on an accurate computation of dose deposition kernels through Monte Carlo calculations, and the convolution of these kernels with distributions of activities measured in the patient. Kernels are usually calculated in homogeneous media, which requires scaling kernels with the radiological distance when applying them to (heterogeneous) patients.

One of the most used programs is the Organ Level Internal Dose Assessment for Exponential Modeling (OLINDA/EXM), based on the Dose Assessment Resource (RADAR) formalism. This platform uses a set of generic phantoms to transform (by scaling factors) the biokinetics of a radiopharmaceutical into absorbed doses to organs. OLINDA/EXM allows the inclusion of some experimental biokinetic measurements, in addition to considering the self‐dose to pathological regions such as tumors. There exists a new generation of phantoms that deform the MIRD/ICRP standards of voxelized computational phantoms to patient characteristics, that is, NURBS‐based phantoms.

Monte Carlo methods can also be used to obtain internal dose distributions. They are considered the gold standard for dosimetry calculations, albeit at the cost of large computational times. There are several MC codes that may be employed for this aim, from general‐purpose, multiapplication codes like Geant4, MCNP, or FLUKA, to codes which are mostly focused on specific applications, like EGSnrc (dosimetry in radiotherapy), or GATE , , (medical imaging, www.opengatecollaboration.org). MC methods have been used to obtain internal dose distributions in medical procedures with radionuclides, both in therapeutic and diagnostic applications. , , Some of them evolved from external radiotherapy platforms, such as MINERVA or DPM. The latter, based on Fortran, presents significant speed advantages over other Monte Carlo platforms in voxelized geometries. Others were developed ad hoc using open‐source Monte Carlo collaborations, like VIDA or RAYDOSE, that implement a patient‐specific dosimetry calculation with Geant4. Recently, the platform RAPID has been introduced, a patient‐specific Geant4‐based dose calculator for molecular radiotherapy, which calculates on the CT of the patient, with sources constructed from PET/SPECT images. The platform was validated by calculating S‐factors in phantoms. All theses codes were benchmarked against phantom‐based methods, however, they have not been extensively used on patients or clinical trials.

In this work, we follow a methodology similar to that of RAPID, but focusing on internal dosimetry in patients undergoing FDG‐PET studies. Our methodology performs individualized MC dose calculations on the CT of the patient, with activity sources constructed from the (single) PET study and a biokinetic model. The code GATE was chosen since it allows easy use of voxelized geometries and sources, in addition to including several tools for internal dosimetry. We present results for 14 patients, and compare these results with those obtained from several nonindividualized methods with generic phantoms: Cristy–Eckerman's (CE), Stabin's formalism, based on ICRP‐89, , and reference computational adults from the ICRP‐133. The uncertainty budget of the calculations is thoroughly studied. Individualized MC dosimetry in patients undergoing PET studies has not been fully addressed, and studies focusing on such computations, and the comparison with phantom‐based methods, may be of interest to the community of medical physicists and nuclear medicine physicians. To our best knowledge, this is the first study presenting MC individualized dose calculations in a cohort of FDG‐PET patients, and comparing those results to phantom‐based, nonindividualized methods.

Materials and Methods

PET/CT equipment

Whole‐body PET/CT acquisitions were obtained with a Ingenuity TF PET/CT (Philips Healthcare, Best, Netherlands). The PET system implements a Time‐of‐Flight technology with an 18 cm axial and a 67 cm transaxial Field‐of‐View. The detectors are based on LYSO crystals with a 4 × 4 × 22 mm3 size. The system has a sensitivity of 7400 cps/MBq (center), 11.7% energy resolution, and 495 ps timing resolution. The CT system has 64 slices with 40 mm coverage, and employs solid‐state GOS detectors.

Study and patient cohort

Patients who were undergoing a PET/CT scan were invited to participate in this study. In the first phase of the study, we recruited 14 patients. We restricted the study to patients who had undergone a full‐body PET/CT, to avoid the introduction of bias in our dosimetry calculations due to missing activity. The main characteristics of our cohort are summarized in the Supplementary Materials (section SM2.1 and Table SM1).

PET/CT protocol

Patients fasted for at least 6 h before the injection of FDG. The activity was weighted by body mass index (BMI) with a mean value of 256.2 MBq (range, 166.5–333.0 MBq) per patient. Patients rested 30 mins preadministration in a warm room, and 45 mins postadministration to ensure correct biodistribution. The PET/CT studies lasted around 35 mins and consisted of a CT scan (helical mode with pitch factor of 0.828 and slice thickness of 3 mm, reconstruction matrix of 512 × 512, 1 mm pixel size, 3 mm slice thickness, 120 kV, current modulation by DoseRight, revolution time 0.5 sec, and slice collimation 64 × 0.625 mm) followed by a PET‐scan (3D acquisition mode, 2 mins per bed, and a 144 × 144 reconstruction matrix).

Image segmentation

The tissues/organs that were deemed important either for dose evaluation or because they factor significantly in the biokinetics of FDG (either directly or because they are well perfused) were segmented in the images. All major organs/tissues were included (a full list is provided in Table SM2, Supplementary Materials). Segmentation was performed on the CT images, by using the Eclipse platform (Varian Medical Systems; Palo Alto, CA) and the “Smart Segmentation” algorithm. This algorithm includes an atlas or model‐based approach for automated contouring of patient anatomy. For the segmentation of the bladder, due to changes in the bladder volume between CT and PET acquisitions, we also used the PET images. The procedure was carried out as follows: firstly, bladder contours were manually drawn to include all the activity in the PET image that was thought to arise from the bladder (including, e.g., spill‐out effects of the bladder activity), these contours are used to model emptying (section 2.5); secondly, the bladder was segmented using the 3D active contours Chan–Vese method. The motivation for and uncertainty analysis of the chosen methodology for bladder segmentation are discussed in the Supplementary Materials (section SM2.2).

Biokinetic model

The acquired PET studies were static, meaning that no kinetic information of the biodistribution of FDG was available. We need to include biodistribution models to obtain time‐activity curves in each voxel, which will then be used to obtain cumulated activities (CA). To determine the biodistribution of the FDG in the body, we used a simplification of the model of Hays and Segal. A diagram of the model is presented in Figure SM1 (Supplementary Materials). The model consists of seven compartments representing blood, liver, brain, heart walls, lungs, bladder contents, and the remaining tissues. The set of differential equations describing variations of activities, y, in each compartment as a function of time can be written as:where the k's are the uptake/release parameters from the i‐th compartment, and λ is the decay constant. Solving the resulting ODE system provides the activity curve as a function of time (TAC) for each compartment present in the biokinetic model. Integration of these TACs in the interval [0, ∞) gives a value of cumulated activity for each compartment. The main advantage of this procedure is that it allows one to obtain a good description of the activity distribution without the need for multiple experimental measurements over time.

The model was optimized for each patient by using a simulated annealing (SA) algorithm. Nine kinetic parameters plus five blood fractions were optimized during the fitting. The cost function was defined as the sum of the weighted square differences between experimental and theoretical data:where E and T are the experimental and theoretical data for a patient, respectively, and σ is the uncertainty associated with experimental data.

Due to the lack of dynamic data, only one experimental point was available (at the study time t) for each organ. When performing such an optimization, there is an intrinsic degeneration of the solution, and the obtained activity curves may differ from reported values. To establish a constraint on the shape of the constructed activity curves and to match the reported experimental activity curves, a second series of data was introduced for the lungs, brain, liver, and blood. Relative activity curves for those organs were obtained from published dynamic PET studies , and used to incorporate activity constraints at t'≈20 min. The new data point was scaled from the published curve according to the experimental activities measured in this study at time t,where A is the time‐activity curve taken from literature, A(t’) is the scaled value of A at time 20 min, and A(t is the activity at time t measured by us.

Different weights were assigned during the optimization of the time‐activity curves, depending on the origin of the data: for our own experimentally measured activities, we used σ = 5%; for the second series (constraints at t’ = 20 min extracted from dynamic studies) we considered σ = 30%; and in the special case of blood we used σ = 50%. Regarding our own data, the 5% that we use is a compromise between large organs with high uptake, which may have uncertainties below that number, and small organs that may have larger uncertainties. On the other hand, for constraint points taken from other experiments we consider much larger uncertainties to allow for interpatient variability and to have flexibility in the fit of time‐activity curves. The addition of data taken from the literature narrowed the set of valid solutions which could be obtained from the optimization, caused by the limitations of experimental time‐points (static PET studies).

As part of our biokinetic model, we assume that patients will empty their bladder at a time t. This will eliminate activity in the bladder and will affect the absorbed dose estimates, especially in its vicinity.

Calculation of cumulated activities

In order to calculate maps of cumulated activities per voxel, firstly, PET files were loaded into a Matlab (The Mathworks, Natwick, MA) 3D array, which was multiplied by the voxel volume to obtain the experimental activity in each voxel, A. The next step involved the application of the biokinetic model described in section 2.5. The resulting system of [eq. (1)] was solved and integrated over a sufficiently longtime. The heterogeneities of activities in organs were conserved during the construction of the source map: the cumulated activity in each voxel belonging to a given organ was computed considering the CA of that organ (biokinetic model), and the activity of the voxel in the PET image:where i runs over the number of pixels, N is the total number of pixels for that organ, and A is the activity of the i‐th pixel in the PET image. The CA of the remaining tissues is shared among all the nonexplicitly considered organs in the biokinetic model according to the same equation.

To account for the deposited dose due to activity in the blood, the activity in the blood compartment has to be distributed among the rest of the structures. The blood CA in a voxel i belonging to a given organ, CAblood,i was calculated as the total CAblood multiplied by an organ weight and divided by the number of voxels in the organ to which the voxel belongs. These organ weights, not to be confused with physical weights, depend on the percentage of blood volume that the organ contains. Reference values were obtained from the literature for different patients’ age and sex (section 7.7.2. in Valentin et al and Table 5 in Wayson et al ). They were modified during the optimization for the liver, lungs, heart, brain and remaining tissues’ blood fractions. The source was finally composed of the sum of blood, CAblood,i, and intrinsic contribution, CAi. Figures SM2 and SM3 show a few example slices of an original PET and a CA source, respectively.

For the sake of completeness, we also investigated the dose distribution using nonindividualized cumulated activities obtained from the ICRP‐106 FDG biokinetic model, which is widely used for internal dose calculations.

Monte Carlo platform

The software vGATE version 8.1. was employed to carry out the MC calculations. GATE is a MC code widely used for imaging applications, which now has tools to compute internal dosimetry features. The methodology involves three main steps to perform the individual dosimetry evaluation: preprocessing, MC simulation, and postprocessing. Pre‐ and postprocessing were performed in Matlab.

The individualized geometry is created from CT scans and introduced as DICOM files in GATE. The source file (cumulated activities per voxel) is constructed as described in the previous section and converted to DICOM format to be used by GATE. The simulations were performed in a virtual machine hosted at CESGA (Galician Supercomputing Centre). For this preliminary stage, CESGA staff deployed a virtual machine with 12 cores and 64 GB of memory running in OpenNebula Cloud. The employed physics package was the emstandard_opt3, with 10 keV cutoffs for γ, e‐, and e+. The simulation was launched for 10 9 disintegrations, which were split into two cores per simulation. This number of disintegrations was chosen because preliminary commissioning of the simulation showed that this number is enough to obtain good statistics. Typically, one full simulation lasted around 36 h (but only 6 h if using the virtual machine's full core capacity). After the simulation, the two jobs were merged following the guidelines in Chetty et al. GATE returns a 3D matrix of absorbed doses, from where mean absorbed doses per organ/tissue were calculated by identifying every voxel belonging to a certain organ/tissue. Detailed information on pre‐ and postprocessing is presented in the Supplementary Materials (section SM2.3).

Validation of MC simulations

The validation of the MC simulations was threefold. Firstly, we compared the dose deposition with GATE to that of EGSnrc, which acts as the gold standard to calculate dose deposition in the energy range employed in medicine. This was done in simple phantoms consisting of parallelepiped water phantoms with cavities of different materials, and in patient geometries (CT). Secondly, we used our code to obtain Specific Absorbed Fractions (SAFs) in the ICRP‐133 male phantom and compared the results to tabulated values. Finally, we measured the doses inside and outside of a Jaszczak phantom filled with a 18F solution using extremity dosimeters, , and compared the measured doses to the GATE simulations. Detailed information can be found in the Supplementary Materials (section SM2.4).

Generic phantom‐based methods

We compared the dose values obtained with our MC‐individualized method to the values obtained with generic phantom‐based methods. In particular, we used three models as a reference: Cristy‐Eckerman's (CE), Stabin's (S89), , and the reference computational adults from ICRP‐133.

To estimate doses, it was necessary to compute Dose Fractions (DFs) from SAFs according to the radionuclide decay scheme. The specific method is well described by Snyder et al. After obtaining the DFs, mean dose for a target organ r was finally computed as:where r is the source organ, and à is the total cumulated activity in each source organ. The array of cumulated activities in each organ was obtained from the source matrix described in section 2.6., that is, we did not use tabulated cumulated activities, but individualized patient activities according to our experimental data and biokinetic model.

The addition of patient‐specific information in generic phantom methods was also studied. Tabulated SAFs for different stylized phantom models were scaled by the organ masses of each patient, as described in Stabin et al.

Uncertainty budget

The uncertainties in the computed doses (both MC and phantom based) may be high. They arise from three main sources: i) statistical uncertainties in the MC calculations; ii) the biokinetic model to include activities during the uptake of the radiopharmaceutical and bladder emptying; and iii) segmentation of organs that can change volume in between CT and PET study (in particular the bladder), as well as fusion of the PET/CT studies. We have performed a thorough study of the contributions of these factors to dose uncertainties.

The statistical uncertainties arising from the MC simulation are proportional to 1/√N, where N is the number of simulated particles. We have investigated this contribution to the simulated number of disintegrations.

Due to the lack of dynamic information, the uncertainties that are associated with the biokinetic model may be high, and it is important to study them. Two different tests were performed for this purpose. A first analysis (Test 1) consisted of measuring how the degeneration of the optimization problem [Eq. (2)] affected cumulated activities, CA: a set of 100 optimizations were run for two different patients and the associated CAs were calculated with the kinetic parameters of each optimization. The second test provided an estimation of the impact of uncertainties in the kinetic parameters on the absorbed doses: the kinetic parameters obtained from one optimization were subjected to Gaussian perturbations with a standard deviation equalling 10% of the mean; the effect of such perturbations on the CAs in each organ was analyzed; this experiment was repeated 500 times. Also, we studied how variations in CAs propagate to variations in absorbed doses. This latter study was limited to one patient.

The time of bladder emptying, if not properly monitored, causes further uncertainties in dose distributions. We have studied the effect of different values of t on internal dose distributions. An in‐depth analysis was conducted regarding the effect of different approaches to the segmentation of the urinary bladder (UB), motivated by critical differences in its volume between CT and PET images. Details about this analysis are included in the Supplementary Materials (section SM2.2).

Statistics

Statistical comparisons have been performed by using a paired‐sample t‐test. We have relied on the statistics toolbox of Matlab. Statistical significance has been defined as P < 0.05. Data dispersion among patients is represented by using boxplots: the central ticks show the median of the data; the boxes span from the 25th to the 75th percentiles; and the whiskers extend to the most extreme data points, excluding outliers (beyond 2.7σ), which are shown as individual points.

Results

Validation

A comparison of the dose distributions obtained with GATE and EGSnrc showed similar results, both in phantom and CT patient geometries. Nevertheless, some inconsistencies were found in low‐density regions (e.g., air), where the GATE doses are lower than the EGS doses for photons (see Figures SM5‐SM8, Supplementary Materials). Regarding the calculation of SAFs in the ICRP‐133 male phantom, the GATE results agree within ± 5% with tabulated values for most target‐organ pairs. In some cases, there are important differences for distant organs (e.g., liver–brain) and low‐energy electrons. The results are shown in Figure SM9 (Supplementary Materials). The observed differences are most likely due to differences in the physics packages in each MC software. The comparison between experimental doses and MC calculations in several locations in and around the Jaszczak phantom filled with a 18F solution generally showed a good agreement (Figure SM10). Detailed results and discussion are presented in the Supplementary Materials (section SM3.1).

Statistical uncertainties arising from the MC simulation are in general low due to the large number of disintegrations. In this work, relative statistical uncertainties of the mean dose deposited in each organ range from ~ 0.02% (large organs receiving important doses, e.g., brain), to ~ 5% for the lenses, which receive low doses and contain few voxels. If necessary, this uncertainty can be reduced by increasing the number of disintegrations, at the cost of increasing calculation times.

The best‐fitting parameters of the biokinetic model applied to each patient are shown in the Supplementary Materials (Tables SM3 and SM4). In Table 1 we show the results of the uncertainty analysis of the biokinetic model. In the degeneration analysis, we found that the existence of multiple minima did not have a heavy impact on the variability of CA for the majority of compartments, being usually less than 5%, but reaching 20% for the liver. On the other hand, the introduction of Gaussian perturbations around optimal kinetic parameters leads to CA variations of around 10% for most compartments. It is important to note that these differences refer to CAs: time‐activity curves could be very different among them, yet result in very similar CAs. In Fig. 1 we show representative time‐activity curves obtained with our fitting method: in general, time‐activity curves match the shape of the curves reported in dynamic studies. In Fig. 2 we show the effect of degeneration and parameter perturbation on representative time‐activity curves.

Table I

Relative variations (std/mean) of cumulated activities in different organs due to perturbations of the biokinetic model for two patients (P1 and P2). Test 1 measured the variability due to the degeneration of minimum of the cost function, while Test 2 quantified the effect of Gaussian perturbations added to optimal biokinetic parameters (see text, section 2.10).

		Blood	Remaining Tissues	Liver	Brain	Heart	Lung	Bladder
Test 1	P1	3.5%	1.2%	20.9%	3.1%	4.1%	1.3%	3.1%
	P2	5.2%	1.2%	19.8%	4.0%	7.1%	5.4%	3.8%
Test 2	P1	7.0%	3.9%	11.4%	9.6%	7.1%	8.8%	11.5%
	P2	6.8%	4.5%	9.9%	10.9%	6.9%	6.9%	11.3%

Fig. 1

Time activity curves obtained from the simulated annealing fit for two patients, (a)‐(c) and (d)‐(f), and different organs: liver, (a) and (d); brain, (b) and (e); lungs, (c) and (f). Cross markers represent the experimental curves which were incorporated from the literature and scaled to the known PET activity at the study time for each patient. Triangle markers show the experimental activity value which was directly obtained from the PET study (at t> 40 min) and the one calculated from scaling the published curve with the PET activity (at t = 20 min). Error bars represent the estimated weight allowed during the SA optimization (5% and 30%, respectively). The solid line shows the fit of the biokinetic model to the experimental data. [Color figure can be viewed at wileyonlinelibrary.com]

Fig. 2

Uncertainty analysis for the biokinetic model. Each column represents a compartment, Liver (a, d), Brain (b, e), Lungs (c, f), and each row a performed test (a‐c for Test 1 and d‐f for Test 2). In the first row, solid gray lines show the result for 75 SA optimizations which were run for a particular patient. In the second row, a Gaussian perturbation in the parameters was applied to the parameters of the biokinetic model with σ = 10%. The dashed black line shows the best fits (with no perturbation). The activity curves are normalized to the injected activity (A0). [Color figure can be viewed at wileyonlinelibrary.com]

Due to the relative ranges of positrons and photons, the major contribution to the absorbed dose is due to positrons. Therefore, we could consider these results on CA uncertainties to be directly correlated with dose uncertainties, avoiding the need to run several simulations. Nonetheless, for one patient we have performed four GATE simulations, one with CAs calculated from optimal kinetic parameters, and three with different sets of parameters which lead to different CAs. In Fig. 3 we show absorbed doses in each organ for each simulation. Variations in the dose, ΔD, are linearly correlated with variations in CA, ΔCA, also shown in Fig. 3. The correlation follows ΔD = 0.7 ΔCA. The slope parameter is not one because part of the dose is due to photons that come from distant voxels/organs, and there is change in the proportions of those two contributions for different k‐values. This slope corresponds to the fit of the whole set of organs: slopes for individual organs could be higher/lower depending on whether the organ's uptake is high/low. From this correlation and the study of uncertainties in CAs due to the biokinetic model, uncertainties in computed doses due to biokinetics should be below 10% for most organs.

Fig. 3

Study of the effect of biokinetic uncertainties: (a) relative differences among MC‐calculated mean doses to different organs/tissues for one patient with four different combinations of biokinetic parameters. D0 refers to doses obtained with baseline parameters for this patient, while Di refers to doses obtained with three different perturbations of such baseline parameters. (b) Correlation between differences in cumulated activities and differences in mean absorbed doses. Different symbols correspond to different simulations. The dashed line shows the best‐fitting curve, ΔD = 0.7 × ΔCA, while the solid line shows the curve ΔD = ΔCA, where ΔCA = CAi/CA0 −1 and ΔD = Di/D0 −1.

We show and discuss the effect of bladder emptying and bladder segmentation in the Supplementary Materials (Figures SM11 and SM4). In short, bladder emptying time and segmentation mostly affect the doses to the bladder wall and nearby organs, and have a negligible effect on the rest of the organs/tissues.

Doses reported from now on correspond to bladder emptying 4h postinjection, and segmentation described as Case C (PET bladder contours are used to calculate cumulated activities, including bladder emptying, and to compute doses; see Supplementary Materials, sections SM2.2 and SM3.2).

Individualized MC dosimetry

In Fig. 4 we present boxplots of MC‐calculated mean doses per unit of activity in different organs/tissues. We also show mean doses to our cohort (notice that patients have been injected different activities, therefore the higher dispersion of these data). Most median dose/activity values are within 0.01–0.03 mGy/MBq, with higher doses delivered to the bladder wall, major veins, and brain (medians of 0.058, 0.060, 0.066 mGy/MBq, respectively, and outliers reaching 0.130, 0.077, 0.088 mGy/MBq). In Fig. 5 we show the dose distribution in representative slices of a patient.

Fig. 4

(a) Boxplots of MC‐calculated mean doses per unit of injected activity to different organs/tissues for the patient cohort. (b) Boxplots of MC‐calculated total mean doses for our patient cohort (notice that patients have been injected different activities, therefore the higher dispersion of these data).

Fig. 5

Dose distribution in representative slices of a patient. Notice that the scale has been saturated to 30 mGy to facilitate differentiation among relevant organs (doses in the bladder reach higher values, up to 100 mGy). [Color figure can be viewed at wileyonlinelibrary.com]

The results for calculated effective dose per administered activity are 0.019 ± 0.003 mSv/MBq, which matches the previously reported value of 0.019 mSv/MBq. We compared the residence time obtained with our biokinetic model to values taken from the biokinetic model published in the ICRP‐106, generally obtaining good agreement for most organs but for the heart (Table SM5).

Comparison with generic phantom‐based dosimetry

In Fig. 6 we present mean doses per unit of activity, calculated with individualized MC and phantom‐based methods for seventeen organs/tissues included in those formalisms. Boxplots of relative differences are presented in Fig. 7 which shows large differences in several cases, ranging from −82% to + 262%.

Fig. 6

Boxplots of mean doses calculated per unit of injected activity with (a) individualized MC, and phantom‐based methods: (b) Cristy–Eckerman, (c) Stabin, and (d) ICRP‐133 for 17 organs/tissues included in those formalisms. [Color figure can be viewed at wileyonlinelibrary.com]

Fig. 7

Boxplots of mean doses calculated per unit of injected activity with individualized MC and phantom‐based methods. For each organ/tissue we report differences between MC and Stabin's (left, light gray), Cristy–Eckerman (center, dark gray), and ICRP‐133 (right, black).

We have studied the significance of the differences between the organ doses evaluated with our Monte Carlo methodology and generic phantom‐based dosimetry. We have limited this comparison to the ICRP‐133 model, as it is the most recent, and, likely, most accurate one. In Fig. 8 we show p‐values for each organ. The differences are significant for several of them, namely, the remaining tissues, adrenals, bladder wall, bones, upper large intestine, heart, pancreas, skin, and stomach wall.

Fig. 8

Significance study of differences between organ doses evaluated with individualized Monte Carlo and generic phantom‐based dosimetry (ICRP‐133 model): p‐values obtained from a paired samples t test. [Color figure can be viewed at wileyonlinelibrary.com]

In Fig. 9 we show the effect of including mass‐scaling corrections in the results obtained with the generic phantoms. In general, this correction does not improve the agreement between Monte Carlo and generic phantom‐based dosimetry. In Figure SM12 we compare MC and generic phantom doses when cumulated activities are calculated by using the biokinetic model presented in this work for the MC (adjusted to each individual patient through an optimization), and the ICRP‐106 model (with and without mass scaling) for the generic phantoms. The use of nonindividualized cumulated activities magnifies the differences between MC and generic phantom calculations, as expected.

Fig. 9

Boxplots of differences in mean dose per unit of injected activity between MC and generic phantom calculations with/without mass‐scaled S‐factors: (a) ICRP‐133, (b) Stabin, (c) Cristy–Eckerman (CE).

Discussion

There is a growing interest in developing methods to compute precise, individualized internal dose distributions in patients undergoing medical procedures with radiopharmaceuticals. Several MC methods have been developed for this aim, but to date, they have been barely used in patients. In this work, we present a MC‐based methodology for individualized internal dose calculations. The calculation method relies on CT images (for organ segmentation and dose deposition), PET images (for organ segmentation and distributions of activities), and a biokinetic model (to obtain cumulated activities). The software vGATE 8.1. was employed to carry out the MC calculations. GATE was chosen due to its flexibility to handle CT images and DICOM source files containing patient‐specific cumulated activities. Dose calculations were performed remotely in a computer cluster. While in this preliminary study we have used modest computing resources, this implementation shall allow us to scale up them by using supercomputers with large numbers of nodes, and to speed up computation times to minutes.

The comparison between dose distributions obtained with GATE and EGSnrc showed similar results both in phantom and in CT patient geometries. Important differences were observed in air regions — the origin of these differences could not be identified and may be related to different implementations of the relevant physics in each code. However, the importance of these differences in the clinical application is negligible, as the dose in air‐filled cavities is not interesting, and it is usually set to zero. Regarding the calculation of SAFs in the ICRP‐133 male phantom, the GATE results agree within ± 5% with tabulated values for most organs. Higher differences arise in far apart organs for low energies electrons. This effect is most likely due to differences in the physics packages in each MC code. A comparison between MC calculated and experimentally measured doses in a Jaszczak phantom filled with an 18F solution also showed good agreement when considering the experimental and simulation uncertainties.

We have applied our methodology to a cohort of patients who underwent FDG‐PET studies and found significant differences between MC‐individualized and phantom‐based nonindividualized doses. The median MC‐calculated dose/activity values are within 0.015–0.03 mGy/MBq, with higher doses delivered especially to the bladder wall, major veins, and brain (medians of 0.059, 0.059, 0.065 mGy/MBq, respectively). Comparison to values obtained with nonindividualized phantom‐based methods has shown important differences in many cases (ranging from approximately −80% to + 250%).

We have performed a statistical test of such differences, finding significance (P < 0.05) for several organs/tissues, namely, the remaining tissues, adrenals, bladder wall, bones, upper large intestine, heart, pancreas, skin, and stomach wall.

The addition of patient‐specific information (i.e., organ mass scaling) did not reduce the differences with GATE for all organs or models. In general, we observed that their improvement was conditioned by several factors. For example, an excessive mass deviation of the segmented organ from that of the phantom tended to enlarge the discrepancies. This is the case for the thyroid, whose mass is greatly above the average for several of our patients. Another factor was the amount of cumulated activity. Organs whose dose depended mainly on the positronic contribution, for example, the liver, tended to give better approximations. Lastly, heterogeneous vs homogeneous activity distribution could explain why easy‐to‐segment high‐uptake organs do not fully match the GATE results, as is the case of the brain.

We have devoted a large part of this study to analyzing sources of uncertainty. It is important to remark that uncertainties in computed dose values are high. These uncertainties arise from two sources mainly: the biokinetic model to include activities during the uptake of the radiopharmaceutical and bladder emptying, and the segmentation of organs that may change volume in between the CT and PET studies, in particular, the bladder. Regarding the former, they could be reduced by performing dynamic PET studies. However, it does not seem clinically feasible to perform dynamic PETs for every single patient to improve the accuracy of internal dosimetry calculations, and users may be forced to rely on biokinetic models (at the cost of increasing the uncertainty of the results). In this work, we have combined a static PET study, with a biokinetic model and constraints on time‐activity curves derived from dynamic PET studies. Our study suggests that while this approach is subject to uncertainties, they may be acceptable when evaluating doses in patients undergoing PET studies. Similar uncertainties should not be acceptable when evaluating doses in molecular radiotherapy, yet it seems important to notice that the large half‐life of those radioisotopes may impact the uncertainties associated with the biokinetic modeling. Regarding uncertainties arising from segmentation, they can be limited if using dual PET/CT studies. In other situations, advanced image fusion tools may prove necessary.

Nonetheless, it is important to remark that these uncertainties are type B, and they apply to both MC and phantom‐based doses, as they affect the calculation of the cumulated activities that are used by both our MC and generic phantom‐based methods. Therefore, even if large, they do not compromise the differences found between doses obtained with our MC method and with phantom‐based methods.

Conclusions

The doses that were obtained with individualized MC calculations were in some cases significantly different from those obtained with generic phantom nonindividualized methods in a cohort of 14 patients undergoing FDG‐PET studies.

The methodology presented in this work is viable and useful to calculate individualized internal dose distributions in patients undergoing medical procedures involving radiopharmaceuticals, with higher accuracy than phantom‐based methods, fulfilling the guidelines provided by the European Council directive 2013/59/Euratom. Remote computation in supercomputers appears to be a promising option to obtain individualized MC dose distributions in short times and would allow for the implementation of these methods in many hospitals which do not host large computing facilities in‐house. The extensive pre‐ and postprocessing required for such calculations could be easily automatized and merged in a software package. The exception remains in the segmentation process which has to be manual and thus is time‐consuming. Nevertheless, promising advances could be made in this field in the coming years.

The methodology could be straightforwardly extended to PET studies involving different radiopharmaceuticals, and also to SPECT/CT studies, provided that there is a good quantification of the activity per voxel in the SPECT image and correction of partial volume effects. While the benefit–cost ratio of demanding computations to determine doses for diagnostic procedures can be debatable, such techniques may serve to evaluate and improve phantom‐based methods.

More importantly, this methodology can be extended to patients undergoing molecular radiotherapy, provided that there are images and models to know the particular biokinetics of the radiolabeled drug. In this regard, it would be important to develop reliable biokinetic models that can assist these calculations.

Ethical Approval and consent to participate

All procedures performed in studies involving human participants were approved by regional (CEIm‐G; study code: JPM‐TRNT‐2019‐01), and national (Agencia Española de Medicamentos y Productos Sanitarios, AEMPS, study code JPM‐IOD‐2019‐01) ethical committees in accordance with the ethical standards of the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Informed consent was obtained from all individual participants included in the study.

Conflict of interest

The authors declare that they have no competing interests.

Table SM1. Characteristics of patient cohort.

Table SM2. Segmented organs in the patient cohort and organ acronyms.

Table SM3. Optimal coefficients of the biokinetic model k [h‐1] for each patient. These are the set of coefficients that were used during the reported dose estimative.

Table SM4. Optimal blood percentage for each patient.

Table SM5. Residence times in the patient cohort and comparison with ICRP‐10643.

Figure SM1. Diagram of the biokinetic model used in our work. This model is a simplification of the more complex model presented in38.

Figure SM2. Static PET of representative slices of a patient. Notice that the color map is saturated, the real maximum at 6.4 × 10−4 MBq.

Figure SM3. CA source of representative slices of a patient.

Figure SM4. Illustration of the different approaches to bladder segmentation (Cases A, B, C, see section “BladderSegmentation” in this document for further information), and effect on the mean doses in each organ (doses relative to Case B). The inner and outer circles represent the bladder as observed in the CT and PET, while the solid and dashed lines represent the contours used to calculate cumulative activities and to compute doses respectively.

Figure SM5. Simulation of depth dose deposition of a 5 × 5 cm2 photon beam incident from z = −15 cm of different energy in a parallelepiped geometry consisting on water, and two layers of air (−10 cm) and lung (−5 cm). The plots represent depth dose profiles calculated with EGSnrc (dashed line) and GATE (solid line). Energies are read from left to right and top to bottom as: (a) 200 keV, (b) 511 keV, (c) 1000 keV, (d) 3000 keV. A marked difference between both calculations is easily observed in air for the lowest energies.

Figure SM6. Simulation of depth dose deposition in a simple setup in GATE against EGSnrc.

Figure SM7. Simulation of depth dose deposition of a 10× 10 cm2 511 keV photon beam in the geometry of a patient: 2D slice and depth dose profile calculated with EGSnrc (dashed line) and GATE (solid line). A marked difference between both calculations is easily observed in low‐density regions (exterior and stomach.

Figure SM8. Simulation of depth dose deposition in a patient's CT geometry in GATE against EGSnrc.

Figure SM9. Specific Absorbed Fractions (SAFs) calculation for the ICRP‐133 male phantom depending on the emission energy (circles) and comparison with the simulations results in GATE (dots). Each panel represents calculated SAFs for a source‐target pair and particle: (a) Liver–Liver, electrons; (b) Liver–Brain, electrons; (c) Liver–Liver, photons; (d) Liver–Brain, photons. Note that x‐ and y‐axis scales are not equal in every plot.

Figure SM10. Experimental validation of GATE MC simulations (red markers) by comparison with measurements in five dosimeters (black markers). The panels show dose readings (a) outside and (b) inside the phantom. (c) Shows a sketch of the location of each dosimeter in the experimental setup. Uncertainties are reported with a coverage factor of k=1.

Figure SM11. Relative effect of different bladder emptying times on organ doses: 2.5 h (leftmost point of each pair, red color) and 4 h (rightmost point of each pair, black color), relative to no emptying.

Figure SM12. Boxplots of differences in mean dose per unit of injected activity between MC (individualized biokinetic model) and generic phantom calculations with cumulated activities calculated by using the ICRP‐106 biokinetic model.

Click here for additional data file.