Optimizing PiB-PET SUVR Change-Over-Time Measurement by a large-scale analysis of Longitudinal Reliability, Plausibility, Separability, and Correlation with MMSE

3.1 Question 1: Is Partial Volume Correction Helpful?

Partial Volume Correction (PVC) attempts to correct for the effect of relative uptake levels in multiple tissue types and/or CSF within each PET voxel. Generally, this uses segmentations from a T1-weighted Magnetic Resonance Image (MRI) in combination with the known resolution of the PET camera used to estimate the fraction of each material expected to lie in each voxel or ROI. In this work we evaluate three options: no PVC; 2-class PVC (PVC2), which attempts to correct for CSF in GM/WM locations; and 3-class PVC (PVC3), which attempts to correct for both CSF and WM in GM locations. Briefly, 2-class PVC uses a binary map of GM+WM locations from MRI segmentation and blurs this by the known resolution of the PET camera (approximately 8mm³ for the scanners in our study (Joshi et al., 2009)) to create an approximate fraction of tissue expected in each voxel. The raw PET scan is then divided by this tissue fraction in each voxel. This algorithm depends on the accuracy of the MRI segmentation and registration between PET and MRI (Meltzer et al., 1990). In 3-class PVC, MRI segmentation is instead used to create separate maps of voxels containing GM and those containing WM. First, an atlas is used to locate the centrum semiovale region and calculate the mean PET signal in this ROI, which is then assumed to reflect the true WM signal in all WM voxels. This mean value is assigned to all voxels segmented as WM, then blurred with the approximate PSF and subtracted from the raw PET scan. This step is designed to remove the signal from each voxel that is due to WM. Next, the image follows the same steps as in two-class PVC, which then attempts to correct further partial voluming by CSF (Müller-Gärtner et al., 1992). Because this method is designed to remove WM signal from each voxel, it is not appropriate when WM voxels are included in either reference or target regions for SUVR calculation. For this reason, we omit pipeline combinations that would use 3-class PVC with WM-containing target or reference regions.

Both algorithms were implemented in-house in C++ using the Insight Toolkit (www.itk.org), and were applied using the MRI scans that corresponded to each subject’s PET scan at each timepoint, in order to allow them to correct for differing amounts of partial volume over time, which is expected due to atrophy and potential pathology. We provide examples in Figure 1.

3.2 Question 2: Which is the optimal reference region?

Most of the debate in this literature has concerned choosing a reference region. We test 55 variations derived from four major regions: cerebellum, pons, midbrain, and supratentorial WM. The 55 variations come from various sub-regions, erosions, and combinations of the above. For supratentorial WM, we include two classes of implementations: those using each subject’s individual segmentation from T1-weighted MRI using SPM12 (Ashburner and Friston, 2005), and those using regions from the Johns Hopkins University “Eve” single-subject WM atlas (Oishi et al., 2009). We show examples in Figure 2.

3.3 Question 3: Which is the optimal GM segmentation within the target region?

In this question, we optimize which voxels to include in the cortical target. We are not comparing different sets of anatomic brain regions, but only how to apply segmentation to determine a subset of the voxels within these regions for each scan. We test a total of six segmentation methods within a fixed target meta-ROI containing those regions primarily affected by β-amyloid deposition: parietal, cingulate precuneus, prefrontal, orbito-frontal, temporal, and anterior cingulate (Figure 3). Within these ROIs, we test three variations with segmentations computed by the Unified Segmentation algorithm in SPM12 (Ashburner and Friston, 2005) with an in-house population-specific set of tissue priors, and three using Freesurfer version 5.3 (Fischl, 2012). For SPM methods, we test three variations: 1) voxels estimated to be either GM or WM together 2) voxels estimated to be GM with at least 50% confidence 3) voxels estimated to be GM with at least 95% confidence. Among Freesurfer methods, we test: 1) voxels segmented as either GM or WM together 2) voxels segmented as GM 3) voxels estimated to lie upon a surface halfway between the GM/WM and GM/CSF (pial) surfaces. Because Freesurfer does not output probabilistic segmentations, we use the output binary segmentations for the GM and GM+WM variations. Freesurfer segmentations are output in the space of its native Talairach template; we resample these back to each subject’s native space using the recommended technique (“How to Convert from FreeSurfer Space Back to Native Anatomical Space,” 2015). The “GM Half Surface” option uses the mri_surf2vol tools to estimate those voxels along a surface halfway between the GM/WM and the GM/CSF (pial) surface. These voxels are theoretically those cortical GM voxels furthest from WM. Segmentations produced by all methods were visually examined for all scans, and none had errors sufficient to warrant manual correction or exclusions. All segmentation and registration steps were performed separately for each scan at each timepoint, with the exception of Freesurfer segmentations, which used the longitudinal stream of Freesurfer.

3.4 Question 4: Which is the optimal analysis space?

Typical T1-weighted MRI scans have a resolution of approximately 1mm³, while PET scans typically have an effective resolution of approximately 6mm³ or larger (Joshi et al., 2009). For this situation, it is necessary to choose sets of voxels defined in some specific voxel space. Thus, resampling one or the other is necessary, but it is unclear which will produce superior results (Figure 4). We test both options (MRI space and PET space) in this work.

3.5 Question 5: SUV or SUVR?

Among discussion of reference regions, one might consider using SUV to normalize instead of a reference region (Table 1). We include this option in the set of reference methods tested, increasing the count to 56.

4.1 Subject Characteristics

We include scans of 129 subjects from the Mayo Clinic Study of Aging (MCSA) and Mayo Clinic Alzheimer’s Disease Research Center (ADRC) studies with three timepoints of serial scans with both PiB and T1-weighted MRI (for a total of 387 scans). Subjects with baseline SUVR > 2.5 were excluded, for reasons discussed later. MCSA is an epidemiological study of cognitive aging in Rochester, Olmsted County, Minnesota (Petersen et al., 2010; Roberts et al., 2008). The ADRC study recruits and follows subjects initially seen as patients in the Mayo Clinic Behavioral Neurology practice. All studies were approved by their respective institutional review boards and all subjects or their surrogates provided informed consent compliant with HIPAA regulations. We provide subject characteristics in Table 2. Freesurfer processing failed to produce an output for one subject, who was therefore excluded from Freesurfer-using pipelines.

4.2 Scan Acquisition Parameters

[11C] Pittsburgh Compound B (PiB) PET/CT studies were acquired using GE scanners (models Discovery 690XT and Discovery RX; GE Healthcare, Waukesha, WI). Subjects were injected with PiB (average 625 MBq, range 256–751 MBq) and a low dose CT scan was acquired. Beginning 40 minutes post-injection, subjects then underwent a 20-minute dynamic PET scan with four five-minute frames. Dynamic PET images were generated (256 matrix, 300 mm field of view, 1.17mm × 1.17mm × 3.27mm voxel size) using an iterative reconstruction algorithm. Standard corrections for attenuation, scatter, randoms and decay were applied as well as a 5 mm Gaussian post filter. The images from the four dynamic frames were averaged to create a single static image.

T1-weighted MRI scans (used for atlas normalization/masking, and for PVC where applicable) were acquired on 3T scanners (models Discovery MR750, Signa HDx, Signa HDxt, and Signa Excite) manufactured by General Electric (GE) using a 3D Sagittal Magnetization Prepared Rapid Acquisition Gradient-Recalled Echo (MP-RAGE) sequence. Repetition time (TR) was ≈ 7ms, echo time (TE) ≈ 3ms, and inversion time (TI) = 900ms. Voxel dimensions were ≈ 1.20mm × 1.015mm × 1.015mm.

4.3 Common Processing

T1-weighted scans were acquired with sagittal-plane gradient distortion correction performed on the scanner. Through-plane correction was performed as part of image processing (Jovicich et al., 2006). Tissue-class segmentation and inhomogeneity (B0 bias-field) correction were performed using the Unified Segmentation algorithm (Ashburner and Friston, 2005) in SPM12 (revision 6225). Several parameters in SPM12 were modified to produce more accurate segmentations for older-adult populations. Firstly, we used an in-house population-specific template and tissue priors that we call MCSA202. This template was created from MRI scans of 202 subjects in the same Mayo Clinic MCSA and ADRC studies from which this study’s subjects were selected. Each subject’s MRI was segmented using SPM12 and our custom template was created from these segmentations using the DARTEL groupwise registration algorithm in SPM12 (Ashburner, 2007). The tissue probability priors were manually edited to correct common segmentation errors, and 122 ROIs were drawn on the anatomic template. In addition to using the MCSA202 template, we also altered the SPM12 segmentation parameters to use two Gaussians to model WM intensities, instead of one, due to the higher prevalence of WM disease in such populations. We also reduced each of the stiffness penalty parameters of the nonlinear normalization of the tissue class priors of half of their defaults, allowing for increased inter-subject variability due to increased prevalence/severity of atrophy and other pathologies. More detailed information about the MCSA202 template and these parameter alterations is beyond the scope of this text, but are forthcoming in future publication. Images were registered to our MCSA202 template using the ANTs SyN registration algorithm (Avants et al., 2008) version 1.9.x with multiple channels: the post-B0-correction T1-weighted image, segmented tissue probabilities, and a mask of total intracranial volume. When each scan was segmented using Freesurfer 5.3, it was directly entered into Freesurfer version 5.3 (Fischl, 2012) using the recon-all pipeline (i.e. not preprocessed with the SPM12 bias-correction described above, because Freesurfer performs its own bias correction step and this would be redundant).

T1 and PiB scans were coregistered using SPM12 with 6 degrees of freedom (DOF). Resampling between MRI and PET resolutions was performed using ANTs software tools with 3^rd-order BSpline interpolation. Atlas ROIs were resampled to subject spaces also using ANTs with nearest-neighbor interpolation.

4.4 1,024 Pipelines

We compare a total of 1,024 different software pipelines for measuring β-amyloid load in PiB scans. In total, we examine 3 methods of PVC (none, 2-class PVC, and 3-class PVC), 56 intensity normalization methods (55 SUVR reference regions + SUV), 6 methods of cortical target segmentations, and 2 potential analysis spaces (MR-Space versus PET-Space). Multiplying all combinations, we implemented a total of 3×56×6×2 = 2016 pipelines. We then excluded pipelines in two classes of theoretically-implausible combinations: (1) 3-class PVC used with any of the 52 references containing WM or with either of the 2 targets containing WM ((1×52×6×2=624) + (1×56×2×2=224) − (1×52×2×2=208 in both) = 640), and (2) target segmentations containing supratentorial WM used with references also containing supratentorial WM (3×44×2×2=528). Class 1 was excluded because images corrected with 3-class PVC have signal in WM removed, making it illogical to then attempt to measure signal in these WM regions. Class 2 was excluded in order to prevent use of the same voxels for both target and reference regions, which would at least partially normalize out any signal of interest. After subtracting both classes, 2016 − 640 − 528 + (1×44×2×2=176 in both exclusion groups) = 1,024 combinations remained, which we analyze in this work.

4.5 Evaluation Criteria and Statistical Methods

In this section we describe the motivations and implementations for each of the four individual criteria, and the weighted, combined metric, used to compare measurement pipelines in this work. Each criterion was designed to address different characteristics desirable in serial PiB measurements. Ideally, one would wish to have an independent, gold-standard measurement of β-amyloid load with which to compare, but the only such measures come from pathological examinations, which are of course impossible for serial measurements. As such, we have chosen criteria based on obtainable data: the plausibility of the serial trajectories produced by each pipeline for each subject. All criteria analyzed are specifically longitudinal, because the goal of this study is to determine the best pipeline for serial measurements, which may or may not also be ideal for cross-sectional studies. Because none of these criterion is flawless nor alone captures all traits desirable in a serial PiB measurement pipeline, and because comparing four separate metrics across different implementations creates a complex array of data, we also present a single metric that is a weighted combination of the four. We describe each below.

6.1 Question 1: Is Partial Volume Correction Helpful?

Our results suggest that under most combinations of other factors, two-class (Meltzer) PVC improved results according to all of our criteria when compared to no-PVC pipelines. Three-class PVC was not applicable to a large majority of pipelines tested because these include WM, which 3-class PVC attempts to remove, in the reference regions. Three-class PVC was often superior among the pipelines with GM-only references, although pipelines using these GM-only references were otherwise not among the top candidates.

6.2 Question 2: Which is the optimal reference region?

This question has received the most attention in prior literature, and also most strongly impacted our analyses. Overall, our results strongly suggest the superiority of reference regions containing both supratentorial WM and whole cerebellum, optionally also including the pons or entire brainstem. Supratentorial WM may be implemented equivalently either by segmenting each individual T1 and eroding the WM voxels with a radius of 3 or 5 voxels, or by using an atlas where the relevant WM ROIs exclude WM near the cortex. Using individual segmentations without such erosion was generally worse. These combined-reference-ROI approaches significantly outperform those with any of their components individually (In Figure 10, the top-performing pipeline, which uses this combined-reference approach, differs significantly from all pipelines using non-WM or supra-/infratentorial regions alone). We visualize these winning variations in Figure 11.

6.3 Question 3: Which is the optimal GM segmentation within the target region?

Among the top-performing pipelines, differences between target segmentation methods were not significant (Figure 10). Across all pipelines, this had a relatively small impact on our quality criteria compared to other methodological factors.

6.4 Question 4: Which is the optimal analysis space?

In almost all cases, pairings of pipelines differing on only this variable trended toward slightly better performance in MRI-space. We omitted showing this data for space reasons.

6.5 Question 5: SUV or SUVR?

Pipelines using SUV, rather than a reference region (SUVR), were by far the worst performing on the straightness and separability criteria, and they were only mid-level performers on the other two. Using the combined criterion, even the worst-performing SUVR pipelines generally outperformed the best-performing SUV pipelines. Therefore, we conclude that use of SUVR with any reasonable reference region is superior to SUV.

6.6 Overall Recommended Pipelines

Our study suggests that an optimal pipeline to measure change in β-amyloid from PiB scans should use two-class PVC, perform calculations in the voxel space of the MRI, and use SUVR with a reference region containing both supratentorial WM and the whole cerebellum, optionally also including pons/brainstem (Figure 11).

6.7 Discussion

Our analysis is unique for its finding that references containing supratentorial WM and whole cerebellum together are significantly superior to either alone for analysis of serial PiB scans. Some previous comparisons using florbetapir have also examined these “composite” reference regions (Landau et al., 2015), but such results could theoretically not be the same for PiB. Some other works favoring supratentorial WM references have hypothesized that the cerebellum is noisier due to its relatively peripheral location in the field of view, where scanner sensitivity is lower, or due to its relatively smaller size (Chen et al., 2015; Landau et al., 2015). Our data supports the hypothesis that cerebellar ROIs are noisier than supratentorial WM references (In Figure 5, the longitudinal reliability criteria, references containing only cerebellar ROIs are among the worst performers). However, the statistical equivalence of including brainstem, and the superiority of some smaller WM variants and some larger WM variants, suggests that larger references are not always superior, perhaps due to diminishing returns, or because larger WM regions tend to include more contamination from adjacent GM.

Although we had no strong reason to believe a priori that our findings would agree with those of previous studies using Fluorine ligands, our findings were consistent with previous studies using Florbetapir in suggesting that references containing supratentorial WM are superior to purely-infratentorial references (Chen et al., 2015; Landau et al., 2015). This agreement with prior studies using differing ligands, populations, and measurement criteria adds to a growing body of evidence that supratentorial WM references are superior to purely-infratentorial references for serial Amyloid PET measurements regardless of tracer used. This agreement across studies using differing criteria also suggests that this finding is strong enough to be seen consistently despite the lack of gold-standard validation measures and despite the limitations of each criterion that has been used in their place.

6.8 Strengths and Limitations of Current Study

The primary strength of this work is its large scope, attempting to test most reasonable combinations of methods across several variables according to four different longitudinal criteria. We use three timepoints for all subjects, giving increased power and enabling more analysis criteria, such as trajectory reliability, versus previous studies using two. This work also addresses a relative lack of such comparisons using PiB, as opposed to ¹⁸F-based amyloid PET ligands. It is important to stress that our work applies only to longitudinal measurements of PiB. We do not assume that our findings necessarily apply to cross-sectional measurements; follow-up studies will be needed to examine this question. We also leave for future work a more exhaustive comparison of PVC methods with more variations, e.g. ROI-based methods (Rousset et al., 1998) and comparisons with methods using data-driven selections of individual voxels (Borghammer et al., 2009; Carbonell et al., 2015; Razifar et al., 2009). We also did not examine PET-only methods (Bilgel et al., 2015; Fripp et al., 2008; Raniga et al., 2007; Zhou et al., 2014) because MRI is generally available in most studies that include amyloid PET.

One limitation of our study is the potential for circularity in some criteria: although we examine different varieties of SUVR measurements, it was necessary to use cross-sectional SUVR measurements for study subject selection, and for their division into subgroups for the group separability criterion. Theoretically, such circularity could bias results toward favoring methods that most resemble the selection criteria. Although we attempted to minimize this potential by using a method that was not among those tested, the selection method was still based on SUVR with two-class PVC and a cerebellar gray reference, and this could have biased results toward similar pipelines. However, methods using cerebellar GM references performed relatively poorly in our analysis, suggesting that this potential bias was not strong enough to incorrectly favor these methods. For the group separability analysis, we also selected thresholded ranges that excluded many subjects between the extremes to minimize the ability for subtle changes in SUVR measurements to impact group selection.

Because external, gold-standard measurements (i.e. autopsy) of change over time in β-amyloid do not exist, our study is primarily limited by the assumptions of each criterion used instead (see Section 4.5). In addition to carefully designing each criterion to minimize the impact of these limitations, we also created the combined criterion that allows each to compensate for the other’s shortcomings, and we weighted each metric in proportion with our confidence in the universality of its assumptions in our dataset. While these weights were chosen a priori, our major findings were largely consistent across all four measurement criteria, and thus across different weightings of the combined criteria (see Supplementary Material). This consistency of our findings across criteria that measured fundamentally different properties of each pipeline suggests that these criteria were each reasonable, despite their varying assumptions and limitations, and adds confidence to our findings overall.

Our study adds to the growing evidence in favor of analysis methods that include WM voxels within reference regions. One potential caveat of WM-containing references for longitudinal analyses is the potential for a confound by changes in WM myelination, which have been shown to affect uptake of β-amyloid PET ligands (Stankoff et al., 2011; Veronese et al., 2015). It is possible that such WM changes over long periods of time could affect the value of a WM-containing reference ROI. However, this effect would likely be reduced by the inclusion of infratentorial voxels in composite references, which are those that performed best in this analysis.