Leveraging multi-echo EPI to enhance BOLD sensitivity in task-based olfactory fMRI

Functional magnetic resonance imaging (fMRI) is one of the most common imaging techniques used to non-invasively study both localized neural activation and whole-brain functional connectivity in humans. Blood-oxygenation-level-dependent (BOLD) fMRI measures brain activity by detecting local magnetic susceptibility fluctuations caused by changes in deoxyhemoglobin concentration. Gradient echo echo-planar imaging (GE-EPI) is the dominant BOLD imaging strategy due to its high temporal signal-to-noise ratio and rapid acquisition capabilities (Budde et al., 2014). However, artifacts caused by the local inhomogeneous magnetic field surrounding high-susceptibility regions present a major technical challenge in utilizing GE-EPI to measure BOLD signals. Olfactory regions, including the orbitofrontal cortex (OFC), amygdala, piriform (olfactory) cortex, and entorhinal cortex, are particularly affected by susceptibility-induced artifacts and signal drop-out due their location near the ethmoid sinuses (air-filled cavities in the bone surrounding the nose). These artifacts limit the applicability of GE-EPI-based BOLD fMRI in pursuit of understanding the human olfactory system (Catani et al., 2013). Moreover, respiratory artifacts, primarily caused by local field shifts during inhalation and exhalation that are highly correlated with event onsets in tasked-based olfactory-related fMRI, further compromise the sensitivity of olfactory fMRI (Glover & Law, 2001).

Several techniques have been proposed to overcome these challenges. Acquiring scans at a tilted angle allows for partial compensation of local susceptibility gradients (Deichmann et al., 2003; Weiskopf et al., 2006). Alternatively, pulse sequences other than GE-EPI, such as T₂-Prepared BOLD fMRI and balanced steady state free precession (bSSFP) (Miao et al., 2021; Parrish et al., 2008), are employed.

Recent work on multi-echo EPI (ME-EPI), which acquires multiple gradient echoes at each imaging timepoint, has shown potential to improve BOLD sensitivity and efficiency in resting-state fMRI (Gonzalez-Castillo et al., 2016; Lynch et al., 2020). In conventional single-echo EPI (1E-EPI), a single volume is acquired after each RF excitation at a specific echo time (TE), which determines both image signal intensity as well as BOLD contrast (Donahue et al., 2011; Gati et al., 1997; van der Zwaag et al., 2009). At a shorter TE, the image is more proton-density weighted, has less BOLD contrast, and less severe signal dropout in regions with high local susceptibility gradients. Conversely, with a longer TE, the image becomes more T₂^*-weighted, has greater BOLD contrast, and more severe signal dropout. The BOLD contrast-to-noise ratio (CNR) is expected to vary with TE, peaking at the T₂^* of each voxel (Donahue et al., 2011; Fera et al., 2004; Gati et al., 1997; Graham et al., 1996; Triantafyllou et al., 2011; van der Zwaag et al., 2009). However, various brain regions exhibit distinct T₂^* values, resulting in varying optimal TEs, making it impossible to optimize TE across the whole brain in a single-echo acquisition.

Unlike 1E-EPI, ME-EPI acquires multiple images across a range of TEs after each RF excitation, resulting in enhanced functional contrast by combining the echo images with different weights (Poser et al., 2006; Posse et al., 1999). The commonly used T₂^*-weighted echo combination method assigns the heaviest weight to the echo that is presumed to have the greatest BOLD CNR, with weighting performed separately for each voxel. By applying a T₂^*-weighted echo combination with ME-EPI, BOLD signal can be maximized at the voxel level while maintaining the temporal resolution needed for fMRI. Moreover, by leveraging the expected exponential signal decay across TEs for true BOLD contrast, noise reduction techniques can be employed to reject artifactual signals and thereby improve the BOLD contrast-to-noise ratio (DuPre et al., 2021; Gonzalez-Castillo et al., 2016; Kundu et al., 2012, 2013, 2017; Van et al., 2023).

In this study, we implemented and tested an ME-EPI acquisition protocol and analysis pipeline specifically designed for task-based olfactory-related paradigms by incorporating shorter TE values in an effort to improve sensitivity in ventral brain regions with high static susceptibility. Using this approach in comparison to a standard 1E-EPI acquisition, we demonstrate superior sensitivity to olfactory task activation, resulting in a significant sample size reduction for detecting group task activation during a representative olfactory three-alternatives forced choice task (3AFC).

Twenty-nine subjects (13 women; aged 22–30 years, mean: 26, SD: 3.15) provided informed consent as approved by the University of Pennsylvania Institutional Review Board (#827217). One subject was excluded from the study due to excessive motion (number of excluded volumes, as defined in Statistical Analysis, exceeded 2% of the number of total volumes) during fMRI scanning session, leaving us with a final sample of n = 28. All subjects reported being right-handed nonsmokers with no history of significant neurological/psychiatric disorder or olfactory dysfunction.

Stimuli consisted of two odorants: lemon oil extract (12.3% v/v, Sigma-Aldrich, MO, USA) and benzaldehyde (0.42% v/v, Sigma-Aldrich, MO, USA). Both odorants were diluted in mineral oil (Sigma Aldrich, MO, USA). Mineral oil alone was used as a control stimulus. The odor stimuli were delivered through a custom-built computer-controlled air-diluted 12-channel olfactometer with a flow rate of 4 L/min controlled using the PsychToolbox package (Brainard, 1997) and MATLAB (Version 2020b, The Mathworks Inc., Natick, MA, USA). Medical-grade air (Airgas, Radnor, PA, USA) was first routed through two computer-controlled 6-channel gradient valve manifolds (NResearch Inc., West Caldwell, NJ, USA), regulated by a mass flow controller (MFC) (Alicat Scientific Inc., Tucson, AZ, USA). Subsequently, it passed through the headspace above 2 mL of the mineral oil-diluted odorant before being delivered to the participant through a fluorinated ethylene-propylene (FEP) plastic tube (U.S. Plastic Corp., Lima, OH, USA) placed directly below the nose.

On Day 0, subjects participated in a behavioral testing session during which odor stimuli were presented outside the MRI environment (Fig. 1a). Participants were introduced to the odors and asked to label them as well as to rate the perceived odor intensity, via the general Labeled Magnitude Scale (gLMS) (Green et al., 1993, 1996). Valence (likeability) was also assessed via the Labeled Hedonic Scale (LHS) (L. V. Jones et al., 1955; Lim et al., 2009; Schutz & Cardello, 2001; Stone et al., 2020), as well as familiarity [Visual Analog Scale (VAS) (Dourado et al., 2021; Flint et al., 2000)]. Ratings were analyzed in MATLAB (Version 2020b, The Mathworks Inc., Natick, MA, USA), using three one-way ANOVAs to evaluate whether there were any differences among stimuli.

Two to 5 days after the behavioral testing session, the same participants took part in an fMRI scanning session, completing six runs, each of which contained 24 trials. On each trial, participants sniffed one of the three possible odorants (lemon oil extract, benzaldehyde, or mineral oil) and were asked to identify the odor using a 3AFC task with the labels they provided earlier during the behavioral testing session. Following the 3AFC task, they were also asked to rate the intensity of the odor. Each odor was presented eight times per run in a random order. The experimental paradigm as well as the sequence of events within a given trial is shown in Figure 1.

All magnetic resonance images were acquired using a 3 T scanner (MAGNETOM Prisma, Siemens Healthineers, Germany) using the vendor’s 64-channel head and neck coil. T₁-weighted structural images were acquired using ME-MPRAGE (TE = 1.69/3.55/5.41/7.27 ms, TR = 2530 ms, TI = 1100 ms, FA = 8°, 1 mm isotropic).

1E-EPI and ME-EPI fMRI data were acquired using the sequence parameters outlined in Table 1. For both protocols, repetition time (TR), field of view (FOV), number of slices, resolution, flip angle (FA), echo spacing, and multiband factors are identical. To optimize the BOLD signals in OFC regions, the field of view was tilted approximately 25° to the anterior commissure-posterior commissure line (rostral > caudal) (Deichmann et al., 2003). The two sequence types were administered in an interleaved fashion in each participant, with the order being counterbalanced across participants. Of note, for 1E-EPI the echo time was set at 22 ms, which is approximately the same as the second echo time out of five echoes in the ME-EPI sequence (see Table 1). This TE aims to balance signal dropout and BOLD sensitivity for olfactory-related regions and is consistent with other olfactory-related fMRI studies (Bao et al., 2019; de Groot et al., 2021; Raithel et al., 2023; Shanahan et al., 2021; Zhou et al., 2019).

Additionally, a pair of spin-echo EPI (SE-EPI) images with two phase-encode blips are acquired for both fMRI acquisition protocols (1E-EPI and ME-EPI) separately with identical FOV, number of slices, resolution, FA, echo spacing, acceleration, and readout bandwidth compared to their respective fMRI acquisition. These pairs of images are then used to estimate inhomogeneities in the static magnetic field (B₀) during preprocessing using a similar method described in (Andersson et al., 2003) and implemented in FSL (Smith et al., 2004) as part of the fMRIPrep package (Esteban et al., 2019; Gorgolewski et al., 2011).

fMRI and structural MRI data were preprocessed in each subject’s native space using fMRIPrep 21.0.0 based on Nipype 1.6.1 (Esteban et al., 2019; Gorgolewski et al., 2011). All images were then registered to ICBM 152 Nonlinear Asymmetrical template version 2009c (MNI-ICBM2009c) (2.4 mm isotropic resolution, consistent with fMRI acquisitions) before further processing (Fonov et al., 2009, 2011). Spatial smoothing was applied to all data using a Gaussian kernel of 8 mm full width at half maximum (FWHM). A detailed description of these preprocessing steps can be found in Supplementary Materials under the Detailed fMRI Data Preprocessing section.

Multi-echo EPI produces separate images for each echo time at each timepoint in the functional imaging run. At each timepoint, these images for the multiple echo times were combined using the T₂^*-weighted combination method (ME-WC) (Posse et al., 1999). For each voxel, T₂^* was estimated by fitting the Eq. (1) using signals from all five echoes,

S₀ is the initial signal magnitude without any T₂^* decay, and S is the signal magnitude acquired at t = TE.

The images from the five echoes are then combined voxel-wise using a weighted average (ME-WC), and the weights for each echo (W_TE) are calculated using Eq. (2) based on estimated T₂^*.

Independent components analysis (ICA) is a common strategy to decompose fMRI signals into distinct independent components, allowing the isolation of BOLD-related neuronal components from sources of noise (Beckmann, 2012; Caballero-Gaudes & Reynolds, 2017; Mckeown et al., 1998, 2003). However, a persistent challenge lies in the manual classification of these independent components, which is labor-intensive and can yield inconsistent results (Kelly et al., 2010). In an attempt to tackle this problem, ME-ICA classifies each component as either a non-BOLD artifact, such as respiratory artifact, or a BOLD-like signal based on their signal characteristics across different echo times (Kundu et al., 2012, 2013). Here, we specifically employed the TE-dependency analysis pipeline (Tedana 0.0.12) (DuPre et al., 2021) to preprocess multi-echo EPI data before registering them to the MNI template.

In the original TEDANA pipeline, a monoexponential model [Eq. (1)] was fit to the data at each voxel using nonlinear model fitting to estimate T₂^* and S₀ maps, using T₂^*/S₀ estimates from a log-linear fit as initial values. This approach generally performs well across most brain regions. However, in regions with significant signal dropout at longer echo times, particularly due to susceptibility artifacts in areas such as the orbitofrontal cortex and entorhinal cortex, this fitting method tends to overestimate T₂^*. Consequently, we replaced the TEDANA estimate of T₂* with our own estimates, generated by fitting the following model voxel-wise:

In this model, x₀ represents the thermal noise floor, which does not decay over TE and whose contribution to the overall signal increases in later echoes as the signal diminishes in high-susceptibility regions. This model was fit to the data by minimizing the 2-norm of the weighted residuals [Eq. (4)] via the Trust-Region-Reflective algorithm (Coleman & Li, 1996).

The weighted residuals are calculated as the difference between the estimated signals $S^i$ and the observed signals $Si$⁠, weighted by the observed signals $Si$ across all n = 1000 observations (200 volumes/run, 5 echoes/volume).

The framewise displacement (FD), calculated during fMRI preprocessing using the method described in (Power et al., 2012), and the onset of odor stimulus for each condition were extracted for each run. The coefficient of determination (R²) was then estimated by calculating the Pearson correlation coefficient between the FD and the odor onset for each stimulus condition convolved with the canonical hemodynamic response function (HRF) for each subject (Pearson, 1895; Siegel, 1956). To compare whether there is a group-level difference in movement parameters for different odor stimuli, two-sided paired t-tests were performed on the R², across all three stimulus conditions ([Lemon vs. Control], [Benzaldehyde vs. Control], and [Lemon vs. Benzaldehyde]) for both 1E-EPI and ME-EPI acquisitions.

The respiratory trace was recorded with a respiratory transducer (breathing belt) using a pressure sensor (BIOPAC System, CA, USA) affixed around the subject’s chest, and then preprocessed before being added as a nuisance regressor: the trace was smoothed using a moving average with a 250 ms window, high-pass filtered at 0.05 Hz with an infinite impulse response (IIR) filter, detrended, and standardized such that the mean of the signals is 0 and the standard deviation is 1 (Liu et al., 2024; Raithel et al., 2023). The preprocessed trace was also downsampled to 1/2.4 Hz to match the TR of the fMRI acquisitions and included as a nuisance regressor. Two additional nuisance regressors were introduced: the trial-by-trial sniff volume was defined as the amplitude difference between the peak and the trough within the cued sniff, and the sniff duration was defined as the time difference between the peak and the trough.

All fMRI data were analyzed in SPM12 using a general linear model (GLM) to estimate the main effects of the odor stimulation. Each experimental condition (control, lemon oil, and benzaldehyde) was modeled using an independent regressor for each run. This results in a 2 (multi-echo and single-echo) by 3 (control, lemon oil, and benzaldehyde) factorial design. The model included the following nuisance regressors: trial-by-trial sniff volume and sniff duration, both convolved with the hemodynamic response function (HRF), as well as the down-sampled breathing trace and 24 movement parameters (six movement parameters derived from spatial realignment, their squares, derivatives, and squares of derivatives). Additional nuisance regressors were introduced to exclude volumes with excessive motion, defined as those with at least one of six movement parameters derived from spatial realignment deviating by more than 6 standard deviations from the mean. In all GLMs, data were high-pass filtered at 1/128 Hz; the temporal autocorrelation was estimated and removed as an autoregressive AR(1) process globally (Friston et al., 2002).

Contrasts of interest were defined as [Odor > Control], [Lemon > Control] and [Benzaldehyde > Control]. Whole-brain analysis was performed to validate that odor-evoked activations were present in the data acquired through both protocols (1E-EPI and ME-EPI). Small volume corrections (SVC) and ROI analysis were separately performed on OFC, amygdala, piriform cortex, and entorhinal cortex, as suggested by previous studies highlighting the involvement of these regions in human olfactory processing (Anderson et al., 2003). Amygdala and piriform cortex were defined anatomically using a human brain atlas (Mai et al., 2016). Delineation of OFC and entorhinal cortex was based on the MNI-ICBM2009c atlas (Manera et al., 2020).

We used our ROI analysis results to estimate the sample sizes required to detect significant olfactory-related group fMRI activation in subsequent studies. A priori power analyses were performed with α = 0.05 and β = 0.80 using G*Power (Faul et al., 2009) for both ME-EPI and 1E-EPI data in the same ROIs described in the previous section.

We first examined the estimated T₂^* values for each olfactory-related region. As expected, signal dropouts are reduced in shorter echo times (Fig. 2a, TE₁ and TE₂), while T₂^*-weighted contrast is greater in later echo times (Fig. 2a, TE₃–TE₅). In particular, the medial and lateral OFC along with the entorhinal cortex exhibited notably smaller T₂^* values due to susceptibility artifacts (Fig. 2a & 2b), and in turn the echo-combined image more heavily weights the shorter echoes (Fig. 2c). Using these weights, the combined image (Fig. 2a, ME-WC) effectively mitigated signal dropout while retaining notable T₂^*-weighted contrast.

We next investigated the regions that were significantly activated when the odors were delivered under ME-EPI acquisition with the ME-ICA denoising method [Odor > Control] (Fig. 3a) (DuPre et al., 2021). Consistent with previous literature, we observed significant activation in the piriform cortex, amygdala, entorhinal cortex, insula, and OFC bilaterally (p < 0.001, uncorrected) (Anderson et al., 2003; Dikecligil & Gottfried, 2024; Gottfried et al., 2003; Howard & Gottfried, 2014; Patin & Pause, 2015; Plailly et al., 2008). We next examined the BOLD effects of individual odor stimuli, namely [Benzaldehyde > Control] and [Lemon > Control]. Despite statistical tests indicating no significant differences in perceived intensity (p = 0.38), valence (p = 0.45), and familiarity (p = 0.11) between benzaldehyde and lemon oil and no significant differences in motion among three experimental conditions ([Lemon vs. Control] p = 0.60, [Benzaldehyde vs. Control] p = 0.33, and [Lemon vs. Benzaldehyde] p = 0.63), our analysis revealed that the main effect [Odor > Control] was driven mainly by the contrast [Lemon > Control]. In comparison, the contrast [Benzaldehyde > Control] produced significant levels of activation only in the bilateral insular region (Fig. 3a).

The difference in BOLD response to these stimuli was not anticipated though not wholly unprecedented, as previous research has demonstrated that different odors can elicit distinct responses within olfactory brain regions even in the absence of differences in low-level stimulus features (Fournel et al., 2016; Howard & Gottfried, 2014; Savic, 2002; Savic et al., 2002). In the single-echo fMRI, a similar trend is observed, with more active regions being detected with lemon oil (p < 0.05, uncorrected, Fig. S3). To our knowledge, there is no literature in humans that directly elicits piriform cortex BOLD activation with a similar sample size for benzaldehyde. However, rodent studies have reported that benzaldehyde demonstrates BOLD activation in the amygdala, as well as the piriform and entorhinal cortices (Kulkarni et al., 2012). Accordingly, to compare and observe the differences between acquisition methods within olfactory regions, we will focus on the activation under lemon oil condition [Lemon > Control] for our further analyses.

Having validated that our ME-EPI acquisition protocol is suitable to capture odor-evoked BOLD responses, we next compared BOLD activation among different acquisition methods (1E-EPI vs. ME-EPI) and analysis pipelines (ME-WC vs. ME-ICA). A whole-brain analysis using either whole-brain multiple comparisons correction with false discovery rate (FDR) (Benjamini & Hochberg, 1995; Benjamini & Yekutieli, 2001; Yekutieli & Benjamini, 1999) or small volume correction (SVC) did not yield any supra-threshold voxels for 1E-EPI. However, the same whole-brain analysis using the T₂^*-weighted combined ME-EPI (ME-WC) highlighted activation in all ROIs (piriform cortex, amygdala, OFC, and entorhinal cortex) (Fig. 3b). Moreover, by using ME-ICA denoising methods (DuPre et al., 2021), we observed a further increase in both the number of activation voxels and peak Z values in all ROIs (Table 2).

We next examined potential reasons for the observed differences. Specifically, we hypothesized the improvement of the BOLD sensitivity in the case of the ME-EPI data was due to larger BOLD signal changes (i.e., greater sensitivity), as we performed a locally weighted combination of all echoes for each individual voxel instead of using one single TE across the brain. BOLD signal changes could be calculated by contrast estimates from the GLM. Therefore, we performed a paired t-test on contrast estimates [Lemon > Control] between 1E- and ME-EPI with ME-ICA [ME-ICA > 1E-EPI] among all ROIs. Only piriform cortex (p = 0.023) and lateral OFC (p = 0.006) showed greater estimates for the contrast [Lemon > Control] for ME-ICA compared to 1E-EPI.

An alternative explanation for the observed increase in BOLD detection for ME-EPI versus 1E-EPI may lie in the reduction of the background noise due to non-BOLD signals (i.e., greater specificity). Accordingly, the signal acquired with ME-EPI exhibits less fluctuation over time, resulting in lower variance and, in turn, increasing Z-scores. To test the hypothesis that ME-EPI signals have lower variance, we performed a Levene’s test (Levene, 1960). This test is commonly employed to assess whether two or more groups have equal variances in a dataset (Chiarello et al., 2009; Farinha et al., 2022; Im et al., 2010; T. B. Jones et al., 2010; Swarup et al., 2013). When using a whole-brain analysis to compare the signal variance between ME-ICA and 1E-EPI, we observed suprathreshold clusters (p < 0.05, FDR corrected) overlapping with the voxels activated for the [Lemon > Control] contrast (p < 0.05, FDR corrected) from ME-ICA (Fig. 3c). We then generated a mask using those activated voxels and found that 18.89% of the masked voxels reached the significance threshold under Levene’s test. This suggests that ME-ICA decreases the signal variance, thereby reducing non-BOLD noise and increasing BOLD specificity. We did not observe any suprathreshold clusters when comparing ME-WC and 1E-EPI.

Taken together, these tests suggest that our observed increase in BOLD detection is not only due to increased sensitivity to BOLD signals, but also the reduction of non-BOLD artifacts enhancing BOLD specificity.

So far, we have established that ME-EPI improved BOLD signal sensitivity by both enhancing BOLD signal and reducing non-BOLD artifacts. An outstanding question is to what degree BOLD sensitivity is improved when using our suggested method in the context of task-based olfactory fMRI. To quantify the change in BOLD sensitivity, we examined the contrast [Lemon > Control] in the brain regions central to olfaction, namely the piriform cortex, amygdala, lateral and medial OFC, and entorhinal cortex for both acquisition methods separately and then compared the statistical outcomes. Whereas activation in medial OFC did not reach statistical significance using data resulting from either acquisition method, all remaining regions show much lower p-values for ME-ICA compared to 1E-EPI (Table 3). Furthermore, to determine the minimum sample size required to observe significance in each region, post-hoc power analyses were carried out (α = 0.05 and β = 0.80) to estimate the effect of each acquisition method on the sample sizes needed to detect effects in similar olfactory task fMRI studies. ME-ICA led to a substantial reduction in the sample size, ranging from by half to one-fourth, compared to 1E-EPI, depending on the specific regions (Fig. 4).

Susceptibility artifacts pose a significant challenge when using GE-EPI for fMRI studies. The primary olfactory cortex not only suffers from some of the most severe static susceptibility artifacts within the human brain but is also affected by respiration-related artifacts. Here, we have developed and validated an ME-EPI protocol and analysis pipeline specifically tailored to olfactory task-based fMRI. In the context of a basic olfactory 3AFC task, we found that ME-EPI, alongside ME-ICA denoising, systemically improves BOLD sensitivity and reduces the required sample sizes needed to infer a significant effect at a group level compared to conventional 1E-EPI. This is achieved by both increasing BOLD signal sensitivity and reducing the contributions of non-BOLD signal changes using the additional denoising features ME-EPI provided.

The echo times used for ME-BOLD provided a range of contrasts across different TEs for each ROI of interest (Fig. 2a). When evaluating the T₂^*-derived weights calculated for the echo-combined image, we confirmed that high-susceptibility regions with short T₂^* values, namely, lateral and medial OFC, and entorhinal cortex, had higher weights assigned to the first two echoes. In contrast, the remaining ROIs had higher weights assigned to the remaining (and later) three echoes. As a result, the T₂^*-weighted combined (ME-WC) images were able to produce more uniform T₂^*-weighted contrast across the whole brain (Fig. 2a, ME-WC).

When investigating the region activated by odor onset using ME-ICA, we detected activation in all major primary olfactory cortex regions (piriform cortex, amygdala, entorhinal cortex, and OFC bilaterally), in addition to the insular region (Fig. 3a). These findings are consistent with previous literature (Anderson et al., 2003; Dikecligil & Gottfried, 2024; Gottfried et al., 2003; Howard & Gottfried, 2014; Patin & Pause, 2015; Plailly et al., 2008). Further analysis revealed that the activation was primarily driven by the lemon extract stimulus in comparison to benzaldehyde. Previous studies have demonstrated that different odor stimuli can activate distinct primary olfactory regions to different degrees (Fournel et al., 2016; Howard & Gottfried, 2014; Savic, 2002; Savic et al., 2002). In this study, we found that that benzaldehyde activates the olfactory system to a lesser extent than lemon oil, such that the observed activation did not reach significance after correcting for multiple comparisons in the whole-brain analysis.

When comparing the results of 1E-EPI and ME-EPI, we observed significantly increased sensitivity in both whole-brain analysis and ROI analysis, which are two commonly used univariate fMRI analysis methods. In the case of whole-brain analysis (Fig. 3b), we found that ME-EPI exhibited higher peak Z values and a greater number of activated voxels compared to 1E-EPI, as shown in Table 2. It is worth noting that even without using the ME-ICA denoising method, ME-EPI still outperformed 1E-EPI by performing a T₂^*-weighted combination of data from all echoes. This improvement was expected, as the weighting of data from each echo time was designed to enhance the sensitivity to BOLD contrast (Posse et al., 1999).

For ROI analysis and sequential power analysis, which were used to estimate minimal sample sizes, ME-EPI with the ICA denoising method demonstrated superiority over 1E-EPI, as shown in Table 3 and Figure 4. We attribute these improvements to two factors: the increase of BOLD signal intensity through T₂^*-weighted echo combination and the reduction of noise by leveraging signal evaluation across the echo times. To support these hypotheses, we conducted a series of statistical tests (Fig. 3c) to show the mean signals are significantly higher and the variance is significantly lower in ME-ICA acquisition compared to 1E-EPI.

While both ME-EPI and ME-WC overwhelmingly increased the statistical power for ROI analysis, we observed an unexpected decrease in statistical power for ME-WC in the entorhinal cortex. Previous studies have demonstrated that local B₀ field shifts affect T₂^* values near the ethmoid sinuses during respiration, resulting in fluctuations of T₂^* values across respiratory cycles (Barry & Menon, 2005; Raj et al., 2000, 2001; Van de Moortele et al., 2002). The T₂^*-weighted combination method used in our analyses adopted a fixed weighting approach throughout a single scan, potentially failing to account for these dynamic T₂^* fluctuations. However, our additional analyses showed that, although statistically significant in most ROIs, T₂^* values exhibit only subtle differences between inhalation and exhalation phases (Table S1). This, in return, demonstrates the robustness of the echo combination method against respiratory phase variations.

We also evaluated whether nonexponential signal decay contributed to suboptimal echo combination weightings, leading to the unexpected result. As shown in Figure S2, signals from most regions fit well with exponential decay, but regions with high susceptibility, such as the entorhinal cortex, medial, and lateral OFC, exhibited more pronounced nonexponential decay. Yet, this observation does not fully explain the unexpected decreases observed exclusively in the entorhinal cortex.

Furthermore, we examined the effects of ME-WC and ME-ICA on the distribution of statistical significance within each ROI (Fig. S1). We have shown that ME-WC increases sensitivity to both positively correlated BOLD activation and anti-correlated signals. Moreover, ME-WC and ME-ICA both increase the peak Z-scores and the number of statistically significant voxels, consistent with voxel-wise whole-brain analysis (Table 2 and Fig. 3b).

Additionally, we investigated the effect of the number of echoes on the main effect, considering whether later echoes contribute significantly. Using the same analysis pipeline, we compared whole-brain analysis results from data including the first 3, 4, and all 5 echoes (Fig. S4). As the number of echoes increased, the size of the clusters also increased.

The current analysis pipeline utilized a specific method of combining multiple echo image series into a single image series: T₂^*-weighted combination, also known as “optimal combination”, is used in most ME-EPI studies as it is designed to maximize BOLD-contrast sensitivity (Baracchini et al., 2021; Gonzalez-Castillo et al., 2016; Kundu et al., 2017; Lynch et al., 2020; Roth et al., 2020). However, an interesting finding emerged when we processed the ME-EPI data separately and treated each individual echo series as a 1E-EPI dataset. Surprisingly, no supra-threshold activation was detected at the first echo (p < 0.05, uncorrected), while the majority of activations were observed at the later echoes (p < 0.05, uncorrected) (Fig. S5). When examining the Z values for the ROI analyses across all 5 individual echoes, we also observe a similar trend (Table S2). However, regions with short T₂^* (such as the lateral and medial OFC and entorhinal cortex in this study) predominantly assigned higher weights to the first two echoes.

Interestingly, in the process of understanding the T₂^*-weighted combination method implemented in the TEDANA package, we observed that directly estimating T₂^* by fitting the Eq. (1) tends to overestimate T₂^* in regions with high susceptibility (Table S3). As described in the Method, we proposed using Eq. (3) as the objective function to estimate T₂^* for each voxel. While this approach improved some estimations, work is needed to improve the accuracy of these estimations in the future.

There have also been recent advances in denoising methods that leverage phase information, such as Noise Reduction with Distribution Corrected (NORDIC) PCA (Dowdle et al., 2023; Moeller et al., 2021; Vizioli et al., 2021). These methods have demonstrated the ability to further enhance BOLD sensitivity in resting-state fMRI. However, our study did not acquire phase data, which prevented us from testing these advanced denoising methods in our data.

In this study, we utilized a limited set of two odor stimuli, along with a control stimulus. Although this approach was adequate to compare the ME-EPI and 1E-EPI acquisition methods, future studies will aim to incorporate a wider range of odor stimuli. By including multiple odor stimuli, we can expand our analysis to encompass multivariate analysis, such as Multivariate Pattern Analysis (MVPA), which has been commonly utilized in other olfaction task-based fMRI studies (Donoshita et al., 2021; Howard et al., 2009; Howard & Gottfried, 2014; Perszyk et al., 2023), though we expect also to see an improvement in multivariate analysis using ME-ICA.

In conclusion, our study demonstrated the utility of an ME-EPI protocol combined with a pre-configured analysis pipeline for fMRI in olfactory cortex. This approach enhanced sensitivity to task-dependent BOLD contrast changes in olfactory regions, and post-hoc power analyses indicate that group effects can be detected with a reduced sample size compared to standard 1E-EPI. Although our protocol and pipeline were specifically designed for task-based fMRI in the presence of large static susceptibility, a longstanding challenge in the field of fMRI, they showcase the superior denoising capabilities of ME-EPI for task-based fMRI more generally. These findings strongly support the adoption of ME-EPI over 1E-EPI for future task-based fMRI studies.

The unthresholded Z-maps showing the whole-brain results for univariate analyses were made available on Neurovault (https://neurovault.org/collections/KGZMOMWR/). The modified TEDANA codes with updated T₂* estimation methods described in Section 2.7 are available on GitHub (https://github.com/luxwig/tedana/tree/altT2s).

L.S.Z.: Conceptualization, Methodology, Software, Formal Analysis, Investigation, Data Curation, Visualization, Writing – Original Draft, and Project administration. C.U.R.: Conceptualization, Methodology, Investigation, and Writing – Review & Editing. M.D.T.: Conceptualization, Methodology, Writing – Review & Editing, Supervision, and Resources. J.A.D.: Conceptualization, Methodology, Writing – Review & Editing, Supervision, and Resources. J.A.G.: Conceptualization, Methodology, Writing – Review & Editing, Supervision, Resources, and Funding acquisition.

This work was supported by the National Institutes of Health award R01DC019405 [awarded to J.A.G.].

The authors have declared that no competing interests exist.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00362.