Resolving laminar activation in human V1 using ultra-high spatial resolution fMRI at 7T

The mesoscopic organization of the human neocortex is of great interest for cognitive neuroscience. However, fMRI in humans typically maps the functional units of cognitive processing on a macroscopic level. With the advent of ultra-high field MRI (≥7T), it has become possible to acquire fMRI data with sub-millimetre resolution, enabling probing the laminar and columnar circuitry in humans. Currently, laminar BOLD responses are not directly observed but inferred via data analysis, due to coarse spatial resolution of fMRI (e.g. 0.7–0.8 mm isotropic) relative to the extent of histological laminae. In this study, we introduce a novel approach for mapping the cortical BOLD response at the spatial scale of cortical layers and columns at 7T (an unprecedented 0.1 mm, either in the laminar or columnar direction). We demonstrate experimentally and using simulations, the superiority of the novel approach compared to standard approaches for human laminar fMRI in terms of effective spatial resolution in either laminar or columnar direction. In addition, we provide evidence that the laminar BOLD signal profile is not homogeneous even over short patches of cortex. In summary, the proposed novel approach affords the ability to directly study the mesoscopic organization of the human cortex, thus, bridging the gap between human cognitive neuroscience and invasive animal studies.


(Supplementary
) 4 . A single session consisted of a total of three anatomical scans (MP2RAGE localizer, MP2RAGE reference, MI-EPI reference) and a total of five functional scans. The sequence parameters, listed in Supplementary Table 1, were optimized in pilot experiments.
In each session, a high-resolution MP2RAGE 5 (0.7 mm isotropic) was acquired first and used as an anatomical localizer for slice planning. The strong tissue contrast of the MP2RAGE was critical for positioning the imaging volume on the ROI for the functional scans. Three functional runs were acquired using an inhouse adapted FLASH 6 with anisotropic spatial resolutions (0.1x1.4x2.0 mm 3 ) (called Anisotropic Voxel FLASH, abbreviated as AVF) along the readout and phase-encoding (matrix size 1200x86) directions. Factor of two GRAPPA undersampling was used, resulting in a total of 43 kspace lines being acquired per slice at TR=41 ms (total slice TR=1.76 s) and a flip angle of 12°. No partial Fourier was used to minimize blurring. At the echo time (TE) of 26 ms, a readout bandwidth of 60 Hz per pixel could be afforded in order to maximize SNR. In addition, RF excitation with high bandwidthtime product (8.5) was applied to reduce slice profile imperfections. The high spatial resolution (0.1 mm) along the readout direction was selected so as to sample cortical depths, and therefore positioned orthogonal to the cortical surface in the ROI. The precise slice placement and orientation was thus adapted on a subject-by-subject basis, based on the curvature of the calcarine sulcus, such that the phase-encoding direction is locally along the cortical distance (see Figure 1). This is critical to prevent any phase-encoding dependent blurring from occurring across the cortical depths in the ROI.
A second MP2RAGE with a high in-plane spatial resolution (0.4x0.4x0.7 mm 3 ) was acquired with the same orientation as the AVF acquisitions. Subsequently, two functional runs were obtained using gradient-echo 3D-EPI 7 with isotropic spatial resolution (0.7mm isotropic) (called Isotropic Voxel EPI, abbreviated as IVE), with the same FoV orientation as the AVF fMRI scans. Lastly, a multiple inversion recovery time EPI (MI-EPI) (0.7mm isotropic) with 64 inversion times 8,9 was acquired for anatomical reference, with the distortions matched to the isotropic 3D-EPI functional scans.
In order to further study the properties of the AVF acquisition, one of the seven subjects participated in an additional fMRI session using the same stimulus paradigm as in the main study. In this session, a high-resolution MP2RAGE was acquired as an anatomical localizer and was followed by an AVF run with highest-resolution sampling the cortical depth (Figure 6a, left). The phase-encoding and frequency-encoding directions were then switched for the second and third functional runs, achieving the highest spatial resolution (0.1mm) along the cortical distance, thereby sampling cortical "columns" (Figure 6a, right). Three additional functional runs were obtained using the AVF with laminar resolution 3 but with increasing number of slices (i.e. increasing coverage), which resulted in consequently longer TRs [Run 1: Slices = 2, TR eff = 3.53 s; Run 2: Slices = 4, TR eff = 7.05 s; Run 3: Slices = 6, TR eff = 10.58 s] ( Figure S2). All other sequence parameters remained the same as in the main study (Table S1).

Functional data
All AVF runs were motion-corrected with Advanced Normalization Tools (ANTs, https://github.com/ANTsX/ANTs) 10,11 using the time-series mean as the reference volume (Translation, antsRegistration). Each motion-corrected functional run was then carefully examined for realignment artefacts, and anatomical features, such as veins (low intensity voxels), were tracked on a volume-byvolume basis as a quality check. Statistical analysis of the fMRI data was done using FEAT 6.0 (FMRI Expert Analysis Tool), part of FSL (FMRIB's Software Library, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) 12 . The grand-mean normalized time-series was temporally filtered (σ high-pass =25xTR, σ low-pass =1.5xTR) and statistical analysis was carried out using FILM with local autocorrelation correction 13 . The GLM model included a double gamma HRF with temporal derivatives and motion parameters as nuisance regressors. The within-subject higher-level analysis was carried out using a fixed effects model by forcing the random effects variance to zero in FLAME (FMRIB's Local Analysis of Mixed Effects) 14,15 .
The quantitative T 1 map from the second MP2RAGE (acquired with the same orientation as the AVF) was registered to the motion-corrected time-series mean AVF fMRI image using the registration tools in ITK-SNAP 3.6 (http://www.itksnap.org) 16 . Note that, in contrast to standard analysis in lowresolution fMRI studies, the anatomy was registered to the functional data in order to keep the functional data in its original space and, hence, avoid any additional smoothing 9 . The images were first aligned using the mutual information cost-function, and further manual fine-tuning of the alignment was performed using ITK-SNAP's interactive registration tool. The transformation matrix was exported in the ITK format and was applied to the anatomical data using ANTs (antsApplyTransforms, 4 th order B-spline interpolation).
All IVE runs were also motion-corrected with ANTs using the time-series mean as the reference volume (Rigid, antsRegistration). No spatial smoothing or distortion-correction was applied (see MI-EPI workflow in Kashyap, et al. 9 ). Further statistical analysis of the fMRI data was done using FEAT 6.0, as described above. A distortion-matched T 1 map was computed by fitting the anatomical MI-EPI data 9 . 4

Segmentation of anatomical images
ITK-SNAP 3.6 was used to semi-automatically segment the subject-specific MP2RAGE and MI-EPI T 1 maps into three tissue classes based on T 1 signal intensity: grey-matter (GM), white-matter (WM) and cerebrospinal fluid (CSF). These binary tissue class masks were manually corrected in ITK-SNAP after co-registration to the functional data, to ensure accurate delineation of the boundaries in the ROI and fix any co-registration artefacts.

Cortical depth analysis
The AVF fMRI time-series data was imported into MATLAB R2015b (Mathworks, USA) for laminar analysis. The average GM thickness in the V1 ROIs across our subjects was ~3 mm, i.e. 30 voxels with 0.1 mm extent in the cortical depth direction. In order to be able to average cortical depth profiles over subjects with varying GM thickness, we linearly interpolated the total number of sampled cortical depths to 30 in each subject's ROI. For the IVE dataset, cortical depth analysis was done using CBS- To check whether similar results are obtained in both acquisition types, if the same ROI is considered, the cortical depth timecourses were averaged along the direction of the cortical distance for all subjects for both AVF and IVE acquisitions. The mean of the positive BOLD and post-stimulus undershoot signals were calculated by averaging over the time windows 12-26 s and 32-42 s following stimulus onset, respectively. The event-related cortical depth timecourses from both datasets were normalized to the mean positive BOLD signal change between 30-70% of cortical depth of each ROI for all subjects, in order to minimize bias to subjects with high BOLD signal changes. The normalized event-related timecourses were then averaged across subjects (Figure 3a). The cortical depth analysis was always carried out using the data across the entire ROI (i.e. no statistical threshold or masking was used).

Variability along cortical distance
Typically, in laminar fMRI studies [21][22][23][24][25] , the sampled depth-specific signals are averaged tangentially along the cortical distance in an ROI. The principle reason for this being, due to the relatively large voxels (~0.7-0.8mm) relative to the histological layer thicknesses, the underlying laminar signals have to be inferred by averaging upsampled fMRI data and accounting for each voxel's partial volumes across the different depths (see section: Effective spatial resolution, below). Due to this limitation, most studies are unable to assess the tangential variability of the acquired signal. In the AVF fMRI data, such upsampling and averaging over an ROI is in principle not necessary, assuming sufficient tSNR. Thus, we assessed the reliability of the AVF fMRI acquisition to probe the tangential variability using Principal Component Analysis (PCA) on the positive BOLD profiles in our ROI with the hypothesis that, if the profiles are similar over trials, then a PCA decomposition over trials should result in most of the variance in the data being explained by the First principal component (PC). For comparison, the PCA was also carried out for the IVE data in its native resolution.

Effective Spatial-Resolution
It is generally assumed that, using the IVE fMRI data, upsampling and averaging over an extensive ROI can in principle yield super-resolved cortical depth profiles with effective spatial resolution only being limited by the extent of the ROI. In order to test this assumption and compare the effective spatial resolution of the IVE with the AVF acquisition, we simulated an anatomically-inspired sample of visual cortex with variable cortical thicknesses (1.7-3.7 mm) at 0.125 mm isotropic resolution 26 ( Supplementary Fig. S3a). Twenty-one equi-volumnar cortical depths were defined, and voxel-partialvolumes were calculated for each cortical depth. Equi-volumnar model of layering was chosen here because of regions of high-curvature in the simulated cortex ( Supplementary Fig. S3b). Two fMRI datasets were generated assuming a Gaussian shaped activation profile in the middle of cortical ribbon with a full-width-at-half-maximum (FWHM) ~1/5 th (Figure 5) of the cortical thickness. The data was then downsampled (by averaging) to a 0.75 mm isotropic resolution. Next, we attempted to reconstruct the original activity profile by sampling the signal with three and twenty-one calculated cortical depths, respectively. Prior to layer sampling, the simulated functional data was upsampled using a nearest-neighbour interpolation method to 0.125 mm resolution and partial volume maps were calculated at this resolution for three and twenty-one numbers of cortical depths, respectively. 6 Signal intensities of small voxels with corresponding partial volume greater than 0.75 within the cortical depth were averaged.
We evaluated (a) the effect of the number of voxels (5-200) within an ROI with randomly assigned location along the cortical ribbon; and (b) the effect of additive noise, with SNR (σ S /σ N ) = 1 (i.e. realistic) and infinite (i.e. ideal) for two sampled number of cortical depth and widths of the activation profile.
The point spread functions (PSF) were first reconstructed and the overlapping area under the true PSF was calculated to quantify blurring induced by subsampling the original signal. All simulations were repeated 200 times for both activation profile widths.

Simulating the large BOLD signal near the CSF/pial boundary
To simulate the effect of leakage of the signal change in the pial vein to the cortical layers due to remaining partial volume effects, we considered activation (amplitude = 5) in the CSF/pial boundary and the very first cortical depth (1 st of 21). The remaining lower 20 cortical depths were considered not to be activated. Two scenarios, with larger (0.75 mm isotropic) and smaller (0.125 mm isotropic) voxel sizes, were compared. This represents a simplified analogy for the acquisition strategies (IVE and AVF) and their differences in laminar spatial resolution. The induced laminar profile (Supplementary Figure S4