# Non-Cartesian non-Fourier fmri imaging for high-resolution retinotopic mapping at 7 Tesla

Zaineb Amor, Pierre-Antoine Comby, Caroline Le Ster, Alexandre Vignaud, Philippe Ciuciu

► To cite this version:

Zaineb Amor, Pierre-Antoine Comby, Caroline Le Ster, Alexandre Vignaud, Philippe Ciuciu. Non-Cartesian non-Fourier fmri imaging for high-resolution retinotopic mapping at 7 Tesla. CAMSAP 2023 - IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, IEEE, Dec 2023, Los Suenos, Costa Rica. hal-04363378

HAL Id: hal-04363378

https://hal.science/hal-04363378v1

Submitted on 24 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# NON-CARTESIAN NON-FOURIER FMRI IMAGING FOR HIGH-RESOLUTION RETINOTOPIC MAPPING AT 7 TESLA

Zaineb Amor¹, Pierre-Antoine Comby¹,², Caroline Le Ster¹, Alexandre Vignaud¹, Philippe Ciuciu¹,²

¹CEA, Joliot, NeuroSpin, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

²Inria, MIND, Université Paris-Saclay, F-91120 Palaiseau, France

## ABSTRACT

Functional MRI (fMRI) is an invaluable tool for neuroscience, but it necessitates a trade-off between spatial and temporal resolution to maintain a reasonable temporal signal-to-noise ratio (tSNR). Non-Cartesian acquisition schemes are more efficient sampling strategies as compared to Cartesian ones and were proposed as a way to push the limits of spatiotemporal resolution further. 3D-SPARKLING is a novel non-Cartesian scheme recently evaluated for fMRI applications. However, 3D-SPARKLING is highly sensitive to $B_0$ field imperfections. Such imperfections are detrimental to fMRI applications, especially at ultra-high magnetic fields. In this work, we collect measurements of the static and dynamic $B_0$ field perturbations concurrently with 3D-SPARKLING fMRI data acquisition and retrospectively correct these perturbations during image reconstruction. The advantages of this image reconstruction strategy are assessed on the statistical sensitivity to the BOLD contrast during a retinotopic mapping fMRI experiment. Importantly, a gain of 44% (resp., 159%) additionally activated voxels was quantified when adopting this brain activity-enhanced image reconstruction technique at a type-I statistical control level of 0.001 without multiple comparisons correction (resp., 0.05 with false discovery rate correction). Additionally, significantly improved retinotopic maps were retrieved on the cortical surface.

**Index Terms—** fMRI, non-Cartesian, Compressed Sensing, ultra-high field MRI, non-Fourier imaging, 7T.

## 1. INTRODUCTION

Functional MRI is today one of the most widespread neuroimaging techniques thanks to its non-invasiveness and high spatial resolution compared to other modalities such as electroencephalography (EEG) or magnetoencephalography (MEG). In contrast, fMRI data has a low temporal resolution on average (around 1s) as the data acquisition, performed in k-space, is relatively slow at least in the Cartesian imaging setting, where the Nyquist criterion must be met. Accelerated Cartesian acquisition schemes such as 2D-SMS EPI [1] or 3D-EPI [2], which sample the k-space using readouts that span parallel planes, combined with parallel imaging [3,4], enable undersampling without aliasing artifacts in the image domain by exploiting the data redundancy from the multiple receiver channels during volume-wise inverse Fourier transform reconstruction [5, 6]. Non-Cartesian fMRI methods such as spiral readouts [7] or TURBINE [8] have been popularized in recent years due to their ability to sample the k-space more efficiently. 3D-SPARKLING [9] is a new encoding scheme that was assessed for fMRI and proved to be competitive with state-of-the-art 3D-EPI to detect evoked brain activity during a visual task [10]. A comparison with 2D-SMS EPI would be more complicated as 2D and 3D acquisitions generate different contrasts.

As an accelerated and massively under-sampled non-Cartesian acquisition method, 3D-SPARKLING requires a nonlinear Compressed Sensing (CS) image reconstruction method involving a non-uniform Fourier operator and sparsity-promoting priors. Proximal gradient methods such as the Fast Iterative Shrinkage-Tresholding Algorithm (FISTA) [11] or the Promixal Optimized Gradient Method (POGM) [12] can be used to solve the associated optimization problem. Furthermore, unlike Cartesian methods where $\Delta B_0$ inhomogeneities result mainly in geometric distortions and can be corrected by means of preprocessing steps such as the TOPUP approach [13], the impact of $B_0$ imperfections (static inhomogeneities in space and temporal fluctuations over time) is more complex in 3D-SPARKLING acquisitions and result in a superposition of blurring, signal loss and geometric distortions. To correct for both static and dynamic imperfections, the model of the collected signal in k-space during each sampling trajectory must be extended from the classic (non-uniform) Fourier transform to a pseudo-Fourier operator.

In this work, we retrospectively perform static and up-to-the-first-order dynamic $B_0$ perturbations correction on resting-state and task-based fMRI data collected during a retinotopic mapping experiment. This correction is achieved by means first of collecting additional external measurements of the field fluctuations using a field Camera [14, 15](Skope Clip-on camera) and second through the definition of a non-Fourier forward operator that encodes an extended signal model taking $B_0$ field imperfections [16] into account. A study has been conducted over three different healthy volunteers, and the impact of $B_0$ imperfections corrections is assessed for image quality, tSNR, and sensitivity to the BOLD contrast.

## 2. THEORY

### 2.1. Extended signal model

3D-SPARKLING is a segmented k-space encoding scheme, meaning that the k-space data is collected over $N_{shot}$ consecutive readouts (also called shots) separated by different excitation RF pulses. All readouts have the same duration denoted $T_{obs}$, and each one of them samples the center of the k-space at the echo time denoted $T_E$. Furthermore, in the case of parallel imaging (multiple receive coils), all coils collect simultaneously their own data in the k-space domain convolved with their own sensitivity profile. During the s-th readout, the NMR signal $\mu_s^\ell(t)$ collected by the $\ell$-th coil for each $t \in [TE - \frac{T_{obs}}{2}, TE + \frac{T_{obs}}{2}]$ can be modeled by the idealized signal model in Eq. (1):

$$ \mu_{\ell}^{s}(t) = \int_{\text{FOV}} \bar{x}_{\ell}(\mathbf{r})e^{-2i\pi k_{s}(t) \cdot \mathbf{r}} \, d\mathbf{r}, \quad (1) $$

where $\bar{x}_{\ell}(\mathbf{r}) = \zeta_{\ell}(\mathbf{r})x(\mathbf{r})$, $x(\mathbf{r})$ and $\zeta_{\ell}(\mathbf{r})$ respectively encode the image and sensitivity map associated with the $\ell$-th coil at the spatial
position $\mathbf{r}$. Note that $\mathbf{k}_s$ ($\text{in m}^{-1}$) is the prescribed (theoretical) k-space position at time $t$ (in s) relative to the s-th readout. Therefore, a non-uniform Fourier transform binds $\mu_{\ell}^s(t)$ to $\bar{x}_{\ell}(\mathbf{r})$.

Eq. (1) supposes a homogeneous $B_0$ and accurate application of the prescribed trajectories ($\tilde{k}=[\mathbf{k}_s]_{1\le s\le N_{\text{shot}}}$) by the MR system. Unfortunately, in the case of highly under-sampled 3D-SPARKLING, strong degradation in the image quality occurs due to $B_0$ imperfections arising from the patient (susceptibility changes at the tissue-air interfaces, physiological noise) and the system (eddy currents, temperature drifts, system instabilities) [17–19]. One way to overcome this issue is to take such imperfections in the signal model into account during MR image reconstruction: Eq. (1) then becomes

$$ \mu_{\ell}^{s}(t) = e^{-2i\pi t \Delta B_{0,s}^{\text{dyn}}} \int_{\text{FOV}} \bar{x}_{\ell}(\mathbf{r}) e^{-2i\pi[\Delta B_{0,s}^{\text{stat}}(\mathbf{r})t + \tilde{k}_{s}(t) \cdot \mathbf{r}]} \, d\mathbf{r} \quad (2) $$

where $\Delta B_{0,s}^{\text{stat}}(\mathbf{r})$ (in Hz) denotes the static inhomogeneities of the $B_0$ field in space and $\Delta B_{0,s}^{\text{dyn}}$ (in Hz) and $\tilde{k}_s = \mathbf{k}_s + \delta\mathbf{k}_s$ (in m$^{-1}$) denote respectively its zeroth order dynamic fluctuation and the measured trajectory (deviated from the prescribed one due to first order fluctuation $\delta\mathbf{k}_s$). $\Delta B_{0,s}^{\text{dyn}}$ is slowly varying and considered constant during a shot. In Eq. (2), the term $\Delta B_{0,s}^{\text{stat}}(\mathbf{r})t$ depends on the image domain making the integral dependent both on the image and k-space domains which is not compatible with the usual Fourier transform model.

## 2.2. Linear approximation of the non-Fourier model

According to Eq. (2), the discretized adjoint operator can be written as follows:

$$ \bar{x}_{\ell}(\mathbf{r}_n) = \sum_t e^{2i\pi t \Delta B_{0,s}^{\text{dyn}}} \mu_{\ell}^{s}(t) e^{2i\pi[\Delta B_{0,s}^{\text{stat}}(\mathbf{r})t + \tilde{k}_{s}(t) \cdot \mathbf{r}]} \quad (3) $$

The mixed term $e^{2i\pi \Delta B_{0,s}^{\text{stat}}(\mathbf{r}_m)t} = \sum_P b_{m,p} c_{p,n}$ is split in a P-rank linear decomposition using a SVD (as implemented in [20] following the methods in [16, 21]) and we then recover a linear combination of P (non-uniform) Fourier transforms as follows:

$$ \bar{x}_{\ell}(\mathbf{r}_n) = \sum_{p=1}^{P} c_{p,n} \sum_{t_m=0}^{T_{\text{obs}}} b_{m,p} e^{\underbrace{\frac{2i\pi t \Delta B_{0,s}^{\text{dyn}}}{\tilde{\mu}_{\ell}^{s}(t_m)}}_{\tilde{\mu}_{\ell}^{s}(t_m)}} e^{\frac{2i\pi \tilde{k}_s(t_m) \cdot \mathbf{r}_n}{\tilde{\mu}_{\ell}^{s}(t_m)}} \quad (4) $$

The higher the number of interpolators P in the sum, the more accurate the approximation, however, at the expense of additional computing time.

Since the term related to $\Delta B_{0,s}^{\text{dyn}}$ is outside of the integral in Eq. (2), the zeroth order dynamic fluctuations can be corrected by simply demodulating each $\mu_{\ell}^s(t)$ by the corresponding $e^{2i\pi t \Delta B_{0,s}^{\text{dyn}}}$ to obtain $\tilde{\mu}_{\ell}^s(t)$. As Eq. (4) holds for all frequencies $\tilde{k}_s(t_m)$ and locations $\mathbf{r}_n$ across all the $N_{\text{shot}}$ readouts, we can summarize the perturbed acquisition in Eq. (4), as a linear combination of adjoint non-uniform Fourier transforms $\tilde{\mathcal{F}}_{\tilde{\Omega}}$, yielding a coil-specific image $\bar{x}_{\ell}$ from the measured frequencies at locations $\tilde{\Omega}$ and associated corrected values ($\tilde{\mu}_{\ell}$):

$$ \bar{x}_{\ell} = \sum_{p=1}^{P} c_p \odot \mathcal{F}_{\tilde{\Omega}}^{*}(b_p \odot \tilde{\mu}_{\ell}) = \tilde{\mathcal{F}}_{P,\tilde{\Omega}}^{*}(\tilde{\mu}_{\ell}), \quad (5) $$

where $\odot$ denotes the element-wise product.

## 2.3. Sparse reconstruction

In this work, we perform a CS-based reconstruction, where each volume is reconstructed independently from the others using L-channel k-space data. Therefore, for each complex-valued volume ($x_j \in \mathbb{C}^N, j = 1, \dots, N_t$), we consider the complex-valued k-space data $\tilde{y}_{j,\ell} = [\tilde{\mu}_{\ell}^1, \dots, \tilde{\mu}_{\ell}^{N_{\text{shot}}}] + n_{j,\ell}$, associated with the set of sampled frequencies points $\tilde{\Omega}_j = \tilde{\Omega} = [\tilde{k}_1, \dots, \tilde{k}_{N_{\text{shot}}}]$ relative to the j-th volume. Future work will address time or scan-varying under-sampling patterns $\tilde{\Omega}_j$. Here $n_{j,\ell}$ is an additive zero-mean white Gaussian noise with variance $\sigma_{\ell}^2$ and $\tilde{y}_{j,\ell} \in \mathbb{C}^M, N$ and $M$ are, respectively, the number of voxels in each volume and that of the k-space measurements collected by each coil. Then the reconstruction problem consists in minimizing the following objective function

$$ \hat{x}_j = \arg \min_{\mathbf{x} \in \mathbb{C}^N} \frac{1}{2} \sum_{\ell=1}^{L} \sigma_{\ell}^{-2} \left\| \tilde{\mathcal{F}}_{P,\tilde{\Omega}} S_{\ell} \mathbf{x}_j - \tilde{y}_{j,\ell} \right\|_2^2 + \lambda g(\Psi(\mathbf{x}_j)) \quad (6) $$

where $S_{\ell}$ denotes the sensitivity profile of the $\ell$-th coil and $\tilde{\mathcal{F}}_{P,\tilde{\Omega}}$ the linearized over P interpolators off-resonance corrected Fourier operator. $\Psi$ is the sparsifying basis (e.g. a tight frame), $g$ is typically chosen as the $l_1$-norm, and $\lambda$ is a regularization parameter. Although more advanced low-rank + sparse regularization could be envisaged easily (see [22] for instance), this requires correcting the number of degrees of freedom during subsequent statistical analysis. Due to this, every volume $\hat{x}_j$ for $j = 1, \dots, N_t$, is reconstructed independently by solving Eq. (6). Practically, in our experiments, the sparsifying basis ($\Psi$) we used was the sym8 wavelet basis and the other parameters were set to $L = 32, P = 30, \lambda = 10^{-8}$.

# 3. MATERIALS AND METHODS

## 3.1. Multi-modal data acquisition

Functional MRI data was collected at 7T (7T Magnetom investigational MRI, Siemens Healthineers, Erlangen, Germany) using a 1Tx-32Rx head coil (Nova Medical, Willmington, CO, USA) from 3 healthy volunteers during resting-state and a retinotopic mapping experiment<fnref n="1" /> at a spatial resolution of 1mm<fnref n="3" />, a temporal resolution ($TR_{vol}$) of 2.4s, $TR_{shot}/TE = 50/20$ms and a 3D field-of-view (FOV) of (192,192,128) mm<fnref n="3" />. Concurrently and for each acquired MR volume, 16 NMR probes (a field camera) were used to monitor and record the zeroth order field fluctuations $\Delta B_{0,s}^{\text{dyn}} = [\Delta B_{0,1}^{\text{dyn}}, \dots, \Delta B_{0,N_{\text{shot}}}]$ and measure the trajectories played by the MR system $\tilde{k} = [\tilde{k}_1, \dots, \tilde{k}_{N_{\text{shot}}}]$. Additionally, an external $\Delta B_{0,s}^{\text{stat}}$ map as well as external sensitivity maps were acquired using a three-echo 3D gradient recalled echo (3D GRE) sequence.

## 3.2. Functional MRI data reconstruction and preprocessing

The fMRI volumes were reconstructed by solving the minimization problem in Eq. (6) using 15 iterations of the POGM algorithm for each volume. This method is implemented in the pysap-mri [23] plugin<fnref n="2" /> of the pySAP package [24]. Motion correction and co-registration of the functional and the $T_1$-weighted anatomical MR images were applied using SPM12³. The

<footnote n="1">¹https://github.com/hbp-brain-charting/public-protocols</footnote>

<footnote n="2">²https://github.com/CEA_COSMIC/pysap-mri</footnote>

<footnote n="3">³https://www.fil.ion.ucl.ac.uk/spm/doc/biblio/</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Comparison of the (A) mean images and (B) tSNR maps yielded by the resting-state fMRI scans sequence collected in the three volunteers (one raw per subject) reconstructed with (A, right col.) and without (A, left col.) correcting off-resonance artifacts due to $B_0$ inhomogeneities. Fully corrected means that static and up-to-the-first-order dynamic terms were taken into account during image reconstruction as in Eq. (6). The orange and green arrows depict examples of recovered signal and better-reconstructed details, respectively. The improved image quality and tSNR maps are systematic across the three volunteers.</figcaption></figure>

segmentation of the $T_1$-weighted MRI scan into the cortical surface meshes was performed using FreeSurfer 7.

### 3.3. Statistical analysis of the task-based fMRI data

Retinotopic mapping fMRI data was analyzed for each participant separately, using a first-level general linear model (GLM) where two paradigm-related regressors (parametric, continuous and sinusoidal) were used to capture the BOLD fluctuations elicited by the stimulus presentation. The design matrix also encompasses 6 motion regressors, a polynomial drift and a baseline. A Fisher test over the task-related regressors was used to estimate the global effect of interest after thresholding the F-statistic maps over the entire brain for two different strategies:

(i) $p < 0.001$ without correcting for multiple comparisons

(ii) $p < 0.05$ with false discovery rate (FDR) control.

The thresholded F-statistic maps were further used to delineate regions of interest (ROIs) in which the retinotopic mapping was computed. The statistical analysis was implemented using Nilearn [25].

## 4. RESULTS

### 4.1. Improved image quality and increased tSNR

Fig. 1-(A) compares, subject-wise, the mean images derived from the resting-state sequence of fMRI volumes reconstructed with and without correcting artifacts due to static and dynamic $B_0$ perturbations.

<table><caption>**Table 1.** Gain in % of median tSNR in corrected images relative to the native tSNR (uncorrected images) computed over the brain mask for the three different volunteers. The highest gains (in bold font) are retrieved when static and up-to-the-first-order field terms are jointly corrected.</caption><thead><tr><th rowspan="2">Terms corrected</th><th colspan="4">Gain in % of median tSNR</th></tr><tr><th>V#1</th><th>V#2</th><th>V#3</th><th>Average</th></tr></thead><tbody><tr><td>ΔB<sub>0</sub><sup>stat</sup></td><td>+5</td><td>+6</td><td>+4</td><td>+5</td></tr><tr><td>ΔB<sub>0</sub><sup>dyn</sup></td><td>+20</td><td>+11</td><td>+22</td><td>+18</td></tr><tr><td>ΔB<sub>0</sub><sup>dyn</sup> & δk</td><td>+26</td><td>+13</td><td>+23</td><td>+21</td></tr><tr><td>ΔB<sub>0</sub><sup>stat</sup> & ΔB<sub>0</sub><sup>dyn</sup></td><td>+28</td><td>+18</td><td><strong>+29</strong></td><td>+25</td></tr><tr><td>ΔB<sub>0</sub><sup>stat</sup> & ΔB<sub>0</sub><sup>dyn</sup> & δk</td><td><strong>+34</strong></td><td><strong>+20</strong></td><td><strong>+29</strong></td><td><strong>+28</strong></td></tr></tbody></table>

We note that the image quality is systematically improved with the proposed correction across the three volunteers. Indeed, performing $B_0$ perturbations correction improves the overall $T_2^*$ contrast. Additionally, the lost signal is recovered, and anatomical details are sharper as depicted, respectively, with orange and green arrows. Fig. 1-(B) shows the boost in tSNR associated with $B_0$ perturbations correction, notably in the anterior and posterior cortex and, along the edges of the brain suggesting that subtle head movement-induced field fluctuations related to breathing were compensated.

Tab. 1 reports the relative gain in % of median tSNR computed over the brain mask when correcting the different field terms during image reconstruction. Although the increase is systematic across volunteers (V#1-3), the relative gain for V#2 is lower than in other
volunteers when dynamic fluctuations are corrected. This is likely because this volunteer has less intense breathing therefore, the native fMRI data is less affected by signal fluctuations. Otherwise, the relative gain reaches a plateau around 30% of gain at maximum.

## 4.2. Increased sensitivity to the BOLD effect and enhanced retinotopic maps

Tab. 2 summarizes the systematic boost in the number of activated voxels and maximum z-score values when static and up-to-the-first-order dynamic $B_0$ fluctuations imperfections were corrected. The number of activated voxels extracted using the thresholding alternative (i) (resp, (ii)) is, on average, 43.3% ± 17.2% (resp, 159.3% ± 38.6%) larger. The reported figures are consistent between the first and third participants. V#2, however, reveals fewer activated voxels, especially after FDR correction. This is likely due to larger (strong) head movement amplitudes.

<table><caption>**Table 2.** Increase in the number of activated voxels detected and the maximum z-score values extracted from task-based fMRI scans with/out correcting $B_0$ imperfections. The highest values (in bold font) are obtained when static and up-to-the-first-order $B_0$ terms - Full Correction- are used. V#2 reveals the lowest statistical significance.</caption><thead><tr><th rowspan="2">Volunteer</th><th colspan="2">p &lt; 0.001</th><th colspan="2">p &lt; 0.05</th><th rowspan="2">Z-score (max)</th></tr><tr><th colspan="2">uncorrected</th><th colspan="2">FDR corrected</th></tr></thead><tbody><tr><td>#1</td><td>6456</td><td><strong>9506</strong></td><td>3253</td><td><strong>7722</strong></td><td>8.02</td></tr><tr><td>#2</td><td>4367</td><td><strong>7405</strong></td><td>1059</td><td><strong>2204</strong></td><td>6.75</td></tr><tr><td>#3</td><td>8529</td><td><strong>10 823</strong></td><td>2503</td><td><strong>7745</strong></td><td>10.83</td></tr><tr><td>Average</td><td>6450</td><td><strong>9244</strong></td><td>2271</td><td><strong>5890</strong></td><td>8.54</td></tr></tbody></table>

Additionally, Fig. 2 shows an example (V#1) of the projection of the retinotopic maps on the cortical surface. We notice clearly that the two visual hemifields project onto the contra-lateral hemispheres in the occipital cortex, a well-known mirroring feature of the primary visual cortex. However, a larger spatial extent of the retinotopic maps is associated with $B_0$ field imperfections correction. Furthermore, the zoomed-in region of interest (ROI) in the left hemisphere illustrates a sharp gain in sensitivity along the cortical surface, suggesting that the activated voxels recovered when performing $B_0$ perturbations correction are highly relevant.

## 5. DISCUSSION AND CONCLUSION

The systematic gain in tSNR on average associated with $B_0$ imperfections correction we observed is generally consistent with previous studies even though the reported figures vary [26, 27]. Such differences arise from the experimental conditions (2D vs 3D acquisition sequence, 3T vs 7T $B_0$ magnetic field, number of scans, volunteer's movement, etc.). The increase in sensitivity associated with $B_0$ field imperfections correction we observed, notably at a more conservative statistical threshold is likely due to the gain in tSNR. Furthermore, the recovery of lost signal plays a major role in retrieving activated voxels. The enhanced image quality is also expected to yield better spatial specificity as blurring is extremely reduced.

In this work, we consider that the static and dynamic field terms evolve independently from each other since such an approximation is easy to implement and remains accurate enough. Nevertheless, it does not reflect the physics of the experiment faithfully. In fact, a truthful model would consider static and dynamic $B_0$ imperfections

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** The projection of the retinotopic maps collected in V#1 and associated with data reconstructed without and with static and up-to-the-first order field terms correction on the cortical surface. These maps were computed for an activation mask corresponding to a p-value of 0.001 and without correcting for multiple comparisons correction (i.e. alternative (i)).</figcaption></figure>

as evolving jointly: A $\Delta B_0$ map could be estimated for each volume in this case analogously to what is proposed in [28,29]. Such a solution should also account for strong head movements to be complete. Furthermore, we performed typical sequential volume-wise reconstructions where the consecutive blocks of k-space data are considered independent from one another. Despite being easy to implement and efficient, such a strategy neglects the temporal structure of the fMRI scans. Global strategies that consider the entire scan data to leverage the temporal structure of the fMRI signal during reconstruction have also been proposed and can be further combined with $B_0$ perturbations correction: Typically, low/fixed rank plus sparse or subspace-constrained reconstructions [22] can be carried out in the same framework, which could also improve the current temporal resolution.

In this work, we demonstrated a systematic and significant benefit in image quality, tSNR, and statistical sensitivity of brain activity detection and localization when correcting static and dynamic $B_0$ inhomogeneities. We can further improve the image reconstruction process by considering a global 4D $\Delta B_0$ map or more advanced priors and even soon deep learning self-supervised strategies, which are still not available for 3D fMRI data at this time.

## 6. COMPLIANCE WITH ETHICAL STANDARDS

The experimental protocol was approved by the national ethics committee under the protocol identifier CPP 100048 issued by the national Comité de Protection des Personnes (CPP Sud Méditerranée 4 number 180913, IDRCB:2018-A011761-53). All participants gave their written informed consent.
7. REFERENCES

[1] S Moeller, E Yacoub, C. A Olman, E Auerbach, and J. H. N Strupp, "Multiband multislice GE-EPI at 7 Tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI," *Magn. Reson. Med.*, vol. 63, no. 5, pp. 1144-1153, 2010.

[2] B. A Poser, P. J Koopmans, T Witzel, L. L Wald, and M Barth, "Three dimensional echo-planar imaging at 7 Tesla," *NeuroImage*, vol. 51, no. 1, pp. 261-6, 2010.

[3] Klaas P Pruessmann, "Encoding and Reconstruction in Parallel MRI," *NMR Biomed.*, vol. 19, no. 3, pp. 288-99, 2006.

[4] A. Deshmane, V. Gulani, M. A. Griswold, and N. Seiberlich, "Parallel MR Imaging," *J. of Magn. Reson. Imag.*, vol. 36, no. 1, pp. 55-72, 2012.

[5] K P Pruessmann, M Weiger, M B Scheidegger, and P Boesiger, "SENSE: Sensitivity encoding for fast MRI," *Magn. Reson. Med.*, vol. 42, pp. 952-962, 1999.

[6] M A Griswold, P M Jakob, R M Heidemann, M Nittka, V Jelus, J Wang, B Kiefer, and A Haase, "Generalized Autocalibrating Partially Parallel Acquisitions (GRAPPA)," *Magn. Reson. Med.*, vol. 47, pp. 1202-1210, 2002.

[7] L Kasper, M Engel, J Heinzle, M Mueller-Schrader, N N Graedel, J Reber, T Schmid, C Barmet, B J Wilm, K E Stephan, and K P Pruessmann, "Advances in spiral fMRI: A high-resolution study with single-shot acquisition," *NeuroImage*, vol. 246, pp. 1187-138, 2022.

[8] N N Graedel, K L Miller, and M Chiew, "Ultrahigh Resolution fMRI at 7T Using Radial-Cartesian TURBINE Sampling," *Magn. Reson. Med.*, vol. 88, no. 5, pp. 2058-2073, 2022.

[9] G R Chaithya, P Weiss, A Massire, A Vignaud, and P Ciuciu, "Optimizing full 3D SPARKLING trajectories for high-resolution Magnetic Resonance imaging," *IEEE Trans. on Med. Imag.*, vol. 41, no. 8, pp. 2105-2117, Aug. 2022.

[10] Z Amor, G R Chaithya, G Daval-Frérot, B Thirion, F Mauconduit, P Ciuciu, and A Vignaud, "Prospects of non-cartesian 3d-sparkling encoding for functional mri: A preliminary case study for retinotopic mapping," *Proceedings of the 31st ISMRM*, 2022.

[11] Amir Beck and Marc Teboulle, "A fast iterative shrinkage-thresholding algorithm for linear inverse problems," *SIAM journal on imaging sciences*, vol. 2, no. 1, pp. 183-202, 2009.

[12] D. Kim and J. A. Fessler, "Adaptive restart of the optimized gradient method for convex optimization," *J Optim Theory Appl.*, vol. 178, pp. 240-263, 2018.

[13] S M Smith, M Jenkinson, M W Woolrich, C F Beckmann, T E J Behrens, H Johansen-Berg, P R Bannister, M De Luca, I Drobnjak, D E Flitney, R K Niazy, J Saunders, J Vickers, Y Zhang, N De Stefano, J. M Brady, and P M Matthews, "Advances in functional and structural MR image analysis and implementation as FSL," *NeuroImage*, vol. 23, pp. S208-S219, 2004.

[14] N. De Zanche, C. Barmet, J. A. Nordmeyer-Massner, and K. P. Pruessmann, "NMR Probes for Measuring Magnetic Fields and Field Dynamics in MR Systems," *Magn. Reson. Med.*, vol. 60, pp. 176-186, 2008.

[15] C. Barmet, N. De Zanche, B. J. Wilm, and K. P. Pruessmann, "A transmit/receive system for magnetic fieldmonitoring of in vivo mri," *Magn. Reson. Med.*, 2009.

[16] B P Sutton, D C Noll, and J A Fessler, "Fast, iterative image reconstruction for mri in the presence of field inhomogeneities," *IEEE Trans. on Med. Imag.*, vol. 22, no. 2, pp. 178-188, 2003.

[17] S J Vannesjo, M Haeberlin, L Kasper, M Pavan, B J Wilm, C Barmet, and K P Pruessmann, "Gradient system characterization by impulse response measurements with a dynamic field camera," *Magn. Reson. Med.*, vol. 69, pp. 583-593, 2013.

[18] L. Kasper, S. Bollmann, S. J. Vannesjo, S. Gross, M. Haeberlin, and K. P. Pruessmann B. E. Dietrich, "Monitoring, analysis, and correction of magnetic field fluctuations in Echo Planar Imaging time series," *Magn. Reson. Med.*, vol. 74, pp. 396-409, 2015.

[19] S. Bollmann and L. Kasper, S. J. Vannesjo, A. O. Diaconescu, B. E. Dietrich, S. Gross, K. E. Stephan, and K. P. Pruessmann, "Analysis and correction of field fluctuations in fMRI data using field monitoring," *NeuroImage*, vol. 154, pp. 92-105, 2017.

[20] G. Daval-Frérot, A/ Massire, B. Mailhe, M. Nadar, A. Vignaud, and P. Ciuciu, "Iterative static field map estimation for off-resonance correction in non-cartesian susceptibility weighted imaging," *Magn. Reson. Med.*, vol. 88, no. 4, pp. 1592-1607, 2022.

[21] J A Fessler, S Lee, V T Olafsson, H R Shi, and D C Noll, "Toeplitz-based iterative image reconstruction for MRI with correction for magnetic field inhomogeneity," *IEEE Trans. on Sig. Process.*, vol. 53, no. 9, pp. 3393-3402, 2005.

[22] Andrii Y Petrov, Michael Herbst, and V Andrew Stenger, "Improving temporal resolution in fMRI using a 3D spiral acquisition and low rank plus sparse (L+S) reconstruction," *NeuroImage*, vol. 157, pp. 660-674, 2017.

[23] L. El Gueddari, C. Giliyar Radhakrishna, Z. Ramzi, S. Farrens, S. Starck, A. Grigis, J.-L. Starck, and P. Ciuciu, "PySAP-MRI: A python package for mr image reconstruction," in *ISMRM workshop on Data Sampling and Image Reconstruction*, 2020.

[24] S. Farrens, A. Grigis, L. El Gueddari, Z. Ramzi, G.R. Chaithya, S. Starck, B. Sarthou, H. Cherkaoui, P. Ciuciu, and J.-L. Starck, "PySAP: Python Sparse Data Analysis Package for multidisciplinary image processing," *Astronomy and Computing*, vol. 32, 2020.

[25] A. Abraham, F. Pedregosa, M. Eickenberg, P. Gervais, A. Mueller, J. Kossaifi, A. Gramfort, B. Thirion, and G. Varoquaux, "Machine learning for neuroimaging with scikit-learn," *Frontiers in Neuroinformatics*, vol. 8, 2014.

[26] S. J. Vannesjo, B. J. Wilm, Y. Duerst, S. Gross, D. O. Brunner, B. E. Dietrich, T. Schmid, C. Barmet, and K. P. Pruessmann, "Retrospective correction of physiological field fluctuations in high-field brain MRI using concurrent field monitoring," *Magn. Reson. Med.*, vol. 73, pp. 1833-1843, 2015.

[27] J. M. Schwarz, R. Stirnberg, P. Ehses, and T. Stocker, "Correction of physiological field fluctuations in high- and low-resolution 3D-EPI acquisitions at 7 Tesla," *Proceedings of the 27th ISMRM*, 2019.

[28] B Dymerska, B A Poser, M Barth, S Trattnig, and S D Robinson, "A method for the dynamic correction of B0-related distortions in single-echo EPI at 7 T," *NeuroImage*, vol. 168, pp. 312-331, 2018.

[29] M W Haskell, A Lahiri, J-F Nielsen, J A Fessler, and D C Noll, "FieldMapNet MRI: Learning-based mapping from single echo time BOLD fMRI data to fieldmaps with model-based reconstruction," *Proceedings of the 31st ISMRM*, 2022.