Visual Abstract
Abstract
The purpose of this study was to examine a nonparametric approach to mapping kinetic parameters and their uncertainties with data from the emerging generation of dynamic whole-body PET/CT scanners. Methods: Dynamic PET 18F-FDG data from a set of 24 cancer patients studied on a long-axial-field-of-view PET/CT scanner were considered. Kinetics were mapped using a nonparametric residue mapping (NPRM) technique. Uncertainties were evaluated using an image-based bootstrapping methodology. Kinetics and bootstrap-derived uncertainties are reported for voxels, maximum-intensity projections, and volumes of interest (VOIs) corresponding to several key organs and lesions. Comparisons between NPRM and standard 2-compartment (2C) modeling of VOI kinetics are carefully examined. Results: NPRM-generated kinetic maps were of good quality and well aligned with vascular and metabolic 18F-FDG patterns, reasonable for the range of VOIs considered. On a single 3.2-GHz processor, the specification of the bootstrapping model took 140 min; individual bootstrap replicates required 80 min each. VOI time-course data were much more accurately represented, particularly in the early time course, by NPRM than by 2C modeling constructs, and improvements in fit were statistically highly significant. Although 18F-FDG flux values evaluated by NPRM and 2C modeling were generally similar, significant deviations between vascular blood and distribution volume estimates were found. The bootstrap enables the assessment of quite complex summaries of mapped kinetics. This is illustrated with maximum-intensity maps of kinetics and their uncertainties. Conclusion: NPRM kinetics combined with image-domain bootstrapping is practical with large whole-body dynamic 18F-FDG datasets. The information provided by bootstrapping could support more sophisticated uses of PET biomarkers used in clinical decision-making for the individual patient.
High-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 18F-FDG using the well-established Huang–Sokoloff 2-compartment (2C) modeling framework (1–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 18F-FDG and their distribution across normal and cancerous tissues has evolved in the years since the Huang–Sokoloff construct was proposed (4–7). The temporal and spatial resolutions of emerging scanners have transformed the ability to objectively assess the accuracy of the 2C framework to represent 18F-FDG time-course data across the diverse tissues encountered in the human body. In this context, the assessment of 18F-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 decision-making for individual patients that is based on complex nonlinear radiomic metrics derived from a kinetic map.
The volume of data produced by a dynamic 18F-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 18F-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 18F-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 18F-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 × 10 s, 30 × 2 s, 4 × 10 s, 8 × 30 s, 4 × 60 s, 5 × 120 s, and 9 × 300 s. Images were reconstructed with a voxel size of 1.65 × 1.65 × 1.65 mm3. 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 × 1.52 × 1.65 mm3.
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, Cp, and the regional tissue residue function. Kinetic parameters are defined in terms of this residue (Fig. 1). Large-vessel vascular blood and distribution volumes (Vb and Vd, respectively) are evaluated as areas under the tissue residue. The apparent rate of retention or flux (Ki) 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 18F-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-of-basis 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 for k of 1, 2,…, K. Here, Rk is the basis element residue and Δ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, , is expressed as Eq. 1Here, δ and (α1, α2,…, αK) are the unknown voxel-level delay and basis-amplitude 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 α 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, , are used to construct an image-domain data generation process (DGP) for bootstrapping. The DGP generates data according to Eq. 2where and the simulated error process, ϵ*, 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 Time-Course Data
Mean VOI time-course data are compared with the corresponding mean VOIs of the fitted voxel-level time courses, , in Equation 2. Mean VOI time-course data are also analyzed using the nonparametric model and the Huang–Sokoloff 2C model including a fractional Vb and delay of the AIF. To facilitate fitting, a wide range of delays of ±5 min is allowed in the NPRM and 2C analysis of the VOI time-course data. The Broyden–Fletcher–Goldfarb–Shanno algorithm as implemented in the optim function in R (R Foundation for Statistical Computing) is used for optimization of the 2C model; more details are shown in the supplemental materials and Supplemental Figure 1.
Results of alternative analyses for a sample case are presented graphically. Formal comparisons are focused on the weighted-residual-sums-of-squares misfit achieved with alternative analyses. The mean relative difference between alternative representations of VOI time-course data and the associated SD is evaluated for each VOI type considered.
VOI Kinetics
Means and SDs of VOI-averaged NPRM kinetics are reported for each VOI type. Kinetics based on nonparametric and 2C analyses of VOI mean time-course data are similarly summarized. Deviations between alternative VOI kinetic values are summarized, and their statistical significance is evaluated using the paired Wilcoxon test.
DGP Model
The bootstrap DGP is expressed in more detail as Eq. 3where the random errors, , are in units of SD and is an overall scale of the model error. In Equation 3, the factors and 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 () should also scale with dose; this is examined graphically. The overall axial pattern variation is described by the scale factor . 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-of-squares 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.
RESULTS
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 18F-FDG patterns expected for key organ structures such as the brain, liver, kidneys, spleen, etc. (2). The uncertainties of Vb, Vd, distribution flow (Kd), and Ki 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 Vd, Kd, and Ki 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 time-course 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 Vb. 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, Ki 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 18F-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 (NB) used and on the accuracy of the representation of the data used in the DGP (20). Computational resources dictate the choice of NB. The results here are based on just an NB of 25, but for the data in Figure 2, a 4-fold increase in NB 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 18F-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 18F-FDG brain data in healthy subjects using much smaller and highly homogeneous VOIs. Similar to what is reported in Table 2 for gray and white matter, the discrepancies primarily impact the accuracy of the initial phase of the 18F-FDG tissue residue—Vb especially—but have much less impact on several other variables including flux and Vd. 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 18F-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 Ki based on nonparametric or 2C analysis would likely yield similarly effective diagnostic values. Indeed, it is well appreciated that even simpler assessments of 18F-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 18F-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 18F-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 Vb and Vd. 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 18F-FDG PET studies.
VOI values of 3 variables, FDG metabolic rate (MRFDG), 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: Ki, Vd, and Vb. Eq. 4Here, 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, 18F-FDG MR is by far the most clinically important of these variables. Note that we do not try to use 18F-FDG as a means to evaluate the glucose MR, as described in Phelps et al. (17). Barrio et al. (6) expressed considerable doubt on the ability to do this in the context of cancer applications. Consideration of the Vb variable is motivated by interest in deriving potentially useful additional diagnostic information related to tissue vascularity from 18F-FDG (1,25,26). There is no intention of questioning PET 15O-H2O as the gold standard for Vb determination. Our Vb 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 18F-FDG MR and Vd, the values are seen to be in the range reported using 2C and Patlak analyses (27). Vb values are compared with those in reports based on PET 15O-H2O 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 Vb variable could be merited. Viability of conducting PET 15O-H2O 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 18F-FDG PET dynamics, the nonparametric technique provides a possibility to also evaluate a summary of the arrival pattern of 18F-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 18F-FDG arrival pattern may be worthwhile.
CONCLUSION
NPRM kinetic analysis together with bootstrap assessment of uncertainty is practically feasible in the context of large-scale long-axial-FOV 18F-FDG PET data. This provides an ability to incorporate patient-specific uncertainty measures of kinetic biomarkers recovered from dynamic PET to support clinical decisions.
DISCLOSURE
This 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 POINTS
QUESTION: Is it feasible to map kinetics together with uncertainty in long-axial-FOV dynamic 18F-FDG PET studies?
PERTINENT FINDINGS: NRPM analysis together with image-domain bootstrapping is a suitable methodology for mapping kinetics.
IMPLICATIONS FOR PATIENT CARE: The ability to derive uncertainties in complex kinetic biomarkers could enhance patient-specific decision-making for guiding treatment of cancer patients.
Footnotes
Published online Apr. 11, 2024.
- © 2024 by the Society of Nuclear Medicine and Molecular Imaging.
Immediate Open Access: Creative Commons Attribution 4.0 International License (CC BY) allows users to share and adapt with attribution, excluding materials credited to previous publications. License: https://creativecommons.org/licenses/by/4.0/. Details: http://jnm.snmjournals.org/site/misc/permission.xhtml.
REFERENCES
- Received for publication September 17, 2023.
- Revision received February 13, 2024.