NORDIC was evaluated on a 3T segmented multi-echo VASO dataset (Huber et al., 2023, here we only used the first echo, TE = 12 ms), covering visual and motor cortices at 0.9 mm isotropic resolution (Figure 2A). The paradigm consisted of six 12-min runs of block-designed simultaneous visual stimulation and bilateral finger tapping (Figure 2B). We tested a 3T dataset rather than a 7T dataset because they are more dominated by thermal noise with a layer-resolution voxel size, making 3T VASO scans more challenging. If NORDIC improves VASO CNR at 3T, then advanced layer-fMRI methods become increasingly applicable, including at sites using 3T MRI systems. Furthermore, it presents a challenging test case given the CNR dependence on the success of low-rank denoising methods (Kay, 2022). Furthermore, it allowed for reasonable ground truth estimates at both the single-voxel level (after averaging across 72 trials, that supposedly mitigates most thermal noise) and the single-trial level (after averaging across an extensive activated area facilitated by combined visuomotor stimulation). Further parameters of this sequence were volume acquisition time for four sets of images (nulled and not-nulled combined, with two echoes each) = 5.2 s, inversion times (TI1/TI2) = 1072/1896 ms, phase partial Fourier = 6/8, matrix size = 216 × 216 × 12 (Huber et al., 2023), and skipped-CAIPI 6 · 1 × 3 z 2 (with 6 being the segmentation factor, 1 being the acceleration factor in the first phase encoding direction, 3 being the acceleration factor in the second phase encoding direction, and 2 being the CAIPI FOV shift in z). This notation is equivalent to segmentation factor 2, with in-plane GRAPPA 3, without GRAPPA z, and a prephased CAIPI in z of 2 (Stirnberg and Stöcker, 2021).

An additional 7T VASO dataset from a different subject was included to explore the effect of g-factor estimation strategies (described in Section 3.7). This dataset was acquired under the NIH-IRB (93-M0170, ClinicalTrials.gov: NCT00001360). Parameters were voxel size = 0.82 mm isotropic, TI1/TI2 = 1232/2386 ms, volume acquisition time (nulled and not-nulled combined) = 3.3 s, TE = 25 ms, k y -GRAPPA = 1, k z -GRAPPA = 3, phase partial Fourier = 6/8, and matrix size = 184 × 184 × 18.

The VASO contrast is cerebral blood volume (CBV)-weighted and is believed to have high spatial specificity toward microvasculature. CBV weighting is obtained by applying an inversion pulse followed by image readout around the blood-nulling time, and the fMRI signal will thereby decrease in proportion to increases in CBV. These nulled images are interleaved with regular BOLD images (not-nulled images), which can be used to correct for positive BOLD contrast counteracting the negative CBV contrast in the nulled images (see details in Huber et al., 2014). As detailed below, it is somewhat unclear how NORDIC should be applied to these alternating contrasts. Furthermore, the VASO contrast has inherent high spatial frequency and temporally dependent phase fluctuations, which is not necessarily compatible with the complex-valued processing in NORDIC (for more background on the basic working principle of NORDIC, see layerfmri.com/nordic; Moeller et al., 2021; Vizioli et al., 2021). Therefore, we first evaluated a comprehensive spectrum of different NORDIC implementation strategies. Questions that we wanted to address were:

We evaluated all 2 3 = 8 implementation combinations of NORDIC and compared them with the non-denoised version (noNORDIC). For all analyses, we used the publicly available MATLAB (Mathworks Inc.) implementation of NORDIC,¹ which comes in two flavors; “NORDIC.m” requires a noise-only volume and a g-factor map as inputs. Here, we focused on the other flavor, “NIFTI_NORDIC.m,” for which the g-factor is estimated from the data itself, and it can be readily applied to a broader range of datasets (but see Section 3.7 “Effect of g-factor estimation method” where the NORDIC.m flavor was explored).

NORDIC was applied to each run individually as the first preprocessing step. Aside from “gfactor_patch_overlap = 6” and “phase_filter_width = 10,” we used default parameters, including the default isotropic patch size: k × k × k, where k = numberOfTRs ⋅ 11 1 / 3 , rounded to the nearest integer. For all versions, time series were then motion-corrected with AFNI (Cox, 1996) using identical transformations estimated from the “combined-withNoiseVolume-complex” version. This was done to avoid interactions between motion estimation and denoising performance (although motion estimates were largely identical across denoising versions, results not shown here). The motion-corrected time series were used to calculate BOLD-corrected VASO using LAYNII (Huber et al., 2021). Trial-wise percent signal changes (PSCs) were then calculated from the ROI-averaged time series as the mean of task blocks minus the mean of rest blocks divided by the latter (the first two volumes were removed from each block to account for intermediate signal values during transition periods). ROIs were manually delineated based on structural contrast in motor and visual brain patches that showed strong BOLD activation (Figure 2D). This approach is chosen to minimize the potential bias of excluding noise-dominated weak voxels in brain areas that are expected to be activated as a whole.

Relevant performance metrics (tSNR, t-values, and PSC) were extracted from a large ROI (>7,000 voxels) including visual and motor cortices (Figure 2D upper panel). Voxel-wise cortical depth estimates were obtained between manually segmented white matter and cerebrospinal fluid boundaries using LAYNII (data upsampled to 0.23 mm isotropic resolution). Laminar profiles reflecting response amplitude as a function of cortical depth were then extracted from a small ROI placed in V1 (Figure 2D lower panel). Note that the response amplitude here was quantified as the raw difference (referred to as delta) between the mean of task blocks and the mean of rest blocks (the first two volumes of each block were still discarded). That is, we refrained from normalizing by the baseline signal to avoid division by values close to 0, which could be a concern for the blood-nulled VASO signal, especially toward CSF. This was not an issue in the non-laminar analyses where the baseline signal was averaged across all ROI voxels before division.

Differences in the average tSNR, t-value, and PSC between versions were evaluated using paired two-sided t-tests. Correction for multiple comparisons was performed using the Holm–Bonferroni method (Holm, 1979). The significance level was set at α = 0.05.

tSNR and t-value maps provide a measure of the degree of noise removal, with tSNR being a measure of temporal variance removed by NORDIC and t-values measuring whether the variance removed included relatively less task-locked signal of interest. As visible in Figure 3, these metrics were significantly increased in all denoised versions relative to noNORDIC. The complex-valued versions with the appended noise volume stood out with the largest gains, whereas the remaining versions were close (Figure 3). Even for these latter versions, NORDIC increased tSNR and the average single-voxel t-value by approximately 50%, which can be translated to, e.g., shorter scan durations, higher resolutions, or lower CNR. Interestingly, the “combined” versions yielded significantly higher tSNR than the “separate” versions, suggesting more noise removal. However, this came at the cost of increased signal loss, as evident from the reversal of this pattern in the t-values (i.e., the combined versions had significantly smaller t-values on average compared with the separate versions). This point is further underscored in the PSC evaluation detailed in “3.2 Response amplitude intactness after denoising.” Thus, the “separate-withNoiseVolume-complex” version seemed to provide the most efficient denoising based on tSNR and t-values alone.

The large tSNR values outside the brain for magnitude-only versions could be an effect of denoising in areas with low SNR where the noise distribution is non-Gaussian for magnitude data. In this Rician noise regime, the thermal noise distribution is a function of signal intensity (Gudbjartsson and Patz, 1995), and the spectral distribution of eigenvalues in the PCA decomposition is less distinct, whereby signal fluctuations are more likely removed. This would mean that the NORDIC assumption of Gaussian-distributed noise is not maintained, resulting in too high noise thresholds. In any case, it does not seem to be an issue inside the brain.

If NORDIC only removes thermal noise, the average PSC across ROI voxels should be unaltered (the ROI has over 7,000 voxels, and thermal noise should thus be practically eliminated after averaging, albeit a small bias will be present if the mean of the noise is non-zero). Figure 4A shows bar plots of PSCs computed from average time series across voxels in the ROI. All versions had a significantly reduced PSC compared with noNORDIC. In line with the conclusion from Figure 3C, this effect was on average significantly larger for combined versions compared with the separate versions (Figure 4B), suggesting that NORDIC should preferably be performed on nulled and not-nulled time series separately. For the separate versions, the effect was quite small, particularly for “separate-withNoiseVolume-magnitudeOnly,” which had a significantly lower PSC reduction than each of the other versions (Figure 4B). Since the patch size depends on the number of timepoints, we considered the possibility of the patch perhaps being over-sized for combined versions. However, the bias prevailed even when the patch size was set to something smaller (28×28×1, N = 784) than for the separate versions (N = 1728) (Figure 4C).

With the above approach to assess bias, signal removal could be masked if the PSC is biased toward 0 for both positively and negatively responding voxels. This motivated a replication of Figure 4 using a split-half approach across runs. Here, we only included within-ROI voxels showing a somewhat clear positive response (t-values >2 based on the first three runs). In line with observations in Figure 4, all versions showed a consistently reduced PSC (across the last three runs) compared with noNORDIC, with the largest bias in the combined versions (Figure 5). For the separate versions, the PSC bias appeared similar to or even reduced compared to Figure 4, suggesting that, also when accounting for this masking effect, NORDIC only minimally biased response amplitudes, particularly in the “separate-withNoiseVolume-magnitudeOnly” version.

In addition to testing for changes in response amplitude, we further evaluated whether NORDIC altered the temporal structure of single-trial responses (Figure 6) by plotting the voxel-averaged PSC as a function of trial number. The versions with the largest bias in PSC based on Figures 4, 5 naturally had a reduced mean. However, the temporal pattern of trial variability seems largely unaltered, as supported by the correlation coefficient between NORDIC and noNORDIC trial time series being >0.95 for all versions. If NORDIC altered the variability across trials (reduced variability is expected at the single-voxel level but not after averaging across >7,000 voxels), it would bias t-values (on top of the bias induced by response amplitude alterations), but that was not the case. The non-thermal “true” trial variance in response amplitude, originating, for example, from differences in attention, movement kinematics, and compliance, thus appears to be well-preserved.

To test whether NORDIC compromised spatial resolution, we computed an image of the mean nulled difference time series (first run), i.e., NORDIC minus noNORDIC, reflecting what is removed by NORDIC (Figure 7A). It should be noted that while all versions seem to have some structure in the periphery, the difference images within the brain look like unstructured random noise scaled by the g-factor as expected if NORDIC did not introduce blurring (Figure 7B).

Finally, we evaluated the spatiotemporal signal intactness after NORDIC based on laminar analysis. Figure 8A shows laminar profiles of each version—if denoising only removed thermal noise, NORDIC profiles should converge with the noNORDIC profile (after averaging 72 trials) in terms of response amplitude and shape, but with smaller error bars. This seems to be fulfilled for all versions (except the “combined, with noise-vol, complex” version), suggesting that robust estimation of laminar profiles could be obtained with fewer trials when using NORDIC. This is illustrated in Figure 8B (for the “separate-withNoiseVolume-magnitudeOnly” version), showing profiles with different degrees of trial subsampling. Subsampled noNORDIC profiles start breaking down already at 24 trial averages, whereas NORDIC profiles seem quite stable even at 12 averages. Furthermore, the observation of preserved profile shapes supports the conclusion in Section 3.4 that denoising did not result in noteworthy spatial blurring of the data.

In summary, the version that seemed to best preserve the signal was the “separate-withNoiseVolume-magnitudeOnly” (wrapper for running this version can be downloaded here: https://github.com/LasseKnudsen1/NORDIC-VASO), closely followed by the “separate-withoutNoiseVolume” versions, none of which were the ones providing the largest noise removal. This highlights that denoising performance should not be evaluated solely from the amount of noise removal; continued careful evaluation of signal intactness is critical. Importantly, more noise removal is often followed by more signal removal, thus representing a trade-off between the two (Kay (2022) for an excellent illustration), which we will look into in the following section.

Based on the analyses above, we prefer the “separate-withNoiseVolume-magnitudeOnly” version as it seemed to be the best performing in terms of removing noise while keeping the signal intact. Accordingly, we chose to move on with this version for the next analysis. Here, we denoised the data with different values of NORDIC’s “factor_error” parameter. The default value is 1, and it simply scales the originally determined singular value threshold. That is, if “factor_error = 1.2,” the threshold is increased by 20% and more components will be removed, resulting in more aggressive denoising. The risk, of course, is a speculated increased likelihood of concurrently removing the signal of interest or imposing temporal signals. Figure 9A illustrates this by showing the ROI-averaged PSC as a function of the scaling parameter. Increasingly more signals might be speculated to be removed as we increase the scaling factor, but the overcompensating noise removal (possibly also including physiological components at higher scale factors) leads to robust sensorimotor activation in the activation maps (t-values), which was not clear at the default scaling of 1 (some of which may also reflect spatial blurring at higher scale factors). To further emphasize the point of adverse signal removal, we computed delta activation maps from the difference time series (noNORDIC minus NORDIC), i.e., the part removed by NORDIC. Figure 9B shows how the difference time series (computed from not-nulled data as explained in the caption) contains increasingly more signals when the scaling factor grows.

With this parameter, we can balance the trade-off between noise and signal removal as desired. For example, it may be useful to push the scaling factor for functional localizers where we do not care so much if the response amplitude is accurate; we just want to find the right region—in this way, the acquisition time spent on localizer runs could be dramatically shortened, and one has more time to do the actual experiment. Another example would be the acquisition of a distortion-matched T1-weighted anatomical reference computed from a VASO time series, where alterations of the signals do not really matter as long as the anatomical structure is intact. The right choice of picking an optimal threshold value in these speculated application cases might require extra care, validation, and training, comparable to the challenges in ICA denoising of fMRI data.

The fact that the PSC is already slightly but significantly reduced at a scaling factor of 0.9 compared with noNORDIC (paired t-test across trials, t (71) = −2.94, p = 0.004) may reflect an important issue for all low-rank denoising methods: in the presence of a noise level that is comparable to the signal magnitude, the estimated principal components will always deviate from the noise-free scenario. That is, estimated components will be “corrupted” to some extent such that thermal noise and signal end up being spread out across all components (Kay, 2022; Olesen et al., 2023). In noisy experimental environments, the effect is always there, and it only becomes a matter of its extent. In many cases, it may be negligible, and a lot of noise removal at the cost of a very small bias will be well worth it, especially in cases where no activation is detectable without denoising.

With respect to the 3T data investigated here, we do find it promising that even for the present dataset where CNR is fairly low (3T VASO at submillimeter resolution), the signal reduction appears to be quite small at the default scaling factor and when the implementation of NORDIC is tailored for VASO. This was also the case in two additional VASO datasets with different field strength (7T), sequence parameters, tasks, number of volumes, participants, etc. The results from one of these datasets are shown in the section below, and the other is from the findings of the study by Huber et al. (2017) (results not shown here, but data can be downloaded at: 10.6084/m9.figshare.26819422).

The MATLAB implementation of NORDIC can either use the g-factor map provided by the scanner (NORDIC.m), which was computed here by the vendor’s reconstruction in IcePat, or the g-factor can be estimated from the input time series (NIFTI_NORDIC.m). The former should be closer to the ground truth and may therefore be preferable. To test this for VASO, we ran a similar analysis as presented in Sections 3.1–3.5, but this time on a 7T dataset (0.8 mm isotropic, six runs of bilateral finger tapping (two of which were discarded due to motion), 12 min each), where scanner-provided g-factor estimates were available.

To test whether it is worthwhile to obtain a g-factor map for each run (which can be impractical and time-consuming) or whether it is sufficient to reconstruct a single map at the beginning of the session, we ran NORDIC either with the same g-factor map for all runs (first run) or with run-specific maps. The current version of NORDIC.m can only be applied to complex data and with appended noise volumes. Thus, to most directly evaluate the effect of actual vs. estimated g-factors, we also ran NIFTI_NORDIC.m with appended noise volumes and on complex data. To test whether these versions were preferable compared with the version recommended based on Sections 3.1–3.5, we additionally ran NIFTI_NORDIC.m with appended noise volumes but on magnitude-only data. All versions were denoised on nulled/not-nulled time series separately as deemed optimal based on Sections 3.1–3.5.

All results (Figure 10) were extracted from a large ROI covering the hand knob of M1 in both hemispheres (ROI for laminar profiles only covered the hand knob in the left hemisphere). First, tSNR and t-value maps, PSC evaluations, and laminar profiles appear highly similar between using the same or run-specific g-factor maps, suggesting that one does not necessarily have to bother with computing g-factor maps for every run. This was the case despite large across-run motion, which resulted in noticeable g-factor differences across runs (Figure 11). Second, the same can be said when comparing NORDIC.m versions with the complex-valued NIFTI_NORDIC.m version, suggesting that denoising is not wildly sensitive to (small) inaccuracies in g-factor estimates (although some local differences can be seen, for example, in the tSNR maps). Third, all these versions had signs of signal alteration, as visible in the bar plots. This problem was largely gone for the NIFTI_NORDIC.m “separate-withNoiseVolume-magnitude-only” version (the “winner” from Sections 3.1–3.5), but that version also had less noise removal (~50% tSNR increase compared to >100% increase for the complex versions).

VASO-specific findings were:

Potential points in future studies to evaluate, optimize, and extend NORDIC: