An automated method for identifying artifact in independent component analysis of resting-state FMRI

Abstract

An enduring issue with data-driven analysis and filtering methods is the interpretation of results. To assist, we present an automatic method for identification of artifact in independent components (ICs) derived from functional MRI (fMRI). The method was designed with the following features: does not require temporal information about an fMRI paradigm; does not require the user to train the algorithm; requires only the fMRI images (additional acquisition of anatomical imaging not required); is able to identify a high proportion of artifact-related ICs without removing components that are likely to be of neuronal origin; can be applied to resting-state fMRI; is automated, requiring minimal or no human intervention. We applied the method to a MELODIC probabilistic ICA of resting-state functional connectivity data acquired in 50 healthy control subjects, and compared the results to a blinded expert manual classification. The method identified between 26 and 72% of the components as artifact (mean 55%). About 0.3% of components identified as artifact were discordant with the manual classification; retrospective examination of these ICs suggested the automated method had correctly identified these as artifact. We have developed an effective automated method which removes a substantial number of unwanted noisy components in ICA analyses of resting-state fMRI data. Source code of our implementation of the method is available.

Keywords: ICA; artifacts; automated classification; automatic; fMRI; functional magnetic resonance imaging; independent component analysis; independent component labeling.

Figure 1

ICA was applied to pre-processed fMRI data yielding spatial component maps with associated time courses and power spectra. SOCK automatically distinguishes between ICs dominated by artifact (Artifact category) and those containing possible neuronal signal (Unlikely Artifact category) by calculating IC features (smoothness measure, edge, CSF, and temporal frequency power).

Figure 2

An overview of identifying “spotty” artifact. The example data is taken from an ICA of a control that underwent a 10 min resting-state fMRI study (see Section 2.5). (A) Discrete Fourier Transform of a single slice of an IC. While a 3D Discrete Fourier Transform is implemented in SOCK, we show a 2D illustration for simplicity. Data in the center of figure contains low spatial frequency information about the image, while data near the periphery represents high spatial frequencies. The sphere with volume V1 is an arbitrary region which contains low spatial frequency information. (B) A plot of the ratio of low to high frequency information [equation (6)] vs. the radius of sphere V1 for all ICs. (C) A spatial map of the top and bottom curves corresponding to the smoothest and least smooth IC respectively. (D) Applying k-means clustering to split the ICs into a set of “smooth” (red curves), “subsmooth” (green curves), and “unsmooth” (blue curves) ICs. Each line in (B,D) represents the ratio for a particular IC calculated over all slices, not just the slice shown in (A).

Figure 3

An overview of identifying motion artifact. The example data is taken from an ICA of a control that underwent a 10 min resting-state fMRI study (see Section 2.5). (A) An edge mask (in green) automatically created by an SPM sub-routine, “New Segment” (single slice shown). This is overlaid onto a mean functional image. In red are an illustration of contiguous clusters which overlap the edge mask. (B) A plot of the edge activity for all ICs. (C) A spatial map of the top and bottom points corresponding to the highest and lowest edge activity ICs respectively. (D) Using k-means clustering, the ICs are automatically divided into a set of “High Edge Activity” components (points colored red in the scatter-plot) and “Low Edge Activity” components (points colored blue in the scatter-plot). Each point in the scatter-plot represents the edge activity for a particular IC summed over all slices.

Figure 4

An overview of identifying CSF artifact. A CSF mask automatically created by an SPM sub-routine, “New Segment” is shown in green (single slice shown). This is overlaid onto a mean functional image. In red are an illustration of contiguous clusters which overlap the CSF mask. If the volume of activity overlapping the CSF mask is 10% or greater, the IC is labeled High CSF Activity (Low CSF Activity otherwise). The CSF activity for a particular IC calculated over all slices, not just the slice shown above.

Figure 5

An overview of identifying ICs with Temporal Frequency Noise (TFN). The example data is taken from an ICA of a control that underwent a 10 min resting-state fMRI study (see Section 2.5). (A) A plot of the TFN activity (sum of power spectrum values from 0.08 Hz to Nyquist frequency) for all ICs. (B) Applying k-means clustering to split the ICs into a set of Low and High TFN ICs. (C) A spatial map of the top and bottom points corresponding to the highest and lowest TFN activity ICs respectively. (D) The associated power spectra for these ICs with the blue and red curves representing the low and high TFN ICs respectively. This is a zoomed in view showing only frequencies beyond 0.08 Hz which is the region of interest.

Figure 6

A selected set of spatial maps from an ICA and corresponding SOCK classification for Subject 4. The numbering of the ICs is based on the order of extraction in the ICA decomposition.

Figure 7

Spatial maps from an ICA and corresponding SOCK classification for all discordant ICs among the first 30 subjects (data sets 1 and 2). All these ICs were classified by an expert as unlikely artifact, whereas SOCK classified them as artifact. The numbering of the ICs is based on the order of extraction in the ICA decomposition. (A,B) Subjects 1 and 8: thresholded spatial maps of IC27 for subject 1 and ICs 7, 23, and 96 for subject 8 overlaid on an edge mask (in green, single slice shown) reveals a significant proportion of activation overlapping the edge mask. (C) Subject 5: a closer look at the unthresholded spatial map of IC90 (shown on the right of the thresholded map), reveals that the IC is not as smooth as it appears; compared to viewing the thresholded spatial map.

Figure 8

Spatial maps from an ICA and corresponding SOCK classification for all discordant ICs among the last 20 subjects (data set 3). All these ICs were classified by an expert as unlikely artifact, whereas SOCK classified them as artifact. The numbering of the ICs is based on the order of extraction in the ICA decomposition. (A) Subject 41, IC70: SOCK rejected this IC as artifact as it was unsmooth. However, an IC with overlapping spatial regions, but with a greater degree of smoothness was observed in the ICA decomposition (compare unthresholded maps of IC70 vs. IC24). IC24 was not classified as artifact by SOCK. (B) Subject 41, IC75: SOCK rejected this IC as artifact as it had substantial temporal frequency noise (see shaded region of power spectrum). However, ICs with overlapping spatial regions, but with less temporal frequency noise were observed in the ICA decomposition (compare power spectra of IC75 vs. ICs 26 and 29). These ICs were not classified as artifact by SOCK.

Figure A1

Receiver Operating Characteristic (ROC) curve for edge and CSF thresholds for 5 test subjects. A total of 525 combinations of edge and CSF thresholds (ranging from 0 to 70%) were tested. Shown in light gray are combinations where the classification in the 5 subjects did not change significantly (within 10%). These represent the ranges, 40–50% for the edge threshold and 20–30% for the CSF threshold. Lower values resulted in at least one neuronal component being misclassified as artifact (i.e., having a sensitivity of less than one which we regard as failure), whilst higher values resulted in fewer artifacts being identified (which we regard as a decrease in performance or smaller specificity). To minimize the chance of failure we chose the highest thresholds before a decrease in performance occurs (i.e., well away from the failure condition); 50 and 30% for the edge and CSF thresholds respectively.

Figure A2

A total of 525 combinations for edge and CSF thresholds ranging from 0 to 70% were tested. Each combination represents an edge and CSF threshold, with combination 0 representing 0% edge and CSF threshold and combination 525 representing 70% edge and CSF threshold. Shown on the graph is combination 367, which corresponds to 50% edge threshold and 30% CSF threshold.