Improved sensitivity and microvascular weighting of 3T laminar fMRI with GE-BOLD using NORDIC and phase regression

Introduction: Functional MRI with spatial resolution in the submillimeter domain enables measurements of activa- tion across cortical layers in humans. This is valuable as diﬀerent types of cortical computations, e.g., feedforward versus feedback related activity, take place in diﬀerent cortical layers. Laminar fMRI studies have almost exclu- sively employed 7T scanners to overcome the reduced signal stability associated with small voxels. However, such systems are relatively rare and only a subset of those are clinically approved. In the present study, we examined if the feasibility of laminar fMRI at 3T could be improved by use of NORDIC denoising and phase regression. Methods: 5 healthy subjects were scanned on a Siemens MAGNETOM Prisma 3T scanner. To assess across-session reliability, each subject was scanned in 3–8 sessions on 3–4 consecutive days. A 3D gradient echo EPI (GE-EPI) sequence was used for BOLD acquisitions (voxel size 0.82 mm isotopic, TR = 2.2 s) using a block design ﬁnger tapping paradigm. NORDIC denoising was applied to the magnitude and phase time series to overcome limitations in temporal signal-to-noise ratio (tSNR) and the denoised phase time series were subsequently used to correct for large vein contamination through phase regression. Results and conclusion: NORDIC denoising resulted in tSNR values comparable to or higher than commonly observed at 7T. Layer-dependent activation proﬁles could thus be extracted robustly, within and across sessions, from regions of interest located in the hand knob of the primary motor cortex (M1). Phase regression led to substantially reduced superﬁcial bias in obtained layer proﬁles, although residual macrovascular contribution remained. We believe the present results support an improved feasibility of laminar fMRI at 3T.


Introduction
Functional MRI (fMRI) with spatial resolutions at the submillimeter scale is a rapidly growing field, motivated, in part, by the ability to resolve cortical layers noninvasively in human subjects. Activation patterns across distinct laminae constitute a fingerprint of feedforward and feedback related information ( Felleman and Van Essen, 1991 ), making the method, which we will refer to as laminar fMRI, a valuable asset for studying hierarchical information flow in the brain. Numerous studies already demonstrated its worth in areas related to, e.g., visual perception and attention ( Aitken et al., 2020 ;Kok et al., 2016 ;Liu et al., 2020 ;Muckli et al., 2015 ), motor control ( Huber et al., 2017 ;Persichetti et al., 2020 ), somatosensation ( Yu et al., 2019( Yu et al., , 2022, auditory perception Moerel et al., 2019 ), and language ( Sharoh et al., 2019 ). Laminar fMRI further holds promise ( Duvernoy et al., 1981 ), and the problem applies in particular to BOLD fMRI using GE-EPI readouts (GE-BOLD) due to its T2 * based contrast ( Yacoub et al., 2003 ). Sequences relying on T2 based contrast such as spin echo EPI (SE-EPI) ( Han et al., 2021 ;Koopmans and Yacoub, 2019 ) and b-SSFP ( Báez-Yánez et al., 2017 ;Liu et al., 2020 ), cerebral blood volume (CBV) based contrast such as VASO ( Huber et al., 2015 ;Lu et al., 2003 ) or cerebral blood flow (CBF) based contrast such as ASL ( Kashyap et al., 2021 ;Shao et al., 2021 ) are less sensitive to macrovasculature. Activation maps of these sequences thus have higher microvascular weighting resulting in improved co-localization with underlying neuronal activity. However, SE-EPI for instance is rarely free of T2 * effects due to long readout times ( Goense and Logothetis, 2006 ) and trade-offs with these alternatives to GE-BOLD generally include lower contrast-to-noise ratio (CNR), lower sampling efficiency and potentially less straightforward implementation. The "best " sequence thus does not exist and which one to choose depends on the specific research goal. GE-BOLD has nevertheless been the most frequently adopted sequence for laminar fMRI hitherto and will be the sequence of focus in the present study.
The vast majority of laminar fMRI studies employed ultrahigh field ( ≥ 7T) scanners. This is motivated by an increased SNR which allows for the high spatial resolution necessary for laminar fMRI studies ( Triantafyllou et al., 2005 ) in addition to an enhanced BOLD effect ( Ugurbil, 2014 ). Moreover, signal contributions from intravascular compartments are reduced as a result of the very short T2 * -value of blood at ultrahigh field ( Jochimsen et al., 2004 ;Uluda ǧ et al., 2009 ). The aforementioned sensitivity and macrovascular challenges are hence both alleviated with ultrahigh field systems ( Chaimow et al., 2018 ;Dumoulin et al., 2018 ;U ğurbil, 2021 ). The downside of a dependence on ultrahigh field is a dramatic decrease in the availability of laminar fMRI, as these systems are still relatively rare, and clinically approved ones, even more so. Successful GE-BOLD 3T implementations of laminar fMRI have been reported ( Kim and Ress, 2017 ;Koopmans et al., 2010 ;Markuerkiaga, Marques, Bains, et al., 2021 ;Puckett et al., 2016 ;Ress et al., 2007 ;Scheeringa et al., 2016Scheeringa et al., , 2022Taso et al., 2021 ;Wu et al., 2018 ), demonstrating that useful layer-dependent activation measures can be obtained outside of ultrahigh field applications. However, sequence and analysis strategies previously adopted to overcome the challenges at lower field strength, such as specialized head coils or large ROIs, may be incompatible in many settings, and widespread use of laminar fMRI at 3T has yet to emerge. To this end, a recently published denoising method, named NORDIC Vizioli et al., 2021 ), has proven effective at increasing the signal stability of magnitude and phase time series data. It is based on principal component analysis (PCA) applied to the full complex valued MRI dataset data with the aim of removing principal components that cannot be distinguished from zero-mean gaussian noise. Thermal noise fits this description and is the dominant source of noise for submillimeter voxels at 3T ( Triantafyllou et al., 2005 ). Initial reports on the effect of NORDIC indicate substantial tSNR increments without any significant spatial blurring being introduced Vizioli et al., 2021 ). It may thus help alleviate the challenge of reduced tSNR at lower field strengths and increase the utility of 3T laminar fMRI.
Accordingly, the main purpose of this study was to examine if the feasibility of 3T laminar fMRI with GE-BOLD could be improved using NORDIC denoising. Furthermore, we examined if increased microvascular weighting could be achieved through phase regression ( Menon, 2002 ), which aims to suppress signal from large veins. Several alternative vein correction approaches have been proposed (discussed in Huang et al., 2021 ;Huber et al., 2021c ;Kay et al., 2019 ;Koopmans and Yacoub, 2019 ;Stanley et al., 2020 ). We decided to test this particular method for three main reasons: (1) it has proved effective for selectively suppressing the response in voxels contaminated by large veins (pial veins and largest ascending veins) while maintaining high sensitivity ( Stanley et al., 2020 ); (2) It is applicable for a large range of laminar fMRI applications as it does not rely on any additional scans or equip-ment, and it does not, in principle, rely on paradigm design (although it has some CNR dependency); (3) its efficiency may elevate with the combination of NORDIC as a result of denoised phase time series. We scanned subjects during a finger tapping task and extracted layer profiles from regions of interest (ROIs) in the hand knob of the primary motor cortex (M1) which enabled comparison with results from the 7T laminar fMRI literature. Finally, we examined the reliability of this 3T setup both within and across sessions by scanning each subject on multiple days. The results could potentially help to increase the utility of laminar fMRI at 3T, which would be valuable for making laminar fMRI available to a much wider community.

Subjects
5 healthy right-handed subjects (1 female) with an age of 25-29 years were included in the study, approved by the regional ethics committee in Region Midt (Study ID 1-10-72-216-21). All subjects were carefully informed about the procedures and provided written consent.

Imaging protocol
Imaging was performed using a Siemens MAGNETOM Prisma 3 T scanner equipped with a standard 32Ch-receive head coil. Anatomical reference images were collected with a MP2RAGE sequence ( Marques et al., 2010 ) and parameters: voxel size = 0.9 mm isotropic, matrix size = 192 × 240 × 256, iPAT = 2, Partial Fourier = 6/8, TE = 2.87 ms, TR = 5000 ms, TI1 = 700 ms, FA1 = 4°, TI2 = 2500 ms, FA2 = 5°, echo spacing = 7.18 ms. Functional images were obtained with a 3D GE-EPI sequence  and parameters: voxel size = 0.82 mm isotropic, TR = 2200 ms, TE = 27 ms, iPAT = 3, Partial Fourier = 6/8 (zero-filling reconstruction), FA = 45 degrees, and matrix size = 176 × 176 × 26 where the 26 axial slices were aligned perpendicularly with respect to the surface of the hand knob in M1 ( Fig. 1 ). Adaptive combine reconstruction of magnitude and phase images was done using the vendor provided Prescan Normalize setting. Additionally, we used a SS-SI-VASO sequence  (voxel size = 0.82 mm isotropic, TR = 4739 ms, TE = 27 ms, iPAT = 3, Partial Fourier = 6/8, FA = 18 degrees, and matrix size = 176 × 176 × 26, volumes = 100) with the same readout strategy and resolution as the functional sequence to obtain T1-weighted images with good anatomical contrast, and distortions similar to the functional images . Matched distortions enabled high quality registration to the functional images without the need for distortion correction, and the T1-weighted contrast further aided nonlinear registration of the MP2RAGE image to EPI volumes. VASO images were only used to get the distortion-matched T1-weighted image and not for functional mapping or similar.

Experimental protocol
Each subject was scanned in 3-8 sessions on 3-4 consecutive days resulting in 23 total sessions, 3 of which were excluded (see Table S1 for detail). Note that one of these sessions was used in a supplementary analysis to assess the effect of motion (Fig. S1). The procedure was identical across all sessions to enable assessment of across-session consistency of laminar responses. The paradigm consisted of a single functional run with 35.2 s trials (17.6 s tapping and 17.6 s rest) repeated 30 times resulting in a total duration of ∼18 min and acquisition of 480 functional volumes (exclusive 4 dummy volumes). Tapping blocks included alternating finger tapping movements between the right index and middle fingers. The frequency of tapping for each finger was ∼2.5 Hz and was guided visually by a blinking marker projected onto a semi-transparent screen placed inside the MRI scanner, visible to the subject via a mirror mounted on the head coil. Before entering the scanner, subjects were The functional paradigm of each session consisted of ∼18 min (including functional localizer) block designed finger tapping (17.6 s tapping alternating with 17.6 s rest, repeated in 30 trials). (C) NORDIC denoising was applied to the functional time series data after DICOM-to-NIFTI conversion. (D) Motion corrected NORDIC and noNORDIC time series were then submitted to phase regression to correct for large vein effects. An example phase time series from a voxel which is located towards CSF, shown in green to the left, has a clear modulation with respect to the task and thus presumably contains a large vein. The linear fit between the phase time series and magnitude time series is subtracted from the original magnitude signal (red) leading to suppression of the paradigm modulation in the resulting "micro " time series (blue). The same is shown to the right for a voxel in gray matter where a clear task modulation can be observed in the magnitude time series only, resulting in a largely unchanged micro time series. E) The ROI of each subject was defined in native space of the first session by first anatomically locating the hand knob ( Yousry et al., 1997 ). The ROI was then manually drawn in a part of the knob where a functional localizer from the first session showed strong deep layer activation. Activation maps from the other sessions were then aligned to the space of the first session and thereby evaluated in the same ROI.
instructed to lie as still as possible, while keeping their gaze on the blinking marker in both tapping and resting blocks. Head motion was further minimized by tape across the forehead for haptic feedback ( Krause et al., 2019 ) and by placing inflatable padding around the head which could be manually adjusted by the subject to optimize comfort. Body motion was reduced by wrapping the subject in a heavy blanket. The duration of each scanning session was kept below 1 h.

Preprocessing
After DICOM-to-NIFTI conversion, functional and VASO time series were denoised with NORDIC  , https: //github.com/SteenMoeller/NORDIC _ Raw , commit 74999d6, downloaded 27,042,022) using the full complex valued data. In NORDIC, singular value decomposition (SVD) is performed on patches of functional time series data (represented as Casorati matrices), where all components with singular values below a certain threshold are removed. This threshold is determined as the average of the largest singular value across Monte Carlo-simulated thermal noise matrices. The variance associated with these random matrices is estimated directly from the measured data, either by an appended noise-only acquisition or from an area outside the brain without signal contributions (see Vizioli et al. 2021 for details). Here, the latter approach was employed. In addition, we employed more aggressive denoising by setting the option ARG.factor_error (FE) equal to 1.15 (default is 1), which scales the aforementioned thresh-old such that more components are removed. The more aggressive denoising comes with a risk of altering the temporal correlations in the data, but as shown in Figs. S2 and S3, additionally removed components appeared to be dominated by thermal noise.
The denoised functional magnitude time series were then motion corrected in SPM12 (Functional Imaging Laboratory, University College London, UK). To optimize the alignment around the hand knob of M1, a spatial weighting mask was applied with largest weights on the hand knob and progressively smaller weights towards the periphery of the brain. The phase images were converted into real and imaginary parts prior to reslicing using realignment parameters estimated from the magnitude images. These resliced (complex) images were converted back to phase values and unwrapped in the temporal direction. For the purpose of comparing scenarios with and without application of NORDIC, the estimated motion parameters were also applied to the same magnitude and phase time series without denoising (noNORDIC).
SS-SI-VASO provides 2 images per TR, one where the blood has been nulled by an inversion pulse, and a second where the blood has not been nulled . These images were first motion corrected separately due to their distinct contrasts, and then registered to the first image of the functional time series in a single interpolation step. The signal variability of the combined nulled and not-nulled time series was then used to compute a T1-weighted image with high anatomical contrast ( Beckett et al., 2020 ;Huber et al., 2017 ). Due to this T1-weighted VASO image theoretically being in the same distorted space as the functional BOLD EPI images, while having similar contrast as MP2RAGE, it was used as a "high-quality reference image " for non-linear registration of MP2RAGE to EPI-space. The reason for using MP2RAGE, even though the T1-weighted VASO image has sufficient contrast for segmentation, was that it enabled an automatic initial segmentation estimate used for guidance (see Methods Section 2.4.4 ), and facilitated across session alignment (see below). MP2RAGE to EPI-space registration was done by first estimating initial parameters for alignment of MP2RAGE to the T1-weighted VASO reference in ITK-SNAP ( Yushkevich et al., 2006 ), followed by rigid, affine and non-linear registration steps ( SyN algorithm) in ANTS ( Avants et al., 2011 ). The same procedure was used to align MP2RAGE images from each session to the MP2RAGE of the first session in each subject. These parameters were then used to align functional maps to the space of the first session in which analyses for all sessions could be performed (see Methods Section 2.4.4 ). We carefully inspected the outputs to ensure all registrations were of sufficient accuracy for laminar analysis.

Implementation of phase regression
Next, phase regression was applied according to the scheme described in previous reports ( Curtis et al., 2014 ;Menon, 2002 ;Stanley et al., 2020 ;Vicente et al., 2021 ) to correct for bias associated with large veins. The main principle is to remove the linear fit between magnitude and phase time series from the magnitude time series. This is based on the physical property that macrovascular directional flow changes result in a change in the phase signal, whereas microvascular flow has no common direction resulting in no net phase change (for details, see Menon 2002 ). In line with the terminology of previous work ( Menon, 2002 ), we will refer to the resulting time series as micro time series. This is not meant to suggest that these will be free of all macrovascular influence, but that they are more weighted towards the micro vasculature. Before phase regression, magnitude and phase time series were filtered on a voxel wise basis by first fitting a general linear model (GLM) with a 16'th order FIR-set modeling the paradigm, and 60 nuisance regressors (24 motion ( Friston et al., 1996 ), 20 RETROICOR ( Glover et al., 2000 ) and 16 high pass) and then subtracting the nuisance fit from the original time series. Prior to nuisance regression the design matrix containing the nuisance regressors was orthogonalized with respect to the paradigm. In the event of stimulus locked motion, this orthogonalization leaves paradigm effects in the phase regressor unchanged even after nuisance filtering. Voxel wise linear fit parameters between filtered magnitude and phase time series were estimated with a deming regression model ( Hall, 2022 ), which is necessary as it accounts for noise on both the x-and y-variables ( Cornbleet and Gochman, 1979 ), whereas standard Ordinary Least Squares (OLS) regression assumes noise free regressors. The fit needs to be conditioned by the relative noise levels between magnitude and phase data, which determines the angle of the residuals to the line of the best fit. To accommodate this, the voxelwise residual variance in the filtered magnitude and phase time series was estimated by further subtracting the fit of the 16'th order FIR-set to also remove the variance explained by the paradigm. The relative noise level was then given as the ratio between the standard deviation of the magnitude residuals and the standard deviation of the phase residuals.

Activation maps
GLMs were then computed using AFNI ( Cox, 1996 ) where design matrices consisted of paradigm regressors convolved with a canonical hemodynamic response function. A separate parameter estimate (beta) was estimated for each of the 30 tapping trials. To avoid "double dipping " ( Kriegeskorte et al., 2009 ), odd trials were used for functional localization of ROIs only (unless stated otherwise), whereas even trials were used in all further analyses. Moreover, all ROIs were defined based on odd trials from the first session only. Separate GLMs were estimated for each of the following time series; filtered magnitude time series with and without NORDIC, and micro time series with and without NORDIC, resulting in 4 statistical maps per session (called NORDICmagn, noNORDICmagn, NORDICmicro and noNORDICmicro, respec-tively). Time series and regressors were scaled to yield beta maps reflecting percent signal change with respect to baseline.

ROI definition, segmentation and generation of layer profiles
ROIs were manually defined in native space of the first session for each subject. Functional maps of all other sessions were aligned to the space of the first session (see Methods Section 2.4.1 ) and thereby evaluated in the same ROI. Each ROI was defined in the hand knob area of Brodmann area 4a (BA4a) similarly to previous layer-fMRI reports focusing on M1 and finger tapping ( Beckett et al., 2020 ;Chai et al., 2020 ;Han et al., 2021 ;Huber et al., 2017 ;Persichetti et al., 2020 ;Shao et al., 2021 ). First, the hand knob was visually identified from its characteristic omega-shape ( Yousry et al., 1997 ) and voxel selection was limited to the upper and lateral part of the knob corresponding to BA4a ( Huber et al., 2017 ). Gray matter (GM) boundaries were manually drawn on an upsampled grid (0.2 mm in plane resolution) in this area guided by an automated segmentation obtained with CAT12 ( http://www.neuro.unijena.de/cat/ ). A depth map comprising the relative cortical depth (values between 0 and 1) of voxels between these boundaries was then computed with LAYNII ( Huber et al., 2021c ) using the equivolume metric ( Waehnert et al., 2014 ). The point at which voxels were estimated to be 50/50 white matter (WM) and GM was located at a relative depth between 0.1 and 0.2. The equivalent point for GM and cerebrospinal fluid (CSF) was between 0.8 and 0.9. The hand knob is functionally divided into topographical digit representations . To identify the index-middle finger representation, which is expected to have relatively strong deep layer activation for the present tapping task, each ROI was defined in columns where the functional localizer (based on odd trials from the first session, see Methods Section 2.4.3 ) had strong activation in deep layers (between depths 0.3 and 0.4 to reduce partial voluming with WM and middle layers, respectively). Note that ROIs were defined as columns going roughly perpendicular onto the surface and without holes to make it somewhat anatomical meaningful. All ROIs spanned multiple slices and the total number of voxels within each ROI were in the range 77-113 (see Fig. S4 which also shows the depth distribution of ROI voxels).
Layer profiles were generated by plotting the beta of each ROI voxel as a function of its cortical depth. The depth of each voxel, that is, the correspondence between its center and the underlying depthmap in upsampled resolution, was inferred by Nearest Neighbor interpolation using the spm_sample_vol -function in SPM12. Beta-estimates from single voxels are quite noisy, so MATLAB's (Mathworks Inc.) movmeanfunction (window size = 0.2) was used to get the final layer profiles for each of the 15 test-betas which were averaged to obtain a single profile per session per data version. In order to run group analyses, without bias from specific depths being sampled differently across subjects, profiles were interpolated at 18 equally spaced depths between 0.075 and 0.925. This approach was chosen, over interpolation of the activation maps, in order to more directly interpolate between voxels sampled at similar depths.

Evaluation of sensitivity and macrovascular contamination
Sensitivity analysis was conducted by submitting the 15 betas to a second level analysis where voxel wise t-values were calculated as = ̄ ∕ √ where ̄ and s is the mean and standard deviation, respectively, across beta estimates (even trials only), and n = 15 beta estimates. This approach is in contrast to the typical calculation of t-values from a first level fMRI analysis where residuals from the GLM is used directly. It was chosen to acknowledge that the degrees of freedom potentially may be affected by NORDIC denoising. By generating t-values from single trial beta-estimates, which are assumed to be independent, this potential concern is largely removed.
A quantitative measure of macrovascular contamination was obtained by fitting a line to the across-session mean layer profiles of each subject. The slopes of these lines represent the degree of bias towards the pial surface, previously referred to as inverse specificity ( Beckett et al., 2020 ). Two-way repeated measures ANOVA with factors NORDICversion (NORDIC/noNORDIC) and veinCorrection (micro/magnitude) was employed to evaluate statistical significance of differences in slopes across the 4 combinations of NORDIC/noNORDIC and magnitude/micro data versions.
Data is presented as mean ± standard deviation unless otherwise stated. The significance level was set to = 0.05.

Data and code availability
The data are available upon request and signed data sharing and data processing agreements as part of Aarhus University data sharing regulations. Analysis code will be made publicly available on figshare: 10.6084/m9.figshare.19925483.v1

Results
The effect of NORDIC on signal stability was quantified by computing an average tSNR-value of all ROI-voxels for each session which was averaged to one tSNR-value per subject. This was done using detrended but otherwise unfiltered motion corrected magnitude time series with and without NORDIC. The mean tSNR was ∼3 times larger with NORDIC (40.46 ± 5.32), compared to without NORDIC (12.32 ± 1.39) (paired ttest, degrees of freedom = 4, p < 0.001). The effect of such increased temporal stability on activation maps is illustrated in an example subject in Fig. 2 A. In the high tSNR case, the hand knob is clearly activated and accompanied by much less false positive activation outside of M1. Furthermore, laminar profiles with and without NORDIC in an example subject are depicted in Fig. 2 B showing that the consistency across trials is clearly facilitated by denoising. Across-trial mean profiles are similar whereas the variability around the mean is notably increased when denoising is omitted. Fig. 2 C depicts absolute single voxel t-values of all sessions in each subject. The first row shows that t-values are generally elevated with NORDIC compared to without NORDIC for both "magn " and "micro ", suggesting improved sensitivity after denoising. Phase regression, on the other hand, has been reported to come along with reduced sensitivity ( Stanley et al., 2020 ) in which case we would expect larger t-values for "magn " compared to "micro ". The second row indeed shows this to be the case when NORDIC is applied. Conversely, noNORDIC t-values seem to be distributed fairly evenly around the identity line. Taken together, this suggests that phase regression reduces sensitivity to a larger degree when NORDIC is applied compared to when it is not. However, as will be addressed in later sections, part of the effect is likely explained by reduced sensitivity towards spatially unspecific draining veins and is thus desired.
Thermal noise is not only reduced by the application of NORDIC, but also as a consequence of averaging across trials and voxels. To investigate the impact of NORDIC on laminar profiles in the presence of averaging, we generated profiles from subsamples of the total number of trials ( N = 30). This is shown in Fig. 3 for different subsampling factors in an example subject. Without NORDIC, notable variability across subsampled profiles already occurs when 15 trials are averaged. Conversely, NORDIC profiles stay similar even when only 5 trials are averaged. This suggests that thermal noise still had a substantial effect on noNORDIC profiles (average of 15 trials) at the present ROI size (90 voxels in total for this subject). Fig. S5 shows results when ROI size is halved by discarding every second voxel (based on depth) to illustrate how the benefit of NORDIC increases when ROI size becomes smaller. To further illustrate the thermal noise contribution at the present ROI size, we computed periodograms of modelled non-white noise and thermal noise. Based on these we found that without NORDIC, thermal noise still seems to have a significant contribution even after averaging across voxels (see Fig. S6). Fig. 4 A shows hand knob activation maps of the 4 different data versions in an example slice from the first session of each subject. Comparison of NORDICmagn and NORDICmicro activation maps suggests that voxels with large response suppression after phase regression are predominantly located where large veins are expected, i.e., towards superficial layers and CSF, supporting the notion that phase regression mainly suppresses macrovascular signal contributions. While the majority of voxels likely to contain large veins (i.e., superficial/CSF voxels with very strong activation in magnitude maps) appear to have strongly suppressed responses after phase regression, there are obvious cases (blue arrows points to examples of this) where the method left residual macrovascular signal, highlighting that although the method reduces macrovascular contributions, it is not perfect. Since the effectiveness of phase regression is CNR-dependent, we suspected its performance to be reduced when NORDIC was omitted, which is supported when comparing noNORDICmagn and noNORDICmicro activation maps, where differences are barely visible and most of the unwanted signal towards large veins remains unsuppressed. This implies that NORDICmicro profiles should have smaller superficial bias compared to noNORDICmicro profiles. This was confirmed using ANOVA with slopes of linear fits to each subject's across-session mean profile as the dependent variable, which revealed a significant NORDICversion by veinCorrection interaction ( p = 0.036). That is, the difference in slopes between magnitudeand micro data versions was larger when NORDIC was applied compared to when it was not (see average slope across subjects for each data version in Fig. 4 B). This effect is further visualized in Fig. 4 C, which illustrates that noNORDICmicro across-subject mean profiles have steeper gradients from WM to CSF compared to NORDICmicro (Bonferroni corrected Post hoc comparison, p = 0.044). Note that the difference in percent signal change between group mean magnitude and micro profiles is on average largest in CSF and decreases progressively towards WM, implying that bias towards superficial layers is reduced in the micro profiles. Fig. 5 A displays NORDICmagn and NORDICmicro profiles for each session plotted together for each subject. It accompanies the acrosssubject mean results in Fig. 4 C by showing how the reduced bias towards the surface replicates across sessions. Fig. 5 A further illustrates the robustness of laminar profiles across sessions. To quantify acrosssession robustness of NORDICmicro laminar profiles, we computed the coefficient of variation (CV), i.e., the standard deviation over sessions divided by the mean (done separately for each layer and then averaged as in  ). The average session-CV across subjects was 26.20 % ± 7.78 %. Similarly, Fig. 5 B depicts the robustness of NORDICmicro profiles across trials for each session and subject. Laminar patterns of activation appear highly consistent across trials, with trial profiles closely following the mean profiles of individual sessions. This was quantified as the trial-CV, i.e., the standard deviation across trials divided by the mean, which was done for each layer and session separately and then averaged. The average trial-CV across subjects was 27.53 % ± 14.67 %.

Discussion
The main purpose of the present study was to test if the feasibility of laminar fMRI at 3T with GE-BOLD could be improved through sensitivity gains obtained by NORDIC denoising. 3T fMRI with submillimeter spatial resolution suffers from poor signal stability which partly explains why nearly all existing laminar fMRI studies were performed at 7T systems. The results presented here demonstrate how this challenge can be overcome with NORDIC denoising, allowing for reliable detection of laminar activation both within and across sessions at 3T. This can be achieved with standard MR hardware, without application of explicit spatial smoothing kernels which degrade spatial specificity, and in a single functional run of 18 min (including functional localizer) making it highly practical. Additionally, we tested phase regression ( Menon, 2002 ) as a method to reduce superficial bias caused by large veins and hence improve microvascular weighting. We found evidence for significant improvements in obtained layer profiles, although some macrovascular contribution likely remained. . Colored profiles represent single trials with across-trial mean profiles plotted on top in black. Mean profiles are similar across the two data versions, whereas variability around the mean is notably higher without NORDIC denoising. Error bars represent SEM across trials ( N = 15). Asterisks denote that half the trials from the first session was used for ROI localization, see Methods Section 2.4.3 C) Absolute t-values of ROI voxels for all sessions in each subject. The first row depicts this for NORDIC versus noNORDIC. The second row depicts this for "magn " versus "micro ".

Sensitivity improvement
A prerequisite for sensitivity, in terms of detecting functional responses, is a stable underlying signal which can be quantified as tSNR.
Here it was computed from motion corrected and detrended NORDICmagn data before motion and RETROICOR regressors were filtered out to avoid uneven comparison with previous ultrahigh field ( ≥ 7T) submillimeter studies. The average tSNR within ROIs was 40.46 ± 5.32 across subjects, which is comparable to, if not higher than what is commonly reported for submillimeter 7T GE-BOLD setups with similar or coarser resolution compared to the present study ( Aitken et al., 2020 ;Beckett et al., 2020 ;Rua et al., 2017 ;Stanley et al., 2020 ;Zaretskaya et al., 2020 ). For example,  , Rua et al. (2017) and Aitken et al. (2020) reported tSNR/ 3 values of 46.4, 25.3, and 24.5, respectively, in motor and visual cortex ROIs (73.6 in the present study). Note that tSNR depends on factors such as EPI-readout strategy (2D versus 3D acquisitions), acquisition parameters (repetition time, flip angle, etc.), image SNR (and thus voxel size), and ROI selection (all brain voxels versus specific subsets). Choices in sequence design, preprocessing and ROI selection may thus partly explain how the present 3T setup has comparable tSNR values to 7T observations. However, the main reason is likely the application of NORDIC de- Fig. 3. Impact of NORDIC on laminar profiles. Green and purple profiles represent NORDICmagn and noNORDICmagn, respectively. Rows and columns represent different sessions and subsampling factors, respectively. Profiles for the first column were generated by averaging across all trials. For the second column, one profile is generated by averaging across trials 1:2:30 (every other trial starting from the first) and the second profile by averaging across trials 2:2:30 (every other trial starting from the second), and similarly for the third and fourth columns. Without NORDIC noticeable variability across subsampled profiles is observable already when 15 trials are averaged. After NORDIC, profiles remain relatively stable also when only 5 trials are averaged, suggesting that denoising still has a considerable effect on laminar profiles even after thermal noise has been reduced by averaging across voxels and trials.
noising, which is apparent from the fact that in our sample it increased tSNR by a factor of ∼3 on average. Importantly, NORDIC increases tSNR without adding noteworthy spatial smoothness to the images (see Fig.  S7 and Dowdle et al. 2021. In one of the original NORDIC articles, the developers found that denoising on average increased tSNR by a factor of more than 2 . The larger observed effect of NORDIC in the present data might partly be explained by the fact that methods for determining the threshold under which components with lower singular values were discarded, differed across studies (see Methods Section 2.4.1 , Figs. S2 and S3). Also, in contrast to Vizioli et al. (2021) , we used the version of NORDIC that works in image space as opposed to raw data space, without the need for an appended noise volume and g-factor map. A complementary explanation is that our runs were significantly longer than those used by Vizioli et al. (2021) resulting in greater data redundancy which theoretically should lead to better separation of signal and noise and thus more effective denoising. Finally, the degree to which thermal noise is dominant scales inversely with field strength ( Triantafyllou et al., 2005 ) which loosely speaking implies that more noise is present to begin with and available for removal at 3T.

Improved microvascular weighting
GE-EPI is the most frequently applied sequence in layer-fMRI studies due to its superior statistical efficiency and ease of implementation. However, its weighting towards macrovasculature is well established ( Menon, 2012 ;Turner, 2002 ;Uluda ǧ et al., 2009 ;Yacoub et al., 2003 ) and needs to be dealt with for laminar purposes where spatial specificity is of the essence. The superficial bias that follows from large vein sensitivity was clearly present in the NORDICmagn profiles which indeed possessed the characteristic positive gradient from WM to CSF ( Figs. 4 C and 5 A) consistently found with GE-BOLD ( Aitken et al., 2020 ;de Hollander et al., 2021 ;Han et al., 2021 ;Huber et al., 2015 ;Kok et al., 2016 ;Shao et al., 2021 ;Stanley et al., 2020 ). Also, activation maps revealed higher percent signal changes towards the surface and in CSF which adds to the evidence for large vein contamination ( Fig. 4 A). Phase regression has recently shown promise as a method to reduce the macrovascular sensitivity of GE-BOLD in 7T submillimeter settings to obtain microvascular specificity comparable to regular SE-BOLD ( Stanley et al., 2020, but see Han et al. 2021 for SE-sequence with improved microvascular weighting). If phase regression predominantly reduces signal orig- Fig. 4. Effect of phase regression on superficial bias (A) Hand knob activation maps for the 4 different data versions in example slices from each subject. With NORDIC, most voxels likely to contain large veins have strongly suppressed responses after phase regression. However, residual macrovascular contributions are still left in NORDICmicro activation maps as pointed out by the blue arrows. Without NORDIC, the difference before and after phase regression appears small and most unwanted signal towards large veins remains unsuppressed. (B) Average of slopes of linear fits to each subject's across-session mean profile for the 4 different data versions (colored profiles in C). Colored dots represent datapoints of each subject. Error bars represent SEM across subjects ( N = 5). (C) Magnitude and micro profiles plotted together for all subjects both with and without NORDIC. Colored profiles represent single subjects (red = magnitude and blue = micro) with across-subject mean profiles plotted on top in black (dashed = magnitude and solid = micro). Straight lines indicate linear fits to each of the 4 mean profiles. The difference in slopes between magnitude and micro profiles is significantly larger with NORDIC, suggesting a more effective phase regression compared to without NORDIC.
inating from large veins, we should see larger signal suppression towards the surface and in CSF. We did indeed find evidence for this in line with previous observations ( Curtis et al., 2014 ;Menon, 2002 ;Stanley et al., 2020 ;Vicente et al., 2021 ); the difference in magnitude between NORDICmagn and NORDICmicro profiles progressively increased from WM to CSF ( Fig. 4 C) and the same effect is observable directly in activation maps ( Fig. 4 A). The degree of superficial bias was estimated as the slope of the linear fit to across-session mean profiles of each subject inspired by previous efforts ( Beckett et al., 2020 ;De Martino et al., 2013 ;Huber et al., 2017 ). The extent of the effect is then reflected in the decrease of the slope from NORDICmagn to NORDICmicro which is shown in Fig. 4 B. Although the slope might be a simplified measure of inverse specificity, as also pointed out by Beckett et al. (2020) , it does illustrate that phase regression substantially reduced the gradient from WM to CSF which is generally believed to originate from large veins. Interestingly, the effect of phase regression, as measured by its Fig. 5. The effect of phase regression on superficial bias and reproducibility across sessions and trials A) NORDICmagn and NORDICmicro layerprofiles plotted together for each subject. Colored profiles represent single sessions (red = magnitude and blue = micro) with across-session mean profiles plotted on top in black (dashed = magnitude and solid = micro). Profiles have seemingly smaller superficial bias after phase regression and appear consistent across sessions. B) NORDICmicro layerprofiles of all sessions and subjects. Colored profiles represent single trials with across-trial mean profiles plotted on top in black. Individual trial profiles seem to closely follow the mean, implying robust laminar patterns of activation within sessions. Error bars = SEM across trials ( N = 15). Asterisks denote that half the trials from the first session was used for ROI localization, see Methods Section 2.4.3 impact on slopes, was significantly larger when NORDIC was applied ( Fig. 4 B and 4 C), illustrating its dependency on CNR of both the magnitude and phase time series. This is supported by a larger reduction in t-values after phase regression for NORDIC compared to noNORDIC ( Fig. 2 C), which we interpret to mainly reflect a reduced sensitivity towards macrovasculature. On the other hand, this also illustrates a limitation of the method, i.e., that superficial bias will deviate across experimental settings with varying noise levels, which needs to be considered when interpreting results.
Having said the above, several lines of evidence indicate that nonnegligible contributions from macrovascular sources were still in play. First, even though the group layer profile for NORDICmicro tends to flatten towards CSF as compared to the NORDICmagn profile which keeps rising, the peak is located at the GM/CSF border, whereas sequences proven to have higher spatial specificity peaks within the GM of M1 during finger tapping (VASO/VAPER: ( Beckett et al., 2020 ;Chai et al., 2020 ;Huber et al., 2017 ;Persichetti et al., 2020 ), SE/GRASE: ( Beckett et al., 2020 ;Han et al., 2021 ), ASL: ( Shao et al., 2021 )). A peak towards CSF after phase regression was also observed for the visual cortex in Stanley et al. (2020) , whereas neural activity is expected to peak in middle layers during visual stimulation ( Liu et al., 2020 ). Second, NORDICmicro activation maps still contain CSF voxels with strong percent signal changes ( Fig. 4 A). Third, the lateral end of the hand knob (BA4a) in M1 has typically been associated with a double peak response when examined with more microvascular weighted sequences at 7T ( Beckett et al., 2020 ;Chai et al., 2020 ;Guidi et al., 2020 ;Huber et al., 2017 ;Shao et al., 2021 ). This feature has been used as a hallmark for spatial specificity and has accordingly been less robustly observed with GE-BOLD where peaks become less distinct as a result of signal drainage caused by intracortical ascending veins and pial veins ( Huber et al., 2017 ;Persichetti et al., 2020 ). In the present study, 3 of 5 subjects consistently had either a double peak response or a superficial peak with a "shoulder " in the deep layers and the group mean profile might best be described as the latter. This is in line with GE-BOLD results from 7T where the double peak feature in M1 during finger tapping was found in 1/3 of subjects  and the group profile had a superficial peak with a "shoulder " in the deep layers ( Huber et al., 2017 , see also ( Beckett et al., 2020 ;Guidi et al., 2016Guidi et al., , 2020Huber et al., 2015 ;Persichetti et al., 2020 ;Shao et al., 2021 ). We find it promising that the results comply with expectations from 7T studies, especially when considering the inherently larger contribution of large vessels at 3T. However, the lack of clear double peak features across subjects, regardless of field strength, highlights the general issue of macrovascular contamination in GE-BOLD, which might be reduced by phase regression but clearly not eliminated fully. It should further be mentioned that although we do believe that the lack of two distinct peaks is explained by macrovascular influence, it may also partly be explained by the fact that we extract signal from relatively large ROIs spanning several slices ( Han et al., 2021 ;Pais-Roldán et al., 2020 ). Moreover, not all subjects have double peak profiles even with non-GE-BOLD sequences ( Beckett et al., 2020 ;Shao et al., 2021 ), thus we cannot rule out that more distinct peaks would emerge in the group profile with different subjects. Finally, it is not completely clear how well the microvascular signal is preserved in GM voxels which also contain large veins. To properly separate micro and macrovascular signal components the method assumes temporal uncoupling between hemodynamic responses of the two compartments ( Stanley et al., 2020 ). Considering that temporal delays between the two types of vasculature are on the order of about 1-3 s ( Kay et al., 2020 ), it may be too subtle to be resolved with TRs commonly employed within the field.
Blood related susceptibility changes around vessels giving rise to extravascular BOLD signal changes are known to be observable even in voxels remote from the vessel, especially for large pial veins Li et al., 2012 ;Moerel et al., 2018 ). It is mainly this so-called "blooming effect " in addition to intravascular BOLD sig-nals from pial veins that are expected to be suppressed with phase regression at high spatial resolution ( Stanley et al., 2020 ). Phase regression might also suppress signal from the largest intracortical ascending veins to some extent ( Klassen and Menon, 2005 ;Stanley et al., 2020 ). However, the majority of these are expected to be too small to produce measurable phase modulations and will thus not be affected by phase regression. Such veins hence represent a likely explanation for some of the remaining macrovascular contamination. Promising methods to reduce the component of superficial bias originating from intracortical ascending veins have been proposed. For example, Markuerkiaga et al. (2016 used spatial deconvolution to unmix the signal leakage from deep to superficial layers (see also Havlicek and Uluda ğ 2020, Heinzle et al. 2016, Marquardt et al. 2018. This method could potentially be combined with phase regression as these techniques have complementary strengths and weaknesses. In summary, the results presented here support that phase regression substantially reduces superficial bias and the associated decrease in sensitivity likely reflects a desired feature of reduced macrovascular weighting. It comes, however, with some limitations that need to be considered when selecting a tool to handle the problem of veins which is mandatory in virtually any laminar fMRI application. We opted to test this particular method due to its practicality and expected synergy with NORDIC (see Introduction). We would like to note that multiple alternative deveining methods have shown great promise and could be considered dependent on study design and aim. A thorough description of such methods is out of the scope of this work, but have been nicely discussed elsewhere (see for example Huang et al. 2021, Huber et al., 2021c, Kay et al. 2019, Koopmans and Yacoub 2019, Stanley et al. 2020 and references therein).

Across-session and across trial reproducibility
Each subject was scanned in multiple sessions on multiple days to assess across-session reproducibility of observed layer profiles after NORDIC and phase regression. To examine the result, we depicted the mean profiles of all sessions for each subject in the same plot ( Fig. 5 A). Visual assessment indicates that individual session profiles followed the across-session mean profiles with reasonable consistency, both regarding the laminar pattern and response magnitude. This is quantified as an average session-CV for NORDICmicro of 26.20 % ± 7.78 % across subjects, which is comparable to the across-days CV-estimate observed for VASO and BOLD (25 % and 30 %, respectively) in a 9.4T laminar fMRI finger tapping study . For reference, nonlaminar supramillimeter studies found CV's of 24 % and ∼28 % (means across subjects) for the variability of response magnitudes during motor tasks in M1 ROIs across 3 sessions at 3T ( Tjandra et al., 2005 ) and 7T , respectively. M1 laminar fMRI studies that scanned the same participants in multiple sessions and where associated layerprofiles of individual sessions were reported are sparse making comparison difficult. However, this was done in Chai et al. (2020) who also used a finger tapping paradigm and focused the analysis onto ROIs from the hand knob, although with a so-called VAPER sequence rather than GE-BOLD. No quantitative estimate was presented, but visual comparison with profiles from that study further supports that an acrosssession reliability in agreement with previous M1 studies, across field strengths and spatial resolutions, could be achieved despite the challenges of submillimeter fMRI at 3T.
NORDICmicro profiles are depicted in Fig. 5 B for each experimental session and shows that also at the level of single trials, profiles seem reproducible. This is quantified as a trial-CV of 27.53 % ± 14.67 % across subjects for NORDICmicro. CV-estimates for single trials are to our knowledge not reported in the laminar fMRI literature, but for reference, the CV across runs was 25 % for both VASO and BOLD in the previously mentioned 9.4T study .

Previous 3T submillimeter studies
The present demonstration of laminar fMRI being feasible at 3T with GE-BOLD is not the first; several previous studies have successfully obtained layer-dependent measures of activation at 3T ( Kim and Ress, 2017 ;Koopmans et al., 2010 ;Puckett et al., 2016 ;Ress et al., 2007 ;Scheeringa et al., 2016Scheeringa et al., , 2022Taso et al., 2021 ;Wu et al., 2018 ). For example,  demonstrated that by averaging responses across patches of cortex to suppress thermal noise, similar CNR profiles could be obtained at low field strength compared to those obtained at 7T. This strategy was for example employed in Scheeringa et al. (2016) to demonstrate how electroencephalographic signals could be related to layer-dependent BOLD responses in early visual cortex. Additional strategies that has been utilized to deal with the reduced signal stability at non-ultrahigh field strength include acquisitions with long TRs ( Koopmans et al., 2010 ), use of specialized coils with improved SNR around a region of interest ( Ress et al., 2007 ;Wu et al., 2018 ), and use of long scan durations ( Puckett et al., 2016 ). The benefit of NORDIC will depend on the degree to which thermal noise is already suppressed by such other means. For example, our results suggest that at the ROI sizes used here, NORDIC is beneficiary in terms of denoising profiles and boosting sensitivity as considerable thermal noise remain after averaging across voxels and trials (see Fig. 3 , Figs. S5 and S6). It will be of interest in future studies to further test how the effect of NORDIC on, e.g., laminar profiles, is altered for larger ROIs, more trials, etc. Notably, even in cases where physiological noise becomes dominant due to averaging, the higher single voxel tSNR afforded by NORDIC may still be beneficiary, e.g., for improving microvascular weighting with phase regression (see Discussion Section 4.2 ). In relation to this, exciting results from two recent conference abstracts by Huber and colleagues indicate that 3T laminar fMRI signals with excellent spatial specificity can be obtained with CBV-weighted VASO albeit with a tradeoff of reduced CNR ( Huber et al., 2021a ;Huber et al., 2021b ).
Based on the above, we believe the present results add to the contributions of previous 3T laminar fMRI studies, by showing that with the implementation of NORDIC, reliable submillimeter activation maps can be obtained with high statistical efficiency without the need for long scan durations, extensive averaging, specialized hardware etc.

Choice of NORDIC parameter
Changing the FE argument to 1.15 (default = 1) in NORDIC, as done here, effectively increases the number of removed time series components (see Methods Section 2.4.1 ). This leads to a dramatic increase in tSNR, but does come with the risk of removing signals of interest or altering temporal correlations. For the present setup, additionally removed components appeared to be dominated by thermal noise as assessed by temporal correlation matrices and Durbin-Watson statistics (Fig. S2). The periodogram of the difference time series (real part of noNORDIC minus NORDIC complex time series) was flat across the frequency range as expected for thermal noise, but it did contain a dip at the task frequency (Fig. S6). This may suggest that some signal containing components were removed, which is supported by an observed tendency towards a reduction in the response magnitude of group-averaged layer profiles after denoising (Fig. S3). However, in our case, the effect on response magnitude appears to be non-meaningful compared with error associated with estimation of group-mean responses (Fig. S3). For these reasons, and given the CNR dependency related to the effectiveness of phase regression, we found it optimal to adjust this NORDIC parameter in the present case. However, the magnitude of the effect likely depends on CNR-related experimental settings, such as sequence, acquisition parameters, field strength, whether the paradigm is block designed versus event related, etc. Alteration of the parameter should thus be done with care and might not even be necessary considering that the average tSNR across sessions and subjects was 22.89 ± 2.74 when NORDIC was applied with default settings (Fig. S2). Further testing of how well NORDIC removes noise and leaves the signal and its temporal structure intact across different CNR levels, input parameters, experimental settings etc. is required, ideally through simulation studies with known ground truths.

Limitations and future work
We opted to prioritize multiple test-retest scans on the same subjects to examine reliability of the current 3T setup. However, we acknowledge that the small sample size of only 5 subjects and 20 total sessions precludes rigorous statistical analysis and the reported reliability measures, as well as the comparison with such estimates from, e.g.,  which were also based on small samples, should be interpreted in that light. Future test-retest studies with larger sample sizes would thus be beneficial to further test our suggestion that layerdependent activation can be reliably measured at 3T. Such studies are to our knowledge currently lacking in the field as a whole, including ultrahigh field, and would be valuable as the validity of any method ultimately depends on its reliability.
Our results support a reduced macrovascular influence through phase regression, i.e., we found a reduced bias of laminar profiles towards CSF and activation maps were associated with suppressed activation of voxels likely to contain large vessels. We take this as one line of evidence for increased microvascular weighting. A limitation in our assessment of the method, however, is that our experimental settings do not allow for quantification of its effect on the functional point spread function ( Engel et al., 1997 ) which is a critical measure for high resolution fMRI. Further studies are thus required for such an evaluation, for example using visual paradigms as outlined previously (see for example, Chaimow et al. 2018, Fracasso et al. 2021, Shmuel et al. 2007 ). Moreover, despite phase regression improvements, the 3T GE-BOLD setup presented here is still associated with macrovascular contamination which is also true for GE-BOLD at higher field strengths. Given that every laminar fMRI study has vastly different requirements regarding sensitivity, spatial specificity, temporal resolution etc., future work at 3T could ideally test alternative sequences known to excel in areas where GE-BOLD is challenged. An example of a candidate sequence is SS-SI-VASO which already showed great promise at 3T ( Huber et al., 2021a ;Huber et al., 2021b ) and might be desired over GE-BOLD when spatial specificity is prioritized over sensitivity. The focus could, along this line, be directed towards alternative postprocessing methods for deveining, with TDM ( Kay et al., 2020 ) being an exciting suitor.
The flip angle of 45 degrees used for functional acquisitions was not optimized for sensitivity, i.e., it was roughly twice the size of the Ernst Angle (as pointed out by one of our reviewers). This means that reported measures of tSNR, stability across profiles, t-scores etc. might be lower than what could be obtained using the Ernst Angle. Conversely, the larger flip angle might have had a beneficiary effect in terms of suppressing CSF signal and thus reducing superficial bias . Future studies are needed to assess how sensitivity and macrovascular influence are modified by an altered flip angle.
Finally, it should be mentioned that certain compromises were necessary in order to achieve submillimeter resolution functional images with a reasonable temporal resolution: (1) The FOV was small (26 axial slices of 0.82 mm thickness) which in some cases led to fold-over artifacts visible in the magnitude images. Although this did not seem to have a big effect on the time series data in the present study, it should be of priority to mitigate the artifact by making sure the imaging slab extends outside the brain which comes at the cost of reduced brain coverage. Alternatively, the FOV can be increased in a tradeoff with longer TRs; (2) Partial Fourier reconstruction was employed for faster coverage of Kspace but comes with a reduction in the effective spatial resolution and it might introduce image artifacts in locations of high frequency phase changes ( Haacke et al., 1991 ). Partial Fourier was implemented with the zero-filling approach in line with previous studies using NORDIC  and phase-regression ( Stanley et al., 2020 ). More advanced methods such as POCS ( Haacke et al., 1991 ) might be beneficiary in order to minimize spatial blurring. Alternatively, Partial Fourier could be omitted entirely at the cost of temporal resolution. In future studies, advanced parallel imaging techniques such as 3D CAIPIRINHA acceleration and multiband excitation may be implemented to increase the coverage and temporal resolution of submillimeter fMRI at 3T.

Conclusion
We used GE-BOLD fMRI to test the feasibility of laminar fMRI at 3T with NORDIC and phase regression. The results suggest that SNR limitations usually accompanied with submillimeter isotropic voxels at 3T can be overcome with NORDIC denoising. As a result, layer-dependent BOLD responses could be detected reliably within and across sessions with high statistical efficiency. Importantly, this was achievable with standard MR hardware, without application of explicit spatial smoothing kernels, and in a single functional run, making it highly practical. As expected, we observed a strong superficial bias in magnitude layer profiles which was substantially reduced after phase regression, although some macrovascular contributions remained. We believe the present 3T results, coupled with the vast amount of research going into improving microvascular weighting, will help making laminar fMRI available to a much wider community.

Data and code availability statement
The data are available upon request and signed data sharing and data processing agreements as part of Aarhus University data sharing regulations. Analysis code will be made publicly available on figshare: 10.6084/m9.figshare.19925483.v1

Declaration of Competing Interest
The authors declare no competing financial interests.

Data availability
The data are available upon request and signed data sharing and data processing agreements as part of Aarhus University data sharing regulations.