EANM practical guidance on uncertainty analysis for molecular radiotherapy absorbed dose calculations

A framework is proposed for modelling the uncertainty in the measurement processes constituting the dosimetry chain that are involved in internal absorbed dose calculations. The starting point is the basic model for absorbed dose in a site of interest as the product of the cumulated activity and a dose factor. In turn, the cumulated activity is given by the area under a time–activity curve derived from a time sequence of activity values. Each activity value is obtained in terms of a count rate, a calibration factor and a recovery coefficient (a correction for partial volume effects). The method to determine the recovery coefficient and the dose factor, both of which are dependent on the size of the volume of interest (VOI), are described. Consideration is given to propagating estimates of the quantities concerned and their associated uncertainties through the dosimetry chain to obtain an estimate of mean absorbed dose in the VOI and its associated uncertainty. This approach is demonstrated in a clinical example.


Introduction
Internal dosimetry following the administration of radio labelled pharmaceuticals for diagnostic and therapeutic purposes is a prerequisite for radiation protection, imaging optimization, patient-specific administrations and treatment planning. The medical internal radiation dose (MIRD) schema [1] has become the most widely accepted formalism for internal dose calculations. The general approach in medicine to determine the validity of a measurement is to compare the accuracy and precision against a "gold standard" measurement. To date, investigations of uncertainties for internal dosimetry have mainly used phantoms or simulated data [2][3][4] as the gold standard comparison. However, due to the diversity of dosimetry data, a subset of phantom experiments cannot necessarily validate the accuracy of measurements made for an in vitro population. It is therefore more appropriate to express the accuracy of a result by characterizing the uncertainty. This involves identification of the major processes and variables within the dose calculation and evaluation of their potential effect on the measurement. An uncertainty estimate should address all systematic and random sources of error and characterize the range of values within which the measured value can be said to lie with a specified level of confidence. The general relevance of Preamble The European Association of Nuclear Medicine (EANM) is a professional nonprofit medical association that facilitates communication worldwide among individuals pursuing clinical and research excellence in nuclear medicine. The EANM was founded in 1985. This guidance document is intended to assist practitioners in providing appropriate nuclear medicine care for patients. The rules provided in the document are not inflexible or requirements of practice and are not intended, nor should they be used, to establish a legal standard of care. The ultimate judgment regarding the propriety of any specific procedure or course of action must be made by medical professionals taking into account the unique circumstances of each case. Thus, there is no implication that an approach differing from the guidance, standing alone, is below the standard of care. To the contrary, a conscientious practitioner may responsibly adopt a course of action different from that set out in the guidance when, in the reasonable judgment of the practitioner, such course of action is indicated by the condition of the patient, limitations of available resources or advances in knowledge or technology subsequent to publication of this guidance document. The practice of medicine involves not only the science but also the art of dealing with the prevention, diagnosis, alleviation and treatment of disease. The variety and complexity of human conditions make it impossible to always reach the most appropriate diagnosis or to predict with certainty a particular response to treatment. Therefore, it should be recognized that adherence to this guidance document will not ensure an accurate diagnosis or a successful outcome. All that should be expected is that the practitioner will follow a reasonable course of action based on current knowledge, available resources and the needs of the patient to deliver effective and safe medical care. The sole purpose of this guidance document is to assist practitioners in achieving this objective.
performing and providing uncertainty information has been discussed in previous guidelines [5]. Flux et al. [6] provided a method to determine the uncertainty of wholebody absorbed doses calculated from external probe measurement data. Whilst whole-body dosimetry is used to predict toxicity in some procedures [7,8], organ and tumour dosimetry are required for treatment planning and in cases where haematotoxicity is not the limiting factor for treatment tolerance.
Specific aspects of uncertainty within the dosimetry chain have been addressed, including the selection of measured time points, [9,10], the chosen fit function [11,12] and uncertainty of model parameters. A comprehensive analysis of propagation of every aspect of the dosimetry calculation chain has yet to be obtained.
Gustafsson et al. [13] adopted a Monte Carlo (MC) approach to investigate the propagation of uncertainties to obtain an uncertainty in estimated kidney absorbed dose in 177 Lu-DOTATATE therapy, using simulated gammacamera images of anthropomorphic computer phantoms. In principle, this approach allows all aspects of the dosimetry process to be taken into account, but the need for multiple samplings from the assigned probability distributions for the quantities involved makes it computationally intensive, and its use for uncertainty assessment in an individual patient basis is challenging.
The Joint Committee for Guides in Metrology (JCGM) Guide to the Expression of Uncertainty in Measurement (GUM) [14] provides a generalized schema for propagating uncertainties. This EANM Dosimetry Committee guidance document provides recommendations on how to determine uncertainties for dosimetry calculations and apply the law of propagation of uncertainty (LPU) given in the GUM to the MIRD schema. This guidance document is presented in the form of an uncertainty propagation schema, and the recommendations are designed to be implemented with the resources available in all nuclear medicine departments offering radionuclide therapy, and are presented using terminology and nomenclature that adhere as far as possible to the GUM.
The uncertainty propagation schema examines each step of the absorbed dose calculation to estimate the standard uncertainty in the mean absorbed dose measured at the organ or tumour level using SPECT imaging. The examples given have been simplified and concentrate only on the mean absorbed dose to a target. However, the approach can be used in different scenarios and expanded to include more complex dose calculations, including cross-dose and multiexponential time-activity curves (TACs). Similarly, aspects of the methodologies described can be implemented in different applications of dosimetry such as those utilizing a hybrid imaging approach or those used to generate 3D dose maps.

Theory
The law of propagation of uncertainty A generic multivariate measurement model is: where is a vector of n generic input quantities X 1 , …, X n and is a vector measurand of m output quantities Y 1 , …, Y m . GUM Supplement 2 [15] gives a generalization of the LPU: where V y is the output covariance matrix associated with y (the estimate of Y). V x is the input covariance matrix associated with the estimate of X, and G x is the sensitivity matrix associated with x, defined as: where ∂f i /∂x j denotes ∂f i /∂X j evaluated at X = x. Element u(x i , x j ) of V x is the covarience associated with x i and x j , and u(x i , x i ) is equal to u 2 (x i ), the squared uncertainty associated with x i . x and V x are obtained from available knowledge, whether statistical (for example, repeated observations) or nonstatistical (for example, expert judgment), about the input quantities.
For a generic scalar measurement model, Eq. 1 becomes Y = f(X), where Y is a scalar quantity and f is a scalar function. Propagation of uncertainty for the estimate y of Y can be achieved using the matrix form of the LPU [15]: where u 2 (y) represents the variance (squared standard uncertainty) associated with the estimate y, and is the gradient matrix in which the ith element denotes the partial derivative of f with respect to the quantity X i evaluated at x. For a two-variable function, Y = f(X 1 , X 2 ), Eq. 8 can be expanded to give: For a two-variable multiplicative function, Y = X 1 X 2 , Eq. 10 can be written in the form: If the two variables X 1 and X 2 are mutually independent, the covariance term of Eq. 11 is zero, and therefore the standard fractional uncertainties u(x 1 )/x 1 and u(x 2 )/x 2 are simply added in quadrature.

Absorbed dose
For situations where the target volume is the source activity volume and the contribution of absorbed dose from other target organs can be considered negligible, a simplified form of the MIRD equation can be used in which mean absorbed dose D is expressed as the product of the cumulated activityÃ and the S-factor (sometimes called the dose factor) S: D ¼ÃS: ð12Þ Following the above notation, D is written as D ¼ fÃ; S À Á ¼ÃS, and the standard uncertainty u D À Á is evaluated at estimates ofÃ and S according to Eq. 11: It follows that the standard uncertainties uÃ À Á and u(S) and the covariance uÃ; S À Á are needed to obtain the standard uncertainty u D À Á associated with D. For the general form of the MIRD equation with meaningful contributions outside the target volume (cross-dose), uncertainties and covariances associated with additional quantities of the formÃ and S should also be considered.
The need for the covariance term of Eq. 13 may not be obvious on first inspection as calculations ofÃ and S are often considered separately, that is, one is derived from scintigraphy data and the other from simulations. However, as can be seen from Fig. 1 (a flow diagram of a typical dosimetry protocol), determination ofÃ and S both rely on a measurement of volume, and therefore a covariance exists between the two parameters. This EANM Dosimetry Committee guidance document describes how the uncertainty in the volume measurement and other confounding factors within the dosimetry protocol can be propagated to estimate an overall uncertainty in absorbed dose.

Volume uncertainty
The volume or mass of an organ or tumour is generally obtained from a volume of interest (VOI) outlined on anatomical or functional imaging data. It is therefore possible to estimate the outlining accuracy by considering factors that affect delineation. The method used will depend largely on the information and resources available at the time of outlining and the method employed by the operator or operators to define the VOI. The following concerns an outline drawn manually by a single operator across all images that comprise the dosimetry dataset.

Operator variability
For any given dosimetry dataset a number of independent VOIs could be drawn by different operators. Ideally, an average VOI boundary of volume v would then be generated and used in the calculation of absorbed dose. In practice, an average VOI boundary cannot be generated and an absorbed dose calculation is based upon the VOI drawn by a single operator. An alternative approach is therefore to estimate operator variability using historical datasets. In this case outlines previously generated by M different operators for N similar VOIs are used. The assumption made is that the VOIs are sufficiently similar with respect to the scanning modality and volume geometry that the fractional uncertainty is the same across all datasets. The VOI outlined by the current operator is then regarded as a random VOI drawn from the populations of outlines represented by the historical data. The standard uncertainty u(v) associated with the drawn volume is then expressed as: where v hist i is the average of M operator volumes for the historical dataset i, and s v hist i ð Þ is the standard deviation of the historical dataset I, The historical datasets should be carefully chosen to match the current study, as differences could lead to an inaccurate estimate of the final standard deviation.

Analytical approach
When historical outlines are not available, it is possible to use an analytical method to determine uncertainty. This approach provides an estimate of the most significant contributions to the uncertainty in the outlining process but is not necessarily exhaustive.
Given that any outlined VOI will be digitized into voxels, the extent of the VOI will be defined, approximately, by the subset of voxels through which the boundary of the VOI passes. The uncertainty of the outlined volume will hence depend on voxel size. Figure 2 shows an outline (Fig. 2a) and the effect of different voxel sizes (Fig. 2b, c).
The mass of a spherical volume may be obtained from an estimate of the diameter of the volume, where the diameter d is measured as the distance between two extreme points, P I and P j , the locations of which can be determined within one voxel dimension (Fig. 3).
Evaluation of the standard uncertainty associated with d can be considered as type B (nonstatistical), as in the GUM  [14]. Given that there is no specific knowledge about the location of point P i on the boundary, other than that it lies within the appropriate boundary voxel, there is a uniform distribution of possible values with variance associated with the mean value given by formula 7 in [14]: where a is one voxel width and u(P i ) is to be interpreted as the standard uncertainty associated with P i when measured on the diametric line between P i and the centre of the sphere. As diametric measurement is the distance between two extreme points, application of the variant of LPU [14] (formula 11a in the GUM) yields (assuming no correlation) the standard uncertainty u vox (d) associated with diameter d due to voxelization: With hybrid imaging it is often possible to use morphological information from CT imaging to aid functional VOI delineation. In this situation the VOI is drawn on the CT dataset and copied to the registered SPECT image. The coordinates of the original boundary will therefore be rounded to the nearest voxel coordinates of the SPECT image. Hence, the SPECT voxel size should be used in Eq. 17.
For many scintigraphic imaging processes the defined voxel size is less than the spatial resolution of the system and therefore the use of Eq. 17 would result in underestimation of the actual uncertainty when the VOI is drawn directly on the SPECT image. To provide a more reliable uncertainty the spatial resolution of the image system must also be considered. Consider a profile through an object approximated as a step function convolved with a Gaussian point spread function (PSF) with the full-width at half-maximum (FWHM) equal to the spatial resolution of the imaging system. The uncertainty in edge definition can be described by the gradient of the convolved step function, where the gradient profile is equal to the Gaussian PSF with standard deviation σ ¼ FWHM 2 ffiffiffiffiffiffi 2ln2 p . As the diametric measurement is the distance between the boundary locations on the profile, application of the variant of LPU [14] (GUM formula 11a) yields the standard uncertainty u res (d) associated with diameter d due to spatial resolution: For situations where both voxelization and resolution contribute significantly to diametric uncertainty (such as outlining directly on a SPECT image), Eqs. 17 and 18 are summed to give the combined uncertainty associated with d: In practice the diameter is not measured and only the volume is reported. However, conceptually the volume is determined through an infinite number of diametric measurements, and the mean diameter of the VOI can therefore be considered. The standard uncertainty u(d) translates into a standard uncertainty u(v) associated with the volume v delineated by the outline. Hence, for some positive constant k, The application of the variant of LPU [14] (GUM formula 12) yields a relationship between the relative volume and diametric uncertainties due to voxelization and resolution: Hence, a fractional standard uncertainty associated with a volume is three times the fractional standard uncertainty associated with the mean diameter of that volume.

Count rate
The total reconstructed count rate, C, within a VOI depends on the VOI delineation, and can be described as a function of volume. Propagation of volume uncertainty into the measurement of count rate is therefore required. As no prior knowledge of the count distribution is generally available, the variation in C within the VOI must be approximated.
A uniformly distributed spherical count rate density H(ρ) with volume v true of radius r, with a total emission count rate of C total , can be described in spherical coordinates as:  Fig. 3 Signal intensity profiles demonstrating that the gradient of a Gaussian blurred function can be described by the Gaussian function where ρ is the radial distance from the centre of the sphere. Due to the spatial limitations of the measuring system the apparent density is described as the spherical volume convolved with a 3D Gaussian function [16]: where σ is the measured standard deviation describing the width of the 3D Gaussian function. Therefore, an observed count rate density distribution can be described as: where * denotes convolution in three dimensions, which can be determined analytically [17] and re-expressed as: The function F(ρ) and that of perfect resolution H(ρ) are shown in Fig. 4a; images of these distributions are shown in Fig. 4b.
Using Eq. 25, the count rate C measured within a VOI of volume v and radius ρ can be expressed as: as shown in Fig. 5, where C is described as the area under the curve. A plot of this function with increasing ρ and v is given in Fig. 5a.
As there is an uncertainty associated with the drawn VOI boundary, the standard uncertainty u(C) associated with C is obtained using the gradient of C at v (Fig. 5d) and the volume standard uncertainty u(v): For a VOI of volume v, where the radius ρ = r, the expression for F(ρ) as given by Eq. 25 can be substituted into Eq. 27: The ratio of the total emission count rate from a source to the count rate measured within a VOI defining the physical boundary is referred to as the recovery coefficient [18] (see section Recovery coefficient): Therefore, the standard uncertainty u(C) associated with C can be rewritten as: Recovery coefficient There are a number of methods to correct for the observed "spill out" of counts from an object due to partial volume effects [16]. The simplest and most widely applied method is to divide the observed count rate by a recovery coefficient determined from phantom data. Dewaraja et al. [19] recommend imaging multiple phantoms of various sizes and geometries using the same acquisition and processing parameters as used for the patient data. An appropriate recovery coefficient is then selected based on the estimated object volume. The recovery coefficient determined from a phantom of volume v true is defined as: where C is the observed count rate measured within a VOI matching the true volume of the phantom and C total is the count rate of all counts originating from the phantom [18].
(a) (b) Fig. 4 a Count density as function of radius for a spherical object with true radius r for a system with ideal resolution (red step function) and realistic system (green curve). b Two-dimensional image planes through the three-dimensional functions H(ρ) and F(ρ) (see text) A dataset comprising a series of volumes and the corresponding factors on the right side of Eq. 32 is fitted by an empirical function of appropriate form. Such a function will have adjustable parameters b = [b 1 ,…, b q ] T and will provide a means for estimating a recovery factor specific to the volume under investigation. The standard uncertainty u(R) associated with the recovery factor can be derived from the fitted parameters b. A covariance matrix V b of dimensions q × q corresponds to the estimate of b determined by ordinary least squares fitting under the assumption of equal uncertainty in all data points that make up the dataset. For a perfectly specified volume v, the squared standard uncertainty associated with R is, using Eq. 8: where g b is the matrix of dimension q × 1 containing the partial derivatives of first order of R with respect to b. The covariance matrix V b can be determined as a by-product of the least squares fitting process. In reality the standard uncertainty u(R) obtained in this manner for a given volume will underestimate the actual uncertainty. To provide a more realistic value for u(R), the standard uncertainty u(v) associated with the clinical outlined volume v has to be taken into consideration. Since this volume is independent of the recovery curve parameters, there will be no covariance associated with v and any of the b j . Accordingly, V b in Eq. 33 is replaced by: where 0 is a matrix of zeros of dimension q × 1 and g b is extended by one element, namely the partial derivative of first order of R(v) with respect to v.

Calibration factor
The final conversion of a partial volume-corrected count rate to activity is achieved by the use of a quantification or calibration factor. The sensitivity of the system is determined by measuring the total count rate C ref of a source of known activity A cal , commonly referred to as a "standard", under the same acquisition and reconstruction conditions as the study data: The quantification factor will therefore depend on the standard activity measurement within the dose calibrator and the reconstructed count rate, and its associated uncertainty accordingly calculated. Methods to determine dose calibrator uncertainty are described by Gadd et al. [20]. Uncertainty in the measurement of reconstructed counts within the standard should be determined statistically from a series of nominally identical observations. The uncertainty of a single measurement can be obtained by calculating the mean and standard deviation of that series. The standard uncertainty u(Q) associated with Q can be determined by Eq. 11, a variant of LPU [14], that combines the fractional uncertainties of the dose . Fig. 5 Activity The expression relevant to the assessment of the uncertainty associated with the activity A i determined from the measured count rate C i in a target VOI at time t i is: Equation 37 corresponds to all measurement times t i , for i = 1, …, n. Using matrix notation: Equation 38 is a multivariate measurement model (section The law of propagation of uncertainty) with n + 2 input quantities Q, R, C 1 , …, C n and n output quantities With no loss of generality, in order to keep the mathematical expressions simpler than they otherwise would be, only two of these activity values, namely, A i and A j , are considered, measured at times t i and t j , respectively: Equation 39 is a bivariate measurement model with four input quantities Q, R, C i and C j , and two output quantities A i and A j . Using Eq. 7, and the input covariance matrix is: Since Q is independent of volume and hence independent of R, C i and C j , the covariance terms u(Q, R), u(Q, C i ) and u(Q, C j ) are zero. Application of Eq. 7 then gives: as the covariance matrix for A i and A j , the elements of which may be expressed as: It follows that the covariances u(R, C i ) and u(C i , C j ) associated with Q, R and the C i must be derived.
Regarding u(C i , C j ) for this situation, Eq. 4 or, more directly, by the use of formula F.2 in the GUM [14], yields: Using Eqs. 28 and 30 gives: Hence: As both the recovery coefficient R and the measured count rate C i depend on the VOI outline they can be expressed as functions of volume v. Again applying the GUM [14] (formula F.2) and using Eqs. 28 and 30 gives: which, after rearrangement, can be expressed as: After substituting the covariance expressions of Eqs. 47 and 49 into Eqs. 43 and 44 it can be seen that: Given the equal fractional uncertainties for all the A i and with perfect covariance between the A i and A j , it is appropriate to treat these uncertainties in a manner similar to a systematic error. Hence the fractional uncertainties in activity can be propagated into a systematic component of uncertainty for cumulated activity u sÃ À Á , where Time-activity curve parameters In addition to the systematic uncertainties associated with quantification and volume determination, uncertainties in the TAC data can arise from other sources such as image noise, patient motion, registration and other imperfect post-acquisition operations such as image reconstruction, including scatter and attenuation corrections. Due to the complexity of these operations, it is assumed that the uncertainties associated with the compensation for effects such as attenuation and scatter are negligible in comparison with the uncertainty associated with the compensation for partial volume effects [13]. It is therefore more appropriate to measure the causality of imperfects in these corrections, and to derive uncertainties in the fit parameters of the TAC. Estimates of the TAC parameters p = [p 1 ,…, p q ] ⊤ can be determined by fitting data points (t i , A i ), i = 1, …, n, where the t i denote the image acquisition times and A i the corresponding measured activities. A least squares approach is recommended, using nonlinear regression techniques to minimize the objective function with respect to p. Note that a weighting term to account for the activity uncertainty is not included due to the covariant nature of the uncertainty. Uncertainties of the fit parameters are then estimated using: where J p is the matrix of first-order partial derivatives of the TAC model with respect to p, evaluated at A. The TAC is generally represented as a sum of exponential functions. For the purpose of presentation, only the case of a single exponential function is described: where A 0 is the activity at time zero and λ is the effective decay constant, for which and

Cumulated activity
The cumulated activity is defined as the integral of the TAC from time t = 0 to ∞, which for a single exponential function is described simply by the ratio of the TAC parameters, that is: Application of Eqs. 8 and 9 to Eq. 57 is used to derive the component of uncertainty associated with random effects: and V p is the covariance matrix for the estimates of the TAC parameters p = [A 0 , λ] ⊤ given in Eq. 53. After expansion of these matrices the component of uncertainty associated with random effects inÃ can be expressed as: Random and systematic components can be combined by considering the general model: where x nom is the nominal value of some parameter x, and r and s are random and systematic effects, respectively. Then, applying LPU: For cumulated activity hence, using Eqs. 50, 51 and 60:

S-factor
Uncertainties associated with S-factors are somewhat less complicated than in the case of cumulated activity, and are predominantly influenced by the uncertainty associated with the volume. The general approach to determining S-factors is to choose a value calculated for a model that closely approximates the organ or region of interest. If a model of the corresponding size does not exist, a scaling can be applied to adjust the S-factor accordingly. Alternatively, an empirical S-factor versus mass representation can be obtained by fitting suitable S-factor data against mass [21,22]. The implicit assumption is that appropriate models exist for the clinical situation. There are uncertainties associated with deviations between the model and reality (for example, a tumour that is not spherical) but these are outside the scope of this framework. Figure 6 shows, on a log-log scale, an example of S-factor data for unit density spheres of different masses [23], empirically fitted by the function: It is possible to apply the same principles as employed in the previous section to determine a covariance matrix for the estimated parameters in the fitting function. However, in this instance the standard uncertainties associated with these fit parameters tend to be very small (<1%) and the mass uncertainty dominates. Therefore, these estimated parameter uncertainties can be ignored and, using Eq. 31, the standard uncertainty in S is: Given that mass is proportional to volume, and assuming a known tissue density with negligible uncertainty, the fractional standard uncertainty associated with S is: The fractional standard uncertainty associated with S is thus proportional to the fractional standard uncertainty associated with volume v, the proportionality constant being the magnitude |c 2 | of the slope of the fitting function on a log-log scale.

Absorbed dose
Having established standard uncertainty expressions forÃ and S, it is evident that both parameters have a dependence on volume. To determine a final uncertainty in the absorbed dose, the covariance between these parameters should therefore be determined. Applying the GUM, [15] (formula F.2) the covariance uÃ; S À Á is evaluated using: Lu-177 Fig. 6 Example plot of S-factor versus mass for the radionuclides indicated for unit density spheres An expression forÃ with respect to volume is difficult to derive. However, using LPU [15] (formula 12) and with only the systematic uncertainty component inÃ having a volume dependence, Using Eq. 51, the fractional standard uncertainty inÃ can be replaced with that of activity to give: Hence, ∂S ∂v is obtained from Eq. 65: with v in place of m, To provide ∂A i ∂v requires a re-expression of activity as a function of volume. Use of Eq. 37 and differentiating with respect to v yields: Using Eq. 26 gives: such that covariance betweenÃ andS can be expressed as: Having established expressions for covariance uncertainty inÃ and S, Eqs. 75, 67 and 64 can be used in Eq. 8 to give a final uncertainty in absorbed dose, given by: where gÃ ;S ½ and VÃ ;S ½ are the respective gradient and covariance matrices associated withÃ and S which, using Eq. 13 for the case of a single exponential TAC, can be written:

Patient example
An example to demonstrate the implementation of the approach described in this paper is given in the following sections, with details of the methodology used to obtain the absorbed dose data and the associated uncertainty analysis. The example given is that of a 47-year-old patient who presented with weight loss, lethargy and upper abdominal cramps. Upper gastrointestinal endoscopy showed a mass in the third part of the duodenum. A subsequent contrast-enhanced CT scan and 68 Ga-DOTATATE PET/CT investigation showed a 6.5-cm mass arising from the pancreatic head and a 3-cm mass within segment 4 of the liver, in keeping with a neuroendocrine tumour arising from the pancreas. The patient underwent 90 Y-DOTATATE radiopeptide therapy in combination with 111 In-DOTATATE for imaging. The administered activity was 4,318 MBq of 90 Y with 111 In given at a ratio of 1:25.

Image acquisition
Absorbed doses for the lesions were calculated using sequential 111 In SPECT acquisitions, performed at 19.7 h, 45.1 h and 66.5 h after administration, acquiring 64 projections in a 128 matrix for 60 s per view. Triple-energy window scatter corrections were applied to the projection data with 20% energy windows centred on the 171 keV and 245 keV photopeaks. The scatter-corrected data for each energy window were then added and reconstructed iteratively with a weighted attenuation coefficient based on the photopeak abundance as described by Seo et al. [24]. The reconstructed SPECT voxel size was 4.67 mm. 111 In-DOTATATE SPECT images are shown in Fig. 7 alongside the 68 Ga-DOTATATE PET/CT images. The primary tumour and the hepatic lesion are indicated on the images from the two modalities.

Volume
VOIs of the two lesions were determined using an adaptive thresholding technique, whereby a threshold for outlining was chosen based on the known threshold required to outline similar sized volumes on phantom data. As the VOI was drawn directly on the SPECT data, the voxelization uncertainty was combined with the spatial resolution element, given in Eq. 19.
The reconstructed system spatial resolution was determined directly from physical measurements of a point source in air. The measured FWHM was 0.9 cm. Volumes and uncertainty components are shown in Table 1.

Recovery coefficient
Recovery data were generated by imaging multiple phantoms with the same acquisition and processing parameters as described for the patient data. The phantoms consisted of a 20 cm × 12 cm cylindrical phantom within which smaller inserts could be placed. Insert volumes ranged from 0.1 ml to 200 ml and were filled with a known concentration of 111 In. For each insert two VOIs were drawn. The first set of VOIs were generated by selecting the appropriate percentage threshold to match the known physical volume of the insert. The second set of VOIs were drawn to encompass all counts (including spill out) that originated from the insert volume. A recovery coefficient for each insert was then determined using Eq. 32. The generated recovery curve is given in Fig. 8. The empirical function fitted to the example data takes the form of a twoparameter logistic function, with respect to volume v [25], namely: Fit parameters of the curve with associated uncertainties were determined using GraphPad Prism fitting software (La Jolla, CA, USA) and are detailed in Table 2. The covariance between the parameters was calculated as 0.0155, which can be expressed as a correlation coefficient, r, defined as: Equations 33 and 34 were used to combine the volume uncertainties with the recovery coefficient uncertainty for both lesions as shown in Table 3. Uncertainty estimates are given with and without the volume component. The importance of propagating the volume uncertainty into the calculation is clearly apparent for the smaller of the two lesions.

Count rate
Count rates for each lesion with each scan time are shown in Table 4 with the associated fractional standard uncertainties. The uncertainty in the VOI count rate is described in Eq. 31.    The covariance between the recovery coefficient and count rates, u(R,C), is defined in Eq. 48, which when combined with the empirical function given in Eq. 78 can be reexpressed as: Substitution of the fit parameters, volumes and count rates into this expression is used to generate the values for covariance which are shown in Table 5. Correlation coefficients relating to these covariance values are 0.99 and 0.92 for the liver and pancreatic lesion, respectively.

Calibration
The system was calibrated by imaging point sources of 111 In at various activities (8 MBq to 30 MBq) in air using the same acquisition parameters as for the patient scans. Images were reconstructed according to the clinical protocol and a spherical VOI was placed over the reconstructed point, ensuring that all counts from the source were contained. A plot of VOI count rate versus activity is given in Fig. 9.
The fractional standard uncertainty in activity was taken as 1.5%, the typical uncertainty for secondary standard calibrators for 111 In as given by Gadd et al. [20]. The statistical uncertainty from repeating the calibration measurement was taken from the standard deviation of the mean. Combining these uncertainties, as shown in Eq. 36, yields: Q ¼ 275 cps=MBq; with a standard error of 8 cps=MBq:

Activity
Calculation of 111 In activity within each lesion is obtained from Eq. 37 using the measured count rate C i , volumespecific recovery coefficient R and calibration factor Q. 111 In activity was converted to 90 Y activity by scaling the ratio of the administered activities and correcting for decay according to the different half-lives of the two isotopes, such that the 90 Y activity is expressed as: The variance associated with the measured activity is given in Eq. 50 and assumes negligible uncertainty in the administered isotope activities. It is therefore a simple case of substituting the relevant variance and covariance values for Q, R and C to form the required uncertainty in A i . These activities and associated uncertainties are shown in Table 6.

TAC fitting
The Gauss-Newton algorithm was used to minimize the objective function described by Eq. 52. A single exponential function was fitted to the data and the uncertainties in the fit parameters A 0 and λ were determined using Eq. 53. TACs for the two lesions are given in Fig. 10 with the fitted exponential   Fig. 9 Activity versus count rate of reconstructed point sources in air functions. Error bars on the data points represent the estimated standard uncertainty in activity. Solution parameters with associated random and systematic components of uncertainty for each TAC are shown in Table 7.

Cumulated activity
The covariance matrix for the solution parameters V p is given in Table 8. The product of the covariance matrix and the gradient matrix g p was used to determine the random component of the variance ofÃ described in Eq. 58.
Random and systematic components of uncertainty in the cumulated activity were determined according to Eqs. 50 and 60. These results and that of combined uncertainty (Eq. 63) are shown in Table 9.

S-factors
The S-factors for the lesions were determined by fitting 90 Y Sfactor data for unit density spheres against mass [21], empirically fitted by the function: The fit parameters of the curve with associated uncertainty were determined using GraphPad fitting software (La Jolla, CA, USA) and are shown in Table 10. As the standard uncertainties associated with these fit parameters are much less than the mass uncertainty of the two lesions the estimated parameter uncertainties can be ignored and Eq. 66 holds. Table 11 shows the determined S-factors for the lesions with associated uncertainties.

Absorbed dose
The uncertainty in the absorbed dose is determined from Eq. 13, for which the covariance uÃ; S À Á is required. Use of Eq. 75 to determine this covariance requires solving the partial derivative ∂R ∂v and substituting the determined parameters S;Ã; R; v; c 1 ; c 2 and the standard uncertainty u(v). For the recovery function defined in Eq. 78, the partial derivative is expressed as: Solutions for uÃ; S À Á are shown in Table 12. It can be seen from the correlation coefficients that the covariance betweenÃ and S is highly significant. In addition, the negative nature of    Table 13. Propagation of uncertainty can be visualized by examination of the fractional uncertainty of each parameter calculated along the dosimetry chain. Figure 11 gives uncertainties for the absorbed doses delivered to lesions and to normal organs. It can be seen that the small volume of the liver lesion has a significant impact on the larger fractional uncertainty compared to the larger lesion and organ volumes.
Using the methodology described, absorbed doses to lesions and normal organs were calculated. In addition, the treatment was repeated for four cycles and an equivalent methodology was employed. Dosimetry results for lesions and normal organs delivered at each fraction are presented in Table 14 and shown graphically in Fig. 12. A significant decrease in absorbed dose to the lesions was observed after the first cycle and an increase in absorbed dose to the kidneys after the fourth cycle.

Propagation of uncertainties from a tracer to a therapy study
A further potential source of uncertainty is the use of a pretherapy or concomitant diagnostic tracer study to predict the absorbed dose that would be delivered from a different therapeutic agent. It has previously been shown that uncertainty in the estimation of the biological half-life of the tracer will have an impact on the uncertainty of the absorbed dose calculated for the therapy procedure as a function of the relative values of the biological and physical half-lives [6]. The uncertainty in the therapeutic effective half-life can then be expressed as; where T eff (th) and T eff (tr) are the effective half-lives of the therapeutic and the diagnostic radionuclide, respectively. In the limiting case of infinite biological retention, the ratio of the uncertainties for the tracer and therapeutic agent will be the ratio of the physical half-lives. An uncertainty in the absorbed dose calculated for 111 In (physical half-life T phys = 67.3 h) will therefore produce a similar uncertainty in the absorbed dose for a 90 Y calculation (T phys = 64.1 h). However, a small uncertainty in, for example, an absorbed dose calculation for 68 Ga (T phys = 1.13 h) would propagate by a factor of~60 to give potentially significant uncertainty in a 90 Y calculation.  It is important to note that nonconformance, for example different administered amounts or affinities between diagnostic and therapeutic radiopharmaceuticals, will also introduce additional uncertainties into the prediction [26,27]. For example, it is assumed that the biokinetics of 111 In-and 90 Y-DOTATATE are equal, whereas the renal uptake of the indium-labelled compound might be higher [28].

Discussion
The methodology presented allows uncertainty analysis to be incorporated into absorbed dose calculations using the MIRD schema [1], the most widely adopted approach for molecular radiotherapy (MRT) dosimetry. The methodology is based on the recommendations described within the GUM [14] and necessarily involves the formation of covariance matrices for several steps of the dosimetry process.
The main objective of this uncertainty propagation schema is to evaluate the standard uncertainty in absorbed dose to a target. The tasks that directly support that objective are the determination of cumulated activity and the S-factor. The cumulated activity, given by the area under a TAC, is obtained from a sequence of quantitative images. Each activity value is expressed in terms of an observed count rate, a calibration factor and a recovery coefficient. The recovery coefficient is based on a recovery curve derived from multiple phantom scans. The presence of a common calibration and recovery factor in all activity values, and the covariance between volume, recovery and measured count rate, can be considered as a systematic uncertainty applied across all TAC data points, and therefore may be applied directly to cumulated activity. For the effects of uncertainties associated with random components of TAC data, a statistical approach using a "goodness of fit" measure is used.
Within the described schema particular functions are used to fit the acquired data, for example for the TAC, recovery and S-factor models. The choice of these functions is not discussed, and an obvious fit function from theory may not always be known. In this case an optimal function can be found, and uncertainties reduced by using model selection criteria and model averaging [10,29,30].
It is suggested that the major factors affecting uncertainty in the absorbed dose originate from the uncertainty in the delineation of the VOI. Two approaches to determine this uncertainty using statistical and analytical methods are presented. In this example an assumption is made that only a single VOI is applied to all datasets. An alternative approach involves the individual delineation of VOIs for each time point, for which the described methods may need to be varied, taking care to account for any commonalities applied across time points.
Propagation of these uncertainties to derive those further along the dosimetry chain requires the covariance between parameters to be evaluated. An understanding of the variation in VOI counts with VOI uncertainty is challenging as there is no prior knowledge of the count distribution. A method for estimating the count distribution is therefore proposed. However, this approach does not model noise or background counts spilling into the VOI. A more rigorous approach would be to determine a function for change in counts versus volume for the dataset being analysed. However, it is considered that   the approach suggested is sufficient since it does not overly complicate the methodology or require additional image processing or analysis, which is not available to the wider nuclear medicine community. An important feature of the schema is that it can be easily implemented using standard nuclear medicine image processing techniques. This feature is demonstrated in the clinical example in which absorbed dose calculations were performed using a standard image processing workstation and a commercial spreadsheet with curve-fitting software. Clinical implementation of this approach clearly demonstrates how different aspects of the dosimetry calculation can influence uncertainty. Uncertainty pertaining to a smaller lesion is clearly affected by the ability to define precisely the lesion volume and can be significant. For larger organs (such as the liver) volume delineation is less significant and the fit to the TAC begins to dominate. The ability to determine the source of larger uncertainties facilitates optimization of dosimetry protocols.
The clinical example given in the appendix demonstrates the importance of uncertainty in reviewing the significance of results. Figure 12 shows the variation in absorbed doses measured in different treatment cycles. With the presence of uncertainties indicated by error bars, it is possible to determine where a significant difference in delivered absorbed dose occurs. If absorbed dose measurements are to be used to aid future treatment (the goal of MRT dosimetry) it is possible that different treatment strategies could be adopted if the absorbed doses delivered are seen to be constant or decrease with sequential cycles. The uncertainty given in the example demonstrates the utility of the guidance to help identify aspects of the calculations that can be addressed to improve accuracy. It is important to note that the scale of uncertainties should be considered in relation to the range of absorbed doses that are delivered from standard administrations.
Whilst the clinical example demonstrates the use of the schema for SPECT-based dosimetry, the methodology can easily be adapted to suit alternative dosimetry protocols (that is, for multiexponential TAC models, external probe counting or 3D dosimetry). However, variations to the proposed schema should always follow the uncertainty guidelines set out by the GUM.
Uncertainty analysis is important for any measured or calculated parameter, whether physical or biological. Such calculations for MRT are rare [5] and leave room for systematic improvement. With the rapid expansion of MRTand an increase in the number of centres performing dosimetry, it is important for adequate interpretation of the data in clinical practice that measurement uncertainties are quoted alongside absorbed dose values. The application of uncertainty analysis may increase the validity of dosimetry results and may become the basis for quality assurance and quality control. Uncertainty analysis may help identify and reduce errors, aiming at an increased likelihood of observing actual dose-response relationships, which in turn would lead to improved treatment regimens.
Acknowledgments This guidance document summarizes the views of the Dosimetry Committee of the EANM and reflects recommendations for which the EANM cannot be held responsible. The recommendations should be taken into context of good practice of nuclear medicine and do not substitute for national and international legal or regulatory provisions.  The guidance document was brought to the attention of all other EANM Committees and of national societies of nuclear medicine. The comments and suggestions from the Radiopharmacy Committee and the Dutch and the Italian national societies are highly appreciated and were considered in the development of this guidance document.

Compliance with ethical standards
Conflicts of interest Jonathan Gear, Katarina Sjögreen Gleisner, and Mark Konijnenberg are members of the EANM Dosimetry Committee, Gerhard Glatting and Glenn Flux are members of the EANM Radiation Protection Committee. All authors declare that they have no conflicts of interest.
Ethical approval All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the principles of the 1964 Declaration of Helsinki and its later amendments or comparable ethical standards. This article does not describe any studies with animals performed by any of the authors.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.