Revealing layer-specific cortical activity in human M1 using high-resolution line-scanning fMRI
Nils Nothnagel, A. Tyler Morgan, Lars Muckli, Jozien Goense
TL;DR
This study uses high-resolution line-scanning fMRI to reveal layer-specific activity in the human motor cortex, showing differences in blood flow responses across cortical layers.
Contribution
The paper introduces a gradient-echo-based line-scanning fMRI method with four saturation regions to achieve high spatiotemporal resolution in human cortical laminae.
Findings
Deeper cortical layers showed faster BOLD response onset and return to baseline compared to superficial layers.
Line-scanning fMRI can capture laminar differences in hemodynamic responses during motor tasks.
The method achieves a narrow FWHM and high resolution while maintaining SAR safety limits.
Abstract
In recent years, ultra-high field functional MRI has allowed researchers to study cortical activity at high spatiotemporal resolution. Advancements in technology have made it possible to perform fMRI of cortical laminae, which is crucial for understanding and mapping of local circuits and overall brain function. Unlike invasive electrophysiology, fMRI provides a non-invasive approach to studying human and animal brain function. However, achieving high spatial resolution has often meant sacrificing temporal resolution. In contrast, line-scanning fMRI maintains both high spatial and temporal resolution, and has been successfully applied to animals to detect laminar differences of the hemodynamic response. Although this method has been extended to human brain imaging in initial studies, staying within SAR safety limits while maintaining a well-defined saturation profile at a short TR is a…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8
Fig. 9| 3T Anatomy | 7T Anatomy | Online GLM | Confirm location | 2D Anat. reference | Functional line scan | |
|---|---|---|---|---|---|---|
| Sequence | MPRAGE | MP2RAGE | EPI | FLASH | FLASH | FLASH |
| Voxel size(mm 3 ) | 1.0 x 1.0 x 1.0 | 0.63 x 0.63 x 0.63 | 1.0 x 1.0 x 3.0 | 1.2 x 1.2 x 3.0 | 0.39 x 0.39 x 3.0 | 0.39 x 3.0 x 3.0 |
| FOV (mm 3 ) read-phase-slice | 256 x 256 x 192 | 240 x 225 x 160 | 150 x 150 x 3 | 224 x 224 x 3 | 224 x 224 x 3 | 224 x 3 x 3 |
| TR (ms) | 2600 | 4680 | 2000 | 250 | 250 | 250 |
| TE (ms) | 40 | 2.09 | 25 | 12 | 14,25,36,47 | 14,25,36,47 |
| FA | 8° | 5°/6° | 68° | 15° | 15° | 15° |
| IPat | 2 | 3 | 3 | 1 | 1 | 1 |
| Has phase encoding? | Yes | Yes | Yes | yes | yes | no |
|
0-20%
|
20-49%
|
49-67%
|
67-100%
| |
| T 0 |
1.1 s
|
0.7 s
|
1.0 s
|
2.3 s
|
- —Medical Research Councilhttp://dx.doi.org/10.13039/501100009187
- —Medical Research Councilhttp://dx.doi.org/10.13039/501100009187
- —Horizon 2020http://dx.doi.org/10.13039/501100007601
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsFunctional Brain Connectivity Studies · Advanced MRI Techniques and Applications · Advanced Neuroimaging Techniques and Applications
Introduction
1
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)andChoi 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^3^captures the same number of protons as an equivalent isotropic voxel of (1.5 mm)^3^. 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.
Methods
2
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 (seesections 2.3and2.4) to ensure that this line passes through the cortex of M1 at a specific, targeted location (Fig. 1C-F).Figure 1Billustrates the amplitude of the collapsed signal andFigure 1C-Fillustrates the laminar arrangement of M1 and S1, providing a visual representation of this process.
Panels (A, C, D) show the location of sensorimotor cortex and panel (B) displays the 1D profile of a line that passes through the layers of the human primary motor cortex (M1) and the primary sensory cortex (S1). Panel (C) shows the location and orientation of a line passing through the hand knob of M1. Panels (E and F) display the underlying microscopy image of the M1/S1 region. The definition of cortical layers in M1 is marked in colors and is based onPalomero-Gallagher and Zilles (2019), showing a thick layer 6, a thinner layer 5, an absence of layer 4, and a thick layer 2/3. Panel (E) shows how the line from Panel B passes radially through the layers of M1. The underlying microscopy image and the atlas are taken from the Allen 34 years old Human Brain Atlas (Allen Institute for Brain Science,https://atlas.brain-map.org/andhttps://atlas.brain-map.org/atlas?atlas=265297126,Ding et al., 2016;Hawrylycz et al., 2012). Panel (F) indicates the cortical depth estimation for the image of Panel (E).
Implementation of the line-scanning sequence
2.1
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 (seeFig. 2). The signal saturation was achieved with four custom Shinnar-Le Roux (SLR) saturation pulses. Two inner SLR pulses were created using thesigpypackage (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 usingsigpywith 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.
The line-scanning sequence. Panel (A) shows the location of a double-oblique slice containing the line. It includes the hand knob area of M1, S1, with the center line of the slice passing perpendicularly through M1 at the target location. Panel (B) illustrates the arrangement of the four saturation pulses (colors) to suppress all signals in the slice except for the center line. Two outer saturation pulses (purple and yellow) suppress the signal from the bulk, while two narrow saturation regions with a sharp edge (blue and pink) create a sharp gap of 3 mm in the center of the slice for the line. Panel (C) shows the pulse diagram with the timing of the RF pulses and gradients, and the readout gradients for an acquisition with four echoes. Panel (D) consists of four images: the first image shows the acquired slice with the line indicated in the center (x-axis: phase encoding direction, y-axis: readout direction). The second image indicates the location of M1 and S1, and the position of the four saturation regions (colored regions). The third image displays the image with the saturation pulses turned on, demonstrating that all the signal is suppressed except for the targeted center line of the slice. The last image overlays the third image (saturation on) onto the first image (saturation off). Panel (E) shows the functional time courses of a 4-min line scan of all voxels in the center line (x-axis: time domain, y-axis: readout direction).
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 (seeFig. 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:
A regular FLASH (Fig. 2D) to provide an anatomical reference that is used for segmenting the gray matter in the area of interest (seesection 2.6.5).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.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 inTable 1.
Volunteer recruitment
2.2
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 (seesection 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.
Pre-session slice planning
2.3
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 (seeFig. 3Aand3B) at 3T using data acquired earlier at 3T (Tim Trio MRI System, Siemens Healthineers, Erlangen, Germany). We acquired an anatomical scan (MPRAGE, seeTable 1) to create a cortical reconstruction using Freesurfer 7.3 (recon-all) (Fischl, 2012).Figure 3shows 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 inMorgan et al. (2020).Figure 3Cindicates 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).
Illustration of the procedure for identification of a suitable location for line scanning in primary motor cortex. Panel (A) indicates the hand area in the anterior primary motor cortex (Brodmann Area 4a, BA4a, pink) on the cortical white matter surface. Panel (B) shows the location of M1/S1 and the target location on the MPRAGE anatomical scan. Panel (C) shows the same cortical surface as in Panel A but with the curvature value for each vertex. We chose as suitable location an area with low curvature (green) within the hand area of BA4a. Panel (D) zooms in on the structural scan of Panel (B) and marks the selected vertex in M1 in voxel space.
After identifying a suitable location in BA4a, we usedMorgan 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).
Setup at the scanner
2.4
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, seeTable 1) at the beginning of the session, we exported the DICOM images to an external computer. We then converted the images to NifTI format usingdicom2niix(Li et al., 2016) and applied LAYNIIsLN_MPRAGE_DNOISEtool 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 (seeTable 1). We then applied the alignment procedure ofMorgan et al. (2020)to obtain the coordinates of the new slice for the current head position (seesection 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 inFigure 4.
Illustration of the procedure for planning the target line for scanning. Panel (A) shows a high-resolution FLASH scan with the line passing through the center from left to right. The outer regions are shaded out since they will be saturated in the line scan. The top row in Panel (B) shows the result of an online GLM of a low-resolution EPI (1.0 x 1.0 x 3.0 mm3) to confirm BOLD activation in the patch. Panels (C and D) zoom in on the target location of a low-resolution (1.2 x 1.2 x 3.0 mm3) (C) and high-resolution FLASH (0.39 x 0.39 x 3.0 mm3) (D). The red line shows the WM surface segmentation (red) around the target spot. This is to visualize how the line passes perpendicularly through the cortex, as indicated by the red line showing the GM/WM boundary. While it appears on the low-resolution FLASH as if the line does not pass perpendicular WM at the identified line-scan position, the smaller voxels and better contrast of the high-resolution FLASH reveal a small, 3–4 mm straight patch between a convex and a concave region that we targeted in our selection process (seesection 2.3).
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 (seeTable 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 (seeTable 1) with a finger-tapping paradigm of the contralateral hand (described insection 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 (seesection 2.3; parameters inTable 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).
Finger-tapping paradigm
2.5
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.
Data processing
2.6
Exporting and reading raw data
2.6.1
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 thepymapVBVDtool, a python port of Philips Ehses’ MATLAB toolmapVBVD(https://github.com/wtclarke/pymapvbvd).
Offline data reconstruction
2.6.2
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’secalibtool (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 theskikitpackage. Both the magnitude and phase data were stored for further analysis to create an anatomical reference (seesection 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 byKundu et al. (2012)with a target T_2_^*^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).
Alignment of functional runs
2.6.3
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 (seesection 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.
Storing data in NifTI format
2.6.4
After the alignment, all functional data were stored in NifTI format using thenibabelpackage (Brett et al., 2023). The phase image of the profile scan was also stored separately for segmentation later on (seesection 2.6.5). The slice orientation information was obtained for all scans from the original DICOM after conversion to NifTI using thedicom2niixpackage (https://github.com/rordenlab/dcm2niix).
Location of the postcentral sulcus in the line and GM segmentation
2.6.5
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 inFigure 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.
shows the definition of the gray matter boundaries in the magnitude (A, C) and phase (B, D) images. Panel (A) shows the magnitude profile of the line passing through M1/S1, with the gray matter voxels marked. The brightest voxel in this profile indicates the CSF border. Panel (C) shows the magnitude of the high-resolution 2D FLASH image. The line is marked in the center. Panel (B) shows the phase of the line, with the gray matter voxels marked. The white matter border is indicated by the location where the phase slightly jumps at the gray/white matter border due to a T2 change, and we used this jump to define the white matter border. Panel (D) shows the high-resolution 2D FLASH phase image.*
Multi-echo ICA denoising
2.6.6
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 usedtedana_workflowwith default options except*-fastica 15andmasktype “decay.”*
Temporal filtering and converting signal to percent change
2.6.7
We used a small temporal filter to reduce thermal noise by employing a running average filter using theLN_TEMPSMOOTHfrom 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.
Definition of laminar profiles
2.6.8
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. 4Aand4D). We used linear interpolation for the gray matter defined insection 2.6.5to resample the time courses of each subject into 24 voxels before assigning them to a cortical layer.
Characterizing time courses and determining onset times
2.6.9
To characterize the shape of the time courses, we conducted a simple numerical search using MATLAB’sfindfunction 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’sfindfunction 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’scurvefittoolbox. We followed the procedure used inAlbers et al. (2018)using an HRF model fromBoynton et al. (1996):
We used the following constraints for the fit: 0 <A< 40; 0.1 < T_0_< 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_0_at each cortical depth.
Results
3
Signal saturation
3.1
To evaluate the performance of the SLR pulses for selecting the line, we compared it with conventional saturation strategies.Figure 6shows 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.
FLASH sequence using different signal saturation strategies. Panel (A): MRI image of the same slice using the following methods: no saturation, standard saturation regions (nominal gap -4 mm), standard saturation regions (nominal gap -8 mm), SLR pulses. Panel (B): Line profile across the center for each image indicating the width of the line. Spatial resolution: 0.39 x 0.39 x 3.0 mm3, TE/TR: 25 ms/250 ms
Representative example of line-scanning data in M1 during finger tapping
3.2
Figure 7shows 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 (seesection 2.6.5), we located nine gray matter voxels.Figure 7D-Gshows 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.
Laminar time courses in a representative participant. Panel (A) displays the magnitude image of the 2D FLASH scan with the saturation pulses turned off. The slice contains the line passing through the center, with a small ROI marking the targeted cortical patch in M1. Panel (B) shows a zoomed-in version of this line where the saturation pulses were turned on. All signals outside the line drop to zero, and the line passes perpendicularly through the WM surface which is indicated as an outline. In Panel (C), the image with saturation turned on (in color) is overlaid on the image without saturation, with the targeted ROI marked in a small box. Panels (D and E) display the trial-averaged functional activation across the cortical layers in the marked region during contralateral (D) and ipsilateral (E) tapping. The y-axis indicates cortical depth (0% at cortical surface and 100% at the white matter border), and the x-axis indicates time from 0 to 40 s. The BOLD activity peaks at the pial surface and reduces in amplitude with increasing cortical depth. Panels (F and G) show the percent signal change of the laminar time courses in the targeted patch during six trials, sorted by cortical depth, together with the stimulus onsets of the experiment (Panel F: contralateral. Panel G: ipsilateral). Next to the plot, we also show the trial-averaged time courses of the data from the same panel. The gray bar above the plot shows the 20 s stimulus cue for the finger tapping. The upper layers show the strongest BOLD signal, with the signal of the deep layers dropping to 0% signal change. Panel (H) shows the betas for all voxels across cortical depth, showing that the betas also drop to zero with increasing cortical depth.
Cortical depth analysis
3.3
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 8shows 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.
Panels show the time courses of the BOLD response to contralateral (A) and ipsilateral (B) finger tapping, averaged across subjects and trials for the layers in M1. The x-axis ranges from 0 to 40 s and the y-axis shows the cortical depth from top to bottom. CSF is on top, and the white matter is on the bottom. Panels (C and D) show the laminar profile with the mean percent signal change across the cortical ribbon (n = 10, mean + sem), 10 s after stimulus onset. The negative BOLD appears inside the gray matter, while the positive BOLD appears mainly at the cortical surface.
Laminar time courses and fitting of onset times at different cortical depths in the primary motor cortex (0–20% depth: CSF, 20–49%: L2/3, 50–67% L5, 67–100%: L6). (A) Trial average response during contralateral (blue) and ipsilateral (orange) finger tapping. The thick line is the group average, and the thin lines are individual data. The stimulus interval is marked with a gray bar (20 s duration, onset and offset marked with dotted lines). (B) Result of the fitting procedure with mean raw data (black) and polynomial fit (red). (C) Timings derived from the average time courses across subjects during contralateral tapping. We could find significant differences (p < 0.05) for peak times between upper and deeper layers (L1/CSF-L6, p = 0.0003; L1/CSF-L5, p = 0.009; L2/3-L6, p = 0.002; L2/3-L5, p = 0.02) when the amplitude reached 50% and 90% (T50, T90) of the maximum amplitude. All p-values were Bonferri corrected for multiple comparisons.
Time course analysis
3.4
Figure 9Ashows the trial-averaged time courses of each cortical depth for both ipsilateral and contralateral tapping.Figure 9Bshows the mean response and the results from the fitting process are described insection 2.6.9.Figure 9Csummarizes 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 inTable 2.
Discussion
4
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, whileBohraus 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)andHe 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, seesection 3.1) that maintains signal within the line (97% of signal coming from the line, seesection 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 toMullinger 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 fromHuber 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, (seeBause 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 asAgrawal (2020)which may reduce pulsation-related artifacts and further improve signal denoising.
Conclusion
5
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.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Agrawal , U. , Brown , E. N. , & Lewis , L. D. ( 2020 ). Model-based physiological noise removal in fast f MRI . Neuro Image , 205 , 116231 . 10.1016/j.neuroimage.2006.01.015 31589991 PMC 6911832 · doi ↗ · pubmed ↗
- 2Albers , F. , Schmid , F. , Wachsmuth , L. , & Faber , C . ( 2018 ). Line scanning f MRI reveals earlier onset of optogenetically evoked BOLD response in rat somatosensory cortex as compared to sensory stimulation . Neuro Image , 164 , 144 – 154 . 10.1016/j.neuroimage.2016.12.059 28012967 · doi ↗ · pubmed ↗
- 3Balasubramanian , M. , Mulkern , R. V. , Neil , J. J. , Maier , S. E. , & Polimeni , J. R . ( 2021 ). Probing in vivo cortical myeloarchitecture in humans via line-scan diffusion acquisitions at 7 T with 250-500 micron radial resolution . Magnetic Resonance in Medicine , 85 ( 1 ), 390 – 403 . 10.1002/mrm.28419 32738088 PMC 7951328 · doi ↗ · pubmed ↗
- 4Bause , J. , Polimeni , J. R. , Stelzer , J. , In , M.-H. , Ehses , P. , Kraemer-Fernandez , P. , Aghaeifar , A. , Lacosse , E. , Pohmann , R. , & Scheffler , K . ( 2020 ). Impact of prospective motion correction, distortion correction methods and large vein bias on the spatial accuracy of cortical laminar f MRI at 9.4 Tesla . Neuro Image , 208 , 116434 . 10.1016/j.neuroimage.2019.116434 31812715 · doi ↗ · pubmed ↗
- 5Bentley , W. J. , Li , J. M. , Snyder , A. Z. , Raichle , M. E. , & Snyder , L. H . ( 2016 ). Oxygen level and LFP in task-positive and task-negative areas: Bridging bold f MRI and electrophysiology . Cerebral Cortex , 26 ( 1 ), 346 – 357 . 10.1093/cercor/bhu 260 25385710 PMC 4677981 · doi ↗ · pubmed ↗
- 6Bergmann , J. , Petro , L. S. , Abbatecola , C. , Li , M. S. , Morgan , A. T. , & Muckli , L . ( 2024 ). Cortical depth profiles in primary visual cortex for illusory and imaginary experiences . Nature Communications , 15 ( 1 ), 1002 . 10.1038/s 41467-024-45065-w PMC 1083744838307834 · doi ↗ · pubmed ↗
- 7Blazejewska , A. I. , Fischl , B. , Wald , L. L. , & Polimeni , J. R . ( 2019 ). Intracortical smoothing of small-voxel f MRI data can provide increased detection power without spatial resolution losses compared to conventional large-voxel f MRI data . Neuro Image , 189 , 601 – 614 . 10.1016/j.neuroimage.2019.01.054 30690157 PMC 6668026 · doi ↗ · pubmed ↗
- 8Bohraus , Y. , Merkle , H. , Logothetis , N. K. , & Goense , J. ( 2023 ). Laminar differences in functional oxygen metabolism in monkey visual cortex measured with calibrated f MRI . Cell Reports , 42 ( 11 ), 113341. 10.1016/j.celrep.2023.113341 37897728 · doi ↗ · pubmed ↗