A Comparison of Partial Volume Correction Techniques for Measuring Change in Serial Amyloid PET SUVR

No PVC

Although partial volume correction of PET images was first proposed decades ago [11, 22] and the methodology has continued to flourish, no strong consensus has emerged for which method is best, and publications using data without partial volume correction are still frequent in the literature. We include no-PVC as an option in this comparison to provide a reference measurement of how much variance is present in serial amyloid PET standardized uptake value ratio (SUVR) measurements without PVC, in the context of serial measurements.

Two-compartment (Meltzer-Style) PVC

Two-compartment (Meltzer-style) PVC is a voxel-based approach where each voxel is modeled as a linear combination of only two classes: tissue and non-tissue. This is one of oldest and most common methods for PVC, in part because it is easy to understand and implement. PET is coregistered to the corresponding MRI and the MRI is segmented to mark each voxel as either tissue or non-tissue. A mask of voxels segmented as tissue is then blurred by the assumed scanner PSF to produce an estimate of the fraction of tissue-originated signal contained in each PET image voxel. Individual PET scan voxels are then each divided by this estimated tissue fraction to correct for partial volume. The signal in non-brain voxels is, implicitly, assumed to be zero. For example, if a PET voxel is estimated to contain 80% tissue and 20% non-tissue according to the PSF-blurred MRI segmentation, its raw intensity is divided by 0.8 to estimate what the signal may have been if the voxel contained 100% tissue [16].

Three-compartment (Müller-Gärtner-style) PVC

Three-compartment PVC is also a voxel-based method that extends the two-compartment model such that uptake in gray matter and uptake in white matter (WM) are modelled separately, rather than together as a single “tissue” class. PET is coregistered to the corresponding MRI and the MRI is segmented to mark each voxel as gray matter, WM, or non-tissue. An atlas is used to locate and measure the mean PET signal inside the centrum semiovale, a large central region of WM. It is assumed that all WM has homogeneous uptake equal to this value, which is then assumed to be the partial volume corrected value for WM (we also examine a less-common, alternative method of estimating WM signal in the Supplementary Material). A binary mask of voxels segmented as WM is blurred by the PSF to produce an estimate of the fraction of WM-originated signal in each PET image voxel. This is then multiplied by the value from the centrum semiovale to estimate the contribution of WM to measured signal in each voxel. The resulting WM-contribution image is subtracted from the original PET image to remove WM signal in each voxel. Finally, the steps from two-compartment PVC (above) are performed on the resultant image to also correct for the effects of non-tissue voxels [12].

Geometric transfer matrix (Rousset-Style) PVC

GTM PVC, also known as region spread function (RSF) PVC, is a region-based method proposed by Rousset et. al. in 1998. Brain PET images are parcellated into many individual regions using MRI segmentation and atlas propagation, and the method attempts to model PSF interactions between entire regions (rather than voxels individually, as in two- and three- compartment PVC). PET signal values within each region are modelled as homogeneous. Consider n regions of interest. Each region (i) is individually segmented and convolved with a 3D Gaussian kernel with the width of the estimated PET PSF. This creates a probabilistic map, called a regional spread function (RSF_i), of how much signal from region i is assumed to be present in each image voxel. The mean RSF value is calculated within each original region of interest (j) to compute an estimate of signal from region i appearing in region j. The GTM model assumes that the observed PET signal in each region is the sum of true (unknown) PET signals in each region combined in a weighted sum where the mean RSF values serve as the weights. This leads to a system of linear equations which may be solved for the unknown, ideal (as if there were no PSF) regional PET values. For an atlas with n ROIs, each of these n equations with n unknowns is inserted into a n×n matrix of equations. Standard matrix methods are then used to solve for the n unknown, PVC-corrected values [13].

Unlike the voxel-based methods, GTM PVC attempts to correct for the effects of PSF causing wash-in and wash-out of signal between neighboring regions, rather than only between broad classes of tissue and non-tissue. A modification of GTM called Symmetric GTM has also been proposed, where the initial uncorrected per-ROI values are measured from the PSF-blurred ROIs (RSFs), rather than the un-blurred binary masks, to improve bias and noise characteristics [23]. This symmetric variant was recently included as the recommended PET PVC method in the latest 6.0 release of the popular FreeSurfer package [24], facilitating its proliferation.

Mayo SUVR pipelines

In this work we refer to a set of in-house software for calculating SUVR values from pairs of MRI and PET images as the Mayo SUVR pipelines. These are based primarily on components from the Statistical Parametric Mapping version 12 (SPM12) [25] and Advanced Normalization Tools (ANTs) [26] packages, and templates/atlases from our own Mayo Clinic Adult Lifespan Template (MCALT) package [27] (https://www.nitrc.org/projects/mcalt/). We have used these components to implement two pipeline variants: a cross-sectional variant where values for each timepoint are computed using only images from that timepoint, and a longitudinal variant where values for each timepoint are computed using images from all timepoints simultaneously. We describe both in this section.

Preprocessing

All input images were converted to nifti format from DICOM. MRI were additionally corrected for through-plane gradient distortion [28] (correction in the sagittal plane was performed on the scanner). These preprocessed images were directly used as inputs for both the Mayo and FreeSurfer pipelines described below.

Cross-sectional variant.

The cross-sectional Mayo PiB SUVR pipeline has been previously described [15, 29]. More recently we have updated it to use the MRI segmentation/normalization steps previously described in [30], and we use this version here. In brief, PET scans were normalized to the corresponding MRI using a 6DOF (rigid) transform computed with SPM12 and resampled using B-spline interpolation. We chose B-spline, rather than linear interpolation, to minimize unnecessary loss in image clarity due to resampling. Although a full comparison of PET image resampling methods is beyond the scope of this paper, we also replicated all of our experiments using linear interpolation and found no significant quantitative differences in measurement reproducibility for any variants (data not presented). The corresponding MRI was corrected for intensity inhomogeneity and segmented with Unified Segmentation [31] in SPM12 using MCALT tissue priors and segmentation parameters. Normalization parameters were calculated between the intensity-corrected MRI and the MCALT T1-weighted template using ANTS symmetric normalization [32], and MCALT atlases were transferred to the MRI space using ANTs’s GenericLabel interpolation. Atlas regions were masked to only include voxels of their corresponding tissue type(s) according to the MRI segmentation, and mean values within each region were computed. The “GlobalPiB” target region is a composite that includes the prefrontal, orbitofrontal, parietal, temporal, anterior cingulate, and posterior cingulate/precuneus regions [29]. A mask of this composite is available in the public MCALT release. Since the best way to normalize PET uptake is still an active area of research, we compared methods across a number of reference regions. The cerebellar crus, cerebellar gray matter, whole cerebellum, and pons reference regions were each defined using MCALT atlases. The supratentorial WM reference region was defined by using the MCALT lobar atlas to remove cerebellar/brainstem regions from a mask of voxels segmented as WM according to the MRI. Eroded supratentorial WM reference regions have been proven effective for serial amyloid PET analysis [33–36], but we included this region without erosion because FreeSurfer/PETSurfer contains no straightforward mechanism to erode WM regions or otherwise separate subcortical from deep WM regions. We do not advocate using non-eroded WM reference regions for amyloid PET due to contamination from gray matter signal, but we included this region in our comparisons because this was the only supratentorial WM region that we could implement similarly with both pipelines. A flowchart of the overall Mayo cross-sectional pipeline is available in Supplementary Figure 2.

Longitudinal variant.

The Mayo longitudinal PET processing pipeline uses the same key components as the cross-sectional pipeline, but critically differs such that all steps are performed by simultaneously using data from all available timepoints in an attempt to stabilize the resulting registrations, segmentations, and parcellations [37]. This approach creates a single-subject template (SST) for each subject’s set of serial PET images (PET-SST), and one for each subject’s set of MRI scans (MRI-SST).

Each MRI-SST was formed by using the buildTemplateParallel iterative group-wise nonlinear registration approach from ANTs [38] to build a mean-space template of all the MRI timepoints. This single image was then segmented and parcellated into regions using the same methods (SPM12, ANTs, MCALT) as in our cross-sectional approach described above, and used for quantifying each PET scan in place of that timepoint’s original MRI. We chose to segment the mean-timepoint image directly and use it for all timepoints, rather than segmenting each timepoint after resampling it to this mean space. We made this choice because conceptually it favors stability (only a single consistent segmentation/parcellation) over accuracy (tissue boundaries in the mean-space image less accurately match any of the individual timepoints). This is more distinct from the cross-sectional approaches and from the Longitudinal FreeSurfer approach, which we also compare. Thus, this choice compares a wider diversity of approaches. During a preliminary phase of this study, we compared a variant with individual-timepoint segmentations/parcellations in the mean space, and the results were not significantly different from the presented variant. For longitudinal studies with a longer time span, where intra-subject differences in atrophy would be larger, we would expect the differences between these approaches to be larger and the re-segmentation approach would be more appropriate. Approaches using similar principles have been shown to improve the reliability and biological plausibility of change measurements using MRI [39].

PET scans from all time points were coregistered to a common space to create a PET-SST using an in-house groupwise registration based on 6-DOF spm_coreg. First, each timepoint was independently 6-DOF-registered to the MCALT_T1 template and resampled in this space. A voxel-wise mean was calculated, forming an initial PET-SST. Each PET scan was then registered to this initial PET-SST, then a new mean was calculated, and the process was repeated for a total of 3 times to form the final PET-SST. A single rigid registration between the MRI-SST and PET-SST was calculated, and each PET scan was resampled (B-spline) to the space of the T1-SST by a composite of the PET⟶PET-SST transform and the PET-SST⟶T1-SST transform. Each of these transformed PET scans was then used with the segmentations/parcellations of the MRI-SST for SUVR quantification and PVC.

To summarize this process, all PET timepoints for a subject were quantified using a single, mean-space MRI of that subject that was derived from all its MRI timepoints. This approach of using a single mean-space MRI for all PET timepoints is designed to minimize variation due to instabilities across timepoints in the MRI processing steps (intensity inhomogeneity correction, segmentation, and atlas registration/parcellations). A mean-space single-subject PET template was also formed, and a single rigid transform was used between the PET and MRI modalities by computing a registration between both templates. This approach of a single PET-MRI registration is designed to minimize the effects of PET-MRI registration instability, which we have previously documented [40]. All steps were performed such that all timepoints were treated equally, to avoid causing bias [41] as in approaches that, e.g., resample all timepoints to match the baseline. We provide a flowchart of all steps in this method in Supplementary Figure 2.

2-compartment and 3-compartment PVC implementations.

For the Mayo pipelines we used in-house implementations of 2- and 3-compartment PVC according to the methods in the original studies [11, 12]. We applied these voxel-based corrections directly prior to calculating the mean values in each atlas region. We used tissue segmentations of the subject MRI from SPM12 as described above, and assumed a PSF of 8 mm full width at half maximum (FWHM). This value was determined by averaging transaxial and axial PSF measurements from an internal study of scans using a F-18 point source at the center of a water phantom, reconstructed identically to the brain imaging protocols.

GTM PVC implementation.

For the Mayo pipelines, we adapted the core GTM PVC logic from a publically-available Matlab implementation of Rousset-style GTM PVC [42]. That implementation was specifically designed for tau PET and to use FreeSurfer brain parcellations as inputs, but it is completely distinct from the PETSurfer implementation later included in FreeSurfer 6.0. Our in-house adaptations only retained the core Rousset-style GTM PVC logic (i.e., RSF estimation and matrix solving) from this implementation. We rewrote the functions around this core to instead use input regions from a combination of our SPM12 MRI segmentations and MCALT atlas regions (instead of FreeSurfer segmentations) to create a dense 132-region atlas where every voxel is included, analogous to but distinct from the gtmseg step in PETSurfer (below). To form these regions, structures in the MCALT_122 and MCALT_Lobar atlases were first parcellated as such, and then classes for cerebrospinal fluid, air, skull+dura, and WM were added using the SPM12 segmentations. We did not include regions/corrections for extra-cortical hotspots from the Baker implementation, because these are specific to tau PET. Unlike FreeSurfer, our atlas splits WM into subcortical and deep WM, but to match FreeSurfer’s regions we combined these afterward to form a total supratentorial WM reference region. As with the voxel-based methods, we assumed a PSF of 8 mm FWHM.

PETSurfer (FreeSurfer 6.0) SUVR pipelines

FreeSurfer is a widely-used software package for surface-based analysis of brain MRI [43]. FreeSurfer analyses can be performed using either the cross-sectional stream or the longitudinal stream, which segments all serial timepoints at once for increased longitudinal stability [44]. In version 6.0, FreeSurfer introduced a new module called PETSurfer, which is used to register a PET scan to a corresponding, previously-segmented MRI and output tables of PET values in each atlas region [45]. Regional values can be output either with or without GTM PVC applied. PETSurfer has no longitudinal pipeline for analyzing serial PET scans using multiple timepoints simultaneously; however, when analyzing PET scans it can use MRI segmentations that were produced using either the cross-sectional or the longitudinal FreeSurfer pipeline, and the latter serves as the recommended approach for serial data [46].

We first preprocessed MRI scans by converting to nifti format and correcting for gradient unwarping as in the Mayo pipelines described above. After this preprocessing, we used PETSurfer with FreeSurfer 6.0 to produce uncorrected and symmetric-GTM-corrected ROI values using standard settings and PETSurfer procedures (https://surfer.nmr.mgh.harvard.edu/fswiki/PetSurfer). We refer to values produced using the cross-sectional MRI segmentations as FreeSurfer_CrossSec, and those using the longitudinal segmentations as FreeSurfer_Longitudinal. As in the Mayo pipelines, we assumed a PSF of 8 mm FWHM. For each reference tested reference region, we used the –rescale option to mri_gtmpvc to provide the appropriate region number(s). To form a composite target region from FreeSurfer’s atlases that is analogous to our GlobalPiB composite region, we used a size-weighted mean of PETSurfer’s output regional SUVR values (with each choice of reference region) of all cortical regions except for those marked as occipital, pre/postcentral, and insula.

Data inspection and preprocessing

We first performed a graphical analysis of outliers and noticed the presence of a small fraction of extreme values due to failures in segmentation or registration. Based on this preliminary analysis, we chose to treat SUVR values that were outside the range of [0,7] as pipeline failures. Although amyloid PET SUVR values of up to 6 would be extremely large by standard SUVR measures, such values are not uncommon when using GTM PVC, and we wanted to use a standard threshold for all methods. Using these thresholds, we identified any subjects for which any method produced an SUVR value that was either outside this range, or not a number (i.e., zero-mean or zero-voxel reference regions). This occurred for only two subjects (1 clinically unimpaired (CU), 1 mild cognitive impairment (MCI)), which we excluded from all further analyses. Next, we identified n = 11 subjects (8 CU, 1 AD, 1 other) that did not have exactly three timepoints of computed SUVR values with every method due to failures in imaging processing pipelines. All of these were failures from FreeSurfer-based methods (no failures were observed with the Mayo methods). Further analyses continued with the remaining n = 278 subjects.

SUVR values from the FreeSurfer pipeline were then plotted versus those from the Mayo pipelines separately by reference region and PVC method in order to validate each of these independent implementations by examining to what extent the estimates agreed (i.e., clustered about the identity line, y = x) and their correlation (i.e., clustered about a regression line).

Estimating the coefficient of variation for each method

We used linear mixed effects regression to model natural logarithm transformed SUVR values over time to estimate the magnitude of within-subject variation for a given measurement method. Even in relatively homogeneous populations PET SUVR values tend to be right skewed due to measurement errors being proportional rather than additive. After a log transformation, SUVR values tend to be approximately normally distributed and have variance that is comparable across the log(SURV) distribution and therefore approximately additive. By a property of the log-normal distribution, the standard deviation (SD) of log-transformed values equals the coefficient of variation (CV) of the original distribution, where CV is the SD divided by the mean. We therefore interpret the residual standard deviation from a mixed of log-transformed SUVR values as providing an estimate of the CV of the method and an indication of the magnitude of within-subject measurement error.

For our primary analysis, our mixed model specification was as follows: We treated time from baseline as a fixed effect and included a random intercept and a correlated random slope for each individual. For each individual, time zero was defined as the midpoint between the first and third scan. We used the arm package’s sim() function [47] to perform 10,000 parametric bootstrap/posterior simulations of the mixed-effects model regression coefficients and variance parameters for each method. We used these simulations when plotting means and 95% confidence intervals.

We assume that change in SUVR over relatively short intervals can be modelled as linear and that deviations from linearity can primarily be attributed to measurement error that an ideal method would minimize. However, fitting a least squares line for each subject results in a high degree of overfitting and would severely underestimate the residual error. By using a mixed model approach, individual intercepts and slopes are shrunk toward the overall population average to reduce overfitting and account for regression to the mean. The resulting model-based residual SD arguably better reflects the true CV of the measurements.

Our justifications for these assumptions have been previously published [36]. In brief, we wish to penalize methods that produce biologically implausible triangle-shaped trajectories (i.e., where values strongly increase and then decrease, or vice versa, over short intervals of time). It is known that amyloid PET trajectories in AD subjects are a roughly sigmoidal shape [15] with an accumulation period of ≈19 years [48]. Having three measurements over a span of ≈3 years, it is reasonable to assume that trajectories should be locally linear. Therefore, we consider significant acceleration or deceleration in measured SUVRs to be much more likely due to measurement error than a true change in amyloid.

Group-wise differences in slopes

To assess how each method might perform for a hypothetical study comparing rates of amyloid accumulation across two groups, we performed a second statistical analysis wherein we separated our subjects into two groups: clinically unimpaired and clinically impaired (MCI + dementia together). We fit the same mixed model as above, but added fixed effect interaction terms to estimate an intercept and slope for each group. We then plotted the difference in slopes (interpreted as a percentage difference due to the log-transformed response variable) between the impaired and unimpaired groups for each method. We recognize that some unimpaired individuals can be expected to be accumulating amyloid while some impaired individuals can be expected to have more or less stable levels of amyloid. However, given that the lowest age in our sample is 51, the proportion of unimpaired subjects who are non-accumulators should outweigh that which might be pre-symptomatic accumulators [49]. Secondly, amyloid continues to accumulate into the symptomatic phase (MCI and early AD), so it is not the case that rates are on average flat in symptomatic individuals [15]. Therefore, it is reasonable to assume that amyloid accumulation should be faster, on average, in impaired subjects than in unimpaired subjects and a correctly-functioning method should detect this difference (we also present an alternative analysis where groups were defined based on amyloid status at visit 3, in the Supplementary Material). This test is additionally important because it would be easy to create a hypothetical measurement method that is highly precise by always returning the same value. Such a method would always win in comparisons based only on precision measures. Therefore, it is important to measure group-wise discrimination for each method to ensure that precision is not eliminating the underlying signal.