Revealing layer-specific cortical activity in human M1 using high-resolution line-scanning fMRI

Ultra-high-field functional MRI (fMRI) in humans has enabled the study of cortical activity at very high spatiotemporal resolutions (Guidi et al., 2016; Huber et al., 2015; Jia et al., 2023; Shao et al., 2021; Shen et al., 2008; Siero et al., 2013). These advances have enabled researchers to perform functional MRI of cortical laminae, which can disentangle local circuits (Finn et al., 2019; Goense, 2018; Huber, Finn, et al., 2021) with implications for studying brain-wide function. Unlike electrophysiology, which is considered the gold standard for measuring neural activity (Self et al., 2019), fMRI is completely non-invasive, making it broadly accessible for studying brain function in humans and animals. As it is based on a hemodynamic response, fMRI is an indirect measure of cortical activity.

The most commonly used sequence in fMRI is the gradient-echo-BOLD (GE-BOLD) sequence due to its ease of application, its strong signal changes, and its sampling efficiency. However, GE-BOLD is sensitive to various physiological parameters such as local changes in cerebral blood flow (CBF), cerebral blood volume (CBV), and oxygen consumption (CMRO2). GE-BOLD shows a bias toward the draining vasculature, which limits its precision in mapping the exact location of cortical activity. Non-BOLD sequences have gained traction because they measure specific aspects of neural activity such as CBF (Kashyap et al., 2021) or CMRO2 (Shao et al., 2021), but due to their complex and often custom implementation, the most commonly used contrast is still the GE-BOLD contrast (Bergmann et al., 2024; Fracasso et al., 2018; Huber et al., 2019; Nothnagel et al., 2023)

The pursuit of high spatial resolution usually requires compromising temporal resolution and/or spatial coverage. To achieve high spatial resolution (voxel size <0.8 mm), the field-of-view (FOV)/imaging slab is reduced to contain the area of interest (10–20 cm), and repetition times (TR) are often long, with reported TRs of up to 9 s (Koiso et al., 2023), although the current standard is 2–3 s.

Simultaneously achieving high spatial and temporal resolution can be achieved via a compromise of reducing coverage. By minimizing both spatial (voxel size <0.8 mm) and temporal (TR <1 s) resolution, dimensions tangential to the cortical surface can be collapsed into a single line. Yu et al. (2014) implemented a line-scanning fMRI method with very high temporal (TR = 50 ms) and spatial (voxel size 50 μm) resolutions and used a sensory stimulation paradigm to elicit functional responses in rat sensory and motor cortex. They detected laminar differences in the BOLD time courses, with the hemodynamic response being initiated in layer IV of sensory cortex and layers II/III and V of motor cortex in rats. Further improvements and cross-validations followed the original investigation. Albers et al. (2018) compared sensory and optogenetically induced BOLD responses using line-scanning fMRI in rats and confirmed that the method can detect time differences in the onsets in the range of 100 ms. Nunes et al. (2019, 2021) introduced diffusion-weighted line-scanning fMRI and demonstrated its ability to detect early onsets (<100 ms) of the functional signal in the middle layers of rat sensory cortex. Choi, Zeng, et al. (2022) and Choi et al. (2023) presented a method for bilateral line scanning and detected laminar patterns of resting-state signals. Recent work has also demonstrated the possibility of functional line scanning using GRASE (Choi, Yu, et al., 2022) and spin-echo-BOLD readouts (Choi et al., 2024).

Line-scanning fMRI has recently been applied to humans (Balasubramanian et al., 2021; Morgan et al., 2020; Raimondo et al., 2021, 2022; Raimondo, Priovoulos, et al., 2023). One of the major challenges of applying line scanning to humans is staying within SAR safety limits while maintaining a well-defined saturation profile at a short TR. Raimondo et al. (2021) reported their initial results and measured GE-BOLD signals with high spatial (200 μm) and temporal (TR = 200 ms) resolution using two saturation bands. However, this approach produced a line profile with a full-width-half-maximum (FWHM) of ~6 mm, which limits the accuracy of laminar BOLD studies when the spatial extent of the activation is smaller than this distance, as signals from adjacent areas would be projected into the line. Heij et al. (2023) highlighted the importance of planning the line location in advance, as the highly convoluted human cortex makes manual line placement less accurate, particularly when aiming for the line to pass exactly radially through the cortical surface. Balasubramanian et al. (2021) presented a line-scanning method for acquiring layer-specific diffusion imaging and demonstrated that they can measure diffusion differences in the cortical layers of M1 and S1 with a 250–500 μm radial resolution.

Here, we demonstrate GE line scanning for layer-specific human brain imaging with a very sharp and thin line profile (3 mm). This makes the method suitable for application to smaller (tangential) activations, such as digit maps in S1 or the hand knob in M1. Line scanning is well suited for cortical depth-specific functional imaging due to its very high resolution in the laminar dimension while maintaining high signal-to-noise ratio (SNR). At a laminar resolution of 0.39 mm, the line scan voxel of 0.39 x 3.0 x 3.0 mm³ captures the same number of protons as an equivalent isotropic voxel of (1.5 mm)³. The line scan voxel integrates signals from the same cortical depth, therefore, reducing partial volume effects and increasing temporal SNR (tSNR) (Blazejewska et al., 2019; Kiebel et al., 2000).

Our GE-line-scanning method involves four saturation regions to create a sharp line profile at acceptable SAR levels. We apply the method to obtain laminar positive and negative BOLD time courses in the primary motor cortex during a finger-tapping task. This stimulus was chosen as it produces a broad activity in the hand knob in M1 and is a common task for human laminar fMRI studies in M1 (Chai et al., 2020; Guidi et al., 2016; Huber et al., 2017; Persichetti et al., 2020; Shao et al., 2022). We find that the time courses of both positive and negative BOLD show laminar specificity. Additionally, we describe a method for accurately segmenting the gray matter from line-scanning acquisitions using both the magnitude and phase information from accompanying high-resolution FLASH scans. We also detail the process of preparing and executing our GE line-scanning sequence and the post-processing, with the overall goal of further advancing the human line-scanning fMRI method.

We implemented a multi-echo GE-based line-scanning sequence and demonstrated its use in non-invasively recording layer-specific functional activity in primary motor cortex (M1). We collapse the entire MRI acquisition into a single one-dimensional readout line through the postcentral gyrus (Fig. 1C). We use a dedicated alignment procedure (see sections 2.3 and 2.4) to ensure that this line passes through the cortex of M1 at a specific, targeted location (Fig. 1C-F). Figure 1B illustrates the amplitude of the collapsed signal and Figure 1C-F illustrates the laminar arrangement of M1 and S1, providing a visual representation of this process.

The line-scanning sequence was based on a multi-echo FLASH (fast low-angle shot) pulse sequence on a 7T Siemens Magnetom Terra (Siemens Healthineers, Germany), incorporating saturation pulses to select the line (see Fig. 2). The signal saturation was achieved with four custom Shinnar-Le Roux (SLR) saturation pulses. Two inner SLR pulses were created using the sigpy package (Ong & Lustig, 2019) with N = 256 samples; passband ripple 0.01%; stopband ripple 0.1%; maximum phase; flip angle 110˚; bandwidth-time product 19; 7800 µs pulse duration; thickness 15 mm. Two outer SLR pulses were created also using sigpy with N = 256 samples; passband ripple 0.01%; stopband ripple 0.1%; maximum phase; flip angle 120˚; bandwidth-time product 2.7; 3800 µs pulse duration; thickness 120 mm. We applied root flipping to lower the B1 peak power of all SLR pulses. We positioned the inner SLR pulse symmetrically around the center to create a well-defined line profile in the center of the slice. We positioned the outer SLR pulse on the rest of the brain tissue to suppress the remaining signal. The minimum TR was 109 ms but we used 250 ms to reduce SAR and increase signal-to-noise (SNR). We expected this TR to capture the essential properties of the dynamics of the hemodynamic response.

We compared this implementation on a regular FLASH image with a conventional saturation strategy by placing vendor-implemented saturation bands to create the line, with otherwise matching parameters. For all approaches, we measured SAR levels as taken from the log files of the scanner. We used MATLAB to calculate the FWHM of each approach taking the average of 10 center profiles on the FLASH images recorded with saturation pulses (see Fig. 5B). We calculated a leakage as percentage as the fraction of signal generated from voxels outside the line, divided by the entire signal of all voxels.

In the functional experiment, our proposed sequence allows for the option to turn the phase-encoding gradient off. The sequence offers several choices:

1. A regular FLASH (Fig. 2D) to provide an anatomical reference that is used for segmenting the gray matter in the area of interest (see section 2.6.5).

2. A FLASH with saturation bands. This shows the target line of 3 mm in the center of the slice (Fig. 2D, third panel), and was used to assess signal suppression quality and to confirm the absence of artifacts that would be projected into the line.

3. Line scanning with saturation pulses on and phase-encoding switched off (Fig. 2E). This turns the sequence into a one-dimensional recording at the center of the slice. We used this option for the functional runs.

The sequence parameters of these scans can be found in Table 1.

We recruited 14 healthy volunteers for this study. We excluded two participants because we could not find a sufficiently flat stretch in their motor cortex to position the line (see section 2.3) and one participant because of head motion. We conducted functional line scanning on 12 healthy volunteers (7/5 m/f, age range 21–35 years, 10 right handed, 1 left handed, 1 ambidextrous). Before their participation, we obtained informed consent from each volunteer. The study was reviewed and approved by the local ethics committee of the College of Medical, Veterinary and Life Sciences at the University of Glasgow (Project no. 302122). We applied standard site-specific inclusion criteria, including the absence of metal in the body, age over 18 years, normal or corrected-to-normal vision, normal hearing, and no diagnosis of neurological or mental health disorders. These criteria were confirmed during the participants’ initial visit to the facility.

Before the 7T session, we identified a flat area in the hand region of the primary motor cortex (M1) on either the left or the right hemisphere (see Fig. 3A and 3B) at 3T using data acquired earlier at 3T (Tim Trio MRI System, Siemens Healthineers, Erlangen, Germany). We acquired an anatomical scan (MPRAGE, see Table 1) to create a cortical reconstruction using Freesurfer 7.3 (recon-all) (Fischl, 2012). Figure 3 shows how we use Freesurfer’s probability atlas to narrow down our search for the anterior primary motor cortex to Brodmann Area 4a (BA4a) on the white matter surface (Huber et al., 2017). We then refined the search within inside BA4a for a sufficiently flat spot following the method described in Morgan et al. (2020). Figure 3C indicates the curvature of the cortex with highlighted area over BA4a. Green areas indicate low curvature (suitable for line scanning) and red and blue areas indicate high curvature (unsuitable for line scanning).

After identifying a suitable location in BA4a, we used Morgan et al.’s (2020) method to calculate coordinates for a slice that contains the line in its center. Given that we could rotate the slice containing the line by 360º around the cortical surface normal, we selected an angle that led to a slice with the least amount of signal outside the line in the image. This procedure led consistently to a double-oblique coronal slice, tilted 20º–40º toward the axial plane and approximately 2º–10º away from the sagittal plane. If no flat area could be found, we did not proceed with the volunteer (we were not able to find a flat spot in two participants based on this criterion).

For the functional line-scanning experiments, we used a 7T Magnetom Terra scanner (Siemens Healthineers, Erlangen, Germany) equipped with a 32-channel head coil (Nova Technologies). After obtaining an anatomical scan (MP2RAGE, see Table 1) at the beginning of the session, we exported the DICOM images to an external computer. We then converted the images to NifTI format using dicom2niix (Li et al., 2016) and applied LAYNIIs LN_MPRAGE_DNOISE tool to denoise the MP2RAGE background (Huber, Poser, et al., 2021). Next, we computed a rigid 6-degree-of-freedom transformation in ITK-SNAP (Yushkevich et al., 2006) to find a transformation between the head position during the pre-session anatomical scan and the current head position in the scanner (see Table 1). We then applied the alignment procedure of Morgan et al. (2020) to obtain the coordinates of the new slice for the current head position (see section 2.3). We then typed in the new slice coordinates as three positions for the location, two angles for the double-oblique slice orientation, and the in-plane slice angle. We confirmed the new slice location using the procedure outlined in Figure 4.

To confirm that the line passes through M1 radially, we acquired a low-resolution 2D FLASH at the new slice location (Fig. 4C). This acquisition had an additional slice excitation pulse that was rotated by 90º around the center line to cast a shadow at the location of the target line (see Table 1). This created a visual guide to confirm that the center line passes radially through M1, and allowed for small adjustments if needed. If the slice did not seem aligned, we would repeat the MP2RAGE and repeat the process. If it still did not work, we would end the session. If we were at the right location, we would then run a single-slice functional EPI scan at the target slice (see Table 1) with a finger-tapping paradigm of the contralateral hand (described in section 2.6) and use an online GLM to confirm functional activity at the target location (Fig. 4B).

We performed an automated B0-based shimming procedure for a region that contains the slice and a minimum thickness of 90 mm (Nassirpour et al., 2018). Next, we would acquire a B1 map and adjust the transmit voltage using a region of interest (ROI) in the motor area. For anatomical reference, we would then acquire a high-resolution 2D FLASH scan at this location (see section 2.3; parameters in Table 1), with and without saturation pulses (Fig. 2D). We then started the functional experiment of four functional line scans (two contralateral tapping and two ipsilateral tapping, with contralateral first).

We used a finger-tapping task using a block design, which elicits a strong BOLD response in M1 and is a standard paradigm for motor tasks (Chai et al., 2020; Guidi et al., 2016; Huber et al., 2017, 2019; Persichetti et al., 2020; Shao et al., 2022). Each block consisted of 20 s of tapping followed by 20 s of rest, starting with a baseline. There were eight trials in total for each run. The participants could see instructions on the screen. During the baseline period, a cross was displayed, and participants were told to relax. During the tapping period, a message instructed them to tap with either their left or right hand. Participants were asked to tap only with their left hand or only with their right hand for the entire run. If the participant could remain still, we recorded at least four runs: two with the contralateral hand, then two with the ipsilateral hand. If there was time left at the end of the session, we would acquire an additional run for contralateral tapping (this was the case for three participants). Before the functional MRI scan, we practised the tapping motion outside of the scanner with the participants. We instructed the participants to tap four fingers simultaneously against the thumb, at a self-chosen pace.

The raw data for all custom line-scanning acquisitions were exported from the scanner using the vendor-specific tool (twix). Subsequently, the data were transferred to an external computer and then converted from the vendor-specific .dat format into regular data arrays using the pymapVBVD tool, a python port of Philips Ehses’ MATLAB tool mapVBVD (https://github.com/wtclarke/pymapvbvd).

The structural 2D data were transformed to the image domain using a 2D FFT along the read and phase dimensions. Then, a coil sensitivity map was calculated for the slice where the line was located using the 2D FLASH scan without signal saturation, employing BART’s ecalib tool (Uecker et al., 2015). Subsequently, the complex coil images of both structural scans (saturated and non-saturated) were combined using a SENSE-1 coil combination. The phase data were unwarped using the skikit package. Both the magnitude and phase data were stored for further analysis to create an anatomical reference (see section 2.6.5).

The 1D functional data were transformed to the image domain using a 1D FFT along the readout dimension. Then, a root-of-sum-of-squares (RSS) coil combination was applied to reduce the array along the channel dimension. The data were then combined along the echo dimension using an optimal weighting suggested by Kundu et al. (2012) with a target T₂^* of 25 ms using the formula

The final data formed a two-dimensional array, with the first dimension covering all voxels along the line and the second dimension covering all TRs in the functional run (Fig. 2E).

We aligned all functional runs individually to the collapsed structural FLASH (with saturation) to ensure consistent activation location across all runs. We chose not to use sub-voxel shifts to avoid potential interpolation-induced blurring. Instead, we calculated the number of voxels to shift each run based on a simple numerical approach to find the highest cross-correlation to the collapsed scan. In our numerical search, we used 30 voxels surrounding the target gray matter in the M1–S1 region and then shifted each functional run by this number of voxels. This ensured that all the gray matter voxels in the functional runs aligned with the collapsed profile where the segmentation was defined (see section 2.6.5).

We used the output of the alignment to detect scans that showed excessive head motion. Since a line scan only provides information about one of the six possible degrees of freedom for motion (the readout dimension), we could not apply retrospective motion correction. Any head motion after the alignment process meant that the signal was recorded from a cortical area different from the intended one. Additionally, the steady state was temporarily affected in the event of head motion which influenced the functional signal until it recovered. Therefore, we chose motion detection over motion correction. Our rule was: If we detected a shift during the alignment of the brightest voxel in the first echo inside the brain of more than three voxels over the entire run, we would exclude the entire run.

Across the 12 participants, we had to reject 12 out of 27 runs for contralateral (44%) and 16 out of 22 runs for ipsilateral (72%) tapping based on the above defined criterion. Most of the rejected scans were at the end of a session. The final subject pool consisted of eight participants with good data from both contralateral and ipsilateral tapping, three participants with good data from contralateral tapping only, and one participant with no acceptable data.

After the alignment, all functional data were stored in NifTI format using the nibabel package (Brett et al., 2023). The phase image of the profile scan was also stored separately for segmentation later on (see section 2.6.5). The slice orientation information was obtained for all scans from the original DICOM after conversion to NifTI using the dicom2niix package (https://github.com/rordenlab/dcm2niix).

First, we started by drawing an initial segmentation of the gray matter on the magnitude image of the profile scan using ITK-SNAP and saved the result in NIfTI format. Next, we collapsed the center line of this segmentation into a 1D line, which should closely label the target gray matter in the line scan data. However, in practice, it may not be perfect. To find a consistent upper border between gray matter and cerebrospinal fluid (CSF), we applied the following method illustrated in Figure 5. We searched for the brightest voxels near one of the edges of the collapsed segmentation and marked this edge as the border between CSF and gray matter. We then used the unwarped (Van Der Walt et al., 2014) phase image of the 2D FLASH (without saturation) to determine the boundary between gray and white matter. Since the phase is sensitive to T2*, we anticipated a slight phase jump between white and gray matter (Fig. 5B). Subsequently, we counted the number of voxels between the CSF and phase jump and marked them as pure gray matter voxels across cortical depth. If a clear phase jump was not found at the white matter border, we used the Freesurfer reconstruction from the MP2RAGE anatomical scan to estimate the gray matter thickness.

In order to improve signal-to-noise and to remove nuisance signals, we applied multi-echo independent component analysis (ME-ICA) denoising. We used the TEDANA v24.0.2 (DuPre et al., 2021) toolbox voxels to remove non-BOLD signals from the data. We provided a binary mask that contains all brain, and used tedana_workflow with default options except -fastica 15 and masktype “decay.”

We used a small temporal filter to reduce thermal noise by employing a running average filter using the LN_TEMPSMOOTH from the LAYNII package (using -box option with a kernel size of 7 TRs) (Huber, Poser, et al., 2021). We then normalized the time courses to percent signal change.

We assigned the voxels at different relative cortical depth to cortical layers based on the histology of human M1 (Palomero-Gallagher & Zilles, 2019). This allowed us to create individual time courses for the cortical surface, LII/III, LV, and LVI for each subject (Fig. 4A and 4D). We used linear interpolation for the gray matter defined in section 2.6.5 to resample the time courses of each subject into 24 voxels before assigning them to a cortical layer.

To characterize the shape of the time courses, we conducted a simple numerical search using MATLAB’s find function to identify the time point at which the percent signal change reached its maximum amplitude at each cortical depth for each subject. We applied a strong temporal filter (gamma = 3 s) on the raw time courses to reduce the effect of high-frequency noise, and used MATLAB’s find function to identify the first index where the signal was above 95%, 50%, and 10% of the signal strength at the time point of the maximum amplitude.

We also fitted a polynomial curve to the raw line-scan data to find laminar onset times of the hemodynamic response (HRF) using MATLAB’s curvefit toolbox. We followed the procedure used in Albers et al. (2018) using an HRF model from Boynton et al. (1996):

We used the following constraints for the fit: 0 <A< 40; 0.1 < T₀< 3.5; 2 <a< 4; 0.5 <b< 3. We used unfiltered, trial-averaged (8 trials), and subject-averaged (n = 10) raw time courses, selecting the interval from the onset to the offset of the tapping stimulus (20 s) to get the best estimation for the onset time T₀ at each cortical depth.

To evaluate the performance of the SLR pulses for selecting the line, we compared it with conventional saturation strategies. Figure 6 shows the signal saturation and line sharpness for several implementations of the line-scanning sequence. The SAR levels were (from left to right) 2%, 14%, 14%, 46%. The FWHM were (from left to right) 10.9 mm, 7.4 mm, 3.9 mm. To achieve the linewidths of 10.9 and 7.4 mm using the conventional saturation approach, we needed to nominally overlap (-4 and -8 mm, respectively) the saturation slabs with respect to the targeted line width (3 mm). This makes it problematic to define signal leakage with respect to the target for the conventional approach. Therefore, we calculated the leakage as the area under the curve using the voxels outside the FWHM, divided by the total area using all voxels. The leakage values for this comparison were (from left to right) 24.2%, 18.2%, 2.1%. If we reference to the target linewidth, leakage is considerably higher, and close to 100% for the conventional saturation pulses. For the conventional approach, further reduction of the linewidth would lead to substantial reductions in signal amplitude within the line.

Figure 7 shows a typical line-scan experiment conducted during finger tapping in a representative participant. The line passes through the M1/S1 region associated with the hand representation. After identifying the gray matter area (see section 2.6.5), we located nine gray matter voxels. Figure 7D-G shows the signal from gray matter and a few extra voxels at the CSF and white matter border. As expected, the strongest amplitude of the BOLD responses was found in the upper layers at the CSF border, as indicated by the time course and beta values for this region.

The mean cortical thickness of our patch in M1 was 3.3 ± 0.2 mm and spanned on average 8.5 ± 0.4 voxels (n = 10). Figure 8 shows the mean time courses from our line-scanning experiment across the cortical depth, comparing tapping with the contralateral and ipsilateral hand. We found that contralateral tapping evoked the strongest BOLD response just above the pial surface (Fig. 8A). We observed BOLD responses with smaller amplitude in the gray matter. The time courses in the gray matter region near the cortical surface (around L2/3) did not reach plateau but rose slowly reaching 95% of their maximum amplitude toward the end of the stimulus period at 22.0 ± 1.9 s after onset (n = 10) (Fig. 9). In contrast, deeper layers plateaued and reached 95% of the maximum amplitude at 9.8 ± 2.4 s after stimulus onset. During ipsilateral tapping, we observed a negative BOLD response that peaked in the middle of the cortex (Fig. 8B) while the surface was inactive or had a slight positive BOLD. We also noted that the deeper layers reached 95% of the maximum negative BOLD earlier (11.2 + 1.5 s) than the gray matter closer to the cortical surface (15.3 ± 1.5 s). The negative BOLD returned to baseline faster after stimulus offset than the positive BOLD in all cortical layers.

Figure 9A shows the trial-averaged time courses of each cortical depth for both ipsilateral and contralateral tapping. Figure 9B shows the mean response and the results from the fitting process are described in section 2.6.9. Figure 9C summarizes the results of the analysis of temporal properties for T10, T50, and T90 when the signal reached 10%, 50%, and 90% of the peak value in each cortical depth.

The results of the fitting process for the onset times are displayed in Table 2.

We present a line-scanning fMRI method for obtaining highly localized laminar BOLD responses in a selected location of the human cortex. Our implementation has low SAR demands and achieves a sharp line profile, which was used to perform layer-specific fMRI in humans with high spatiotemporal resolution. By placing the line on the primary motor cortex (M1) during a finger-tapping task, we demonstrate its ability to sample distinct BOLD time courses across cortical depth. This makes the method a promising complementary tool for neuroscientific and clinical applications of high-resolution fMRI.

Several studies have found differences in the time courses between layers and paradigms, although mostly in primary sensory cortex. Jin and Kim (2008) reported laminar differences in onset times and time-to-peak of BOLD and CBV-weighted recordings. Tian et al. (2010) showed that the hemodynamic response has the earliest onset in layer IV, suggesting it begins in layer IV and propagates to the superficial layers. Yacoub et. al. (2006) also showed different time courses in superficial vessels and gray matter in the cat, while Bohraus et al. (2023) found similar laminar differences in the time courses between BOLD, CBF, and CMRO2 in macaques. Insights from animal studies using electrophysiological and oxygen measurements (Bentley et al., 2016) suggest that the shape of BOLD time courses may be modulated by various underlying mechanisms, which may also differ between cortical areas and layers. Shmuel et al. (2006) and He et al. (2022) studied HRFs of negative BOLD signals and found that they are significantly different from positive BOLD signals, suggesting that they could be caused by different underlying neural or vascular dynamics. de la Rosa et al. (2021) also noted that the HRF of negative BOLD significantly differs from the HRF of positive BOLD, and is not just its inverted version. Goense et al. (2012) reported that negative BOLD in macaque V1 has different laminar profiles, suggesting that different neurovascular coupling mechanisms operate at different cortical depths. Kashyap et al. (2018) used localized fMRI at a very high spatial resolution in humans and found the BOLD response within gray matter had a stronger post-stimulus undershoot than closer to the pial surface. In summary, differences in the time courses between layers and between paradigms are common. Although most of the studies focussed on primary sensory or visual cortex, we recorded time courses in M1 and also found different laminar time courses between positive and negative BOLD. Our method provides a tool to study these differences in humans with high laminar and tangential precision.

Line-scanning fMRI of cortical layers was first applied in animals (Yu et al., 2014). They made use of its high spatiotemporal resolution and measured laminar specificity of the BOLD response in the rat cortex (Choi et al., 2023, 2024). In human studies, line-scanning fMRI is a relatively new frontier (Morgan et al., 2020; Raimondo et al., 2021). The SAR limitations in humans are more stringent than those in animals, which makes it challenging to achieve GE-based line scanning with saturation RF pulses that have both a high bandwidth-time product and a very short TR. Raimondo et al. (2021) introduced a line-scanning GE sequence using a standard saturation scheme, with a line with FWHM of ~6 mm, which works in V1 or where activation is extensive, but may cause blurring by including BOLD and nuisance signals from non-target areas when the tangential spread of the activation is small. When we reduced the linewidth to ~4 mm using the conventional saturation pulses, there was a high signal loss in our target line, indicating the challenge of GE-based line scanning of small target areas or where sharp lines are needed. Spin-echo line scanning dramatically improves the sharpness of the line (Raimondo, Heij, et al., 2023), but at the cost of SNR and temporal resolution. They demonstrated that SE line scanning can create sharp line profiles, making this a promising approach with a clear line selection.

The approach we present here involves slightly longer TR due to the extra saturation pulses, but achieves a very sharp line profile (FWHM of 3.9 mm, see section 3.1) that maintains signal within the line (97% of signal coming from the line, see section 3.1) and reduces the contribution of signals from outside the line (2.2% leakage signal from outside the line). Fat signals did not tend to leak into our line in M1, but that could be different for other brain regions. Although there is no risk that fat produces functional signals, it could introduce signal instabilities. We recommend to apply fat saturation in those cases even though it may increase TR. The smallest possible TR with our GE-based sequence was 109 ms. This TR leads, however, to high SAR and reduced SNR. We, therefore, chose a slightly longer TR than other studies (Raimondo et al., 2021). Sequence improvements such as ramp sampling, stronger gradients or lower gradient spoil moments, and/or RF pulse schemes with reduced SAR demand could potentially lower the TR in the future. Regarding the spatial resolution limit, the sequence is capable of voxel sizes down to around 0.1 mm, but such voxel sizes may be more susceptible to head motion.

We demonstrate the potential of GE-line scanning by using it to extract laminar fMRI signals in the hand area of M1 while participants were performing a motor task. Our method successfully detected positive and negative BOLD time courses, and revealed differences across cortical depth. We chose this stimulus as it is widely used and known for creating strong BOLD amplitudes, and because we expected to see differences in the time courses (Schäfer et al., 2012). For the positive BOLD response, we observed a peak at the cortical surface and a faint “shoulder” in the deeper layers, which is in agreement with previous high spatial resolution studies at lower temporal resolution (Huber et al., 2017; Knudsen et al., 2023; Shao et al., 2021). For the negative BOLD response, we see a trend of a negative “double peak” with the signal peaking inside the cortex twice, with a clear negative peak in the upper layers (around L1/2, 0–20% cortical depth) and a weaker peak in the deeper layers (L6, 67–100%). This could correspond to neural inhibition (Stefanovic et al., 2004) with the upper “input” layers and the deeper “output” layers of M1 receiving transcallosal inhibitory input (Persichetti et al., 2020). There is also a hint of a post-stimulus laminar pattern for the negative BOLD which could relate to Mullinger et al.’s (2017) findings that suggest underlying inhibitory neural process. We think this could be explored in a future study using a stimulus paradigm tailored to elicit post-stimulus responses, perhaps in combination with a CBF-weighted reference. The pattern we observe may be in line with findings from Huber et al. (2017), considering that their cortical activation shows large variations across M1 and not all areas in BA4a seem to have a “double-peak” pattern. However, other studies report a negative BOLD response during ipsilateral tapping (Mullinger et al., 2014; Newton et al., 2005; Stefanovic et al., 2004; Yuan et al., 2011), although at lower spatial resolution. We speculate that the difference could be explained by some variability of the experimental design, assumptions of the HRF shape using a GLM, age dependency, or by individual differences. Further studies, especially a CBV-weighted or CBF-weighted line-scanning implementation such as VASO (Huber et al., 2017) or VAPER (Chai et al., 2020) may confirm this. We also note a high percent signal change with a few participants toward the cortical surface. We explain this by proximity to superficial pial veins or larger draining veins in these participants, (see Bause et al., 2020). We did not specifically select areas without large draining veins for our line positioning, and included all voxels in our analysis. In future studies, we aim to choose alternative locations during line positioning further away from large veins to reduce this bias.

The main limitation for GE line scanning is head motion (Raimondo, Priovoulos, et al., 2023), and our method is affected by this problem too. When participants move their heads even slightly (>0.3 mm) after the alignment process, the assignment of a time course to a cortical depth at a specific location is lost and cannot be recovered from the data. We, therefore, relied on recruiting experienced subjects, strong head fixation with cushions, training of participants to keep their head in the same position, and running the experiment as quickly as possible. This left us a window of maximum 20 min to collect the functional data while more standardized high-resolution scannings involve typically >1 h of scanning. We note that this could be seen also as an advantage as subjects can keep their attention high for the task. In addition, the demand for storing and processing the data is negligible in comparison with submillimeter 3D scanning with millions of voxels. We had to exclude data due to motion, especially in the last 10 min where we collected data for the negative BOLD. Raimondo, Priovoulos, et al. (2023) suggested prospective motion correction to tackle this problem, but the authors reported unwanted T1-transient signals due to their navigators. Another limitation of our method is that it reaches the limits of SAR for human scanners at very short TRs. Our current setup reaches the SAR limit around TR = 200 ms depending on the participant’s body height/weight. Future optimizations would need to find ways to reduce TR, which would be beneficial for experiments that study laminar onset times as they would be in the expected range of 100–150 ms across the layers based on animal research (Yu et al., 2014). We also noted that not all of our participants were right handed, presenting a potential cofound for the interpretation of our data.

Another improvement would be a more optimizing saturation strategies as suggested by Cloos et al. (2024) to achieve shorter TRs and improved SNR, or to explore different contrasts (spin echo, MT weighted, SSFP fMRI, diffusion fMRI, VASO/VAPER). Other optimizations will focus on ways to reduce motion problems by head fixation strategies, external motion cameras, or more optimal sequence-based navigators to decide which data or trials to reject due to motion. In addition, line-scanning data could benefit from model-based physiological noise removal techniques such as Agrawal (2020) which may reduce pulsation-related artifacts and further improve signal denoising.

We demonstrate the implementation of a GE-based line-scanning method for human fMRI with a narrow FWHM (3.9 mm) at high spatiotemporal resolution (voxel size 0.39 x 3.0 x 3.0 mm, TR = 250 ms). Our findings demonstrate the suitability of the method for laminar fMRI by extracting time courses across the cortical depth of the hand knob in the primary motor cortex during a finger-tapping task. The results reveal distinct temporal properties across cortical depths and paradigms. The ability to detect subtle differences in the hemodynamic response across cortical depths opens up new possibilities for non-invasively studying neural dynamics in humans.

We will make code and data available upon reasonable request.

N.N., A.T.M., and J.G. conceptualized and designed the study. N.N. wrote the MRI sequence and collected the data. A.T.M. created the software for aligning the line location to the head position. N.N., A.T.M., and J.G. performed data analysis including creating data analysis tools. N.N. and J.G. drafted the original manuscript. N.N., A.T.M., L.M., and J.G. reviewed and edited the manuscript. L.M. and J.G. provided resources and supervised the project.

This work was supported by the Medical Research Council (MR/R005745/1 awarded to J.G.; MR/N008537/1, awarded to L.M. and J.G.) and by the European Union’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement Nos. 785907 and 945539 (Human Brain Project SGA2 and SGA3) awarded to L.M.

We have no competing interests to declare.

We would like to thank Frances Crabbe for helping with the acquisition of anatomicals at the 3T. During preparation of the draft, N.N. used Grammarly to correct grammatical errors and spelling mistakes. After using these tools, we reviewed and edited the content as needed and we take full responsibility for the content of the publication.