Decoding of columnar-level 1 organization across cortical depth 2 using BOLD-and CBV-fMRI at 7 T 3

16 Multivariate pattern analysis (MVPA) methods are a versatile tool to retrieve information from 17 neurophysiological data obtained with functional magnetic resonance imaging (fMRI) techniques. 18 Since fMRI is based on measuring the hemodynamic response following neural activation, the 19 spatial speciﬁcity of the fMRI signal is inherently limited by contributions of macrovascular 20 compartments that drain the signal from the actual location of neural activation, making it 21 challenging to image cortical structures at the spatial scale of cortical columns and layers. By 22 relying on information from multiple voxels, MVPA has shown promising results in retrieving 23 information encoded in ﬁne-grained spatial patterns. We examined the spatial speciﬁcity of the 24 signal exploited by MVPA. Over multiple sessions, we measured ocular dominance columns 25 (ODCs) in human primary visual cortex (V1) with different acquisition techniques at 7 T. For 26 measurements with blood oxygenation level dependent contrast (BOLD), we included both 27 gradient echo-(GE-BOLD) and spin echo-based (SE-BOLD) sequences. Furthermore, we acquired 28 data using the vascular-space-occupancy (VASO) fMRI technique, which is sensitive to cerebral 29 blood volume (CBV) changes. We used the data to decode the eye-of-origin across cortical depth. 30 While ocularity information can be decoded with all imaging techniques, the macrovascular 31 contributions in GE-and SE-BOLD limit their use for discriminating signals between cortical layers. 32 However, the cortical proﬁle of decoded ocularity information from VASO measurements better 33 reﬂects the expected proﬁle of neural activity, suggesting the combination of VASO and MVPA to 34 be a promising approach for investigating the mesoscopic circuitry of the human cerebral cortex. 35


Introduction
In the cerebral cortex, neurons tend to cluster into functional units across cortical depths (Mountcastle, 1957;Hubel and Wiesel, 1962), which are usually called cortical columns and often denoted as the fundamental building blocks of the cortex (Mountcastle, 1997).A prominent example is Haenelt et al. | bioR𝜒iv | September 28, 2023 | 1-33 found in the primary visual cortex (V1).V1 mainly receives thalamocortical projections from the lateral geniculate nucleus (LGN) (Wandell, 1995), which contains monocular neurons that are segregated into eye-specific layers (Andrews, Halpern, and Purves, 1997).The monocular information is preserved when entering V1, and projections from the left and right eye are sent to segregated cortical columns, widely known as ocular dominance columns (ODCs) (Hubel and Wiesel, 1969;Tootell, Hamilton, et al., 1988;Dougherty et al., 2019), which form a repeating stripes pattern of alternating eye preference (Adams, Sincich, and Horton, 2007).
Functional magnetic resonance imaging (fMRI) is a versatile neuroimaging technique for noninvasive measuring and mapping of brain activity by assessing the hemodynamic response following neural activation (Buxton, 2013).However, due to the limited spatial resolution, conventional fMRI techniques only allow the detection of relatively large pieces of cortex involved in the execution of a specific task (Glover, 2011).Therefore, ODCs with an approximate column width of around 1 mm in humans (Adams, Sincich, and Horton, 2007) and other cortical columns were out of reach for usual fMRI applications.
With the development of MR scanners with higher magnetic field strengths and more sophisticated radiofrequency (RF) coils providing higher signal-to-noise ratio (SNR), mesoscopic structures like ODCs became accessible in humans at the expense of prolonged acquisition times and usage of anisotropic voxels (Menon, Ogawa, et al., 1997;Menon and Goodyear, 1999;Dechent and Frahm, 2000;Goodyear and Menon, 2001;Cheng, Waggoner, and Tanaka, 2001;Yacoub et al., 2007).Only with the emergence of ultra-high field MRI at a field strength of 7 Tesla and above, it became possible to measure ODCs with isotropic voxels at sub-millimeter resolution (Nasr, Polimeni, and Tootell, 2016;Hollander et al., 2021;Zaretskaya, Bause, et al., 2020;Akbari et al., 2023).
Given the average cortical thickness of 2-4 mm (Fischl and Dale, 2000) and its convoluted structure, the use of isotropic voxels at sub-millimeter resolution is necessary for the reliable sampling of data at different cortical depths (Turner and Geyer, 2014).This recent possibility is intriguing since the cerebral cortex is known to be composed of several layers, e.g., in terms of cytoarchitecture (Brodmann, 1909), myeloarchitecture (Vogt and Vogt, 1919), and vasculature (Duvernoy, Delon, and Vannson, 1981).Furthermore, cortical layers generally differ in their connectivity profile within and to other cortical areas, e.g., feedforward and feedback signaling between cortical areas in a hierarchically organized cortical system (Felleman and Van Essen, 1991).Thus, the mapping of cortical columns at different cortical depths with fMRI enables studying the local microcircuitry of the cerebral cortex in vivo (Yang et al., 2021).
The monocular feedforward signal from the LGN enters V1 in layer 4C of corresponding ODCs (Kennedy et al., 1976;Dougherty et al., 2019).Layer 4C is located directly below layer 4B, which contains the highly myelinated external band of Baillanger, also called stria of Gennari (Trampel, Ott, and Turner, 2011).Typically, layer 4C is further divided into layers 4Cα and 4Cβ, which receive color-selective parvocellular and "color-blind" magnocellular input from corresponding LGN layers, respectively (Nieuwenhuys, Voogd, and Huijzen, 2008).Above and below layer 4C, the signals from the two eyes converge onto single neurons, which lead to a variable degree of ocularity across cortical depth.However, individual neurons of the same column still tend to receive input predominantly from either the left or right eye, respectively (Wandell, 1995).In this regard, V1 is the first main stage of binocular integration, which is important, for example, for the processing of stereopsis (Poggio, 1995).
However, fMRI applications most often rely on the blood oxygenation level dependent (BOLD) signal sampled with a gradient echo-based sequence (GE-BOLD) that is known to be most sensitive to macrovascular compartments of the cerebral cortex (Turner, 2002), specifically draining veins that carry the deoxygenated blood back to the cortical surface (Polimeni, Fischl, et al., 2010;Markuerkiaga, Barth, and Norris, 2016).This usually leads to a signal accumulation toward the pial surface, limiting the ability to associate the BOLD response with a specific cortical layer.Alternatively, a spin echo-based sequence (SE-BOLD) can be used at high magnetic field strengths (Yacoub et al., 2007).SE-BOLD promises a more specific signal due to the refocusing of extravascular signal contributions from around larger veins (Boxerman et al., 1995).This has the advantage of increasing signal weighting to the microvasculature, which is believed to be closer to the actual location of neural activation.Furthermore, recent advances of imaging approaches with contrast weighted by cerebral blood volume (CBV) using vascular-space-occupancy (VASO) fMRI at higher magnetic fields show promising results in terms of increased laminar specificity (Huber, Handwerker, et al., 2017;Huber, Finn, et al., 2021) at the expense of overall sensitivity.
Next to the choice of the proper acquisition technique, multivariate pattern analysis (MVPA) (Haxby, 2012) methods have been shown to retrieve information manifested in spatial patterns of fMRI activity, which promise increased sensitivity compared to univariate methods (Kriegeskorte and Bandettini, 2007;Formisano and Kriegeskorte, 2012;Vizioli et al., 2020), for example, for the dissociation of bottom-up and top-down processing into different cortical layers (Muckli et al., 2015;Kok et al., 2016;Iamshchinina et al., 2021).However, though the presence of pattern information provides strong evidence for neuronal effects, the spatial scale of the exploited information remains unknown (Formisano and Kriegeskorte, 2012).Interestingly, already at a conventional resolution of 3 × 3 × 3 mm 3 using GE-BOLD at 3 T, decoding of orientation information is possible from responses in V1 (Haynes and Rees, 2005a;Kamitani and Tong, 2005), which is known to be encoded at a much finer spatial scale at the level of cortical columns (Obermayer and Blasdel, 1993).In the same year, the eye-of-origin could also be decoded from V1 voxels based on a binocular rivalry stimulus (Haynes and Rees, 2005b).These studies started a controversy several years ago (Boynton, 2005;Op de Beeck, 2010;Swisher et al., 2010;Gardner, 2010;Shmuel et al., 2010;Kriegeskorte, Cusack, and Bandettini, 2010;Chaimow et al., 2011;Misaki, Luh, and Bandettini, 2013) about the source of the exploited information.Possible mechanisms were suggested like the aliasing of high spatial frequency information encoded above the Nyquist frequency of the MRI sampling process (Boynton, 2005) (but see (Chaimow et al., 2011)), the contributions from random irregularities of the fine-scale columnar pattern, which leads to information at low spatial frequencies (Haynes and Rees, 2005a;Kamitani and Tong, 2005;Kriegeskorte and Bandettini, 2007) or the exploitation of large-scale information that are not related to the fine-scale columnar pattern (Op de Beeck, 2010).Growing evidence showed that functional biases can also be introduced by large vessels (Turner, 2002;Gardner, 2010;Shmuel et al., 2010;Sengupta et al., 2017), which might be thought of as a sort of lowpass filtering of the neural pattern into a coarser venous spatial pattern (Formisano and Kriegeskorte, 2012).Therefore, neural patterns can generally be represented at multiple spatial scales in the fMRI signal however (Swisher et al., 2010;Sengupta et al., 2017).
To study the microcircuitry of the cerebral cortex, it is of importance to know the source of the decoded information, e.g., by relating the decoded information to specific cortical layers.In this regard, it might be appealing to use fMRI acquisition techniques that are less sensitive to large vessels in combination with MVPA methods to benefit from the increased sensitivity of multivariate methods, while keeping a high spatial specificity of the exploited signal.However, most decoding studies use the GE-BOLD technique, which is known to be inherently limited by macrovascular contributions, reducing the potential benefits.
In our study, we acquired ODC data from five participants using GE-BOLD, SE-BOLD, and VASO in different sessions to study the laminar specificity of the respective acquisition technique in combination with MVPA to decode the signal of the stimulated eye in V1.Functional data were acquired with nominal isotropic voxel size of 0.8 mm allowing data sampling at different cortical depths.
From the perspective of neural processing, we expected highest eye-of-origin discriminability in deeper cortical layers since eye-specific segregation is most preserved in the input layer 4C.However, due to the drainage of deoxygenated blood toward the pial surface, macrovascular contributions to the fMRI signal were expected to change or "smear" the discriminability across cortical depth.Therefore, studying decoding performances of a feedforward signal between acquisition techniques across cortical depth enables the analysis of their different sensitivities to draining vein contributions.We believe that this study gives insights into the capabilities and limitations of using

Participants
A total of five healthy volunteers participated in this study, of which two were female (age = 28.00 ± 2.61, mean ± standard deviation in years).Written informed consent was obtained from all participants, and the study received ethical approval from the local ethics committee of the University of Leipzig.All participants had normal or corrected-to-normal visual acuity.We performed the Miles Test (Miles, 1929) with each participant to determine eye dominance, which is stated in Supplementary Figure 1-Supplementary Figure 5 for single participants.

General procedure
Each participant was invited to multiple scanning sessions on different days with an ultra-high field MR scanner operating at 7 T.The first session was used for reference measurements, during which a high-resolution anatomical reference scan and retinotopy data (Sereno et al., 1995;Engel, Glover, and Wandell, 1997) were acquired.In addition, a high-resolution functional time series without task (GE-BOLD) was obtained using the same parameters as the functional measurements in the subsequent sessions, in order to aid with between-session registration.The remaining six sessions were exclusively devoted to ODC mapping (2× GE-BOLD, 2× SE-BOLD, 2× VASO).A subset of the retinotopy data had previously been utilized in another experiment (Movahedian Attar et al., 2020), but it underwent independent processing for this study.

Visual stimulation
For the purpose of visual stimulation, an LCD projector (Sanyo PLC-XT20L) with custom-built focusing objective lens was used (refresh rate: 60 Hz, pixel resolution: 1024 × 768) that was positioned in the magnet room.To prevent interferences with the MR scanner, the projector was housed within a custom-built Faraday cage.The stimuli were projected onto a rear projection screen, mounted above the participants' chest within the bore.Participants viewed the stimuli by means of a mirror attached to the head coil.In order to minimize scattered light reaching the participants' eyes, the projection screen was surrounded by black felt, and all ambient lighting was turned off during data acquisition.This setup allowed visual stimulation within an approximate visual angle of 22 • × 13 • (width × height).Stimulus generation and presentation were carried out using the Psychophysics Toolbox (3.0.14, http://psychtoolbox.org/) (Brainard, 1997;Pelli, 1997;Kleiner et al., 2007) with GNU Octave (4.0.0, http://www.gnu.org/software/octave/).

ODC mapping
We used a block design with two experimental conditions that was previously reported in detail (Nasr, Polimeni, and Tootell, 2016;Haenelt et al., 2023), with the following minimal modifications for the current study.Every scanning session comprised ten runs, each lasting for 270 s.Within each run, a baseline period of 15 s was placed at the beginning and end, during which participants were presented with a uniform black background.The experimental protocol consisted of eight blocks, each lasting for 30 s, allowing four distinct stimulation periods targeting the left and right eye, respectively.The ordering of blocks was pseudorandomized.Throughout the runs, participants were instructed to maintain fixation on a central point (0.2 • ×0.2 • ) and respond on a keypad when the fixation point changed its form (square or circle).Presented stimuli consisted of red or green random dot stereograms (RDS) (Julesz, 1971) shown on a black background (dot size: 0.1 • , dot density: ∼ 17%) that were viewed through custom-built anaglyph spectacles using Kodak Wratten filters No. 25 (red) and 44A (cyan) and enabled the stimulation of either the left or right eye in separate blocks, see Figure 1.RDSs performed a horizontal sinusoidal movement (temporal frequency: 0.25 Hz, amplitude: 0.11 • ), and phases of dots were initialized to create the appearance of an 8 × 6 checkerboard with independent movement of squares.To reduce cross-talk between the eyes, the luminance of the dots was maintained at a low level (red through red filter: 3.1 cd/m 2 , red through cyan filter: 0.07 cd/m 2 , green through cyan filter: 5.7 cd/m 2 , green through red filter: 0.09 cd/m 2 ).It is worth noting that the luminance of the green dots was approximately doubled relative to red to ensure a similar excitation of cone photoreceptors for both colors (Dobkins, Thiele, and Albright, 2000).

Imaging
We used a whole-body MR scanner operating at 7 T (MAGNETOM 7 T, Siemens Healthineers, Erlangen, Germany) for measurements.The scanner was equipped with SC72 body gradients (maximum gradient strength: 70 mT/m; maximum slew rate: 200 mT/m/s).We used a single-channel transmit/32-channel receive head coil (Nova Medical, Wilmington, DE, USA) for RF signal transmission and reception.To optimize the transmit voltage over the occipital lobe, we always acquired a low-resolution transmit field map at the beginning of each scanning session using a sequence that exploits the ratio of consecutive spin and stimulated echoes, which was provided by the scanner vendor.
For ODC mapping measurements, we acquired functional data with GE-BOLD, SE-BOLD, and VASO in different sessions.GE-and SE-BOLD data were acquired using a single-shot sequence with 2D echo planar imaging (EPI) readout (Feinberg, Moeller, et al., 2010;Moeller et al., 2010).For VASO measurements, we used a single-shot slice-selective slab-inversion (SS-SI) VASO sequence (Huber, Ivanov, et al., 2014) with a 3D EPI readout (Poser et al., 2010).An oblique-coronal slab was imaged positioned over the occipital pole.For all acquisition techniques, we used the following parameters: In C, the surfaces before (depicted in red) and after (depicted in green) alignment with the GBB method are presented.This technique is implemented in the GBB package (0.1.6,https://pypi.org/project/gbb/).Before alignment, surfaces were transformed into functional space via a rigid registration.The method improves spatial correspondence between the surfaces and the GM/WM boundaries observed in the functional images.GM: gray matter, WM: white matter.
parameters were the following: TE = 25 ms, TI = 1370 ms for the blood-nulling point, FA = 26 • , 7.7% slice oversampling.50 slices were acquired in GE-BOLD measurements that covered the whole stimulated area of V1.Due to specific absorption rate (SAR) constraints, the number of slices was limited for SE-BOLD and VASO measurements.For VASO, we acquired 26 slices.For SE-BOLD, we used the maximum number of allowed slices that varied between subjects and sessions and was between 16 and 29 slices.

Data preprocessing
Functional time series from individual ODC mapping sessions were first subjected to motion correction to address within-run and between-run motion using SPM12 (v6906, https://www.fil.ion. ucl.ac.uk/spm/) with Matlab R2019b (MathWorks, Natick, MA, USA).In the case of VASO measurements, the time series were initially separated into individual time series for nulled and not-nulled images prior to motion correction.Motion correction was then independently applied to each of these time series.Final VASO time series were obtained by correcting the nulled time series for residual BOLD contamination.To achieve this, the motion-corrected nulled and not-nulled VASO time series were temporally upsampled onto a common grid using 3drefit from AFNI (19.1.05,https://afni.nimh.nih.gov/)(Cox, 1996), matching the effective temporal resolution of GE-and SE-BOLD measurements.Subsequently, the nulled time points were divided by the not-nulled time points to perform BOLD correction (Huber, Ivanov, et al., 2014).All time series underwent then highpass filtering (cutoff frequency: 1∕270 Hz), and a voxel-wise statistical analysis was performed for each session using a fixed-effects general linear model (GLM) as implemented in SPM12 with both experimental conditions as regressors convolved with the canonical hemodynamic response function (HRF).Note that GLM results were only used to visualize statistical maps and for the repeatability analysis (see Consistency of ocular dominance maps), while the main analysis was based on raw fMRI time series.
The functional time series obtained from retinotopy measurements underwent similar preprocessing steps.However, prior to motion correction, each time series was corrected for different slice timings by voxel-wise temporal interpolation to a common time grid using 3drefit.Following motion correction, the time series were subjected to highpass filtering (cutoff frequency polar angle runs: 1∕192 Hz, cutoff frequency eccentricity runs: 1∕96 Hz), and the data from the first quarter stimulus cycle was discarded from further analysis.A voxel-wise Fourier transform was computed, and the signal at stimulus frequency was averaged from runs with opposite stimulus directions to compensate for the hemodynamic lag.The phase at stimulus frequency from polar angle runs was used to delineate the borders of V1.
The MP2RAGE (UNI) image was used for surface reconstruction of the cerebral cortex.Initially, the UNI image underwent bias field correction using SPM12.The corrected image was then fed into the recon-all pipeline in FreeSurfer (6.0.0, http://surfer.nmr.mgh.harvard.edu/)(Dale, Fischl, and Sereno, 1999;Fischl, Sereno, and Dale, 1999) with the hires flag to perform segmentation at the original voxel resolution (Zaretskaya, Fischl, et al., 2018).The brain mask was separately created based on the second inversion image of the MP2RAGE by using the SPM12 segmentation algorithm and excluding voxels in a binary mask that exceeded the tissue class threshold of 10% in all nonwhite matter (WM) and non-gray matter (GM) tissue classes.Subsequently, generated boundary enhance the robustness of this method, we increased the GM/WM contrast in functional images following the method suggested in Fracasso, Petridou, and Dumoulin, 2016 that weights the magnitude image by its phase (both provided by the online reconstruction on the sanner) as conventionally practiced in susceptibility-weighted imaging methods.For this purpose, the magnitude time series was corrected for motion using AFNI.Each image of the phase time series was individually phase unwrapped using the method by Abdul-Rahman et al., 2005 implemented in Nighres (1.2.0, https://pypi.org/project/nighres/)(Huntenburg, Steele, and Bazin, 2018), and computed motion parameters were subsequently applied to the unwrapped phase time series.The temporal mean of both magnitude and phase data was calculated, and the phase data underwent thresholding and normalization.Finally, the contrast of the magnitude data was enhanced by assigning weights to each voxel based on the contrast-reversed phase data.
Nine equidistant surfaces were computed and positioned between boundary surfaces.This resulted in 11 cortical layers for subsequent analysis.
To achieve registration between the whole-brain anatomy and the functional time series without task, the anatomical image underwent an initial transformation to align with the functional space based on the scanner coordinate system.Only for registration, the mean functional image was bias field corrected (Tustison et al., 2010).Both images were then brain-masked and rigidly registered using ANTs (2.3.1, http://stnava.github.io/ANTs/).A similar procedure was employed for registering functional images from other sessions to the functional time series without task, except that a nonlinear registration was performed using the Symmetric Normalization (SyN) algorithm (Avants et al., 2008) implemented in ANTs.
For sampling data onto reconstructed surfaces, surfaces were first moved into the space of individual functional sessions based on the computed registration.Subsequently, the functional data were sampled onto the surface mesh using linear interpolation.

Pattern classification
We used a linear support vector machine (SVM) algorithm for pattern classification from single time points of motion-corrected and detrended raw functional time series.Each ODC mapping session and each cortical depth was analyzed independently.For classification, functional time series were first sampled onto a cortical layer.One run contained 90 time points, and 10 runs were acquired per session.All time points from the baseline condition were discarded.Additionally, the first two time points from each experimental condition were discarded from further analysis to omit contamination from transient effects of the hemodynamic response during classification.This correlation were chosen for further analysis.We did not include more vertices for the principal analysis to decrease the chances of overfitting.
For classification, we used the SVM implementation sklearn.svm.SVC with fixed regularization term  = 1 that is based on the libsvm library (Chang and Lin, 2011).This method was repeated for all possible splittings of training and test data sets using a leave-one-run-out cross-validation procedure to estimate mean prediction accuracies., Hubel, and Wiesel, 1975).This is the expected topography as depicted in (Adams, Sincich, and Horton, 2007;Adams and Horton, 2009).

Topography of ocular dominance columns
The blind spot is a further distinctive monocular region of V1 (Tootell, Hadjikhani, et al., 1998).
Due to the lack of photoreceptor cells on the optic disc of the retina where the optic nerve bundles and passes through, there is an oval area in V1 on the contralateral hemisphere that is solely "filled" by the response from the ipsilateral eye.In Figure 3A, there is a spatially extended response from the ipsilateral eye at the anterior end of the stimulated area (see cyan arrow in Figure 3A), which could be the blind spot representation on this hemisphere.However, due to the limited visual field in our experiment, we did not expect to have covered the blind spot region, which should be found at around 15 • eccentricity (Tootell, Hadjikhani, et al., 1998).Therefore, we assume that this response is of vascular origin or a response that was elicited by the border of the stimulus.).Only data points (n = 200) were used that were also selected for the decoding analysis.Note that we inverted the y-axis in C for consistency with A and B. Mean percent signal changes across cortical depth with all V1 data can be found in Supplementary Figure 6.Percent signal change curves from single participants can be found in Supplementary Figure 7.

GE-BOLD SE-BOLD VASO
We cannot exclude the possibility that some columns merged due to idiosyncrasies in local A double peak can be identified in C that is not visible in A and B, respectively.Red solid and dashed lines show the mean across participants from the first and second session, respectively.Black lines indicate the mean across participants and scanning sessions.The gray area demarcates the bootstrap 95% confidence interval (n = 1,000).In A-C, data were significantly different (p < 0.05) from a 50% chance level at each cortical depth.The p-value was determined by bootstrapping (n = 1,000) and corrected for multiple comparisons of individual layers (FDR correction using the Benjamini and Hochberg procedure).Prediction accuracy curves from single participants can be found in Supplementary Figure 8.
distribution  2 = (1−) and adding an upper bound of 3 to the number of samples exceeding the 387 test statistics (Burt et al., 2020).A corrected -value of < 0.05 was considered statistically significant.

388
Table 1 shows the results of the correlation analysis from all participants.

394
We used the same vertices that were included in the classification analysis after feature selec-395 tion.As expected, GE-BOLD signal changes were overall larger than SE-BOLD and VASO.Note that 396 signal changes for VASO, which has a negative relationship with CBV changes, were inverted for 397 visual purposes.

398
Across cortical depth, both GE-and SE-BOLD showed a steady increase toward the pial sur-399 face, most likely reflecting draining vein contributions to the signal (Polimeni, Fischl, et al., 2010; 400 Markuerkiaga, Barth, and Norris, 2016).The VASO signal profile was more restricted to GM and 401 shows a peak within GM.But an overall trend toward the pial surface could be seen as well.In Sup- analysis.An independent classification was performed for each cortical depth.Black lines indicate the mean across participants and sessions with the corresponding 95% bootstrap confidence interval.Red lines depict mean prediction accuracies from single sessions.Supplementary Figure 8 further illustrates prediction accuracies from single participants.
With all acquisition techniques, the eye-of-origin could be decoded with statistical significance at all cortical depths (chance level: 50% -value determined by bootstrapping).Among acquisition techniques, GE-BOLD showed the highest prediction accuracies.Furthermore, prediction accuracies increased toward the pial surface, mirroring the increase of univariate responses across cortical depth as shown in the previous section (Univariate contrasts across cortical depth).However, prediction accuracies did not show a steady increase compared to signal change profiles but saturated around mid-cortical depth.A similar behavior could be seen for SE-BOLD with an overall reduced level of prediction accuracies.
Compared to GE-and SE-BOLD, VASO showed a different cortical profile with a peak in deeper cortical layers.Note that this peak was a reliable finding reproduced across sessions (see red lines).
From a neuronal perspective, this is an expected finding since thalamocortical projections from the LGN primarily enter in layer 4C of V1 (Nieuwenhuys, Voogd, and Huijzen, 2008), which is located slightly below mid-cortical depth, see (Weber et al., 2008;Oga, Okamoto, and Fujita, 2016).The profile also showed a second peak in upper cortical layers.This might not be explainable from a neuronal point of view.However, acknowledging that VASO exploits a hemodynamic response and showed a modest signal increase in the univariate response of upper cortical layers depicted in Fig- ure 5C, the second peak in upper layers may also reflect remaining dependencies on macrovascular signal contributions.

Discussion
In this study, we mapped ODCs in V1 and decoded the eye-of-origin from raw fMRI time series in a group of human participants.While high-resolution fMRI was already used to analyze ODCs across cortical depth based on GE-BOLD (Hollander et al., 2021) and VASO measurements (Akbari et al., 2023), we provide a detailed analysis of the cortical depth dependence of monocular feedforward information captured with MVPA, which were measured with GE-BOLD, SE-BOLD, and VASO.
Our ODC in vivo maps (see Figure 3) robustly show a similar topography as found in histological studies (Adams, Sincich, and Horton, 2007;Adams and Horton, 2009).As expected, overall signal changes were highest for measurements with GE-BOLD, which led to increased repeatability of activation maps compared to SE-BOLD and VASO, see Figure 4 and Table 1.
The MVPA analysis showed that decoding of the eye-of-origin from sampled raw fMRI time series volumes is possible with all acquisition techniques at all cortical depths, see Figure 6.Our analysis revealed the best decoding performance for GE-BOLD, followed by SE-BOLD and VASO.
Interestingly, a similar increase of prediction accuracies toward the pial surface could be seen in GE-and SE-BOLD profiles, mirroring the increase of BOLD percent signal changes shown in Figure 5.
Similar profiles were derived for GE-BOLD (Vizioli et al., 2020;Hollander et al., 2021), which most probably do not reflect the profile of neural activity but are a sign of macrovascular draining vein contributions in the decoded BOLD signal.
VASO shows a distinct decoding profile compared to GE-and SE-BOLD with a peak in deeper cortical layers.This matches the expected pattern from neural processing, where the monocular response can be best discriminated in cortical layer 4C since it contains the highest segregation of left and right eye input.In contrast, layer 4B is targeted already after intracortical processing with increased combination of information from both eyes (binocular integration) (Wandell, 1995).In our analysis, we divided the cerebral cortex into 11 equidistantly spaced cortical depths.Regarding VASO signal.
Since MRI is sensitive to myelin content (Stüber et al., 2014;Trampel, Bazin, et al., 2019;Weiskopf et al., 2021), the stria of Gennari can be depicted in high-resolution anatomical images (Trampel, Ott, and Turner, 2011;Fracasso, van Veluw, et al., 2016).A precise anatomical localization of layer 4B is out of the scope of the current study but might be considered in future studies to verify the location of the found peak.Interestingly, high-resolution fMRI studies with GE-BOLD (Koopmans, Barth, and Norris, 2010) and VASO (Huber, Finn, et al., 2021) found that the highest cortical activation coincided with the stria of Gennari.However, while these studies considered binocular feedforward stimulation, our study focused on the processing of monocular information in V1, where layer 4C has a more prominent role than layer 4B.
It should be noted that we did not observe a peak in deeper layers in univariate VASO profiles as shown in Figure 5C, but see (Akbari et al., 2023).Possible explanations are that the laminar pattern is lost in univariate cortical profiles due to regional averaging of voxel responses.Furthermore, the higher sensitivity to noise contributions at the pial surface due to motion and brain pulsatility that lead to variation in voxel time series introduced by partial volume effects (Polimeni, Greve, et al., 2010;Pfaffenrot et al., 2021) might have further obscured the cortical profile at the single-voxel level.Higher resolution should decrease this effect, and interestingly, a recent study by Feinberg, Torrisi, et al., 2022 showed the activation pattern of feedforward binocular stimulated across cortical depth for GE-BOLD and VASO acquired with a nominal isotropic voxel size of around 0.4 mm, i.e., an 8-times smaller voxel volumes compared to the current study, which strongly resembled our cortical profile with a peak in deeper layers at the location, where layer 4C is expected.Therefore, MVPA methods in combination with VASO fMRI might be a promising technique to further increase the laminar specificity of VASO fMRI by being more robust to the aforementioned noise sources.
The results from Feinberg, Torrisi, et al., 2022 also showed a second peak in the upper layers, which was also seen in our decoding profile, see Figure 6C.This is not an expected finding from a neural perspective, when considering the monocular feedforward input.But as for Figure 6A and Figure 6B, this could also reflect information carried by draining veins.Therefore, it might be possible that both peaks have different origin, e.g., high spatial frequency information encoded for GE-BOLD are shown for V2 (A) and V3 (B) across cortical depth, respectively.Both areas were split in half based on retinotopy, with V2αand V3αbeing the half closer to V1.In A-B, data were significantly different ( < 0.05) from a 50% chance level at each cortical depth.The -value was determined by bootstrapping ( = 1, 000) and corrected for multiple comparisons of individual layers (FDR correction using the Benjamini and Hochberg procedure).These decoding performances in areas V2 and V3 cannot be attributed to responses at the columnar level and indicate that also decoding performances in V1 may not be exclusively due to responses at the columnar level.V2: secondary visual cortex, V3: tertiary visual cortex.at columnar level in deeper cortical layers and low frequency information due to the underlying vasculature in upper cortical layers.However, even if we cannot resolve the source of information that is used by the classifier, it is apparent that the VASO decoding profile does not follow the steady increase toward the pial surface as seen in GE-and SE-BOLD, which most probably is a reflection of macrovascular signal contributions.Future studies may show the origin of both peaks that were observed in this study.
The main decoding analysis was restricted to 200 features (vertices).It might be argued that the peak in deeper layers, though repeatable across sessions, might depend on the selected vertices.
Figure 7 presents cortical profiles of prediction accuracies, which were independently derived for different numbers of vertices [1, 2, … , 500].It can be seen that only a few voxels were necessary to decode the eye-of-origin, which was similarly found for orientation decoding (Haynes and Rees, 2005a).GE-and SE-BOLD showed a similar increase toward the pial surface for different numbers of included vertices.The peak for VASO measurements remained stable for various numbers of vertices.However, with a further increasing number of vertices, the cortical profiles became more similar to the profiles of GE-and SE-BOLD.We note that including more voxel can lead to overfitting and might have overshadowed the weak effect in deeper layers.This can also be seen in the univariate profiles shown in Supplementary Figure 6, which show overall decreased signal changes when more voxels are included due to the inclusion of noise.
The interpretation of the laminar profile is built on the assumption that the monocular feedforward information is exploited in V1, which is encoded at the fine-grained level of ODCs.Note that the larger monocular regions in V1, like the blind spot (Tootell, Hadjikhani, et al., 1998) and the temporal monocular crescent (Nasr, LaPierre, et al., 2020), were not covered in our experiment due to the limited field of view.However, we cannot exclude that other features besides ocularity might have contributed to the successful eye-of-origin decoding.Therefore, we conducted an additional analysis, in which we decoded the stimulated eye from cortical areas outside of V1 that are known not to be driven by monocular input.
Figure 8 shows cortical profiles of prediction accuracies from GE-BOLD data (200 vertices) sampled in secondary visual cortex (V2) and tertiary visual cortex (V3), respectively.V2 and V3 were further divided into two halves (α: half closer to V1, β: half further away from V1).The stimulated eye could be decoded in both V2 and V3 across cortical depth, but with overall decreased decoding performances compared to Figure 6A.Furthermore, a similar increase toward the pial surface was visible.Since no information about ocularity is expected from extrastriate cortex, the exploited fMRI signal also needs to contain other information that enables classification.Interestingly, decoding performance decreased in V2 when using data further away from V1.
That could be due to parts of V1 still being included when sampled close to the V1/V2 border.However, the opposite was seen in V3, where decoding performance increased again with growing distance to V1.Because V1/V2 and V2/V3 show mirror symmetry around areal borders representing the vertical meridian and horizontal meridian, respectively, this indicates that prediction accuracies became more prominent when data around the representation of the vertical meridian were included.It should be noted that parts of the lower vertical meridian might have been occluded due to the participant's nose when seeing through one eye (R. B. H. Tootell, personal communication).
This can lead to differences in induced activation in some parts of the visual field when stimulating the left or right eye, respectively, which is therefore not related to ODCs.This feature could also be seen in other ODC studies (Feinberg, Vu, and Beckett, 2018), which showed an extended area around the representation of the lower vertical meridian biased for one eye.
Next to differences caused by the variability in visibility of some parts of the visual field, we note that the used stimulus also differed in color and luminance between eyes that was not explicitly controlled for.This might have led to decodable information along the parvo-and magnocellular streams inside but also outside of V1 (Tootell and Nasr, 2017).For example, Supplementary Figure 1-Supplementary Figure 5 illustrate ODC maps from single participants, which generally show higher responses for the left eye, irrespective of eye dominance of single participants (eye dominance is stated in corresponding figure captions), which might be caused by remaining luminance differences between colors and therefore between eyes.Similar observations were made in an early fMRI decoding study, in which the eye-of-origin was decoded from a binocular rivalry stimulus (Haynes and Rees, 2005b).In binocular rivalry, the left and right eye receives incongruent stimuli, which were presented via anaglyph goggles.In that study, color filters were swapped between successive fMRI scanning runs in a control experiment.This resulted in decreased decoding performances in V1, whereas in extrastriate area V3 it stayed above chance level.From these results, it was concluded that performances in V1 were mostly based on ocularity information, while extrastriate areas V2 and V3 exploited more the color information in the stimulus.While not having the data to confirm these results in our experiment, we hypothesize that a similar effect drives the decodability in extrastriate areas as seen in Figure 8.
Our fMRI signal is therefore influenced by several biases that are not related to ocularity information.These biases will also lead to differences in the expected laminar profile.However, we emphasize that, compared to other decoding studies exploiting information encoded at the columnar level with a conventional resolution, we could map and visualize ODCs in all single participants.That means that fine-grained information at the spatial scale of ODCs was present and the dominant pattern in univariate activation maps (see Figure 3), which potentially was exploited by the linear classifier.Although it is probable that the classifier used multiple sources of information at different spatial scales, we are confident that the induced signal changes were mainly based on ocularity information, which is encoded at the level of ODCs with a peculiar laminar profile.
In VASO measurements, we exploit a CBV-weighted contrast that has a different temporal evolution compared to the BOLD response (Buxton, Wong, and Frank, 1998;Silva, Koretsky, and Duyn, 2007).More specifically, the CBV response has no initial dip, a shorter time-to-peak after stimulus onset, no poststimulus undershoot after stimulus offset, and needs more time to return to baseline.However, for the univariate analysis and the repeatability analysis, we processed data from all acquisition types with the same canonical HRF.As a control, we also analyzed the VASO data with a modified HRF that more closely resembled the CBV response's time evolution (data not shown), which only resulted in minor differences to the presented results.Note that we did not use an HRF model for the multivariate analysis, since analysis was based on the raw fMRI time series.
Our study showed that GE-BOLD, SE-BOLD, and VASO have different properties regarding the laminar specificity for the decoding of cortical feedforward information.For the first time, we used VASO in combination with MVPA to retrieve information from fine-grained cortical structures at the level of cortical layers.GE-BOLD is a very time-efficient acquisition method with larger SNR compared to SE-BOLD and VASO.This enables GE-BOLD to decode columnar information with high accuracy.However, the signal is weighted toward macrovascular signal contributions, limiting its capabilities to resolve information at the level of cortical layers.In comparison, VASO is a two-shot method, which limits its time efficiency.In addition, the BOLD correction in VASO is performed by a division operation, which enhances noise in the time series.This manifests in overall lower decoding accuracies for VASO.However, since the signal is more sensitive to microvascular contributions, information at the laminar level can be better localized.

Figure 1 .
Figure 1.Illustration of used stimuli.For visual stimulation, we used red (A) and green (B) random dot stereograms (RDSs) that were viewed through anaglyph goggles by participants, respectively.Stimuli were based on Nasr, Polimeni, and Tootell, 2016 and enabled full field of view visual stimulation of the left or right eye in separate experimental blocks.RDSs formed the percept of an 8 x 6 checkerboard during stimulation with independent sinusoidal movements in the horizontal direction of individual squares.
Visual stimuli consisted of a flickering (4 Hz) black-and-white radial checkerboard restricted to a clockwise/anticlockwise rotating wedge (angle: 30 • , temporal frequency: 1∕64 Hz) or expanding/contracting ring (temporal frequency: 1∕32 Hz) shown in separate runs.Each run presented 8.25 cycles of stimulation, with a baseline block of 12 s at the beginning and end of each run, in which a uniform gray background was shown.The mean luminance of the stimuli was set to 44 cd/m 2 .Throughout the run, participants were instructed to maintain fixation on a central point.No explicit task was given.

Figure 2 .
Figure2.Illustration of the GBB method.The method is used to enhance the alignment of cortical boundary surfaces based on an undistorted whole-brain anatomy to the cortical borders found in distorted functional images.A shows the temporal mean of the functional time series without task (GE-BOLD, 200 time points, subject 3) in coronal view that was acquired in the first session.B To enhance the GM/WM border and thereby increase the robustness of the proposed method, we weighted the temporal mean by its phase (see to Data preprocessing for detailed information) as usually done in susceptibility-weighted imaging methods.In C, the surfaces before (depicted in red) and after (depicted in green) alignment with the GBB method are presented.This technique is implemented in the GBB package (0.1.6,https://pypi.org/project/gbb/).Before alignment, surfaces were transformed into functional space via a rigid registration.The method improves spatial correspondence between the surfaces and the GM/WM boundaries observed in the functional images.GM: gray matter, WM: white matter.

Figure 3 .Figure 4 .
Figure 3. Representative maps of ocular dominance columns (ODCs).Thresholded activation maps (contrast: left eye > right eye) are shown for the left hemisphere of one representative participant (subject 1).Data were sampled at mid-cortical depth.In A, the contrast from GE-BOLD sessions (average across two sessions) is shown on the inflated surface.Several columns confined to V1 can be identified.The green arrow points to a small location (gray area) outside of the imaging field of view.B-C show the contrast from single GE-BOLD sessions on the flattened surface.The similar appearance of both maps illustrates the consistency of the columnar pattern across sessions conducted on different days.D-E show the contrast from SE-BOLD and VASO sessions (average across sessions).Due to the reduced number of slices, the area around the foveal representation was not covered (see the white dotted line that outlines the covered area).A similar ocular dominance pattern can be seen in all maps (see black dots with white outline for reference).Note that VASO has an inverted contrast compared to BOLD.The white line shows the representation of the vertical meridian (V1/V2 border) that was based on a separate retinotopy measurement.White asterisks indicate the location of the foveal representation.The black line in E shows a scale bar (5 mm).Maps from all participants can be found in Supplementary Figure 1-Supplementary Figure 5.

Figure 3
Figure3shows ocular dominance column maps (contrast: left eye > right eye) for a representative participant sampled at mid-cortical depth.Maps from single participants can be found in Supplementary Figure1-Supplementary Figure5.Figure3Ashows the average activation map across two GE-BOLD sessions.Some features can be seen that are expected from ODCs: (1) V1 shows a fine-scale pattern.(2) The pattern is constrained to area V1. (3) Around the approximate location of the horizontal meridian, columns are oriented more in parallel to both vertical meridians (V1/V2 border)(LeVay, Hubel, and Wiesel, 1975).This is the expected topography as depicted in(Adams,

Figure 5 .
Figure 5. Percent signal changes across cortical depth.Mean percent signal changes (contrast: left eye and right eye > baseline) for GE-BOLD (A), SE-BOLD (B), and VASO (C) are shown across cortical depth.Red solid and dashed lines show the mean across participants from the firstand second session, respectively.Black lines indicate the mean across participants and scanning sessions.The gray area demarcates the bootstrap 95% confidence interval (n = 1,000).Only data points (n = 200) were used that were also selected for the decoding analysis.Note that we inverted the y-axis in C for consistency with A and B. Mean percent signal changes across cortical depth with all V1 data can be found in Supplementary Figure6.Percent signal change curves from single participants can be found in Supplementary Figure7.

FigureFigure 4 Figure 6 .
Figure 3Cand Figure3Dshow the average activation maps across sessions for SE-BOLD and

389Univariate contrasts across cortical depth 390 Figure 5
Figure5shows the strength of cortical responses by plotting the percent signal changes of left and 391

402 plementary Figure 6 ,Figure 6
Figure6shows mean prediction accuracies across cortical depth from the pattern classification 413

Figure 6C ,Figure 7 .
Figure 6C, the peak in deeper layers corresponds to a relative cortical depth of around 73%.When compared to an exemplary laminar cytoarchitecture of V1 (Palomero-Gallagher and Zilles, 2019), this approximately corresponds to layer 4C, located below layer 4B, the location of the stria of Gennari.It is therefore intriguing to argue in favor of increased laminar specificity found in the

Table 1 . Repeatability of ODC maps across scanning sessions for single participants.
Spearman's rank correlation coefficients and corresponding p-values are shown to illustrate the consistency of activation maps (contrast: left eye > right eye) between scanning sessions for single participants.Only data from V1 sampled at mid-cortical depth were used.Statistical significance was determined by permutation testing (n = 10,000).Due to the spatial covariance of data from neighboring vertices, only randomly selected 10% of all data points were used for significance testing.
(Pedregosa et al., 2011)rformed by only considering time series data from locations within V1 that were present in the FOV of all functional sessions.Based on the training data, we further used an  -test implemented in the scikit-learn library (1.2.0, https://scikit-learn.org/)(Pedregosa et al., 2011), specifically sklearn.feature_selection.f_classif, to select the vertices whose time series strongest correlated with the experimental paradigm.The top 200 vertices with the highest