Alpha-180 spin-echo-based line-scanning method for high-resolution laminar-specific fMRI in animals

Line-scanning fMRI has been successfully applied to investigate circuit-specific neuronal activity by measuring dynamic hemodynamic responses across cortical layers with high spatiotemporal resolution (Albers et al., 2018; S. Choi, Chen, et al., 2023; S. Choi et al., 2018; S. Choi, Xie, et al., 2023; S. Choi, Yu, et al., 2022; S. Choi, Zeng, et al., 2022; Raimondo, Priovoulos, et al., 2023; Raimondo et al., 2021; X. Yu et al., 2014). This is initially originated from Mansfield’s line-profile mapping studies in early 1970s (Mansfield & Maudsley, 1976; Mansfield et al., 1976). The advantage of the current line-scanning fMRI method is to sample cortical layers with ultra-high spatial resolution. Meanwhile, the line-scanning method only acquires a single k-space line per timepoint, enabling an ultrafast sampling rate. This high spatiotemporal laminar fMRI sampling scheme has been being utilized for bottom-up and top-down blood-oxygenation-level-dependent (BOLD) fMRI mappings in both animal and human fMRI studies. Previously, X. Yu et al. (2014) developed a line-scanning fMRI method to delineate laminar fMRI onset time with distinct laminar-specific neural inputs such as thalamocortical input and corticocortical input in the rat brain with high spatial (50 μm) and temporal resolution (50 ms). Line-scanning fMRI has also been combined with optogenetic control to further investigate the temporal features of the fast neural inputs across cortical layers in rodents (Albers et al., 2018). Beyond preclinical fMRI studies, line-scanning fMRI for human brain mapping has demonstrated a good correspondence with BOLD responses of 2D echo planar imaging (EPI) at the same temporal scale (200 ms) (Raimondo et al., 2021). This line-scanning fMRI also motivated the cortical depth-dependent diffusion-based fMRI mapping schemes (Balasubramanian et al., 2021). Lately, the ultra-fast line-scanning fMRI with k-t space reshuffling scheme has even provoked some interesting investigations of direct neuronal activity measurements (Tan Toi et al., 2022), albeit it has not been replicated in animals and humans (S.-H. Choi et al., 2023; Hodono et al., 2023).

The typical gradient echo (GRE)-based line-scanning fMRI (GELINE) method needs to dampen signals outside of the region of interest (ROI) to avoid aliasing artifacts along the phase encoding direction (X. Yu et al., 2014; Albers et al., 2018; Raimondo et al., 2021; S. Choi, Zeng, et al., 2022; Raimondo et al., 2023; S. Choi, Chen, et al., 2023; S. Choi, Xie, et al., 2023). Two saturation slices with additional RF exposure are applied for this purpose. However, two issues should be further investigated. One is the imperfect elimination of the aliasing artifacts (including inflow effects) due to imperfect RF performance and inhomogeneous B0 field. The other is the specific absorption rate (SAR) problem stemming from high duty cycle sequences. To alleviate large draining vein contribution to non-specific GELINE responses, large-tip-angle spin-echo based line-scanning fMRI method was previously proposed for diffusion-based fMRI (dfMRI) studies in animals by combining both a typical spin-echo sequence and saturation RF pulses at a short TR (Nunes et al., 2021). This study provided the first direct link between ultra-fast dfMRI signals upon forepaw stimulation and intrinsic optical signals upon optogenetic stimulation while demonstrating ultrafast dfMRI could reflect neuromorphological coupling. Even though the concept of the proposed spin-echo line-scanning fMRI method was implemented in humans, laminar-specific fMRI signals were not observed possibly due to low temporal signal-to-noise ratio (tSNR), potential motion artifacts (no motion correction applied), and/or low functional sensitivity (Raimondo, Heij, et al., 2023). Based on our initial line-scanning fMRI study (S. Choi et al., 2018), here, we further developed an α (alpha)-180 line-scanning fMRI method in animals to solve these problems while achieving high spatiotemporal resolution. We modified a spin-echo (SE) sequence by altering the refocusing 180° RF pulse perpendicular to the excitation slice (S. Choi et al., 2018; Mansfield et al., 1976; Mansfield & Maudsley, 1976). This adjustment allows to only highlight a specific line-profile across the cortical layers without the need for additional saturation RF pulses. Nevertheless, there is an inevitable trade-off between T2*-weighted GELINE and T2-weighted SELINE: GELINE has low specificity but high sensitivity to BOLD whereas SELINE has low sensitivity but high specificity to BOLD. As reported in previous works (Boxerman et al., 1995; Budde et al., 2014; Han et al., 2019; Norris, 2012; Siero et al., 2013; Zhao et al., 2004), GRE-based BOLD responses are dominated by the macrovasculature (e.g., large draining veins) and SE-based BOLD responses have greater microvascular sensitivity (e.g., capillary vessels) than GRE-based BOLD responses, specifically indicating extravascular and intravascular contributions to BOLD responses should be taken into account. In high-field MRI, extravascular signal changes likely predominate while intravascular signal changes mostly diminish (Boxerman et al., 1995). In contrast to the GELINE method, the SE-based line-scanning fMRI (SELINE) method thus has the potential to effectively exclude the surface draining vein effects. However, it should be noted that the laminar patterns of BOLD signals in SELINE can still be highly varied across different cortical layers in anesthetized rats. Furthermore, we can also shorten the repetition time (TR) to 200 ms for the SELINE method to sample the high-resolution T2-weighted fMRI signals, demonstrating the feasibility of the fast sampling of laminar fMRI with effective ROI selectivity in rodents.

The study was performed in accordance with the German Animal Welfare Act (TierSchG) and Animal Welfare Laboratory Animal Ordinance (TierSchVersV). This is in full compliance with the guidelines of the EU Directive on the protection of animals used for scientific purposes (2010/63/EU) and the MGH Guide for the Care and Use of Laboratory Animals. The study was reviewed by the ethics commission (§15 TierSchG) and approved by the state authority (Regierungspräsidium, Tübingen, Baden-Württemberg, Germany) and the MGH Institutional Animal Care and Use Committee (Charlestown, MA, USA). A 12-12 h on/off lighting cycle was maintained to assure undisturbed circadian rhythm. Food and water were available ad libitum. A total of four male Sprague–Dawley rats were used in this study.

Anesthesia was first induced in the animal with 5% isoflurane in the chamber. The anesthetized rat was intubated using a tracheal tube, and a mechanical ventilator (SAR-830, CWE, USA) was used to ventilate animals throughout the whole experiment. Femoral arterial and venous catheterization was performed with polyethylene tubing for blood sampling, drug administration, and constant blood pressure measurements. After the surgery, isoflurane was switched off, and a bolus of the anesthetic alpha-chloralose (80 mg/kg) was infused intravenously. After the animal was transferred to the MRI scanner, a mixture of alpha-chloralose (26.5 mg/kg/h) and pancuronium (2 mg/kg/h) was constantly infused to maintain the anesthesia and reduce motion artifacts.

All data sets from rats were acquired using a 14.1T/26 cm (Magnex, Oxford) horizontal bore magnet with an Avance III console (Bruker, Ettlingen) and a 12 cm diameter gradient system (100 G/cm, 150 μs rising time). A home-made RF transceiver surface coil with a 10 mm diameter was used on the rat brain. For the functional map of BOLD activation (Fig. 1A), a 3D gradient-echo EPI sequence was acquired with the following parameters: TR/TE 1500/11.5 ms, FOV 1.92 × 1.92 × 1.92 cm³, matrix size 48 × 48 × 48, spatial resolution 0.4 × 0.4 × 0.4 mm³, and readout bandwidth 133928 Hz. A high order (e.g., 2^nd or 3^rd order) shimming was applied to reduce the main magnetic field (B0) inhomogeneities at the region-of-interest. For anatomical reference of the activated BOLD map, a RARE sequence was applied to acquire 48 coronal images with the same geometry as that of the EPI images. The fMRI design paradigm for each trial comprised 200 dummy scans to reach steady-state, 10 pre-stimulation scans, 3 scans during stimulation, and 12 post-stimulation scans with a total of 8 epochs.

GELINE datasets (9 trials of 4 rats) were acquired with a 6-mm diameter home-made transceiver surface coil in anesthetized rats for evoked fMRI. GELINE was applied by using two saturation slices to avoid aliasing artifacts in the reduced field-of-view along the phase encoding (i.e., from left to right) direction (Fig. 1B and 1C). 2D line profiles were acquired to evaluate saturation RF pulses performance (Fig. 1D). The details of the saturation RF pulse were as follows: pulse shape sech.exc (adiabatic) installed on Bruker PV 5.1, length 1 ms, bandwidth 20250 Hz, FA 90°, bandwidth factor of the pulse 20250 Hz·ms, normalized shape integral 0.106428, rephasing factor 50 %, and derived power 0.2048 W. Laminar fMRI responses were acquired along the frequency-encoding direction (Fig. 1I and 1J). The following acquisition parameters were used: TR/TE 100/12.5 ms, TA 10 min 40 s, FA 50°, slice thickness 1.2 mm, FOV 6.4 × 1.2 mm², 1D readout matrix 128 (for the 2D line profiles, FOV 6.4 × 12.8 mm², matrix 64 × 128), and readout bandwidth 9014 Hz. The fMRI design paradigm for each epoch consisted of 1 s pre-stimulation, 4 s stimulation, and 15 s post-stimulation within a total of 20 s. A total of 6400 lines (i.e., 10 min 40 s) in each cortex were acquired every single trial in evoked fMRI. Evoked BOLD activation was induced by performing electrical stimulation to the left forepaw (300 µs duration at 2.5 mA repeated at 3 Hz for 4 s). A GELINE 2D image was additionally acquired with and without outer volume suppression (Figs. S2 and S3A). For this 2D image, 14T/13 cm (Magnex Scientific, horizontal bore) and a 6 cm diameter gradient system (100 G/cm, 150 μs rising time) was additionally used. A home-made RF transceiver surface coil with a 25 mm diameter was used on the rat brain. The following acquisition parameters were used: for without outer volume suppression, TR/TE 100/12.5 ms, TA 2 min 33 s, Average 4, FA 30°, slice thickness 1.2 mm, FOV 19.2 × 19.2 mm², and matrix 384 × 384, and readout bandwidth 17857 Hz; for with outer volume suppression, TR/TE 100/12.5 ms, TA 2 min 33 s, Average 4, FA 30°, slice thickness 1.2 mm, FOV 25.6 × 25.6 mm², and matrix 128 × 128, and readout bandwidth 9091 Hz.

SELINE datasets (18 trials of 4 rats) were acquired in anesthetized rats for evoked fMRI. SELINE was applied by the 180˚ RF pulse oriented perpendicular to the α˚ excitation RF pulse as moving the refocusing gradient to phase encoding gradient in order to obtain high spatial resolution without reduced FOV aliasing problem along the phase encoding (i.e., from left to right) direction (Fig. 1E and 1F). 2D line profiles were also acquired to evaluate the refocusing RF pulses performance (Fig. 1G). Laminar fMRI responses were acquired along the frequency-encoding direction (Fig. 1K and 1L). The following acquisition parameters were used: for 1000 ms SELINE acquisition, TR/TE/FA 1000/20 ms/90°, TA 10 min 40 s, slice thickness 1.2 mm, FOV 3.2 × 1.2 mm², and 1D readout matrix 64 (for the 2D line profiles, FOV 6.4 × 12.8 mm², matrix 64 × 128), and readout bandwidth 5000 Hz; for 200 ms SELINE acquisition, TR/TE/FA 200/10 ms/100° or 130° or 150°, TA 10 min 40 s, slice thickness 1.2 mm, FOV 6.4 × 1.2 mm², 1D readout matrix 64, and readout bandwidth 9014 Hz. A SELINE 2D image was also acquired with and without inner volume suppression (Figs. S2 and S3B) on the 14T MRI (Magnex Scientific) used to acquire the 2D GELINE (Figs. S2 and S3A). The following acquisition parameters were used: for without outer volume suppression, TR/TE 200/10 ms, TA 2 min 33 s, Average 4, FA 150°, slice thickness 1.2 mm, FOV 19.2 × 19.2 mm², and matrix 192 × 192, readout bandwidth 26455 Hz; for with outer volume suppression; TR/TE 200/10 ms, TA 2 min 33 s, Average 4, FA 150°, slice thickness 1.2 mm, FOV 12.8 × 12.8 mm², and matrix 128 × 128, and readout bandwidth 10000 Hz. The fMRI experiment set-up was identical to those of the GELINE in evoked fMRI.

All signal processing and analyses were implemented in MATLAB software (Mathworks, Natick, MA) and Analysis of Functional NeuroImages software (Cox, 1996) (AFNI, NIH, USA). For evoked fMRI analysis for Figure 1A, the hemodynamic response function (HRF) used was the default of the block function of the linear program 3dDeconvolve in AFNI. BLOCK (L, 1) computes a convolution of a square wave of duration L and makes a peak amplitude of block response = 1, with $g(t)=t4e−t/[44e−4]$⁠. Each beta weight represents the peak height of the corresponding BLOCK curve for that class. The HRF model was defined as follows:

Cortical surfaces were determined based on signal intensities of fMRI line profiles as described in the previous work. The detailed processing was conducted as provided in the previous line-scanning studies (S. Choi, Chen, et al., 2023; S. Choi, Zeng, et al., 2022; X. Yu et al., 2014). For quantitative comparison of background signals between GELINE and SELINE (Fig. 1D and 1G), the background signals were calculated as follows: N_bkg/S_roi × 100 %, where N_bkg was the mean of the outside-of-ROI signals and S_roi was the mean of the ROI signals at the cortical regions. For Figure 1I and 1K, demeaned fMRI time courses were used as follows: (x - μ), where x was the original fMRI time courses, and μ was the mean of the time courses. The line profile map concatenated with the multiple fMRI signals was normalized by a maximum intensity (Fig. 1I, 1K, 3A, and 3C; bottom). For Figure 3F, laminar-specific BOLD signals were normalized by a maximum intensity (i.e., a maximum mean value plus its standard deviation) for both GELINE and SELINE. The Z-score normalized time courses were calculated as follows: (x - μ)/σ, where x was original fMRI time courses and μ, σ were the mean and the standard deviation of the time courses, respectively (zscore function in MATLAB). Average BOLD time series and percentage changes were defined as (S-S0)/S0 × 100 %, where S was the BOLD signal and S0 was the baseline. S0 was obtained by averaging the fluctuation signal in the 1-s pre-stimulation window in evoked fMRI that was repeated every 20 s with the whole time series (640 s). The BOLD time series in each ROI were detrended (“polyfit” function in Matlab, order: 3) and bandpass filtered (0.01–0.1 Hz, FIR filter, order: 4096). The bandpass filtering was performed as a zero-phase filter by “fir1” and “filter” functions in Matlab, compensating a group delay (“grpdelay” and “circshift” functions in Matlab) introduced by the FIR filter. Temporal signal-to-noise ratio (tSNR) values were calculated across the cortical depths to compare tSNR differences between GELINE and SELINE. Note that σ was calculated as the standard deviation of the whole time series. We did not use the standard deviation of the time courses from repeated baseline periods because, given our evoked fMRI design paradigm, the baseline period may not be long enough to be considered as resting state (4 s stimulation vs. 16 s) (Chen & Pike, 2009; Mandeville et al., 1999) and thus could be influenced by the post-stimulus undershoot fluctuations. This calculation for the standard deviation may alter the Z-score normalized time courses and tSNR values. Student t-test was performed with the tSNR values of GELINE and SELINE (Fig. 1H). The p-values <0.05 were considered statistically significant.

To optimize the α˚ (alpha) excitation flip angle with short TR (i.e., 200 ms) in SELINE, signal intensities were calculated as a function of an excitation flip angle by solving the Bloch equation (Blenman et al., 2006; Diiokio et al., 1995), by employing the refocusing 180˚ RF pulse. The maximum signal intensity occurred at the optimal angles which was defined as follows:

where $α,β$ indicate excitation and refocusing flip angles, respectively, and T1, T2 indicate longitudinal and transverse magnetization parameters, respectively. For GELINE, a steady-state signal was calculated as a function of an excitation flip angle by solving the Bloch equation (Ernst & Anderson, 1966) and defined as follows:

where $θ$ indicates an excitation flip angle. T1 and T2 values (i.e., 2211 and 24 ms of somatosensory cortex at 16.4T, respectively) were estimated from the previous study (Pohmann et al., 2011) because the relaxation values should not be much different from those at 14T.

We developed the SELINE method to map laminar-specific BOLD responses across cortical layers at the primary forepaw somatosensory cortex (FP-S1) of anesthetized rats, which could be compared with the conventional GELINE method (X. Yu et al., 2014). First, unilateral electrical stimulation of the left forepaw of rats showed robust BOLD responses in the right FP-S1 using EPI-fMRI method (Grandjean et al., 2023) (Fig. 1A). Using the GELINE method, the selected FOV was defined by two saturation slices to avoid the aliasing problem along the phase encoding direction (Fig. 1B). In contrast, the same FOV could be selected by applying a refocusing 180° RF pulse perpendicular to the excitation slice with SELINE (Fig. 1E). To compare ROI selectivity between GELINE and SELINE, 2D in-plane images were acquired by turning on a phase encoding gradient (Fig. 1C and 1F) and 1D profiles were plotted by averaging all readout voxels of the 2D image (Fig. 1D and 1G). Full width at half maximum (FWHM) of the 1D profiles was estimated: For GELINE, trial #1) and #2) ~1.8 mm, and for SELINE: trial #1) and #2) ~1.5 mm. Background signals were estimated from the areas outside of the FOV (for details, see the Method section): For GELINE, trial #1) 10.6 %, #2) 15.8 %, and for SELINE: trial #1) 4.9 %, #2) 5.0 %. The variation in the background signal suppression of GELINE was potentially caused by imperfect saturation RF pulses and B0 field inhomogeneity in the outside of the ROI (Fig. S2). This result indicated the efficiency of the SELINE method to produce sharper 2D slice profiles and lower background signals.

To study the laminar fMRI characteristics of GELINE and SELINE across the cortical layers, we calculated tSNR with 1D line-profiles which were acquired by turning off the phase encoding gradient. The tSNR of SELINE was higher than those of GELINE (Fig. 1H). The tSNR graph of SELINE showed a gradually decreasing trend across the cortical depth while those of GELINE showed a gradually increasing trend. To predict tSNR difference by solving the Bloch equation (see the Methods), we calculated theoretical tSNR with steady-state signals given T1 (~2200 ms), T2 (~24 ms), and T2* (~20 ms) values and scan parameters such as TRs (1000 ms vs. 100 ms), TEs (20 ms vs. 12.5 ms), flip angles (90° vs. 50°), readout bandwidths (5000 Hz vs. 9014 Hz), and readout FOVs (3.2 mm vs. 6.4 mm), assuming tSNR increased linearly with SNR (Han et al., 2019) and noise contribution to BOLD signals was identical in both acquisitions. The ratio of the simulated tSNR of SELINE vs. GELINE was ~2.2 (~7.0 x 10^-6/3.2 x 10^-6) while the ratio of the experimental tSNR was ~1.6 (~50.8/31.0). In this simulation, we did not consider different effects of noise sources (e.g., thermal noises, physiological noises, and background noises from outside of ROI) (Khatamian et al., 2016; Krüger & Glover, 2001; Ragot & Chen, 2019) which was also likely to lead to the difference not only between the simulated and experimental results but also between GELINE and SELINE. It should be noted that the transceiver surface coil caused a non-uniform B1 field and possibly contributed to the tSNR difference between GELINE and SELINE across cortical depths.

As shown in Figure 1I-L, we demonstrated dynamic BOLD responses across different cortical layers of FP-S1 from a representative trial in individual GELINE (Fig. 1I and 1J) and SELINE (Fig. 1K and 1L) studies. Figure 1I demonstrated periodic evoked BOLD signals upon left forepaw electrical stimulation with the T2*-weighted GELINE method, showing the dynamic laminar-specific BOLD responses as a function of time peaked around the superficial layer in the FP-S1 (4 s on/16 s off for each 20 s epoch, total 32 epochs). Average BOLD time series and laminar-specific BOLD maps illustrated that the peak BOLD response is located at L1, highlighting large draining vein effects at the cortical surface (J. Goense et al., 2016; Goense & Logothetis, 2006; Han et al., 2019, 2021; Polimeni et al., 2010; Zhao et al., 2004, 2006) (Fig. 1J). In comparison to GELINE, SELINE also detected robust FP-S1 BOLD signals across different cortical layers (Fig. 1K), but showed the peak BOLD signal located at L4, presenting improved spatial specificity to deeper cortical layers (J. Goense et al., 2016; Goense & Logothetis, 2006; Han et al., 2019, 2021; Harel et al., 2006; Zhao et al., 2004, 2006) and similar BOLD signal change with a longer TR (1000 ms) (Fig. 1L).

We further investigated the reproducibility of laminar-specific peak BOLD responses, as well as the variability of laminar-specific BOLD response patterns, between the two methods (14 trials from 3 animals). The GELINE method detected peak BOLD signals primarily located at L1, but the peak BOLD signal detected by the SELINE method was much deeper (Fig. 2). In animal #3, the strong BOLD signal detected in the superficial voxel indicates a large draining vein dominating the voxel BOLD signal (Fig. 2G, 2H, and 2I). A similar BOLD response was also detected by the SELINE method, which is most likely to stem from non-negligible intravascular effects of the large draining vein in the superficial voxels with only 50 μm thickness or increased motional narrowing effects due to the fast diffusion rate of spins from cerebrospinal fluid (CSF) (Pfaffenrot et al., 2021; Uludaǧ et al., 2009).

Interestingly, the layer-specific BOLD signal varied largely across animals in both GELINE and SELINE maps. Besides the primary BOLD peak in L1 of GELINE, a second peak appeared in L4 in some animal (Fig. 2D and 2F). And for the SELINE method, the primary peak also varied at L4 and L2/3 (Fig. 2B, 2C, 2E, and 2F), which presented highly different laminar patterns from GELINE even when acquired from the same animal with interleaved trials during experiments (Fig. 2A and 2D). These results have suggested that the profile of laminar-specific BOLD signals can vary largely across animals, which may present varied dynamic patterns of BOLD responses due to the altered neurovascular coupling across different cortical layers.

We performed BOLD fMRI experiments with a 200 ms TR by applying optimized flip angles based on the Bloch equation (Blenman et al., 2006; Diiokio et al., 1995) (see the Method section). For comparison, we also performed the GELINE method in the same anesthetized rat. As shown in Figure 3A-D, we demonstrated the evoked BOLD responses across the cortical layers upon the periodic electrical stimulation with the GELINE (Fig. 3A and 3B) and SELINE (Fig. 3C and 3D) methods, showing the average BOLD time series and percentage changes peaked at L1 in both GELINE and SELINE. To characterize the laminar-specific BOLD responses, the normalized BOLD signals were plotted across the cortical layers. As shown in Figure 3F, the GELINE method had a steep signal drop from L1 to L2/3, while the SELINE method had a gradual signal drop across the cortical depth. Although the laminar-specific BOLD responses of the 200 ms SELINE maintained large vessels sensitivity in the superficial layers (Fig. 3D), the peak BOLD signal of the SELINE were much lower than that of the GELINE: ~10 % vs. ~30 % and the slopes of the normalized BOLD signal plot illustrated that the SELINE method had less bias to large draining veins than the GELINE method at the superficial layers (slope at L1 and L2/3: SELINE; -0.31 vs. GELINE; -0.49). Taken together, these results indicate that the high temporal SELINE method reduces the large vessel contribution to the BOLD responses by minimizing magnetic susceptibility effects at the superficial layers (i.e., L1 and L2/3).

To select an optimized flip angle, the tSNR of different flip angles was plotted (Fig. 3G). Even though the optimal flip angle for TR 200 ms was ~150° and had the highest tSNR, the difference of the tSNR change was relatively small in multiple trials with the different flip angles regardless of B1^- inhomogeneity correction (Fig. 3G and Fig. S1) (Delgado et al., 2020). This result was possibly caused by the long T1 effect (~2200 ms) of the cortex in the SELINE acquisition with a short TR (200 ms) (Pohmann et al., 2011). Same as the theoretical predictions based on the Bloch equation (Blenman et al., 2006; Diiokio et al., 1995), $e−TR/T1$ was almost close to 1 and thus, the maximum intensity at the optimal flip angle did not change a lot. Average tSNR values for GELINE and SELINE were 15.3 and 27.0 while the tSNR efficiency of those was 48.4 and 60.4, respectively. It is noteworthy that the average tSNR of SELINE was higher at the superficial and middle layers than that of GELINE (Fig. 3E and 3G) due to a larger flip angle (100-150° vs. 50°) and longer TR (200 ms vs. 100 ms). In summary, these results not only demonstrated less magnetic susceptibility effects at the superficial layer, but also highlighted laminar specificity enhancement in SELINE with high temporal resolution.

In this study, we applied the SELINE method to investigate laminar-specific evoked BOLD responses across cortical layers with high spatial and temporal resolution. The SELINE method has sharper and better ROI-selectivity than the GELINE method, employing the refocusing 180° RF pulse perpendicular to the excitation plane. It should be noted that a part of tissue on the upper side of the 2D SELINE image (Fig. 1F) remained outside the ROI and was severely distorted due to nonlinear gradient around the air-tissue boundary (i.e., skull) and chemical shift of fat signals (Sakurai et al., 1992) at an ultra-high magnetic field (14.1T). Since our analyses focused on cortical layers (0–2 mm) and the varying displacement resulted in geometric distortion occurred outside of the cortex, the influence of the imperfect suppression was negligible. As ascertained in the 2D GELINE image (Fig. 1C), a major signal source of the pile-up displacement artifact was fat tissue. This issue thus can be alleviated by applying fat suppression pulses and high readout bandwidth (Raimondo, Heij, et al., 2023; Sakurai et al., 1992). Our results show that the peak signal of SELINE is spread across the cortical layers while that of GELINE is at the superficial layer (Han et al., 2019; Krüger & Glover, 2001). By pushing the temporal resolution of SELINE to 200 ms, we also demonstrate the feasibility to map laminar-specific BOLD responses with less large draining vein effects (Boxerman et al., 1995; J. Goense et al., 2016; Goense & Logothetis, 2006; Han et al., 2019, 2021; Weisskoff et al., 1994; Yacoub et al., 2003, 2005; Zhao et al., 2004, 2006), in comparison to the GELINE method.

Significant effort with high-field fMRI has been made to explore laminar fMRI responses corresponding to distinct information flows (e.g., top-down/bottom-up or feedforward/feedback) at high spatial and temporal scales in both animals and humans. Among these efforts detecting BOLD, cerebral blood volume (CBV), and cerebral blood flow (CBF) signals with both SE and GRE methods, cortical depth-dependent fMRI has identified hemodynamic regulation, blood volume distribution, circuit-specific laminar responses, and hierarchical information streams across cortical layers in animal (Albers et al., 2018; S. Choi, Chen, et al., 2023; S. Choi, Zeng, et al., 2022; J. B. Goense & Logothetis, 2006; Jung et al., 2021; Lu et al., 2004; Shen & Duong, 2016; Silva & Koretsky, 2002; X. Yu et al., 2014; Zhao et al., 2004, 2006) and human brains (Finn et al., 2019; Huber et al., 2017; Kashyap et al., 2018; Sharoh et al., 2019; Y. Yu et al., 2019). In particular, high-resolution CBV-fMRI, based on the VASO mapping scheme, has been used to measure layer-specific directional functional connectivity across human motor cortex and somatosensory and premotor regions (Huber et al., 2017). It should be noted that the cortical thickness of human brains is in the range of 1–4 mm, which is highly comparable to that of rodent brains in the range of 1–2 mm (Fischl & Dale, 2000). Given the limited spatial resolution of the high field laminar-fMRI method (~600-700 um), the thoroughly counted voxels across different cortical depth regions are in the single digit number. This could be much better improved by the developed line-scanning fMRI method, as well as with ultra-fast sampling rates that enable the detection of fast hemodynamic responses across cortical layers without compounding artifacts (Caballero-Gaudes & Reynolds, 2017).

Recently, the GRE-based line-scanning BOLD mapping scheme has been implemented to investigate BOLD signals across cortical layers in human fMRI studies (Morgan et al., 2020; Raimondo et al., 2021; Raimondo, Priovoulos, et al., 2023). Nevertheless, since SAR is proportional to the square of the magnetic field (B0) and the duty cycle of the sequence, high temporal resolution in ultra-high-field fMRI studies can be constrained by SAR limits especially representing a safety-related limit in high-field human MRI system (e.g., 7T and 9.4T). Compared to the SELINE method, the required saturation RF pulses of the GELINE method theoretically result in higher SAR and total RF power limits with short TRs, inducing more complicated aliasing problems. For the SELINE method, the beam-like line-scan projection has been previously applied for probing myeloarchitecture across cortical layers in the primary somatosensory cortex (S1) and the primary motor cortex (M1) of the human brain (Balasubramanian et al., 2021) and mapping irreversible and reversible transverse relaxation rates (i.e., R2 and R2´) in primary visual cortex (V1), S1, and M1 of human brains (Balasubramanian et al., 2022). The sharper line-profile has been also demonstrated in human fMRI studies by employing the SELINE method at a cost of compromising tSNR and BOLD sensitivity (Raimondo, Heij, et al., 2023). Although pre- and post-data processing steps (i.e., NORDIC (Vizioli et al., 2021), SNR-optimized coil combination (S. Choi et al., 2016; Raimondo et al., 2021; Roemer et al., 1990), and independent component analysis (Mckeown et al., 2003)) were applied to enhance tSNR and functional sensitivity, no task-driven activation was observed in this initial human fMRI study (Raimondo, Heij, et al., 2023). This suggests that prospective motion correction (Bause et al., 2020), averaging of more runs (Huettel & Mccarthy, 2001), and localized surface coils (Kashyap et al., 2018), which were not incorporated into this study, are required for the future work. In contrast, the previous large-tip-angle spin-echo line-scanning fMRI study in anesthetized rats elucidated that early onset of laminar-specific BOLD responses occurred at the middle layers for comparison of diffusion fMRI onset while potential aliasing issues of saturation RF pulses remained (Nunes et al., 2021). We thus applied this SELINE method to better characterize layer-specific fMRI features across cortical depths at FP-S1 of rodent brains without the need for additional saturation RF pulses. The SELINE method employed the spin-echo scheme to reduce the large draining vein effect, which could be further distinguished from the deeper cortical layer responses given the high spatial resolution (Fig. 1I–L).

As reported in previous studies (Duong et al., 2003; J. Goense et al., 2016; J. B. Goense & Logothetis, 2006; Yacoub et al., 2003, 2005; Zhao et al., 2004, 2006), GELINE is more sensitive to large veins at the pial surface but has poor specificity across different cortical depths, whereas SELINE is less vulnerable to superficial large draining veins but has good sensitivity to micro-vessel across cortical layers. However, the largely varied laminar patterns of the BOLD responses were observed in both methods (Fig. 2). This may suggest that the varied patterns of laminar-specific BOLD signals pertain to microvascular biases and baseline blood volume distribution across cortical layers (Hartung et al., 2022a, 2022b). In addition, TE is crucial for optimizing the SELINE BOLD responses because the relative micro- and macro-vascular contributions to BOLD signals can be changed by altering the TE (Boxerman et al., 1995; Uludaǧ et al., 2009). Here, we chose a TE of 20 ms (for a TR of 1000 ms) which was analogous to the tissue T2 value (~20 ms) to minimize macrovascular contribution to the BOLD responses (Pohmann et al., 2011). For high temporal SELINE acquisition (TR 200 ms), a TE of 10 ms was chosen to optimize macrovessel suppression and tSNR preservation. Also, the TEs were relatively long compared to the venous blood T2 value at 14T (<10 ms) (Lin et al., 2012), resulting in small intravascular contribution to the SELINE BOLD responses. Note that, nevertheless, intravascular SE-based BOLD signals do not completely disappear even at high-field MRI when the partial volume contribution of vessels to the given voxel is not negligible (Uludaǧ et al., 2009). In some cases of SELINE, strong BOLD responses at the superficial layer existed (Fig. 2H and 3D). As accepting that both intra- and extra-vascular signals contribute to the total BOLD responses, it is most likely to originate from the following factors: 1) The extravascular effect in combination with non-negligible cerebral blood volume effect from large vessels can be dominant in the GRE-based acquisition scheme. Depending on the surface draining vein and diving artery localization, the superficial voxels can be heavily influenced by the extravascular effects, as well as active or passive vessel dilation effect. The combined signal changes could dramatically influence the signal changes in the GRE scheme. 2) The non-negligible intravascular effects can remain in the SE-based acquisition scheme for superficial voxels with large partial volume contribution from vessels (Boxerman et al., 1995). In these superficial voxels, the vessel volume contribution is far higher than 2-4% as reported for conventional fMRI studies (Ji et al., 2021; Kim & Ogawa, 2012; Weber et al., 2008), and the direct intravascular effects caused by oxy- and deoxy-hemoglobin ratio changes could contribute to the BOLD signal from these voxels. In parallel, the SE scheme would reduce the extravascular effect, which further makes the intravascular effect from vessels less negligible. 3) Varied inflow effects (Axel, 1984) can influence laminar fMRI signals given the vascular distribution. The blood flow would alter the spin-echo effect given the TE selected for 200 ms and 1000 ms TRs. The spins in the flowing blood would experience less rephasing effect based on the SE scheme and show weaker signals. Nevertheless, if spins excited by a 90º (or αº) RF pulse, originally outside of the line profile, flow into the 180º RF excitation slice, it would contribute to the line-profile signals erroneously and impact the laminar BOLD responses. For GRE-acquisition scheme, if flowing fresh blood enters into the image plane, new spins would experience an excited RF pulse. With increasing blood flow velocity, the fresh blood signal increases since the total number of experienced RF pulses decrease. 4) CSF partial volume effect (Pfaffenrot et al., 2021; Uludaǧ et al., 2009) might occur at superficial voxels due to slow CSF flow at a TR of 1000 ms. This is similar to the inflow effects discussed in the previous section. Meanwhile, for the superficial voxels when CSF took a non-negligible partial volume contribution similar to blood vessels, it could also contribute to the altered laminar BOLD responses. Since the spatial resolution of the fast SELINE method was lower than that of the GELINE method, the layer boundaries of the laminar response patterns might be vague and blurred in the fast SELINE results (Fig. 3). Whereas the varied peak profiles of BOLD responses exist across different cortical layers presumably due to the confounding factors, these results illustrate the feasibility of the line-scanning method to detect distinct laminar BOLD responses. It provides a high-resolution mapping scheme when investigating altered neurovascular coupling events and functional connectivity across cortical layers.

The main limitation of SELINE is the slow sampling rate. The sampling rate (TR) is determined by two times of the TE. TE is generally limited by the duration of the excitation and refocusing RFs, the duration of the slice rephasing gradient, and the duration of the frequency encoding gradient (i.e., readout matrix size and bandwidth). We attempted to shorten the TR by adjusting the excitation flip angle (α). Based on the Bloch equation (Blenman et al., 2006; Diiokio et al., 1995), we have estimated the appropriate angles with a short TR (i.e., 200 ms). Our results show the feasibility of the fast SELINE method which has a good sampling capacity capturing dynamic BOLD signals from superficial to deeper layers. For future work, the fast SELINE method should be optimized in terms of spoiling and phase cycling schemes to enhance tSNR. Furthermore, simultaneous GRE- and SE-type fMRI acquisitions can be applied to better characterize laminar-specific fMRI patterns and minimize time dependency of dynamic fMRI responses by employing GRASE (Oshio & Feinberg, 1991)-based line-scanning in rodents as already suggested for the human fMRI mapping (S. Choi, Yu, et al., 2022). For laminar human fMRI studies, a fat suppression RF (e.g., SPIR, SPAIR, STIR, etc.) should be applied to avoid fat aliasing artifacts in cortical areas at relatively low magnetic fields (e.g., 7T) (Raimondo et al., 2021; Raimondo, Heij, et al., 2023). Imperfect fat suppression can also be alleviated by adjusting acquisition parameters (e.g., readout matrix size and bandwidth).

All other data generated during this study are available from the corresponding author upon reasonable request. The related image processing codes are available from the corresponding author upon reasonable request.

X.Y. designed and supervised the research; X.Y. and K.S. conceptualized the research; L.G., W.M., and X.Y. performed animal experiments; S.C. and X.Y. acquired data; S.C. analyzed data; D.H., L.G., W.M., R.P., N.A., and X.Y. provided technical support; S.C. and X.Y. drafted the manuscript; and S.C., D.H., R.P., L.G., K.S., and X.Y. reviewed and edited the manuscript.

The authors declare no competing interests.

This research was funded by NIH funding (RF1NS113278, R01NS124778, R01NS122904, R01NS120594, R21NS121642), NSF grant 2123971, and the S10 instrument grant (S10 MH124733–01) to Martino’s Center, German Research Foundation (DFG) Yu215/2-1, 3-1, DGF SPP2041 SCHE658/17 and SCHE658/12, BMBF 01GQ1702, and the internal funding from Max Planck Society. We thank Dr. J. Engelmann, and Ms. H. Schulz for technical support, Dr. P. Douay, Ms. R. König, and Ms. M. Pitscheider for animal support, the AFNI team for the software support.

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