LayNii: A software suite for layer-fMRI

1Functional Magnetic Resonance Imaging (fMRI) aims to measure correlates of brain activity at the level of individual voxels. Recent advances in fMRI hardware, readout strategies, and fMRI contrasts provide estimates of brain activation with voxel-sizes in the sub-millimeter domain. Thus, modern fMRI data from ultra-high fields (UHF) are approaching the spatial scale of cortical layers and columns. While the first decade of sub-millimeter fMRI research focused on data acquisition methods, the methodological research questions of the layer-fMRI field have since shifted towards addressing analysis challenges. As such, a recent survey of the ISMRM study group Current Issues in Brain Function showed that most high-resolution fMRI researchers consider analysis challenges to be more relevant than acquisition challenges.

2The specific shortcomings of common analysis approaches for high-resolution fMRI are listed in multiple review articles, and there are many fMRI analysis software packages that are able to minimize these shortcomings to some degree: AFNI/SUMA (https://afni.nimh.nih.gov,), ANTs (http://stnava.github.io/ANTs/,), BrainVoyager (http://www.brainvoyager.com,, CBSTools/Nighres (http://www.nitrc.org/projects/cbs-tools/,, https://nighres.readthedocs.io/,), FreeSurfer (https://surfer.nmr.mgh.harvard.edu,, FSL (http://fsl.fmrib.ox.ac.uk/fsl,, LIPSIA (https://www.cbs.mpg.de/institute/software/lipsia,), and SPM (http://www.fil.ion.ucl.ac.uk/spm). Some of these packages, however, were originally developed, validated and established for the application for fMRI data that have resolutions in the range of 1.5-5mm, large brain coverage and T2*-weighted fMRI contrasts. Such fMRI data constitute the vast majority of the fMRI literature and are hereafter referred to as “conventional fMRI”. However, layer-fMRI data are different from conventional fMRI data and layer-fMRI data are limited by many unique high-resolution constraints. Thus, layer-fMRI data cannot be straightforwardly analyzed with conventional analysis pipelines (Fig. 1). The specific layer-fMRI data constraints are listed below: Low SNR: Layer-fMRI data are noisier compared to conventional fMRI data. In layer-fMRI studies, tSNR values in the single-digit regime are not uncommon. This indicates that layer-fMRI data are usually more limited by thermal noise than physiological noise. Therefore, many quality assessment (QA) tools that have originally been developed for conventional fMRI data are not applicable to layer-fMRI data. For example, in the thermal noise-dominated regime, the QA metric of tSNR cannot capture many sources of MR-sequence-related and physiological-noise-related signal clutter. These sources of signal instability are below the thermal noise floor (for more discussions of inappropriate QA metrics see Section 2.5). Thus, for layer-fMRI data, new analysis tools with additional (and optimized) QA metrics are needed.

3Low SNR: Layer-fMRI data are noisier compared to conventional fMRI data. In layer-fMRI studies, tSNR values in the single-digit regime are not uncommon. This indicates that layer-fMRI data are usually more limited by thermal noise than physiological noise. Therefore, many quality assessment (QA) tools that have originally been developed for conventional fMRI data are not applicable to layer-fMRI data. For example, in the thermal noise-dominated regime, the QA metric of tSNR cannot capture many sources of MR-sequence-related and physiological-noise-related signal clutter. These sources of signal instability are below the thermal noise floor (for more discussions of inappropriate QA metrics see Section 2.5). Thus, for layer-fMRI data, new analysis tools with additional (and optimized) QA metrics are needed.

4Due to the low tSNR in layer-fMRI, it can become necessary to average multiple data points across space (e.g. spatial smoothing, or voxel pooling). In many layer-fMRI application studies, it can be assumed that neighboring columnar structures are performing the same neural task (e.g. top-down attention) with the identical layer-specific processing signature (albeit potentially different feature representations). In those cases, it can be beneficial to apply local smoothing within the layer direction and without signal blurring across layers to enhance the layer-specific processing signature. Alternatively, it can be beneficial to apply intracortical non-isotropic smoothing based on the directional functional activation or based on the nonisotropic MRI signal intensity across the cortical depth. Thus, for layer-fMRI analyses, new software tools with a larger variety of anatomically informed spatial smoothing methods are needed. Necessity to manually intervene: Accuracy and precision wins over streamlining. In layer-fMRI, the data acquisition procedure is usually pushed to its limits. This means that every acquisition protocol -and image data from each participant- may suffer from individual image artifacts and other experimental shortcomings. While minimal manual tuning of conventional analysis pipelines provides more than sufficient quality for conventional fMRI resolutions, it is not uncommon in layer-fMRI applications to invest 8-12 work hours per dataset for manual high-quality corrections of the segmented borders. Furthermore, it is not uncommon that the experimenter needs to spend more than an entire work week per participant to complete an analysis pipeline. This is due to the layer-fMRI specific requirement to have exceptional accuracy of the analysis and high precision of gray matter (GM) borders. This user-dependent layer-fMRI analysis paradigm is in direct opposition to the demands of large scale population studies of conventional fMRI (UKBiobank, HCP, etc.). In those large scale studies, the analysis pipelines need to be executable with minimal user interventions; and in those population studies, analysis robustness can outrank small accuracy tradeoffs. While these standardized conventional fMRI analyses have gone through several decades of protocol optimizations and aims for stable and robust streamlined analysis that can process thousands of participants without necessary user-interventions, this is neither possible, nor currently demanded for layer-fMRI analyses. In layer-fMRI, the researcher needs to have the flexibility to specifically adjust the analysis pipeline for each individual data set. Thus, a timely layer-fMRI analysis toolbox should be modular and optimized to allow straightforward user-based manipulations of a wide range of algorithm parameters. This would aid the user to iteratively tweak the parameters for every dataset to end up with fine-tuned pipelines that provide the highest possible quality. Such an approach is necessary to optimize layer-fMRI analysis strategies globally in the long term until the data quality will allow streamlined analyses.Restricted coverage: Layer-fMRI acquisition protocols are usually pushed to the limits of what the scanner can achieve. Thus, high spatial resolutions are often achieved by significant tradeoffs in coverage and even single-slice protocols are accepted as a compromise for resolution. These coverage constraints in high resolution fMRI are usually due to the limitations of the MR hardware (mainly gradient performance) and the need for a TR in the range of seconds to allow statistical analysis of the fMRI signal. Furthermore, layer-fMRI data are often acquired with small receiver arrays that have limited coverage and do not facilitate straightforward whole-slice analyses. While it is without question that small fMRI coverage can be troublesome for many analyses, small coverage data play a crucial part in this rapidly developing field. However, many conventional fMRI analysis pipelines are not designed for (or are incompatible with) small coverage data. Therefore, new layer-fMRI analysis tools that are more accommodating to small coverage data are needed.No topology requirements of the cortical sheet: Almost all current layer-fMRI studies use reduced field of views and do not exceed a coverage volume beyond a thin slab. However, many of the current topography-related analysis packages have been designed for whole brain analyses that need to fulfill a number of topological requirements. Namely, fMRI data are used in brain models of closed, continuous GM sheets, without holes or crossing surfaces. However, even though the biological structure of the cerebral cortex fulfils this requirement, common layer-fMRI EPI slab-data do not fulfill this requirement. It is not uncommon for some parts of the functional data to suffer from holes in the cortical sheet and not fulfill the original topology requirements. There are discontinuities, resolution constraints, as well as constraints of brain coverage (commonly 0.5-2.5 cm slabs) that violate the topology constraints in virtually all layer-fMRI data. Therefore, there is a need for analysis software that allows topographical analysis of layers and columns with fewer topological requirements.Necessity to work in the data type of voxels: There are various philosophical approaches to layer-fMRI analyses. Many brain researchers prefer to analyze, depict and interpret their data in a data format that best resembles the object of interest. Since the brain consists of a folded sheet of stacked layers, one could analyze and depict brain data in surface space (using mesh vertices). Alternatively, many MRI and reconstruction methodologists prefer to keep their data as raw (and original) as possible and thus prefer to analyze, depict and interpret their data in the format they are acquired (in voxel space) and without additional necessary data conversions. Since layer-fMRI is still largely limited by the data acquisition strategies, it might be appropriate to favor data processing tools that operate in voxel space to avoid additional confounds that come from transitioning to surface space. Thus, in layer-fMRI, there is a need for purely voxel-based analysis software in addition to already existing surface-based analyses. If researchers had access to layering analyses in voxel space, the data acquisition and analysis could be brought closer together. This would allow researchers to optimize acquisition and analysis concurrently and would enable a new form of fMRI studies that were still unthinkable until today. For example, with layer estimates directly accessible in the raw scanner voxel space, online analysis at the scanner console are doable and layer-dependent neural-feedback studies can become possible.No anatomical reference requirements: Since layer-fMRI data are limited by a high noise level, a successful study often depends on extensive averaging of countless task trials. In light of limited research funding for expensive scan time and ethically appropriate finite scan durations, it is often challenging for the experimenter to acquire additional non-fMRI auxiliary data that might aid the analysis. As such, it can be challenging to obtain high-quality high-resolution whole brain anatomical reference data, B0 and B1 maps, as well as reliable distortion inverted EPI reference data. It is also not uncommon that a researcher is confronted with the decision: either a) to obtain many auxiliary analysis-facilitating reference data without remaining scan time for sufficient functional averages, or b) to obtain sufficient functional averages for decent functional interpretability but without additional reference data. Furthermore, the authors believe that even in cases when these auxiliary data are available, their utilization at submillimeter accuracy level is rarely perfect.While the lack of reference data is not ideal, there is a need for software analyses that can still extract useful information from layer-fMRI without reference data. Thus, there is a need for a layer-fMRI analysis software that is able to perform layer-specific analyses directly in the distorted native EPI space without the requirement of additional non-fMRI data.Non-BOLD fMRI contrasts: At high spatial resolutions, conventional GE-BOLD fMRI sequences are limited by spatially unspecific draining veins that obscure the underlying layer-specific neural activation. Thus, it is getting popular to use alternative fMRI contrasts that do not suffer this shortcoming. One example is the blood volume sensitive VASO (vascular-space-occupancy) contrast that is used in approximately 30% of current layer-fMRI studies. The VASO sequence, however, does not directly provide a single contrast time course like conventional fMRI GE-BOLD data do. Instead, raw VASO time series data consist of multiple interleaved contrasts with variable repetition times (TR) and require additional preprocessing steps to extract pure BOLD and pure bloodvolume weighted time courses, respectively. Thus, there is a need for an analysis software that performs respective VASO processing analyses including temporal image resorting, nifti header manipulations and dynamic contrast divisions.Draining vein problem: In cases where layer-fMRI data have already been acquired with GE-BOLD, the unwanted venous signal can no longer be accounted for with advanced acquisition methods. Instead, one could model the directional blood flow in the cortical vasculature to estimate the venous signal contamination and remove it as much as possible. Thus, a timely layer-fMRI analysis software should contain corresponding tools for model-based venous signal removal.

5Necessity to manually intervene: Accuracy and precision wins over streamlining. In layer-fMRI, the data acquisition procedure is usually pushed to its limits. This means that every acquisition protocol -and image data from each participant- may suffer from individual image artifacts and other experimental shortcomings. While minimal manual tuning of conventional analysis pipelines provides more than sufficient quality for conventional fMRI resolutions, it is not uncommon in layer-fMRI applications to invest 8-12 work hours per dataset for manual high-quality corrections of the segmented borders. Furthermore, it is not uncommon that the experimenter needs to spend more than an entire work week per participant to complete an analysis pipeline. This is due to the layer-fMRI specific requirement to have exceptional accuracy of the analysis and high precision of gray matter (GM) borders. This user-dependent layer-fMRI analysis paradigm is in direct opposition to the demands of large scale population studies of conventional fMRI (UKBiobank, HCP, etc.). In those large scale studies, the analysis pipelines need to be executable with minimal user interventions; and in those population studies, analysis robustness can outrank small accuracy tradeoffs. While these standardized conventional fMRI analyses have gone through several decades of protocol optimizations and aims for stable and robust streamlined analysis that can process thousands of participants without necessary user-interventions, this is neither possible, nor currently demanded for layer-fMRI analyses. In layer-fMRI, the researcher needs to have the flexibility to specifically adjust the analysis pipeline for each individual data set. Thus, a timely layer-fMRI analysis toolbox should be modular and optimized to allow straightforward user-based manipulations of a wide range of algorithm parameters. This would aid the user to iteratively tweak the parameters for every dataset to end up with fine-tuned pipelines that provide the highest possible quality. Such an approach is necessary to optimize layer-fMRI analysis strategies globally in the long term until the data quality will allow streamlined analyses.

6Restricted coverage: Layer-fMRI acquisition protocols are usually pushed to the limits of what the scanner can achieve. Thus, high spatial resolutions are often achieved by significant tradeoffs in coverage and even single-slice protocols are accepted as a compromise for resolution. These coverage constraints in high resolution fMRI are usually due to the limitations of the MR hardware (mainly gradient performance) and the need for a TR in the range of seconds to allow statistical analysis of the fMRI signal. Furthermore, layer-fMRI data are often acquired with small receiver arrays that have limited coverage and do not facilitate straightforward whole-slice analyses. While it is without question that small fMRI coverage can be troublesome for many analyses, small coverage data play a crucial part in this rapidly developing field. However, many conventional fMRI analysis pipelines are not designed for (or are incompatible with) small coverage data. Therefore, new layer-fMRI analysis tools that are more accommodating to small coverage data are needed.

7No topology requirements of the cortical sheet: Almost all current layer-fMRI studies use reduced field of views and do not exceed a coverage volume beyond a thin slab. However, many of the current topography-related analysis packages have been designed for whole brain analyses that need to fulfill a number of topological requirements. Namely, fMRI data are used in brain models of closed, continuous GM sheets, without holes or crossing surfaces. However, even though the biological structure of the cerebral cortex fulfils this requirement, common layer-fMRI EPI slab-data do not fulfill this requirement. It is not uncommon for some parts of the functional data to suffer from holes in the cortical sheet and not fulfill the original topology requirements. There are discontinuities, resolution constraints, as well as constraints of brain coverage (commonly 0.5-2.5 cm slabs) that violate the topology constraints in virtually all layer-fMRI data. Therefore, there is a need for analysis software that allows topographical analysis of layers and columns with fewer topological requirements.

8Necessity to work in the data type of voxels: There are various philosophical approaches to layer-fMRI analyses. Many brain researchers prefer to analyze, depict and interpret their data in a data format that best resembles the object of interest. Since the brain consists of a folded sheet of stacked layers, one could analyze and depict brain data in surface space (using mesh vertices). Alternatively, many MRI and reconstruction methodologists prefer to keep their data as raw (and original) as possible and thus prefer to analyze, depict and interpret their data in the format they are acquired (in voxel space) and without additional necessary data conversions. Since layer-fMRI is still largely limited by the data acquisition strategies, it might be appropriate to favor data processing tools that operate in voxel space to avoid additional confounds that come from transitioning to surface space. Thus, in layer-fMRI, there is a need for purely voxel-based analysis software in addition to already existing surface-based analyses. If researchers had access to layering analyses in voxel space, the data acquisition and analysis could be brought closer together. This would allow researchers to optimize acquisition and analysis concurrently and would enable a new form of fMRI studies that were still unthinkable until today. For example, with layer estimates directly accessible in the raw scanner voxel space, online analysis at the scanner console are doable and layer-dependent neural-feedback studies can become possible.

9No anatomical reference requirements: Since layer-fMRI data are limited by a high noise level, a successful study often depends on extensive averaging of countless task trials. In light of limited research funding for expensive scan time and ethically appropriate finite scan durations, it is often challenging for the experimenter to acquire additional non-fMRI auxiliary data that might aid the analysis. As such, it can be challenging to obtain high-quality high-resolution whole brain anatomical reference data, B0 and B1 maps, as well as reliable distortion inverted EPI reference data. It is also not uncommon that a researcher is confronted with the decision: either a) to obtain many auxiliary analysis-facilitating reference data without remaining scan time for sufficient functional averages, or b) to obtain sufficient functional averages for decent functional interpretability but without additional reference data. Furthermore, the authors believe that even in cases when these auxiliary data are available, their utilization at submillimeter accuracy level is rarely perfect.

10While the lack of reference data is not ideal, there is a need for software analyses that can still extract useful information from layer-fMRI without reference data. Thus, there is a need for a layer-fMRI analysis software that is able to perform layer-specific analyses directly in the distorted native EPI space without the requirement of additional non-fMRI data.

11Non-BOLD fMRI contrasts: At high spatial resolutions, conventional GE-BOLD fMRI sequences are limited by spatially unspecific draining veins that obscure the underlying layer-specific neural activation. Thus, it is getting popular to use alternative fMRI contrasts that do not suffer this shortcoming. One example is the blood volume sensitive VASO (vascular-space-occupancy) contrast that is used in approximately 30% of current layer-fMRI studies. The VASO sequence, however, does not directly provide a single contrast time course like conventional fMRI GE-BOLD data do. Instead, raw VASO time series data consist of multiple interleaved contrasts with variable repetition times (TR) and require additional preprocessing steps to extract pure BOLD and pure bloodvolume weighted time courses, respectively. Thus, there is a need for an analysis software that performs respective VASO processing analyses including temporal image resorting, nifti header manipulations and dynamic contrast divisions.

12Draining vein problem: In cases where layer-fMRI data have already been acquired with GE-BOLD, the unwanted venous signal can no longer be accounted for with advanced acquisition methods. Instead, one could model the directional blood flow in the cortical vasculature to estimate the venous signal contamination and remove it as much as possible. Thus, a timely layer-fMRI analysis software should contain corresponding tools for model-based venous signal removal.

13In short, layer-fMRI data often suffer from multiple shortcomings that hinder a straightforward analysis of layers and columns. While the acquisition methodology needs to be further advanced to remove those shortcomings, there is a current need for a layer-analysis software that works despite each and all of the current shortcomings (Fig. 1).

14Here, we present an fMRI analysis software suite LayNii that is specifically designed for layer-fMRI data and that addresses all of the above listed challenges. LayNii is designed to perform layer analyses entirely in voxel space and consists of many modular programs that perform layer-fMRI specific tasks: Cortical depth and thickness measurements with corresponding assignments of layer tags to each GM voxel (a.k.a. layerification) in any segmented volume data set,Estimation of columnar structures across voxels in any volume data set,Performing various ways of layer-specific smoothing,Providing a set of appropriate QA metrics for thermal noise-limited layer-fMRI data,Pre-processing blood-volume based VASO fMRI contrasts that are optimized for the layer-specific microvascular responses,Applying model-based deveining algorithms of layer-fMRI GE-BOLD data.

15Cortical depth and thickness measurements with corresponding assignments of layer tags to each GM voxel (a.k.a. layerification) in any segmented volume data set,

16Estimation of columnar structures across voxels in any volume data set,

17Performing various ways of layer-specific smoothing,

18Providing a set of appropriate QA metrics for thermal noise-limited layer-fMRI data,

19Pre-processing blood-volume based VASO fMRI contrasts that are optimized for the layer-specific microvascular responses,

20Applying model-based deveining algorithms of layer-fMRI GE-BOLD data.

21The LayNii software suite is available as stand-alone C++ source code without external dependencies or needing libraries. It is distributed as pre-compiled (binary) installation packages for all major operating systems (Linux, Windows, macOS). LayNii can read and write nifti files (Neuroimaging Informatics Technology; https://nifti.nimh.nih.gov/) in uncompressed NII and compressed NII.GZ formats.

22The name LayNii is derived from the two words “Layer-fMRI” and “Nifti”. The prefix LAY stresses that this suite is particularly built for layer-fMRI (which includes sub-millimeter fMRI and columnar fMRI). The suffix NII emphasizes that this software suite operates in the voxel space of nifti data.

23The structure of LayNii is largely inspired by the philosophy of alternative software packages and is aimed to be as modular as possible. It is designed to consist of multiple individual lightweight programs that are executable from the command line and can be combined in pipeline scripts that may or may not include additional elements of alternative software suites like, AFNI, fslmaths, ANTs, etc. LayNii’s modularity will also help to straightforwardly include LayNii’s programs into already established software suites in the future. We follow the AFNI principle of providing mechanisms, not policies. We aim to give the user the power to assemble computing pieces in different ways to make customized analysis, which in turn means that it is the user’s responsibility to know what the individual programs with all its tuning parameters do. This might give LayNii a high flexibility to increase the precision and accuracy for each given data set at the cost of straightforward streamlining. Adjusting LayNii’s tuning parameters can have a large influence on how the final results look. For the sake of reproducibility, we encourage users of LayNii to report the version number (each LayNii release has a unique DOI using Zenodo) and all specifics of the executed LayNii commands (input parameters) in their respective publications. Currently, the aim of the LayNii package is not to provide a one-shot analysis solution for entire pipelines, as this is already well-covered by other available packages. There is version control of LayNii, so researchers can reproduce results from older studies even after several updates.

24Some of the most essential LayNii programs are discussed below. An exhaustive list of additional programs and tutorials of all (>30) programs is available on Github. Each individual program has a -help option that describes the main functionality and how to use it.

25Layer-fMRI necessitates determining the cortical depth of each gray matter voxel. This information can be used to pool together every voxel within a certain cortical depth range to generate “layers”. Since the cortex varies in thickness across regions, absolute (as opposed to relative) cortical depth measurements are often not very useful to consistently determine layers. To account for this variation across regions, local cortical thickness measurements are used to normalize the cortical depths. We refer to the resulting normalized cortical depth measurement as “equi-distant metric”. The equi-distant metric yields equi-distant layers, when quantized.

26In addition to cortical thickness, the location and size of the cortical neurobiological layers also vary with regards to the cortical curvature. The effect of cortical curvature on layering was originally demonstrated by Bok in 1929, where he subdivided the cortex into columns to show that deeper layers get thicker and superficial layers get thinner around gyri. Meanwhile, deeper layers get thinner and superficial layers get thicker around sulci. This observation is known as the equi-volume principle.

27In LayNii, we provide both equi-distant and equi-volume layering options to the user. While solutions already exist, our implementation differs in a unique way. It is computed completely on the discrete lattice (as opposed to by means of surface approximations). To give a description of our layering algorithm, first we define the terms that will be referred to later (also see Fig. 2)): Layer(s): In neuroscience, the term “cortical layers” often refers to neurobiological layers (e.g. histologically defined). Here we are not using layers to exclusively mean neurobiological layers, but instead to mean “a thickness of some voxels laid within cortical gray matter” which may or may not correspond to neurobiological layers. Due to the thickness of neurobiological layers and the typical spatial resolution used in depth-dependent fMRI, an actual matching of artificial and neurobiological layers can only be achieved in case of layerification with higher resolutions than the native spatial resolution of the fMRI sequence used.Inner gray matter surface: Portion of gray matter that mostly faces white matter.Outer gray matter surface: Portion of gray matter that mostly faces cerebrospinal fluid.Streamline: A line that connects inner and outer gray matter surfaces based on some principle (not necessarily the shortest Euclidean distance).Cortical thickness: Shortest streamline distance between outer and inner gray matter borders.Cortical depth: Distance from inner or outer gray matter surface, a portion of cortical thickness. Slightly different from layer because this term indicates a quantitative measure which can be used to define the neurobiological layers.Column: A group of voxels that are successively penetrated by one or multiple streamlines. Not to be confused with neurobiological cortical columns which indicates a group of neurons, not voxels.Unit column: Smallest units of columns that can be defined in a discrete lattice (ordered set of voxels) given the image resolution and topology of the cortex. These unit columns are solely used for the purpose of facilitating subsequent layerification. Their definition is chosen to be algorithmically convenient and not intended to be physiologically plausible in all instances. Their width is chosen to be as thin as possible while touching at least one outer gray matter voxel and at least one inner gray matter voxel. This definition of unit columns should not be confused with the definition of cortical distances in Section 2.3 (below).Curvature: A measure of ‘how much a shape bends’. More specifically, here we are defining it as a measure of ‘how much a unit column bends’. E.g. bananashaped or pyramidal-shaped.Rim: A volume of cortical gray matter. Here, we specifically use this term to indicate our input image, which consists of four integer values (0=irrelevant voxels, 1=inner gray matter surface voxels, 2=outer gray matter surface voxels, 3=pure gray matter voxels; see panels of Fig. 2).

28Layer(s): In neuroscience, the term “cortical layers” often refers to neurobiological layers (e.g. histologically defined). Here we are not using layers to exclusively mean neurobiological layers, but instead to mean “a thickness of some voxels laid within cortical gray matter” which may or may not correspond to neurobiological layers. Due to the thickness of neurobiological layers and the typical spatial resolution used in depth-dependent fMRI, an actual matching of artificial and neurobiological layers can only be achieved in case of layerification with higher resolutions than the native spatial resolution of the fMRI sequence used.

29Inner gray matter surface: Portion of gray matter that mostly faces white matter.

30Outer gray matter surface: Portion of gray matter that mostly faces cerebrospinal fluid.

31Streamline: A line that connects inner and outer gray matter surfaces based on some principle (not necessarily the shortest Euclidean distance).

32Cortical thickness: Shortest streamline distance between outer and inner gray matter borders.

33Cortical depth: Distance from inner or outer gray matter surface, a portion of cortical thickness. Slightly different from layer because this term indicates a quantitative measure which can be used to define the neurobiological layers.

34Column: A group of voxels that are successively penetrated by one or multiple streamlines. Not to be confused with neurobiological cortical columns which indicates a group of neurons, not voxels.

35Unit column: Smallest units of columns that can be defined in a discrete lattice (ordered set of voxels) given the image resolution and topology of the cortex. These unit columns are solely used for the purpose of facilitating subsequent layerification. Their definition is chosen to be algorithmically convenient and not intended to be physiologically plausible in all instances. Their width is chosen to be as thin as possible while touching at least one outer gray matter voxel and at least one inner gray matter voxel. This definition of unit columns should not be confused with the definition of cortical distances in Section 2.3 (below).

36Curvature: A measure of ‘how much a shape bends’. More specifically, here we are defining it as a measure of ‘how much a unit column bends’. E.g. bananashaped or pyramidal-shaped.

37Rim: A volume of cortical gray matter. Here, we specifically use this term to indicate our input image, which consists of four integer values (0=irrelevant voxels, 1=inner gray matter surface voxels, 2=outer gray matter surface voxels, 3=pure gray matter voxels; see panels of Fig. 2).

38In what follows, we detail the implementation steps of layerification in LayNii. The layerification algorithm in LayNii starts from a segmented rim image (see Fig. 2) Panel C). This segmented rim image contains the segmentation of GM and its borders. The segmentation itself is not done in LayNii (see Section 3.1). The LayNii program LN_RIMIFY can convert the segmentation output of alternative software packages to a LayNii-optimized rim file.For each voxel, we measure two cortical depths. One relative to the inner gray matter surface voxel and another relative to the outer gray matter surface voxels. We compute these cortical depths by using an iterative growing algorithm.By adding relative distances to inner and outer gray matter surfaces, we compute the cortical thickness per voxel.Cortical depth relative to inner gray matter surface voxels divided by thickness provides an equi-distant metric (ranges between 0 to 1).One streamline per gray matter voxel (green) is defined by the closest outer (red) and the inner (blue) gray matter surface voxels. We use the streamlines to describe the curvature type (as either gyrus, sulcus or straight wall). This is done by means of counting the unique inner and outer gray matter surface voxels each streamline of a unit column connects to. Note that this method avoids traditional curvature computations while still yielding the useful curvature type information.Then we define our unit columns. A unit column has to touch at least one outer gray matter voxel and at least one inner gray matter voxel. A unit column might correspond exactly to one streamline where the cortex curvature is mostly straight (and thin). However, unit columns will often correspond to multiple streamlines when the cortex is curved (and thick). Such unit columns exhaustively cover all gray matter voxels (see illustrative examples in Fig. 2E). The width of the unit column is as small as the voxel resolution allows. They are defined to have at least one voxel of the inner and at least one voxel of the outer gray matter.By exploiting the equi-distant metric, we count the voxels that fall close to the superficial side versus the deeper side of the middle gray matter surface.In order to ultimately estimate equi-volume layers, we compute how much and to which direction the middle gray matter should be pushed to in order to balance the number of voxels on each side. We call the resulting numbers equi-volume factors.Then we exploit the equi-distant metric as a constrained vector space (adding up to always one when the normalized cortical depth relative to inner and outer gray matter surfaces are added together). This constrained vector space can be recognized as a simplex space of two dimensions. The n-dimensional simplex space was first described by August Ferdinand Möbius in 1827 under the name “barycentric coordinates” and brought back to attention in the modern era by John Aitchison in 1986 under the name compositional data. We perform a mathematical translation operation that is defined as a linear operator within the simplex space (called perturbation in ) by using the equivolume factors. This translation balances the metric space in a way to yield the equi-volume metric. Note that this operation is done on a per unit-column basis.Since the unit-columns have sharp transitions in between, it can be advantageous to smooth these transitions to have a more biologically plausible equivolume metric distribution over gray matter voxels. We use iterative smoothing with a small kernel to prevent leakage around kissing gyri and tight sulci.As the final step, we quantize the equi-distant and equivolume metrics to provide the desired number of layers to the user. However, note that we are also providing the equi-distant and equi-volume metrics as default outputs. This is advantageous to allow the users to choose their own quantization methods.

39The layerification algorithm in LayNii starts from a segmented rim image (see Fig. 2) Panel C). This segmented rim image contains the segmentation of GM and its borders. The segmentation itself is not done in LayNii (see Section 3.1). The LayNii program LN_RIMIFY can convert the segmentation output of alternative software packages to a LayNii-optimized rim file.

40For each voxel, we measure two cortical depths. One relative to the inner gray matter surface voxel and another relative to the outer gray matter surface voxels. We compute these cortical depths by using an iterative growing algorithm.

41By adding relative distances to inner and outer gray matter surfaces, we compute the cortical thickness per voxel.

42Cortical depth relative to inner gray matter surface voxels divided by thickness provides an equi-distant metric (ranges between 0 to 1).

43One streamline per gray matter voxel (green) is defined by the closest outer (red) and the inner (blue) gray matter surface voxels. We use the streamlines to describe the curvature type (as either gyrus, sulcus or straight wall). This is done by means of counting the unique inner and outer gray matter surface voxels each streamline of a unit column connects to. Note that this method avoids traditional curvature computations while still yielding the useful curvature type information.

44Then we define our unit columns. A unit column has to touch at least one outer gray matter voxel and at least one inner gray matter voxel. A unit column might correspond exactly to one streamline where the cortex curvature is mostly straight (and thin). However, unit columns will often correspond to multiple streamlines when the cortex is curved (and thick). Such unit columns exhaustively cover all gray matter voxels (see illustrative examples in Fig. 2E). The width of the unit column is as small as the voxel resolution allows. They are defined to have at least one voxel of the inner and at least one voxel of the outer gray matter.

45By exploiting the equi-distant metric, we count the voxels that fall close to the superficial side versus the deeper side of the middle gray matter surface.

46In order to ultimately estimate equi-volume layers, we compute how much and to which direction the middle gray matter should be pushed to in order to balance the number of voxels on each side. We call the resulting numbers equi-volume factors.

47Then we exploit the equi-distant metric as a constrained vector space (adding up to always one when the normalized cortical depth relative to inner and outer gray matter surfaces are added together). This constrained vector space can be recognized as a simplex space of two dimensions. The n-dimensional simplex space was first described by August Ferdinand Möbius in 1827 under the name “barycentric coordinates” and brought back to attention in the modern era by John Aitchison in 1986 under the name compositional data. We perform a mathematical translation operation that is defined as a linear operator within the simplex space (called perturbation in ) by using the equivolume factors. This translation balances the metric space in a way to yield the equi-volume metric. Note that this operation is done on a per unit-column basis.

48Since the unit-columns have sharp transitions in between, it can be advantageous to smooth these transitions to have a more biologically plausible equivolume metric distribution over gray matter voxels. We use iterative smoothing with a small kernel to prevent leakage around kissing gyri and tight sulci.

49As the final step, we quantize the equi-distant and equivolume metrics to provide the desired number of layers to the user. However, note that we are also providing the equi-distant and equi-volume metrics as default outputs. This is advantageous to allow the users to choose their own quantization methods.

50This layering procedure works in isotropic and asymmetric voxels without directionality biases.

51Layering in LayNii works directly on 2D and 3D images. In addition, it works equally well for partial coverage and whole brain images. As an example, see top row of Fig. 3 for layering applied on synthetic/simulated a 2D image. For reallive examples of 2D-layerification on individual histological slices and pictures thereof see https://layerfmri.com/2dlayers/. The middle row of Fig. 3 depicts an example of a 3D partial coverage image (cut from respective whole brain data) and bottom row depicts a whole brain image.

52Since the LayNii algorithms are operating directly in voxel space, the layerification can be significantly faster than conventional approaches. For instance, the duration for wholebrain 0.5 mm iso. data processing on a single CPU takes approximately 8 sec for equi-distancing and 48 sec for equivoluming. The corresponding computation times for 0.2 mm whole brain data (BigBrain) takes 2 and 12 min respectively. For more information on the computation times, see: https://thingsonthings.org/ln2_layers/.

53For a practical (video-) tutorial of performing layerification directly on layer-fMRI EPI data see here: https://layerfmri.com/analysispipeline/and here: https://layerfmri.com/quick-layering/. For an in-depth discussion and comparison of equidistance layering, equi-volume layering and additional layering algorithms that are implemented in LayNii see: https://layerfmri.com/equivol/.

54Fig. 4 depicts examples of the layerification in native distorted EPI space. While it is often advantageous to perform the layerification directly in the native distorted EPI space, it is not necessary (nor advised) to perform the layerification in the same spatial grid size. Depending on the number of extracted layers, it can be advantageous to perform the layerification of a finer spatial voxel grid. This can be done, while leaving the raw fMRI data in the original spatial grid and is comparable to vertex coordinates that are defined with finder precision as the voxels. For more hand on examples of how to do this with LayNii, see https://layerfmri.com/how-many-layers-should-i-reconstruct/.

55While voxels with the same relative distance to the GM surfaces are collected into bins of layers, it can also be possible to group voxels into bins of the same columnar distances, in the orthogonal direction to layers, where the distances are measured from specified landmarks. Such estimates of columnar structures that span across all layers of the cortical depth are also often needed for many steps of layer-fMRI analysis. That is, in studies where layer-specific activity cannot be assumed to be identical across large patches of GM, it becomes necessary to perform simultaneous laminar and columnar analyses. As such, estimates of columnar distances are vital for research questions addressing the topographical distributions of layer-dependent brain activation. Past examples are a) depictions of somatotopically aligned body-part representations in the primary sensory cortex, b) layer-dependent signal distribution across visual eccentricity, c) for columnar-specific functional hierarchy mapping, d) for cortical unfolding, or e) the analysis of tapping induced activation patches that do not extend across a cortical patch of few millimeters in the primary motor cortex. In LayNii, columnar distances are calculated in a six-step algorithm that is schematically illustrated in Fig. 5A: In contrast to layer estimates, there is usually no clear physiologically defined coordinate system origin. While layer distances are inherently normalized to a coordinate system between WM and CSF borders, the origin of columnar estimates is highly dependent on the specific study and research question. Thus in LayNii, columnar distance estimates are generated based on a manually set landmark (Fig. 5A). This is the starting point of the algorithm (white arrow). In past applications of this algorithm, LayNii users placed the landmarks at anatomical reference points (e.g. position of smallest curvature radius in hand knob) or functional reference points (e.g. position of V1-V2 border, border of thumb finger representation in S1).This landmark is used as the origin of a subsequence growing algorithm. In order to account for the fact that the GM of a neighboring sulcus can be closer than the distance of a column within a sulcus of kissing gyri, the distance is estimated here with a local grow-algorithm that extends the local patch of columnar distances iteratively voxel-by-voxel (step-size=1 voxel).Since the growing algorithm works in voxel space, it can only estimate distances in units of integer multiples of voxel distances. This means that a diagonal voxel neighbour is estimated as being two steps away (one step in each orthogonal direction), whereas the Euclidean distance is actually smaller (square root of two). The corresponding Pythagorean errors are mitigated by the application of an additional processing step that applies local smoothing along the layer direction only.In the next step, the columnar distances are extrapolated across layers. The columnar distance value of every voxel is simply inherited from the next closest voxel that has a determined distance estimate already. This is done within the GM ribbon only (to avoid leakage from neighboring sulci).Again, this step comes along with Pythagorean errors (mismatch of growing iteration steps and Euclidean distance) that can be accounted for with within-layer smoothing.Finally, the user can choose the desired column thickness. This is done with the optional -Ncolumns flag of the LayNii program. When this parameter is not specified by the user, the maximum number of columns is used, and the resulting columnar estimates are one voxel thick. For smaller values of -Ncolumns, the columns become thicker. The width of the columnar distances is defined to be equal in the middle layer along the cortical ribbon. This means that at locations of strong cortical curvature, the columnar distances have a shape of a frustum of a pyramid exhibiting different widths across the cortical depth.

56In contrast to layer estimates, there is usually no clear physiologically defined coordinate system origin. While layer distances are inherently normalized to a coordinate system between WM and CSF borders, the origin of columnar estimates is highly dependent on the specific study and research question. Thus in LayNii, columnar distance estimates are generated based on a manually set landmark (Fig. 5A). This is the starting point of the algorithm (white arrow). In past applications of this algorithm, LayNii users placed the landmarks at anatomical reference points (e.g. position of smallest curvature radius in hand knob) or functional reference points (e.g. position of V1-V2 border, border of thumb finger representation in S1).

57This landmark is used as the origin of a subsequence growing algorithm. In order to account for the fact that the GM of a neighboring sulcus can be closer than the distance of a column within a sulcus of kissing gyri, the distance is estimated here with a local grow-algorithm that extends the local patch of columnar distances iteratively voxel-by-voxel (step-size=1 voxel).

58Since the growing algorithm works in voxel space, it can only estimate distances in units of integer multiples of voxel distances. This means that a diagonal voxel neighbour is estimated as being two steps away (one step in each orthogonal direction), whereas the Euclidean distance is actually smaller (square root of two). The corresponding Pythagorean errors are mitigated by the application of an additional processing step that applies local smoothing along the layer direction only.

59In the next step, the columnar distances are extrapolated across layers. The columnar distance value of every voxel is simply inherited from the next closest voxel that has a determined distance estimate already. This is done within the GM ribbon only (to avoid leakage from neighboring sulci).

60Again, this step comes along with Pythagorean errors (mismatch of growing iteration steps and Euclidean distance) that can be accounted for with within-layer smoothing.

61Finally, the user can choose the desired column thickness. This is done with the optional -Ncolumns flag of the LayNii program. When this parameter is not specified by the user, the maximum number of columns is used, and the resulting columnar estimates are one voxel thick. For smaller values of -Ncolumns, the columns become thicker. The width of the columnar distances is defined to be equal in the middle layer along the cortical ribbon. This means that at locations of strong cortical curvature, the columnar distances have a shape of a frustum of a pyramid exhibiting different widths across the cortical depth.

62Note that the columnar distances can be bent in areas of strong cortical folding (Fig. 5A, step 6). While there are cytoarchitectonic and angioarchitectonic indications that this might be physiologically plausible, a general and quantitative parametrization of the columns bending and shape has not been described in the literature so far. Thus, further research is needed to investigate the physiological accuracy of the columnar bending in LayNii. A more detailed discussion of the columnar distance shape in LayNii is given in: https://layerfmri.com/equivol/.

63These columnar distance estimates can then be further used for cortical unfolding Fig. 5B.

64Fig. 5C depicts an example application of LayNii’s columnar distance estimation and subsequent cortical flattening in three dimensional voxel space. Here the mirrored finger representations in the primary motor system are used as a toy-model to exemplify the purpose of the columnar distance analyses in LayNii. The left panels depict the estimates of the columnar coordinates in the native volume space of T1-weighted functional EPI. Here, two orthogonal axes of columnar distances are used. One axis goes along the medial-lateral direction (colors left-right), and one axis goes across the direction of the depth of the central sulcus (upwards-downwards). This spans a locally confined orthogonal grid of columnar distances across the entire anterior bank of the central sulcus and can be used to extract functional signals across all layers (third dimension). These three units of distances (medial-lateral, upwards-downwards, and layers). Can be used to re-grid the signal into a new nifti file with the LayNii program LN_IMAGIRO, by simple signal value copying of each voxel to another nifti file based on the three new coordinates. The result is shown on the right panel of Fig. 5C as flattened cortical patches of functional finger dominance maps in 3D-nifti pace. It is shown how the flattened volume representation of the data facilitates straightforward analysis across columnar profiles projected across cortical depths.

65A more detailed explanation of the columnarization in the LayNii program LN_3DCOLUMNS, and the unfolding program LN_IMAGIRO is given here: https://layerfmri.com/columns/.

66Due to the low tSNR in layer-fMRI, it can become advantageous to apply spatial smoothing. In order not to compromise the locally specific layer-fMRI signal of interest, the spatial smoothing kernel needs to be exclusively applied in anatomically informed directions only. Spatial smoothing in specific layer directions has been originally proposed for vertex-based analyses in, re-implemented for voxel-based layer smoothing in, and it has been ultimately described in full depth later. This form of layer-specific smoothing is commonly applied locally in a way that the signal leakage drops off based on a Gaussian kernel of a given full-width-half-maximum (FWHM), analogous to conventional isotropic smoothing. In contrast to conventional isotropic smoothing, however, the layer-specific smoothing has an additional signal leakage drop off penalty perpendicularly to the layer-direction. Thus, for layer-specific smoothing, the signal leakage is restricted across voxels in different layers, despite the case that they might be very close to each other in Euclidean space. In the LayNii program, LN_LAYER_SMOOTH this is implemented as volume smoothing in voxel space by restricting the smoothing kernel across voxels whose centroids are located in the same cortical depth. The smoothing is thus applied within masks of previously defined layers only.

67This form of layer-specific smoothing, can be seen as isotropic smoothing within a collection of binary masks (comparable to AFNIs 3dBlurInMask). So if the user chooses to have very few (but thick) layers, smoothing would be isotropic within each of those layers, with hard boundaries in between layers where no smoothing occurs. On the other hand, when the user would choose to have as many as 100 (very thin) layers, the smoothing algorithm solely smoothes across voxels that are in the same percent-bin of cortical depth. Using more layers, reduces signal leakage across cortical depth within the same depth-bin.

68When the smoothing kernel sizes are particularly large and approach the spatial scale of the distance between touching gyri, it becomes necessary for the algorithm to consider signal leakage across sulci in very close proximity (Fig. 6A). While the problem of signal leakage across neighboring gyri is inherently avoided in surface-based smoothing approaches, it needs to be explicitly taken care of in layer-specific volume smoothing programs of LayNii. In the LayNii implementation of layer-specific smoothing, a specific -NoKissing flag ensures that smoothing does not occur blindly across all voxels in close proximity, but only in those voxels that are connected to the same GM area.

69While this form of intracortical smoothing requires accurate knowledge of the local layer geometry within the cortex, there are alternative approaches to apply layer-specific smoothing in preferred layer-directions without the need to have predefined layer labels available. For example, it has been suggested to apply intracortical non-isotropic smoothing based on the local functional activation or based on the non-isotropic MRI signal intensity across the cortical depth. In LayNii, this form of smoothing is implemented in the program LN_GRADSMOOTH by taking advantage of the fact that highresolution EPI intensities contain informative anatomical information too. This form of smoothing estimates the signal leakage kernel with two penalty functions; a) Gaussian weighted Euclidean distance of seed and target voxel, and b) local entropy weighted signal difference of the seed and target voxel. Just like for isotropic and layer-smoothing, the size of the Gaussian kernel can be adjusted by means of a user-defined -FWHM parameter in units of millimeter. The entropy weighted signal difference can be adjusted by the user by means of a -selectivity parameter. E.g., for a common - selectivity parameter of 0.08, the signal leakage is restricted to the local voxels that have a similar signal intensity, which differ from the central voxel by less than 8% of the local signal variance. For a more conservative -selectivity parameter of 0.04, the signal leakage is restricted to the local voxels that have a similar signal intensity, which differ from the central voxel by less than 4% of the local signal variance (see examples in Fig. 6B). The optimal selectivity value depends on the image’s contrast, the acceptable blurring across voxels with similar contrast, and the desired smoothness parameter (FWHM). Since the smoothing kernel is determined in a data-driven way, its application is significantly more convenient for the researcher to use. Labor intensive segmentation is not necessary to be able to perform this way of smoothing.

70This algorithm of intensity-based smoothing is very similar to a large number of previously described edge-preserving smoothing methods, which are implemented in AFNI, FSL, and CBS tools (which is also inherited by the Nighres wraper) as 3danisosmooth, LaminarIterativeSmoothing, or SUSAN, and some of them are also referred to as diffusion filter. The LayNii implementation LN_GRADSMOOTH differs from those algorithms in the sense that it is optimized for applications in layer-fMRI. Since layer-fMRI data are usually very noisy, the smoothing kernel should not be generated on the actual time series data or the statistical activation data themselves. Instead, the LayNii implementation estimates the smoothing kernel on an independent optimized input file (given to LN_GRADSMOOTH with the -gradfile parameter), e.g., a mean EPI image with anatomical contrast (T1-EPI). This is readily doable in layer-fMRI data because of the exquisite structural details visible in sub-millimeter EPI. Furthermore, the LayNii implementation can avoid signal leakage of unconnected islands of similar signals intensities across kissing gyri (Fig. 6A).

71This form of signal intensity based smoothing can also be used in combination with layer masks as described for LN_LAYER_SMOOTH. In this case, the originally hard borders of smoothing across layer-bins can be softened with larger -selectivity parameter values.

72These forms of anatomically informed gradient smoothing and layer smoothing are particularly beneficial in layer-fMRI studies where it can be assumed that neighboring columns perform the same neural task with the identical layer-dependent activation signature. And they are counter-productive, when signal differences of neighboring columnar structures are of interest.

73Fig. 6C-G shows representative example results of anatomically informed smoothing applications. It can be seen that the characteristic double peak in the primary motor cortex is preserved across a wide range of algorithm parameters. Furthermore, it can be seen that for larger smoothing kernels, it becomes important to prevent signal leakage across kissing gyri.

74Additional descriptions of the smoothing algorithms and their respective application in LayNii is given here: https://layerfmri.com/anatomically-informed-spatial-smoothing/.

75Rigorous quality assessment of fMRI time series has become a fundamental pillar of any good code of conduct in fMRI research. In conventional fMRI, a consensus about good practices has been found and these practices are implemented in the major standardized analysis streams. In layer-fMRI, however, conventional QA metrics are not sufficient or can even be completely inapplicable. Examples of metrics insufficient for layer-fMRI QA include a.) tSNR, b.) motion displacement, or c.) variance explained of neural ICA components. tSNR: Since tSNR is calculated by means and the temporal standard deviation (STDEV), it cannot easily capture multiple sources of variance that are adding non-linearly (sum of squares). Since layer-fMRI is usually dominated by thermal noise, layer-fMRI specific artifacts are not well captured in tSNR maps. They are rather hidden in the thermal noise.Motion displacement: While it is critically important in layer-fMRI to manually check displacement estimates like in conventional fMRI, it is not sufficient. In layer-fMRI motion usually comes along with higher-order EPI-Phase inconsistencies too. Namely, even when the rigid head motion can be corrected for, motion-related distortion changes and motion related B0-changes cause intermittent artifacts.Variance explained in neural ICA components: While the relative variance explained by neural ICA (independent component analysis) components is a useful measure for conventional fMRI, it is limited in layer-fMRI applications. ICA works best, if Gaussian and non-Gaussian fluctuation sources can be separated from each other. Unfortunately however, in layer-fMRI, the neurally-driven temporal fMRI signal fluctuations are deeply hidden in the dominating Gaus-sian thermal noise and it is not uncommon that there isn’t a single neural ICA component selectable in layer-fMRI analyses.

76tSNR: Since tSNR is calculated by means and the temporal standard deviation (STDEV), it cannot easily capture multiple sources of variance that are adding non-linearly (sum of squares). Since layer-fMRI is usually dominated by thermal noise, layer-fMRI specific artifacts are not well captured in tSNR maps. They are rather hidden in the thermal noise.

77Motion displacement: While it is critically important in layer-fMRI to manually check displacement estimates like in conventional fMRI, it is not sufficient. In layer-fMRI motion usually comes along with higher-order EPI-Phase inconsistencies too. Namely, even when the rigid head motion can be corrected for, motion-related distortion changes and motion related B0-changes cause intermittent artifacts.

78Variance explained in neural ICA components: While the relative variance explained by neural ICA (independent component analysis) components is a useful measure for conventional fMRI, it is limited in layer-fMRI applications. ICA works best, if Gaussian and non-Gaussian fluctuation sources can be separated from each other. Unfortunately however, in layer-fMRI, the neurally-driven temporal fMRI signal fluctuations are deeply hidden in the dominating Gaus-sian thermal noise and it is not uncommon that there isn’t a single neural ICA component selectable in layer-fMRI analyses.

79Due to the limits discussed above, conventional QA metrics in layer-fMRI data, additional and adapted metrics are needed. Since most layer-fMRI data are acquired in the thermal noise-dominated regime, the noise can often be considered to be Gaussian. This means that non-thermal noise (e.g. artifact-induced signal clutter) can be identified by its non-Gaussian characteristics. Two useful measures to quantify how non-Gaussian a noise distribution is, are the measures of kurtosis and skewness. Fig. 7E,H depict a layer-fMRI example, where skew and kurtosis maps can depict artifact affected areas that are not clearly detectable in conventional tSNR maps.

80While thermal noise is not coupled across space and time, some forms of layer-fMRI artifacts result in strong noise coupling. This noise coupling can be used to identify artifacts and depict them despite the presence of dominating Gaussian noise. Informative measures that can be estimated on a voxel-wise level are: a) auto-correlation and b) the correlation with the global signal (GS). Fig. 7F,G depicts layer-fMRI examples where these measures can isolate and depict phase inconsistencies that are particularly limiting in layer-fMRI acquisition protocols.

81While the above QA metrics are very useful to characterize the time course quality of individual voxels, they cannot easily capture spatiotemporal interactions. E.g. if neighboring voxels share the same sources of noise (e.g. due to spatial signal leakage), the above mentioned QA metrics are ignorant to this. This is particularly problematic because the EPI readout for layer-fMR is usually very long and layer fMRI data are, thus, expected to suffer from unwanted T2*-related temporal noise coupling between neighboring voxels.

82To characterize the strengths and the pattern of the temporal correlation of neighboring voxels, a so-called noise kernel calculation is implemented in LayNii as described in Fig. 8A. The LayNii program estimates the average temporal correlation of neighboring voxels. In contrast to functional smoothness estimations in alternative software suites, the LayNii implementation of the noise kernel estimation does not simplify the noise coupling with individual Gaussian FWHM values and/or exponential decay terms. Instead, LayNii also writes out the 3-dimensional noise kernel as a volume nii file. This can be advantageous because a full noise kernel can capture negative side lobes, anticorrelations, as well as common saw-tooth patterns between odd and even lines. Furthermore, it also provides estimates along the diagonal directions, which can be helpful for 3D acquisition sequences (3D-EPI or 3D-GRASE). This form of three dimensional noise-correlation kernel can in the future be used for data-driven deconvolution-based image sharpening filters (e.g. see the ’-laurenzian’ flag in the LayNii program LN_DIRECT_SMOOTH). Furthermore, such analyses could be applied across layers and columns to estimate the topographical features of functional connectivity within the cortex (see Fig. S8 in).

83Additional descriptions of the QA algorithms, the limits of their interpretability and their respective application in LayNii is given here: https://layerfmri.com/qa/.

84The VASO contrast is believed to be inversely proportional to CBV changes and can non-invasively capture microvascular layer-specific signal responses through selective detection of signal changes in the extravascular volume compartment. The specific UHF-optimized Slab-selective slice inversion (SS-SI)-VASO approach acquires bloodnulled images and GE-BOLD images concomitantly. To directly provide a CBV and BOLD contrast time course comparable to conventional fMRI signal contrast analyses, the raw VASO time series need to be temporally resorted and orthogonalized into clean T2*-dependent and T1-dependent functional contrasts. The corresponding preprocessing steps are illustrated in Fig. 9A. To remove unwanted T2*-related signal in CBV-time series, LayNii assumes a simple magnetization model that obeys the Bloch equations (Fig. 9C). According to the Bloch equations, T1 and T2* decays behave like two independent relaxation effects that can be described in a multiplicative fashion. This means that unwanted T2* contrast in the CBV-weighted signal can be canceled out by means of dynamic division of MR images with and without preceding T1 -weighting (Fig. 9C).

85Fig. 9D-E depicts a representative example of VASO and BOLD data as they are generated from LayNii. Note that the VASO output is inversely correlated with dynamic CBV changes.

86For a hands-on example of how the application of LayNii works for layer-fMRI VASO applications see here: https://layerfmri.com/analysispipeline/.

87Even though high-resolution GE-BOLD is known to suffer from unwanted locally-nonspecific signals of large draining veins and the location of large draining veins on the cortical surface, due to its high sensitivity it is the most popular acquisition method for layer-fMRI data. Due to the blood drainage towards the cortical surface, the unwanted GE-BOLD signal of large veins is known to increase towards the cortical surface. Based on the angioarchitecture, it has been proposed that the unwanted layer-specific venous signal bias in GE-BOLD can be accounted for in retrospective model-based correction. A number of models are popular in the field (Fig. 10A): The linear-offset model is based on the assumption that the layer-specific microvascular signal is added on top of a task-independent macrovascular signal. Thus, it is assumed that a simple task-contrast subtraction or a simple linear-trend-removal would get rid of the unwanted macrovascular component of the BOLD signal. The linear offset model has been explicitly discussed in the context of linear correlation analyses and in the context of layer-dependent shape parameters. This linear model is furthermore implicitly used for in countless task contrast subtraction analyses of previous layer-fMRI studies.The scaling model is based on the assumption that the layer-specific bias in GE-BOLD is driven by layer-dependent variations of vein density. As such, it is assumed that the superficial signal is larger than the signal in deeper layers, simply because the amount of venous blood volume is higher in superficial layers compared to deeper layers. Or in other words, the macrovasculature acts like a layer-specific signal amplification (gain). Thus, it is assumed that a simple multiplicative (or divisory) normalization can correct for the macrovascular signal bias. This layer-dependent signal normalization has been proposed as part of several layer-fMRI analysis strategies, including: scaling the layer-dependent signal with estimates of layer-dependent venous CBV,normalizing the signal difference between different task responses by the mean signal response of all involved task responses,refraining to infer neuroscience conclusions solely based on task response differences in favor of focusing on task response ratios that are normalized by the presumably vascularly driven signal fluctuations.The leakage model is based on the assumption that the BOLD signal in each layer constitutes a signal mixture of activity originating in multiple layers. Namely, it is assumed that a voxel in a given layer contains some signal from the layer of interest plus an additional unwanted integrated signal of all the layers beneath. E.g. it is assumed that the BOLD signal in a voxel of a superficial layer contains both the desired superficial signal plus unwanted signal from the deeper layers. The BOLD signal model in each layer is then usually parameterized as a weighted sum of the layer itself and the signal from deeper layers. To use this model to correct for unwanted leaked signals, a spatial signal deconvolution approach is applied. (Note that some implementations of such deconvolution models refrain from the term “deconvolution” in favor of the terms “matrix inversion” or “consecutive signal subtraction”). Very often, the leakage model is combined with the scaling model. Layer-fMRI focused implementations of BOLD signal models that assume spatial signal leakage across layers have been described by a large number of research labs, including models from,,,,, and. All of these signal leakage models estimate the signal in superficial layers as a sum of the microvascular response within the given layer plus a weighted sum of the signal from all other layers below. The main difference between the various leakage models comes from the procedure by which the respective summation weights are estimated. The model in describes the venous signal draining with travelling wave equations. The models described in perform forward simulations without the explicit aim to invert the model for venous signal removal. The models described in, propose such a model inversion for the removal of unwanted venous signal in layer-fMRI signal processing. The newer model in entails more physiologically-informed constraints than the model in. In most of those models, the summation weights of layer-dependent signal leakage are derived based on more-or-less appropriate assumptions of the layer-dependent blood vessel architecture and/or experimental data. The weights usually depend on resting CBVv and CBF assumptions across cortical depth.

88The linear-offset model is based on the assumption that the layer-specific microvascular signal is added on top of a task-independent macrovascular signal. Thus, it is assumed that a simple task-contrast subtraction or a simple linear-trend-removal would get rid of the unwanted macrovascular component of the BOLD signal. The linear offset model has been explicitly discussed in the context of linear correlation analyses and in the context of layer-dependent shape parameters. This linear model is furthermore implicitly used for in countless task contrast subtraction analyses of previous layer-fMRI studies.

89The scaling model is based on the assumption that the layer-specific bias in GE-BOLD is driven by layer-dependent variations of vein density. As such, it is assumed that the superficial signal is larger than the signal in deeper layers, simply because the amount of venous blood volume is higher in superficial layers compared to deeper layers. Or in other words, the macrovasculature acts like a layer-specific signal amplification (gain). Thus, it is assumed that a simple multiplicative (or divisory) normalization can correct for the macrovascular signal bias. This layer-dependent signal normalization has been proposed as part of several layer-fMRI analysis strategies, including: scaling the layer-dependent signal with estimates of layer-dependent venous CBV,normalizing the signal difference between different task responses by the mean signal response of all involved task responses,refraining to infer neuroscience conclusions solely based on task response differences in favor of focusing on task response ratios that are normalized by the presumably vascularly driven signal fluctuations.

90scaling the layer-dependent signal with estimates of layer-dependent venous CBV,

91normalizing the signal difference between different task responses by the mean signal response of all involved task responses,

92refraining to infer neuroscience conclusions solely based on task response differences in favor of focusing on task response ratios that are normalized by the presumably vascularly driven signal fluctuations.

93The leakage model is based on the assumption that the BOLD signal in each layer constitutes a signal mixture of activity originating in multiple layers. Namely, it is assumed that a voxel in a given layer contains some signal from the layer of interest plus an additional unwanted integrated signal of all the layers beneath. E.g. it is assumed that the BOLD signal in a voxel of a superficial layer contains both the desired superficial signal plus unwanted signal from the deeper layers. The BOLD signal model in each layer is then usually parameterized as a weighted sum of the layer itself and the signal from deeper layers. To use this model to correct for unwanted leaked signals, a spatial signal deconvolution approach is applied. (Note that some implementations of such deconvolution models refrain from the term “deconvolution” in favor of the terms “matrix inversion” or “consecutive signal subtraction”). Very often, the leakage model is combined with the scaling model. Layer-fMRI focused implementations of BOLD signal models that assume spatial signal leakage across layers have been described by a large number of research labs, including models from,,,,, and. All of these signal leakage models estimate the signal in superficial layers as a sum of the microvascular response within the given layer plus a weighted sum of the signal from all other layers below. The main difference between the various leakage models comes from the procedure by which the respective summation weights are estimated. The model in describes the venous signal draining with travelling wave equations. The models described in perform forward simulations without the explicit aim to invert the model for venous signal removal. The models described in, propose such a model inversion for the removal of unwanted venous signal in layer-fMRI signal processing. The newer model in entails more physiologically-informed constraints than the model in. In most of those models, the summation weights of layer-dependent signal leakage are derived based on more-or-less appropriate assumptions of the layer-dependent blood vessel architecture and/or experimental data. The weights usually depend on resting CBVv and CBF assumptions across cortical depth.

94In current layer-dependent applications that use model-based deveining, the models are applied area-wide (e.g.). I.e., fMRI signals are pooled from the layers of large cortical patches and their effect is solely investigated in the form of one-dimensional layer profiles. In LayNii, however, the various vein mitigation algorithms work on a voxel-by-voxel level. This allows the researcher to appreciate the signal distribution along the topographical space across layers and columns, simultaneously. Furthermore, the resulting signal maps provide an intuitive understanding of the noise amplification that comes along with the various deveining algorithms of the respective models (See Fig. 10B for examples).

95All three model-based deveining methods are included in LayNii. The linear offset model is used by executing LN2_DEVEIN with the flag -linear. The CBV scaling model is used by executing LN2_DEVEIN with the flag -CBV. The leakage model is used if LN2_DEVEIN is executed without either of those flags. The formalism of implementation of the leakage model is largely based on. Despite the different approach to estimate leakage weights, the formalism implemented in LayNii that describes the signal leakage as a spatial weighted signal summation is identical for all published leakage models mentioned above. Since CBVv can be obtained from the GLM residuals, no additional measurements are needed. Additional descriptions of the model-based deveining algorithms and their respective application in LayNii is given here: https://layerfmri.com/devein/. Note that the model-based deveining described here should not be confused with the auxiliary method of (semi-)automatic detection of large veins and masking of vein-dominated voxels.

96While the authors are excited about the emergence of model-based deconvolution methods in many research labs, these model-based vein mitigation methods are just starting to be validated, their shortcomings are just starting to be understood, and their applications are not widely established except for exploratory use. Currently, it is still more established in the field to avoid unwanted venous signals by means of advanced acquisition strategies, rather than adding additional signal analysis (deconvolution) steps. Specific extensions to current models that are currently being investigated are: So-called blooming effects (cross-voxel dipole fields) of very large draining veins into neighboring voxels and layers that are currently not accounted for in model-based deveining methods.More detailed knowledge of the signal distribution between intra- and extravascular BOLD effect across field strengths and echo times for various readout trajectories during T2*-decay needs to be obtained, validated, and incorporated in future models. Due to the sheer numbers of model parameters and experimental setups, elaborate and comprehensive validation studies are needed.Orientation effects of pial and ascending veins with respect to the external magnetic field that can lead to over-estimation or underestimation of the contribution in model-based deveinig, if not accounted for.A crucial requirement of model-based deveining is the assumption of a good baseline signal. Thus, impractical inter-stimulus-intervals and task-duration dependencies might hamper straightforward application of model-based deveining in neuroscience studies.

97So-called blooming effects (cross-voxel dipole fields) of very large draining veins into neighboring voxels and layers that are currently not accounted for in model-based deveining methods.

98More detailed knowledge of the signal distribution between intra- and extravascular BOLD effect across field strengths and echo times for various readout trajectories during T2*-decay needs to be obtained, validated, and incorporated in future models. Due to the sheer numbers of model parameters and experimental setups, elaborate and comprehensive validation studies are needed.

99Orientation effects of pial and ascending veins with respect to the external magnetic field that can lead to over-estimation or underestimation of the contribution in model-based deveinig, if not accounted for.

100A crucial requirement of model-based deveining is the assumption of a good baseline signal. Thus, impractical inter-stimulus-intervals and task-duration dependencies might hamper straightforward application of model-based deveining in neuroscience studies.

101Here, we present an fMRI analysis toolbox that is specifically developed for layer-fMRI applications. It has minimal requirements on the input data and works in the native data format of voxels as they come off the scanner. Its application ranges from layer/column assignment, over anatomically informed smoothing and QA, up to model-based deveining and VASO analyses.

102LayNii is not designed to duplicate analysis features that have already been implemented in alternative fMRI software packages with the same capabilities. Even though there are specific layer-fMRI demands for partial coverage image alignment, LayNii does not perform it itself. Since LayNii also works in the native EPI-space, image alignment is not absolutely necessary to begin with. If a layer-fMRI user is interested in alignment constraints that are specific to layer-fMRI data, we refer to the instructions of how to perform image alignment of layer-fMRI data with ITK-snap and ANTs (https://layerfmri.com/high-quality-registration/), or with AFNI (https://blog.cogneurostats.com/aligning-partial-volume-fmri-in-afni/) and. Furthermore, even though tissue-type segmentation is particularly important in layer-fMRI analyses, LayNii does not contain its own segmentation tools. In layer-fMRI, GM borders often need to be manually adjusted and corrected anyway. Even though this is a time consuming and tedious task of several work days, in layer-fMRI, it is often most efficient to fully manually segment the areas of interest within the limited EPI coverage without the necessity of automated segmentation tools. For large coverages and/or a rough first estimate of GM borders as an input to LayNii, alternative software packages can be used. E.g. the interested reader is referred to an instruction of how to use FreeSurfer, nighres and SUMA to generate an input rim file for LayNii here: (https://layerfmri.com/getting-layers-in-epi-space/), this segmentation can be further corrected with the semi-manual segmentation tool Segmentator (https://github.com/ofgulban/segmentator). The LayNii program LN_RIMIFY can be used to convert the segmentation output of common third-party software tools to a LayNii-readable rim file.

103A collection of non-LayNii scripts and programs that are designed to work in tandem with LayNii available to download here: https://github.com/ofgulban/LayNii_extras.

104While it is a dedicated aim of the developers to include it in the LayNii software suite, the current version of LayNii does not contain a program that estimates the geometric angle between the cortical surface and the voxel grid or the main magnetic field. (Note by author during the article production stage: cortial orientation vectors have been included into LN2_LAYERS with the option -streamline in the LayNii release v2.0).

105Here, ‘layer-fMRI’ is used as an umbrella term for depth-dependent fMRI, intra-cortical fMRI, and sub-millimeter fMRI in general1. This means that LayNii is also explicitly intended for signal analyses of columnar structures (e.g. see Fig. 6)). LayNii is originally intended for -but not limited to- functional MRI. In fact, LayNii is applicable and has been applied to high-resolution structural MRI, as well as histology data with and without restricted field of views (E.g. Fig. 3 and). We believe that LayNii might also be specifically suited for in-vivo sub-millimeter diffusion-weighted data. While LayNii is time-efficiently applicable to conventional high-resolution whole brain anatomical data, it is not specifically advertised for this purpose. While LayNii can be applied on conventional structural data with distortion-corrected and aligned functional data, there are alternative analysis software packages that are more extensively tested, documented and supported with exhaustive educational tutorials for this purpose.

106While LayNii is specifically designed for high-resolution applications, individual programs can be applied to lower spatial resolution fMRI data too. Namely, the VASO-specific analyses can -and have been-applied to lower-resolution 3 mm data too. Similarly, the quality assessment metrics can also be beneficial at lower resolutions.

107Since LayNii works directly in voxel space without imposing topological constraints, it can straightforwardly generate layers and column estimates in 2D data. This allows applications in common single slice data of preclinical MRI in rodents, cats and monkeys, and it allows layer profile extraction of figures from any electronic publication including the seminal microscopy images of an entire century of cytoarchitecture, myelo-architecture and angioarchitecture research.

108Note, however, that in 2D, the equi-volume layering approach converges with an equi-area approach. For more information on how to generate layer-profiles of literature figures see: https://layerfmri.com/2dlayers/.

109The LayNii software suite has been used in multiple highresolution fMRI studies across the world. It’s Github site welcomes >1400 unique visitors per year and its code is being cloned >100 times per month. It found particular application for layer-fMRI studies in the motor cortex, in sensory cortex, DLPFC, across association cortices for columnar imaging in the motor cortex, in layer-specific functional connectivity mapping, for mental imaginary layer-fMRI, for methods development of new sequences, for visual layer-fMRI, for model-based removal of vein effects in layer-fMRI-EEG and for methods debugging of human 9.4T layer-fMRI. In the early days of LayNii, its programs were shaped and optimized by continuous interactions with its users. Current and future requests for new features and support are encouraged via LayNiis github issue page.

110All the discussed LayNii programs (BOLD correction, equidistance and equi-volume layering, etc.) have been initially implemented in the ODIN (Object-oriented Development Interface for NMR) environment. This LayNii predecessor software suite was originally used in many initial layer-fMRI studies. However, since this software package had more than 30 dependencies, and required files to be located in certain root-locations, it was not straight-forward to use across various operating platforms, nor was it installable on servers with conservative user access. Thus, with guidance from the NIH data science and sharing team (https://cmn.nimh.nih.gov/) and the AFNI team (https://afni.nimh.nih.gov/), the layer-fMRI algorithms were reimplemented outside of ODIN with a data I/O that is inherited from the original NIFTI release (https://nifti.nimh.nih.gov/). Thus, LayNii does not have any residual external dependencies on third-party software, nor does it rely on libraries.

111Here we introduce a new dedicated software toolbox for layer-specific (functional) MRI. It is open source and available for installation via source code or pre-compiled binaries for Linux, Windows and macOS and provides a comprehensive set of advanced techniques for layer-fMRI analyses. The current functionality largely focuses on applications where the challenges of layer-fMRI data do not allow the application of standard analysis pipelines of the major software packages. Though, we hope that the modular package structure facilitated augmentation of analysis pipelines that also contain other software packages.