Submillimeter balanced SSFP BOLD–functional MRI accelerated with 3D stack‐of‐spirals at 9.4 T

This work aims to improve the speed of balanced SSFP (bSSFP) acquisition with segmented 3D stack‐of‐spirals for functional brain studies at ultrahigh field.


INTRODUCTION
The speed, sensitivity, and availability of gradient-echo (GRE) EPI makes it an ideal candidate for most BOLD functional experiments.As a result of shorter T * 2 at ultrahigh fields (>7 T), submillimeter functional MRI (fMRI) more often relies on a segmented acquisition.The longer T 1 relaxation decreases either the acquisition speed or the SNR.In addition, longer echo trains result in signal loss, blurring, and distortion depending on the encoding strategy.Therefore, in terms of SNR, there is an advantage for unspoiled balanced SSFP (bSSFP) acquisition where a high ADC duty cycle (i.e., percentage of total acquisition time spent acquiring data) acquisition is possible even with shorter echo trains.Furthermore, the spin-echo behavior of the bSSFP BOLD contrast is believed to be more selective to the smaller vessels and therefore closer to neuronal activation. 1 Among other spin-echo acquisitions, 2 bSSFP is particularly well suited for high-field applications because of its flexibility and lower specific absorption rate (SAR) requirements.The main challenge of bSSFP at high field is to maintain the B 0 homogeneity within the imaging volume to avoid banding artifacts.Ironically, the bSSFP acquisition was originally proposed to increase the sensitivity of BOLD fMRI with its large signal change in the transition band. 3,4However, the present study was performed at a higher field with a relatively longer TR.Therefore, a more stable and less sensitive technique to detect BOLD signal changes in the passband is adopted. 5lthough spiral acquisition has been proposed for fMRI since its inception, 6,7 its adoption has been slow because it is more susceptible to system imperfections and requires a more complex reconstruction due to its non-Cartesian sampling pattern.It is particularly challenging at high fields with increased B 0 inhomogeneity and signal instabilities.With the advances in correction algorithms, [8][9][10] field monitoring, 11 open source tools, 12,13 and inexpensive computing, the image quality and speed of non-Cartesian acquisition has improved significantly.Recently, a high-resolution 2D gradient-spoiled single-shot spiral fMRI was demonstrated at 7 T with field monitoring and improved correction techniques for system imperfections. 14Spiral acquisition with bSSFP is popular for dynamic cardiac imaging at lower fields. 15owever, ultrahigh-field spiral bSSFP fMRI remains an unexplored domain, with only one documented study conducted at 3 T. 16 In previous functional studies, highly segmented EPI 17,18 and spiral readouts 16 have been proposed to improve the speed of bSSFP acquisition.In this study and our reported preliminary findings, 19,20 we aim to further improve speed and image quality by combining parallel imaging and correction techniques for system imperfections.Several functional studies with a visual task were performed to investigate the contrast and speed improvement of a 3D stack-of-spirals bSSFP acquisition.The dependence of bSSFP BOLD contrast on flip angle (FA), TR, and TE has been investigated in simulation 1,21,22 and in vivo. 17,18,23However, the isolated effect of TE is difficult to assess in a conventional bSSFP acquisition, because TE is usually placed at or close to TR/2 for sampling efficiency reasons.As the rapid T * 2 decay is quite significant at ultrahigh fields, we attempt to measure the effect of TE with minimal bias within a range relevant to our experiments (2-6 ms).We analyze relative signal change, significance of activation, and thermal and temporal signal stability to explain the effect.Similar experiments and analyses for TR dependence (6-10 ms) are performed to compare with the previously reported results. 18Finally, two submillimeter resolution protocols (0.6 mm isotropic and 0.8-mm isotropic) and a whole-brain coverage protocol at 1.2-mm isotropic are measured to demonstrate the speed, flexibility, and robustness of the proposed method.

Experimental setup
All experiments were carried out using a Magnetom 9.4T Plus whole-body MRI scanner(Siemens Healthineers, Erlangen, Germany) equipped with a SC72 gradient system (70-mT maximum gradient strength and 200-mT/m/ms slew rate) and an in-house-built 16-transmit/31-receive RF coil array. 24ight healthy subjects (5 male, mean age 34.8 ± 8.9 years) participated in the measurement after providing informed consent as per guidelines set by our local ethics committee.Six subjects were scanned for the TR/TE dependence and high-resolution protocols in a single session to reduce the interscan variability.Each session consisted of seven functional scans, a B 0 field mapping scan, and a T 1 -weighted anatomical scan.Two more volunteers took part in the whole-brain functional studies.To maximize the coverage of the visual cortex, the acquisition volume was oriented parallel to the calcarine sulcus in all protocols except the whole-brain protocol.
A full visual field checkerboard stimulus was presented in a simple block design starting with 20 s of resting period and followed by 10 s of stimulus presentation.Twelve cycles of the stimulus were presented during the individual functional experiments lasting about 6 min.The subjects were instructed to count the occurrences of a particular fixation point color change to sustain and check their attention level.The visual stimulus was designed in Psychopy 3.0. 25The trigger timings and stimulus state were logged and converted into a paradigm file to reduce timing and human errors.

Acquisition
A slab-selective segmented 3D stack-of-spirals sequence was used for all measurements.A uniform-density Archimedean spiral was generated according to the FOV, resolution, segmentation, and readout time.Then, the time optimal trajectory was calculated using the optimal control algorithm. 26The gradient constraints are the maximum available gradient strength per channel of 42 mT/m and a slew rate limit of 150-170 mT/m/ms, depending on the peripheral nerve stimulation limit.The important acquisition parameters of all the protocols are summarized in Table 1.All TR/TE scans were performed with a short 2.4-ms readout to reduce the T * 2 bias.The other protocols were optimized to maximize spatio-temporal resolution.In case of whole-brain functional experiment, the TR was reduced to 5 ms to increase the usable volume coverage (i.e., to reduce banding artifacts).Even at short TR, an ADC duty cycle (i.e., percentage of total acquisition time spent acquiring data) of about 63% was maintained using the shorter first-order water excitation.Because of the thick slab size, an additional 2-fold undersampling was possible in the slice direction resulting in a total 8-fold parallel-imaging acceleration.We chose a regular undersampling pattern to reduce the additional eddy current effects due to reordering. 27Additionally, we performed a whole-brain T 1 weighted MPRAGE 28 scan at 0.8-mm isotropic for anatomical reference (other sequence parameters: segment TR = 3.36 s, TR = 6.16 ms, TE = 3.03 ms, matrix size = 264 × 264 × 224, GRAPPA = 2 × 2, TI = 1.34 s, FA = 9 • ) and a 3D saturated TurboFLASH scan 29,30 at 3.5-mm isotropic resolution for FA mapping.Further information about the measurement protocols and the stimulus can be found in the shared software repository.

Water excitation
For rapid water excitation, two half FA slab-selective pulses (binomial-11 31 ) with bipolar gradients were used as illustrated in Figure 1.The shorter first-order water excitation, in which the fat processes  rad between the pulses, was used for the whole head protocol.The second-order water excitation with longer interpulse duration (1050 ), in which the fat precesses 3 rad, was used for the submillimeter protocols with thin slabs to accommodate the increased gradient switching needed for higher resolutions.

Gradient correction
The gradient impulse response function (GIRF) was measured using a field camera with 16 field probes (Clip-on camera; Skope, Zurich, Switzerland) arranged optimally in a frame to measure dynamic field changes up to second-order spherical harmonics. 32A protocol of 20 triangular pulses of duration 20 − 210 , slew rate of 200 mT/m/s, TR of 2 s, and 50 averages was used.Along with the GIRF measurement, various spiral waveform of different length and gradient axis were acquired for validating the measured GIRF.The calculated self-responses and cross-responses were low-pass-filtered with 60-kHz and 40-kHz FWHM to reduce noise amplifications.The B 0 temporal modulation of the MR system eddy current compensation was estimated and removed retrospectively.
The GIRF prediction up to first-order spherical harmonics including the cross-terms was used for reconstruction.The estimation of vendor B 0 modulations and gradient predictions was based solely on the spiral readout gradients, with no inclusion of other gradients in the prediction process.
In a separate measurement, four field probes were mounted outside the RF coil, and a frequency-modulated RF pulse from the field probes was recorded by both the MR system and the field camera system to measure the data delay (a combination of filter group delay 33 and propagation delay) between them.This synchronization scan was repeated for different ADC dwell times, and a simple linear regression model of data delay as a function of ADC dwell time was calculated.This allows proper synchronization of the data and the corrected trajectories.

Image reconstruction
The image reconstruction pipeline shown in Figure 1 was developed using MATLAB (MathWorks, Natick, MA, USA) and is similar to the established model-based iterative reconstruction. 9,34Before reconstruction, the dependencies for the reconstruction like B 0 map and coil sensitivity maps were calculated.The B 0 map was calculated from an accelerated 3D multi-echo gradient-echo (GRE) acquisition using a spatial unwrap algorithm 35 and a regularized field map estimator 36 for denoising after reconstruction.The ESPIRiT implementation of the BART toolbox was used to calculate 3D coil maps from the noise-decorrelated calibration data. 12Computationally efficient forward and adjoint operators for non-Cartesian gridding were formulated using the nonuniform fast Fourier transform algorithm using min-max interpolation 13 and the time segmentation algorithm for fast off-resonance correction. 9,37hree-dimensional nonuniform fast Fourier transform  was used, as the trajectory in slice direction does not follow Cartesian sampling after GIRF correction.An iterative conjugate-gradient SENSE method was used to invert the signal model to obtain the reconstructed images. 10With the density compensation function 38 preconditioning, an earlier convergence was achieved after, at most, 10 iterations for all protocols.Additionally, Tikhonov regularization with a heuristically chosen regularization factor based on visual inspection was used to improve the stability of the ill-posed model inversion.The individual functional scans were split into multiple parts and reconstructed simultaneously across various compute nodes of a high-performance compute cluster to increase data throughput.The reconstruction time was approximately 4.8, 20, 18, and 98 min for 1.0, 0.8, 0.6, and 1.2 mm data sets, respectively, and showed high variability depending on the number of time segments required for B 0 correction.
In addition to the offline reconstruction, a simplified real-time reconstruction was implemented in the Gadgetron framework 39 to assess image quality during a scan session.In case of banding artifacts in the region of interest (ROI), additional shimming with reduced shim volume and additional system frequency adjustments were performed.

Data and fMRI analysis
The reconstructed volumes were motion-corrected using the AFNI toolbox 40 (motion parameters can be found in Figure S1), and a general linear modeling was performed using the FSL package. 41Only the minimal preprocessing steps like high-pass temporal filter (0.03 Hz) for baseline correction and a mild spatial smoothing (FWHM = nominal resolution) for robust activations were performed.Other default settings, including temporal derivative to accommodate for timing inaccuracies and cluster thresholding for multiple comparison correction, were left untouched.No other nuisance regressors like motion parameters or physiological traces were used.The anatomical MPRAGE volume was coregistered and resliced to the functional volume space using the SPM package. 42he tSNR was calculated by dividing average signal intensity by the SD of the motion-corrected volumes.To reduce BOLD contamination, only the volumes acquired during the last 14 s of the resting period were used here.The regression coefficients from the general linear modeling analysis were normalized with the average signal time course to obtain the relative signal change (s) in percent.Additionally, SNR 0 (thermal signal to noise ratio) and g-factor maps were calculated for all protocols using the pseudo-replica method. 43All reported data metrics like tSNR efficiency (tSNR divided by the square root of the acquisition time), B 0 inhomogeneity, z-scores, and signal change were calculated within a brain mask after motion correction and removing 20% of the outermost slices corrupted by a poor slice profile (i.e.20% slice oversampling was applied).For z-score and signal change, only the activated voxels (z-score > 4) within the brain mask were used.tSNR efficiency metrics were estimated from the occipital region to eliminate bias from regions with banding artifacts (ROI used for tSNR efficiency estimation is highlighted in Figure S2).

RESULTS
The GIRF correction was validated by comparing the corrected trajectory to the trajectory measured with the field probe system.The average RMS error of the delay-corrected nominal and measured gradient waveforms was approximately 1% of the maximum gradient amplitude.After GIRF correction, the average RMS error was reduced to about 0.1% of the maximum gradient amplitude, regardless of the gradient rotation or readout length (Figure S3).The resulting corrected trajectory showed only a deviation of less than 5 rad/m for the short readouts.

Water excitation
The spectral selectivity of short and long water excitation pulses along with their performance in whole brain and thin slab acquisition are exhibited in Figure 2. In both cases, binomial excitation significantly reduces fat artifacts.The effect is more pronounced in the long water excitation case applied in the high-resolution protocols, in which the extended readout necessitated by the higher resolution significantly blurs the fat signal.The poor B 0 homogeneity results in suboptimal water excitation near the frontal regions.After transmitter adjustments, the actual FA in the visual cortex is about 70%-80% of the nominal FA.

Image quality and reconstruction
All images shown in this work were reconstructed using the described image reconstruction pipeline with all the corrections for system imperfections.The isolated improvements in image quality from the gradient and B 0 correction are illustrated in Figure S4.The reconstructed volumes show typical bSSFP contrast with minimal gray-matter and white-matter contrast because of their similar T 2 /T 1 ratio at high field (see Figure S5).The SNR 0 and g-factor for all protocols from a single subject (Subject 7 for whole head and Subject 2 for the rest of the protocols) were calculated from 128 pseudo replicas.Additionally, tSNR of the same data set was calculated and is shown in Figure S2.To compare protocols with different acquisition time, both SNR 0 and tSNR are displayed as efficiency maps (i.e., SNR 0 and tSNR divided by the square root of the acquisition time).displayed along with key statistics such as median signal change, number of activated voxels, median z-scores, and tSNR to elucidate the relationship.There is a significant increase in median signal change of 0.89 percentage points for a 4-ms TE increment from 3.69% at 2-ms TE.We observed a slight drop in tSNR and a fairly constant median z-scores and number of activated voxels with increasing TE.The shim volume was restricted to the visual cortex for some subjects to reduce bandings in the ROI, resulting in severe bandings in the frontal regions.Similarly, TR dependence of bSSFP BOLD among the same subject group is shown in Figure 4.In the rather small range of 6-10-ms TR (TE = 2 ms, 2.4-ms readout duration), the median signal change increased by 0.78 percentage points, from 2.91% at 6 ms TR.As a result of the fixed readout length, the ADC duty cycle dropped from 40% to 24% with increasing TR (6-10 ms).We observed a 25% drop in tSNR and a mild decrease in the extent and statistical significance of the functional activation in the same TR range.

Optimized acquisition protocols
The speed, coverage, and flexibility of this acquisition method was explored using three optimized protocols with maximum possible ADC duty cycle (Table 1).In all cases, distortion-free and highly significant activation along the gray-matter gyri were observed.Figure 5 presents the functional activation for full field visual stimulus at 0.8-mm isotropic on the average functional volume and coregistered anatomical MPRAGE volume.The geometric accuracy of the functional volume is highlighted with the overlay of gray-matter outlines.The orthogonal sections of activation maps from other subjects are also shown.In a similar manner, the results of the first 0.6-mm isotropic functional bSSFP experiment is shown in Figure 6.Statistics such as median signal change, z-score, tSNR, and B 0 homogeneity across all subjects for both high-resolution protocols are displayed in Figure 7 (additional summary metrics are provided in the Table S1).The result of a single-subject, whole-brain acquisition at 1.2-mm isotropic is shown in Figure 8.In this case, a median signal change of 2.4% was observed among 8957 activated voxels (z-score > 4).A much higher mean tSNR efficiency of 20.03 s −1/2 was obtained within the ROI.

DISCUSSION
We have demonstrated the feasibility of passband bSSFP with an accelerated stack-of-spirals readout at 9.4 T. This study achieved at least comparable and often higher spatio-temporal resolution than many of the GRE functional studies (see Table S2) for the following reasons: (1) high ADC duty cycle combined with spiral encoding; (2) the weaker TE dependence leading to more flexible and efficient encoding; (3) no need for a time-consuming fat saturation module due to the use of fast water excitation; and (4) no need for navigator modules for intershot phase correction due to bSSFP's high signal stability toward physiological noise compared with a 3D-GRE acquisition. 17

Water excitation
Despite the relatively short spiral readout (≤ 6.5 ms), the large chemical shift of fat at 9.4 T (≈1400 Hz) significantly degrades the image quality.The off-resonant fat signal blurrs in the in-plane direction due to low-bandwidth spiral readouts and translates along the slice direction during slice selection.This mislocalization of fat was not accounted for in the signal model, and hence affects the convergence of the reconstruction.The water excitation Visual activation maps overlaid onto the orthogonal slices of the 0.8-mm isotropic mean spiral volume acquired in 1.96 s and the anatomical MPRAGE volume of Subject 4. Gray-matter boundaries from the structural volume superimposed on the functional volume highlight the lack of distortion.The bottom panel shows the zoomed-in anatomical scan with activations from 4 additional subjects.

F I G U R E 6
Visual activation maps overlaid onto the orthogonal sections of the 0.6-mm isotropic mean spiral volume acquired in 2.16 s and the anatomical volume of Subject 4. Gray-matter boundaries from the structural volume superimposed on the functional volume highlight the lack of distortion.The bottom panel shows the zoomed-in anatomical scan with activations from 4 additional subjects.
module that was included in the sequence mitigated these issues without disturbing the steady state, unlike conventional off-resonant fat saturation.An image-degrading effect of the narrower spectral selection of the long water excitation compared with short water excitation was not observed in the images due to the better shim achievable in thinner slabs.The use of faster bipolar gradients with the binomial pulse results in poor water excitation in the slab boundaries.Therefore, a relatively high bandwidth pulse to reduce the fat displacement in the slice direction was preferred over a longer pulse to reduce SAR.Alternatively, fat can be suppressed using an alternating TR bSSFP acquisition 44 or through joint fat-water estimation 45 in combination with a multi-echo readout, although this comes at the cost of acquisition efficiency.

Image quality
The model-based reconstruction improved the image quality by accounting for system imperfections such as trajectory deviation and spatial B 0 inhomogeneity.The GIRF correction primarily reduces the blurring and distortion of the high-resolution features, because larger trajectory  S1. deviations occur primarily in outer k-space.Additionally, the GIRF correction improves the registration of the spiral image with the Cartesian field map and coil sensitivity maps for better field-inhomogeneity correction and parallel-imaging performance.Further improvements in GIRF predictions can be achieved by considering gradients before the spiral readout gradients (see Figure S6).This requires additional sequence simulations or exporting all gradient definitions during the measurement.
From the point spread function analysis, the limited k-space support increases the FWHM to 1.21 and 1.40 times the nominal resolution in slice and in-plane directions, respectively. 46The T * 2 decay along the short readout (6.5 ms) minimally increases in-plane FWHM by approximately 5%.Static spatial B 0 inhomogeneity has minimal influence after B 0 correction even in the presence of small errors in the field-map estimation.Because the spiral readouts are relatively short (≤ 6.5 ms), the images show almost no blurring in the passband region.The geometric accuracy improvement due to the B 0 correction is primarily visible in the higher-order (off-resonant) passband regions with large B 0 inhomogeneity.Therefore, the residual blurring of the functional volumes could be the result of unbalanced higher-order gradient moments (see Figure S7) disturbing the geometric congruence of the echo pathways. 27These imperfections cannot be corrected retrospectively but must be minimized during acquisition.Furthermore, factors such as subject motion, the nonlinear nature of global B 0 fluctuation, 8 and higher-order spatial dynamics also degrade image quality.

TE/TR dependence
The influence of TE, TR, and FA on bSSFP-BOLD contrast in terms of signal change and vessel-size specificity has been investigated extensively in simulations 1,21,22,47 and partially in measurements. 17,18,23In general, BOLD signal change can be increased to some extent by increasing any of these three sequence parameters.However, because of the SAR limitation at 9.4 T and to maintain acquisition speed, the optimal FA for bSSFP BOLD for a given TR is almost never reached.Therefore, more often the TR dictates the maximum possible FA.Also, with the TE ≈ TR/2 constraint in conventional Cartesian bSSFP studies, TR becomes the primary factor determining the in vivo contrast at higher fields.Previously, an increase in group mean signal change of approximately 2.1 percentage points at 9.4 T was reported for a similar parameter range, representing the combined effects of TR (4-10 ms), TE (2-5 ms), and echo train length (0.8-7 ms). 18In this study, we attempted to disentangle the influence of TE and TR separately with minimal T * 2 bias.We observed a group average increase in signal change of 0.78 percentage points for increasing TR (6-10 ms) and 0.9 percentage points for increasing TE (2-6 ms) for a constant readout length of 2.4 ms.Although the TE dependence of the bSSFP BOLD is much weaker when compared with GRE BOLD, the effect is quite significant at 9.4 T due to the short T * 2 (≈18 ms).The sensitivity of the functional experiments also depends on the temporal signal stability.Passband bSSFP generally achieves higher tSNR than GRE methods because of higher SNR 0 and being less sensitive to physiological noise at short TR. 17 It is worth noting that the SNR gain due to unspoiled transverse magnetization in bSSFP acquisition at high fields is minimal due to the diminishing T 2 /T 1 ratios compared with lower fields.Nevertheless, the high duty-cycle readout and shorter TEs improve SNR 0 across all resolutions (see Figure S2) when compared with previously reported values. 18,48The dependence of noise characteristics on the details of the regularized iterative reconstruction makes a quantitative SNR comparison difficult.The lower values in the g-factor maps (see Figure S2), which predict spatially varying noise amplification due to parallel imaging, suggest optimal performance of our receiver array.However, a corresponding increase in the tSNR efficiency maps was not observed, possibly due to considerable physiological noise at higher fields and stronger eddy current-related signal instabilities.
The sensitivity gain from signal change at large TE was lost with decreasing tSNR, as indicated by the flat response of the median z-score and activation extent with respect to TE (see Figure 3).Because of the fixed readout length of the TR dependency scans in this study, the experiments with large TR have an unfair SNR 0 and tSNR disadvantage.Therefore, even with the increasing signal change of 0.8%, there was a mild decrease in the extent and statistical significance of the functional activation.This can be attributed to dropping tSNR efficiency at higher TR because of the lower acquisition efficiency and possibly to the decreasing signal stability.Generally, increasing TR allows both higher FA and acquisition efficiency, which should synergistically improve the sensitivity at the expense of coverage loss due to bandings, as shown in previous studies. 17,18egarding specificity, a TE in the middle of the RF pulses (i.e., TE = TR/2) and a lower TR are more desirable for large-vessel suppression. 1,22The refocusing component responsible for vessel-size specificity reaches its maximum amplitude near the middle of the RF pulses. 22,49In the simulations similar to Scheffler et al. 1 (see Figure S8), no significant degradation of "vessel size" specificity with decreasing TE was observed.Additionally, the spin-echo behavior of bSSFP acquisition is only guaranteed if the intravoxel spectrum is within the banding width (1/TR) and the T 2 relaxation within the repetitions is minimal (TR ≪ T 2 ).All measurements were acquired with a TR ≤ 10 ms, which should be comfortably in the bSSFP regime (TR ≤ 15 ms at 9.4 T). 48

High-resolution functional scans
In this work, a spiral-out trajectory with small TE was adopted rather than a spiral-in/out trajectories (with TE close to TR/2) 19 for better acquisition speed and tSNR.In all the used high-acquisition-efficiency protocols, an ADC duty cycle of 60%-65% was maintained regardless of resolution and TR, resulting in high thermal SNR and spatiotemporal resolution.The minimal overhead of the bSSFP acquisition makes it possible to obtain high acquisition efficiency with a short readout.Given the hardware limits, the acquisition time shows an asymptotic trend with respect to the TR (see Figure S9).For resolutions up to 0.8-mm isotropic, there is only a marginal improvement in acquisition speed by prolonging the TR (and thus the readout) beyond 10 ms with our gradient system.
To achieve the same functional sensitivity as the GRE techniques, the bSSFP acquisition requires approximately twice the tSNR, 50 because the bSSFP-BOLD signal change is approximately half that of the GRE BOLD at 9.4 T. 48 We obtained a slightly higher average tSNR efficiency of 10 s −1/2 for the 0.8-mm isotropic protocol across all subjects compared with a mean tSNR efficiency of 8.4 s −1/2 for 2D-GRE spiral at 7 T for the same resolution. 14Previous EPI-bSSFP experiments at 9.4 T reported a much higher tSNR efficiency of approximately 26 s −1/2 for 0.9-mm isotropic resolution. 18The absence of expected tSNR improvements due to higher SNR 0 appears to be due to the reconstruction and/or the spiral acquisition and must be investigated further.The observed tSNR together with mild spatial smoothing (FWHM = resolution) was sufficient to robustly detect activations in the ROI across subjects.The detected functional activation was primarily localized on the gray-matter ribbon without any visible distortion (compare Figures 5 and 6).Due to the generic nature of the stimulus and the mild spatial smoothing, it is difficult to quantify the spatial specificity in the high-resolution functional experiments.Further investigation of the underlying bSSFP-BOLD contrast and its applicability for cortical layer studies is warranted in a dedicated multirun functional experiment with a layer-specific paradigm.For the voxels with large blood fraction (i.e., large veins), the signal change is dominated by apparent T 2 relaxation due to water exchange pathways rather than susceptibility-related contrast. 47As a result, highly significant activations due to intravascular effects were observed near the large veins.
In the whole-brain experiment, TR was reduced to 5 ms to increase the banding-free volume.This relaxes the B 0 homogeneity requirement to approximately 140 Hz (70% of the bSSFP band spacing).However, we observed banding artifacts near the temporal and frontal regions with a conventional second-order spherical shim.Additionally, the combination of increased proximity and poor signal stability in the banding regions resulted in strong parallel-imaging artifacts near the ear canals.Nevertheless, a higher tSNR was achieved in the ROI facilitating the detection of small signal change activations.
Despite all the challenges of both passband bSSFP and an ultrahigh-field spiral acquisition, the bSSFP-BOLD contrast shows great potential at high field because of the increased extravascular effects of smaller vessels. 48ith its flexible 3D acquisition, this method can provide isotropic high-resolution functional images free from distortions for investigating cortical layers or columns, and high temporal resolution for event-related functional studies.Moreover, when compared with T * 2 -based approaches, bSSFP BOLD potentially shows less signal bias (due to extravascular effects) from large vessels 51 in the superficial layer, thereby enhancing the clarity of layer activations.

CONCLUSIONS
We have demonstrated the capability of passband bSSFP fMRI when combined with accelerated 3D stack-of-spirals readouts at ultrahigh field.

s). (B)
The effects of these trajectory deviations with and without k 0 on the reconstruction quality along with their difference (in %).In this study, gradient impulse response function (GIRF) predictions performed with only the spiral gradients active during the readout are used for reconstruction (including k 0 ).The accelerated 0.8-mm isotropic protocol is used for sequence simulation, and the reconstructed images are from the functional scan performed on Subject 4. Figure S7.Additional phase errors in degrees due to higher-order eddy currents (second and third orders only) predicted by gradient impulse response function (GIRF) for the 0.8-mm isotropic protocol.The phase accrual between the RF pulses at eight different spatial locations (see colored markers on the image) for all the excitations (7 × 28) of a single volume.The distance of the corresponding voxel from the gradient isocenter is shown in the title of panels.All gradient events within 500 ms are considered for the GIRF prediction.The following parameters were used: G max = 42 mT/m, S max = 180 mT/m/ms, FOV = 210 × 210 mm 2 , and ADC dead times between 1.5 and 4.5 ms for resolutions between 2 mm and 0.5 mm, respectively.These requirements can be further reduced by about 3-4 times with parallel imaging.Table S1.Subject-specific measurements to assess the experimental conditions and the functional sensitivity of the high-resolution functional MRI experiments.Abbreviation: P 95% ΔS, 95th percentile signal change.Table S2.Comparison of acquisition efficiency (matrix size/[volume TR * undersampling factor]) of previous passband balanced SSFP (bSSFP) functional studies and some gradient-echo (GRE) functional studies.Acceleration factors due to parallel imaging and partial Fourier are not included in the acquisition efficiency calculation to compare acquisition with different coverage.This table is adopted for bSSFP functional studies from a similar table comparing acquisition efficiency of 2D-GRE spiral studies presented in the supporting information of Kasper et Abbreviations: Dur., RF pulse duration; FA, flip angle; N PE1 , segmentation in first phase encoding directions (Interleaves); N PE2 , segmentation in second phase encoding directions (partitions); N vol , number of acquired volumes; Res., isotropic nominal resolution; R in-plane , in-plane under sampling factor; RO, spiral readout length; R slice , undersampling factor along slice direction; TBWP, time bandwidth product; TR vol ,

Figure 3 2
Figure3shows the influence of TE on bSSFP BOLD contrast across 6 subjects.The center slices of all subjects for TE = 2/4/6 ms (TR = 10 ms, 2.4-ms readout duration) are

3
Results of TE dependence scans acquired at 1-mm isotropic resolution (pulse TR fixed at 10 ms).(A-C) The middle slices with activations from 6 subjects are presented for three TEs of 2, 4, and 6 ms.(D) The Venn diagrams show the unique and shared activated voxels (z-score > 4) between different TEs for all subjects.(E) The median relative signal change (%), median z-score, number of activated voxels (z-score > 4), and the mean temporal SNR (tSNR) efficiency ( tSNR/ √ volume TR ) with respect to TE for each subject are plotted together with their mean values.
C) Results of pulse TR dependence scans acquired at 1-mm isotropic resolution (TE fixed at 2 ms).The center slices with activations of 6 subjects are presented for three TRs: 6, 8, and 10 ms.(D) Venn diagrams show the unique and shared activated voxels (z-score > 4) between different TRs for all subjects.(E) The median relative signal change (%), median z-score, number of activated voxels (z-score > 4), and mean temporal SNR (tSNR) efficiency ( tSNR/ √ volume TR ) with respect to TR for each subject are plotted together with their mean values.

7
Subject-specific measures to assess experimental conditions and functional sensitivity of the high-resolution functional MRI (fMRI) experiments.The temporal SNR (tSNR) efficiency and spatial B 0 inhomogeneity distributions in the region of interest of each subject for the 0.8-mm isotropic protocol are shown in (A) and (B), respectively.(C,D) Relative signal change and z-scores of the activated voxel (z-score > 4) of each subject for the 0.8-mm isotropic protocol.(E,F) The corresponding measures for the 0.6-mm isotropic protocol.The SD of the distribution is marked in (B) and (F).The number of activated voxels is mentioned in (D) and (H).The mean (black bar), median (red bar), and 95th percentile (blue plus) are marked in the distribution along with the corresponding group-level mean (dotted line) and SD.The 95th percentile specifically captures the top 5% of voxels with the largest signal change.Additional measures of the high-resolution fMRI experiments are provided in Table

F I G U R E 8
Exemplary transversal, coronal, and sagittal slices of the 1.2-mm isotropic volume acquired in 2.88 s with functional activations.The lower panel shows the alignment of activations along the gray-matter ribbon in a coregistered MPRAGE volume.
L) of all measurements from Subject 2 and the whole-brain measurement from Subject 7 are shown.The regions of interest (ROIs) used for temporal SNR (tSNR) estimation are shown as blue contours.The mean and SD within the mask are displayed in the top-right corner of the image.The lower g-factor for high-resolution protocols with the same undersampling factor is a result of early stopping.

Figure S3 .
The measured gradient impulse response function (GIRF; A) and data delay with respect to the ADC dwell time (B) of our MRI system.(C) GIRF correction of a spiral-out gradient waveform.(D) Gradient amplitude deviations in mT/m before and after GIRF corrections.The solid and dashed lines are the x and y gradient axis, respectively.FigureS4.Reconstruction with gradient correction and off-resonance correction of fully sampled 0.6 mm isotropic from subject 4. The reconstruction with nominal (delay corrected trajectory), gradient impulse response function (GIRF) corrected trajectory, time segmentation (MTI), and frequency segmentation (MFI) deblurring algorithms.The amplified difference maps and the spatial B 0 map are shown in the bottom panel.The calculation of interpolation weights for combining the frequency segments is highly computation-intensive and memory-intensive with our current MFI implementation.Therefore, we used MTI in the reconstruction pipeline.

Figure S5 .
The mean high-resolution (0.8-mm isotropic and 0.6-mm isotropic) functional volumes from all subjects.Outer slices corrupted by poor slice selection and motion correction are not shown.FigureS6.(A) Trajectory deviations when considering only the spiral readout gradient and additional trajectory deviations due to the past gradient events (up to 0.5

Figure S8 .
(A,C,E) Simulated BOLD contrast with respect to sphere/cylinder radius for various TRs at a TE of 2 ms.(B,D,F) BOLD contrast for various TEs at a fixed TR of 10 ms.(A,B) Simulated with spheres and measured microsphere probe relaxation times (T 1 ∕T 2 ≈ 1726∕649 ms) at 3 T. (C,D) Simulated with cylinders and relaxation times (T 1 ∕T 2 ≈ 1310∕71 ms) at 3 T. (E,F) Simulated with cylinders and relaxation times (T 1 ∕T 2 ≈ 2000∕35 ms) at 9.4 T. The flip angle for all the simulations was 30 • .Figure S9.(A) The 3D stack-of-spirals sequence timings for different resolutions without parallel imaging.The acquisition efficiency (number of voxels/volume TR).(B) Acquisition time per slice (second phase encoding).(C) Number of excitations per slice for some isotropic resolutions with respect to TR.The typical acquisition efficiency range of conventional gradient-echo (GRE)-EPI sequences is shown as a green region for comparison.