Layer-dependent functional connectivity methods.

CBV Imaging Method with VASOPanel 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 (Huber et al. 2017)).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 CBV-weighted laminar fMRI with whole-brain coverage, at sub-second sampling rates, and with high resolution of 0.5 mm, it is not yet possible to achieve all of this at the same time.

Popularity of human layer-fMRI VASO across recent years and around the globePanel 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.

Comparison of VASO with various GE-EPI and SE-based acquisition method.Panels A-F) depict the MRI sequences that are compared: VASO, GE-BOLD EPI, SE-BOLD EPI, T1ρ-prep TFE, T2-prep TFE, diffusion-weighted T2-prep TFE.Panels G-L) depict the raw image with functional activity elicited by finger tapping overlaid (12 min experiment). VASO (panel G) and SE-EPI (panel I) show indications of a double- layer response (black arrows).Panel M) depicts the respective layer profiles of the compared imaging contrasts.Panel N) summarizes the respective sensitivity and localized specificity of all contrast mechanisms. Here, specificity and sensitivity are approximated by means of the profile slope and activity z-score. For depiction of alternative approximations, see (Huber et al. 2017). Panel N shows that the fMRI contrasts typically exhibit either high sensitivity or high specificity, but not both (dotted line). In this form of comparing the different methods, VASO does not fall on this line. VASO shows a compromise of moderate sensitivity and moderate specificity.Panel O) depicts the expected vascular origin of the respective methods.

High resolution methods without large vein biasesThe 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 while not in the GE-BOLD signal. The underlay is a 0.2 × 0.2 × 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.

Functional connectivity data with large coverage VASO methodsWith 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.

Resting-state directional connectivity in the visual systemThis figure depicts one proposed way of using resting-state time series analyses to investigate layer-dependent functional connectivity.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 Polimeni (2010a).

Mapping the columnar-specific layer hubness across brain areasPanel 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.

Hierarchy mapping procedure by means of a seed-based layer-dependent clustering analysisPanel 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.

Whole-brain layer-dependent connectome mappingThis 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 (Huber et al. 2017), 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).

Network decomposition and sub-decomposition with iterative ICAThis 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). The collection of four maps on the top of panel C) depict ICA maps as separate figures. The map at the bottom depicts two of those independent components in different colors superimposed on each other. The blue and red arrows show that a red-blue pattern is visible across sulci and gyri of BA1, BA 3b, BA 4, and BA6. This pattern looks consistent with finger representation maps of tapping induced activity (inset).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). Seven individual ICA components are depicted as individual maps at the right. The location of each component in the hand knob is highlighted with colorful arrows. These very components are also shown superimposed to each other in an enlarged panel on the left with the same arrows. 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.

Potential connectivity procedure to investigate the layer-dependent source of individual differencesPanel 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).