Fast Detection and Reduction of Local Transient Artifacts in Resting-State fMRI

2.1. Subjects

Fifteen right-handed males (mean age 20±2.14 years) without a known history of neurological or psychiatric disorders participated in this study. Written, informed assent and consent were obtained from all participants and/or their parent/guardian. The National Institutes of Health Institutional Review Board approved the study.

2.2. Image acquisition

Subjects were sequentially scanned using an 8-channel head coil array equipped GE 3 Tesla (3T) scanner. For six of the total fifteen scans, the scanner had a hardware system malfunction identified to be a defect in one channel of the head coil. Subjects N01-N09 comprised the normal data set and were scanned in the absence of the head coil defect. Subjects A01-A06 comprised the abnormal data set and were scanned using the same head coil, while the head coil defect was present. Subject order and scan dates are shown in Table 1. During scanning sessions, all subjects were instructed to fixate on a cross in the center of the screen. Resting-state (RS) fMRI time series were acquired for 490 s (140 volumes) by a T2*-weighted gradient echo pulse sequence with high spatial resolution (1.72×1.72×3.00 mm³, TR=3.5 s, TE=27 ms, flip angle = 90 °). Respiration volume and heart rhythm were monitored with a respiration belt and pulse oximeter. Belt diameter and pulse oximeter readings were sampled at 50 Hz during scanning. Two high-resolution (0.94×0.94×1.20 mm³) T1-weighted anatomical images were also acquired with an MPRAGE sequence.

The data sets from the 15 subjects were assessed by three QC tests assessing local white matter signal: local white matter regressor (WMe_LOCAL) [3], global correlation (GCOR) [15], and coefficients of variance of EPI data (CVAR_EPI) [16] (see Table 1). The GCOR measure is computed as the brain-wide average correlation over all possible combinations of voxel timeseries. The calculation of GCOR is simplified by calculating the average dot product of the average unit-variance time series. The CVAR measure is computed the coefficient of variation of input voxels, which is the standard deviation of timeseries divided by the absolute value of mean timeseries.

2.3. Defining the target artifact

Unlike head motion artifacts, which can be visualized in EPI data during scanning, hardware artifacts may be revealed when correlation maps show marked biases that violate neurophysiological norms. To target a hardware artifact, one might consider a persistent or transient, local, abnormal correlation between brain regions. In a previous study, we found that the temporal signature of a local artifact could be found across the time series of multiple tissue voxels within a region imaged by a defective array component, not exclusively in deep WM regions that could have less partial volume effect with gray matter voxels [3]. Consequently, there should be high correlations between the time series in a voxel in the artifact region and in the time series in the surrounding area, assuming no involvement of other, stronger artifacts such as large head motion [17].

All procedures in this study used AFNI packages (http://afni.nimh.nih.gov) [18]. The high-resolution anatomical image was aligned to the fifth EPI volume of the EPI time series [19], and then processed with texture-based tissue classification to get the WM masks of individual subjects [20]. All WM masks were then resampled to EPI resolution. To reduce partial volume effects across different tissue types, we eroded WM masks by one voxel along each of the three axes. The first 4 volumes of the RS time series were removed to ensure that all remaining volumes in the time series were at magnetization steady state. Despiking was done to suppress the spiky artifact from rapid head motion and to enhance the accuracy of motion correction [17]. Rigid-body registration was used to estimate subject movement during RS EPI scans and to correct for slice timing [9]. All MRI images presented in this work were spatially normalized to the N27 brain template for ease of interpretation by brain region [21].

To precisely identify the target artifact, we began with a basic regression model for RS fMRI, defined as follows:

Here, yⁱ is the EPI time series vector at voxel i. X_MO is the matrix of second order Legendre polynomials to fit quadratic temporal trends and six regressors containing rigid-body motion parameter estimates (three translations and three rotations) intended to model any residual effects of movement after motion correction. X_RI is the matrix of eight slice-wise regressors for retrospective image correction (RETROICOR) that models the effect of respiration and heart cycle [1], and a respiration volume per time (RVT) regressor that models slow blood oxygenation level fluctuations [2]. In Eq. (1) and subsequent equations, $r_1^i$ is the residual signal referring to the RS fMRI signal with unwanted nuisance components projected out.

The model for ANATICOR is defined in Eq. (2), where $X_{WM{e}_{LOCAL}}^i$ is the average local WM signal for each voxel i:

Using the spatial pattern of correlation maps with the seed point at right middle frontal gyrus (RMFG) and marginal explained variance (R²) maps for the regressor $X_{WMeLOCAL}^i$, the subjects were categorized into two groups; normal cases (N01-N09) and cases with hardware artifacts (A01-A06).

2.4. Outlier tests

We used QC measures in the AFNI package with three criteria: the outlier test, head motion estimation, and GCOR [15]. In the outlier test, at each time point of EPI data, the inter-TR image differences were calculated as the number of outlier time points. In each time series, time points that differed from the trend were deemed outliers. The threshold level for the outlier time point (e_thr) was defined by:

Where, N is the length of time series, M is the median absolute deviation of the time series, and D⁻¹ is the inverse of the reversed Gaussian cumulative distribution function. For all subjects, the number of outlier voxels was less than 1 percent of the number of voxels in the brain mask, and the Euclidean L²-norms (“enorm” in the AFNI package) of head movement were less than 0.25 (a.u.). Therefore all subjects met a practical threshold for minimizing artificial signal bias by head motion [17]. Additionally, GCOR of the whole brain was calculated to detect a concentrated correlation pattern that might be a signature of artifactual time courses [15], with the temporal noise-to-signal-ratio (tNSR) of EPI time series as a typical QC measure.

2.5. Artifact removal for fast denoising

For fast denoising, the nuisance model utilizing WMe_LOCAL (ANATICOR) has been limited by calculation time. For this data set, calculating the localized WM signal (WMe_LOCAL) for a single volume took over fifty TRs on average. To secure fast processing within 1 TR time, we replaced the localized regressor with isotropic Gaussian-distribution weights by distance from the voxel i:

Where, Wⁱ is a vector for the distance weight by isotropic Gaussian kernel of full-width-at-half-maximum (FWHM) of 30 mm at each voxel i. The processing time was remarkably reduced in this model, in part by replacing the voxel loop in the pipeline with a weighted white matter regressor (wWMe_LOCAL). Diagrams of the kernel weighting methods are presented in Fig. 1.

The local WM regressor considers the strength of covarying artifacts from each head coil array, which will differ by distance from the voxel i (see Figs. 1-C1 and 1-C2).

We replaced X_RI (9 components) with respiration volume per time X_RVT (1 component) to decrease processing time and to secure degrees of freedom for residual time series at the 30-TR time window.

The models BASE, ANATICOR, Weighted ANATICOR, and Fast ANATICOR (FANATICOR) contain 18, 19, 19, and 11 regressors, respectively.

2.6. Artifact indicators for image quality control

To demonstrate the transient nature of the local artifact, we used the R² and tNSR maps of the wWMe_LOCAL regressor to visualize the spatial distribution of the hardware artifact over time. For time windows of 30 TRs, we calculated the tNSR for two seed points, one in the LMFG and one in the artifactual region, the RMFG. To evaluate fast QC in those windows, we directly compared R² maps of the nuisance regressors used in the denoising regression models: NULL, MO, MO+RVT, and Fast ANATICOR.

2.7. Spatial patterns of captured artifact and seed-based correlation maps within the time window

To visually inspect the spatiotemporal characteristics of the transient artifact for the entire time-windowed series (105 sets of 30 TR time windows), we obtained the R² maps of the nuisance regressor used in the denoising regression models (NULL, MO, MO+RVT, and Fast ANATICOR) for normal and abnormal scans. To determine if each nuisance model could appropriately remove the targeted transient artifact over the full time series, we obtained the seed-based correlation maps for the RMFG seed point within the artifactual region for the artifactual group. The correlation and R² maps from the artifactual subjects were compared with those of the normal subjects.

2.8. Variances of sliding-time-windowed correlations

Typically, real-time fMRI applications consider regions-of-interest (ROI) instead of individual voxels in order to shorten computation time and improve ease of interpretation (see [22], 2011, for review). To determine ROIs for each subject, EPI data were first split into 105 sets by 30-TR time windows. The residual time series were then obtained by the preprocessing pipeline (see Fig. 2) with the four denoising models (NULL, MO, MO+RVT, and Fast ANATICOR) and averaged within 90 ROIs using the 90-fROIs atlas from Shirer and colleagues [23]. For each denoising model and each individual subject, 105 time-windowed Pearson correlation matrices (90×90 matrix size) were calculated, and then transformed to Fisher’s z scores. To observe variance changes in time-windowed correlations for the different denoising models, the CVAR of time-windowed correlation (CVAR_ρ) was calculated.

2.9. Non-stationarity of sliding-time-windowed correlations

We also calculated the non-stationarity of the correlation matrices along the time axis by the augmented Dickey-Fuller unit root test [24]. Non-stationarity varied between the different preprocessing methods, with Fast ANATICOR showing the most non-stationarity. For individual time-windowed correlation matrices from an individual subject, we obtained a binary matrix representing significant non-stationarity at a threshold p < 0.05 from the unit root test. A frequency (the number of subjects) map showing significant non-stationary correlations was then calculated for the 15 individual subjects with each denoising model.

3.1. Time costs of denoising and image quality control

The total elapsed time for nuisance regression by the Fast ANATICOR model was about 3.8 s for a single central processing unit (CPU), which is a marked improvement upon the previous method requiring 161.9 s per CPU on average. The calculation of the wWMe_LOCAL regressor in the Fast ANATICOR model was about 52 times faster than the typical WMe_LOCAL calculation in the ANATICOR model (see Table 2 for details). For QC with the wWMe_LOCAL regressor, mapping tNSR_wWMe could be accomplished at about 1 s per CPU, but mapping R² took much longer, 5.8 s for a single CPU, on average (see Table 2).

3.2. Transient artifacts detected by weighted local white matter regressor

We tested the ability of nuisance models BASE, ANATICOR, Weighted ANATICOR, and of Fast ANATICOR to capture patterns of hardware artifacts present in our dataset. Fig. 3 shows how these indicators analyzed data from subject N01, who had a normal scan, and from subjects A01 and A06 whose scans had known hardware issues. In Fig. 3-A, there is no visual artifact pattern in the EPI data for both normal and abnormal scans, although abnormal scans did show abnormal correlation patterns between RMFG and other brain areas (see Fig. 3-B). Interestingly, we found that the tNSR of wWMe_LOCAL and WMe_LOCAL could detect the hardware artifact (see Fig. 3-F).

In our dataset, the ANATICOR and Fast ANATICOR models differed significantly in processing time. The captured artifact patterns in R² maps and in correlation maps with a seed point within the artifactual area (RMFG) were highly similar in the two models (see Fig. 4).

The tNSR changes in wWMe_LOCAL regressors (tNSR_wWMe) at two seed locations within and outside of the coil artifact area, at the right and left middle frontal gyri (RMFG and LMFG), respectively, are presented in Fig. 5-A. When the transient artifact occurred in the RMFG, the tNSR_wWMe at the RMFG increased. The tNSR_wWMe at the LMFG, which was outside the artifactual area, did not change. During the time window when the artifact was absent [5, 34], the tNSR_wWMe map showed no artifactual spatial patterns and the correlation maps for the RMFG seed calculated by the four denoising models were nearly indistinguishable (see Fig. 5-B) (see Fig. 5-B). Notably, during the time window [51, 80] when the artifact was present, the artifactual pattern was detected in the tNSR_wWMe map, and the Fast ANATICOR model could clearly reduce its effect on the correlation map (see Fig. 5-C).

3.3. Explained variances of nuisance regressors

We created R² maps with a seed point in the RMFG for the normal and artifactual scans across the entire time series and compared them using each of the four denoising models. Fig. 6 shows this comparison for subject N01, who had no coil issues during scanning. As expected for the normal scan, no spatially concentrated signals were captured by the nuisance regressors of all the denoising models and the time-windowed correlation maps varied along time intervals with varying degrees and ranges in each model. Subjects N02-N09 showed similar results to N01.

For the abnormal scan of subject A01, at time intervals with high tNSR_wWMe values at the seed location, correlation patterns from the NULL, MO, and MO+RVT models showed the hardware artifact (see Figs. 7-A, B-1, and C-1), but R² maps did not. Interestingly, R² maps of wWMe_LOCAL captured the hardware artifact at the appropriate time intervals (see Fig. 7-D4). The Fast ANATICOR model not only consistently identified the artifact, but effectively removed it, showing correlation patterns similar to those of the normal subjects at all time intervals (see Fig. 7-D1).

3.4. Coefficients of variances for sliding-time-windowed correlations

For all subjects and nuisance models used in this study, CVAR_ρ matrices of time-windowed correlation matrices across the functional ROIs are shown in Fig. 8. The NULL and Fast ANATICOR models showed the fewest and most CVAR_ρ values, respectively.

3.5. Non-stationarity of sliding-time-windowed correlations

For all subjects and nuisance models in this study, the frequency matrices of subjects that showed significant non-stationary correlations are shown in Fig. 9. The NULL and Fast ANATICOR models showed the least and most non-stationarity, respectively. In general, the non-stationarity of time-windowed correlations was greater in correlations related to specific functional networks: auditory, basal ganglia, sensorimotor, visuospatial, anterior salience, and primary visual.

4.1. Possibility of real-time denoising for sub-second repetition time

Calculation time limits the application of correction methods in real-time fMRI. The bottleneck created by the voxel-loop in ANATICOR’s pipeline has contributed significantly to processing time. We demonstrated marked temporal improvements to the ANATICOR model by replacing the voxel-loop with a weighted white matter regressor and using tNSR as an QC measure, thereby creating Fast ANATICOR (see Table 2, Fig. S1). With calculation speeds that permit QC during image acquisition, real-time implementation of Fast ANATICOR will depend on the processing stream and data transfer speeds. Sub-second processing for denoising can be accomplished by multi-core computing systems containing more than 12 computing cores or CPUs with parallel computing frameworks such as OpenMP libraries (http://openmp.org). Regarding system integration for real-time fMRI applications, however, more computing power might be required to operate accurate QC by tNSR_wWMe, with real-time denoising by Fast ANATICOR. Current data transferring methods may need improvement in order to permit simultaneous transmission of data to both QC and denoising/analyzing systems that have been developed separately. Our suggested pipeline for real-time imaging applications is described in Fig. 2.

4.2. Variance and non-stationarity change in time-windowed correlations by denoising models

As more nuisance regressors were removed in the interest of calculation time, the CVAR_ρ and non-stationarity in time-windowed correlations increased within the ROIs (Figs. 8 and 9). This trend was more evident in the Fast ANATICOR model than the MO or MO+RVT models, even though the head-coil artifact partially overlapped with the visuospatial, right executive control, and anterior salience networks. It could be suggested that nuisance signals operate like a drift process in the time course of time-windowed correlations and that the localized artifact could be a strong inducer of stationarity between the time-windowed correlations.

4.3. Fast ANATICOR and data quality

In this study, we used relatively high-resolution EPI data. Notably, the WM erosion step typically requires EPI data of a resolution comparable to or higher than that of the data used here. At lower-resolutions, such as 3.5 mm isocubic voxel, little remains of WM voxels after the erosion operation. Fast ANATICOR uses a Gaussian-curved weighting kernel to compute each voxel-wise regressor, giving more weight to the white matter that is closest to a given voxel. If there is no WM voxel close to a given voxel (due to a poor segmentation or stringent erosion operation in a low-resolution data grid, for example), the Gaussian curve will likely find WM voxels farther away instead of creating an empty regressor.

4.4. Physiological nuisance regressors

One physiological regressor, RVT, was included in this Fast ANATICOR model due to processing time limitations. In an interesting contrast to previous studies using typical, longer series of EPI data [3], the MO+RVT model in this study suppressed more artifactual correlations than NULL and MO models during the artifactual time interval (see Fig. 5-C). It is possible that physiological nuisance signals have a stronger effect in the input data with smaller time points in a time window analysis than the full-length EPI data.