The minimal preprocessing pipelines for the Human Connectome Project

Abstract

The Human Connectome Project (HCP) faces the challenging task of bringing multiple magnetic resonance imaging (MRI) modalities together in a common automated preprocessing framework across a large cohort of subjects. The MRI data acquired by the HCP differ in many ways from data acquired on conventional 3 Tesla scanners and often require newly developed preprocessing methods. We describe the minimal preprocessing pipelines for structural, functional, and diffusion MRI that were developed by the HCP to accomplish many low level tasks, including spatial artifact/distortion removal, surface generation, cross-modal registration, and alignment to standard space. These pipelines are specially designed to capitalize on the high quality data offered by the HCP. The final standard space makes use of a recently introduced CIFTI file format and the associated grayordinate spatial coordinate system. This allows for combined cortical surface and subcortical volume analyses while reducing the storage and processing requirements for high spatial and temporal resolution data. Here, we provide the minimum image acquisition requirements for the HCP minimal preprocessing pipelines and additional advice for investigators interested in replicating the HCP's acquisition protocols or using these pipelines. Finally, we discuss some potential future improvements to the pipelines.

Keywords: CIFTI; Grayordinates; Human Connectome Project; Image analysis pipeline; Multi-modal data integration; Surface-based analysis.

Figure 1

Components of the CIFTI grayordinates spatial dimension and standard space. Left and right cerebral cortices contribute about 30k surface vertices each (less than the 32k vertex standard meshes, because the non-cortical medial wall is masked with standard ROIs). Additionally, 19 subcortical grey matter structures combine to contribute about 30k volume voxels. In total, there are 91,282 grayordinates corresponding to all of the grey matter sampled at a 2mm average vertex spacing on the surface and as 2mm voxels subcortically. The minimal preprocessing pipelines guarantee that each subject will have 91,282 maximally aligned grayordinates, facilitating cross-subject comparisons of multiple data modalities.

Figure 2

A CIFTI dense timeseries. The rows are grayordinates and represent the spatial dimension of the timeseries, whereas the columns are timepoints. Reading along a row produces a timecourse (green) for a given spatial point, and a different spatial point produces a different timecourse (red). The columns represent the spatial pattern across the entire brain at a given timepoint (yellow) and can be also visualized. A different spatial pattern at a second timepoint is also shown (cyan).

Figure 3

An example of a CIFTI dense connectome. Both dimensions of this 2D matrix are grayordinates. Rows of the matrix represent seed points for functional connectivity, and the columns represent the spatial map of the connectivity. If one clicks on a point on the surface (green), Connectome Workbench finds the appropriate row in the dense connectome matrix, and displays the connectivity values on the appropriate cortical and subcortical structures. A different point (cyan) shows a different connectivity pattern.

Figure 4

A CIFTI dense scalar file. CIFTI dense scalar files can contain any kind of scalar data, but in this case, two task fMRI contrast maps are shown (z-statistics). Each contrast corresponds to a different column in the schematic matrix. The first task contrast, Places – AVG Categories (yellow), is from a working memory task that used a number of categories of visual stimuli, including places. The second task contrast, Story – Math (cyan), is from a language task that contrasted listening to a story to doing math problems. For more details, see Barch et al THIS ISSUE.

Figure 5

A comparison of the geometrical effects of various fMRI acquisition resolutions on folded cortex averaged across the Q1 20 unrelated HCP subjects all on the same color scale. A, 4mm isotropic, B, 3mm isotropic, C, 2.5mm isotropic, D, 2mm isotropic. When mapping fMRI data of a given resolution to the midthickness surface, larger voxels can artificially increase the correlation between geodesically distant surface vertices, as larger voxels can span both sides of a sulcus or a thin gyral blade of white matter. The measure being displayed is calculated for each vertex by finding the set of vertices that could potentially share a voxel's data with the selected vertex during volume to surface mapping with trilinear interpolation. For each of the vertices in this set, we compute the expected correlation to the selected vertex that would result from mapping a white noise timeseries by trilinear interpolation, aligned such that a voxel center is located midway between the two vertices. We then multiply by the geodesic distance between the vertices in order to upweight the measure for the vertices that are far apart across the surface. Finally the maximum of the measure is taken over the set of vertices found to be in range. This model is a significant simplification of the problem, however, as it does not account for the variable pointspread function of BOLD at 3T or the effects of larger draining veins. The model also uses a simple trilinear interpolation volume to surface mapping method instead of the more complicated ribbon-based method used in the pipelines below. This model does, however, point to the potential benefits of higher acquisition resolution in reducing the induced cross-sulcal and cross gyral correlations between vertices that are far apart along the surface (and thus potentially lying in different cortical areas), but close together in the volume.

Figure 6

Illustration of the spatial specificity of the high resolution HCP fMRI acquisition and minimal preprocessing pipelines in the left (top) and right (bottom) hemispheres of one HCP subject. Two resting state ICA components are shown, one involving the POS2 and retrosplenial cortex on the left and the other involving the peripheral visual cortex on the right. These two components exist on opposite sides of the parieto-occipital sulcus and have very distinct functional connectivity patterns. With the 2mm isotropic acquisition used here, the component spatial maps fall nicely between the white (green) and pial (cyan) surfaces that define the grey matter ribbon on either side of the sulcus.

Figure 7

An overview of the HCP Minimal Preprocessing Pipelines. The HCP Structural Pipelines are run first, and then either the HCP Functional Pipelines or the Diffusion Preprocessing Pipeline can be run. Afterwards, resting state, task, and diffusion analysis can proceed.

Figure 8

An overview of the volume and surface spaces produced by the HCP pipelines. Each space has specific uses, which are listed.

Figure 9

The steps of the PreFreeSurfer pipeline. T1w and T2w images are processed through the same steps in parallel up until cross-modal registration. A fieldmap is used to remove readout distortion in the T1w and T2w images, which slightly improves their alignment and resulting myelin maps. The PreFreeSurfer pipeline also produces the Native and MNI Volume Spaces, which will be used in later pipelines and analyses.

Figure 10

Effect of gradient nonlinearity distortion on the T1w image. Panel A is the T1w image before gradient distortion correction, whereas Panel B is the T1w image after gradient distortion correction. The red lines in both panels are from a rough tissue segmentation generated with FSL's FAST on the gradient distortion corrected image shown in Panel B. The ellipse highlights the region of maximum distortion, which is as much as the width of the grey matter ribbon in some places.

Figure 11

Precise alignment of T2w (A) and T1w (B) images after correction for readout distortion. The coronal slice above is in the same anatomical region as van der Kouwe et al (2008) Figure 4, which illustrates the effects of readout distortion. White matter (green) and pial (blue) surface contours from the FreeSurfer pipeline are in the same locations relative to the anatomy in both T2w and T1w images. The ellipses indicate the locations where readout distortion would be present in uncorrected images.

Figure 12

The steps of the FreeSurfer Pipeline. FreeSurfer's recon-all forms the basis of most of the pipeline, but it is interrupted at certain steps to improve the robustness of the brain extraction, to fine tune the T2w to T1w registration, and to more accurately place the white and pial surfaces with high-resolution 0.7mm T1w and T2w data. These inputs are fed back into recon-all at each stage.

Figure 13

The benefits (stars) of a high resolution (0.7mm isotropic) T1w image for reconstructing the white matter surface in thin, heavily myelinated cortical regions such as the early visual cortex in the calcarine sulcus (A) and the early somatosensory cortex in the central sulcus (B). The red surface contour is the output of standard recon-all on the downsampled 1mm isotropic T1w data and the green surface contour is after adjustment with the original 0.7mm isotropic T1w data. In general, the white matter surface is placed too superficially in these regions (Type A error, Glasser and Van Essen 2011).

Figure 14

The benefits of using a high resolution (0.7mm) T2w image to assist in pial surface placement. The blood vessels and dura are close to isointense with grey matter in the T1w image (A), however they are clearly of different intensity in the T2w image (B), allowing them to be easily excluded from the grey matter ribbon. The red surface is from the T1w image whereas the cyan surface is after cleanup with the T2w image.

Figure 15

The Conte69 (Glasser and Van Essen 2011) myelin maps before (A) and after (B) correcting a systematic Type B artifact where the pial surface was being placed too deep. In this artifact, lightly myelinated superficial cortical layers were being excluded in many subjects by the default recon-all algorithm (in FreeSurfer versions 4.5, 5.0, and 5.1) leading to artifactually higher myelin map values within the ellipse. Use of grey-matter specific normalization and the T2w surface eliminate the problem.

Figure 16

An overview of the PostFreeSurfer pipeline. In addition to generating a final brain mask and myelin maps, it converts data to NIFTI and GIFTI formats and creates Connectome Workbench spec files in the surface spaces shown in Figure 8 above.

Figure 17

The unsmoothed myelin maps from an individual HCP subject.

Figure 18

The effect of residual bias field correction on a particularly inhomogeneous HCP subject. A) Original Myelin Map. B) Reference group average myelin map (Conte69, from Glasser and Van Essen 2011). C) Estimated residual bias field. D) Bias corrected myelin map. The artifactual anterior-posterior very low spatial frequency gradient of T1w/T2w intensity is removed and the corrected myelin map now resembles the reference. All maps are displayed on the subject’s inflated surface.

Figure 19

The steps of the fMRIVolume Pipeline. All images undergo gradient distortion correction. Motion correction is to the single-band reference (SBRef) image. The SBRef image is used for EPI distortion correction and is registered to the T1w image. The spin echo EPI images have reversed phase encoding directions and are used to calculate deviations from the B0 field caused by susceptibility. All transformations, including the MNI space transformation are concatenated and applied to the original fMRI timeseries in a single spline interpolation step. Finally, this MNI space timeseries is masked and its intensity is normalized to a 4D global mean.

Figure 20

The steps of the fMRISurface pipeline, which maps cortical grey matter voxels onto standard surface vertices and subcortical grey matter voxels onto standard subcortical parcels. These combine to form a standard number of grayordinates that have been aligned according to the nonlinear surface and volume registrations respectively. This allows precise cross subject correspondence between data in the grayordinates standard space.

Figure 21

The coefficient of variation of the surface timeseries without (A) and with (B) exclusion of locally noisy voxels. While certain parts of cortex are very noisy due to signal loss from susceptibility effects, many other parts of cortex, particularly some gyral crowns, have locally large amounts of noise that is effectively removed with this technique.

Figure 22

The steps of the diffusion preprocessing pipeline. After b0 intensity normalization, the b0 images of both phase encoding directions are used to calculate the susceptibility-induced B₀ field deviations. Then the full timeseries from both phase encoding directions is used in the “eddy” tool for modeling of eddy current distortions and subject motion. Gradient distortion is corrected and the b0 image is registered to the T1w image using BBR. The diffusion data output from eddy are then resampled into 1.25mm native structural space and masked. Though not shown in Figure 22, the diffusion directions and the gradient deviation estimates are also appropriately rotated and registered into structural space.