Mapping 18F-FDG Kinetics Together with Patient-Specific Bootstrap Assessment of Uncertainties: An Illustration with Data from a PET/CT Scanner with a Long Axial Field of View

Visual Abstract

Hi gh-resolution dynamic whole-body PET scanning enhances the ability to map metabolic characteristics of tissue, particularly in the context of cancer.The current focus has been on dynamic PET studies with 18 F-FDG using the well-established Huang-Sokoloff 2-compartment (2C) modeling framework (1)(2)(3).Although 2C modeling has had widespread application in PET imaging, far beyond the brain setting in which it was developed, the biochemical understanding of the transporters involved in the metabolism of 18 F-FDG and their distribution across normal and cancerous tissues has evolved in the years since the Huang-Sokoloff construct was proposed (4)(5)(6)(7).The temporal and spatial resolutions of emerging scanners have transformed the ability to objectively assess the accuracy of the 2C framework to represent 18 F-FDG time-course data across the diverse tissues encountered in the human body.In this context, the assessment of 18 F-FDG kinetics based on more flexible nonparametric analysis approaches (8,9) may be necessary.The most recent implementation of the nonparametric voxel-level analysis scheme (9) is particularly efficient, largely because of an extensive reliance on quadratic programming techniques, and its nonparametric aspect provides an ability to apply an image-domain bootstrapping process for evaluation of statistical uncertainties in derived kinetic maps and associated biomarkers (10,11).Uncertainties in diagnostic information recovered from PET scans could augment decisionmaking for individual patients that is based on complex nonlinear radiomic metrics derived from a kinetic map.
The volume of data produced by a dynamic 18 F-FDG PET study on a state-of-the-art scanner with a long axial field of view (FOV) is a practical computational challenge for voxel-level analysis of kinetics.The bootstrap uncertainty assessment requires that comprehensive voxel-level analyses be applied to multiple simulated datasets, each created to match the full character and extent of the original data.This significantly adds to the computational challenge involved.
The work here uses a series of dynamic 18 F-FDG data acquired on a long-axial-FOV scanner (2) to investigate the approach.Apart from the demonstration of the practical feasibility of kinetic mapping with uncertainty evaluation, the analysis allows regional comparisons between nonparametric and 2C modeling results in terms of both derived kinetics and accuracy of data representation.

MATERIALS AND METHODS
An extended materials and methods description is provided in the supplemental materials (supplemental materials are available at http:// jnm.snmjournals.org)(12).

Patient Scans and Volumes of Interest (VOIs)
The data considered arise from a set of 24 patients with different types of cancer who participated in an institutionally approved 18 F-FDG PET/CT study at Bern University Hospital (KEK 2019-02,193).Details of the study were reported previously (2).In summary, PET scanning was conducted on a Biograph Vision Quadra scanner (Siemens) with a 106-cm axial FOV and a nominal in-plane resolution of 3.3 mm in full width at half maximum (13).Data were acquired in list mode starting 15 s before an intravenous bolus injection of 18  Automated segmentation algorithms based on CT and PET were used to define VOIs corresponding to several tissue structures, including gray and white matter in the brain, liver, lungs, kidneys, spleen, and bones (2).A further set of 49 VOIs corresponding to tumor tissue was identified by an experienced nuclear medicine physician.Finally, a VOI placed in the descending aorta was used to define the whole-blood arterial input function (AIF) used for kinetic analyses (2).Further scanning and study protocol details are available in the supplemental materials.

Parametric Imaging Techniques
Tissue Residue and Kinetic Parameters.When the Meier-Zierler (14) formalism is followed, the analysis assumes the PET-measured time course for a tissue region is represented as a convolution between the local AIF, C p , and the regional tissue residue function.Kinetic parameters are defined in terms of this residue (Fig. 1).Large-vessel vascular blood and distribution volumes (V b and V d , respectively) are evaluated as areas under the tissue residue.The apparent rate of retention or flux (K i ) of the tracer, measurable by PET over the scan duration, is the height of the residue at the end of the acquisition period.Also, the mean transit time of the tracer in the tissue and extraction fraction are defined as ratios of amplitude and integral measurements.A variety of approaches might be used to approximate the residue: a nonparametric method is used here.Patlak analysis uses a constant residue (15).Compartmental model forms, for example, the 1-compartment Kety-Schmidt ( 16) model for water and the 2C Huang-Sokoloff (17) model for 18 F-FDG in the brain, represent residues by positive linear combination of exponentials.In the 6-parameter 2C model, there is additive adjustment for an arterial signal.By adding a sharp residue element to the 2-exponential form, a Meier-Zierler residue is also available for this model.This allows residue-defined metabolic parameters for the extended compartmental model to be evaluated via the decomposition shown in Figure 1 (18).Supplemental materials provide a review of how Meier-Zierler residue parameters link with rate constants in the 2C model.

Nonparametric Residue Mapping (NPRM) of Kinetics. NPRM
approximates the voxel-level residue by the positive linear sum-ofbasis elements that have been selected by a cross-validation-guided analysis of a comprehensive collection of time courses produced by segmentation of all the available data in the study (10,18).Individual basis elements are of the form m k ðtÞ5 Here, R k is the basis element residue and D k is its associated delay factor.Note that cross validation is used to select the number of basis elements (K).Given the basis set, PET-measured voxel-level time-course data over the available set of J time frames, fzðt j Þ, j51, 2, . . ., Jg, is expressed as Eq. 1 Here, d and (a 1 , a 2 , … , a K ) are the unknown voxel-level delay and basisamplitude parameters, respectively, and є (t) represents (random) model error.A weighted least-squares criterion, with weights proportional to the product of the frame duration and the decay correction factor used to convert raw counts to decay-corrected tracer activity, is used for optimization of the unknown parameters.For any delay, the optimal set of a coefficients is found by quadratic programming.A crude grid search is used to optimize delay (10).
Bootstrap Assessment of Uncertainty.Model residuals across N voxels and J time frames, fz i ðt j Þ2ẑ i ðt j Þ, i51, . . ., N; j51, 2, . . ., Jg, are used to construct an image-domain data generation process (DGP) for bootstrapping.The DGP generates data according to where ẑðt j Þ5â 1 m 1 ðt j 2 dÞ1 . . .1â K m K ðt j 2 dÞ and the simulated error process, e*, mimic the stochastic character of analysis residuals.Analysis of bootstrapped datasets arising from the DGP leads to a set of bootstrapped kinetic parameter values at each voxel.The SD of these values estimates the voxel-level SE of the kinetic parameter.Similarly, the SEs for more complex quantities, such as the maximum-intensity projection (MIP) for a kinetic map, are created as the SD of the bootstrapped MIPs of the kinetic parameter (Fig. 2).Numeric studies (10,11) have shown that image-domain DGP bootstrapping matches the accuracy of the much more computationally intensive list-mode bootstrapping approach of Haynor and Woods (19).
The number of bootstrap simulations impacts the accuracy of the SEs it produces (20); this is discussed in the supplemental materials.

Statistical Analysis
NPRM kinetic analysis with 25 bootstrapped simulations is evaluated for each of the studies in the series.Results are examined in 4 separate ways.Technical details with formulas are in the supplemental materials.

Representation of VOI
where the random errors, є Ã i ðtÞ, are in units of SD and ŝe is an overall scale of the model error.In Equation 3, the factors ĉi and ŵj are scale-free quantities representing the relative uncertainty across voxels (i) and time frames (j).As the PET-measured activity scales with dose, the DGP error scale (ŝ e ) should also scale with dose; this is examined graphically.The overall axial pattern variation is described by the scale factor ĉi .In a uniform cylindric phantom, this has a familiar U-shaped pattern related to scanner sensitivity (10).With a patient in the scanner, the distribution of activity and attenuation is far from uniform.Physiologic patient motions, such as breathing, may also impact axial variation.Skewness is a key feature of iteratively reconstructed PET data.A histogram of scaled residuals shows how the DGP captures this aspect.After adjustment for spatial scale factors, the 3-dimensional power spectrum of the normalized residual process provides insight into the effective resolution of the scanning.Coordinatewise autocorrelation functions associated with the spectrum give insight into the actual resolution of the scanner.Again, physiologic movements may well lead to the actual resolution's deviating from what might be predicted on the basis of static phantom measurements.
SEs of VOI Kinetics.In theory, uncertainty in parameters recovered by kinetic model fitting should be proportional to the scale of the residual model error, but it may also be a function of the relevant sensitivity matrix for the model.We examine the relation between the bootstrap assessment of mean VOI kinetic SEs and suitable explanatory factors including the weighted-residual-sums-ofsquares fit of the VOI and the mean VOI kinetic values.For each kinetic parameter, linear regression analysis on a logarithmic SE scale is applied.Adjustment of this regression analysis based on the VOI type and the kinetics are explored.Regression predictions of SEs are graphically compared with the true.Correlation values are also summarized.

Illustration
Sample kinetic MIP maps with associated SEs obtained using the NPRM technique and bootstrapping are shown in Figure 2. A video of all coronal MIP maps is provided as Supplemental Video 1.Note that the dataset is the same as that used in a previous report (2).The results are of high quality and are well aligned with the vascular and metabolic 18 F-FDG patterns expected for key organ structures such as the brain, liver, kidneys, spleen, etc. ( 2).The uncertainties of V b , V d , distribution flow (K d ), and K i are generally higher for regions with larger magnitudes for the kinetic variable.This is perhaps related to the fact that these parameters, which are linear functions of the fitted voxel-level residue, ultimately scale with the magnitude of the time-course data.Mean transit time and extraction fractions deviate somewhat from this pattern.This is likely to be related to the fact that both the mean transit time and extraction fraction are defined in terms of ratios of the V d , K d , and K i variables and, as a result, do not necessarily scale with the scale of the voxel time course.The large blood vessels are seen to impact the structure of the MIP uncertainty for several parameters.The algorithms developed allow kinetic mapping, including the bootstrapping process, to be achieved in a timely fashion.On a single 3.2-GHz processor, the compute time for the NPRM kinetic analysis including the definition of the DGP is 140 min; each bootstrap replicate took 80 min.

Statistical Analysis
Representation of VOI Time-Course Data.The full time course as well as the time course over the first minute of data acquisition are shown in Figure 3. Average VOI time-course data are fit directly using the nonparametric and 2C models; averages of voxel-level fits are also provided.This gives a reference to the results reported previously (2).Although the 2C fitting of some VOIs is reasonable, for example, gray and white matter, there are clearly some VOIs where 2C modeling is substantially inferior (e.g., kidney, liver, bone, and bladder).The data fit achieved by the VOI averaging of the voxel-level nonparametric fit is quite good overall and especially over the first 1 min of acquisition.However, it is important to appreciate that almost half of the total number of frames occur in the first 80 s.For this example, over the first minute, differences between the VOI average of the voxelwise 2C fits and the fit of the 2C model to the mean of the VOI time-course data are quite pronounced.In contrast, differences between the corresponding nonparametric fits are much smaller.
Quantitative summaries of the nonparametric fitting of VOI timecourse data and comparisons with direct analysis of the mean VOI time-course data using nonparametric and 2C analysis are presented in Table 1.Although values from the weighted-residual-sums-of-squares fit for VOIs are similar based on the VOI average of voxel-level nonparametric fits or by direct fitting of the VOI time-course data, there is a marked increase in weighted-residual-sums-of-squares fits when the VOI time course is approximated using the best-fitting 2C model.VOI time-course fitting by the nonparametric model is consistently improved by averaging voxel-level nonparametric fits; the percent improvement is a modest 50%.VOI time-course fitting by the 2C model is substantially worse than the nonparametric fitting.The mean percent improvement here is almost 390%.
VOI Kinetics.VOI kinetics are reported in Table 2. Statistically significant deviations between the kinetics recovered by alternative methods are largely linked to early time-course parameters (Fig. 1), particularly for V b .Deviations between voxel-averaged parameters and values recovered from nonparametric and 2C analysis of the VOI time course are much smaller for nonparametric analysis than for 2C analysis.However, it is noteworthy that, for most VOIs, K i is quite similar in magnitude across all 3 analyses.This might be because flux is a late-time-course parameter (Fig. 1), and alternative methods fit the late time course quite similarly (Fig. 3).DGP Model. Figure 4 and Supplemental Figure 2 show an expected linear relation between the scale of the DGP and study dose; the linear correlation of 0.68 is highly significant.The axially averaged spatial scale of the DGP increases toward the top and bottom of the patient in the FOV.As expected, the increased scale is not just a function of the nominal sensitivity but is clearly impacted by patient-specific factors including the varying uptake, attenuation, and perhaps any impacts of small patient movements.The skewed nature of random fluctuations in the DGP model, which vary on the basis of the data coefficient of variation, are fully consistent with patterns for iteratively reconstructed PET data (10,21).The full width at half maximum of the autocorrelation functions in each direction is on the order of 2-3 mm.The coordinatewise autocorrelation functions show greater spatial persistence in the x (perpendicular to scanning bed) and z (axially) directions (Supplemental Fig. 3).This could align with involuntary patient movements during scanning.
SEs of VOI Kinetics.SEs of VOI kinetics (voxel-level nonparametric) are well approximated using a log-linear model that accounts for the VOI type, the VOI mean kinetics, and the residual weighted root mean square error of the voxel-level nonparametric fit of the VOI time course (Fig. 5; Supplemental Fig. 4).The overall correlation between the bootstrap-measured SE and the SE values predicted by log-linear modeling is 0.96, which is also quite high for individual kinetic parameters.

DISCUSSION
This work demonstrates the practicality of using image-domain bootstrapping for the construction of patient-specific uncertainty assessment in kinetics variables for voxel, VOI, and more complex derived quantities such as MIPs from a whole-body dynamic 18 F-FDG PET study.This development creates an opportunity to   incorporate uncertainty about a PET-derived kinetic biomarker that might be used to guide a clinical decision for a patient.This could be particularly helpful in cases where the biomarker value is close to a boundary between alternative treatment options.Bootstrap reliability depends both on the number of bootstrap simulations (N B ) used and on the accuracy of the representation of the data used in the DGP (20).Computational resources dictate the choice of N B .The results here are based on just an N B of 25, but for the data in Figure 2, a 4-fold increase in N B leads to little qualitative change in derived voxel-level SE (Supplemental Fig. 5).Table 1 clearly demonstrates the benefit of using a nonparametric methodology in the DGP.Relative to the well-established 2C 18 F-FDG model, substantial and highly significant improvements in data representation are achieved using the nonparametric approach.These benefits are mostly associated with the ability of the nonparametric technique to capture the highly resolved early time-course pattern of data from the current generation of PET scanners.The generally more modest deviations between nonparametric and 2C fits beyond the early time period, say after 1 min, suggest that the deficiencies in the 2C model may primarily relate to the lack of sophistication in the representation of the vascular components of blood-tissue exchange (22).The high temporal resolution of the scans here, as well as the use of a bolus injection, contributes to the ability to scrutinize the 2C model in ways that have likely not been possible in the past.The VOIs here are large and heterogeneous-far from the assumption of homogeneous well-mixed compartments that underly the 2C model.However, it is notable that our previous work (23) reported significant discrepancies between 2C and nonparametric representation of dynamic 18 F-FDG brain data in healthy subjects using much smaller and highly homogeneous  2 for gray and white matter, the discrepancies primarily impact the accuracy of the initial phase of the 18 F-FDG tissue residue-V b especially-but have much less impact on several other variables including flux and V d .However, statistically significant differences between voxel nonparametric and VOI 2C parameters do not imply that parameters are unrelated.For example, Figure 6 and Supplemental Figure 6 show pairwise plots and summary correlations for the 18 F-FDG metabolic rate (MR) flux scaled by the plasma glucose in Equation 4. The strong linear dependence in Figure 6 emphasizes the importance of differentiating statistical and practical significance.Calculated K i based on nonparametric or 2C analysis would likely yield similarly effective diagnostic values.Indeed, it is well appreciated that even simpler assessments of 18 F-FDG flux by Patlak analysis and SUV are also highly effective.
The nonparametric technique here uses a linear basis, but the structure and number of elements involved are adapted to the full 4-dimensional dynamic data and guided by cross validation to prevent overfitting (10).The accuracy and stability of a kinetic mapping procedure are best evaluated numerically, which was reported previously (24)-studies based on a 2-min constant infusion injection of 18 F-FDG and a temporal sampling protocol in which the shortest time frames were 20 s in duration, providing mean-square-error performance characteristics of NPRM and 2C kinetic mapping of 18 F-FDG PET data as a function of the study dose and as a function of whether the underlying ground truth is governed by a compartmental model or not.In this study, the accuracy of the flux is largely unaffected by whether a 2C or an NPRM mapping technique is used.Across other kinetic variables, when the ground truth is noncompartmental, the NPRM approach is much better.Remarkably, when the ground truth is a 2C model, the NPRM continues to outperform the 2C approach, especially for variables such as V b and V d .Further study of the mean-square-error performance would clearly be useful, particularly in settings where the ground truth, study protocol, and scanning methods are similar to those encountered with the current generation of whole-body 18 F-FDG PET studies.
VOI values of 3 variables, FDG metabolic rate (MR FDG ), distribution volume (DV), and vascular blood flow (BF), are compared with literature reports.Each variable is directly obtained by simple scaling of our summary kinetic values: K i , V d , and V b .
Here, m glc is the plasma glucose concentration and t* is the value used to define the vascular component in the decomposition of the Meier-Zierler residue in Figure 1.In a cancer setting, 18 F-FDG MR is by far the most clinically important of these variables.Note that we do not try to use 18  interest in deriving potentially useful additional diagnostic information related to tissue vascularity from 18 F-FDG (1,25,26).There is no intention of questioning PET 15 O-H 2 O as the gold standard for V b determination.Our V b formula is an application of the central volume theorem (14) based on an assumed mean transit time in the vasculature of t*/2 (here, 7.5 s) for the collection of tracer atoms whose tissue transit time in the local voxel is less than 15 s.Table 3 compares the VOI averages of 3 variables to those in the literature.For 18 F-FDG MR and V d , the values are seen to be in the range reported using 2C and Patlak analyses (27).V b values are compared with those in reports based on PET 15 O-H 2 O and dynamic susceptibility contrast MR techniques.The results for the NPRM approach are remarkably similar to those in the literature, particularly given that the study group here is older and unhealthy (28).Further examination of the V b variable could be merited.Viability of conducting PET 15 O-H 2 O on this scanner was previously demonstrated (29).Note that some of the deviation in Table 3 may be related to scaling differences between the use of whole-blood activity as an AIF (like ours) and other analyses that used the arterial plasma activity time course as an AIF.
Although our focus has been on parameters that have traditionally been used to quantify 18 F-FDG PET dynamics, the nonparametric technique provides a possibility to also evaluate a summary of the arrival pattern of 18 F-FDG at the voxel level.A sample amplitude-weighted average of voxel-level basis element delay as shown in Equation 1 is shown in Figure 7.There is early arrival of the signal to the lung and much more delayed arrival to the bladder and more peripheral regions (1).More detailed consideration of the 18 F-FDG arrival pattern may be worthwhile.PERTINENT FINDINGS: NRPM analysis together with imagedomain bootstrapping is a suitable methodology for mapping kinetics.

IMPLICATIONS FOR PATIENT CARE:
The ability to derive uncertainties in complex kinetic biomarkers could enhance patientspecific decision-making for guiding treatment of cancer patients.
F-FDG (with activity of $3 MBq/kg of patient weight) to the left or right arm, followed by flushing with 50 mL of saline solution.The plasma glucose level was measured for each patient.Emission data were acquired for 65 min and binned into 62 contiguous time frames with durations of 2 3 10 s, 30 3 2 s, 4 3 10 s, 8 3 30 s, 4 3 60 s, 5 3 120 s, and 9 3 300 s.Images were reconstructed with a voxel size of 1.65 3 1.65 3 1.65 mm 3 .Low-dose CT scans (voltage, 120 kV; tube current, 25 mA; CARE Dose4D and CARE kV [Siemens]) were acquired as part of the examinations.The CT images were reconstructed with a voxel size of 1.52 3 1.52 3 1.65 mm 3 .
Time-Course Data.Mean VOI timecourse data are compared with the corresponding mean VOIs of the fitted voxel-level time courses, ẑðt j Þ, in Equation 2. Mean VOI timecourse data are also analyzed using the nonparametric model and the Huang-Sokoloff 2C model including a fractional V b and delay of the

FIGURE 1 .
FIGURE 1. Meier-Zierler tissue residue (R) with decomposition into vascular (R b ), in-distribution (R d ), and extracted (R e ) components.Decomposition was used to define indicated metabolic parameters.MTT 5 mean transit time; Ext 5 extraction fraction.

FIGURE 2 .
FIGURE 2. MIP maps of NPRM kinetic parameters and associated SE.SEs are based on SD of MIP results for each of 25 bootstrap replications.Top row shows CT images for selected cross sections through volume and PET MIP maps at indicated times.K 5 9 basis elements were determined for data by NPRM methodology.MTT 5 mean transit time; Ext 5 extraction fraction.

FIGURE 3 .
FIGURE 3. Results of alternative fitting of VOI data used in Figure 2. Data are points, and line colors correspond to methods used.Full time course is on left; first minute is on right.GM 5 gray matter; NP 5 nonparametric; WM 5 white matter.

CONCLUSION
NPRM kinetic analysis together with bootstrap assessment of uncertainty is practically feasible in the context of large-scale longaxial-FOV18 F-FDG PET data.This provides an ability to incorporate patient-specific uncertainty measures of kinetic biomarkers recovered from dynamic PET to support clinical decisions.DISCLOSUREThis research is supported by Science Foundation Ireland grant PI-11/1027 and by the National Cancer Institute USA grant R33-CA225310.Kuangyu Shi and Axel Rominger are funded by Siemens Healthineers and Novartis.Hasan Sari is an employee of Siemens Healthineers.No other potential conflict of interest relevant to this article was reported.KEY POINTSQUESTION: Is it feasible to map kinetics together with uncertainty in long-axial-FOV dynamic18 F-FDG PET studies?

FIGURE 7 .
FIGURE 7. (A) Delay image corresponding to coronal CT slice in Figure 2; values are amplitude-weighted delay values, fd1D k , k51, 2, . . ., Kg, in Equation 1. Data are centered so that mean delay in descending aorta is 0. (B) Box plots of distribution of mapped delay values in seconds by VOI.GM 5 gray matter; WM 5 white matter.

TABLE 2
VOI Kinetics Recovered Using Different Methodologies VOIs.Similar to what is reported in Table