Considerations and code for partial volume correcting [18F]-AV-1451 tau PET data

[18F]-AV-1451 is a leading tracer used with positron emission tomography (PET) to quantify tau pathology. However, [18F]-AV-1451 shows “off target” or non-specific binding, which we define as binding of the tracer in unexpected areas unlikely to harbor aggregated tau based on autopsy literature [1]. Along with caudate, putamen, pallidum and thalamus non-specific binding [2], [3], we have found binding in the superior portion of the cerebellar gray matter, leading us to use inferior cerebellar gray as the reference region. We also addressed binding in the posterior portion of the choroid plexus. PET signal unlikely to be associated with tau also occurs in skull, meninges and soft tissue (see e.g. [4]). We refer to [18F]-AV-1451 binding in the skull and meninges as extra-cortical hotspots (ECH) and find them near lateral and medial orbitofrontal, lateral occipital, inferior and middle temporal, superior and inferior parietal, and inferior cerebellar gray matter. Lastly, the choroid plexus also shows non-specific binding that bleeds into hippocampus. We are providing the code (http://www.runmycode.org/companion/view/2798) used to create different regions of interest (ROIs) that we then used to perform Partial Volume Correction (PVC) using the Rousset geometric transfer matrix method (GTM, [5]). This method was used in the companion article, “Comparison of multiple tau-PET measures as biomarkers in aging and Alzheimer's Disease” ([6], DOI 10.1016/j.neuroimage.2017.05.058).


Specifications
Value of the data -[ 18 F]-AV1451, although useful, displays some nuisance off-target binding that, due to partial volume effects, affects regions of the brain (temporal cortex, hippocampus, reference region) that are of great interest in studying tau as it relates to aging and disease. -Fundamentally we are providing the code (uploaded with the article) that will create an ideal set of ROIs for studying tau deposition, along with a set of ROIs that must be accounted for if proper partial volume correction is to be performed. -This code addresses the reference region in two ways: removing superior portion of the cerebellar gray which typically shows [ 18 F]-AV-1451 binding thereby contaminating the reference region, and also scans the remaining reference region (inferior cerebellar gray) for the rare occasion of isolated increased [ 18 F]-AV-1451 binding. -The PVC code specifically addresses uptake in the choroid plexus, splitting choroid plexus into 2 clusters of higher and lower uptake. It also addresses off-target binding of tracer in the skull creating clusters of extra cortical hotspots (ECH) used as ROIs in the partial volume correction algorithm.

Data
We are sharing the code (http://www.runmycode.org/companion/view/2798) used to create regions of interest and deploy the GTM approach for partial volume correction [5] of [ 18 F]-AV-1451 tau PET data. This code creates an inferior cerebellar gray reference region, splits the choroid plexus region into high and low binding, creates subject-specific extra-cortical hotspot regions, creates high and low binding CSF and skull þmeninges regions.
The purpose of the code is to thoughtfully create 1. a set of ROIs that explore interesting scientific questions, and 2. a set of ROIs where nuisance off-target binding occurs resulting in partial volume effects in the former ROIs. We will focus more on 2. In order to test the effect of specific steps in the ROI creation for the PVC, we tested 10 different ROI configurations: 1. FreeSurfer ROIs grouped as shown in Table 2 plus the whole choroid plexus. Any voxel not defined by FreeSurfer (cerebrospinal fluid (CSF), skull, meninges) was assumed to be 0. 2. Same as 1, but choroid plexus was split into high and low choroid plexus ROIs. We found choroid plexus uptake was higher in the ventral portion that was next to the hippocampus, and lower in the dorsal portion and wanted to explore the effect of segmenting this ROI into high and low on the PVC. 3. Same as 2, but also included extra-cortical hotspots (ECH). ECHs were determined on a subject-bysubject basis. An ECH was defined as a cluster of 4500 voxels with SUVR 41.6, each ECH was added as its own ROI for PVC. 4. Same as 2, but included an ROI for SPM12 c3 (CSF) and another ROI for c4 þ c5 (skull þ meninges), no ECHs. 5. Same as 3, ECHs with threshold of 1.6, but included an ROI for SPM12 c3 and another ROI for c4þc5 (as described in step 4). 6. Same as 4, but the c3 mask was divided into a c3 high (SUVRZ 1) and c3 low (SUVR o1), and c4þc5 mask was divided into c4þc5 high (SUVR Z1) and c4 þc5 low (SUVR o1) masks. Did not include ECHs. 7. Same as 5, but the c3 mask was divided into a c3 high (SUVR Z1) and c3 low (SUVR o1), and the c4þc5 mask was divided into c4 þc5 high (SUVR Z1) and c4þ c5 low (SUVR o1) masks. Included ECHs using 1.6 as threshold. 8. Same as 7 with ECH threshold was 1.6, and added a search to remove voxels 4 1.6 from inferior cerebellar gray. 9. Same as 7, but ECH threshold was 1.3 and added a search to remove voxels 4 1.6 from inferior cerebellar gray. 10. Same as 7, but ECH threshold was 1.9 and added a search to remove voxels 4 1.6 from inferior cerebellar gray.
In order to quantify the success of PVC with the 10 different ROI configurations, within a subject and ROI configuration, a PVC image was created (example: Fig. 1C), in which the PV-corrected SUVR value of each ROI was assigned to all voxels within that ROI. This image was then smoothed by the  resolution of the scanner to create a calculated pre-PVC image (example: Fig. 1D). If the ROI configuration and smoothing kernel of the scanner explained the original SUVR perfectly, the calculated pre-PVC image would look exactly the same as the original SUVR, although estimation of the smoothing kernel, subject motion, inhomogeneity of ROIs, coregistration between PET and MRI, and imperfect segmentation of the MRI are a few possible sources of error. The difference between the original SUVR and the calculated pre-PVC image was calculated (example: Fig. 1E), these are the residuals. The mean within ROIs of the difference between original SUVR and calculated pre-PVC is 0, or close to 0, and does not offer any information as to how well the PVC ROI configuration fits the original data. However, the standard deviation within an ROI is useful in quantifying the performance of an PVC ROI configuration; it is lower when the PVC configuration fits the data better and higher for worse fits. We normalized the within ROI standard deviation by the mean of the original SUVR (PVCnstd) in order to quantify how well an ROI configuration is performing. Fig. 2 shows the mean and standard deviation of PVCnstd within each subject group averaged over ROIs 1-77 from Table 2. Including c3 and c4þc5 masks without ECHs in the PVC (steps 4 and 6) causes an increase in the PVCnstd. Therefore, not including ECHs in the ROI configuration results in a worse fit of the data. If steps 4 and 6 are ignored, there was a decrease of PVCnstd from ROI group 2 to 3 to 5 to 7, reflecting an improvement of the fit to the data when ECHs, c3 and c4þc5 masks and splitting those masks into high and low uptake are included. The addition of scanning inferior cerebellar gray for voxels 4 1.6 showed no improvement in the model because this rarely occurs. Setting the threshold for defining ECH at 1.3, 1.6, and 1.9 (steps 9, 10 and 8, respectively) also had no effect overall on the success of the PVC from a global quantitative perspective.
ROI configuration 8 was chosen as the final version [6], although we have written the code so that it can flexible. Information about the ROIs added to the ROI configuration for step 8 is found in the description of the code below.  Table 1 þ choroid plexus. ROI configuration 2, choroid plexus was divided into low and high ROIs. ROI configuration 3 was the same as 2 plus ECH (SUVR threshold¼1.6). ROI configuration 4 was the same as 2, no ECHs, and the c3 mask and c4 þ c5 mask were added. ROI configuration 5 included ECHs (threshold ¼ 1.6) and c3 ROI and a c4 þ c5 ROI. ROI configuration 6 was the same as ROI configuration 4 except the c3 mask and c4 þc5 mask were both divided into high and low masks (threshold¼1). ROI configuration 7 included ECHs (threshold ¼ 1.6), c3 low and high (threshold¼ 1) and c4 þc5 low and high (threshold ¼ 1) ROIs. ROI configuration 8 was the same as 7 except the inferior cerebellar gray was scanned for high voxels. ROI configuration 9 was the same as 8 except a threshold of 1.3 was used for ECH. ROI configuration 10 was the same as 8 and 9 except a threshold of 1.9 was used for ECH. The ROIs used for calculating the normalized standard deviation corresponded to those in Table 1.

The code
The main function is TAUPVC_RUNME_Create_ROIs_For_Rousset. The inputs into this Matlab function are full paths and filenames for: 1. aparcþ aseg.nii from FreeSurfer, 2. c1.nii, c2.nii, c3.nii, c4. nii, and c5.nii from SPM12, 3. Reverse normalized Cerebellar SUIT template (resliced to the  dimensions of the MRI), 4. Mean or sum of realigned [ 18 F]-AV-1451 frames coregistered and resliced to the MRI, and 5. Approximate PET scanner resolution. The outputs are 1. [ 18 F]-AV-1451 normalized to inferior cerebellar gray, and 2. an edited version of aparcþ aseg file in which voxels from each ROI used for PVC are assigned a different integer value, 3. an SUVR image normalized by the inferior cerebellar gray cortex with no partial volume correction, and 4. PV-corrected values are saved in variable roigroups in a final FINAL0_roigroups.mat.
The first step in the code was to assign new index values and group ROIs (as described in Table 2) using the aparcþ aseg.nii file from the FreeSurfer segmentaion. Ventricles were unassigned and therefore assumed to have a 0 PVC value.
The second step was to create the inferior cerebellar gray ROI from the reverse-normalized cerebellar SUIT template (rnCereSUIT). The inferior cerebellar gray was used as the reference region due to frequent binding seen in the dorsal cerebellum as well as bleeding in from adjacent cortical regions (Fig. 3). The rnCereSUIT was divided into binary masks for the inferior portion (indices 6, 8-28, 33 and 34) and the superior portion (1-5 and 7) of the cerebellar gray. These masks were smoothed by 8mm 3 , ensuring that all voxels defined as cerebellar gray in the original aparcþaseg.nii had a nonzero value in either the smoothed inferior or smoothed superior mask. This was necessary because the original aparcþ aseg cerebellar gray region did not perfectly overlap with rnCereSUIT. If a voxel was defined as cerebellar gray in the original aparcþaseg.nii and the smoothed inferior mask 4 smoothed superior mask, it was added to the final inferior cerebellar gray ROI. If a voxel was defined as aparcþaseg cerebellar gray and the smoothed superior mask 4 smoothed inferior mask, it was added to the final superior cerebellar gray ROI. A new SUVR was calculated by normalizing the [ 18 F]-AV-1451 SUVR by the mean inferior cerebellar gray; this newly-normalized SUVR was used for the rest of the code.  5. young healthy control subject with ECH near the occipital lobe. A is the original SUVR image. B has the calculated pre-PVC for the ROI group with an ECH threshold of 1.3. C is the same as B but ECH threshold was 1.6, and D had an ECH threshold of 1.9. The third step was to take the choroid plexus, as segmented by FreeSurfer, and divide it into low and high choroid plexus ROIs (Fig. 4). The choroid plexus was not a homogenous ROI but showed a bimodal distribution of tracer binding and therefore had to be divided into choroid plexus high ROI and low ROI. The ventral portion of the choroid plexus that is next to the hippocampus shows offtarget binding (Fig. 4A, red arrows). In order to realistically quantify the activity in the hippocampus, this had to be addressed with PVC. Bimodal distribution cutoff, kurtosis and skewness is shown in Supplementary figure 1 for the whole choroid distribution across all subjects. We defined high and low choroid plexus ROIs based on if voxels had more or less binding than the reference region. Within the choroid plexus, the low choroid plexus ROI contained voxels where the SUVRo ¼1 and had at least 100 contiguous voxels; the high choroid plexus ROI contained voxels where SUVR 41 and had at least 100 contiguous voxels. This did not account of all choroid plexus voxels. So the high and low choroid plexus masks were smoothed (8 mm 3 ). Any unassigned voxels within aparcþaseg choroid plexus were assigned to high choroid plexus if smoothed high choroid plexus 4 smoothed low choroid plexus, and low choroid plexus if smoothed low choroid plexus 4 smoothed high choroid plexus. In the edited_aparc þaseg.nii, high choroid plexus voxels ¼79 and low choroid plexus voxels ¼ 80.
The fourth step was to define ECHs. ECHs are defined in each subject's native space, and are defined on an individual basis. This was done in order to make no assumptions regarding the location and spread of ECHs. As mentioned above, we explored 3 possible thresholds (1.3, 1.6, 1.9) for ECH before finalizing our PVC code. They resulted in similar values for PVCnstd for cortical ROIs, so we relied on qualitative inspection of the calculated pre-PVC images in comparison to the original SUVR. Fig. 5 is an example of such an image. A threshold of 1.3 included too many voxels and 1.9 at times resulted in too small ECHs, 1.6 seemed to be the best compromise.
In the implementation, we wanted to look for ECHs in c4 and c5 (skull þ meninges), but also some of c3 that was not close to cortex. c3 is defined as CSF in SPM12 but this probability mask also contains parts of dura and soft tissue where ECHs have occurred, unrelated to cortical activity. To define the ECH-searchable ROI, we used the SPM12 segmentation c3, c4, and c5 probability maps, a brain mask, and the SUVR. The brain mask was a binary mask of the non-zero voxels in the aparcþ aseg, this smoothed to the resolution of the scanner. The c4 and c5 probability maps were added together, a binary mask of any voxels in c4þc5 4 0.5 was created and smoothed to the resolution of the scanner. We then searched in the c3 probability mask: any voxels with a probability 40.5 of being c3 (CSF) and that were closer to the c4þ c5 mask than to the brain mask (smoothed Fig. 6. Age versus volume of ECH; Spearman correlation shows age is correlated with the volume of ECH in healthy controls but is not significant in the AD/MCI subgroup within that subsample. c4þ c5 mask 4 smoothed brainmask) were added to the ECH-searchable ROI. Any voxels where c4þ c540.3 were also added to the ECH-searchable ROI. Within the ECH-searchable ROI, voxels with SUVR4 1.6 were considered probable ECHs. Contiguous probable ECH voxels were separated into individual clusters, and considered a full-fledged ECH if the cluster contained 4 500 voxels. Each ECH was added to the edited_aparcþaseg with a unique integer value starting at 85 and counting up for each subsequent ECH. We found that the number of overall ECH voxels was negatively correlated to age in healthy controls (p o 0.001, Fig. 6). However, the correlation with age to number of ECH voxels was not significant when looking only in healthy female controls (n¼51, age range ¼ 53-93.5 years) but was significant (p o0.001) in age-matched healthy male controls (n¼34, age range ¼ 56.5-93.8).
The fifth step was to assign remaining c1, c3 and c4þc5 voxels that were not defined in the original aparcþ aseg.nii (aparcþaseg ¼0) and were not classified as ECHs. Voxels were defined as CSF low (or c3 low) if 1. CSF (c3) or gray matter (c1) probability 40.3, 2. CSF (c3) probability 4 0, 3. CSF (c3) probability 4 bone þ meninges probability (c4 þc5) and 4. SUVR was 4 0.1 and o 1. CSF high met the same criteria except SUVR was 4 ¼1. For bone þ meninges low (c4 þc5 low), c4 þc5 had to have a probability 4 0.3, c4 þc54c3 and SUVR 4 0.1 and o 1. Bone þ meninges high met the same criteria as bone þ meninges low except SUVR was 4 ¼1. In the edited_aparcþaseg file, CSF low was assigned a value of 81, CSF high ¼82, bone þ meninges low ¼ 83, and bone þ meninges high ¼ 84. The SPM12 gray matter mask (c1) was included in the search because we found there was not 100% overlap between FreeSurfer defined gray matter voxels and SPM12 defined gray matter voxels. In this case, if a voxel was not defined as gray matter in FreeSurfer, but SPM12's probability of the voxel being gray matter was 30% or higher, we assumed SPM12 misclassified the voxel and it should be classified as CSF.
Next, voxels with gray matter mask (c1) 4 0.3, 0% probability of being c3, c4, or c5, SUVR 40.1, and unassigned in the edited_aparcþaseg (¼0) were addressed. For each of these voxels, a matlab function looked at the integer value of the surrounding voxels and assigned the stray voxel the integer value most common to its neighbors. We wanted to assign unassigned voxels because most of these voxels were gray matter missed by the FreeSurfer segmentation that showed some tracer uptake, and if left unassigned the value would assume to be 0 during PVC.
Lastly, the inferior portion of the cerebellar gray was searched for clusters of voxels with SUVR4 1.6 and more than 500 contiguous voxels in a cluster (Fig. 3). If any existed, each was individually assigned a unique integer value (counting up after the highest ECH integer value) in the edited_aparcþaseg. A new mean of the inferior cerebellar gray was calculated with these nuisance reference region hotspots removed, the SUVR was normalized one last time to this new mean of the inferior cerebellar gray. Across all subjects in [6] (n ¼216) inferior cerebellar ECH was only found in three subjects.