Resting-state fMRI in the Human Connectome Project

Abstract

Resting-state functional magnetic resonance imaging (rfMRI) allows one to study functional connectivity in the brain by acquiring fMRI data while subjects lie inactive in the MRI scanner, and taking advantage of the fact that functionally related brain regions spontaneously co-activate. rfMRI is one of the two primary data modalities being acquired for the Human Connectome Project (the other being diffusion MRI). A key objective is to generate a detailed in vivo mapping of functional connectivity in a large cohort of healthy adults (over 1000 subjects), and to make these datasets freely available for use by the neuroimaging community. In each subject we acquire a total of 1h of whole-brain rfMRI data at 3 T, with a spatial resolution of 2×2×2 mm and a temporal resolution of 0.7s, capitalizing on recent developments in slice-accelerated echo-planar imaging. We will also scan a subset of the cohort at higher field strength and resolution. In this paper we outline the work behind, and rationale for, decisions taken regarding the rfMRI data acquisition protocol and pre-processing pipelines, and present some initial results showing data quality and example functional connectivity analyses.

Fig 1

A) Temporal-SNR estimations for various acquisition protocols from the FMRIB multiband motion piloting, using just the “normal head motion” runs. For each resting-fMRI run, and for each protocol, the raw temporal-SNR image was formed and its median value found. This was also calculated from the data after ICA-based artefact cleanup. The boxplots show distributions over the 6 subjects. It is clear that both decreasing the voxel size and increasing the acceleration result in much lower SNR. B) However, if the increased number of timepoints is taken into account, in terms of its effect on simple timeseries statistics, it is clear that the acceleration is of great statistical value, and approximately counters the loss in SNR caused by the increase in spatial resolution.

Fig 2

Simplified graphical overview of HCP rfMRI data organization and analysis flow, including example generation of dense and parcellated connectomes at the group level. Each subject's 15-minute rfMRI dataset is spatially and temporally pre-processed, resulting in two versions of the pre-processed timeseries data – volumetric (in 3D MNI152 space) and grayordinates (surface vertices plus subcortical and cerebellar grey matter voxels). Either of these standard-space versions of the data can be combined across runs and subjects, but the grayordinate version is more compact (it contains only grey-matter data) and should provide better alignment across subjects than the volumetric version. A simple group-level analysis might concatenate all subjects’ datasets in time; from this, the “dense connectome” of average correlations could be estimated. This can be fed into a parcellation, and from this (via the parcels’ associated mean timeseries estimated from the group timeseries data) one can estimate the “parcellated connectome”.

Fig 3

Example SBRef and single-timepoint fMRI images, before and after corrections for distortion. Top: SBRef images acquired with L-R phase encoding (left) and R-L phase encoding (right). The central columns are the raw images in native 2mm space; the outer columns are the images after distortion corrections. Middle: equivalent example single-timepoint images from the 4D fMRI timeseries; these are in the same space as the SBRef images, having the same dropout and distortion, but with lower SNR and tissue contrast. The red dilated-brain-edges are the same in all cases, being derived from the average L-R and R-L corrected SBRef images, to allow for easier comparisons across different images. The orange arrows indicate example areas of signal dropout, which are different in L-R and R-L, and are not removed by the distortion corrections. The blue arrows indicate example areas of distortion, which are well corrected by the distortion corrections. Bottom: distortion-corrected SBRef images after alignment to the structural images, with the white-grey boundary (estimated by FreeSurfer from the structural images) overlaid in green. The left-right and right-left phase encoding direction SBRef images are shown separately; there is no obvious residual distortion that is different between these two images. The same corrections are applied to the fMRI timeseries data, which have the same distortions as the SBRef images, but we show just the latter here as the better SNR and tissue contrast make it easier to see the high quality of the alignment to the structural data. The cortical surface overlay views and all surface renderings in following figures were created using the Connectome Workbench display tool (humanconnectome.org/connectome/connectome-workbench.html).

Fig 4

Structural and functional images, averaged across 20 HCP subjects (i.e., 20 structural images and 80 rfMRI runs); left is right. Images are in MNI152 space, after all spatial pre-processing, including distortion correction (gradient corrections + Topup), head motion correction (FLIRT), affine registration of functional to structural (FLIRT+BBR) and nonlinear registration of structural to MNI152 space (FLIRT+FNIRT). The mean structural image is shown both in native 0.7mm resolution and after resampling to 2mm; all other images are shown in 2mm resolution. The SBRef images are reference-EPI images with no multiband acceleration and no T₁ saturation. Overlaid, in colour, are edges derived from the mean structural image (different colours indicate different edge gradient strengths); these show excellent alignment and lack of distortion in the EPI data. SBRef cross-subject averages are also shown separately for the L-R vs R-L phase-encoding directions; the asymmetry in dropout can be seen, but there is very little residual distortion. The multiband-accelerated EPI mean-timeseries images are also shown; these have the same (well-corrected) distortion and dropouts as the SBRef images, but much poorer tissue contrast. The asymmetry in dropout is quantified by subtracting the mean R-L image from the mean L-R, dividing by the maximum of the two, and multiplying by 100; this is shown in the colour overlay, having a maximum difference of approximately 60%. Finally, mean tSNR images are shown, as well as mean of 1/tSNR. These are shown for data after temporal highpass filtering in both cases, and comparing without vs. with artefact cleanup (removal of bad ICA components and motion confounds). The histograms show the distributions of tSNR values in the two cases, with a tSNR of 40 marked on the x axis. The maximum display intensity for the mean tSNR images is set to 50 in both cases, and the maximum display intensity for the mean 1/tSNR images is set to 0.1 in both cases.

Fig 5

The multiband-accelerated EPI mean-timeseries images, averaged across 20 subjects (80 rfMRI runs). This is similar to the mean EPI images shown in the previous figure, except that the timeseries are averaged across subjects in grayordinate space and only the cortical surface data are shown. The intensity display is arbitrary units; the highest mean intensity is 5 times greater than the lowest.

Fig 6

Examples of “bad” (above) and “good” (below) components from ICA applied to a single 15-minute resting-state fMRI session. The component spatial maps are shown (in red/blue, overlaid on the raw fMRI data), as well as the component timecourse and its power spectrum. Insets show expanded views of 3 slices from the spatial maps. The classifications of these particular components are quite clear from looking at the spatial maps and the timecourses, whereas in this case the spectrum of the bad component is not obviously artefactual. This is an example of why different features are important for accurate classification in different components. The display tool (“Melview”, an in-house program developed specifically for this purpose) is a convenient way to visualise and hand-label components, for feeding into the FIX classifier training.

Fig 7

The hand classification of the 25 HCP subjects was carried out using Connectome Workbench with combined surface and volume visualisation. The volume views looked very similar to those of the “Melview” program. Two examples of Connectome Workbench surface views of two single-run ICA components’ spatial maps are shown. One component is clearly artefactual (above), and the other non-artefactual (below). The combined volume and surface approach was found to be useful for several reasons. 1) Cortical signal of interest is generated from the cortical grey matter ribbon and thus always maps onto the surface as strong and distinct patches of activation. 2) Artifacts often don't map onto the surface or do so irregularly (e.g., more on gyri than sulci). Thus comparison of volume and surface maps is often enough to distinguish good signal from bad. 3) Projection to the surface can make certain artefactual patterns easier to spot (such as artefacts in the axial slice plane).

Fig 8

Example RSN (parts of the default mode network) from a single 15-minute run from a single subject. ICA was run on the volumetric data, and the resulting RSN's map is shown at the top, overlaid onto the single-band reference scan. On the bottom is shown the corresponding spatial map in grayordinate space. Both views are thresholded at abs(Z)>3. The black dot corresponds to a local maximum, and is used as the seed location for the dense connectomes (correlation maps) shown in the following figure.

Fig 9

Functional connectivity (full correlation converted to Z-statistics, using the FSLNets package) between a seed point (single grayordinate seed) in the default mode network and the rest of the cortical grayordinates. The correlation map is shown for a single-run, a single subject (4 runs concatenated) and 20 subjects (80 runs concatenated). Positive correlations are thresholded at Z>5 and negative correlations are thresholded at Z<-2.5 (negative correlations tend to be weaker than positive, so we use a more liberal threshold here (but still consistent across the 3 datasets) in order to show interesting anti-correlated structure.

Fig 10

Functional connectivity maps in two nearby seed locations. The top row shows individual (left) and group (right) connectivity maps from seeds in the retrosplenial cortex (the larger marker is the seed in each case). The middle row shows individual (left) and group (right) maps from seeds in the immediately adjacent posterior cingulate cortex. The bottom row shows individual (left) and group (right) functional connectivity gradients that highlight the location of this change in functional connectivity. The functional connectivity colour palette is scaled so that the 98^th percentile is yellow (or cyan if negative) and the 2^nd percentile is black. The gradients are scaled between 96% (red) and 4% (black).

Fig 11

Spatial maps of the ICA components found by MELODIC group-ICA, run on all 20 subjects (80 runs), with an ICA dimensionality of 30. 8 components were excluded (being judged as either artefactual or highly inconsistent across subjects), leaving the 22 components shown. All colour overlays show Z-statistic versions of the spatial maps, thresholded at Z>5. (A) shows the original spatial maps from the group-ICA, carried out on the grayordinate versions of the fMRI datasets. (B) shows representative axial slices (3 per component) from the spatial maps in volumetric MNI152 space.

Fig 11

Spatial maps of the ICA components found by MELODIC group-ICA, run on all 20 subjects (80 runs), with an ICA dimensionality of 30. 8 components were excluded (being judged as either artefactual or highly inconsistent across subjects), leaving the 22 components shown. All colour overlays show Z-statistic versions of the spatial maps, thresholded at Z>5. (A) shows the original spatial maps from the group-ICA, carried out on the grayordinate versions of the fMRI datasets. (B) shows representative axial slices (3 per component) from the spatial maps in volumetric MNI152 space.

Fig 12

Mixed-effects cross-subject group-level spatial maps estimated in volumetric (top) and in grayordinates space (middle and bottom), in both cases thresholded at Z>7. The grayordinates map is shown in the bottom row overlaid on the “midthickness” cortical surface that runs halfway between the outer and inner grey matter boundaries. The same map is shown in the middle row on the “very inflated” cortical surface.

Fig 13

The increase in peak mixed-effects Z-statistics when carrying out cross-subject group-level analysis in grayordinate space instead of in volumetric MNI standard space. For each RSN from the 100-dimensional group-ICA, the maximum Z-statistic (across space) was computed for both analyses, and the difference between the two maxima computed. The histogram shows the distribution of this difference across RSNs. There is no RSN having a larger mixed-effects Z-statistic peak value in volumetric space. The mean difference (reflecting the extent to which grayordinate-based cross-subject modelling is superior) is 2.8.

Fig 14

22×22 correlation matrices derived from the timeseries associated with the 22 group-ICA components. Below the diagonal is shown the full correlation; above the diagonal is shown the partial correlation matrix. Each row or column is the set of correlations between a single network “node” and all other nodes; the nodes have been reordered from the original ordering, according to a hierarchical clustering algorithm applied to the full correlation matrix (depicted at the top), that attempts to form clusters of nodes, seen as blocks along the diagonal of the full correlation matrix (more clearly seen with the larger number of components in the following figure). See the main text for discussion of the network edges marked as A-E. The figure is generated using the FSLNets package (fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSLNets).

Fig 15

Correlation matrices and node clustering derived from the 100-dimensional group-ICA, resulting in 78 non-artefactual components; see previous figure caption for further details.

Fig 16

Spectra and amplitudes of the RSN timeseries found from regression of group-ICA spatial maps into individual runs (in grayordinates). The boxplots show the distribution of the timeseries amplitudes across all runs and all components (red crosses mark outlier values). The spectra are averaged across all runs and all components.

Fig 17

Group Z-statistic maps for one component from the 30-dimensional group-ICA, without and with artefact cleanup. Top: Slices concentrating on the superior cortical parts of the RSN, thresholded at abs(Z)>7 and brightest yellow/blue corresponding to Z=±40. Bottom: The same RSN, but with slices intersecting the cerebellum; as the signal is weaker here than in cortex, the maps are thresholded at abs(Z)>5, with brightest yellow/blue corresponding to Z=±15.