Multiband multi-echo simultaneous ASL/BOLD for task-induced functional MRI

Typical simultaneous blood oxygenation-level dependent (BOLD) and arterial spin labeling (ASL) sequences acquire two echoes, one perfusion-sensitive and one BOLD-sensitive. However, for ASL, spatial resolution and brain coverage are limited due to the T1 decay of the labeled blood. This study applies a sequence combining a multiband acquisition with four echoes for simultaneous BOLD and pseudo-continuous ASL (pCASL) echo planar imaging (MBME ASL/BOLD) for block-design task-fMRI. A multiband acceleration of four was employed to increase brain coverage and reduce slice-timing effects on the ASL signal. Multi-echo independent component analysis (MEICA) was implemented to automatically denoise the BOLD signal by regressing non-BOLD components. This technique led to increased temporal signal-to-noise ratio (tSNR) and BOLD sensitivity. The MEICA technique was also modified to denoise the ASL signal by regressing artifact and BOLD signals from the first echo time-series. The MBME ASL/BOLD sequence was applied to a finger-tapping task functional MRI (fMRI) experiment. Signal characteristics and activation were evaluated using single echo BOLD, combined ME BOLD, combined ME BOLD after MEICA denoising, perfusion-weighted (PW), and perfusion-weighted after MEICA denoising time-series. The PW data was extracted using both surround subtraction and high-pass filtering followed by demodulation. In addition, the CBF/BOLD response ratio and CBF/BOLD coupling were analyzed. Results showed that the MEICA denoising procedure significantly improved the BOLD signal, leading to increased BOLD sensitivity, tSNR, and activation statistics compared to conventional single echo BOLD data. At the same time, the denoised PW data showed increased tSNR and activation statistics compared to the non-denoised PW data. CBF/BOLD coupling was also increased using the denoised ASL and BOLD data. Our preliminary data suggest that the MBME ASL/BOLD sequence can be employed to collect whole-brain task-fMRI with improved data quality for both BOLD and PW time series, thus improving the results of block-design task fMRI.


Introduction
Functional MRI (fMRI) is a powerful, noninvasive tool to measure brain function. Two major contrasts used for fMRI are blood oxygenation-level dependent (BOLD) and arterial spin a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 labeling (ASL). While the BOLD fMRI signal is sensitive to magnetic susceptibility fluctuations caused by changes in blood oxygenation, it is also related to changes in cerebral blood flow (CBF), cerebral blood volume (CBV), and the cerebral metabolic rate of oxygenation (CMRO2). In contrast, ASL fMRI measures blood flow changes directly by magnetically tagging blood flowing into the brain. Therefore, CBF and BOLD are commonly used to study the hemodynamic response to neuronal activity and extract information regarding the role of the microvasculature in the brain [1][2][3][4][5].
Sequences have been developed to obtain ASL and BOLD contrast simultaneously by collecting ASL-and BOLD-sensitive echoes in one acquisition [1][2][3][4][5][6][7][8]. These sequences have been used to assess the contributions of CBF to the BOLD response [3,4], obtain calibrated BOLD and cerebrovascular reactivity (CVR) measurements [9,10], and reduce total imaging time by acquiring two image contrasts in one acquisition [5]. Furthermore, simultaneous acquisition of ASL and BOLD is an important tool for neurovascular coupling measures [6] and has been applied in research of aging and stroke [7,8]. However, ASL requires a tagging module and post-labeling delay (PLD) to allow tagged blood to flow into the brain. For pseudo-continuous ASL (pCASL), the recommended approach for ASL imaging [11], the suggested tagging time and PLD are each more than 1.5 s [11]. This approach results in long TRs and total ASL acquisition readout times that are severely limited by the short T1 relaxation of the tagged blood. In turn, it reduces the signal-to-noise ratio (SNR) and restricts the image resolution and total number of slices that can be acquired. Slices also experience varying PLD, which can affect accurate CBF quantification. Thus, most simultaneous ASL/BOLD studies acquire a 64×64 matrix, slice thicknesses >5 mm, and less than 20 slices, which can be inadequate when studying finer structures.
To address these issues, our group has recently developed a multiband (MB), multi-echo (ME) simultaneous ASL/BOLD sequence [12]. This sequence leverages two advanced techniques, MB imaging and ME imaging, to concurrently acquire whole-brain ASL and BOLD time-series with high tSNR. MB imaging, where multiple slices are excited and acquired simultaneously, can be used to increase spatial coverage and/or temporal resolution [13,14]. MB imaging has been developed and validated for task fMRI [15] and resting-state fMRI [14,16] and has been combined with ASL to acquire high-resolution blood flow images [17][18][19]. One major advantage of the addition of MB to ASL is a reduction in interslice labeling delay time differences. This is because MB imaging allows more slices to be acquired in fewer excitations compared with single band acquisitions. Fewer excitations also means T1-relaxation of the labeled blood is reduced. This leads to an increased SNR and more accurate CBF estimations [17][18][19]. MB-ASL has shown similar CBF estimation compared with single band ASL [19].
The MBME ASL/BOLD sequence also utilizes a multi-echo approach. In contrast to typical simultaneous ASL/BOLD sequences, which acquire only two echoes, our MBME ASL/BOLD sequence acquires four total echoes. A growing body of research indicates the image SNR, temporal SNR (tSNR), and BOLD sensitivity can be increased using ME echo-planar imaging (EPI) approaches [20][21][22][23][24][25]. The BOLD contrast is highest when the TE equals the T2 Ã relaxation of the tissue of interest. However, T2 Ã varies across the brain. This can be corrected for by acquiring several (>2) echoes and averaging them, weighted by the voxelwise T2 Ã . Combining echoes in this way has resulted in improved BOLD sensitivity [20,21,26]. Images can be further denoised using an automated ME independent component analysis (MEICA) approach [23,24,27], which can robustly separate BOLD and non-BOLD components. The MEICA technique classifies independent components as either BOLD or non-BOLD depending on whether their amplitudes are linearly dependent on TE or not dependent on TE, respectively. Components deemed unrelated to BOLD fluctuations are removed from the data.
In a previous study, the MBME ASL/BOLD sequence was used to acquire resting state fMRI data in a cohort of healthy control subjects [12]. Functional connectivity strength and network size was significantly increased following MEICA denoising of the BOLD signal. Robust connectivity was also observed in well-known brain networks using ASL-derived perfusion-weighted time-series. In the current study, we applied the MBME ASL/BOLD sequence to task fMRI. Four total echoes were acquired and combined to increase the image tSNR and BOLD sensitivity and remove non-BOLD signal from the BOLD data. Furthermore, we have modified the MEICA algorithm and introduce an automated ASL denoising technique to increase the tSNR of ASL data. Task fMRI data were acquired and compared between perfusion-weighted (PW), PW denoised (PWDN), second echo (E2), multi-echo combined (MEC), and multi-echo combined, denoised (MECDN) time-series.

Subjects
This study was approved by the Medical College of Wisconsin Institutional Review Board, and all subjects provided written informed consent before participating. 13 healthy adult volunteers (six males, seven females; mean age = 30.2 +/-8.5 years, ranging from 20-50 years) were recruited for this study. All subjects were right-handed. Subjects were asked to refrain from intake of caffeine before the MRI exam.

MBME ASL/BOLD sequence
A sequence diagram for the MBME ASL/BOLD sequence is shown in Fig 1. The sequence is described in detail in [12]. In short, the sequence consists of an unbalanced pCASL tagging module [28,29], followed by a PLD, a MB excitation, and ME gradient-echo EPI readout. Blipped-CAIPI [16] was also applied to reduce g-factor noise amplification caused by the sliceunaliasing in MB imaging, and is indicated by the z-gradient blips. In-plane acceleration was also implemented. Calibration repetitions were acquired at the start of the MBME ASL/BOLD acquisition using an interleaved approach and an MB excitation with phase-cycling [12,30]. These repetitions were used for the subsequent unaliasing. The last repetition in each acquisition was not tagged in order to obtain an M0 image, which is the equilibrium brain tissue magnetization used to normalize the subtracted PW maps for CBF quantification.
During the task-based fMRI scan, subjects performed a finger-tapping task. A blockdesign paradigm, consisting of four alternating periods of rest and bilateral finger tapping, was applied. Rest and tapping periods lasted 10 TRs (40 sec) each. During tapping periods, subjects were instructed to tap their thumb to the other four digits sequentially at their own pace.

Reconstruction
All image reconstruction was performed in Matlab (The MathWorks, Inc., Natick, MA, USA) and described in detail in [12]. First, Nyquist ghosting correction was performed using navigator echoes collected at the beginning of each excitation followed by echo separation. The calibration repetitions, acquired at the start of each MBME ASL/BOLD scan, were unaliased using a discrete Fourier transform and used to generate kernels for slice and in-plane unaliasing. A slice-GRAPPA algorithm [16] was implemented for slice unaliasing and applied separately for each echo. A traditional 1D-GRAPPA approach [31] was used to perform in-plane unaliasing following the slice-unaliasing procedure. Partial k-space was reconstructed using a homodyne method [32].

Preprocessing
The fMRI data processing pipeline is shown in Fig 2. Preprocessing was performed on each echo separately using AFNI [33,34] (https://afni.nimh.nih.gov/afni) and FSL (http://fsl.fmrib. ox.ac.uk/fsl/fslwiki) [35]. Anatomical segmentation and registration was performed in SPM. First, the MPRAGE image was segmented into gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF). The gray matter and white matter segmentations were then combined resulting in a brain-only MPRAGE image. This image was then transformed to Montreal Neurological Institute (MNI) space using a nonlinear registration algorithm. Motion was estimated for the first echo and those estimates were applied to the subsequent echoes using The sequence consists of an unbalanced pCASL labeling train, followed by a PLD, and finally an ME EPI readout. The first echo train was acquired after the acquisition of navigator echoes through the center of kspace. The phase was then rewound to the start of k-space, and the next echo train was acquired. This was repeated three times for a total of four echoes. MB imaging was implemented by replacing the single-band excitation pulse with a MB excitation pulse. Finally, blipped-CAIPI was utilized to shift the FOV of aliased slices and reduce g-factor penalties associated with MB imaging. Reprinted from [31] under a CC BY license, original copyright 2017. 3dvolreg in AFNI. Echoes were then coregistered to the anatomical MPRAGE image using epi_reg in FSL [36] with an affine registration with 12 degrees of freedom. Next the functional data was transformed to MNI space via applywarp in FSL using the transformation matrix output from the MPRAGE-MNI registration and normalization.

Multi-echo combination and denoising
All four echoes were combined using the T2 Ã -weighted approach [26,37]. First, the voxelwise mean across time of each individual echo dataset was used to estimate the signal immediately after excitation, S 0 , and the voxelwise T Ã 2 , T Ã 2ðfitÞ using log linear regression (Eq 1). The voxelwise T Ã 2ðfitÞ was then used to determine the weights, wðT Ã 2 Þ (Eq 2), which were used in a weighted summation of the echoes. TE n represents the n th echo time.
Following the echo combination, the data were denoised using the automated MEICA technique in afni and the meica.py plugin (v2.5). Default parameters were used except for the "daw" parameter, a weight use to control ICA dimensionality, which was raised from 10 to 20 Schematic showing the processing pipeline for the ASL and BOLD echoes. The first and second echoes were processed separately to yield the PW none and E2 data, respectively. Echoes were combined using a T2 Ã -weighted approach to generate the MEC dataset. This dataset was further denoised using MEICA, resulting in the MECDN dataset. No additional regression was performed in the GLM for the PW and MECDN datasets. Example activation curves and model fits are shown for the different datasets. to account for the low number of timepoints. This method is described in detail elsewhere [23,24,27] and classifies independent components as BOLD or non-BOLD based on whether their amplitudes are linearly dependent on TE [23,24,27]. Components deemed non-BOLD were removed from the data. For all cases, the tag-control modulation was identified as a non-BOLD component and removed from the data.

ASL processing
BOLD signal is known to contaminate the PW signal [38][39][40]. Liu et al showed there are two components of this BOLD contamination [38]. The first component is a multiplicative term related to the non-zero echo time of the sequence. The second component is a spurious component that can be reduced by a low pass filter. These filtering techniques however, are not perfect [38,40].
Perfusion-weighted denoising. Therefore, we have further modified the MEICA algorithm in order to denoise the first-echo images by removing both the artifact and BOLD components from the signal prior to generating the PW time-series (i.e. surround subtraction and high-pass filtering and demodulation). First, data from all echoes was low-pass filtered using a band-pass filter (3dTproject in afni) with a maximum frequency of 0.09 Hz. This frequency is below the label-control oscillation frequency (for TR = 4s, f = 1/(2 Ã TR) = 0.125 Hz), the main perfusion frequency (0.1125 Hz), and the second harmonic of the perfusion signal (0.1 Hz). The perfusion frequencies were determined by subtracting the task frequency (1/80s = 0.0125 Hz) from the label-control oscillation frequency. Filtering had the effect of removing PW signal oscillations from the data. Without low-pass filtering, one or more ASL components, characterized by these label/control oscillations, were identified. These components tended to be classified as artifacts and would be removed from the data, significantly reducing the tSNR of the subsequently calculated PW time-series. Next, the original MEICA algorithm was run, as described above, using the low-pass filtered echoes. This procedure combined the echoes using the T2 Ã weighted approach [26,37], extracted independent components, and classified those components as BOLD, artifact, and indeterminate. The BOLD and artifact components were then regressed out of the original, unfiltered first-echo data. A flow chart showing the PW denoising algorithm is shown in Fig 3. The PW denoising procedure relies heavily on the original MEICA algorithm; however, there are two key differences between the techniques. First, a temporal filtering procedure is applied prior to MEICA to remove perfusion-related signal oscillations. The original MEICA algorithm is then run without modification on this data. Second, whereas the original MEICA algorithm removes non-BOLD components from the combined ME time-series and keeps BOLD components, here we remove artifact and BOLD signals from the unfiltered, heavily perfusion weighted first-echo.
An additional denoising procedure was employed where only the artifact components were removed from the first-echo time-series. This was compared to the data where both artifact and BOLD components were regressed to determine whether regressing the BOLD signal provided an added benefit over regressing only the artifact components. Furthermore, to verify PW signal was not being removed by the denoising procedure, the mean and temporal standard deviation of the PW signal in gray matter were analyzed.
Perfusion-weighted time-series. A PW time-series was generated for both the nondenoised and denoised first-echo data using two common, related, approaches: by surround subtracting label and control images (SS) [41] and by high-pass filtering the data followed by demodulation (HD) [4,42]. Both these approaches are equivalent to demodulation followed by low-pass filtering, but differ in the filtering frequency. HD data was high-pass filtered with Echoes were first despiked, volume registered, and coregistered to MNI space. Each echo was then individually low-pass filtered at f 0.09 Hz. Echoes were then combined using a T2 Ã -weighted approach. This low-pass filtered, multi-echo combined dataset was fed into the MEICA algorithm, which extracted independent components and classified them as artifact, BOLD, or indeterminate. The BOLD and artifact components were regressed from the unfiltered first-echo data resulting in a denoised first-echo dataset. Surround subtraction and high-pass filtering followed by demodulation were performed on this data leading to denoised PW datasets. A GLM was employed on this data to determine activation.
https://doi.org/10.1371/journal.pone.0190427.g003 a minimum frequency of 0.09 Hz for the reasons described above. The modulated ASL component is present only above this frequency. Thus, high-pass filtering above this frequency separates the modulated ASL signal from the BOLD component present below 0.09Hz [4].

fMRI processing
The above procedures resulted in three BOLD and six PW datasets for each scan that underwent further processing for fMRI analyses: E2 (second echo, TE = 25ms), combined multiecho (MEC), MEC after MEICA denoising (MECDN), PW without denoising for the SS and HD data (PW ss,none and PW hd,none ), PW denoised with artifact components regressed for the SS and HD data (PWDN ss,art and PWDN hd,art ), and PW denoised with artifact and BOLD components regressed for the SS and HD data (PWDN ss,art+BOLD and PWDN hd,art+BOLD ). All data were blurred with a 4.5 mm full width at half maximum (FWHM) Gaussian kernel. For the E2 and MEC data, the six rigid-body motion parameters derived from the motion correction processing were regressed out of the data, and label-control oscillations were regressed out of the data by including a column of alternating −1s and 1s in the design matrix.

Statistical analysis
A general linear model (GLM) was used for the task-based fMRI analysis using 3dDeconvolve in afni. For the E2 and MEC data, in addition to the nuisance regression of motion parameters and the label-control sequence, the model included detrending with a third-degree polynomial. This was not necessary for the MECDN data where MEICA removed the low frequency drifts and motion-related artifacts. Individual activation maps were thresholded at an uncorrected p<0.01. For the MEICA denoised BOLD and ASL data, the number of degrees of freedom was reduced by the number of components regressed from the data on an individual subject basis. For the MECDN data, an average of 43.1 +/-7.7 total components were identified, and the number of regressed components ranged from 10 to 24, with a mean of 14.9 +/-4.6. For the PW data, an average of 32.9 +/-5.6 total components were found. For the PW data with artifact components removed, the number of regressed components ranged from ranged from 4 to 18, with a mean of 11.2 +/-3.9. Finally, for the PW data with artifact and BOLD components removed, the number of regressed components ranged from 13 to 33, with a mean of 22.8 +/-5.3. For the BOLD data, following 3dDeconvolve, a restricted maximum likelihood (REML) model (3dREMLfit) was used to model temporal autocorrelations in the data. This program uses an ARMA(1,1) to model the time-series noise correlation in each voxel. ASL data have been shown to not have significant temporal autocorrelation, so this model was not used for the PW data [40].
Group maps were calculated for each dataset (E2, MEC, MECDN, PW ss,none , PW hd,none , PWDN ss,art , PWDN hd,art , PWDN ss,art+BOLD , and PWDN hd,art+BOLD ) using a one-sample t-test. Group maps were thresholded at p<0.005 with cluster-size thresholding to correct for multiple comparisons. The minimum cluster size was determined using the recommended approach from Cox et al to reduce the family wise error rate [43]. First, spherical autocorrelation parameters (acf) were estimated using the residuals from the t-test and 3dFWHMx in afni. These parameters were then averaged for the BOLD and ASL data separately to create one set of acf parameters for the BOLD and ASL data. Parameters were then fed into afni's 3dClustSim program to determine minimum cluster sizes for α<0.05. The minimum cluster size for the BOLD data was 182 and for the ASL data was 131.
For each subject and dataset, the tSNR was computed on a voxelwise basis and defined as the mean signal divided by the standard deviation of the noise across the time-series. For the BOLD data, the mean tSNR was extracted from a whole-brain mask. For the PW data, the mean tSNR was extracted from the GM segmentation. In addition, the mean and maximum t-scores of the fMRI task were computed in activated voxels using an overlap mask that was created for each subject from voxels active in all BOLD and all PW datasets separately. Noise signal was defined as the residual between each voxel's best fit to the model and the signal itself. Statistical comparisons between mean values were made using a Bonferroni-corrected paired t-test to compare across datasets with significance set at p < 0.05.

CBF/BOLD relationship
For most simultaneous ASL/BOLD studies, the ASL signal is used to provide a better interpretation of the activation, instead of another way to compute activation. To that end, the relationship between the ASL and BOLD activation responses was evaluated. Mean signal was extracted for E2, MECDN, PW ss,none, PW hd,none , PWDN ss,art+BOLD , and PWDN hd,art+BOLD data from an overlap mask consisting of voxel active for all six datasets (uncorrected p < 0.01). The time series were converted to percent signal change by dividing by the mean of the residual time-series from the GLM. The four activation blocks for each dataset were then averaged. The activation magnitude was computed as the mean of the middle five time points from the activation. The ratio between CBF and BOLD activation magnitude was computed for the E2 and MECDN data vs. the PW and PWDN datasets respectively. The relationship between these ratios and baseline CBF was also computed.
Finally, CBF/BOLD coupling was computed for the same datasets. CBF/BOLD coupling was assessed by correlating the signals from the E2 and MECDN datasets to the PW and PWDN data respectively on a voxelwise basis using Pearson correlation. The correlation maps were thresholded at an uncorrected p<0.01 and the mean correlation values were extracted. The correlation maps were then averaged across subjects to create group CBF/BOLD coupling maps.

Data quality
Representative individual echo, MEC, MECDN, and mean PW ss,none and PW ss,art+BOLD images are shown in Fig 4. Image quality was better for the MEC and MECDN data compared to any individual echo. The PW ss,art+BOLD data was less noisy compared to the PW ss,none data. This is most noticeable in the white matter. Group-averaged tSNR is shown in Fig 5. Quantitative results are shown in Table 1 and mirrored the results from the previously published resting state MBME ASL/BOLD study [12,30]. Specifically, whole-brain tSNR significantly increased for the MEC data compared to the E2 data (p<0.001) and for the MECDN data compared to the MEC and E2 data (p<0.001). In addition, tSNR was significantly increased for the PWDN ss,art+BOLD data compared to the PWDN ss,art and PW ss,none data (p<0.001) and for the PWDN ss,art vs. PW ss,none data (p<0.001). The same trend was observed for the HD data with PWDN hd,art+BOLD tSNR increased compared to the PWDN hd,art and PW hd,none data (p<0.001) and for the PWDN hd,art vs. PW hd,none data (p<0.001).

Task fMRI
BOLD results from the group finger tapping fMRI analysis are shown in Fig 6. Fig 6A shows the results of the group analysis across subjects for the E2, MEC, and MECDN. The MECDN dataset detected the most activation, followed by the MEC and E2 datasets. Fig 6B shows mean signal curves from one representative subject in active voxels for the E2, MEC, and MECDN datasets. An ANOVA comparing datasets did not show significant differences.
PW results are shown in Fig 7. Additional activation was observed for the PWDN ss,art+BOLD data compared to the and PW ss,none data (Fig 7A, left); however, no appreciable differences were observed for the PWDN hd,art+BOLD data compared to the PW hd,none data (Fig 7A, right). In addition, increased activation was seen for the HD data compared to the SS data. Fig 7B  shows mean activation time-series from a representative subject for SS (left) and HD data (right). The SS denoised time-series appears markedly cleaner with less variance than the non- Quantitative results are shown in Table 1. In summary, the mean t-score in active voxels was significantly increased for the MECDN and MEC data compared to the E2 data and for the MEC data compared to the E2 data. Mean t-score was also increased for the PWDN ss,art+BOLD data compared to the PWDN ss,art and PW ss,none data and for the PWDN hd,art+BOLD data compared to the PW hd,none data. Maximum t-score in active voxels was not significantly different for the BOLD data, but increased with denoising for the PW data.

CBF/BOLD relationship
No differences were seen between processing schemes for the CBF/BOLD response ratio ( Fig  8A). In addition, a significant negative correlation was observed between the CBF/BOLD ratio and mean baseline CBF for all processing schemes except for HD, PW/E2 which trended toward significance (p = 0.08) (Fig 8B). For the SS technique, CBF/BOLD coupling significantly increased for the PWDN and MECDN datasets compared to the E2 and PW datasets. For the HD technique, CBF/BOLD significantly increased for the MECDN dataset compared to the E2 data (Fig 9). There was a slight, but non significant increase in coupling for the PWDN data for the HD method. These results are also observed in the group maps.

Discussion
The purpose of this work was to evaluate our recently developed MBME ASL/BOLD sequence for task-fMRI. To accomplish this, a finger-tapping task fMRI study was performed using this sequence. Four total echoes were collected. Combining the echoes resulted in increased BOLD sensitivity. The BOLD data was further cleaned using the MEICA denoising procedure. The MEICA technique was also modified to denoise the PW data. Motor activation was detected from the finger-tapping task using non-denoised and denoised PW and BOLD time-series. Data denoising resulted in increased sensitivity, tSNR, activation strength, and activation volume for the PW and BOLD datasets, and increased CBF/BOLD coupling.
One major advantage of our MBME ASL/BOLD implementation, compared with other typical simultaneous ASL/BOLD sequences, is that it gives us the ability to collect more than two echoes. A T2 Ã weighted combination of the echoes resulted in increased tSNR for the BOLD echoes. MEICA could also be used to denoise the data [23,24,27]. The automated MEICA denoising process relies on the TE dependence of BOLD-related ICA components to remove non-BOLD signals from the data. Using MEICA, the tSNR increased approximately 2.5 times  compared to E2, resulting in more activation and increased mean t-statistics. Previous work showed increased BOLD sensitivity for resting state fMRI using the MBME ASL/BOLD sequence by combining echoes and incorporating MEICA denoising [12,30]. Here, we confirm these increases in BOLD sensitivity extend to task-based fMRI as well.
The MECDN data showed increased activation, both in terms of the activation volume and activation strength, compared with the MEC and E2 data, as well as a higher overall tSNR. Though the temporal resolution of the MBME ASL/BOLD sequence was relatively low (TR = 4.0s), the increased tSNR from the additional echoes outweighed the drawbacks of the slight TR increase caused by these echoes. Murphy et al. showed that with increased tSNR, much shorter time-series are needed to detect activations for a certain effect size [44]. In line with this notion, our results also suggest the scan duration could be much shorter by using MECDN or PWDN art+BOLD for fMRI to detect activations with similar effect size as Robust bilateral activation was seen in the motor cortex, including the pre and postcentral gyrus, medial frontal gyrus, and cerebellum for all datasets. Increased activation was observed for the MEC data compared to E2 data and for the MECDN data compared to the MEC and E2 data. Activation was also observed in subcortical areas for the MECDN data. All maps were thresholded at p<0.005 and cluster corrected with a minimum cluster size of 182 voxels (α<0.05). (B) Average BOLD time-series extracted from a mask of voxels active for all BOLD datasets.
https://doi.org/10.1371/journal.pone.0190427.g006 conventional E2 or PW scans. The MBME ASL/BOLD sequence and subsequent denoising strategies provide a trade-off between scan time and activation strength. This is an important implication for future clinical applications.
Moreover, an increased area of activation was present in MECDN data in the group analysis. This included activation in subcortical areas not seen with the MEC and E2 data. Extensive bilateral activation was found in the insula and thalamus. The insula and thalamus have been shown to be part of the somatomotor network [45,46]. The thalamus is also known to be part of somatosensory pathways, which control our perception of touch. Finger-tapping necessarily involves touch, so activation in these deep gray matter regions should occur. Gradient echo EPI has typically shown poor performance in subcortical regions and requires higher sensitivity to detect significant signal changes [47]. Our results indicate MEICA denoising may be useful for identifying task-based subcortical activation due to the increased tSNR and resulting BOLD sensitivity. Kundu et al. found robust resting-state functional connectivity between subcortical and cortical structures following MEICA denoising, while no clear connectivity patterns could be detected by using standard denoising [24].
In this study, the MEICA algorithm was modified to denoise the PW data by regressing artifact and BOLD components from the first echo. BOLD signal can contaminate the PW signal, however, the BOLD and CBF time courses are strongly temporally coupled. Thus, there is the potential for true PW signals to be removed by removing BOLD components. To avoid this, all echoes were low-pass filtered below the label/control oscillation frequency and perfusion frequencies. This removed the majority of the perfusion contribution to the signal while , bilateral activation was observed in the motor cortex for the PW ss,none and PWDN ss,art+BOLD data. An increased activation area was seen for the PWDN ss,art+BOLD data compared to the PW ss,none data. The HD data (right) showed increased activation compared to the SS data, however no differences were seen between the denoised and non-denoised data the. All maps were thresholded at p<0.005 and cluster corrected with a minimum cluster size of 131 voxels (α<0.05). (B) Average SS PW signal from one representative subject (left) and average HD PW signal from the same subject (right). All PW signal was extracted from a mask of voxels active for all PW datasets. The denoised SS time-series appear less noisy with less variance compared to the non-denoised time-series. This effect is less apparent for the HD data. keeping the BOLD contributions. BOLD and artifact components were extracted using this filtered signal and then removed from the unfiltered first echo. To verify PW signal was not being removed, mean PW signal was analyzed. No differences were seen between the denoised and non-denoised data.
It is important to note this denoising technique is limited to block design fMRI studies where the PW signal is relegated to high frequencies, and thus a wideband lowpass filter can be used to remove the PW components. For event-related and resting state studies, this type of filter could potentially include PW components, which will then be removed from the data by the MEICA process.
In this study, both SS and HD techniques were used to extract the PW time-series. While these are two common approaches, they are fundamentally the same. Surround subtraction has been shown to be equivalent to demodulation followed by lowpass filtering, which itself is equivalent to highpass filtering followed by demodulation. We showed PW data quality benefits from the regression of BOLD components in addition to artifact components for both the The mean ratio of CBF to BOLD signal is plotted across subjects in active voxels (CBF/BOLD response ratio) for SS (left) and HD data (right). The response ratio was examined for non-denoised (PW/E2) and denoised data (PWDN/MECDN). No significant difference was observed between the non-denoised and denoised response ratios averaged across the middle five activation TRs (TR #s [8][9][10][11][12]. (B) CBF/BOLD response ratio plotted against baseline CBF for SS (left) and HD data (right). A significant negative correlation was observed between the CBF/BOLD ratio and mean baseline CBF for all processing schemes except for HD, PW/E2 which trended toward significance.
https://doi.org/10.1371/journal.pone.0190427.g008 SS and HD time series. PWDN ss,art+BOLD and PWDN hd,art+BOLD data had increased tSNR and mean t-score compared to PWDN ss,art and PWDN hd,art data respectively (the HD t-score trended toward significance). Both of these techniques have been shown to minimize the contribution of BOLD signal to the ASL signal [4,38,42]. These results indicate there may still be BOLD contamination of the PW signal using these methods.
High-pass filtering followed by demodulation was used as a stricter filtering approach to examine if the MEICA denoising technique was still effective. To provide a fair comparison between denoised and non-denoised HD data, the highpass filter frequency was the same as the lowpass filter frequency used for the MEICA denoising approach. By removing low-frequency BOLD and artifact components from the signal, MEICA denoising acts as an imperfect high-pass filter. Thus, in theory, highpass filtering and demodulating the data without denoising should provide similar, if not better, results compared to MEICA denoising as all frequency components below 0.09 Hz, including the inderminate components left by the denoising procedure, are removed. In fact, the HD method without denoising resulted in increased activation strength and volume compared to the SS method with denoising. However, highpass filtering and demodulating the denoised signal did lead to slight, but significant increases in mean t-score and tSNR at the individual level compared to highpass filtering and demodulating the non-denoised signal. SS filtering is not as strict as HD filtering (i.e. the frequency cutoff is lower). Therefore, MEICA denoising has a larger effect as less noise is inherently removed by the SS process, and there is increased BOLD contamination in the signal. Future studies could comprehensively examine the effect of filter cutoff frequency on PW data.
Finally, the relationship between CBF and BOLD activation was investigated. CBF/BOLD response ratios and CBF/BOLD coupling were compared between the denoised data (MECDN and PWDN) and non-denoised data (E2 and PW). The ratio between CBF and BOLD activation magnitude did not change with denoising indicating denoising does not change the quantitative relationship between CBF and BOLD. This has implications for techniques such as BOLD calibration, which uses ASL and BOLD data to estimate CMRO2 [48]. We also found a significantly negative relationship between the CBF/BOLD activation magnitude ratio and baseline CBF that did not change with denoising. This was likely driven by the increased baseline CBF leading to a reduced percent signal change for CBF activation as a larger CBF increase was needed to produce the same percent signal change as a lower baseline CBF. For the SS data, both BOLD and PW denoising with MEICA contributed to increased CBF/BOLD coupling. For the HD data, increases in CBF/BOLD coupling were driven by BOLD denoising with a limited effect from PW denoising.
This study was not without limitations. First, the number of subjects was relatively small. The purpose of this study was to determine the feasibility of using the MBME ASL/BOLD sequence to detect brain activation. Thus, the small subject size was justified. Additionally, only one scan was collected per subject, so a repeatability analysis could not be conducted. Future studies should look at the repeatability/reproducibility of MBME ASL/BOLD scans.
A relatively short PLD (1.5s) was employed for this study compared to the recommended PLD of 1.8 s for pCASL [11]. Shorter PLDs have the benefit of increased SNR, but can lead to intravascular artifacts if blood does not have adequate time to reach capillary-feeding small arteries. In addition, when the PLD is too short, the CBF response to activation can be overestimated. For example, if the arterial transit time is reduced in the activated state, the proportion of tagged blood in the voxel will be higher creating an artificially increased change in CBF [49]. This could impact studies where quantitative accuracy of CBF is necessary and should be examined in future studies. The PLD chosen here was a trade-off between TR, SNR, and the arterial transit time.
The issue of short PLD is especially important for interleaved MB acquisitions with MB slice packets that simultaneously cover inferior and superior slices. Intravascular artifacts were not detected in this study, and we were also able to consistently detect PW activation. Looking into locations for the simultaneously excited slices in MB pCASL may be worthwhile. Increasing PLD necessarily increases TR. The TR for ASL acquisitions is already long, owing to the long tagging and PLD segments of the sequence. However, we have shown that additional echoes collected in the MBME ASL/BOLD sequence can compensate for the longer TR.
Of note, most recently developed ASL sequences have utilized a segmented 3D readout, such as 3D GRASE, or stack of spirals approaches [11]. These techniques can also acquire whole-brain, high resolution ASL data; however, they are not compatible with multi-echo (ME) time-series acquisitions as they lead to very long TEs (for the 3D GRASE case) or TRs (for the segmented stack of spirals case). Background suppression (BS) was not used for the ASL scans in this study. There remains debate whether BS is appropriate for 2D approaches since BS can only be optimized for one slice. Furthermore, BS necessarily reduces BOLD SNR and tSNR. Future studies should examine the effect of BS on MBME ASL/BOLD.
In conclusion, we applied the MBME ASL/BOLD sequence to task fMRI. Motor activation was robustly detected using ASL and BOLD. In addition, the collection of more than two echoes allowed MEICA denoising to be applied to both the BOLD and PW data, resulting in increased tSNR, activation strength and volume, and CBF/BOLD coupling.