The optimized combination of aCompCor and ICA-AROMA to reduce motion and physiologic noise in task fMRI data

The data for this study was obtained from a fMRI study investigating the effects of transcranial direct current stimulation (tDCS) on episodic future thinking (EFT) (results yet to be published). The included data sets were from 24 females and 6 males (aged 25 ± 3 years) without any neural or psychiatric disorder. All participants were Dutch speaking and right-handed and signed an informed consent form. The current retrospective and the original study were approved by the ethical committee of UZ Brussels.

All scans were performed using a GE 3 T MR750W scanner equipped with a 24-channel head coil. For each participant, we selected the T1 BRAVO (TI/TR/TE = 450/8.48/3.23 ms, flip angle = 12°, FOV = 256 × 256 mm², matrix = 320 × 320, slice thickness = 0.8 mm, 230 slices, NEX = 2), field map (2D gradient echo, TR/TE1/TE2 = 500/4.9/7.4 ms, flip angle = 60°, FOV = 240 × 240 mm², matrix = 64 × 64, slice thickness/gap = 3.0/0.2 mm, 48 slices, magnitude and phase images), and multi-echo fMRI (ME-fMRI) (gradient echo EPI, TR/TE₁/TE₂/TE₃ = 2000/17/39.2/61.4 ms, flip angle = 70°, FOV = 224 × 224 mm², matrix = 64 × 64, slice thickness/gap = 3.5/0.0 mm, 48 slices, 210 dynamics) data.

To avoid the effects of the tDCS session on the neural activity, we selected only the data taken prior to the tDCS session. However, since the tDCS session was performed inside the MRI scanner immediately after the selected fMRI scan, all subjects wore MRI-compatible tDCS electrodes (Soteric Medical) attached to the head, with the anode placed at the left dorsolateral prefrontal cortex (DLPFC) and the cathode placed at the right supraorbital. No artifacts were observed in the fMRI images owing to the presence of these electrodes.

The fMRI task was performed in E-Prime ⁵ and synchronized with the fMRI scan using a scanner trigger. To allow the scanner to reach steady-state signals, the task started 6 s after the trigger.

During the fMRI scan, the participants performed an EFT task. EFT refers to the ability to imagine events that occur in one's personal future (Schacter et al 2018). EFT is closely related to episodic memory (thinking of one's own past) and uses a neural network comprising the default mode network (DMN): the medial temporal lobe (MTL) including the hippocampus, parahippocampal regions and amygdala, the posterior cingulate cortex (PCC), medial prefrontal cortex (MPFC), and lateral temporal (LTC) and parietal cortex (LPC). Additional EFT activations have been reported in the DLPFC, inferior parietal lobe (IPL), and ventromedial prefrontal cortex (vmPFC).

During the EFT task, a series of 18 common Dutch words was presented on a screen, and the participants had to imagine an associated future event that could happen in their own lives. The participants were instructed that the imagined scenarios must be new, non-repetitive, realistic, limited in time and space and imagined from a first-person perspective (e.g. presented word: car; imagined event: driving to a friend next weekend). Once an event was imagined, the participants had to push a button with their right hand. Each trial lasted 20 s with 4 s rest in between. The presented words were chosen randomly from a list of 72 words per subject.

All fMRI scans were visually checked for artifacts and manually oriented according to the anterior-posterior commissure line with the image origin set in the anterior commissure.

The data were preprocessed using a Nipype pipeline (Gorgolewski et al 2011) based on SPM12 ⁶ and FSL (Smith et al 2004).

The three TE images were combined using a T2*-weighted weighting factor (w_i for TE_i with i∈[1:3]) (Kundu et al 2017):

$$\begin{eqnarray}&&{w}_{i}=\displaystyle \frac{T{E}_{i}.{e}^{(-T{E}_{i}/T{2}^{* })}}{\displaystyle {\sum }_{j=1}^{3}T{E}_{j}.{e}^{(-T{E}_{j}/T{2}^{* })}}\end{eqnarray} \tag{ 1 }$$

with T2* determined by an exponential fit to the TE data from the first dynamic.

Second, the first 3 dynamics were removed as dummy scans. Using the magnitude and phase images from the field map scan, a voxel displacement map (VDM) was created (SPM: FieldMap) and used during the realignment and unwarping step (SPM: Realign & Unwarp) to correct susceptibility distortions (Andersson et al 2003). The SPM ArtRepair tool was used to create movement and frame-wise displacement (FD) plots to determine volume outliers. Based on a visual check of these plots, the data of 2 subjects were discarded from further processing owing to excessive head motion. To avoid adverse effects on the autocorrelation of the time series, no censoring of the volume outliers was performed (Siegel et al 2014, Liu 2016, Parkes et al 2018, Mascali et al 2021). The realigned time series were slice time corrected (SPM Slice Timing) using the slice times found in the json header file.

The T1 anatomical scan was cropped to the brain volume (FSL: robustfov) and the brain was extracted (FSL: BET). White matter (WM), gray matter (GM) and CSF maps were created during the segmentation step (FSL: FAST). The cropped anatomical scan and the tissue maps were coregistrated and resampled to the MNI template (SPM: Coregistration) and the fMRI scan was coregistrated to the resampled anatomical scan. Thereafter, the resampled anatomical scan was normalized (SPM: Normalization) and the normalization parameters were reused to normalize the tissue maps and the fMRI series. Finally, the normalized fMRI data were smoothed using a 6 mm Gaussian filter.

As motion regressors, we used 6 realignment parameters (3 translations and 3 rotations), their temporal derivatives, and squared regressors (24 regressors in total). Given the BOLD signals reported in WM (Gawryluk et al 2014, Grajauskas et al 2019), we used the temporal signals of 5 PCA components determined in the CSF using aCompCor, as physiological noise regressors.

For the ICA-AROMA analyses, the ICA components were determined using FSL Melodic. Using the script ⁷ of Pruim et al (2015), the noise components were identified as noise based on their location at the brain edges and in the CSF, the amount of high-frequency (> 1 Hz) content, and their correlation with noise regressors. In the first ICA-AROMA (ICA-AROMA-1) analysis only the 24 motion regressors were used as noise regressors. In the second ICA-AROMA analysis (ICA-AROMA-2), the 24 motion and 5 aCompCor signals were used as noise regressors.

As reference from the current practice, in the first analysis the temporal signals of the noise components determined by ICA-AROMA-1 and the five regressors from the aCompCor analysis were regressed out (nilearn: clean_img) from the fMRI signals (ICA-AROMA-1/regression). As alternative strategies, only the temporal signals of the noise components identified by ICA-AROMA-2 were regressed out (nilearn: clean_img) from the fMRI signals (ICA-AROMA-2/regression) in the second analysis and in the last analysis, the ICA components identified as noise were non-aggressively removed (FSL: fsl_regfilt) from the fMRI data (ICA-AROMA-2/removal).

The fMRI data were processed using SPM12. For each subject, a design matrix was built based on the EFT trial onsets and durations. The duration of each trial was defined as the time between the onset of the presented word and the button response. In the case of a missing response, the average response time for that participant was used.

The individual activation (positive model fit: EFT > 0) and deactivation (negative model fit: EFT < 0) contrast maps were thresholded using a voxel significance threshold of p(uncorrected)<0.001. A weighted overlap map (WOM) per contrast was created (Seghier and Price 2016) as

$$\begin{eqnarray}&&\begin{array}{l}WOM=\displaystyle \frac{1}{28}.\displaystyle \sum _{i=1}^{28}{w}_{i}\,with\,{w}_{i}\\ \,=\left\{\begin{array}{lll}1 & if & {p}_{i}\lt 0.001\\ {e}^{-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{p}_{i}-0.001}{0.0245}\right)}^{2}} & if & {p}_{i}\geqslant 0.001\end{array}\right.\end{array}\end{eqnarray} \tag{ 2 }$$

(w_i: weighting factor depending on the significance p_i of the contrast for subject i (i∈[1:28]))

One-sample t-tests were performed based on the individual contrast maps. To reduce the chance of false positive results, the group statistical maps were thresholded using a voxel significance threshold p(uncorrected)<0.001 combined with a WOM mask (WOM > 10%). This threshold was based on the Monte Carlo simulations of 1000 studies done on 28 subjects using RepSim ⁸ . Per study, 28 contrast maps were simulated and thresholded at p(uncorrected)<0.001, a WOM was determined and a 1-sample t-test was performed to simulate the group analysis. The resulting statistical map was thresholded at p(uncorrected)<0.001. The Monte Carlo simulations revealed that the chance of having a repeated false result in 3 subjects (∼10% of 28) in combination with a false positive group effect in at least 1 voxel is 0.18. On average, this corresponded to 1 ± 2 voxels per simulation.

To compare the fMRI results after applying the various denoising strategies, we calculated the global fMRI signal as the averaged signal over the entire brain. For this global signal we calculated the temporal signal-to-noise ratio (tSNR) and the fractional amplitude of low-frequency fluctuations (fALFF). The fALFF is defined as the power within the BOLD frequency range (0.08–0.1 Hz) over the total frequency power. The tSNR and fALFF values were compared pairwise between the three denoising strategies.

At last the overlap between the resulting group activation and deactivation maps were determined.

The fMRI series were split into 47 ± 18 components by Melodic. ICA-AROMA-1 identified 21 ± 15 components as noise components. ICA-AROMA-2 identified 22 ± 15 components as noise components. In 16 subjects, no differences were found in the classification of the components. In the others, 1 to 6 components were identified as noise by ICA-AROMA-2 but not by ICA-AROMA-1. Visual inspection of these latter components confirmed their classification as noise.

As can be seen in table 1, tSNR was significantly higher in ICA-AROMA-1/regression than in ICA-AROMA-2/regression and ICA-AROMA-2/removal. The tSNR in ICA-AROMA-2/regression was significantly higher than in ICA-AROMA-2/removal. The fALFF in ICA-AROMA-2/regression and ICA-AROMA-2/removal was significantly higher than in ICA-AROMA-1/regression.

In 1528 voxels, a significant activation was found for ICA-AROMA-1/regression, of which 97% overlapped with the activations found in ICA-AROMA-2/regression and 95% overlapped with those of ICA-AROMA-2/removal. 8723 voxels showed a significant activation in ICA-AROMA-2/regression, of which 17% overlapped with the activations found in ICA-AROMA-1/regression and 79% with those found in ICA-AROMA-2/removal. In ICA-AROMA-2/removal, 14890 voxels showed a significant activation, of which 10% overlapped with the activations found in ICA-AROMA-1/regression and 47% with those found in ICA-AROMA-2/regression. ICA-AROMA-1/regression revealed a significant deactivation in 864 voxels, while only in 4 voxels a significant deactivation was found in ICA-AROMA-2/regression and none in ICA-AROMA-2/removal. The 4 voxels of ICA-AROMA-2/regression overlapped with the deactivations of ICA-AROMA-1/regression. None of the deactivations found in ICA-AROMA-1/regression overlapped with the activations found in ICA-AROMA-2/regression or ICA-AROMA-2/removal. The overlap between the activation and deactivations maps is shown in figure 1.

Zoom In Zoom Out Reset image size