Layer-dependent functional connectivity methods

Recent methodological advances in fMRI contrast and readout strategies have allowed researchers to approach the mesoscopic spatial regime of cortical layers. This has revolutionized the ability to map cortical information processing within and across brain systems. However, until recently, most layer-fMRI studies have been conﬁned to primary cortices using basic block-design tasks and macro-vascular-contaminated sequence contrasts. To become an established method for user-friendly applicability in neuroscience practice, layer-fMRI acquisition and analysis methods need to be extended to more ﬂexible connectivity-based experiment designs; they must be able to capture subtle changes in brain networks of higher-order cognitive areas, and they should not be spatially biased with unwanted vein signals. In this article, we review the most pressing challenges of layer-dependent fMRI for large-scale neuroscientiﬁc applicability and describe recently developed acquisition methodologies that can resolve them. In doing so, we review technical tradeoffs and capabilities of modern MR-sequence approaches to achieve mea-surements that are free of locally unspeciﬁc vein signal, with whole-brain coverage, sub-second sampling, very high resolutions, and with a combination of those capabilities. The presented approaches provide whole-brain layer-dependent connectivity data that open a new window to investigate brain network connections. We exemplify this by reviewing a number of candidate tools for connectivity analyses that will allow future studies to address new questions in network neuroscience. The considered network analysis tools include: hierarchy mapping, directional connectomics, source-speciﬁc connectivity mapping, and network sub-compartmentalization. We conclude: Whole-brain layer-fMRI without large-vessel contamination is applicable for human neuroscience and opens the door to investigate biological mechanisms behind any number of psychological and psychiatric phenomena, such as selective attention, hallucinations and delusions, and even conscious-ness itself.


Abstract
Recent methodological advances in fMRI contrast and readout strategies have allowed researchers to approach the mesoscopic spatial regime of cortical layers. This has revolutionized the ability to map cortical information processing within and across brain systems. However, until recently, most layer-fMRI studies have been confined to primary cortices using basic block-design tasks and macro-vascular-contaminated sequence contrasts. To become an established method for user-friendly applicability in neuroscience practice, layer-fMRI acquisition and analysis methods need to be extended to more flexible connectivity-based experiment designs; they must be able to capture subtle changes in brain networks of higher-order cognitive areas, and they should not be spatially biased with unwanted vein signals. In this article, we review the most pressing challenges of layerdependent fMRI for large-scale neuroscientific applicability and describe recently developed acquisition methodologies that can resolve them. In doing so, we review technical tradeoffs and capabilities of modern MR-sequence approaches to achieve measurements that are free of locally unspecific vein signal, with whole-brain coverage, sub-second sampling, very high resolutions, and with a combination of those capabilities. The presented approaches provide whole-brain layer-dependent connectivity data that open a new window to investigate brain network connections. We exemplify this by reviewing a number of candidate tools for connectivity analyses that will allow future studies to address new questions in network neuroscience. The considered network analysis tools include: hierarchy mapping, directional connectomics, source-specific connectivity mapping, and network sub-compartmentalization. We conclude: Whole-brain layer-fMRI without large-vessel contamination is applicable for human neuroscience and opens the door to investigate biological mechanisms behind any number of psychological and psychiatric phenomena, such as selective attention, hallucinations and delusions, and even consciousness itself.

Introduction
Methodological advancements of functional Magnetic Resonance Imaging (fMRI) contrasts and readout strategies in recent years have allowed researchers to approach the mesoscopic spatial regime of cortical layers Huber et al. 2017;Kok et al. 2016;Koopmans et al. 2010;Polimeni et al. 2010b). This has allowed investigations of the connectivity and function in neural microcircuits across cortical layers within cortical areas (Callaway 1998;Lund 1988). Non-human primate studies have shown that in hierarchically organized brain systems, inter-area neural input arrives in different layers for bottom-up or top-down connections. Specifically, feed-forward, bottom-up activity (e.g. V1 → V5/hMT) terminates predominantly in the middle granular layer, while feedback, top-down activity (e.g. V1 ← V5/hMT) terminates predominantly in superficial and/or deeper layers (Angelucci et al. 2002;Felleman and Van Essen 1991). Thus, the ability to map brain activity across cortical layers revolutionizes the ability to tackle cortical information processing within brain systems and might become an important tool for cognitive network neuroscience (Lawrence et al. 2019;Stephan et al. 2019). Until now, however, most human layer-fMRI studies have been confined to primary motor and sensory cortex and have used basic fMRI task block-designs. There are several reasons why layer-fMRI methods are not applied more widely. One of the major limitations is that large draining veins in conventional GE-BOLD fMRI methods spatially blur the fMRI activity, imposing a cortical depth specific contrast weighting across the layers -with superficial layers being most heavily weighted (Uludag and Blinder 2018). Thus, locally specific interpretations of individual layer activity are not straightforward with GE-BOLD (Kay et al. 2019). Alternative non-BOLD fMRI sequences (reviewed in Huber et al. (2019)), such as the CBV-weighted VASO (Lu et al. 2003) method, have been shown to be locally more specific and having contrast that is more evenly weighted across the cortical depths (Huber et al. 2015). Until now, addressing layer-specific neuroscientific research questions with any layer-fMRI imaging tool (BOLD and non-BOLD) has been limited by: A) low sensitivity, B) asymmetric voxel sizes, C) restricted brain coverage, D) unclear potential for network analyses.
A) Low sensitivity: The reduced sensitivity of the early initial layer-fMRI methods limited its applicability to strong block-design tasks that evoked massive brain activations in primary brain areas, e.g. strong visual tasks to evoke activity in primary visual cortex (Kashyap et al. 2018;Kok et al. 2016;Koopmans et al. 2010;Marquardt et al. 2018;Polimeni et al. 2010b), strong auditory tasks to evoke activity in primary auditory cortex , and strong sensorimotor tasks to engage the primary somatosensory and primary motor cortex (Huber et al. 2015;Yu et al. 2019). If these methods are to contribute to our understanding of cognitive neuroscience, studies of higher-order regions with more sophisticated tasks and task designs will be necessary.
B) Asymmetric voxel size: Increasing the sensitivity of the acquisition methods by means of using asymmetric (non-isotropic) voxel sizes has allowed more advanced investigations, however it has limited the addressable neuroscience questions to small patches of cortex (Finn et al. 2019;Guidi et al. 2016;Huber et al. 2015;Kashyap et al. 2018;Koopmans et al. 2010;Menon and Goodyear 1999;Yu et al. 2019). Particularly in brain areas with stable folding patterns across many individual people (e.g. primary motor cortex), the approach of thicker MRI-slices can help improve sensitivity. However, to fully grasp the neural representation of laminar activation across the folded cortical ribbon in higher-order brain areas with variable folding patterns, isotropic submillimeter resolutions are vital.
C) Small coverage: Restricted field-of-view imaging allows investigations of high-resolution brain connections in individual brain areas with great reliability (reviewed in Schluppeck et al. (2018)). The reduced coverage allows fast data sampling rates and short scan times. However, the slice positioning and protocol setup at the scanner is more difficult for the experimenter. And furthermore, it does not take advantage of the ultimate potential in layer-fMRI to apply network analysis tools for investigating the information processing and the interplay within and across brain areas.
D) Unclear potential for network analyses: Until now, the field of layer-fMRI has been mostly focusing on acquisition and analysis methodologies. Thus, other aspects of the experimental design have been kept as conservative as possible. This means that in virtually all layer-fMRI studies published so far, the investigators focused on 1-3 specific predefined brain areas of healthy participants and investigated the layerdependent activity changes in these areas with a small number of discrete task conditions. The potential of less conservative experimental setups and areaindependent analysis methods in layer-fMRI is unclear and has not been systematically applied so far. Thus, layer-fMRI has not yet become an established tool for use in whole-brain connectivity studies, which are of increasing importance to human neuroscience research.
To increase the utility of VASO-based layer-fMRI acquisition and analysis methods for the application in human cognitive neuroscience, all of the above limitations must be addressed. Only, then these methods can be used in more sophisticated task designs, in higher-order brain areas, and while capturing subtle changes in brain networks. In this special issue article, we argue that with the methodological advancements in the past years, the aforementioned limitations are largely resolved. Layer-fMRI connectivity measures are now available without venous contamination and with whole-brain coverage. This paves the road for a new class of neuroscience questions to be investigated. We aim to exemplify how layerdependent whole-brain fMRI datasets might be used for a novel family of connectivity analyses by discussing multiple examples of layer-dependent connectivity analysis procedures that will shed new light into the internal working principles and functional organization of the human brain 1 . In this paper, we show several examples of layer-specific fMRI.
1 Note to the reviewers: As you can see, this manuscript contains aspects of a review article, as well as some new exemplary data. General layer-fMRI methodology is reviewed by means of previously unpublished example data. This manuscript format is chosen based on the special issue invitation and after consultation with the editors. Furthermore, this manuscript is on purpose right at the border between an acquisition methods paper and analysis paper for neuroscientific applications. Therefore, we are starting with a methods section that describes in detail the approaches taken.

Methods
Data from N = 31 participants were used for the examples shown in this paper participating in a movie watching study (N = 12), in a visual study (N = 12) and in a sequence development study (N = 5). All experiments were conducted in accordance with the Belmont Report after participants granted informed consent. Data were acquired in FMRIF/NIMH/NIMH in Bethesda, USA and at Scannexus in Maastricht, The Netherlands. Experiments at NIH were conducted according to the US Federal Regulations that protect human subjects and approved by the Combined Neuroscience Institutional Review Board protocol 93-M-0170 (ClinicalTrials.gov identifier: NCT00001360). All experiments of this study were conducted on SIEMENS 7T scanners (Siemens Healthineers, Erlangen, Germany) using the vendor-provided IDEA environment (VB17A/VE11U/VE12K-UHF). The scanners were equipped with the standard body gradient coils (maximum effective gradient strength used here: 48 mT/m; maximum slew rate used: 198 T/m/s). For RF transmission and reception, volume-transmit and 32-channel receive head coils (Nova Medical, Wilmington, MA, USA) were used. All functional data were acquired with VASO contrast (Hua et al. 2013;Huber et al. 2014b;Lu et al. 2003) and a 3D-EPI readout (Poser et al. 2010). Slab-oversampling was applied in the second phase direction. 3D slice aliasing was minimized by using a sharp slab-excitation pulse profile with a bandwidth-time-product of 25. The readout parameters were: in-plane FLASH GRAPPA 3 (Talagala et al. 2016), with regularization strengths chi = 1-5000 (for instructions how to do this, see here: https://layerfmri.com/grapparegularization/), TE = 25 ms, partial Fourier = 6/8 with POCS reconstruction of 8 iterations. Online reconstruction was performed using a GRAPPA kernel size ×6 × 5 × 3. Resolutions of fMRI protocols varies between 0.5 × 0.5 × 0.6 and 0.8 × 0.8 × 0.8 mm 3 . For anatomical reference, 3 averages of 0.5 mm iso MP2RAGE data with partial brain coverage were acquired. For optimal use in layer-fMRI (https://layerfmri.com/mp2rage), background noise (outside the brain) was minimized with a regularization term lambda = 0.1, fold-over-artifacts were minimized by using an excitation pulse BWTP of 25, and no partial Fourier was applied in the second phase encoding direction. Physiological traces of respiration and heartbeat were recorded for the first 12 participants that underwent the visual task. Further recordings were discontinued because in the thermal noisedominated regime of 0.8mm voxel size here, we did not observe any improvements by conducting RETROICOR (Li et al. 2000) physiological noise correction ). All fMRI images were motion corrected using SPM12 (UCL, UK) (Penny et al. 2007). Volume realignment and interpolation were performed with a 4th-order spline. In order to minimize effects of variable distortions (non-rigid motion), the motion was estimated in a manually drawn mask. The 8 outermost slices were excluded from the analysis to minimize any residual 3D-EPI related slab fold-over artifacts and reslicing artifacts, when necessary. Raw VASO time series usually consist of interleaved acquisitions of MR signal with and without blood nulling. Thus, the analysis pipeline treats odd and even time points separately for the motion correction and subsequently divides the images by each other in order to correct for BOLD contaminations in VASO (Huber et al. 2014b). Boundary lines of the gray matter (GM) ribbon to cerebrospinal fluid (CSF) and white matter (WM) were estimated with FreeSurfer (v6.0) (Fischl 2012) from wholebrain 0.7mm MP2RAGE data that had been acquired in previous studies and were registered here to the functional EPI space using ANTS (Avants et al. 2008). These boundaries were used as visual guidance to manually segment GM in EPI space slice-by-slice, as described in the manual ( https://layerfmri.com/getting-layersin-epi-space). A coordinate system across cortical layers and columnar structures was estimated in LAYNII (https://github.com/layerfMRI/LAYNII). LAYNII is an open source C++ software suite for computing layer and column functions. We estimated the depth of equidistantly distributed layers 2 . With the resolution of 0.8 mm, we obtained 4-6 independent data points across the thickness of the cortex. Across these data points, we created 20 layers 3 across the thickness of the cortex on a 4-fold finer grid than the effective resolution. The entire pipeline of this analysis is described here: https://layerfmri.com/analysispipeline. Cortical 'columnar' structures were also determined in LAYNII's LN_3DCOLUMN program with the following algorithm: For each voxel, the next closest (Euclidean distance) voxel in each other layer is determined. This group of voxels is considered a 'columnar' structure. This means that every voxel has a unique column assignment. For more information on the algorithm see here: https://layerfmri.com/columns. Note that the terminology of 'columns' refers to the geometric columnar shape of fMRI voxel groups. It does not necessarily refer to groups of neurons with the identical receptive field (a.k.a. hypercolumn, macro column, functional column, cortical module).

MR-sequences for layer-activity CBV mapping with large coverage
Previous layer-fMRI CBV-weighted sequences could not capture more than 5-18 slices resulting in small field of views (6 -12 cm), even when using anisotropic voxels (Guidi et al. 2016;Huber et al. 2020b;2018a;. Recent developments in VASO readout strategies now allow substantial coverage improvements to capture functional connectivity in larger brain networks. In the following sub-sections, we describe these approaches, while also mentioning their caveats: 2.1.1. SS-SI-VASO covering large brain networks. The CBV weighting in VASO is based on T 1 -selective nulling of the blood signal, while keeping the extravascular signal with its different T 1 preserved for image creation (Hua et al. 2013;Lu et al. 2003). Thus, VASO is inherently a single-slice method, as multi-slice acquisitions would suffer from inconsistent or insufficient blood nulling across slices. With 3Dreadout methodologies, the readout window can be substantially increased up to a time window of 0.5 -1 seconds ( T 1 ) without compromising blood nulling across the entire imaging volume (Hua et al. 2013;Poser and Norris 2009;. Extending the readout acquisition window even further by another 1 -4 seconds (≈ T 1 ) is technically possible, however, it comes along with some challenges. Namely, the longer TRs mean that the timing of the inversion-recovery relaxation is faster than the k-space acquisition. Compared to previous VASO sequences, these limitations have been recently mitigated with a series of new approaches: 1) Dead-times between VASO and BOLD images have been removed. This increases the imaging duty cycle of the sequence, at the cost of an asymmetric readout. This means that, in post-processing, an additional temporal interpolation step needs to be included.
2) More efficient variable-flip angle regimes were designed to manipulate the signal strength along the second phase-encoding direction in 3D-EPI. It was proposed using an exponentially increasing flip angle for the first half of k-space and keeping a relatively high flip angle constant for the second half of k-space. The flip angle increase is chosen to match the T 1 -decay of GM. The level of asymmetry is chosen such that the negative side lobes of the complex point-spread function are compensated for with the level of overflipping in the second half of k-space. While this approach allows substantial increases in the readout duration, it can come along with edge-enhancement artifacts in other T 1 -compartments, such as at WM/CSF borders.
3) The readout duration can be extended such that only the k-space center is acquired at the blood nulling time. This means that the outer k-space segments will have a contribution of residual intravascular signal, which can result in CBF contributions for the high spatial frequencies. Since CBF is believed to be dominated from capillary water exchange only, this will not compromise the layer-specificity of the sequence. However, this form of VASO signal amplification comes along with a reduced quantifiability of determining CBV in units of ml per 100 ml of tissue.
These optimizations allow increases in coverage up to 28 slices (2 slices oversampling) with isotropic 0.8 mm resolution and a sampling rate of 1.7s and 2.8s. See Fig. 1A for a graphical depiction of the sequence. Note that as long as there is a reference image acquired without the inversion pulse, but with the identical readout, the T 1 -related signal changes can be separated from T * 2 -related signal changes to ultimately extract CBV signal traces without BOLD contaminations. This BOLD correction is based on a dynamic division analysis under the assumption that, at 7T, the GE-BOLD effect is dominated by extravascular T * 2changes, as previously validated in (Huber et al. 2014a).

MAGEC-VASO for whole brain CBV imaging.
While VASO was originally proposed as a blood-nulling method (Lu et al. 2003), it has in the last 15 years been generalized to a general T 1 -contrast without specific blood-nulling requirements. Early VASO versions without blood nulling used a T 1 -selective GM-nulling procedure to estimate an inverse VASO contrast (Shen et al. 2009;Wu et al. 2008). Later on, the VASO formalism was further generalized to extract CBV changes at any inversion time (Ciris et al. 2014;W.C. et al. 2007). This literature has shown that the experimental trick of blood-nulling is not the only way of obtaining a CBV-weighting. In fact, as long as there is a different T 1 -weighting between the extravascular signal and intravascular signal, any volume redistribution between these pools of z-magnetization, will result in a VASO signal change. Thus, instead of using an inversion pulse, T 1 -weighting can also be introduced by variable flip angles that create a dynamic steady-state across k-space segments along the 3D-EPI trajectory. This approach has the advantage that the T 1weighting can be maintained in a dynamic equilibrium for as long as needed. Since the blood z-magnetization is not completely nulled here, the MAGEC approach can contain small CBF-dependent VASO signal amplifications. Since CBF is believed to be dominated from capillary water exchange only, this will not compromise the layer-specificity. It rather improves the sensitivity. As long as there is a reference image acquired without this T 1 -weighting, the T 1 -related signal changes can be separated from T * 2 -related signal changes to ultimately extract CBV signal traces without BOLD contaminations. Analogously to the MAGIC VASO method (Lu et al. 2004) with multiple inversion pulses, the variable flip angle approach is called Multiple Acquisitions with Global Excitation Cycling (MAGEC) VASO. Since MAGEC VASO does not rely on a given inversion time, the readout can be prolonged as much as needed (at the cost of TR). This allows for increased coverage with up to 72-104 slices (and more) at 0.8 mm isotropic resolution and with TR of 6.5 -8s. See Fig. 1A for a graphical depiction of the MAGEC sequence discussed here. For further discussions of the MAGEC sequence, see https://layerfmri.page.link/WholeBrain_ISMRM20.

In-plane segmented 3D-EPI.
Using conventional 7T systems with 32-channel receive arrays and body gradients, the effective resolution of single-shot EPI is encoding-limited by the readout speed relative to the T * 2 -decay rate. Given the typical size of a human head, and the constraints set by peripheral nerve stimulation, this limits the resolution of single-shot EPI to approximately 0.8mm. Segmented EPI approaches can overcome this limitation at the cost of longer TRs (Jin and Kim 2008; Menon and Goodyear 1999) (which, however, might limit some neuroscience applications whose task design requires short TRs). The correspondingly higher achievable spatial resolution can be used to provide better understanding of the underlying localization specificity and contrast-generating mechanisms. Here, we discuss such a segmented 3D-EPI approach (Stirnberg et al. 2016; to exemplify the laminar and columnar specificity and vessel sensitivity in conventional GE-BOLD and VASO. A matrix of 316 × 316 × 26 was used with 2 shots per kz-segment (TR = 3.9s).

Tasks to illustrate functional connectivity methods
Three different functional task classes are discussed here.
1) Method validation tasks: For validation and illustration of recent fMRI sequence developments, we used basic visual and sensorimotor tasks. In visual experiments, we used a block-design moving star field task (Huk et al. 2002) alternating between 30s rest, 30 sec static stars and 30s swirly stars. This task was chosen because it evokes robust and strong signal changes in early visual areas (including V1 and hMT+) across participants. Participants were instructed to fixate on the center of the screen. Each run lasted 12 min and was repeated 2-4 times.
In sensorimotor experiments, we used a block-design finger-tapping task. Participants were instructed to mimic a video of a moving hand performing this movement. Tapping and rest periods were alternated every 30s. Each run lasted 12min and was repeated 1-2 times (totaling 2-3). Such basic visual and sensorimotor tasks are established in the field of fMRI method developments as testbeds to validate the performance of new fMRI sequences.
2) Resting-state tasks: In order to illustrate the schematics of layer-fMRI functional connectivity analyses, we conducted resting-state experiments. During these experiments, we instructed the participants to keep their heads still and not fall asleep. All participants, that underwent the visual stimulation tasks, also participated in two 12 min resting-state experiments (in the same scan session), while imaging the same part of the brain.
3) Movie watching: In order to illustrate the applicability of layer-fMRI acquisition methods with 'naturalistic tasks', we conducted 12 sessions of movie watching. It has been suggested that movie watching is beneficial over alternative block-designed or resting-state tasks for various reasons: a) Movie watching can engage multiple brain systems at the same time and thus allows the researcher to extract more information in shorter time (Huth et al. 2012).
b) Participants move less during movie watching compared to resting-state (Vanderwal et al. 2015).
c) Watching the same movie multiple times results in highly reproducible time-locked brain activity patterns (Mandelkow et al. 2017). Thus, movie watching tasks allows the experimenter to repeat the data acquisition. This can be highly beneficial in layer-fMRI, e.g. to average datasets before conducting functional connectivity analyses, or to consecutively acquire different portions of the brain and then retrospectively combine the data for synthetic whole-brain connectivity analyses. d) Movie watching can be advantageous to capture individual differences compared to resting-state (Vanderwal et al. 2017).
e) Movie watching-induced activity is expected to be time-locked across participants (Kauppi et al. 2010). Thus, inter-participant correlation analyses can be performed to extract functional connectivity results without biases of participantspecific global physiological noise.
Here, we used an already established collection of 5 short video clips (https://layerfmri.page.link/7Tmovie). This collection of video clips has been used in the 7T HCP project (https://www.humanconnectome.org/study/hcpyoung-adult/article/first-release-of-7t-mr-image-data). More information about the movie task can be found at the HCP documentation under this link.

Proposed analysis algorithms to investigate directional connectivity in layer-dependent connectivity data
To show the type of new information about functional connectivity that can be extracted with layer-fMRI acquisition methods, we discuss a series of proposed analysis algorithms. Fig. 3). To illustrate the data quality and sensitivity of whole-brain layer-fMRI VASO during movie watching tasks, we extracted functional activation time courses of the most common functional networks with the following inter-subject correlation procedure: a) We used 98 datasets (not including siblings) from the 7T HCP database. These data were acquired with GE-BOLD at 1.6mm with whole brain coverage.

Network extraction with naturalistic tasks (in relation to
b) Signal traces were extracted from the functional network masks provided in Smith et al. (2009) and averaged across participants.
c) These data were temporally downsampled to match the layer-fMRI TR.
d) The resulting signal traces were then used as regressors (design matrix) in a conventional GLM-analysis in FSL-FEAT ).

2.3.2.
ROIs-based directional functional connectivity (in relation to Fig. 4). To demonstrate how seed-based analyses can be used to confirm layer-specific inter-area connections using resting-state fluctuations, we focus on the early visual system here. First, ROIs of the brain areas LGN, V1 and hMT+ were defined based on functional and structural localizers (see Fig.  4). Then, time series of the seed ROIs were extracted from all layers and orthogonalized to the time series of a corresponding control region (see Fig. 4). This means that the resulting time series solely contained temporal events that were unique to the seed of the ROI, without global temporal events (e.g. physiological noise). These time series were used as regressors in a conventional GLM analysis (FSL-FEAT ). The resulting activation and connectivity estimates were extracted as beta values and z-scores.

Feed-forward vs. feedback classification (in relation to
Figs. 5 and 6). According to the expectations of a canonical microcircuit (Felleman and Van Essen 1991), feed-forward driven input terminates mostly in middle cortical layers, while feedback input terminates in superficial and deeper layers. With the new whole-brain layer-fMRI VASO connectivity data available nowadays, this simplified layer-model allows researchers to develop binary classification algorithms that determine whether any given columnar set of voxels is better described as predominantly feed-forward driver vs. predominantly feedback driven.
Here, two templates of feedback vs. were predetermined feed-forward driven layer-profiles were predetermined (blue and red in Fig. 6). In the feed-forward driven case, the template profile has one peak close to the center of GM. In the feedback driven case, the template profile has two peaks. One peak is located in the superficial layers and a secondary peak is located in the middle and deeper layers. Any columnar unit can then be classified as being feed-forward or feedback driven depending on whether the correlation of the profile is higher for either one of the template profiles. Each column's feedback or feed-forward dominance can be represented by color maps that code for the relative correlation strength to either of the two template profiles. Here, each column's layer profile was estimated in two ways: a) For an ROI independent approach to determining the feed-forward vs. feedback nature of a column, we estimated the hubness of every layer within every columnar unit. Similar to the implementation in AFNIs 3dT-corrMap (Cox 1996), we calculated the correlation of the time series of each individual layer with all other layers. When the correlation of any given layer is high, this suggests that this layer represents the overall ongoing correlation in the entire columnar unit. When the correlation of any given layer is low, it indicates that this layer is not contributing a lot to the overall ongoing fluctuations in that columnar unit. The resulting feed-forward vs. feedback driven map of this ROI independent classification can be used to indicate which areas have resting-state fluctuation that are predominantly coming from feed-forward or feedback input.
For a graphical depiction of this algorithm, see Fig. 5. b) Alternatively, the feed-forward vs. feedback classification was also estimated in a seed-based approach. The time courses of manually selected seed regions were extracted and used as a regressor in a GLM analysis. The resulting beta-value maps were used to extract activity layer-profiles for every columnar structure. Finally, those beta-value layer-profiles were used to classify the feed-forward vs. feedback dominance. This procedure was repeated for 8 manually selected seed regions along the visual processing stream. The resulting feed-forward vs. feedback driven map of this seedbased classification algorithm can be used to determine which areas receive feed-forward input from the seed and which area receives feedback input from the seed. For a graphical depiction of this algorithm, see Fig. 6. Fig. 7). At conventional 1.5-3mm fMRI resolutions, with whole brain coverage, functional connectivity is often investigated by means of connectivity matrices, also described as functional connectomes. Here, we show how layer-dependent functional connectivity data can add additional dimensionalities and valuable directionality information to connectivitymatrix-analyses. Here, the Shen atlas (Shen et al. 2013) of 268 parcels was used to define approximate masks of brain areas. First, the parcels were transformed from MNI space to the individual participants EPI space with ANTs (Avants et al. 2008) using the non-linear warping SyN algorithm. Next, the time series of every layer was extracted individually for every brain area. Finally, Pearson correlation values of every time course with every other time course were estimated and depicted in connectivity matrix style. Fig. 8). Functional brain networks usually incorporate multiple local and distant brain areas at a macroscopic spatial scale. It has been shown, however, that functional networks can be further separated into smaller sub-networks (Braga and Buckner 2017;Heinzle et al. 2011;Smith et al. 2009). With VASO-based, submillimeter fMRI data, it becomes possible to investigate the topographical sub-division patterns of larger functional networks into smaller and smaller units, without unwanted signal leakage from macro veins. Here, we used iterative FSL Melodic (Multivariate Exploratory Linear Optimized Decomposition into Independent Components) ICA decomposition  to extract functional networks across multiple macroscopic and mesoscopic spatial scales. We focused on individual manually selected components and iteratively decompose them into smaller and smaller sub-components.

Iterative ICA (in relation to
a) First, we used a conventional FSL-Melodic ICA decomposition analysis pipeline with 30 subcomponents.
b) We individually selected components of interest (e.g. the sensorimotor network) and temporally regressed out all other components (including non-neural noise components). c) We repeated the ICA pipeline that subdivided the selected components into further sub-components. d) We repeated steps b-c multiple times.
2.3.6. Investigating individual differences with layer-dependent inter-subject correlation analyses (in relation to Fig.   9). To investigate participant-dependent layer-specific signal variations in movie watching tasks, we conducted the following procedure: a) fMRI time courses during movie watching tasks were extracted from all 98 (non-sibling) 7T-HCP participants in all ROIs (Shen parcels, (Shen et al. 2013)).
b) These data were temporally downsampled to match the layer-fMRI TR.
c) The downsampled time courses of all 98 7T-HCP participants were used as regressors in a FSL-FEAT ) GLM analysis of our highresolution layer-fMRI VASO data.
d) Beta-values were extracted for every brain area and plotted as layer-profiles. e) These participant-specific layer-profiles were then used to determine peak locations and draw conclusions about the feed-forward vs. feedback nature of selected areas within the network.
For a graphical depiction of this analysis, see Fig. 9. Figure 1 illustrates the acquisition challenges in layerdependent fMRI. For large coverage protocols, it can be beneficial to interleave the contrast generation periods with readout periods. In MAGEC VASO (Fig. 1A), the CBVweighting is introduced between every 3D-EPI readout block instead of a sole single inversion pulse. This allows an increase of the acquisition window and, thus enables more slices to be acquired. Similarly, for in-plane segmented readout strategies, multiple contrast-generating inversion pulses are applied per image. One of the most critical and most investigated challenges of layer-fMRI acquisition is the unwanted vascular bias of large draining veins. Fig. 1B conceptualizes the corresponding localization specificity in conventional GE-BOLD and blood volume sensitive VASO. It can be seen that VASO has a higher localization specificity, without strong signal leakage effects of large pial veins. However, it has a lower detection threshold with fewer significantly activated voxels. Two stripes of layer-dependent activity are visible in VASO, representing the input and output activity of layers II/III and layer Vb/VI, respectively. To make layer-fMRI connectivity mapping more attractive to neuroscience application research, it is desired to achieve a series of quality features including: high-sensitivity, whole-brain coverage, 0.5 mm resolutions, sub-second TRs, and minimal venous biases. While advanced layer-fMRI acquisition strategies can achieve each of those quality features individually, constraints can arise when all of these features are tried to be achieved simultaneously. Fig. 1C illustrates the layer-fMRI parameter space as a triangle of pushing coverage, resolution and TR. Fig. 2 highlights the advantage of CBV-based fMRI for layer-fMRI applications with an example of the primary visual cortex. It can be seen how the CBV signal change is dominated from the middle cortical layers at the location of the Stria of Genari (green arrow). BOLD signal, on the other hand, shows the largest signal peaks in the superficial layers, most probably due to unwanted signal spillage from deeper layers. The signal variations along the columnar dimension, however, does not seem to be that much affected by signal spillage. This might be coming from the fact that most of the larger intra-cortical veins are orthogonally-oriented to the cortical surface. These results suggest that the CBV-contrast is most vital to interpret single-contrast variations across cortical layers. For the interpretation of single-contrast columnar variations, a corresponding specificity advantage of CBV is less clear. Fig. 3 illustrates the capability of whole-brain layer-fMRI sequences to reliably capture connectivity measures across the brain. It can be seen how CBV-fMRI allows the extraction of conventional functional networks without unwanted biases of the pial vasculature.

Characteristics of layer-dependent connectivity analysis approaches
High-quality layer-fMRI data, as shown in Figs. 1-3, open the door to a new class of connectivity analyses that have not been possible without layer-fMRI specificity. These data are still emerging and are just starting to become available. Nevertheless, we review a number of these new connectivity analysis methods, which we believe might become relevant for a larger research field in the near future (see Fig. 4-9). Layerdependent functional connectivity data can be used for: a) Mapping and confirming layer-dependent connectivity patterns with seed-based correlation analysis (see Fig.  4 and Polimeni et al. 2010a;Wu et al. 2018)). Such analyses were previously only accessible in invasive animal research Jung et al. 2019).
b) Data-driven hierarchy mapping with iterative seedbased connectivity clustering along cortical processing streams (see Fig. 5 with cation).
c) Layer-dependent hubness mapping to characterize the functional embedding of brain areas in larger networks (see Fig. 6 with caption).

Figure 1. CBV Imaging Method with VASO
Panel A) schematically depicts the working principle of the CBV-fMRI methods reviewed here. The CBV sensitivity of VASO is based on volume distributions between different T1-compartments of intravascular and extravascular space in a T1-weighted sequence. For layer-fMRI applications, it is customary that VASO is acquired interleaved with BOLD data (see different contrast in brain insets). For SS-SI-VASO, the T1-contrast is generated with an initial 180°-inversion pulse. For MAGEC-VASO, the T1-contrast is facilitated with variable flip angles interleaved with readout modules. Panel B) depicts the expected vascular origin of VASO, compared to GE-BOLD. Most of the oxygenation changes that GE-BOLD is sensitive to are happening in large draining veins (blue) that unidirectionally smear the functional contrast across layers. CBV changes, however, are believed to be dominated by micro-vessels that are closely located at the laminar-aligned neurons and synapses. Because of the high localization specificity of VASO, its functional activity maps can reveal laminar stripe patterns or superficial and deeper layers (here in human M1 ). Panel C) depicts data acquisition trade-offs of imaging coverage (a.k.a. FOV), isotropic resolution, and sampling efficiency (a.k.a. TR). While it is possible to perform CBVweighted laminar fMRI with whole-brain coverage, at sub-second sampling rates, and with very high resolution of 0.5 mm, it is not yet possible to achieve all of this at the same time.

Figure 2. High-resolution methods without large vein biases
The purpose of this figure is to illustrate the layer-dependent localization specificity of CBV-weighted fMRI methods. Panel A) illustrates the location of CBV change during a flickering checkerboard task at nominal 0.5 mm resolution. The underlay is the mean VASO signal with its inherent T1-contrast, which facilitates straightforward layerification in EPI space. The zoomed section shows that CBV changes are confined to the middle layers of the calcarine sulcus, without unwanted sensitivity to pial veins. Panel B) shows the lack of sensitivity to local veins in VASO and GE-BOLD signal. The underlay is a 0.2 x 0.2 x 0.5 mm FLASH image (3 averages) to illustrate the location of the Stria of Genari and the location of large intracortical veins (arrows in shades of gray). Layer-profiles reveal that the VASO signal change is dominated by the Stria of Genari (green arrow), the location of expected input from thalamus for this task. GE-BOLD signal, on the other hand, is dominated by superficial layers. This is presumably due to unwanted sensitivity to large draining veins. These data suggest that CBV fMRI can be advantageous to reveal laminar-specific correlates of neural activity. In contrast to layer-profiles, the columnar profiles show that VASO and BOLD exhibit very similar activity distributions. At the location of large diving veins (gray arrows), BOLD signal seems to be slightly larger. This effect is small compared to the overall variance along the cortical ribbon. These data suggest that the vein bias in GE-BOLD signal is not as severe in columnar fMRI as in layer-fMRI. d) Layer-dependent mapping of the functional human connectome to obtain valuable information for directional graph theory analyses (see Fig. 7 with caption). e) Iterative decompartmentalization of cortical systems to investigate the topographical working principles of macroscopic functional networks (see Fig. 8 with caption).
f) Investigating the layer-dependent source of individual differences with time-locked naturalistic tasks (see Fig.  9 with caption). g) Generalized psychophysiological interaction analysis (gPPI) to investigate interactions between signals originating from different depths (see (Sharoh et al. 2019)).

Open questions in layer-fMRI connectivity analyses
Figures 2-9 strongly suggest that it is technically possible to measure layer-dependent functional connectivity straightforwardly and reliably across the entire brain in living humans.
To use such data for routine system neuroscience interpretations, additional aspects need to be considered that are discussed below.

Figure 3. Functional connectivity data with large coverage VASO methods
With the variable flip angle MAGEC-VASO approach, the number of slices can be increased as desired. This comes at the cost of TR. Here, 104 slices were acquired in 8.3s at 0.8 mm isotropic resolution, without z-GRAPPA acceleration. The panels show characteristic functional networks during a movie watching task (average of four repetitions of 14 min each). Activation was estimated in three steps: 1.) extracting signal traces from 98 participants of the 7T-HCP movie dataset, 2.) resampling them to the same TR used for layer-fMRI acquisitions, 3.) using the HCP signal traces as regressors (in the design matrix) of a conventional GLM-analysis. Panel A) depicts the 'parietal network', panel B) depicts the 'default mode network', panel C) depicts the 'fronto-parietal network', and panel D) depicts the 'visual network'. The activation maps reveal clear layer-specific structures. For example, the zoomed section of the intraparietal sulcus exhibits a double-stripe pattern (green arrows) following the cortical ribbon.

Shared sources of signal fluctuations
As in conventional functional connectivity-based analyses at lower resolutions, 'connectivity' is estimated based on the temporal similarity of functional time series. This method is prone to overestimating the connectivity strength if multiple areas are affected by the same source of erroneous signal fluctuations. One prominent cause of unwanted signal fluctuation can be respiration induced signal changes, which can introduce biases of false positive connectivity between neurally unconnected brain areas. In layer-fMRI analyses, this can have a higher effect in the superficial layers with a larger vascular density compared to deeper layers with reduced vascular density. Another source of unwanted signal fluctuations, when estimating the connectivity of two brain areas of interest, is a common input from a third unknown area. In this case, the common input would induce the same fluctuations in both ROIs and make them look more connected. Figs. 4-7 depict a number of potential approaches to minimize the effect of these unwanted sources of shared variance. One approach would be to restrict connectivity interpretations to differential analyses. Unwanted sources of fluctuations can be removed, by orthogonalizing the seed region's time course with an appropriate control region (Fig. 4B). Another approach of differential layer-dependent connectivity analysis is depicted in Figs. 5 and 6. Using a template matching approach refrains from interpretations of absolute connectivity strengths. Instead, it is based on the differential similarity, i.e. it quantifies which template correlates stronger with the profile. Yet another approach to minimize the effect of shared sources of unwanted variance would be to restrict interpretations of layer-dependent connectivity to asymmetric off-diagonal elements in the connectivity matrices (Fig. 7C). Common sources of physiological noise would be symmetric Panel A) depicts a toy model of expected directional connectivity in the early visual system. The primary visual cortex V1 receives feed-forward input from the thalamus mostly in the middle/deep layer IV and it receives feedback input from V5/hMT+ mostly in superficial and deeper layers. With 0.8 mm fMRI resolution, input to superficial layers (II/III) is expected to be separable from input to layers IV/V/VI. However, layers IV and V/VI might be too close together to be separable with 0.8 mm resolution. Panel B) schematically illustrates the procedure for layer-dependent time course analysis. First, the time series of the seed ROI is extracted from all layers and is orthogonalized to the time series of a control region. This means that the resulting time series solely contains temporal events that are unique to the seed of the ROI, without global temporal events (e.g. physiological noise). This time series is used as a regressor in a conventional GLM analysis. The resulting activation and connectivity score are extracted as beta values or z-scores. Panel C) illustrates how the ROIs were defined in this study. While the thalamus and hMT+ are defined based on functional localizers, V1 is defined based on the occurrence and borders of the Stria of Gennari (black arrows). Panel D) depicts the resulting connectivity profiles across layers. As expected from the model in panel A), feed-forward connectivity is dominantly terminating in middle and deeper layers (black arrow), while feedback input has additional connectivity in the superficial layers (brown arrow). The blue profile refers to a seed region in the contra-lateral V1 and can be interpreted as a measure of overall ongoing fluctuations arising from global physiological noise, thalamic input, V5/hMT+ input, and other input, for comparison. The layer-dependent functional connectivity analyses conducted here are inspired by previous work from Jonathan Polimeni (2010a).
to the diagonal axis in both areas and can thus be identified and removed.

Inherent connectivity across layers
The brain is highly interconnected (Schuez and Braitenberg 2002). Each neuron receives input from up to 10000 other neurons. Most brain areas are connected to most other brain areas and all layers are connected to all other layers (Constantinople and Bruno 2013;Felleman and Van Essen 1991;Harris et al. 2019). Thus, as far as the time scales of fMRI concerns, it can be challenging to treat different parcels of the brain (brain areas, layers, columns) as independent entities. In fact, even before the onset of the fMRI signal (500 ms -600 ms after stimulus onset), the first input into a specific layer of a given brain area can spread the signal across all cortical depth and brain areas. Before the fMRI signal reaches its peak amplitude (6s-10s after the stimulus onset) the neural signal might have travelled to multiple other brain areas back and forth multiple times already. While this fast interconnectivity theoretically challenges the interpretability of layer-dependent functional connectivity, it does not seem to be a limiting factor in practice. Several potential mecha-nisms have been proposed, why this might be the case. They are discussed in the following sections.

Mind the magnitude.
While electrophysiology studies show fast neural activity propagation in the time scales of 30 ms -40 ms across different layers within the same brain area (Godlove et al. 2014;Ninomiya et al. 2015), the magnitude of the neural activity decays rapidly. The initial neural activity in the feed-forward input layers IV of the primary visual cortex is by far the strongest activity. The magnitude of the further propagated activity in superficial and deeper layers is significantly lower as about 25%-40% (Ninomiya et al. 2015). Even though inter-layer connectivity can contribute to the overall signal, it is often negligible in magnitude to the original initial layer-specific input.

The neuro-vascular coupling favors the first input.
Previous time-resolved layer-fMRI studies in animals have shown that the onset of the fMRI signal to secondary connected layers has a surprisingly small amplitude and is surprisingly late. It has been shown with line scanning in the rat barrel cortex that the feed-forward induced fMRI signal starts rising in input layer IV ≈ 500ms after the stimulus onset.

Figure 5. Mapping the columnar-specific layer hubness across brain areas
Panel A) depicts how the cortex is parceled into columnar structures. The resting-state time course of every columnar unit is extracted as a mean value. Panel B) illustrates how the layer-specific fMRI fluctuations are used to determine a functional measure of hubness. The term 'hub' is used here to describe nodes (e.g. layers) with exceptionally higher functional connectivity compared to other nodes. These nodes are thought to play a major role in the coordination of information flow within brain networks (Bullmore and Sporns 2009;Bullmore and Bassett 2011;Sporns et al. 2007). Here, hubness is defined as the correlation between the layer-specific time course and the mean time course of all remaining layers. Calculating the hubness of every layer in a column allows the generation of hubness layer-profiles. These profiles can be characterized based on their respective peak location in granular vs. agranular layers. Panel C) depicts an example of clustering the brain into feed-forward driven areas with largest hubness measures in the granular layer versus feedback driven areas with largest hubness measures in infra-granular or supra-granular layers. A bilateral gradient between frontal and parietal areas can be seen. Panel D) illustrates similarities of this anterior-posterior pattern with measures of cortical thickness and myelination.

Figure 6. Hierarchy mapping procedure by means of a seed-based layer-dependent clustering analysis
Panel A) First, characteristic layer-dependent profiles are determined. At 0.8mm resolution, feedback activity in superficial layers (II/III) and deeper layers (IV/VI) can be separated as two separate peaks (red). Feed-forward activity in the deep layer IV can be seen as a single peak in middle/deeper cortical depth (blue). At 0.8mm resolution, the layer IV peak cannot be separated from the layer V/VI peak with Nyquist sampling. Thus, the feed-forward hump looks very similar to the deeper feedback hump. Here, the feed-forward and feedback profiles are used in a differential analysis. This means that even though any given layer profile usually contains a superposition of feed-forward and feedback peaks, ultimately it only matters to which of the two templates the profile is more similar to. Panel B) For a given seed region, the layer-profile is determined for all column in the field of view. Here 'columns' are considered as smooth 1mm patches of the cortex. Each column's layer-profile can then be clustered into one of the predefined classes based on the highest correlation similarity (relative correlation strength). Columns with layer-profiles that are dominated from superficial and deeper layers can be considered to mostly receive feedback input from the seed region. Columns with layer-profiles that are solely dominated from middle layers can be considered to mostly receive feed-forward input. In an iterative approach, the seed region can be picked based on the location of the clusters from the previous step. Panel C) Example clusters for a number of seed-regions (indicated with green arrows). It can be seen that clusters of feed-forward and feedback dominance are bilaterally organized along the geodesic distance. Note that the separation into two cluster groups results in an algorithm-enforced simplified view of the brain hierarchy. In fact, it is not expected that single columns are either 100% feed-forward driven or 100% feedback driven. Instead, it is expected that most of the columns exhibit a superposition of the two. The algorithm, however, enforces two binary clusters solely based on the maximum similarity to the templates. This is also visible in the color scale of the two clusters.

Figure 7. Whole-brain layer-dependent connectome mapping
This figure shows a possible analysis approach and representative example data to exemplify what kind of information layer-fMRI can contribute to interpret the brain's connectome. Panel A) depicts the raw VASO EPI data quality for whole-brain layer-dependent connectomics. Panel B) illustrates how functional connectome matrices are commonly generated: First, the brain is parcelated into a number of brain areas (colored masks overlaid on brain refer to the Shen (2013) atlas). Then, the average time courses within each brain area is correlated against all other brain area's time courses. The combinations of all correlation values are summarized in a functional connectivity matrix. Any value refers to one edge of the brain connectome and represents the functional connectivity strength between two brain areas. Panel C) shows that the resolution of layer-fMRI can add an additional dimension in connectome analyses. Since each brain area can be subdivided into multiple layers (colored masks overlaid on the brain), each node in the whole-brain connectivity matrix represents a layer-to-layer connectivity matrix in itself. One example node is highlighted (cyan). Here, rows and columns refer to layers. Superficial layers are depicted at the top and on the left, while the deeper layers are depicted on the bottom and on the right. Off-diagonal elements can be used to interpret directional connectivity. High connectivity values on the bottom left suggest that the connectivity is dominated from connections between middle/deeper layers of area 2 and superficial layers in area 1. Area 1 sends input into feed-forward layers of area 2, while area 2 send feedback input to area 1 in the superficial layers. Panel D-G) depict representative layer-dependent connectivities of common large networks. Panel D) depicts the 'visual network'. Selected correlation diagrams between V1 and V5/hMT+ confirm data from Fig. 4D. Namely, V1 receives top-down feedback in superficial layers from V5, while V5 receives bottom-up input in the middle/deeper layers (red circles). Panel E) depicts the 'sensory motor network'. As expected from previous layer-fMRI studies , the primary motor cortex receives input from the sensory areas solely in superficial layers (dark blue ellipses). Panel F) shows an example of the 'default mode network'. Cyan ellipses highlight that the PCC is the only middle-layer dominated ROI. The other ROIs seem to be more feedback driven. This can be taken as an indication that the PCC is the major hub of the 'default mode network', while the other areas are being passively driven perhaps by PCC activity. Panel G) depicts the 'fronto-parietal network'. Orange squares depict how the superficial and deeper layers have strong within-region connectivity and weak connectivity between each other. They almost look like two independent brain areas. This is consistent with electrophysiology data previously presented in monkeys (Maier et al. 2010).

Figure 8. Network decomposition and sub-decomposition with iterative ICA
This figure aims to illustrate a connectivity algorithm to investigate the submillimeter topology of common macroscopic networks. Panel A) depicts representative ICA components in axial slices covering the sensory motor system. Panel B) depicts the manually selected network for further decomposition. Panel C) depicts how ICA can further decompose the selected network from panel B). The sub-components do not separate into different Brodmann areas (E.g., BA1, BA2, BA3a/b, BA4, BA6), instead they rather separate into body part representations that span across multiple involved brain areas (see blue vs. red arrows). This is consistent with previously shown results (Kuehn et al. 2017). Panel D) depicts yet another iteration of sub-decomposition of the blue network from panel C). The sub-components are now at the spatial scale of the voxel size (0.8 mm).
The sub-networks do not appear to separate into different layers (E.g. input and output layers in M1 that are ≈ 1.2 mm apart). Instead, they separate into columnar like subnetworks of 0.8-1.2 mm distance that are consistent across layers.
fMRI activity in secondary connected layers II/III is delayed by another ≈ 300 ms -400 ms and has a much smaller magnitude (Yu et al. 2014). The transition time that it takes until the secondary connected layers show fMRI signal change is longer than what is expected from the electrophysiology (Godlove et al. 2014;Ninomiya et al. 2015). Similarly, recent layer-dependent work across multiple areas in the sensory system has revealed that the fMRI onset of secondarilyconnected areas can be delayed by up to several seconds (Jung et al. 2019). Based on this empirical evidence, it has been hypothesized that the fMRI response might be more susceptible to the initial input of a brain area. Secondary activity in inter-connected layers might need more time to accumulate sufficient neural activity until the fMRI signal reaches a significant increase.

Functional connectivity is not the same thing as struc-
tural connectivity. Only because the layers are structurally connected, doesn't guarantee that these connections functionally are engaged. It usually depends on the content of the neural activity fluctuation. In the context of conventional fMRI resolutions, the difference between structural and functional connectivity has been extensively discussed and widely established (Bullmore and Sporns 2009;Friston 2011). While functional connectivity depends on structural connectivity, structural connectivity is not a guarantee of functional connectivity. For example, at the spatial scale of large brain areas, it is well established that V1 and V5 are highly interconnected. However, the engagement of this connectivity depends on the motion energy of the stimulus. Thus, when the visual stimulus contains a lot of motion components, the connection is more engaged as opposed to static stimuli. Thus, despite the high structural connectivity between V1 and V5, Figure 9. Potential connectivity procedure to investigate the layer-dependent source of individual differences Panel A-B) illustrates the same movie task was used across all participants (98 HCP participants at 1.6 mm and 12 participants with layer-resolutions). This means that movie related brain activity changes are expected to be synchronized across participants. Some synchronized activity events can be attributed to specific semantic labels (Huth et al. 2016). Panel C) depicts how this experimental setup is used to extract layer-dependent information of individual-differences. Time courses of 98 non-sibling HCP participants are extracted across ROIs. Here, an example of the interparietal sulcus is depicted. These participant specific signal traces are then used as regressors in a GLM analysis of the submillimeter data. Panel D) depicts individual differences of these time courses for three clusters (k-means) of participants. There are time frames, when all participants have very similar synchronized fMRI fluctuations (pink arrow) and time frames when they are less similar (green arrow). Layer-profiles reveal that these inter-participant differences are solely caused by superficial feedback layers. Middle feed-forward layers are more consistent across participants. Since all participants are looking at the same movie, their retina activity is expected to be identical. In addition, the extraction of the low-level visual features in the early visual brain areas are expected to be similar. It is not surprising that the personal experience differences of the movie must be contributing to the brain activity further along the visual hierarchy as feedback input. Here, we focus on the interparietal sulcus because it is robustly detectable across participants and because it is positioned relatively high in the cortical hierarchy. Thus, it presumably does not only represent low-level visual features that are expected to be independent of participants. Instead, it is expected to be 'cognitive' enough (i.e. receiving high-level input or performing high level output) to also represent participant-specific components of personal movie experiences. As part of the Fronto-Parietal-Network, this area has been suggested to be most affected by individual differences during movie watching tasks (Finn et al. 2015;Vanderwal et al. 2017).
in resting-state analyses, these areas can be investigated as two separate entities. This reasoning can be extended from the scale of Brodman areas to the scale of layers too (Sotero et al. 2010). Despite the fact that individual layers in every column have many structural connections, layer-specific functional connections can still be extracted from isolated cortical depths (Fig. 3-7).

Sampling speed of fMRI fluctuations compared to timing of neural fluctuations
Meaningful neural activation modulations are happening across a wide range of temporal frequencies. While depthdependent electrophysiology studies often focus on the modulation of neural activity, connectivity and phase amplitude coupling changes in the range of 50 ms -300 ms (Godlove et al. 2014;Sotero et al. 2015), optical imaging studies examine meaningful resting-state connectivity across from the regime of 100 ms up to the 10 s regime (Ma et al. 2016). Due to the hemodynamic delay of the vascular response, conventional resting-state fMRI focuses on signal fluctuations in the time frame of 6-10 s 4 , implying that fMRI is usually only sensitive to a small frequency window of a wide spectrum of neural fluctuations. The acquisition approaches for whole brain layer-dependent connectivity analyses discussed in Figs. 1C, 3, and 7, are optimized for this temporal frequency window of ≥10s. Since resting-state fMRI fluctuations follow the pattern of scale free dynamics (He 2011), the focus on this frequency window is expected to be largely representative of functional connections at any temporal scales. Future work in combining the MAGEC-VASO approach with acceleration in both phase-encoding directions (Huber et al. 2020a) will become important to confirm this temporal invariance.

Lack of universal layer-dependent models
The seminal meta-study from Felleman and van Essen (Felleman and Van Essen 1991) summarizes the layer-dependent feed-forward connections in layer IV vs. feedback connections mostly in superficial and deeper layers. This simplified model is often considered to be canonical and evident across the entire neocortex (Douglas and Martin 2004), and thus, provides the basis of many layer-fMRI studies to date. The universal applicability of this simplified layer-dependent model, however, has been recently called into question (Constantinople and Bruno 2013; Harris et al. 2019). Recent neuroanatomical studies point to a more complicated layout of layer-dependent hierarchy-defining connections than previously assumed. Future insight of appropriate layer-dependent models needs to be taken into account (Markov and Kennedy 2013) when interpreting layer-dependent functional connectivity data. In any case, the discovery of more complex layerspecific pathway-models of inter-layer communication further underscores the importance of high-resolution laminar fMRI for elucidating principles of the cortical connectome.
For more comprehensive discussions about the lack of accurate neural models for layer-fMRI interpretation, please see the discussion section of the layer-fMRI review article by Stephan et al. (2019).

Unconfirmed scalability, optimizability, and efficiency of the discussed connectivity methods
The field of layer-fMRI connectivity is still emerging, and most methods-focused UHF research labs are still in the process of streamlining layer-fMRI connectivity tools for their neuroscience colleagues. The working principle of the connectivity analysis strategies that are conceptually illustrated in Figs. 2-9 are reviewed here as a set of hypothetical research tools that will possibly shed new light into the layerspecific organization of the cortex in future research. Different to conventional fMRI analysis tools, such as blockdesigned GLM and resting-state correlation analyses, these tools did not yet experience a decade-long evolution of optimizations and validations.

Biases of microvascular density
VASO fMRI is just one method of a large zoo of fMRI acquisition sequences that have been proposed for application of layer-fMRI, including: GE-BOLD, SE-BOLD, GRASE, ASL, etc. All of them (including VASO) are based on neurovascular coupling and can only capture fMRI signal changes in vascularized tissue. This means that, whenever a voxel has a higher vascularization and a higher microvascular density compared to other voxels, it will exhibit a higher fMRI responsiveness (a.k.a. vascular reactivity). As a result, layers with a higher vascular density are expected to be biased toward a larger fMRI signal change, for any task conditions. Thus, independent of the localization benefits of VASO compared to other methods, the layer-dependent microvascular density can introduce biases in layer-profiles. One extreme example would be layer I. Since layer I is almost free of microvascularization (Duvernoy 1981;Weber et al. 2008), none of the layer-fMRI sequences (including VASO) is expected to capture meaningful layer-specific functional signal from layer I. In other layers (II-VI), the VASO bias of layerdependent vascular density is expected to be negligible for the majority of cases. In contrast to the large diving veins that bias the GE-BOLD signal, the VASO-relevant microvascular density is rather homogeneous across cortical depth in most brain areas. As such, in the extra striate cortex, variations of the vascular density are usually smaller than 15% across layers II-VI (Kennel et al. 2017;Weber et al. 2008). There are exceptions, however. In V1 and S1, the microvascular density can be as much as 25% higher in layer IV compared to other layers (II, III, V, and VI), which can eventually introduce biases in VASO profiles. Due to these variations in layer-dependent microvascular density, it can be challenging to interpret single-condition layer-profiles without reference conditions. This challenge can be mitigated by introducing the control-analysis steps described in the section shared sources of signal fluctuations. Alternatively, this bias can be mitigated by explicitly measuring the layer-dependent vascular reactivity, e.g. by means of gas calibrations, by means of resting-state, or by means of the residuals of the task design 2020).

Limits of layer resolvability
Cytoarchitectonically defined cortical layers can be considerably thinner than the resolutions of any of the data reviewed here. Individual cortical layers can have thicknesses between 100 µm and 800 µm. Until now, it is not fully understood what the ultimate limit of the localization specificity in bloodvolume-based layer-fMRI will be. While animal studies have shown that layer-dependent CBV responses can be locally precise to an accuracy level of 200 µm in the rat olfactory bulb (Poplawsky et al. 2019), the best current spatial resolution of human layer-fMRI with VASO is ≈ 500 µm (Fig. 2). Thus, it is important to note that the layer-specific connectivity tools discussed here, do not represent activity of individual neuron clusters in isolated layers with Nyquist sampling across cortical depth. Instead, the findings on directional connectivity reviewed here are based on depth-specific signals variations that represent variable super-positions of neural activity from multiple cytoarchitectonically-defined cortical layers.

Accessibility and usability of layerdependent CBV-fMRI
In the past 5 years, the number of layer-fMRI applications with VASO has multiplied faster than any other layer-fMRI acquisition methodology (Fig. 10A). Since the advent of human layer-fMRI VASO in 2014-2015, layer-fMRI has become the sole driver of VASO-fMRI in humans. While the number of layer-fMRI VASO studies is still building up, there are 25 published peer-reviewed journal articles on layer-fMRI VASO until today (see Fig 10A). This is a significant portion of the entire field of layer-fMRI (Fig. 10B). While most of the VASO sequence development has been conducted in a handful of research labs only (Johns Hopkins, Max-Planck Leipzig, SFIM at NIH, MBIC in Maastricht, and DZNE in Bonn), the layer-fMRI user-base has significantly extended beyond these sites. Fig. 10C depicts that there are currently more than 30 labs around the globe that are using layer-fMRI VASO (effective January 2020). Most of these layer-fMRI VASO labs are in Europe, followed by Asia. Despite the high 7T scanner density in North America, the number of layer-fMRI VASO labs in USA is still emerging.

Alternative fMRI methods to measure layerdependent connectivity
The connectivity analysis methods reviewed here, happen to be exclusively illustrated with example data using VASO. We believe, however, that the discussed connectivity tools are generally applicable to any fMRI acquisition contrast that can provide whole brain layer-specific signals without the sensitivity of macro-veins. Until today, VASO just happens to the only acquisition sequence that can provide such signals non-invasively in humans. Other promising candidate methods that might also be optimized in the future research to be usable for layer-dependent connectivity analyses are briefly discussed below 5 : a) Model-based GE-BOLD signal deconvolution: Having reasonable assumptions of the biophysical properties of the layer-dependent vascular architecture, the unwanted venous drainage effect in GE-BOLD can be accurately estimated and cancelled out from ROI-specific GE-BOLD layer-profiles (Havlicek and Uludag 2019;Heinzle et al. 2016;Markuerkiaga et al. 2016). As soon as this approach will be advanced to also work on a voxel-by-voxel level, it might also be able to provide layer-specific data for the connectivity analyses too.
b) SE-BOLD: Spin-echo EPI has been shown to have an improved localization specificity compared to GE-EPI for applications in layer-fMRI (Goense and Logothetis 2006;Zhao et al. 2006), but are limited in coverage, energy deposition and acquisition speed (Koopmans and Yacoub 2019). Recent advancements of those methods using multisection excitation by simultaneous spin-echo interleaving (MESSI) with complexencoded generalized slice dithered enhanced resolution (cgSlider) simultaneous multislice echo-planar imaging (Han et al. 2019) can mitigate these limitations and bring SE-BOLD methods closer to high-resolution, whole brain protocols. And thus, they might also become usable for whole brain connectivity analyses.
c) 3D-GRASE: 3D-GRASE (Oshio and Feinberg 1991) has been shown to have an improved layer-dependent localization specificity compared to other BOLD contrasts (De Martino et al. 2013;Kemper et al. 2015). In most layer-dependent applications, however the imaging coverage in the segment direction did not exceed more than 12-18 slices. Recent advancements with using variable flip angles (Kemper et al. 2016) and compressed sensing ) strategies can mitigate the coverage limitations to some extent.

Conclusion
This article provides an overview of submillimeter fMRI methodology to map layer-dependent functional connectivity across brain areas. With recent sequence advancements in fMRI contrast generation and readout strategies, it is possible to overcome previous limits of resolution, coverage, and venous contaminations. The layer-fMRI tools reviewed here provide a starting point for mapping layer-specific connections across the entire brain, in the context of cognitive neuroscience. Layer-fMRI using VASO represents a paradigm shift that promises to transform the methods driven field of layer-fMRI to the application domain and will make layer-fMRI a reliable tool of neuroscience in two important aspects. First, the ability to image the entire brain with sufficient sensitivity on UHF scanners around the world makes layer-fMRI tools applicable for neuroscience cortical circuitry questions Figure 10. Popularity of human layer-fMRI VASO across recent years and around the globe Panel A) depicts the frequency of peer-reviewed journal publications using VASO. In the decade following its discovery in 2003, many researchers focused on low-resolution VASO studies to fully characterize the working principle of its functional contrast. Only with the advent of submillimeter imaging protocols in 2014/2015, VASO found it's 'killer application': layer-fMRI. 2019 was the first year, when layer-fMRI became the sole driver of VASO fMRI. An itemized list of all layer-fMRI VASO studies can be found at https://layerfmri.com/VASOworldwide. Panel B) depicts the overall trend of layer-fMRI for reference, including all fMRI contrasts (VASO, GE-BOLD, SE-BOLD, GRASE, ASL, etc.). Note the different scaling of the y-axis compared to panel A. Panel C) depicts the distribution of research labs that use layer-fMRI VASO. The vast majority of layer-fMRI VASO research is being conducted in Europe, followed by Asia.
For an itemized list of all layer-fMRI VASO users with references and example data, see https://layerfmri.com/VASOworldwide. without the need of extensive MR-physics expertise. Secondly, many influential network theories of brain function that posit distinct connection across individual cortical layers (e.g., predictive coding (Stephan et al. 2019), hierarchical processing streams) may now be directly tested in humans. This will allow a large field of researchers to investigate mechanisms behind any number of neuropsychological and psychiatric phenomena, such as selective attention, impulse control, learning, adaptation, hallucinations, cognitive control, and multi-sensory integration (Lawrence et al. 2019), to name a few. from the Institute for Basic Science (IBS-R015-D1). We thank Dimo Ivanov for his help with implementing the 3D-EPI ASL/VASO in previous versions of the sequence. We thank our collaborators who are also investigating the visual system at high resolutions for helpful discussions: Elisha Merriam, Insub Kim, Robert Trampel, Saskia Bollmann, Zvi Roth, Markus Barth, Marcello Venzi, and Kevin Murphy. We thank Sriranga Kashyap for his ANTs registration script that was used to display the background in Fig. 2B. Structural data in the background of Fig. 2B have been acquired at SFIM with a customized FLASH sequence that was optimized in MBIC with friendly advice from Masaki Fukunaga. We thank Jonathan Polimeni and Laura Lewis for helpful comments on how to select LNG/Pulvinar ROI used in data of Fig. 4. We thank Sunil Patil and Erik van den Bergh from SIEMENS Healthineers for support with sequence code sharing via C2Ps. We thank Dave Jangraw for the stimulation presentation code that was developed in previous studies and was also used here. We thank Insub Kim for optimizing our SPM motion correction scripts that were used here. We thank Eli Merriam for suggesting temporal interpolation in the BOLD correction scheme used here. We thank Daniel Glen and Rick Reynolds for sharing C++ code that is used for the nii I/O of the LAYNII-analysis of this study. We thank David Feinberg for discussions that lead to the final implementation of Fig. 5C. We thank Federico DeMartino for discussions and contributions to previous versions of Fig.  6. We thank FMRIF (especially Andy Derbyshire), NMRF (especially Joelle Sarlls and Lalith Talagala), and Scannexus (especially Chris Wiggins), where the parts of the data were acquired. Data were also provided in part by the Human

Data and Software availability
Anonymized MRI data that are presented in this article can be anonymously downloaded from OpenNeuro. Volume maps of data presented in Fig. 2 can be downloaded at http://doi.org/10.18112/openneuro.ds002274.v1.0.2 (go to actualfiles/V1_LAYERING). The raw and processed data of multiple participants of the study shown in Fig. 4 and 6 are available here: www.doi.org/10.18112/openneuro.ds001547 .v1.1.0. Data shown in Fig. 5 and 8 are publicly available from previous published studies ). Underlying tables and lists used for Fig. 10 are available on https://layerfmri.com/VASOworldwide. All custom written software (source code) and evaluation scripts are available on Github (https://github.com/layerfMRI/repository). The authors are happy to share the MAGEC-VASO 3D-EPI MR sequence upon request via a SIEMENS C2P agreement. A complete list of scan parameters used in this study is available on Github (https://github.com/layerfMRI/Sequence_Github). The source code of the layer-specific analysis software is available on Github https://github.com/layerfMRI/LAYNII with analysis pipelines explained on www.layerfmri.com. All stimulation presentations were implemented in Psychopy and can be downloaded here: https://github.com /layerfMRI/Phychopy_git. The collection of movie clips is available on https://layerfmri.page.link/7Tmovie.