Improved kinetic analysis of dynamic PET data with optimized HYPR-LR

INTRODUCTION

HighlY constrained backPRojection (HYPR) is a family of image reconstruction and post processing algorithms that have made a large impact on magnetic resonance angiography (MRA), allowing for undersampling factors on the order of several hundred fold and dramatic signal-to-noise ratio (SNR) improvements in dynamic datasets.^{1, 2, 3} The central idea of HYPR is to estimate individual frames in a dynamic study by weighting a temporally summed composite image. In the ideal case, individual frames will take on the noise properties of the composite image. It is the potential improvement in the SNR that makes HYPR an attractive tool for other imaging modalities, and we have recently applied the post-processing version, HYPR-LR (for local reconstruction), to positron emission tomography (PET) with promising results.⁴

One of the primary motivations for applying HYPR-LR to PET data is to improve the kinetic analysis methods used to analyze dynamic PET data, all of which are sensitive to noise to some degree. In our previous work, we demonstrated this with simple parameters estimated from phantom studies (e.g., nonlinear least squares fitting to estimate radionuclide decay constants).⁴ In these studies, we summed all the frames of the study to form the composite image, as this will provide the greatest enhancement in the SNR. However, it is well understood that differences between the composite image and the frame of interest can distort the temporal signals of neighboring regions with different temporal behaviors.¹ In the context of quantitative kinetic analysis, the bias caused by HYPR-LR processing may have an impact on the analysis, which should be accounted for and controlled as best as possible. Traditionally, bias in HYPR processing has been minimized by using a sliding composite window, so that composite images are more similar to their corresponding frames.¹ While such an approach may potentially be effective for dynamic PET data, it is arbitrary, and the optimal implementation of HYPR-LR will be dependent on the specific imaging task. In addition, as the noise reduction provided by HYPR-LR will be dependent on the number of counts in the composite (i.e., its temporal duration) relative to the frame of interest,^{1, 2} it will provide the greatest benefit when the composites are made using as long an integration time as possible without introducing any bias. If composite images are shorter than this, they will provide no increased benefit in terms of reduced bias, and if they are made larger than this the additional noise improvements gained will come at the cost of bias.

In this work, we aim to develop an optimized implementation of HYPR-LR that can be tailored to individual PET studies such that it will provide the greatest reduction in variance while minimizing bias introduced by the algorithm, and to closely evaluate the performance of HYPR-LR processing in the context of quantitative kinetic analysis. We attempt to fulfill these aims using a combination of simulated and real human data of a tracer of interest to us, [¹¹C]-Pittsburgh compound-B (PIB). [¹¹C]-PIB is amenable to analysis using a number of techniques, including Logan graphical analysis,⁵ a data driven method with a well-known noise-dependent bias,^{6, 7, 8} and simplified reference tissue model (SRTM) methods,⁹ which are model based and not sensitive to bias.^{10, 11, 12} The basis function implementation of SRTM, so called receptor parametric mapping (RPM or RPM2), is regarded as more robust than the Logan method with respect to both variance and bias.⁹ We examine the effects of our proposed implementation of HYPR-LR processing on both the bias and precision of nondisplaceable binding potential (BP_ND) estimates obtained with both Logan graphical analysis and RPM2. We do this in the context of parametric image generation, as single voxel analysis will be the most sensitive to noise. Simulated data are used to illustrate the trade-off between bias and variance in absolute terms where the truth is known, and human data are used to illustrate the algorithm’s properties in the context of a real imaging task. While, we focus on two analysis methods for [¹¹C]-PIB data here, we also intend to provide a framework for determining the optimal implementation of HYPR-LR for a given imaging task that can be generalized to other tracers and analysis methods.

THEORY

All of the formulations of HYPR make use of a temporally summed composite image in a dynamic set of images to provide a low-noise estimate of the true image at an individual time frame. This composite is weighted by a low resolution spatial comparison of each frame and the composite image. In the case of HYPR-LR, the spatial comparison is achieved by convolving the frame of interest and the composite image with a filtering kernel

where H_t is the HYPR-LR estimate of an individual frame at time t, C_t is the composite image used for time t, W_t is the weighting image, I_t is an initial reconstruction of the frame, F is the filtering kernel, and ⊗ represents the convolution operation. When images are sparse, as they are in angiography, so long as two objects with different temporal behaviors do not overlap in the filtering process, they will be perfectly reproduced.² Unfortunately, in non sparse images, such as those in PET, regions with different temporal behavior will invariably overlap in the filtering process, potentially creating a bias. This bias will manifest itself at high-frequency boundaries since the weighting image is formed by a low-frequency comparison.

The only way for the HYPR-LR estimate of a PET frame to be perfectly accurate is if the contrast between objects present in an individual frame is the same as the contrast between those objects in the frame’s composite image. Failing this, errors at high-frequency boundaries will depend on the difference in contrast between neighboring regions and the size of the kernel used (Fig. 1).

Quantitatively accurate application of HYPR-LR to PET data may still be possible because of the predictable behavior of PET tracers. This predictable behavior can be used to form composite images that more closely fulfill our criteria of having the same contrast between objects as their corresponding frames, and our proposed optimized implementation of HYPR-LR for PET data attempts to achieve this. We divide a typical tracer’s behavior into three phases: uptake, specific retention, and equilibrium. Immediately after a tracer’s injection, during the brief uptake phase, the activity in all regions will be increasing at a relatively constant rate. As a tracer is retained in regions of specific avidity, these regions begin to distinguish themselves from the background in a predictable manner. Finally, the tracer will reach a time of equilibrium, known as t*, the attainment of which is required for application of graphical analysis methods and after which the image changes very little.^{6, 7, 13, 14} Figure 2a illustrates these behaviors for [¹¹C]-PIB for different regions of the brain. A composite image could thus be formed for the uptake and equilibrium phases by summing all the frames in these phases. A sliding composite could be used during the retention phase such that the change in contrast is nearly linear for the frames used, creating a composite very similar to the individual frame being processed.

The formation of the composite images can thus be described as:

Note that t denotes frame, not time. The phases of the study are demarcated here by t_uptake, the frame at which the uptake phase ends, t*, the frame at which the steady state begins, and t_max, the final frame of the study. The parameter α determines the number of frames to be used in the sliding composite during the specific retention phase, and Δ_t is the duration of frame t. It is also important to note that data should not be decay corrected, so that each frame is weighted based on counts in its contribution to the composite image. Decay correction is applied after HYPR-LR processing.

The temporal parameters in Eq. 2a, 2b, 2c must be determined from the data to be analyzed. We propose doing this on a region of interest (ROI) level. Specifically, ROIs should be drawn over structures in the image, and the ratios of neighboring structures must be examined. More emphasis can be put on structures of relevance. We illustrate this here with the parietal cortex and its surrounding structures in the case of a real [¹¹C]-PIB dataset [Fig. 2b]. The parameters t_uptake and t* can be determined by identifying the end of the flat uptake phase and the start of the flat equilibrium phase, respectively. α is more arbitrary, and we propose selecting a value such that the composite image is 10–20 min (5–7 frames) in duration for a typical PET imaging study. The frames at the beginning and end of the retention phase will largely determine the size of α, and α should be made as large as possible while keeping the contrast in the composite images for these frames as close as possible to the contrast in the frames themselves. This process could be performed for individual patients, or for a general population for whom a tracer’s temporal behavior is relatively consistent.

The size of the kernel, F, used in the formation of the weighting image is the other parameter which must be determined. The size of the kernel will determine the degree to which regions with disparate temporal behaviors overlap and contribute to bias, and also determine the spatial frequencies of individual frames that effectively take on the noise properties of their composite images. If our composite scheme is implemented correctly and the contrast between objects in individual frames and their respective composites is nearly the same, overlap between regions should contribute little to temporal errors. A larger kernel should thus be desirable, as it will suppress noise at a wide range of spatial frequencies. In this work, the largest kernel used was a 3D Gaussian with a FWHM approximately double the resolution of the scanner, 9 mm in this case. If the composite images do not match their respective frames, a relatively small kernel should be used to limit bias. The smallest kernel use in this work was a 3D Gaussian with a FWHM of 3 mm, close to the thickness of the cortex, the primary object of interest for [¹¹C]-PIB. We evaluate the effects of different kernel sizes as outlined in our methods.

Our proposed method of forming multiple composite images will itself be limited in a number of respects. For one, it will likely not be possible to perfectly fulfill our criterion of forming composite images having the same contrast between objects as in their respective frames. However, so long as there is not substantial bias introduced into multiple frames, having a few frames slightly biased by imperfect composite images should not affect kinetic analysis a great deal. Second, not all dynamic PET studies will be amenable to our particular method of forming composites described in Eq. 2a, 2b, 2c, for example tracers with a washout phase. Finally, using shorter composite images will increase the variance in individual frames and parametric images. Nevertheless, the ratio of the activities of neighboring regions should at least be examined to provide insight into what time-series data will produce composite images those objects in the respective frames. We believe this will allow for the maximum benefit in terms of variance reduction while minimizing bias for a given PET tracer and imaging task.

METHODS

RESULTS

DISCUSSION

HYPR-LR is a promising denoising technique for PET, and here we have attempted to illustrate its ability to improve the variance and bias of parametric images derived from kinetic analysis techniques. In order to maximize the improvement in variance and bias while minimizing the error introduced by HYPR-LR processing, we have proposed a method of forming multiple, time dependent, composite images (HYPR-LR-MC) that minimizes the difference of the contrast between objects in each frame of the study and the frame’s respective composite image. Such an approach will provide less of an improvement in noise and introduce less bias than using all the frames in the study to form the composite image. It is also less arbitrary than the traditional approach of using a sliding composite window.¹

Bias introduced by HYPR-LR processing will be dependent on differences in contrast that do exist between frames of interest and their respective composites. In the case of the simulated and human [¹¹C]-PIB presented here, HYPR-LR-FC processing distorts TACs [Fig. 3g] and thus inappropriately increases the BP_ND values obtained with both Logan graphical analysis and RPM2 in regions of high uptake where they border regions of lower uptake, and decreases the BP_ND values of the surrounding voxels (Fig. 3). This bias can also be seen in the overestimated BP_ND values of the frontal and parietal cortices in the parametric images generated from the noisy simulated data processed with HYPR-LR-FC (Fig. 4), and in the BP_ND values of these same regions in the parametric images of the human datasets (Fig. 6). HYPR-LR-MC processing virtually eliminates the bias caused by HYPR-LR (Fig. 3), as it almost completely removes the noise-dependent Logan bias without overestimating BP_ND values, and introduces little bias into parameters obtained with RPM2 (Figs. 46). Bias can also be limited with HYPR-LR-FC if a smaller kernel size is used. In this work, the smallest kernel used was a Gaussian with a FWHM of 3 × 3 × 3 mm³, only slightly larger than the thickness of the cerebral cortex, our primary object of interest. For the simulated data studied here, HYPR-LR-FC with this size kernel performed particularly well when the data were analyzed with RPM2, only introducing slightly more bias than HYPR-LR-MC with a large 9 mm FWHM filtering kernel (5.39% versus 3.73%). HYPR-LR-FC with the smaller kernel size did not provide as much of a benefit to the data analyzed with the Logan graphical method, as BP_ND values are still relatively biased (Fig. 4).

Both methods of implementing HYPR-LR demonstrate an ability to reduce variance in the parametric images presented here. HYPR-LR-FC processing with a larger filtering kernel in particular substantially reduces the variance in parametric images generated with both the Logan graphical method and RPM2 (Figs. 457). This is not surprising, as each individual frame in a study is given the noise properties of the fully summed composite image. Processing with our proposed HYPR-LR-MC method reduces variance more substantially in parametric images generated with the Logan graphical method (Figs. 457). HYPR-LR-MC also reduces variance in the parametric images created with RPM2, but more moderately, particularly in the simulated data (Fig. 4). The results from the real data are more promising and show that HYPR-LR-MC processing visually reduces the variance seen in the example parametric image (Fig. 5), and increases the correlation between BP_ND values obtained with the reference Logan method and RPM2 (Fig. 7 and Table TABLE II.). As expected, a larger kernel size results in better variance reduction with the HYPR-LR-MC algorithm while introducing little or no additional bias. A larger kernel likewise improves the variance reduction with the HYPR-LR-FC algorithm, but more bias is introduced as a result (Fig. 4). Using HYPR-LR-FC with a small filtering kernel may therefore be desirable in some applications. For example, when the simulated data are analyzed with RPM2, HYPR-LR-FC with a smaller kernel still achieves good variance reduction (from 16.6% to 8.7%) while introducing little bias, as discussed above.

In this work, we have compared HYPR-LR processing to simple spatial smoothing. While many other denoising methods have been developed, spatial smoothing, either within image reconstruction or afterward, remains one of the most common ways to control noise in PET data. In addition, the processing time required for HYPR-LR will not be substantially longer than the time required to spatially smooth each frame, as processing time will largely be determined by the number of convolution operations that must be performed. In this work, the number of convolutions required for HYPR-LR-MC processing is 1.6 times greater than spatially smoothing each frame. HYPR-LR processing compares very favorably to smoothing when the data are analyzed with the reference Logan graphical method. HYPR-LR-MC and HYPR-LR-FC control both bias and variance very well, whereas smoothing with the 3 × 3 × 3 mm³ FWHM Gaussian only slightly improves the bias and variance, and smoothing with the 6 × 6 × 6 mm³ Gaussian controls variance very well but introduces a substantial bias of its own because of excessive blurring (Figs. 456). The improvements offered by HYPR-LR are less dramatic compared to simple smoothing when parametric images are created with RPM2. In the simulated data, smoothing with a 3 × 3 × 3 mm³ FWHM Gaussian improves variance more than HYPR-LR-MC (10.5% versus 14.0%), though the smoothed BP_ND values are more biased relative to the unprocessed data as there is a loss of spatial resolution (Fig. 4). HYPR-LR-FC using a 3 × 3 × 3 mm³ Gaussian kernel does provide more variance reduction than simply smoothing with a 3 × 3 × 3 mm³ Gaussian (8.7% versus 10.5%) for a comparable amount of bias relative to the original data (Fig. 4), although the bias caused by HYPR-LR-FC is not due to any loss of spatial resolution. The human data are consistent with this, demonstrating that simple smoothing with a 3 × 3 × 3 mm³ FWHM Gaussian reduces the variance in the RPM2 images to a greater degree than HYPR-LR-MC while only modestly blurring the data (Figs. 56). The lower correlation between the RPM2 and Logan BP_ND values following smoothing with the 3 × 3 × 3 mm³ FWHM Gaussian is due to the high variance in the Logan images (Fig. 7). HYPR-LR processing can also be done in addition to simple smoothing to provide further reductions in variance in the parametric images without introducing any additional loss of spatial resolution (Figs. 45).

In the future, HYPR-LR denoising must be considered more fully in the context of other denoising methods, for example wavelet denoising and iterative image reconstruction algorithms, including the numerous proposed approaches to four-dimensional PET reconstruction, which have previously been explored as a means of improving kinetic analysis and parametric image generation.^{19, 20, 21, 22, 23, 24, 25} As HYPR-LR uses temporally integrated data to reduce noise, comparing it to denoising methods that likewise utilize the time domain, such 4D reconstructions, will be particularly important. Iterative reconstructions, wavelet denoising, and HYPR-LR have all demonstrated an ability to substantially reduce noise, but each also has drawbacks. A full comparison between HYPR-LR and these other denoising processes, explicitly examining the pros and cons of each, is beyond the scope of the current study and will likely depend on the specific imaging task. Furthermore, as demonstrated here with spatial smoothing, as a post processing technique that is fast and relatively simple to implement, HYPR-LR could easily complement these other denoising techniques. This is particularly relevant for performing HYPR-LR on OSEM reconstructions, which are now routinely done on both clinical and research scanners and provide some noise control themselves.

In this work, we have shown that following HYPR-LR processing, Logan and RPM2 analysis perform quite comparably. After denoising the simulated data with HYPR-LR-MC and HYPR-LR-FC processing with larger filtering kernels, both the coefficient of variation and the bias are similar between the Logan method and RPM2 (Fig. 4). If anything, after HYPR-LR processing the greater bias is seen with RPM2 with little additional benefit in variance reduction. The human data confirm this, showing little bias in the BP_ND values obtained with either Logan analysis or RPM2 after HYPR-LR-MC processing, and greater bias in the RPM2 data following HYPR-LR-FC processing (Fig. 6). The level of variance in the parametric images generated from the human data with the two analysis methods is also comparable following HYPR-LR processing (Figs. 57).

In addition, we found the performance of RPM2 to be more variable and sensitive to user selected parameters. In particular, we selected a minimum k₂ slightly greater than that predicted by some of the non-linear SRTM fits to the ROI TACs, and greater than the value previously reported for [¹¹C]-PIB analysis.⁹ We also see a slight positive bias in the BP_ND values obtained with RPM2 in the simulations (Fig. 4). This might be because the assumption of a simplified reference tissue model does not adequately describe the kinetics of [¹¹C]-PIB in this case.⁵

While our objective here is not to determine the best way to analyze [¹¹C]-PIB data, the increased reliability and greater robustness seen with Logan analysis following HYPR-LR processing may prove valuable. The Logan method is easy to implement, does not assume any particular model, and does not require imposing limits on the data, for example the range of k₂ in RPM2. HYPR-LR processing does also improve the parametric images created with RPM2, albeit not as dramatically relative to simple smoothing.

Presumably, HYPR-LR processing should also improve other data driven methods, such as the multi-linear regression methods that are not susceptible to a noise-sensitive bias but do exhibit greater variance than the Logan graphical method.^{26, 27} Likewise, it should improve other model based methods like compartmental analysis.

The results presented here demonstrate that both the proposed method for forming multiple time dependent composites and simply using all the frames in the composite formation may have value in different contexts. HYPR-LR-MC introduces minimal bias but the bias introduced by HYPR-LR-FC is not that great in the data studied here, particularly when a smaller filtering kernel is used. And while HYPR-LR-MC certainly reduces variance, HYPR-LR-FC does so to a greater degree. When the focus is on relatively large regions of high uptake in [¹¹C]-PIB data, HYPR-LR-FC appears to perform relatively well. However, we still urge caution in using all frames of the study for forming the composite image. In the case of [¹¹C]-PIB, it appears as if the bias caused by HYPR-LR-FC results in greater contrast between structures in the brain. While such an image may be appealing to look at, it may be a misleading result. Bias caused by HYPR-LR will also likely be greater in studies that have greater contrast between areas of interest and their surrounding background, for example [¹¹C]-raclopride or any of a number of tracers used to study tumor biology. It is clear that using HYPR-LR with composite images that have been formed in a way that accounts for the kinetic behavior of the tracer being studied can certainly provide an improvement in the variance of kinetic parameters while introducing very little bias.

While we have focused on the application of HYPR-LR processing to [¹¹C]-PIB data, a tracer of interest in the neuroscience community, the approach that we have presented here should be generalizable to other tracers and other applications outside the brain. Indeed, the fairly detailed structure of the brain, and of the cortex in particular, provides a good means of testing HYPR-LR processing. HYPR-LR processing may prove particularly valuable for providing more detailed physiologic information on a smaller anatomical scale for cardiac and oncology applications of dynamic PET imaging. Not all tracers will follow the exact uptake pattern we have described, and the optimal composite scheme, including simply using a fully summed composite, will depend on the application. Nevertheless, examining the ratio of the activities of neighboring regions over time will still provide valuable insight as to how HYPR-LR can be applied to maximize variance reduction while minimizing introduced bias.

CONCLUSION

HYPR-LR is a promising denoising technique for a number of medical imaging modalities, and we have previously demonstrated its potential for denoising dynamic PET data. In this work, we have shown that HYPR-LR processing can improve kinetic analysis techniques used for processing dynamic PET data, and we have introduced a method for forming multiple time-dependent composite images that minimizes the bias the HYPR-LR algorithm causes while maximizing the improvement in variance it provides. This implementation of HYPR-LR could thus improve the kinetic analysis of dynamic PET data without sacrificing accuracy. In addition, our comparison between Logan graphical analysis and RPM2 before and after HYPR-LR processing provides a framework for testing the validity of HYPR-LR processing in the context of a given tracer and two methods of analysis. We envision HYPR-LR being particularly valuable in PET applications that suffer from high noise, such as PET scans requiring high spatial or temporal resolution, vulnerable patient populations who require less radiation dose, tracers used to screen large populations, and tracers utilizing unique but dosimetry limited radionuclides such as ¹²⁴I and ⁶⁴Cu. HYPR-LR is a simple denoising tool, and our proposed method for its quantitatively accurate implementation could easily be implemented for different tracer behaviors, or for individual studies.