Experimental Procedure of the main experiment.

1.2 s prior to the stimulus onset, the fixation indicator would turn from green to red, indicating the participant to avoid blinking. After a period of 0.4 s (0.8 s before stimulus onset) the last fMRI volume of the previous trial finished recording. For a period of 3 s no fMRI data was recorded in order to avoid gradient artifacts in the EEG data. The stimulus presented is a left or right (± 45°from vertical axis) oriented grating. 16.7% of stimuli (2 × 8.3% for left and right respectively) were modified to have a slight wavy pattern (see example grating). Participants were asked to respond to those oddball trials with a button press. If no oddball trial was presented, the stimulus would remain on the screen for 1.6 s, followed by a 0.6 s period where only the fixation indicator is shown. In case of a response, corresponding feedback (“correct”, “false alarm” or “miss”) is displayed instead. After that, the fixation indicator turns back to green, indicating to the participant that the period to avoid blinking has ended. Now, three consecutive 3D EPI volumes (TR : 3.3 s) are recorded, before the next trial starts. Overall, 240 trials (four blocks with 60 trials each) have been recorded. Additionally, a high resolution (0.8 mm iso voxel size) T1 weighted full brain image has been recorded before the main experiment to obtain the individual participant’s anatomy. Furthermore, three blocks of pRF mapping (128 trials each) have been performed after the main experiment using the same fMRI sequence as in the main experiment (without gaps).

Data preparation and combined EEG-fMRI analysis for three cortical layers.

A) Layer specificity. Cortical layers are constructed as shell-like meshes, taking the gray and white matter boundaries as reference. Between those two reference shells two additional shells divide the space between gray and white matter into three layers. The area outwards relative to the gray matter boundary is assigned to the CSF layer, whereas the area inwards from the white matter boundary is assigned to the white matter layer. Each voxel contributes a fraction of its signal to those layers, depending of the proportional volume of a voxel within each shell like mesh of each layer (see example). Respective fractions are later used as weights to split the β coefficients resulting from GLM into different layer contributions. B) EEG based regressors. After transforming the EEG data into virtual channel data for each grid point within the gray matter (using LCMV beamforming), virtual channels are time-frequency transformed. For each hemisphere and frequency band separately (2 Hz to 32 Hz for α and 20 Hz to 120 Hz for γ) the virtual channel with the highest γ power increase or α power decrease after stimulus onset is selected respectively (∑ = 4 virtual channels). Regressors are built for each time-frequency bin separately. Time bins are averaged to boost SNR. Power values over trials are convolved with the HRF as built into SPM12 resulting in one parameter modulation regressor for each TF bin. C) Voxel selection. Voxels are selected based on t-maps resulting from first level contrasts. Thereby, the contrast between both stimulus orientations, and the response to any or both stimulus orientations alone are considered. D) Statistical inference. EEG based regressors are entered into a general linear model (GLM) as predictors for the BOLD signal in each selected voxel. The resulting β coefficients for each voxel are multiplied with the respective layer weights (excluding white matter and CSF layers) in each voxel and averaged for the respective time window of interest (0.1 s to 0.8 s after stimulus onset) to obtain the final depth by frequency resolved data. This data was tested against the hypothesis that there was no significant relationship between the EEG and fMRI data (β coefficients do not differ from zero) for each respective condition, using a cluster permutation test (Maris and Oostenveld, 2007). Resulting clusters were averaged in each layer over frequencies for the widest possible window selected across layers. The layer profiles of the averaged clusters were tested against the hypothesis that the layer profile is as likely as any other layer profile - under the assumption of interchangeability of the data - using an auto-regressive rank order permutation (aros) test (Clausner and Gentili, 2022).

Intermediate results for the EEG and fMRI data.

A) Full brain DICS beamformer results. Participant average of log-ratios between stimulus and baseline for 11 Hz (α) and 60 Hz (γ) for both sets of separately filtered EEG data. This serves illustrative purposes only, since the virtual channels of interest have been selected from time-frequency transformed virtual channels obtained using LCMV beamforming. Here, the 5% vertices with the strongest decrease (top) or increase (bottom) are shown. B) Time-frequency representation of virtual EEG channels. Participant average of log-ratios between stimulus and baseline of time-frequency transformed virtual channels, obtained using LCMV beamforming (2 × 2 virtual channels for low and high frequencies and both hemispheres separately). Only the right hemispheric channels are shown. The white empty square indicates the data points that were included in the combined EEG-fMRI analyses. Average reaction time and stimulus onset are indicated by a continuous or dashed white line, respectively. C) Average t-value distribution. Surface projection of average t-map of the first level contrast for the general fMRI activation (top) and the contrast between left and right stimulus orientation (bottom) for illustrative purposes.

V1 feature-unspecific relationship between laminar BOLD signal and EEG.

Average coefficients from a GLM predicting laminar BOLD signals in V1 from trial-by-trial EEG spectral power, shown separately for voxels with a strictly positive (Pos) or strictly negative (Neg) response to both stimulus orientations. 10% most extreme t-values were selected. coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. (A): α-band (8 – 14 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). (B): average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V1 relationship between laminar BOLD signal and EEG for the feature-specific contrast and feature specific BOLD increase.

Average coefficients from a GLM predicting laminar BOLD signals in V1 from trial-by-trial EEG spectral power. Voxels where selected based on the first level fMRI contrast of left – right oriented stimuli (feature contrast; A–C) or with an additional constraint, selecting only voxels with a positive response to any stimulus orientation alone (feature contrast BOLD increase; D–F). 10% most extreme t-values were selected. EEG regressors were computed from trials of the Congruent orientation (matching the preferred voxel orientation) or the Incongruent orientation (opposite). Furthermore, the difference between congruent and incongruent models has been computed (Congruent - Incongruent). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. A, D: α-band (8 – 14 Hz) layer × frequency plots. B, E: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. C, F: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995). See Figures S4 and S5 in Supplementary Figures for results for the 5%and 25% thresholds.

V1 feature-unspecific relationship between laminar BOLD signal and EEG.

Average coefficients from a GLM predicting laminar BOLD signals in V1 from trial-by-trial EEG spectral power, shown separately for voxels with a strictly positive (Pos) or strictly negative (Neg) response to both stimulus orientations. Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V2 feature-unspecific relationship between laminar BOLD signal and EEG.

Average coefficients from a GLM predicting laminar BOLD signals in V2 from trial-by-trial EEG spectral power, shown separately for voxels with a strictly positive (Pos) or strictly negative (Neg) response to both stimulus orientations. Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V3 feature-unspecific relationship between laminar BOLD signal and EEG.

Average coefficients from a GLM predicting laminar BOLD signals in V3 from trial-by-trial EEG spectral power, shown separately for voxels with a strictly positive (Pos) or strictly negative (Neg) response to both stimulus orientations. Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V1 relationship between laminar BOLD signal and EEG for the feature-specific contrast.

Average coefficients from a GLM predicting laminar BOLD signals in V1 from trial-by-trial EEG spectral power. Voxels where selected based on the first level fMRI contrast of left – right oriented stimuli (feature contrast). Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). The GLM was computed using voxels with a stronger response to one orientation over the other. EEG regressors were computed from trials of the Congruent orientation (matching the preferred voxel orientation), the Incongruent orientation (opposite). Furthermore, the difference between congruent and incongruent models has been computed (Congruent - Incongruent). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V1 relationship between laminar BOLD signal and EEG for the feature-specific contrast.

Average coefficients from a GLM predicting laminar BOLD signals in V1 from trial-by-trial EEG spectral power. Voxels where selected based on the first level fMRI contrast of left – right oriented stimuli and their (positive) response to either stimulus orientation (feature contrast BOLD increase). Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). The GLM was computed using voxels with a stronger response to one orientation over the other. EEG regressors were computed from trials of the Congruent orientation (matching the preferred voxel orientation), the Incongruent orientation (opposite). Furthermore, the difference between congruent and incongruent models has been computed (Congruent - Incongruent). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V2 relationship between laminar BOLD signal and EEG for the feature-specific contrast.

Average coefficients from a GLM predicting laminar BOLD signals in V2 from trial-by-trial EEG spectral power. Voxels where selected based on the first level fMRI contrast of left – right oriented stimuli (feature contrast). Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). The GLM was computed using voxels with a stronger response to one orientation over the other. EEG regressors were computed from trials of the Congruent orientation (matching the preferred voxel orientation), the Incongruent orientation (opposite). Furthermore, the difference between congruent and incongruent models has been computed (Congruent - Incongruent). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V2 relationship between laminar BOLD signal and EEG for the feature-specific contrast.

Average coefficients from a GLM predicting laminar BOLD signals in V2 from trial-by-trial EEG spectral power. Voxels where selected based on the first level fMRI contrast of left – right oriented stimuli and their (positive) response to either stimulus orientation (feature contrast BOLD increase). Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). The GLM was computed using voxels with a stronger response to one orientation over the other. EEG regressors were computed from trials of the Congruent orientation (matching the preferred voxel orientation), the Incongruent orientation (opposite). Furthermore, the difference between congruent and incongruent models has been computed (Congruent - Incongruent). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V3 relationship between laminar BOLD signal and EEG for the feature-specific contrast.

Average coefficients from a GLM predicting laminar BOLD signals in V3 from trial-by-trial EEG spectral power. Voxels where selected based on the first level fMRI contrast of left – right oriented stimuli (feature contrast). Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). The GLM was computed using voxels with a stronger response to one orientation over the other. EEG regressors were computed from trials of the Congruent orientation (matching the preferred voxel orientation), the Incongruent orientation (opposite). Furthermore, the difference between congruent and incongruent models has been computed (Congruent - Incongruent). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

V3 relationship between laminar BOLD signal and EEG for the feature-specific contrast.

Average coefficients from a GLM predicting laminar BOLD signals in V3 from trial-by-trial EEG spectral power. Voxels where selected based on the first level fMRI contrast of left right oriented stimuli and their (positive) response to either stimulus orientation (feature contrast BOLD increase). Columns correspond to voxel-selection thresholds (top 5%, 10%, 25% most extreme t-values). The GLM was computed using voxels with a stronger response to one orientation over the other. EEG regressors were computed from trials of the Congruent orientation (matching the preferred voxel orientation), the Incongruent orientation (opposite). Furthermore, the difference between congruent and incongruent models has been computed (Congruent - Incongruent). coefficients were weighted by laminar contribution weights, yielding estimates for superficial (S), middle (M), and deep (D) layers. a–c: α-band (8 – 14 Hz) layer × frequency plots. d–f: average (± SEM) for lower (8 – 10 Hz) and upper (11 – 13 Hz) α sub-bands per voxel selection, tested with a linear mixed-effects model. g–i: γ-band (50 – 70 Hz) layer × frequency plots. Black rectangles indicate significant clusters (Maris and Oostenveld, 2007); where significant, lollipops to the right show the layer rank order of the effect (Clausner and Gentili, 2022). All p-values are FDR-corrected (Benjamini and Hochberg, 1995).

Signal-to-Noise Ratio.

Several measures have been computed to estimate the temporal signal-to-noise ratio across layers, as indicated by the title of the respective sub-plot. All sub-plots except the rightmost sub-plot have been computed based on the used layering algorithm. We have included the tSNR estimate for the laminar segmentation as provided by LayNii (Huber et al., 2021), to allow for a comparison between our fraction-based approach and LayNii’s approach, which assigns entire voxels to specified layers. In general, the laminar profile of tSNR estimates was unexpected, because the vascular draining effect (Markuerkiaga et al., 2016) typically leads to higher signal changes in superficial layers. However, since both functionally relevant signal components, as well as physiological noise drains towards the cortical surface, a higher tSNR for deep layers is not entirely surprising. Importantly however, the tSNR profile is not reflected in our main results, which strengthens the overall validity.

EPI example with layer segmentation.

Left: a segment of a slice of a single participant’s raw fMRI volume during the task. Centre: The same segment, overlaid with pial and white matter boundaries that were used for reconstructing cortical depth. The white box indicates an example voxel selection for which the respective layer attributions are shown on the Right: Each voxel was assigned a fraction by how much this voxel belongs to a specific layer.

Correlation of α Frequency and Task Performance.

Task performance for each participant has been defined as the average reaction time (RT) to correctly identified oddball stimuli, as well as d’ as a measure for average response accuracy. Individual α frequencies have been determined for each participant based on the average over all correctly identified non-oddball trials. Within the time of interest (0.1 to 0.8 s after stimulus onset) and frequencies of interest (8 to 14 Hz) the frequency with the largest decrease was chosen. Afterwards, per-participant frequencies have been correlated with per-participant average RTs and d’ values. We found that participants with the strongest decrease of α power in higher frequencies, also respond more accurately on average. No such correlation was observed between α frequency and RT. However, RT and d’ have been found to be negatively correlated as well.

Average Motion of Participants.

Average across participants for average and maximum frame-wise displacement (Power et al., 2012) using a radius of 50 mm, total translation and total rotation within one experimental block.