LayNii: A software suite for layer-fMRI

Abstract

High-resolution fMRI in the sub-millimeter regime allows researchers to resolve brain activity across cortical layers and columns non-invasively. While these high-resolution data make it possible to address novel questions of directional information flow within and across brain circuits, the corresponding data analyses are challenged by MRI artifacts, including image blurring, image distortions, low SNR, and restricted coverage. These challenges often result in insufficient spatial accuracy of conventional analysis pipelines. Here we introduce a new software suite that is specifically designed for layer-specific functional MRI: LayNii. This toolbox is a collection of command-line executable programs written in C/C++ and is distributed opensource and as pre-compiled binaries for Linux, Windows, and macOS. LayNii is designed for layer-fMRI data that suffer from SNR and coverage constraints and thus cannot be straightforwardly analyzed in alternative software packages. Some of the most popular programs of LayNii contain 'layerification' and columnarization in the native voxel space of functional data as well as many other layer-fMRI specific analysis tasks: layer-specific smoothing, model-based vein mitigation of GE-BOLD data, quality assessment of artifact dominated sub-millimeter fMRI, as well as analyses of VASO data.

Fig. 1. There is a need for a new software suite that is explicitly developed for layer-fMRI.

Panel A) illustrates that because layer-fMRI data suffer from countless shortcomings, the corresponding processing software needs to give the user the flexibility and responsibility to harmonize analysis algorithm parameters for their specific data’s features. Panel B) illustrates that for layer-fMRI data, the possibility and necessity of user-dependent algorithm tweaking can result in reproducible results of higher quality as opposed to enforcing standardized homogeneous user-independent black-box analysis pipelines without any user interventions necessary. This case-specific analysis tweaking required broad knowledge of the user about the algorithm and all its tuning parameters.

Fig. 2. Description of the layerification algorithm.

Panel A) shows an ideal cortex model consisting of sharp tissue boundaries together with terms exemplified on this model. Panel B) shows zoomed-in regions of the model together with a discretized version. This discretized image is analogous to CSF nulled T₁ - weighted MRI data. Panel C) shows the segmentation required by LayNiis layering algorithm. Black voxels indicate the irrelevant voxels, blue voxels indicate the inner gray matter surface, red voxels indicate the outer gray matter surface voxels, and green voxels indicate the pure gray matter voxels. Panel D) visualizes the spatial intuition behind measuring the cortical depth together with showing the measured cortical depths. Note that here the cortical depths are only measured from the outer gray matter surface voxels (red), however LayNii also measures the cortical depth relative to the inner gray matter surface voxels (blue). Therefore, each pure gray matter voxel (green) is described by two cortical depths. Also note that when added, these two cortical depths measure the cortical thickness. Panel E) shows the unit columns defined by LayNii’s algorithm. These unit columns are used to compute the equi-volume metric.

Fig. 3. Layering metrics generated in LayNii.

The top row shows an application with a synthetic 2D image. The middle row shows the empirical layers from Ding et al. (2016) (0.2 mm iso.). The bottom row shows BigBrain (0.1 mm iso., native space) (Amunts et al. 2013) with cortical borders provided in Wagstyl et al. (2020). The equi-distant metric is shown in the middle column and equi-volume metric is shown in the right column for each image type. To better appreciate the difference between the equi-volume and equi-distance layers on the BigBrain data, see the gif animation in Fig. 6) online: https://thingsonthings.org/ln2_layers/. The arrows highlight areas where the equi-distant and the equi-volume metric differ considerably.

Fig. 4. Examples of the layerification in native distorted EPI space.

Example of performing layerification on the whole brain fMRI at 0.8 mm isotropic resolutions. Here, an example is shown with three extracted layers. Typically, at 0.8 mm resolutions, not more than 2-3 layers can be extracted without losing smoothness along the three dimensional cortical folding (bottom left panel). However, when the layerification is performed on a finer spatial grid, the smoothness is improved (top panels).

Fig. 5. Estimating columnar distances in voxel space with LayNii.

Panel A) schematically describes the underlying algorithm of LayNii’s columnar distance estimation. Panel B) depicts the corresponding MRI signal in two independent coordinate systems: a) the scanner coordinate system with folded GM and b) the unfolded cortical ribbon with orthogonalized depth and column distances. The data presented in panels A-B are acquired with an 8 weeks old female cat, Varian 9.4T at CMRR, resolution: 0.125 × 0.125 × 0.5 mm³, Gradient Echo MultiSlice imaging sequence (GEMS, Agilent technology, Inc.) sequence. Panel C) depicts a potential application study of the columnar coordinate system for topographic mapping of functional movement representations. The data presented in panel C are acquired with VASO at a SIEMENS magnetom 7T at FMRIF/NIH with 0.8 mm³ resolution and have been previously described in (Huber et al. 2020b).

Fig. 6

Layer-specific smoothing. Panels A-B) describe two different algorithms of layer-specific smoothing. Panel A) depicts that for spatial smoothing in the folded cortex, an Euclidean distance metric alone might not be a good enough estimate of columnar distances and cannot prohibit signal leakage across kissing gyri. Panel B) depicts example smoothing kernels in contrast-specific smoothing (anisotropic smoothing a.k.a. diffusion filter). Here the mean VASO EPI signal intensity with its inherent T₁ -weighting is used as a reference to generate the voxel-wise smoothing kernel. The elongated shapes of the smoothing kernels depict the weighted signal leakage to voxels with similar anatomical contrast for two example voxels in green and red. The length of the kernel is determined by the FWHM parameter, while the width is determined by the selectivity parameter. Note that those Kernels are highly anisotropic. Panels C-G) exemplify the different results of the respective smoothing methods for an example dataset that comes with the software package. Here the characteristic double-stripe pattern is used as a layer-specific feature that is aimed to be preserved. Panel G) depicts that for extensive signal smoothing with FWHM of 4mm, layer-specific activation can leak across kissing gyri, if not explicitly prevented. Panel G) furthermore shows that the layer-specific double peak is only preserved for very conservative selectivity values. Note that the contrast specific smoothing algorithm with conservative selectivity values preserves the double peak layer activation pattern even for FWHM values of 4mm.

Fig. 7. Conventional and layer-fMRI-specific quality metrics of time series data.

Panels A-D) depict maps of conventional fMRI QA metrics of Mean, tSNR, and activation scores. Panels E-H) depict higher-order QA metrics that can be informative in layer-fMRI. While the QA measures in the top row suggest that the underlying time series is of high quality and that it is not severely limited by artifacts, higher-order QA measures in the bottom row show that there are indeed typical phase errors present in this example time series. For example, the turquoise and red marked areas highlight locations of non-Gaussian signal clutter that are likely referring to EPI phase problems during the readout. While these artifacts are clearly visible in most quality metrics of E-H), they are hardly visible in first order quality metrics A-D).

Fig. 8. Spatiotemporal noise kernel.

Panel A) depicts the major algorithm steps to estimate the noise kernel. Panel B) depicts representative results of the noise kernel in the whole-brain VASO layer-fMRI. It can be seen that the PSF in the second phase encoding direction has negative sidelobes, which suggests that the PSF is not well characterizable with FWHM estimates. The LayNii program LN_NOISE_KERNEL estimates the noise kernel for any time series dataset and writes it out as a three dimensional nifti file. These noise kernel files can then be used to characterize the quality of the time series. E.g., the user can then browse through several noise kernels from multiple layer-fMRI protocols and make judgements about which one is the best. For example, the user can ask questions about which protocol has the least T₂* blurring (pink line in step 6 of panel A), and/or which flip angle scheme in 3D-EPI would result in the smallest negative point-spread function side lobes (blue dotted outline in panel B).

Fig. 9. VASO processing at high magnetic fields with the LayNii program LN_BOCO.

Due to the short T₂* at UHF and the long readout of high-resolution EPI, all VASO data are contaminated with unwanted BOLD contrast. In layer-fMRI application of VASO, this is accounted for by concomitant acquisition of BOLD and VASO in the SS-SI-VASO approach, analogously to ASL-fMRI. The purpose of this figure is to illustrate the relevant processing steps in the LayNii program LN_BOCO to correct for unwanted BOLD contamination in VASO. Panels A-B) depict how the original time series of alternating images with and without blood nulling are temporally interpolated and sorted into two parallel time series of respective contrasts only. Panel C) depicts that in SS-SI-VASO, the BOLD contamination is believed to be a multiplicative factor and be taken care of with a division operation. Panels D-E) depict a representative functional dataset of a flickering checkerboard experiment. It can be seen in the time courses that VASO is anti-correlated to BOLD. While the BOLD signal shows a signal increase during activation, upon BOLD correction, VASO shows a signal decrease. The data presented here have been acquired on a SIEMENS Terra in Glasgow with a segmented-EPI sequence from Ruediger Stirnberg (Stirnberg and Stöcker, 2020).

Fig. 10. Layer-dependent model-based deveining strategies.

The three most often used strategies of layer-dependent vein mitigation are based on a linear offset model, a scaling model, or a leakage model. The respective models are illustrated in panel A). While all models can be used to predict the increasing GE-BOLD signal towards the cortical surface, their assumed physiological signal origin and the corresponding vein mitigation algorithm is fundamentally different. Panel B) exemplifies the application of layer-dependent vein mitigation in LayNii and depicts representative results.