Experimental Procedure of the main experiment

1.2s prior to the stimulus onset, the fixation indicator would turn from green to red, indicating the subject to avoid blinking. After a period of 0.4s (0.8s before stimulus onset) the last fMRI volume of the previous trial finished recording. For a period of 3s no fMRI data was collected 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.6s 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 subject 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 collected. Additionally, a high resolution (0.8 mm iso voxel size) T1 weighted full brain image has been collected before the main experiment to obtain the individual subject’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).

Intermediate results for the EEG and fMRI data.

A) Full brain DICS beamformer results. Subject 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. Subject 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.

False discovery rate (FDR) adjusted p-values, computed using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995) for the respective voxel sub-selections in V1.

FDR corrections where computed on the p-values resulting from either the cluster permutation tests (Maris and Oostenveld, 2007) or the corresponding aros tests Clausner and Gentili, (2022) across the respective selection thresholds (5%, 10% and 25%). The circle symbol (°) indicates a p-value smaller than 0.05 before multiple comparison correction was applied. Conditions are abbreviated as follows: fu act.: feature unspecific activation; fs cont.: feature specific contrast; fs incr.: feature specific BOLD increase. See also Table 1 (supplementary material) for more information and Figure 4 and Figure 3 for a visual representation.

Average β coefficients

(average contribution per voxel) for feature unspecific activation. For each frequency bin, a separate regressor has been used and results were weighted with the respective layer contribution weights. Frequency bins (0.5 Hz steps for α and 2.5 Hz steps for γ) are indicated at the out-most half circle. Layers are indicated by coloured abbreviations for superficial (S), middle (M) and deep (D) layers. Percentages indicate the threshold at which voxels were selected and refer to the respective sub-selection (e.g. 5% would refer to the 5% most positive or negative t-values given a respective subselection). More saturated areas indicate a significant cluster after a cluster permutation test (Maris and Oostenveld, 2007) and FDR correction (Benjamini and Hochberg, 1995). For each cted significant cluster an aros test (Clausner and Gentili, 2022) across all layers for the widest possible frequency range that the cluster spans. In case of a significant result after FDR correction, white circles indicate in which layer the respective effect was strongest (3 circles), medium strong (2 circles) and weakest (1 circle). The strongest effect in the α band corresponds to the most negative β values, whereas in the γ band the strongest effect refers to the most positive values. Corresponding cluster and aros p-values are shown in Table 1.

Average β coefficients for V1

(average contribution per voxel) resulting from GLM using EEG power regressors to predict the BOLD signals for various voxel sub-selections. For each frequency bin a separate model has been computed and results were weighted with the respective layer contribution weights. Frequency bins (0.5 Hz steps for α and 2.5 Hz steps for γ) are indicated at the out-most half circle. Layers are indicated by coloured abbreviations for superficial (S), middle (M) and deep (D) layers. Voxel were selected based on the first level contrast as follows: (A) Voxel that respond stronger for one orientation over the other irrespective of the response to one condition alone; (B) Voxel that respond strictly positive for either orientation and stronger for one orientation over the other; See also Table 1 for an overview of how the conditions are constructed. Percentages indicate the threshold at which voxels were selected and refer to the respective sub-selection (e.g. 5% would refer to the 5% most positive or negative t-values given a respective sub-selection). Congruent thereby means that if a voxel e.g. responds stronger to left oriented gratings over right oriented gratings, it was paired with EEG based predictors for the left orientation trials (and vice versa). In turn, incongruent means that if a voxel e.g. responds stronger to left oriented gratings over right oriented gratings it was paired with EEG based predictors for the right orientation trials (and vice versa). More saturated areas indicate a significant cluster after a cluster permutation test Maris and Oostenveld, (2007) and FDR correction (Benjamini and Hochberg, 1995). Black dots indicate a cluster p-value of > .05 after, but < .05 before the correction was applied. For each significant cluster, an aros test (Clausner and Gentili, 2022) has been conducted across all layers for the widest possible frequency range that the cluster spans. In case of a significant result after FDR correction, white circles indicate in which layers the respective effect was strongest (3 circles), medium strong (2 circles) and weakest (1 circle). If the aros p-value was < .05 before, but > .05 after the correction, aros shapes are indicated using gray circles. The strongest effect in the α band corresponds to the most negative β values, whereas in the γ band the strongest effect refers to the most positive. Corresponding cluster and aros p-values are shown in Table 1. The right sub-panel indicates the average upper (11 − 13 Hz) and lower (8 − 10 Hz) α power. The p-value was calculated based on a linear mixed effect model controlling for individual layer activation and has been FDR corrected (Benjamini and Hochberg, 1995).

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 build 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, both stimulus orientations vs baseline, each stimulus orientations vs baseline or the contrast between both orientations 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.1s to 0.8s 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).