Fast layer-fMRI VASO with short stimuli and event-related designs at 7T

Layers and columns are the dominant processing units in the human (neo)cortex at the meso-scopic scale. While the blood oxygenation dependent (BOLD) signal has a high detection sensitivity, it is biased towards unwanted signals from large draining veins at the cortical sur-face. The additional fMRI contrast of vascular space occupancy (VASO) has the potential to augment the neuroscientific interpretability of layer-fMRI results by means of capturing com-plementary information of locally specific changes in cerebral blood volume (CBV). Specifically, VASO is not subject to unwanted sensitivity amplifications of large draining veins. Because of constrained sampling efficiency, it has been mainly applied in combination with efficient block task designs and long trial durations. However, to study cognitive processes in neuroscientific contexts, or probe vascular reactivity, short stimulation periods are often necessary. Here, we developed a VASO acquisition procedure with a short acquisition period (895 ms volume acqui-sition) and sub-millimetre resolution. During visual event-related stimulation, we show reliable responses in visual cortices within a reasonable number of trials (∼20). Furthermore, the short TR and high spatial specificity of our VASO implementation enabled us to show differences in laminar reactivity and onset times. Finally, we explore the generalizability to a different stimu-lus modality (somatosensation). With this, we showed that CBV-sensitive VASO provides the means to capture layer-specific haemodynamic responses with high spatio-temporal resolution and is able to be used with event-related paradigms.


Introduction 20
Cortical layers give rise to fundamental processing units like the cortical microcircuit (Bastos 21 et al., 2012;Douglas and Martin, 2004) and may inform about feedforward and feedback char-22 acteristics of brain regions in human and animal brains (Felleman and Van Essen, 1991;Markov 23 et al., 2014). Due to technological advancements in the last decade, neuroscientists are now 24 able to study these structures non-invasively in humans using functional magnetic resonance 25 imaging (fMRI) at high resolutions (<1µL voxel volume, Polimeni et al., 2010). However, when 26 approaching this spatial scale, several challenges remain. Most notably, the so-called "draining 27 vein" effect is known to have a major impact on the fMRI signals measured across cortical 28 depths when using the popular gradient echo (GE) blood-oxygenation level dependent (BOLD) 29 sequences (Turner, 2002). Specifically, the venous blood is drained from deep to superficial 30 layers within the cortex leading to a spatial displacement of neuronal responses measured indi-31 rectly by the BOLD signal. This decreases the effective spatial specificity, despite having small 32 voxel sizes (De Martino et al., 2013;Menon and Goodyear, 1999;Self et al., 2019). 33 To counter the neurovascular confounds in the BOLD signal, additional measurements of 34 cerebral blood flow (CBF) and cerebral blood volume (CBV) have been proposed. These 35 approaches promise higher specificity to the neural site of activation by being less affected by the 36 draining veins (L. Huber et al., 2019). One of the most widely used non-invasive measurements 37 of CBV is vascular space occupancy (VASO, Hua et al., 2013;Lu et al., 2003). Slice-selective 38 slab-inversion (SS-SI) VASO (L. Huber et al., 2015, henceforth referred to as VASO) has been 39 used to study fundamental neuroscientific mechanisms on the mesoscopic scale (L. Huber et 40 al., 2017;Persichetti et al., 2020;Y. Yu et al., 2019). However, the challenging contrast-to-41 noise-ratio of sub-millimeter fMRI and the constrained sampling efficiency of VASO have so 42 far hindered wide adaptation of its application in neuroscientific application studies with fast 43 stimulus designs. While previous studies suggest that slow event-related paradigms with 6 s 44 long task durations are feasible with layer-fMRI VASO (Persichetti et al., 2020), it is not clear 45 if VASO can be pushed to even faster designs. Furthermore, it has not been clear until now if 46 and how much the detection sensitivity is compromised in layer-fMRI VASO for event-related 47 design tasks compared to commonly used block-design tasks. 48 Event-related task designs are widely-adopted in fMRI to conduct cognitive neuroscience 49 research (Huettel, 2012). This is because event-related task designs provide flexibility in stim-50 ulus presentation to control for a wide range of confounds like psychological effects (carryover 51 effects, habituation, anticipation, etc., Rosen et al., 1998). Other uses of event-related task 52 designs revolve around studying the individual trials e.g. to characterize learning, correlating 53 responses to behavioral variables (e.g. reaction time), or post-hoc sorting of trials based on 54 e.g. subjective perceptions [e.g. Formisano et al., 2002;Heynckes et al., 2023;Schneider et al., 55 2019). Although some of these aspects are also achievable with longer stimulus durations (e.g. 6 56 second long randomized stimuli as implemented by Persichetti et al. (2020), most other benefits 57 of event-related task designs can only be exploited by using shorter stimulus durations. 58 In this study, we implemented, tested, and validated a VASO sequence protocol with a short 59 repetition time (TR, 2.57 s) that provides the means to acquire layer-specific haemodynamic 60 responses with high temporal resolution and sufficient coverage to capture multiple distant 61 brain areas. For this, we used high resolution (0.88 mm3 isotropic nominal voxel size) fMRI at 62 ultra-high field strength (7 Tesla). We further characterized the sub-millimeter CBV responses 63 and their detection sensitivity to block and event-wise stimulation. Third, we explored the 64 generalizability to different stimulus modalities (vision and somatosensation). As the VASO 65 sequence yields VASO and BOLD images, we made use of both contrasts to present our results, 66 which offers complementary information. Using our new sequence protocol, we showed that it is 67 possible to capture CBV-VASO responses to stimuli with short durations in a fast event-related 68 design within a reasonable amount of trials. Our results form the foundation of a new avenue 69 for further adoption of VASO in cognitive neuroscience and physiology studies.   Nakamura et al., 2016), reconstruction with 8 iterations, FLASH GRAPPA 3 (Talagala et al.,84 2016), varying flip angles between 26 and 90°, read bandwidth = 1266 Hz/Px, first phase 85 encoding direction along anterior-posterior axis with axial-like tilted slice orientation, matrix 86 size = 152 × 148, field of view (FOV) = 133 × 129.542 mm in read and first phase en-87 coding directions respectively. We used the vendor-provided 'improved GRAPPA WS' algo-88 rithm with at least 1000 fold over regularization and small GRAPPA kernels of 2×3 (phase 89 × read). Parameters are summarized in Figure 1B, a scan protocol is available on github 90 (https://github.com/sdres/sequences/blob/master/20211022 Seb.pdf) and the sequence is avail-91 able for VB17B via SIEMENS' C2P app store Teamplay.

92
Two participants were invited for a second session, in which we compared the above protocol 93 to a version of the sequence with a longer TR. Specifically, we inserted 500 ms delays between 94 the image readouts with and without blood-nulling, as well as between the readout of the BOLD 95 (also referred to as 'not-nulled') image and the consecutive inversion pulse. Thus, the pair TR 96 was 3570 ms, while all other parameters remained equal.

97
Slice position and orientation were chosen individually for each participant. We covered 98 bilateral calcarine sulci and the hand area of the right postcentral gyrus (both indicated in 99 Figure 1A). We localized the area of interest in the right somatosensory cortex based on 100 the position of the 'hand-knob' area on the precentral gyrus, opposite of which the finger 101 representations are located on the postcentral sulcus. Note that, targeting both the visual and 102 the somatosensory cortex turned out to be challenging without significant fold-over artifacts 103 due to the small FOV and limited number of slices. In cases where coverage of both visual and 104 somatosensory cortices was not possible (n = 5), we prioritized the visual cortex.

105
Finally, we either acquired high-resolution whole-brain anatomical images (0.7 mm isotropic, 106 240 slices, TI1/TI2/TR/TE = 900/2750/5000/2.47 ms, FA1/FA2 = 5°/3°, bandwidth = 250 107 Hz/Px, acceleration factor = 3, FoV = 224 × 224 mm) using MP2RAGE (Marques et al., 2010), 108 or images were available from previous scans with similar parameters. For one participant, we 109 neither acquired MP2RAGE data due to time constraints, nor was there previous data available. 110 Stimuli 111 Stimuli consisted of flickering checkerboards presented at 16 Hz and vibrotactile stimulation of 112 all 5 digit-tips (left hand) by means of a piezoelectric stimulator at 25 Hz (mini PTS system, 113 Dancer Design, UK). Both means of stimulation are displayed in Figure 1D. Stimuli were pre-114 sented either in blocks or as events while the run duration was kept to 12.5 minutes in both cases 115 ( Figure 1C). We acquired 1-2 blocked-and 2-4 event wise stimulation runs per participant for 116 the main experiment. For the comparison between long and short TR protocols, we acquired 4 117 runs of block-wise stimulation runs (2 per TR duration). During block-wise stimulation runs, 118 we presented stimuli for 30 s with 30 s of resting fixation in between blocks, resulting in 12 119 blocks per run with an additional rest period before the first stimulation block. Event-wise 120 stimulation runs started with a 19 s baseline period and then stimuli were presented for 2 s per 121 trial with inter-trial intervals between 3 and 10 s chosen from a uniform distribution, resulting 122 in 84 trials per run (12.5 minutes). These highly irregular inter-trial intervals were chosen to 123 accurately capture the haemodynamic responses of BOLD and VASO. The stimulation pattern 124 was generated once and then used for all runs in order to allow averaging. The stimulation was 125 controlled via Psychopy 3 v2020.2.4 (J. Peirce et al., 2019;J. W. Peirce, 2007) scripts, which can 126 be found at: <https://github.com/sdres/eventRelatedVASO/tree/main/code/stimulation>. 127 Figure 1: fMRI data acquisition and stimulation protocol. A Functional EPI data (warm colors) overlaid on MP2RAGE UNI image (sub-07). The fMRI slab covered bilateral calcarine sulci (lower right) and when possible the hand area of the right postcentral gyrus (upper right). B Functional EPI data acquisition details giving a schematic overview of relevant imaging parameters. (Note: nulled and BOLD images are acquired in an interleaved fashion. Both volumes are later used to derive the 'VASO' images.) C Illustration of block and eventwise stimulation runs. Block-wise stimulation consisted of 30s on-off periods. Event-wise stimuli were 2 seconds in duration separated by inter-trial intervals of 3-10 seconds which were randomly chosen from a uniform distribution. D Means of stimulation. Left: Flickering checkerboard. Contrast reversals were displayed at 16 Hz. Right: A single piezoelectric stimulation device. The silver disc vibrates at 25 Hz and stimulates the fingertips. On the right, 5 devices are attached to all fingers of the left hand with elastic tape.
In the scenario described above, visual and tactile stimulation was always presented con-128 currently. For five participants, we deviated from this stimulation pattern. During one session 129 of one participant, the vibrotactile stimulation device was not available. Therefore, we only 130 presented visual stimuli during that session. For four other participants, we randomized vi-131 sual only and visuo-tactile stimulation during the event-related runs. This way, we aimed to 132 disentangle potential multisensory and physiological effects of laminar reactivity. The analysis code is available at <https://github.com/sdres/eventRelatedVASO>. Briefly, 135 initial motion-correction was performed within runs for nulled and not-nulled time-series sepa-136 rately using ANTsPY v0.2.7 (Avants et al., 2011). We used an automatically generated brain 137 mask (3dAutomask in AFNI 22.1.13, Cox, 1996) to calculate cost functions in brain tissue only. 138 We then computed T1w images in EPI-space derived from the combined nulled and BOLD im-139 ages for each run (3dTstat with -cvarinv) using AFNI, performed bias field correction using 140 ANTS v2.3.5.dev238-g1759e (N4BiasFieldCorrection, Tustison et al., 2010) and registered the 141 resulting image to the first event-wise stimulation run using ANTsPY. We then performed the 142 motion correction again from scratch, while concatenating the transformation matrices from 143 within and between run registration. Due to the small imaging slab, registering different ses-144 sions of the same participant was challenging. Therefore, we treated them as independent 145 datasets. However, as they were only acquired to assess differences between long and short TR 146 flavors of the protocol, this does not influence the results of our main experiment.

147
The VASO sequence acquires images with BOLD and CBV weighting concomitantly in 148 an interleaved fashion. Therefore, we temporally upsampled the motion corrected data with 149 a factor of 2 (3dUpsample in AFNI with 7th order polynomial interpolation). We then per-150 formed BOLD correction and computed standard measures for quality control (tSNR,skew,151 and kurtosis) using LayNii (v2.2.1, LN BOCO and LN SKEW respectively; L. ( Huber et al.,152 2021). General linear model (GLM) analyses were performed in FSL (Woolrich et al., 2001, 153 v6.0.5.2, contrast: stimulation vs rest). MP2RAGE UNI images were registered to the T1w 154 EPI image of the first event-wise stimulation run, first, using manual alignment in ITK-SNAP 155 (v3.8.0, Yushkevich et al., 2006, then an automatic alignment (affine transformation, mutual 156 information cost metric) while using the brain mask generated for the motion correction to 157 guide registration. Finally, anatomical images were registered nonlinearly to the T1 weighted 158 image in EPI space using ANTS. contrast. For the participant, where no MP2RAGE data was available, we drew the borders 163 based on residual contrast of the T1-weighted images directly in EPI space. All manually 164 drawn masks are available alongside the raw data. Finally, layering was performed on spatially 165 upsampled data (3dresample factor 5 in-plane) using the equidistant approach as implemented 166 in LayNii (LN2 LAYERS). For the estimation of cortical depth profiles of GLM analyses, we 167 defined 11 depth bins, while for the finite impulse response (FIR) analysis across cortical depth, 168 we estimated three depth bins (deep, middle, superficial) to allow for appropriate averaging. 169 For the analysis of the block-wise stimulation runs, we ran a GLM with one predictor for 170 the stimulation times, convolved with a standard gamma haemodynamic response function 171 without temporal derivative (mean lag: 6 s, std. dev: 3 s). Here, we applied a high-pass filter 172 (cutoff = 0.01 Hz) and no additional smoothing. To estimate the temporal response for blocked 173 stimulation, we averaged epochs with stimulation with 4 volumes before stimulus onset and 8 174 volumes after cessation. The percent signal change was computed with respect to the 30 s 175 rest-periods in between stimulation blocks.

176
Responses to event-wise stimulation runs were first estimated with a GLM, similar to the 177 blocked stimulation for better comparison. We then ran a GLM analysis using finite impulse 178 response models on the event-wise stimulation data. Here, we modeled 10 timepoints after 179 stimulus onset, resulting in a window of 13.08 seconds.

180
To estimate the number of trials necessary to obtain a reliable signal, we extracted time-181 courses from the ROIs, epoched the data with a window from stimulus onset until 8 s after 182 cessation. We then computed the overall epoch activity across all participants. For each run, 183 we then computed a cumulative average for 1 up to 84 trials, where 84 was the number of trials 184 in a run. For each number of averaged trials, we computed the sum of squares between the 185 values for each timepoint separately and the average response across all participants. Based 186 on anatomical landmarks, we will refer to ROIs in the visual cortex as V1 and ROIs on the 187 postcentral gyrus as S1.

189
In the following, we will focus on the results in visual cortices as this was the primary area of 190 interest of this study. The exploratory results for the somatosensory cortex will be described 191 briefly afterwards. Finally, one participant (sub-10) consistently showed excessive head motion 192 well above the voxel size and was therefore excluded from further analyses (Supplementary 193 Figure S1).

194
Regarding the main experiment, we first examined the results from block-wise stimulation for 196 two reasons. First, we compared our results to previous research to see whether the results 197 obtained with our fast protocol match the expected patterns. Second, we used these data to 198 determine ROIs for further analyses.

199
As expected, block-wise stimulation elicited strong signal changes in the visual cortices for 200 both BOLD and VASO in all participants. Figure 2A shows GLM activation in response to 201 stimulation (z-scores, contrast: stimulation vs. rest) overlayed on an MP2RAGE UNI image 202 for one individual. As expected, BOLD activation scores are generally higher than for VASO. 203 On the other hand, the active voxels are mostly constrained to cortical GM for VASO, while 204 BOLD shows largest z-scores in CSF. Based on this, we manually delineated ROIs from which 205 we extracted signal changes over time and z-scores across cortical depth. As expected, the 206 temporal profiles of the BOLD and VASO responses show sustained activity, peaking around 207 10s with an additional off-response after stimulus cessation ( Figure 2B). Interestingly, we 208 also observed a strong post-stimulus undershoot for BOLD, while this was less pronounced in 209 VASO. Finally, Figure 2C shows the activation profile across cortical depths. Here, we see the 210 well reproduced result of increasing activity towards the surface for BOLD and the inverted 211 U-shaped activity pattern for VASO. These likely reflect the drainage of deoxygenated blood 212 towards the pial surface and the thalamocortical input to layer 4, respectively (Felleman and 213 Van Essen, 1991;Turner, 2002). Taken together, these results assured us that our protocol 214 gives robust results with paradigms.

215
Event-wise stimulation 216 Event-wise stimulation also elicited signal changes in the visual cortices for both BOLD and 217 VASO in all participants, albeit to a lesser degree. Figure 3A shows GLM activation in 218 response to stimulation (z-scores, contrast: stimulation vs. rest) for two runs of one individual 219 overlayed on an MP2RAGE UNI image. Despite the lower z-scores compared to block-wise 220 stimulation in both BOLD and VASO, clear patterns can be seen without the need for additional 221 processing (e.g. denoising/ smoothing). Specifically, activated VASO voxels follow the cortical 222 ribbon as expected, whereas BOLD responses are strongest at the cortical surface and in CSF. 223 Figure 2: Block-wise stimulation results show robust activation. A GLM activation in response to stimulation (z-scores, contrast: stimulation vs. rest) overlayed on MP2RAGE UNI image. Note that here we show results from a single run (12.5 minutes). B Group-level averages showing the BOLD and VASO signal change across all layers in response to stimuli with a duration of 30s (shaded gray area). Colored shaded areas represent 95% confidence intervals across blocks. C Group-level layer profiles showing GLM activation in response to stimulation (z-scores, contrast: stimulation vs. rest) for BOLD and VASO in V1 across cortical depths. Colored shaded areas represent 95% confidence intervals across participants.
To investigate the haemodynamic response over time, we performed a deconvolution analysis 224 using a FIR model as implemented in FSL with three layer bins ( Figure 3B). Responses for 225 both BOLD and VASO followed the expected haemodynamic response. As expected, BOLD 226 signal changes increased from deep to superficial layers. This was less apparent in VASO. 227 To our surprise, group-level VASO responses showed comparable signal changes for middle 228 and superficial layers as opposed to the strongest response in middle layers for block-wise 229 stimulation. This was highly variable across participants (Supplementary Figure S2). Most 230 notably however, responses in superficial layers were delayed with respect to deep and middle 231 layers (white arrow in the right panel of Figure 3B). This finding was highly consistent across 232 participants, irrespective of superficial layer response amplitude (Supplementary Figure S2). 233 Also note that the response in deep layers does not differ strongly between BOLD and VASO. 234 This might indicate similar detection sensitivity between BOLD and VASO for GM close to 235 the WM/GM border.

236
The stimulus consisted of visual and tactile stimulation and so, we tested whether the lam-237 inar timing differences might be explained by effects of multisensory integration. For this, we 238 randomly presented visual or visuo-tactile stimuli in an event-related fashion in 4 participants. 239 The results are shown in Supplementary Figure S3. Briefly, visual and visuo-tactile stimu-240 lation evoked very similar responses in superficial and middle layers for BOLD and VASO, thus 241 ruling out multisensory effects as an explanation for the delayed peak in superficial layers for 242 VASO. In deep layers, visuo-tactile stimulation showed a slightly prolonged response compared 243 to visual stimulation only. This effect was present in the BOLD and VASO data, while being 244 more pronounced in the latter. However, this effect is rather small (within error bars) and has 245 to be interpreted with caution.

246
Finally Figure 3C shows the activation profile across cortical depths. Also here, the 247 increasing activity towards the surface for BOLD is clearly visible. The inverted U-shaped 248 activity pattern for VASO is also present albeit less clear than for the block-wise stimulation 249 results. Taken together, these results show that it is possible to capture event-related responses 250 to short stimuli using CBV-based VASO measurements. 251 Figure 3: Event-wise stimulation results also show robust activation. A GLM activation in response to event-wise stimulation (z-scores, contrast: stimulation vs. rest) in V1 of a representative participant (sub-14) overlayed on an MP2RAGE UNI image. Note that here we show data from two runs (25 minutes). B Group-level average responses showing the BOLD (left) and VASO (right) model fit (FIR model with 10 predictors) for three layers (deep, middle, superficial) in response to stimuli with a duration of 2 seconds (shaded gray area). The white arrow indicates the delayed peak for the superficial cortical depth in VASO. The colored shaded areas represent 95% confidence intervals across runs. Supplementary Figure S2 shows the same plot for 3 individuals. C Group-level layer profiles showing GLM activation in response to stimulation of 2 seconds (similar to Figure 3A, z-scores, contrast: stimulation vs. rest) for BOLD and VASO in V1 across cortical depths. Colored shaded areas represent 95% confidence intervals across participants. Figure 4A shows the response stabilization when cumulatively averaging trials. Here, the 253 average of a given number of trials is compared to the response across all participants. The 254 number of average trials ranges from 1 to 84 (a run had a maximum of 84 trials). Thus, the 255 last datapoint on the x-axis constitutes the difference between the average response of a run 256 and the overall mean. While the difference between a few trials and the response across all 257 participants is large, the response quickly approaches the overall average within a few more 258 trials. Across participants, the 5% of the error dynamic range for VASO is reached when 259 averaging 20 trials. The 5% of the error dynamic range for BOLD is reached when averaging 260 23 trials. The slightly higher number of averages needed for BOLD is due to the smaller 261 dynamic range. Still, this shows that the noise level of our data is acceptable to efficiently 262 capture the response to event-related stimuli without the need to average over long periods 263 of time. Note that we used strong stimuli (flickering checkerboards) and did not perform any 264 experimental manipulation. Most neuroscientific investigations use weaker tasks with control 265 conditions that are very similar to the experimental condition. Therefore, in these experiments 266 the number of presented stimuli will most likely have to be higher in order to reliably separate 267 responses to the conditions. In Figure 4B we directly compared GLM activation in response 268 to block-and event-wise stimulation for BOLD and VASO separately. In our BOLD data, we 269 see lower signal intensity for event-than for block-wise stimulation. This decrease in activation 270 scores is similar across cortical depths. In VASO, the pattern is more complex. While overall 271 signal intensity does not differ much between block-and event-wise stimulation for deep and 272 superficial layers, the group-level response peak is more biased towards the surface for event-wise 273 stimulation. As discussed before, this response peak shows large variability across participants 274 (see previous section and Supplementary Figure S2). Further, we aimed to estimate the 275 change in detection sensitivity when employing block-vs. event-wise stimulation designs. From 276 the data plotted in Supplementary Figure 4B, we approximated the area under the curve 277 for block-and event-wise stimulation separately for BOLD and VASO, by summing the z-scores 278 of all the layer-bins for block-and event-wise stimulation. We then subtracted the sum of the 279 event-wise from the block-wise results and normalized them to the latter. Across subjects, 280 detection sensitivity decreased 36 % in BOLD and 22 % in VASO. Note, we performed 1-2 281 block-wise and 2-4 event-wise stimulation runs per participant. Still, the direct comparison of 282 block-and event-wise stimulation shows that it is possible to use event-related designs within 283 an acceptable acquisition time of 25 minutes (see Figure 4A). Figure 4: VASO responses stabilize after 20 averaged trials and event-wise stimulation yields lower but plausible activation compared to block-wise stimulation. A Response stabilization curves for BOLD and VASO. Note, a run had a maximum of 84 trials. Thus, the last datapoint on the x-axis constitutes the difference between the average response of a run and the overall mean. Dashed vertical line shows the 5% of the error dynamic range for VASO. B Data from Figure 2C and Figure 3C are replotted to compare block-with eventwise stimulation directly for BOLD (left) and VASO (right) separately. Note, we performed 1-2 block-wise and 2-4 event-wise stimulation runs per participant.

284
Exploratory results in the somatosensory cortex 285 In 5 sessions we succeeded in including both the visual and right somatosensory cortex in the 286 fMRI imaging slab. Therefore, we can explore the applicability of event-related stimulation 287 using VASO in S1. Here, we obtained similar results as in V1, albeit less clear.

288
In general, fingertips are represented sparsely along the cortical sheet. Therefore, even if 289 we managed to include S1 in some participants, we did not necessarily include all five finger 290 representations for each of them. Therefore, we limited our analysis to the likely site of one 291 fingertip. Figure 5A shows GLM activation in S1 in response to block-wise stimulation (z-292 scores, contrast: stimulation vs. rest) overlayed on an MP2RAGE UNI image. For the BOLD 293 data, we can see clear activation. For VASO, only a few voxels exceed the detection threshold. 294 Nevertheless, we were able to observe sustained periods of elevated signal in both imaging 295 modalities Figure 5B) and group-level layer profiles showing GLM activation in response 296 to block-wise stimulation ( Figure 5C). Also the event-wise stimulation results show similar 297 qualities of the responses in V1, both with higher noise levels ( Figure 5D&E). Namely, just 298 like in V1, in S1 VASO responses in superficial layers show indications of a later peak compared 299 to deeper layers. These results point towards the applicability of event-wise stimulation using 300 VASO in S1. 301 Figure 5: Exploratory results point towards feasibility of event-related designs in S1 using VASO. A GLM activation in response to block-wise stimulation (z-scores, contrast: stimulation vs. rest) overlayed on MP2RAGE UNI image. Note that here we show data from two runs (25 minutes). B Group-level (n = 5) averages showing the BOLD and VASO signal change across all layers in response to stimuli with a duration of 30s (shaded gray area). Colored shaded areas represent 95% confidence intervals across blocks. C Group-level (n = 5) layer profiles showing GLM activation in response to block-wise stimulation (z-scores, contrast: stimulation vs. rest) for BOLD and VASO in S1 across cortical depths. Colored shaded areas show 95% confidence intervals across participants. D Group-level (n = 5) average responses showing the BOLD (left) and VASO (right) model fit (FIR model with 10 predictors) for three layers (deep, middle, superficial) in response to stimuli with a duration of 2 seconds (shaded gray area). The colored shaded areas represent 95% confidence intervals across runs. E Group-level (n = 5) layer profiles showing GLM activation in response to stimulation of 2 seconds (similar to Figure 3C, z-scores, contrast: stimulation vs. rest) for BOLD and VASO in S1 across cortical depths. Colored shaded areas represent 95% confidence intervals across participants.
Because we stepped into uncharted territory with our short acquisition times, we wanted to test 303 whether the activation profiles across cortical depth are influenced by the short TR compared to 304 conventional acquisitions with longer intervals between readouts. SS-SI VASO is based on the 305 assumption that the intravascular water magnetization in the imaging slab has been inverted 306 only once (L. Huber et al., 2014). This means that blood should be refilled in the period between 307 two consecutive inversion pulses (here 2.57s). While this condition is expected to be fulfilled 308 for the thin imaging slab used here, this has not been validated for the unprecedented fast 309 imaging of this study. Therefore, we directly compared data from our new protocol and a more 310 conservative approach. In the latter, we added 500ms delays between the nulled and BOLD 311 acquisitions and between the BOLD and the consecutive inversion pulse. As this decreased 312 the efficiency of our temporal sampling, we tested this approach using block-wise stimulation 313 (30s on-off) with our visuo-tactile paradigm. We found several differences between the two 314 types of acquisitions. First, we observed a reduction of T1w contrast for short compared to 315 long TR acquisitions ( Figure 6A). While we can clearly see the borders between GM and 316 CSF/WM, for long TRs, the borders between GM and WM are less visible for our short TR 317 protocol. This different structural T1-contrast across TRs is expected due to the different z-318 magnetization steady-state in absence of a relaxation delay before each inversion pulse. Because 319 of the decreased contrast in short TR protocols, we used the whole brain anatomical MP2RAGE 320 images and registered them to the functional data in order to delineate GM/WM and GM/CSF 321 borders as commonly done in GE-BOLD type acquisitions. Second, the tSNR of the BOLD 322 time-course was slightly higher for the long TR, while the VASO tSNR was largely unaffected 323 ( Figure 6B). The higher tSNR for BOLD for long TRs is likely due to the further signal 324 relaxation of the longitudinal tissue magnetization subsequent to the VASO readout. Finally, 325 the statistical GLM activation in response to stimulation (z-scores, contrast: stimulation vs. 326 rest) was slightly higher for long-than short TRs ( Figure 6C). However, this difference was 327 rather uniform across cortical depths. The slightly lower z-scores are to be expected due to 328 the lower tSNR in BOLD. Evidently, this hampered the fit of the GLM for the BOLD data. 329 The differences we see for the VASO data can be attributed to the dynamic division of BOLD 330 and nulled data during the BOLD-correction step of the VASO data. Here, a higher noise level 331 of the divisor will also enhance noise in the resulting signal. In conclusion, as can be seen 332 in Figure 6C, the layer-profiles for long and short TRs were qualitatively not substantially 333 different thus ensuring comparable signal origin and applicability of the VASO assumptions 334 (L. Huber et al., 2014). 335 Figure 6: Comparison between short and long TR fMRI data in two participants. A T1w image in EPI-space for long (upper) and short (lower) TR fMRI data. Note the high tissue contrast between white and gray matter in the long TR data. B tSNR comparison between short and long TRs. Right: tSNR maps for long and short TR data. Left: Voxel-wise tSNR kernel density estimation plots (similar to smoothed histograms) for both participants together. Values were taken from V1 gray matter (example indicated in the uppermost map on the right). C Comparison between short and long TR data activation (z-scores, contrast: stimulation vs. rest) across cortical depths for both participants together. Values were taken from the same V1 gray matter regions as tSNR.
In this work, we characterized the applicability of VASO for event-related paradigms. With 338 a short 895 ms volume acquisition TR, we were able to capture the haemodynamic response 339 for VASO and BOLD within a few (∼20) average trials. Furthermore, the short TR and 340 high specificity of VASO enabled us to show subtle laminar-dependent timing differences of 341 neural processing and vascular reactivity features. Specifically, superficial layers showed delayed 342 responses in VASO but not in BOLD. Besides general advantages of event-related designs, we 343 believe that this result demonstrates the added value of employing short stimulus paradigms 344 using VASO, which greatly increases the range of possible research questions to be addressed 345 with CBV-based laminar fMRI.

346
Response peak location and timing across cortical depth 347 A distinctive feature of VASO is its microvascular weighting, which leads to strong activation 348 of middle cortical layers during bottom-up processing (Akbari et al., 2022;Jin and Kim, 2008;349 Y. Yu et al., 2019). In our data, microvascular weighting of VASO seems to vary depending on 350 stimulus duration. For block-wise stimulation we observe a peak in middle layers ( Figure 2C). 351 For event-wise stimulation we observe a peak between middle and superficial layers (  ours, the authors find highest CBV amplitudes in superficial layers. In contrast to Hirano et al. 363 (2011), Shen et al. (2008) report strongest CBV responses to short stimuli in middle layers. 364 Note however, that for Shen et al. (2008) the middle layer response is only significantly different 365 from deep but not superficial layers. Therefore, the results from Shen et al. (2008) also match 366 our data, as we observe a peak between middle to superficial layers for event-wise stimulation. 367 In these preclinical studies, the relative microvascular weighting of CBV responses for short 368 compared to long stimuli has been discussed in the context of laminar differences in the vascular 369 architecture. Specifically, microvascular blood compartments tend to have the earliest onset 370 times (<1 s, Silva and Koretsky, 2002;Tian et al., 2010;X. Yu et al., 2014), but take long to 371 peak (>15 s) (Mandeville, 2012). On the other hand, larger, actively muscle-controlled arterial 372 compartments closer to the cortical surface, tend to have later onset times (>1 s, Tian et al.,373 2010) but shorter times to peaks (Kennerley et al., 2012). For similar investigations of response 374 times across CBV compartments in humans, time-resolved CBV measurements are necessary, 375 but were not available so far.

376
Our event-related results indicate, for the first time, similar reactivity of CBV-compartments 377 in humans. Specifically, we report earliest responses in deep and middle layers, paralleling the 378 early responses in microvascular compartments found in rodent models (X. Yu et al., 2014). 379 Furthermore, we find a delayed response-peak in superficial layers for VASO, supporting later 380 macrovascular responses (Kennerley et al., 2012;Tian et al., 2010). To our knowledge, this is 381 the first demonstration of laminar response time differences in humans using non-invasive CBV 382 measurements. On the other hand, we do not find the longer time to peak for superficial layers 383 in the BOLD data, as reported by Siero and colleagues (Siero et al., 2011;Siero et al., 2013). 384 However, this can likely be explained by two factors. Firstly, Siero et al. report a smaller effect 385 size. In their data, the difference in the time to peak is ∼0.23 s/mm (Siero et al., 2011). In 386 the visual cortex, this would result in peak time differences between deep and superficial layers 387 of ∼0.6 s. This effect is much smaller than in our VASO data (peak time differences of about 388 1.3 s). Secondly, the higher temporal resolution and shorter stimulus durations in their studies 389 compared to ours. Here, we invested image encoding time in higher spatial resolutions and used 390 longer stimuli. This allowed us to differentiate layer effects with less partial voluming of the 391 temporal evolution of pial vessels. The subtle depth-dependent timing of BOLD-and CBV-392 haemodynamic responses across stimulus durations in humans could be investigated further by 393 using systematic jitters to increase effective sampling rate.

395
In VASO protocols with longer acquisition times (e.g. >1.5 s volume acquisition time), the 396 inherent T1 contrast is sufficient to delineate GM/WM and GM/CSF borders in EPI space 397 (L. Huber et al., 2015). However, in VASO protocols with shorter acquisition times (e.g. <1 s 398 volume acquisition time), such as the one used here, the T1 contrast of the EPI images is weaker 399 due to reduced relaxation times ( Figure 6A). Thus, acquisition and registration of anatomical 400 reference images is necessary. Another option for future studies might be to acquire additional 401 VASO run(s) with longer acquisition times. For example, the acquisition of run(s) with a 402 long TR could be straightforwardly implemented in the paradigm when acquiring independent 403 localizers with strong tasks using block-wise stimulation. Furthermore, we saw slightly lower 404 tSNR values in the BOLD data for short TR acquisitions ( Figure 6B). This also translated 405 into slightly decreased activation scores ( Figure 6C). We estimated the decrease in detection 406 sensitivity when employing event-related paradigms with 2 s stimulation as opposed to 30 s 407 on/off block-designs and found a 36 and 22 % decrease in detection sensitivity for BOLD and 408 VASO, respectively. Therefore, if event-related paradigms are crucial, and detection sensitivity 409 is required to remain comparable to block-wise stimulation, up to 2 times longer experiments 410 may become necessary in future event-related layer-fMRI VASO studies.

411
Limitations & Future directions 412 We see the present study as a proof of concept that event-related stimulation is feasible using 413 VASO. As a result, there are several aspects that can be improved in future, more extensive 414 investigations.

415
In the design of our stimulation protocol, we opted for efficiency. Therefore, the inter-trial 416 intervals we chose are rather short (3-10 s), which lead to highly overlapping haemodynamic 417 responses between events. While this might be appropriate for neuroscience experiments in 418 order to investigate depth-dependent differences between haemodynamic responses to different 419 stimuli (Dale and Buckner, 1997;Glover, 1999), investigations of the haemodynamic response 420 per se might be compromised. Rather, longer inter-trial intervals could be chosen to let the 421 response return to baseline (van Dijk et al., 2021). This would allow researchers to account 422 for higher-order non-linear HRF effects that are not commonly considered in the linear super-423 position assumption of GLM/ FIR analysis models. However, this would render the paradigm 424 less efficient. Another aspect that is compromised by fast acquisition times at high resolutions 425 is brain coverage. The coverage of our high temporo-spatial VASO protocol (133 × 129.542 426 × 114.08 mm) was sufficient to image V1 in all, and V1 and S1 in some participants. This is 427 common for most layer-fMRI applications, which focus on single areas of interest (Schluppeck 428 et al., 2018). However, if the goal is to study brain-wide distributed networks, this is not 429 sufficient. Therefore, sacrificing imaging speed might be necessary to investigate connectivity 430 between distant brain regions across cortical depth (Koiso et al., 2022).

431
In general, we see various routes for future applications. The most immediate benefit would 432 be for neuroscientific applications that necessitate short TRs for efficient sampling.  (Goebel et al., 2003) or dynamic causal modeling 442 (Friston, 2012) may further corroborate investigations of the directional connectivity, however, 443 these methods are highly reliant on fast sampling of responses which, to date, was not available 444 with VASO. With our advances in short TRs in VASO, we are approaching the feasibility of 445 these methods in future applications.

446
Finally, measuring different aspects of the haemodynamic response (CBF, CBV, CMRO2 447 and BOLD) is pivotal for its complete characterization (Uludag et al., 2009).  leagues have investigated the haemodynamic response to short stimuli across cortical depths in 449 humans using BOLD fMRI (Siero et al., 2015;Siero et al., 2011;Siero et al., 2013). However, the 450 invasive nature, constrained sampling efficiency, and low SNR of CBF and CBV measurements 451 have so far hindered thorough investigations of the CBF-and CBV-haemodynamic response 452 at the laminar level in humans. Here, we demonstrate laminar response time differences in 453 humans using non-invasive CBV measurements. We believe that with our implementation of 454 fast VASO acquisition, we have enabled crucial investigations of the cortical depth-dependent 455 evolution of biophysical hemodynamic (Havlicek and Uludag, 2020;Puckett et al., 2016) and 456 neurophysiological processes (Petridou and Siero, 2019).

457
Conclusion 458 High-resolution VASO has proven to be a robust method to study brain function across cortical 459 depth. For an even wider application in neuroscientific research, fast acquisition schemes and 460 stimulation protocols are crucial. Huettel (2012) commented on the impact of event-related 461 designs as follows: "No other advance-not stronger magnetic fields, not improved pulse se-462 quences, nor even sophisticated new analyses-has contributed more to the popularization of 463 fMRI than event-related approaches to experimental design" (Huettel, 2012). Indeed, layer-464 fMRI VASO is highly dependent on strong magnetic fields, improved pulse sequences and 465 sophisticated analyses. Now, with this study, we hope to have contributed to the further 466 popularization of VASO by showing that fast event-related designs are possible and provide 467 meaningful insights.

469
Data was acquired at Scannexus (Maastricht, the Netherlands). We thank Chris Wiggins for 470 providing the 3rd order shimming tools used here. We thank Federico DeMartino for discussions 471 on the analysis of event-related fMRI data. We thank Vojtěch Smekal for helpful discussions 472 on the manuscript. Special thanks goes to the remaining "Maastricht layer-seminar" members 473 Lonike Faes, Miriam Heynckes, Kenshu Koiso, Alessandra Pizzuti and Yawen Wang for count-474 less discussions on topics related to layer-fMRI. We thank Benedikt Poser for kindly sharing 475 his 3D-EPI sequence code used here and helpful advice on layer-fMRI VASO optimization and 476 reconstructions. Finally, we thank Dr. Amanda Kaas, who was supervising Sebastian Dresbach 477 in the beginning of this project and provided the piezo-electric stimulation setup.  Figure S1: Participant with excessive head motion. Framewise displacements (FDs) for each run of the only participant consistently showing motion peaks greater than a voxel size (0.88mm, indicated by white dashed line). This participant was therefore excluded from further analyses. Figure S2: FIR-model results show high inter-participant variability with respect to peak layer. Same as Figure 3B but for 3 individual participants (upper: sub-05, middle: sub-08, lower: sub-14). Note that the superficial layer VASO activity varies across participants. While sub-08 resembles the group results, sub-05 shows greatest activity in superficial layers and sub-14 shows lowest responses in superficial layers. Nonetheless, the delayed peak for superficial layers is preserved in all participants. Figure S3: Comparison between visual and visuo-tactile stimulation in V1. BOLD (upper) and VASO (lower) activation (n = 4) in response to visual and visuo-tactile stimulation separately. Data is plotted separately for the three layers (deep, middle, superficial). Visual and visuo-tactile stimulation evoke very similar responses in superficial and middle layers for BOLD and VASO. In deep layers, visuo-tactile stimulation shows a slightly prolonged response compared to visual stimulation only. This effect was present in the BOLD and VASO data, while being more pronounced in the latter. However, this effect is rather small (within error bars) and has to be interpreted with caution. Still, we believe that this might indicate a multisensory integration effect taking place in deep layers of V1 in response to additional tactile input. Future studies on neuroscientific event-related applications might investigate this further.