Impact of B0$$ {\mathrm{B}}_0 $$ field imperfections correction on BOLD sensitivity in 3D‐SPARKLING fMRI data

Static and dynamic B0$$ {\mathrm{B}}_0 $$ field imperfections are detrimental to functional MRI (fMRI) applications, especially at ultra‐high magnetic fields (UHF). In this work, a field camera is used to assess the benefits of retrospectively correcting B0$$ {\mathrm{B}}_0 $$ field perturbations on Blood Oxygen Level Dependent (BOLD) sensitivity in non‐Cartesian three‐dimensional (3D)‐SPARKLING fMRI acquisitions.


INTRODUCTION
The race toward higher resolution in functional MRI (fMRI) has motivated a significant amount of technological development in both hardware (higher fields and better coil architectures) and software (acquisition, reconstruction, and preprocessing tools). 1 In this collective endeavor to reach subsecond and/or submillimeter resolution, accelerated acquisition schemes play a major role: Cartesian Echo-Planar-Imaging (EPI) schemes such as two-dimensional SMS-EPI 2 and three-dimensional (3D)-EPI 3 combined with parallel imaging are today the proven norm in high-resolution fMRI studies due to their efficient use of the gradient capabilities.They are currently challenged by spiral readouts 4 which are today one of the most promising and investigated non-Cartesian sampling patterns for high-resolution fMRI.3D-SPARKLING 5 is a novel non-Cartesian acquisition method based on compressed sensing (CS) theories.It generates optimization-driven multi-shot sampling patterns that fit a target sampling density in k-space and meet the hardware constraints regarding maximum gradient amplitude and slew rate.This multishot readout performs a variable density sampling while maintaining local uniformity in the sampling, and results in arbitrary shots.
It produces low-level noise-like artifacts after compressed sensing-based MR image reconstruction.Upon comparison with 3D-EPI at 1-mm isotropic spatial resolution and 2.4s-TR, 6 3D-SPARKLING has recently been assessed as a competitive acquisition technique for high spatial resolution task-based fMRI.
Like EPI and spiral readouts, 3D-SPARKLING is sensitive to static and dynamic B 0 perturbations, which can be detrimental for fMRI applications as they degrade the sensitivity to the Blood Oxygen Level Dependent (BOLD) contrast and the temporal signal-to-noise ratio (tSNR). 7On one hand, moving to ultra-high magnetic field ensures a higher tSNR and temporal Contrast-to-Noise Ratio and, therefore, improved sensitivity to the BOLD effect. 8Additionally, at ultra-high magnetic fields, the boost in the contribution of capillaries in the BOLD contrast is more pronounced as compared to large veins, 9 which enhances spatial specificity.On the other hand, the influence of B 0 imperfections becomes more prominent at ultra-high magnetic fields (7T and beyond).If uncorrected, such imperfections can cancel out the gain in sensitivity and specificity expected by moving to higher fields: Correcting off-resonance effects due to B 0 inhomogeneities is therefore one way to secure the advantages of high fields without being subject to their disadvantages.
The static contribution of spatially varying B 0 field offsets arises mainly from the susceptibility differences at tissue-air interfaces and is, to some extent, routinely corrected in MRI in general.In brain MRI, this issue is most prominent near the sinuses, buccal cavity and ear canals.Static B 0 inhomogeneities can cause geometric distortions, blurring, shading as well as strong signal loss.They are usually minimized by means of shimming: Optimizing shimming for brain MRI, especially at high fields (3T and beyond) is an active field of research. 10,113][14] Preprocessing steps such as the TOPUP approach 15,16 can also be applied in the case of EPI fMRI to further correct geometric ΔB 0 distortions.
Dynamic perturbations of the B 0 field can stem from the system (instabilities in the gradient system, temperature drifts, eddy currents 17 ) or from the volunteer (e.g., breathing and heart rate).They are typically divided into a global zeroth order term, first-order deviations from the nominal trajectories, and some higher-order terms. 18,19ynamic field fluctuations are time-varying and constitute a potential confound for fMRI in regions where a significant BOLD effect is sought 20 : They may act as nuisance variables 21 by causing intensity changes that are unrelated to the neuro-vascular coupling and degrade the tSNR.Similar to static field inhomogeneities, dynamic field fluctuations can be compensated prospectively during the acquisition or corrected retrospectively during image reconstruction using time-varying estimates 14 : Navigators are typically implemented into sequences to monitor dynamic field fluctuations. 22However, navigators require alteration of the pulse sequences and sometimes lengthen the acquisition time.Physiological probes such as breathing belts can also be used to correct physiology-related field fluctuations during acquisition but only give an indirect measure of the induced field perturbations.
A method that uses NMR probes 23,24 to monitor the field fluctuations has been proposed. 25Unlike navigators, using NMR probes does not require alteration of the pulse sequence and allows for field monitoring in a time-locked fashion instead of having access only to average measurements of the temporal field fluctuations.The system presented in References 23 and 24, uses transmit and receive probes that contain an NMR active product excited by applying a radiofrequency pulse at its Larmor frequency.It later evolved into a field camera using doped Fluorine ( 19 F) NMR probes and was used for field monitoring with anatomical MRI and fMRI. 19,26,27Given the difference in Larmor frequency of doped Fluorine and Hydrogen, such probes enable monitoring the field concurrently with the imaging process in the case of proton MR imaging.However, the T 2 relaxation time of the probes is around 45ms at 7T, implying that a relatively long TR probe is needed to eliminate residual transverse magnetization between the shots and achieve a proper steady-state at the level of the probes.Using only radiofrequency spoiling, the field camera requires a constraining minimum TR probe of roughly 110 ms.Such a situation is particularly challenging for 3D fMRI applications as all shots are, ideally, monitored in real-time using TR probe = TR shot ≤50 ms.
In Reference 27, the authors managed to use this system to acquire realistic 3D-EPI fMRI data with a short TR shot and a long TR probe : They assumed repeatable readouts between the shots, skipped monitoring some shots and interpolated the missing NMR probe data.A similar strategy was used in Reference 4 where the same system was used to correct up-to-the-first-order dynamic field fluctuations in single-shot high-resolution spiral fMRI data: Field fluctuations were recorded for every third shot using a TR probe = 270 ms and a TR shot = 90 ms.Such strategies are impractical for 3D-SPARKLING applications given the pseudo-random nature of the sampling pattern, that is, the variability across shots and thus the impossibility to interpolate NMR probe data across shots: In Reference 28, we studied the benefit of correcting ΔB 0 perturbations on dynamic 3D-SPARKLING acquisitions without extending the study to realistic fMRI data because of the long TR probe constraint.
In this paper, we alleviate this issue using an external spoiling gradient and adapting the experimental protocol to enable the use of Skope's Clip-on field Camera with TR probe = TR shot = 50 ms challenging its long TR probe constraint and evaluate it for 1mm 3 3D-SPARKLING retinotopic mapping and resting-state fMRI acquisitions.We demonstrate the feasibility of our experimental protocol and study the impact of static and dynamic B 0 field imperfections correction during image reconstruction on fMRI data in terms of image quality, tSNR, sensitivity to the BOLD contrast and quality of the retinotopic maps.To the best of our knowledge, it is the first time such a thorough comparison has been performed using the Skope Clip-on field Camera for fMRI data at 7T and non-Fourier reconstruction.dr (1)

Extended signal model
where x  (r) =   (r)x(r) is the image x multiplied by the sensitivity profile of the th coil   at the spatial position r.Note that ΔB 0,stat (r) (in Hz) denotes the static inhomogeneities of the B 0 field in space whereas ΔB s 0,dyn (in Hz) and ks = k s + k s (in m −1 ) denote respectively its zeroth-order dynamic fluctuation and the measured trajectory deviated from the prescribed one (k s ) due to first-order fluctuation k s .ΔB s 0,dyn is slowly varying and considered constant during a shot but a more temporally resolved measure can be written as k s 0 (t) ⋍ 2tΔB s 0,dyn (in rad).In Equation (1), the term ΔB 0,stat (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 (nonuniform) Fourier transform model.

Linear approximation of the non-Fourier model
According to Equation (1), the discretized adjoint operator can be written as follows: The mixed term e 2ΔB 0,stat (r n )t m = ∑ P p=1 b m,p c p,n can be split in a rank-P linear decomposition using a SVD 13,29 to recover a sum of P (nonuniform) Fourier transform as follows: Since the term related to ΔB s 0,dyn is outside of the integral in Equation (1), the zeroth-order dynamic fluctuations can be corrected by simply demodulating each  s  (t) with the corresponding e 2tΔB s 0,dyn ⋍ e k s 0 (t) to obtain μs  (t).As Equation (3) holds for all frequencies ks and locations r n across the N shot readouts, we can summarize the perturbed acquisition in Equation ( 4), as a sum of adjoint nonuniform Fourier transforms  Ω, yielding a coil-specific image x  from the measured frequencies locations Ω and associated corrected values (μ  ): where ⊙ denotes the element-wise product.

Experimental protocol and data acquisition
The study was conducted at 7T MRI (7T Magnetom investigational device, Siemens Healthineers) on three healthy volunteers (two males) aged between 20 and 40 years with normal-to-corrected vision using a 1Tx-32Rx head coil (Nova Medical).The experimental protocol was approved by the national ethics committee (Comité de Protection des Personnes) under the protocol identifier CPP 100048 (CPP Sud Méditerranée 4 number 180913, IDRCB:2018-A011761-53).All participants gave their written informed consent.
Task-based fMRI was performed using two consecutive runs (one clockwise and one counter-clockwise) of a classical retinotopic mapping paradigm using the stimulation code available online * .Retinotopic mapping and resting-state data with normal breathing (NB) was collected from the three participants.Additional data with forced breathing (FB) and while performing a hand-to-chin movement (HC) was collected at resting-state from one volunteer.The functional data was collected with T * 2 -weighted 3D-SPARKLING acquisitions at a spatial resolution of 1 mm 3 , a temporal resolution of 2.4 s, TR shot /TE = 50/20 ms and a 3D field-of-view (FOV) of (192, 192, 128) mm 3 .A three-echo 3D gradient recalled echo (3D GRE) sequence was used to obtain both a ΔB 0 (ΔB 0,stat ) and external sensitivity maps.Table 1 details the parameters of the sequences.
Concurrently and for each acquired fMRI volume, 16 NMR probes from the field Camera, Cranberries edition (Skope Magnetic Resonance Technologies AG) were used to monitor and record the zeroth order field fluctuations over the acquisition window ] and measure the trajectories played by the MR system k = [ k1 , … , kN shot ] for a TR probe = TR shot = 50 ms.The residual magnetization resulting from the use of such a short TR probe was destroyed using the strongest spoiling gradient (470 mT*ms/m) implementable within a 50 ms-TR shot 3D-SPARKLING sequence.

Image reconstruction and preprocessing
The fMRI volumes were reconstructed independently from each other following the extended recorded signal model in Equation (1) and by solving the minimization problem in Equation (S1) in Section S1 using the Proximal Optimized Gradient Method (POGM) algorithm.This method is implemented in the pysap-mri 30 plugin † of the pySAP package. 31otion correction and co-registration of the functional and the T 1 -weighted anatomical data were applied using SPM12 ‡ .Except for a very specific step when estimating BOLD phase maps, no spatial smoothing was applied to maintain the native spatial resolution.The segmentation of the cortical surface into a pial, white, inflated, and sulcus meshes was performed using the T 1 -weighted anatomical scan and FreeSurfer 7.

Statistical analysis
Retinotopic mapping data was analyzed for each participant separately, using a two-session first-level general linear model (GLM) and the following block-diagonal design-matrix X: where N vol and Q∕2 are respectively, the number of volumes and the number of regressors within a single retinotopic session.The non-zero diagonal blocks X 1 and Specifications of the different sequences used in each participant: 3D-SPARKLING was used to acquire the functional data at 1 mm 3 , 2.4s-TR v resolution and with a readout duration of 26.88ms.

TE (ms)
TR Notes: 120 volumes were collected for each retinotopic run (either clockwise or counter-clockwise), whereas one resting-state run consisted of 125 volumes.The GRE sequence used had three echoes: Raw data of the full dataset (3 echoes) and raw data from only the first echo were used to estimate the ΔB 0 and sensitivity maps, respectively.The estimated maps were interpolated to match the 3D FOV and resolution of the fMRI scans.The XFL sequence was used for B + 1 mapping and calibration.The MP2RAGE sequence was used to acquire an anatomical T 1 w scan.
X 2 are, respectively, associated with the experimental paradigm that is carried out during the clockwise and counter-clockwise sessions.Each block X s is defined as follows: where two paradigm-related parametric, continuous and sinusoidal regressors (x s,1 and x s,2 ) serve to capture the BOLD fluctuations elicited by the stimulus presentation.The motion regressors are denoted x mot s,1 to x mot s,6 and x pol s,1 and x bas s,1 model a polynomial drift and the baseline, respectively.
First, a Fisher test over the task-related regressors (x task 1,1 , x task 1,2 , x task 2,1 , and x task 2,2 ) 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. 32e thresholded F-statistic maps were further used to create activation masks over the regions where the activations are statistically significant.
Second, a Student's t-test was performed over each task-related regressor separately and the corresponding z-scores were used to estimate the BOLD phase maps  Clock and  CClock as the voxel-wise arctangent of the ratio of the corresponding task-regressors (x task 1,1 and x task 1,2 for the clockwise session and x task 2,1 and x task 2,2 for the counter-clockwise session).
Then, after compensating for the recorded BOLD response delay ( ) due to the haemodynamic delay in  Clock and  CClock , we derived the overall retinotopic phase estimate as follows: The implementation was done through the Nilearn § package.

Evaluation
The mean images computed from the resting-state time series were used to inspect the impact of B 0 perturbations correction on image quality qualitatively.In order to gain a deeper insight into the impact of each field term, six distinct reconstruction strategies accounting for different combinations of field terms were used (cf, Table S2): The resting-state fMRI scans were reconstructed using strategies (a)-(f) was and further used to compute the tSNR maps in order to rank the impact of each contribution both visually on the image quality and quantitatively.
In order to gain insight on the differences between the three volunteers and several physiological noise scenarios, the power spectra of the recorded ΔB 0,dyn measurements were computed over the whole resting-state fMRI acquisition at an individual level and visualized for specific frequency intervals.
To go one step further and investigate the influence of B 0 imperfections correction on the statistical performances when detecting evoked brain activity, the retinotopic mapping fMRI data was used to estimate activation maps by computing z-score maps from the global effects of interest at the subject level.The impact on the sensitivity to the BOLD effect as well as the prevalence of true versus false positives was assessed according to the following qualitative (q) and quantitative (Q) criteria: (1) qActiv1: The statistically significant z-score maps obtained from the fMRI volumes reconstructed with strategies (a) and (f) and thresholded according to strategies (i) and (ii) as explained in Section 3.3 were compared one another subjectwise.
(2) QActiv1: The Number of activated voxels (defined using thresholding strategies (i) and (ii)) and maximum z-score values were compared at the subject level once the fMRI volumes were reconstructed using strategies (a) and (f).
Finally, the impact on the quality of the BOLD phase maps was evaluated according to the following criterion: (3) qPhase1: The volumetric BOLD phase map derived from the fMRI scans reconstructed using strategies (a) and (f) were compared subjectwise.Furthermore, the BOLD phase maps corresponding to V#1 were projected onto the cortical surface for both hemispheres and visually assessed.

Achieving steady state signal of the NMR probes
Figure 1 depicts S raw (t = 0) (the raw signal measured by the field camera before any estimation of field fluctuations) for the first 96 consecutive shots of the multirepetition 3D-SPARKLING sequence used to acquire in vivo fMRI data without/with the application of an external spoiling gradient through the 3D-SPARKLING sequence and using a TR shot = TR probe = 50 ms.The data plotted in this figure is acquired in vitro (phantom data) and comes from a single probe.It illustrates the temporal evolution of the 19 F NMR signal in the absence of participant-induced perturbations: Stimulated echoes due to the very short TR probe prevent steady-state at the level of the probes in the absence of a spoiling gradient(blue trace) whereas an external 470 mT*ms/m spoiling gradient overcomes this issue (orange trace).Therefore, It becomes possible to use the field monitoring system in this alternative setting to challenge its long TR probe constraint by simply applying an external spoiling gradient that is strong enough to eradicate the residual transverse magnetization.
However, according to the signal equation at steady state, even at the Ernst angle, using a smaller TR induces a lower NMR signal and, therefore, a degraded SNR at the level of the probes.The following subsections will demonstrate that even with a degraded SNR, the field camera manages to estimate field fluctuations that are accurate enough to yield beneficial correction for fMRI data.

Enhanced image quality
Figure 2 shows a comparison of the mean images computed from the resting-state fMRI data collected in volunteer #1 (V#1) and reconstructed according to the strategies (a) to (f) mentioned in Section 3.4.It is noteworthy that the image quality associated with strategy (a) is very degraded.Additionally, the figure demonstrates that the overall T * 2 contrast is significantly enhanced when correction of both static and dynamic fluctuations of the B 0 field is performed: The blue arrows illustrate how the correction of either the dynamic or static contribution alone (strategy (c) or (d)) resulted in partial improvement whereas correcting both terms (strategy (f)) induced a more significant enhancement.Moreover, the lost signal in strategy (a) is now recovered, and anatomical details are more finely reconstructed as depicted by the orange and green arrows respectively.Again, it is worth noting that strategies (b)-(e) yield limited improvement when compared to strategy (f) suggesting that none of the field terms has a negligible influence on image quality.
Figure 3A displays a subjectwise comparison between the mean images from the resting-state sequence of fMRI scans reconstructed without and with static and up-to-the-first-order dynamic field terms correction: The gain in image quality is systematic across the three volunteers and for three physiological movement scenarios namely NB, FB, and HC.

Increased tSNR
Figure 3B qualitatively depicts the boost in tSNR obtained when correcting static and up-to-the-first-order dynamic field terms: A notable increase is observed across the three volunteers for the NB scenario, 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.V#2, however, yielded slightly superior tSNR maps compared to the two other participants.This is likely due to the fact that they had less-energetic breathing fluctuations, as demonstrated in Figure S1A (in Section S3) where the power spectra of the breathing-induced field fluctuations during the NB scenario are shown for each subject.Despite a discernible increase in tSNR in the data collected in V#3 during FB and when performing the HC, we did not recover levels that are comparable with the NB scenario.This is likely due to the involvement of additional large-amplitude head Example of the NMR signal decay from one probe of the field camera over 96 free induction decays with/out an external spoiling gradient for a TR shot = TR probe = 50 ms: Stimulated echoes prevent steady-state in the absence of an external spoiling gradient whereas a 470 mT*ms/m spoiling gradient ensures a steady state at the level of the probes.

F I G U R E 2
Comparison of the mean image of resting-state functional MRI scans collected from V#1 reconstructed using strategies (a)-(f).From left to right, the top row (resp., bottom row) depicts the mean images yielded by the images reconstructed using strategies (a)-(c) (resp., (d)-(f)).The overall contrast is enhanced, the lost signal is better recovered, and anatomical details are better reconstructed as illustrated by the blue, orange, and green arrows, respectively.movements that we are currently unable to correct using our current protocol (cf, Figure S2 in Section S4).
These qualitative findings are quantitatively supported by Table 2: It 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 and scenarios, the relative gain for V#2 is approximately one-half (resp., one-third) lower than for the two other volunteers when strategies (b) and (c) (resp., strategies (e) and (f)) are used.Otherwise, the relative gain in median tSNR reaches a plateau around 30% at maximum.The limited improvement in the tSNR maps corresponding to the FB and HC scenarios (Figure 3B) highlights the limits of the correction in extreme cases.

Increased sensitivity to the BOLD contrast
Figure 4 illustrates a subjectwise comparison of the statistically significant activation patterns thresholded using alternatives (i) and (ii) and yielded by the uncorrected (strategy (a)) and fully corrected (strategy (f)) data (qActiv1) on selected axial slices.
First, a larger effect size is associated with B 0 imperfections correction for both statistical thresholding strategies.This is reproducible across the three volunteers but most visible in V#1 where, in addition to an increased statistical significance, the activation pattern seems to better delineate the gray matter in the visual cortex suggesting a finer fit to the cortical surface.
Second, comparing the results at two statistical control levels (i) and (ii), we observed a smaller discrepancy between the activation maps obtained when static and up-to-the-first-order dynamic field terms are corrected.This suggests that at a given false positive rate the sensitivity to the BOLD contrast is increased and that the results are not merely biased by false positives.
These qualitative findings are supported by the figures reported in Table 3, which summarizes the systematic gain in the number of activated voxels and maximum z-score values when B 0 field imperfections are corrected (strategy (f)): 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 greater boost in effect size when an FDR control is applied suggests once again a higher prevalence of true positives over false positives in the corrected fMRI scans.The reported figures are consistent between the first and third Comparison of the (A) mean images and (B) temporal signal-to-noise ratio (tSNR) maps yielded by the resting-state functional MRI scans sequence collected in the three volunteers and using the different physiological movement scenarios reconstructed using strategies (a) and (f).The improved image quality when strategy (f) is used is reproducible across volunteers.
participants.V#2, however, reveals fewer activated voxels, notably after applying FDR control.This is likely due to larger head movement amplitudes as shown in Figure S3 in Section S5 where the time course of motion regressors of the translation over the z-axis are depicted.

More accurate BOLD phase maps
Figure 5 depicts the BOLD phase maps derived from uncorrected and corrected fMRI volumes on selected axial views.Firstly, we noticed a larger spatial extent of the Gain in % of median tSNR in corrected data (strategies (b)-(f)) relative to the native temporal signal-to-noise ratio (tSNR) (uncorrected data, i.e. strategy (a)) computed over the brain mask.

Terms corrected volunteer
Note: The highest figures (in bold) are retrieved when strategy (f) is used.
BOLD phase maps when B 0 field imperfections are corrected.Second, the maps shown are in better agreement with the prior knowledge about the projection of the visual field onto the occipital cortex.A zoomed-in example is shown for V#1 as she/he was the most compliant volunteer: We clearly observe 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.Furthermore, and without ambiguity, the top (resp., bottom) parts of the visual field project onto the bottom (resp.top) parts of the occipital cortex.In comparison, the BOLD phase map yielded by the uncorrected data seems noisy and matches the expected gradient only partially.
To gain a better insight, the BOLD phase maps ought to be projected on the cortical surface as shown in Figure 6A (qPhase1).In the first column, we observe significantly enhanced BOLD phase maps when B 0 imperfections are corrected: A poor retinotopic organization is retrieved when no field term is corrected whereas an improved one can be inferred when they are.Furthermore, estimating the BOLD phase maps using an ROI defined using a p-value of 0.05 with FDR control (alternative (ii)) has a considerable demeaning effect on the phase maps yielded by the uncorrected fMRI volumes.However, the degradation is almost imperceptible on BOLD phase maps obtained from the corrected fMRI volumes.
Additionally, we compared in Figure 6B the BOLD phase maps derived from raw corrected (strategy(f)) and smoothed uncorrected (strategy(a)) fMRI volumes.Two distinct isotropic 3D Gaussian kernels (FWHM=1.1mmand FWHM = 1.5 mm) were used to smooth the raw fMRI volumes before performing a general linear model analysis.FWHM = 1.1 mm was tested because, all else being equal, it corresponds to an increase in tSNR of about 33%, which roughly corresponds to the maximum gain in tSNR we claim is possible to achieve using our protocol.Firstly, we observed that the smoothed uncorrected volumes yield BOLD phase maps with lower quality than the raw corrected scans.Secondly, we noticed that smoothing the images at FWHM = 1.5 mm does not recover the lost signal, whereas correcting B 0 field imperfections does (zoomed-in region in the left hemisphere).Furthermore, while it is true that we recover more spatially extended BOLD phase maps on the right hemisphere when smoothing the images, this map remains less spatially specific than in the raw corrected volumes.We conclude that the increase in median tSNR observed when correcting B 0 imperfections yields an increase in sensitivity that exceeds the expected improvement if we choose to degrade the spatial resolution to increase the tSNR.

DISCUSSION
In this work, the impact of static B 0 field inhomogeneities and dynamic B 0 field fluctuations retrospective correction on fMRI volumes has been assessed.

Main findings
Through the application of a suitable external spoiling gradient, we have first demonstrated that it is possible to use the field camera outside its standard experimental setting, which requires a quite long TR probe with respect to the TR shot in fMRI, and accurately estimate up-to-the-first-order dynamic field fluctuations.Second, on resting-state fMRI data, we have proved that B 0 field imperfections correction has a hugely beneficial impact on 3D-SPARKLING image quality as well as tSNR: Up to 30% increase in median tSNR was quantified.Finally, a rather exhaustive evaluation of the impact on the fMRI volumes demonstrated a significant improvement in the sensitivity and quality of the BOLD phase maps.

F I G U R E 5
Comparison of the blood oxygen level dependent (BOLD) phase maps over an ROI defined using the thresholding alternative (i) yielded by retinotopic funtionalMRI scans reconstructed without and with static and up-to-the-first order dynamic field terms correction.over the brain when performing FB (resp., HC) for 1mm 3 3D-EPI fMRI (20 repetitions) data at 7T. Upon comparison with 2 mm 3 fMRI data, they showed that the gain in tSNR is more significant at higher spatial resolution, that is, 1 mm 3 .The differences in the reported numbers arise from the experimental conditions (two dimensional vs. 3D, 3T vs. 7T, number of repetitions, volunteer's movement).
To the best of our knowledge, no existing study in the literature using a field camera pushed the comparison up to the statistical analysis.The increase in sensitivity associated with B 0 field imperfections correction we observed remains notable at stricter statistical thresholding levels as demonstrated by criteria QActiv1: The recovery of lost signal plays a major role in retrieving true positives.The enhanced image quality is also expected to yield better spatial specificity as blurring is extremely reduced.Such a claim is supported by the projections of the retinotopic maps onto the cortical surface for V#1.
The BOLD phase maps obtained from the corrected volumes were not merely better than those produced by the uncorrected ones, but their quality also exceeded that of those yielded by smoothing (using either FWHM = 1.1 mm or FWHM=1.5 mm) uncorrected volumes: In fact, when B 0 field imperfections correction is performed, confounding factors due to the participant's physiological movement and the system's instabilities and drifts are (A) Blood oxygen level dependent (BOLD) phase maps yielded by the data collected from V#1 and reconstructed without/with B 0 field imperfections correction and their projections on the cortical surface.(B) Projected BOLD phase maps yielded by the data collected from V#1 and yielded by the smoothed (using a Gaussian kernel with FWHM = 1.1 mm and FWHM = 1.5 mm) B 0 -uncorrected data and the raw B 0 -corrected data.eliminated prior to preprocessing and statistical analysis while larger voxel sizes only help in aggregating more signal within a voxel without effectively minimizing the effects of such confounding factors.

Is it relevant to go beyond first-order dynamic field fluctuations correction?
It worth noting that our protocol is feasible and has proved without ambiguity that static and up-to-the-first-order dynamic B 0 field imperfections correction is beneficial for 3D-SPARKLING fMRI statistical analysis by illustrating the gain in image quality, tSNR, and accuracy of the retonotopic mapping.Nevertheless, it is essential to keep in mind that data collected from 16 probes was used to estimate four unknown terms for each time point t, namely one zeroth-order term (k 0 (t)), and three terms translating the first-order field fluctuations (k(t) = (k x (t), k y (t), k z (t))).This means that the lower SNR of the probes expected from using a shorter TR probe is compensated with data redundancy.
A further step would be to correct the higher-order field terms 33 during MR image reconstruction.However, it involves an additional and heavy computational burden during image reconstruction, especially in the non-Cartesian setting.Depending on the encoding scheme, it would be necessary to consider the higher order terms 34 or not. 4,26In our case, the compromise between the possible benefit and computational load of such a strategy could be further investigated.Furthermore, additional investigation will be required to prove that the SNR at the levels of the probes remains high enough to estimate higher-order terms accurately as well.

Limitations and perspectives
In the literature, several Cartesian and non-Cartesian methods 3,4,27,35 seem to yield improved image quality compared to that reported in this work.However, it is actually challenging to perform relevant and fair comparisons as different experimental setups are used in the competing studies 3,4,27,35 in terms of brain coverage, number of shots, two-dimensional versus 3D acquisitions, etc.In a former work, 36 we performed an extensive comparison in fMRI between 3D-EPI and 3D-SPARKLING for similar acquisition parameters (same spatiotemporal resolution-1 mm 3 and 2.4 s-TR-and FOV = 192 × 192 × 128 mm 3 , same TE/TR shot = 20/50 ms and similar T obs ).
In this previous study, the dynamic field fluctuations were not corrected in 3D-SPARKLING fMRI volumes, and we found a poorer image quality as compared to 3D-EPI.
This is likely due to the fact that 3D-SPARKLING is more sensitive to the acquisition's imperfections than 3D-EPI.Indeed, in contrast to 3D-EPI where the same readout is repeated along the k-space planes, the random nature of 3D-SPARKLING trajectories make imperfections due to B 0 or motion accumulated along different orientations (cf, Figure S4 in Section S6).This effect directly translates into more complex artifacts in the image domain.However, we demonstrated in Reference 36 that despite a poorer image quality the sensitivity/specificity statistical trade-off for the detection of evoked brain activity based on the BOLD contrast is, on average over six volunteers, similar for both techniques (3D-SPARKLING vs. 3D-EPI).With the addition of dynamic field perturbations correction, we believe that we further improve the image quality for 3D-SPARKLING.
As demonstrated by our findings on resting-state fMRI data collected for the FB and HC scenarios, the fact that head movements are not corrected hinders the efficiency of the correction performed for the most extreme cases.In contrast, in typical cases of unvoluntary moderate movement such as for V#2, the findings showcase undeniable benefits.In Reference 37, the authors present an experimental protocol for using the field camera for motion monitoring in anatomical imaging.Such a solution should work for fMRI as well.However, it is beyond the scope of this paper.
In this work, we resort to external measurements of the dynamic field fluctuations, overlooking the fact that each shot crosses the center of the k-space in 3D-SPARKLING.In fact, such a feature enables self-navigation and the estimation of an average zeroth-order dynamic term ΔB 0,dyn per shot, similarly to the solution proposed for TUR-BINE in Reference 35.Nevertheless, such estimates may not be as accurate as the field camera measurements.In any case, it would be interesting to implement such a strategy and to compare its estimates with external measurements.
Furthermore, 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 MR physics of the experiment faithfully.In fact, a truthful model would consider static and dynamic B 0 imperfections as evolving jointly: A ΔB 0 map would be estimated for each volume in this case analogously to what is proposed in References 38 and 39.As demonstrated in Reference 29, it would actually be possible to estimate a volume-wise ΔB 0 map from 3D-SPARKLING data.Such a solution would also avoid any possible inconsistencies between an external ΔB 0 map and the fMRI volumes in case of patient movement for instance.However, the pipeline in Reference 29 was conceptualized for anatomical scans and would need further tuning for fMRI data since, to date, it does not produce good estimates of ΔB 0 maps from highly accelerated fMRI data.

CONCLUSIONS
In this work, we demonstrated a systematic and significant benefit in image quality, tSNR, and in terms of effect size and accuracy of the brain activity detection and localization when static and dynamic B 0 field imperfections are corrected retrospectively during MR image reconstruction: In fact, using retinotopic fMRI data, we have noted an increase in sensitivity, notably at a stricter false positive rate control level, and more accurate BOLD phase maps.This study was conducted using a field camera in an alternative setting challenging its TR probe constraint to monitor up-to-the-first-order dynamic B 0 field fluctuations.

( a )
No correction: None of the field terms were taken into account.(b) ΔB 0,dyn : Only the zeroth-order dynamic contribution was accounted for by the correcting the data itself.(c) ΔB 0,dyn & k: The zeroth-and first-order dynamic field fluctuations were corrected.(d) ΔB 0,stat : Only the static inhomogeneities were corrected.(e) ΔB 0,stat & ΔB 0,dyn : Static inhomogeneities and zeroth-order dynamic fluctuations were corrected.(f) ΔB 0,stat & ΔB 0,dyn & k: Static and up-to-the-first-order dynamic contributions were included in the signal model for reconstruction.

Notes:
The activated voxels are defined using two distinct statistical significance levels: a p-value of 0.001 without multiple comparisons correction and a p-value of 0.05 with FDR control.The highest figures (in bold) are obtained when strategy (f) -Full Correction-is used.V#2 reveals the lowest statistical significance.
Number of activated voxels and the maximum z-score values extracted from task-based fMRI volumes with/out correcting ΔB 0 imperfections.