Cortical laminar resting‐state signal fluctuations scale with the hypercapnic blood oxygenation level‐dependent response

Abstract Calibrated functional magnetic resonance imaging can remove unwanted sources of signal variability in the blood oxygenation level‐dependent (BOLD) response. This is achieved by scaling, using information from a perfusion‐sensitive scan during a purely vascular challenge, typically induced by a gas manipulation or a breath‐hold task. In this work, we seek for a validation of the use of the resting‐state fluctuation amplitude (RSFA) as a scaling factor to remove vascular contributions from the BOLD response. Given the peculiarity of depth‐dependent vascularization in gray matter, BOLD and vascular space occupancy (VASO) data were acquired at submillimeter resolution and averaged across cortical laminae. RSFA from the primary motor cortex was, thus, compared to the amplitude of hypercapnia‐induced signal changes (tSDhc) and with the M factor of the Davis model on a laminar level. High linear correlations were observed for RSFA and tSDhc (R2 = 0.92 ± 0.06) and somewhat reduced for RSFA and M (R2 = 0.62 ± 0.19). Laminar profiles of RSFA‐normalized BOLD signal changes yielded good agreement with corresponding VASO profiles. Overall, this suggests that RSFA contains strong vascular components and is also modulated by baseline quantities contained in the M factor. We conclude that RSFA may replace the scaling factor tSDhc for normalizing the laminar BOLD response.

conclude that RSFA may replace the scaling factor tSD hc for normalizing the laminar BOLD response.

K E Y W O R D S
7 T-fMRI, calibrated fMRI, hypercapnia, laminar fMRI, resting-state fMRI, VASO

| INTRODUCTION
Functional magnetic resonance imaging (fMRI) based on echo-planar imaging (EPI) with gradient-recalled echoes (GRE) is currently the most widespread technique to look at human brain function. However, its baseline-dependent nature makes data comparisons across participants, brain areas, and brain states challenging. Apart from true differences in neuronal activity, significant differences in the blood oxygenation leveldependent (BOLD) responses could be attributable to factors that are difficult to control for, such as caffeine intake, time of the day, age, baseline oxygenation, vascularization density, or baseline cerebral blood flow (CBF) (Krieger et al., 2014;Whittaker, Driver, Bright, & Murphy, 2016). The spatial signature of these effects becomes particularly clear at laminar resolutions (Polimeni, Fischl, Greve, & Wald, 2010;Yen, Zhao, & Kim, 2012), such that the so-called physiological noise (Krüger & Glover, 2001) is strongly linked to the baseline cerebral blood volume (CBV) (Koopmans, Barth, Orzada, & Norris, 2011).
Calibrated fMRI attempts to solve this issue by scaling the taskrelated BOLD response by one obtained with a purely vascular challenge, typically hypercapnia. Combined with a biophysical model, this has been used to extract information on (relative) changes of the cerebral metabolic rate of oxygen consumption (CMR O2 ) (Davis, Kwong, Weisskoff, & Rosen, 1998;Gauthier & Hoge, 2013;Hoge et al., 1999). This normalization is believed to result in a quantity that is more directly coupled to the neuronal response than unscaled BOLD signal changes. The underlying rationale is that the oxygen uptake needed for adenosine triphosphate (ATP) production happens at the closest capillary bed. Consequently, it should be more spatially specific to the site of neural activity than the BOLD signal changes, which are driven by the venous vasculature (Huber, Uluda g, & Möller, 2019). The method has proven to account for a high degree of variability in the BOLD response and, therefore, found some application (Blockley, Griffeth, Simon, & Buxton, 2013;Pike, 2012).
Most calibrated fMRI studies suffer from two main drawbacks: (a) They rely on a gas manipulation scan (apart from the functional paradigm). (b) They require the acquisition of an additional non-BOLD functional contrast for quantifying the vascular response. The vascular challenge has been traditionally taken to be a gas manipulation inducing hypercapnia, under the assumption that it does not induce changes in oxidative metabolism. This isometabolic assumption has not yet been robustly proven, and previous studies have reported a decrease Zappe, Uluda g, Oeltermann, U gurbil, & Logothetis, 2008), an increase (Jones, Berwick, Hewson-Stoate, Gias, & Mayhew, 2005), or no change (Chen & Pike, 2010;Jain et al., 2011) in CMR O2 . The metabolic effect of hypercapnia is believed to be dependent on the carbon dioxide (CO 2 ) content of the inspired gas, with higher concentration inducing higher oxidative metabolism changes (Jones et al., 2005;Zappe et al., 2008). The setup needed for the gas manipulation usually involves a facial mask or a mouth piece and a nasal clip, which causes discomfort to the participant. Additionally, the inflow of gas may be perceived as unpleasant (e.g., a dry mouth and throat) and may lead to dizziness if the CO 2 concentration of the gas mixture is high (e.g., above 5%). Such setup-specific problems can be circumvented by using a breath-hold task (Kastrup, Li, Glover, & Moseley, 1999). However, this comes at the expense of reduced signal quality (e.g., due to enhanced task-related head motion) and reproducibility (because of a reduced response amplitude). Moreover, respiratory manipulations might be inapplicable in noncompliant subject populations, such as children, or might be too demanding for patients with pulmonary or cardiac disease (Moreton, Dani, Goutcher, O'Hare, & Muir, 2016) or for elderlies. Finally, the additionally required functional contrast, such as CBF recorded with arterial spin labeling (ASL) (Alsop et al., 2015;Detre, Leigh, Williams, & Koretsky, 1992;Lorenz, Mildner, Schlumm, & Möller, 2018;Mildner et al., 2003) or CBV recorded with vascular space occupancy (VASO) techniques (Hua, Jones, Qin, & van Zijl, 2013;Huber et al., 2014b;Lu, Golay, Pekar, & van Zijl, 2003;Lu, Hua, & van Zijl, 2013), is penalized by a reduced sensitivity compared to the GRE-BOLD response. Consequently, this has a crucial impact on the quality of the CMR O2 estimation (Huber et al., 2019).
Given these obstacles, measures obtained with so-called restingstate fMRI have been proposed as alternative normalization approaches.
Nevertheless, most of the calibrated BOLD studies still rely on gas manipulation or breath-hold scans rather than on resting statebased approaches. This might be due to the fact that several different scaling factors have been proposed and most of them are missing a clear validation. Given that the band of interest is the low-frequency one, which is considered to be more closely coupled to spontaneous neural activity and, thus, used for connectivity studies (Birn, 2012;Murphy, Birn, & Bandettini, 2013), a mixture of neuronal and vascular contributions is expected. This "confounding" neural contribution represents one of the biggest impediments to the use of resting-state fluctuations as a scaling factor (Lipp, Murphy, Caseras, & Wise, 2015;Liu, 2013). Moreover, the choice of the repetition time, TR, impacts the amount of aliased cardiac and respiration frequencies into the low-frequency band (Viessmann, Möller, & Jezzard, 2017;Viessmann, Möller, & Jezzard, 2019;Wise, Ide, Poulin, & Tracey, 2004).
So far, these methods have only been applied at relatively low spatial resolution. Their applicability to high-resolution fMRI, as required for obtaining depth-dependent information, has not yet been investigated. An increase in resolution is limited, in general, by the achieved temporal signal-to-noise ratio (tSNR). A more specific potential problem is related to the interpretation of the BOLD response because traditional assumptions about neurovascular coupling mechanisms may not be valid locally on a submillimeter spatial scale. Voxels of a dimension on the order of the cortical thickness (or above) in a given brain region contain similar mixtures of arterial, capillary, and venous blood. On this spatial scale, it can be assumed that BOLD signal changes are modulated by consistent dynamic changes in CBF, CBV, and CMR O2 (Buxton, Uluda g, Dubowitz, & Liu, 2004). On a submillimeter scale, however, different voxels of the same brain region may have different total CBV fractions or different portions of arterial, capillary, and venous blood. For example, voxels located at the pial surface should contain relatively large fractions of arterial and venous blood but should be mostly devoid of capillary blood. In such voxels, significant changes in the deoxyhemoglobin concentration, [dHb], may be observed if they contain veins draining an (upstream) activated area (Markuerkiaga, Barth, & Norris, 2016). However, this local BOLD signal change may not have a colocalized CMR O2 change.
One goal of the current work was to investigate the use of RSFA indices for normalizing BOLD signal changes recorded at submillimeter spatial resolution on a laminar basis. This was done by comparing the obtained results to other metrics recorded with gas-calibration experiments in a separate session. Given that the human neocortex is organized into six layers, each of which can be treated-to some extent-as a functional unit, averaging was performed along laminae rather than across the full patch of cortex. Compared to standard-resolution scans, laminar fMRI yields additional information about an approximate location within the cortical ribbon where BOLD signal changes occur. This helps to better disentangle contributions from the capillary bed and those from the downstream vasculature (Huber et al., 2015;Kim & Kim, 2010;Koopmans, Barth, & Norris, 2010). For further validation, a second goal was a direct comparison of laminar CBV changes and simultaneously recorded BOLD signal changes with and without RSFA-based normalization. Ten right-handed healthy volunteers (6 males, mean age: 256 ± 4 years) with no history of neurological disorders participated in the first part of the study after giving informed written consent. The experimental procedures had been previously approved by the Ethics Committee at the Medical Faculty of the University of Leipzig (Reg.-No. 273-14-25082014).
Participants were asked to refrain from coffee and alcohol intake on the day of the experiment. A physician was present during each session to monitor physiological recordings during the breathing manipulation task.
The gas mixture used to induce hypercapnia contained 5% CO 2 , 21% oxygen (O 2 ), and 74% nitrogen (N 2 ). It was delivered to the participant through a mouthpiece connected by a tube to a gas wall socket.
The mouthpiece was adjusted to fit tightly to the participant's face in order to avoid inflow of room air. A nose clip was used to prevent breathing from the nasal cavity. The inflow of gas was manually adjusted to a rate of approximately 18 L/min (standard temperature and pressure). Normocapnia was restored by disconnecting the tube from the wall socket, thereby letting room air flow into the tube. A second tube connected to the mouth piece was used for directing the exhaled air outside the scanner bore. Inhaled and exhaled gases were prevented from mixing via valves placed in the interior of the mouth piece. Heart beat and respiration timecourses were recorded using the scanner's physiological monitoring unit (pulse sensor and respiratory belt). The end-tidal partial pressures of CO 2 (P ET CO 2 ) and O 2 (P ET O 2 ) were recorded using an MP150 system (BIOPAC Systems, Goleta, CA). Written verbal instructions were projected onto a screen, which could be viewed by the participant via a mirror mounted on the RF coil. In case of visual deficiencies, MRI-compatible glasses were provided to the participant.

| Functional paradigm
The paradigm consisted of two sessions, each 15 min long. The first session was a block-design hypercapnia task, followed by a resting-state session. The hypercapnia task consisted of an alternation of breathing room air and the CO 2 -enriched gas mixture. After an initial 2-min block of room air, the gas mixture was administered for two 3-min blocks, separated by an equally long block of room air. A final 4-min block of room air breathing was added for baseline localization. For the resting-state scan, a fixation image was projected onto the screen, and the participants were asked either to fixate it or to stay with the eyes closed for the whole duration of the experiment and not to fall asleep.

| Imaging pulse sequence
Previously acquired MP2RAGE (magnetization prepared two rapid gradient echoes) (Marques et al., 2010) images were inspected in order to localize the primary motor cortex (M1) of each individual participant. Positioning the three-dimensional (3D) slab for fMRI followed previously established procedures (Guidi, Huber, Lampe, Gauthier, & Möller, 2016;Huber et al., 2015;Huber et al., 2018). In particular, a slice orientation perpendicular to the cortical surface achieves a sufficiently high resolution along the cortical depth while the slice thickness can be increased depending on the individual anatomy to improve the SNR. However, even with a sufficiently wide slab, it is typically not possible to align the slices perpendicularly to M1 on both sides at the same time. Therefore, slice positioning was always optimized for the hand knob area of the left M1, which was exclusively selected for the subsequent data analysis. A slice-saturation slabinversion (SS-SI) VASO sequence (Huber et al., 2014b) with a hybrid 3D EPI readout (Poser, Koopmans, Witzel, Wald, & Barth, 2010) was used for the acquisition. Previous work has shown that the 3D readout is beneficial for sub-millimeter applications (Huber et al., 2018).

| Image preprocessing and layering
All timeseries were corrected for rigid-body motion using the function "Realign: Estimate and Reslice" in SPM12 (http://www.fil.ion.ucl.ac. uk/spm/). The recorded respiratory and cardiac traces were used for denoising the resting-state timecourses using the Analysis of Functional NeuroImages (AFNI) (Cox, 1996) implementation of RETRO-ICOR (Glover, Li, & Ress, 2000). All fMRI timeseries were linearly detrended to remove low-frequency signal drifts.
Anatomical references for layering were obtained from resized T 1weighted EPI maps generated from the functional gas-manipulation and resting-state timeseries. The T 1 -weighting was derived from the original SS-SI-VASO timeseries consisting of alternating EPI slabs acquired with and without a preceding inversion pulse, that is, prior to the CBV/BOLD contrast splitting (Guidi et al., 2016). Since the signal difference of the interleaved T 1 -weighted VASO and BOLD acquisitions is dominated by longitudinal relaxation, this "anatomical" gray matter (GM)/cerebrospinal fluid (CSF) contrast in EPI space is comparable to the contrast on T 1 -weighted MP2RAGE uniform ("UNI") images. A representative example is shown in Figure 1.
For layering purposes, the original in-plane matrix size of 132 × 44 was resized to 528 × 176. Consistent with most segmentation software packages, this upsampling was done to obtain smooth laminae without angularity limitations in voxel space. Gray matter/ CSF and GM/white matter (WM) borders were then manually drawn on the slice with the best through-plane orientation. The region of interest (ROI) was defined as the portion of GM that was (a) well aligned to the slice orientation across the cortical depth and that (b) additionally showed sufficient activation on BOLD and VASO maps. In this ROI, 15 laminae were grown as in previous work (Guidi et al., 2016;Huber et al., 2015) employing C++ and Object-oriented Development Interface for NMR (ODIN) libraries (Jochimsen & von Mengershausen, 2004). Briefly, the algorithm follows an equivolume approach that preserves the volume fraction of each lamina in cortical segments ( Figure 1) and accurately resembles the arrangement of anatomical layers (Waehnert et al., 2014). Throughout this work, we use the term "laminae" rather than "layers" to stress the fact that they correspond to certain depths from the cortical surface, whereas they do not a priori correspond to single histological layers.

| Quantitative analyses
Following Kannurpatti, Rypma, and Biswal (2012), the resting-state fluctuation amplitudes of the unfiltered timeseries (RSFA full ) as well as of their low-frequency (RSFA low ) and high-frequency (RSFA high ) portions were taken to be the corresponding temporal standard deviation (tSD rs ). The low-frequency and high-frequency timeseries included frequencies in the bands between 0.01 and 0.1 Hz and F I G U R E 1 Schematic illustrating the four steps involved in the calculation of depth-dependent profiles. The slice with the best through-plane orientation was chosen (1), and the gray matter (GM)/ cerebrospinal fluid (CSF) surface and the gray matter (GM)/white matter (WM) surface were manually delineated (2). Equivolume laminae were then generated between these two borders (3). Finally, the quantities of interest were averaged along each individual lamina and plotted as a function of depth (4). Here, the laminar profile of tSD hc measured in a representative participant is shown as an example between 0.1 and 0.15 Hz, respectively. The unfiltered timeseries included frequencies in the entire range between 0.01 and 0.15 Hz as defined by scan duration and TR. All RSFA values were calculated for each lamina (Figure 2). Similarly, the amplitude of hypercapniainduced signal changes was taken to be the temporal standard deviation of the corresponding timeseries (tSD hc ) (Kannurpatti & Biswal, 2008).
The calculation of the calibration factor, M, followed the approach presented in previous work (Guidi et al., 2016). In brief, the Davis model (Davis et al., 1998;Hoge et al., 1999) was modified to express BOLD signal changes in terms of CBV and CMR O2 and to account for the CBF term using Grubb's relationship (Grubb, Raichle, Eichling, & Ter-Pogossian, 1974). The original Davis model can be written as ΔS = (S − S 0 )/S 0 is the relative BOLD signal change from the baseline level (indicated by an index "0"); f = CBF/CBF 0 and r = CMR O2 /CMR O2, 0 are the normalized cerebral blood flow and oxidative metabolic rate, respectively; and α t = 0.38 is the Grubb exponent relating the total blood volume (CBV t ) to CBF. β is a magnetic field-specific constant describing the coupling between the reversible transverse relaxation rate, R 0 2 , induced by inhomogeneous external fields and [dHb]. In general, it depends on the vessel size but may be approximated as β ≈ 1 at 7 T, where intravascular signal contributions become negligible at typical TE values (Bright, Croal, Blockley, & Bulte, 2019;Kida, Kennan, Rothman, Behar, & Hyder, 2000;Martindale, Kennerley, Johnston, Zheng, & Mayhew, 2008). Consequently, the change in R 0 2 is approximately linear in [dHb]. M represents the maximum possible BOLD signal change and depends on TE and the baseline levels of venous blood volume (CBV v ) and [dHb] in addition to further parameters (e.g., vessel size, brain region, magnetic field strength), which may be lumped together into a proportionality constant κ: In the modified Davis model, the dependency on f in Equation (1) is replaced by a dependency on v t = CBV t /CBV t, 0 , yielding (Guidi et al., 2016): Here, α v = 0.2 is a modified Grubb exponent relating CBV v to CBF (Chen & Pike, 2009). For mild hypercapnia assumed to be isometabolic, r ≈ 1, and the evoked BOLD signal change can be expressed as: If laminar BOLD and CBV changes are recorded during a hypercapnic challenge, M is easily obtained with Equation (4) (Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012) and AFNI (Cox, 1996). Subsequently, The same type of comparison as in Figure 5 in a participant without a clear double-peak pattern in the laminar blood oxygenation level-dependent (BOLD) response. with R 2 = 0.92 ± 0.06 (mean ± 1 SD across participants) for RSFA full , 0.94 ± 0.06 for RSFA low , and 0.92 ± 0.08 for RSFA high (Table 1). For the correlations of RSFA and M, R 2 was 0.62 ± 0.19 (p < .05 in 8 out of 9 participants) for RSFA full , 0.60 ± 0.26 (p < .05 in 7 out of 9 participants) for RSFA low , and 0.55 ± 0.26 (p < .05 in 7 out of 9 participants) for RSFA high (Table 2). Of note, applying RETROICOR correction did not substantially affect these correlations. This is consistent with previous research demonstrating that corrections methods attempting to minimize contributions of physiological noise become less relevant at high (e.g., submillimeter) resolution, where noise is increasingly dominated by thermal fluctuations Murphy, Bodurka, & Bandettini, 2007).

| Comparison of RSFA-normalized BOLD and CBV laminar profiles
Selected maps of tapping-induced percent BOLD and VASO signal changes and of the RSFA extracted from the resting-state and from the residuals of the tapping data as well as corresponding laminar profiles in M1 are shown in Figures 5 and 6. The laminar CBV profile in Figure 5b shows a clear indication of a double-peak pattern, consistent with previous reports (Huber et al., 2015;Huber et al., 2017).
The corresponding laminar BOLD profile demonstrates signal amplification in superficial laminae due to venous drainage (Turner, 2002). A pronounced amplitude increase towards the pial surface is also evident in both RSFA profiles. Figure 5c shows RSFA-normalized laminar BOLD profiles in direct comparison to the CBV profile from Figure 5b indicating a slightly improved depiction of the peak in upper layers.
We note, however, that laminar fMRI responses can be quite variable across participants. As such, while the double peak was robustly visible in previous VASO data, only one third of the participants showed the same feature in the BOLD response . To further exemplify the working principle of the RSFA method, we illustrate a corresponding case without a double-peak BOLD pattern in  Figure 6c).

| DISCUSSION
Nearly perfect correlations of RSFA and tSD hc were found for all examined frequency bands. This confirms previous results obtained on a voxel-wise level at standard fMRI resolution (Kannurpatti & Biswal, 2008). Our correlation strengths were somewhat higher than those reported by Kannurpatti et al. (2012). This may be attributed to the fact that linear regression was performed on a laminar basis in our work, that is, on more local level and with substantial blurring across laminae due to the relatively high degree of interpolation in this direction. There were no significant differences in the correlation strengths when comparing the high-and low-frequency bands, indicating similar scaling along the cortical depth. This assumption is further corroborated by the good agreement of laminar RSFA-normalized BOLD profiles and VASO profiles (Figure 5c).
Significant positive correlations were also found between RSFA and M. As for the correlations with tSD hc , similar correlation strengths were obtained for the different frequency bands in the majority of the participants. M can be considered in our work to be only modulated by subject-specific parameters, since sequence-specific parameters were identical in all acquisitions (see Equation (2)). This is not the case for RSFA, which does not entirely depend on the baseline levels but also on dynamic changes of CBF and CBV. Thus, we cannot exclude that lower correlations (relative to the above results for tSD hc ) might be partially due to the different physiological modulators of RSFA with respect to M. However, reduced correlations are also expected due to the variability of the M parameter, which is computed by a division of relatively subtle BOLD and CBV changes and is, hence, inherently noisier than a pure contrast.
The low-and high-frequency bands as defined here showed similar scaling as a function of depth. The low-frequency band (0.01-0.1 Hz) is usually considered to be the one reflecting the relevant neuronal activity (Biswal, Yetkin, Haughton, & Hyde, 1995;Margulies et al., 2010). However, it is obvious from our laminar profiles ( Figure 2) that even this band is dominated by contributions from the superficial draining vasculature, which does not correspond to the precise location of neuronal activity (Turner, 2002). While this does not rule out the possibility that the neuronal component is present in this band, it is likely to add some confounds to it (Markuerkiaga et al., 2016). The high frequency band (0.1-0.15 Hz) selected in our analysis is narrower than the traditional high-frequency band (0.1-0.25 Hz), whose definition is based on typical fMRI studies with TR ≈ 2 s. Similar to the results observed for the low-frequency band, the highfrequency band is also dominated by the pial surface, therefore, yielding almost identical correlations. Given the spatial nature of pial vein contaminations, it is reasonable to expect that the correlation would be disrupted as soon as its relative contribution becomes negligible, and the shape of the depth-dependent RSFA starts flattening out. This effect was visible at frequencies above 0.15 Hz (Guidi et al., 2017) and is also consistent with the profiles shown in Figures 5b and 6b.
Based on our results, RSFA seems to provide a good substitute for the amplitude of hypercapnia-induced signal changes to scale the BOLD response. Regarding the assumption of an isometabolic hypercapnia challenge, RSFA and tSD hc are likely to be similarly driven by reactivity features of the underlying vasculature (both following a relation similar to Equation (1)) and, therefore, look like scaled versions of each other.
Despite the good agreement of RSFA and M, a similar statement cannot be made in this case. This is due to their different physiological origins, with M depending solely on baseline parameters as long as it is accurately measured (Griffeth & Buxton, 2011). In fact, a strict correlation of RSFA and M might be not as straightforward as it seems (Lu, Hutchison, Xu, & Rypma, 2011). The calculation of M relies on several assumptions, whose validity is called into question, especially at high resolution. Such potential pitfalls include (among others) the validity of a Grubb-like relation between CBF and CBV with a suitable coupling exponent α, the assumption of an isometabolic hypercapnia challenge, or the validity of Fick's principle on a laminar basis. More extended discussions of these issues have been published recently (Guidi et al., 2016;Hua et al., 2019;Huber et al., 2019). For example, simulations assuming different values for α and β indicate that these variations primarily produce a scaling effect, whereas the laminar profiles were largely preserved (see Supplementary Information in [Guidi et al., 2016]). This seems to indicate that these factors should not critically affect our main conclusions.
The application of fMRI signal normalization along the cortical depth may further be associated with limited interpretability in a more general way: • We have shown in earlier work (Huber et al., 2015) that-in the process of normalization using depth-dependent physiological parameters-activity in some cortical layers may be undesirably underestimated. In particular, higher values of, for example, cerebrovascular reactivity (CVR), CBV 0 , or CBV v,0 in upper cortical layers may lead to artificially reduced normalized fMRI signal changes in these layers, despite the fact that the microvascularrelated response has a similar magnitude (see also figure 8 from Huber et al. (2015)).
• When the fMRI signal is normalized, the resulting activity measure depends on a number of estimated parameters as opposed to nonnormalized activity measures. Due to nonlinear error propagation of multiple parameters during the normalization, this might result in a noisier normalized activity measure.
• In the current work, we investigated the possibility of fMRI signal To further evaluate the quality of laminar profiles obtained after RFSA-based scaling, we used results obtained with VASO (i.e., a measure of the CBV response) as a reference of high spatial specificity. In general, the RFSA-normalized laminar BOLD profiles agreed well with the VASO results achieving a subtle improvement in the depiction of the expected double-peak contour (Figures 5c and 6c). This suggests that the normalization accounts for venous bias inherent to the BOLD response with respect to signal amplification. It is to note that quite similar results were obtained with both RSFA estimates (i.e., from the task-based time series and from an additional resting-state time series) indicating that sufficient information can be extracted directly from the task-based fMRI time series. However, some deviation remained, which was most pronounced at the pial surface pointing to the fact that voxel-wise RSFA-normalization does not account for local signal leakage.
Recent modeling results indicate that the baseline CBV variation due to the arrangement of ascending veins rather than venous draining alone contributes to the amplitude increase of the laminar BOLD response (Havlicek & Uluda g, 2020;Markuerkiaga et al., 2016).
The RSFA normalization as described here is a linear scaling approach that may be considered to correct-at first order-for local signal amplification in layers with higher CVR, However, such scaling does not account for signal leakage due to venous drainage from the location of distant layers. In particular, the laminar point-spread function of the BOLD response is not entirely defined by the local baseline CBV but depends nonlinearly on multiple physiological factors and on the level of activity (Havlicek & Uluda g, 2020;Huber et al., 2014a).
Future development will be necessary to incorporate higher-order normalization of signal leakage. As such, an integration of RSFA factors as proposed here into laminar-deconvolution models (Havlicek & Uluda g, 2020;Heinzle, Koopmans, den Ouden, Raman, & Stephan, 2016;Markuerkiaga et al., 2016;Merola & Weiskopf, 2018) might be an interesting topic of future research. Such models currently rely on assumptions on cortical depth-dependent fMRI reactivity (e.g., from ex-vivo data), which are not easily generalizable across brain areas beyond the primary visual cortex (V1) (Marquardt, Schneider, Gulban, Ivanov, & Uluda g, 2018). In this context, the RSFA method might provide a data-driven means for venous bias correction in laminar fMRI across brain areas and participants.
Other researchers tried to account for venous bias through linear postprocessing strategies, for example, by regressing out a linear slope of the laminar GRE-BOLD profile as a zeroth-order correction (Fracasso, Petridou, & Dumoulin, 2014) or by a decomposition of the profile into a constant and a linear term (Gau, Bazin, Trampel, Turner, & Noppeney, 2020). Similar to the RSFA method, these concepts are based on some form of linear scaling. Unlike the RSFA method, however, they are limited by the fact that the correction factors are not determined form a temporally orthogonal fMRI signal, which might introduce some level of circularity into the analysis.
Moreover, these methods cannot account for scaling effects that have nonconstant slopes. We may, therefore, speculate that the RSFA approach could be useful for accounting for additional sources of venous scaling compared to linear slope methods.
As another strategy, Muckli et al. (2015) proposed to account for venous scaling bias by means of constraining the fMRI activity interpretation to statistical measures of "classification accuracy." Unlike the BOLD signal magnitude, this measure is inversely proportional to the signal variance, which in turn is dependent on CVR. Thus, the measure of 'classification accuracy' should be inherently weighted by a CBVdependent scale factor, which is comparable to RSFA-based normalization. Lawrence et al. (2018) used a similar statistical scaling by means of layer-dependent t-scoring. It is argued that this normalizes the BOLD signal change by the vein-dependent signal variance (similar to assumptions underlying RSFA normalization). However, these methods are not exclusively specific to CVR-related signal variance but might be also affected by high-frequency signal fluctuations, including variance from thermal noise and CSF motion. In this sense, the RSFA approach might be a more accurate measure of vasculature-driven signal variance.

| CONCLUSIONS
We have shown that the laminar amplitude of spontaneous brain fluctuations resembles almost perfectly the scaling of the laminar amplitude of hypercapnia-induced signal changes, which are frequently used as a biomarker of CVR. This result points to the fact that RSFA measures can be used to replace tSD hc as a scaling factor for normalizing the BOLD response. Despite its different physiological origin, the calibration parameter M also showed remarkable similarities with the RSFA profile. The shape of laminar BOLD signal changes reflects spatial variations in baseline CBV v , which explains the similarities observed for the laminar profiles of M, tSD hc , and RSFA, and the consequently strong correlations.