Characterisation of laminar and vascular spatiotemporal dynamics of CBV and BOLD signals using VASO and ME-GRE at 7T in humans

Cortical layers are a key feature of mesoscopic scale brain organisation. Recent developments in functional magnetic resonance imaging (fMRI) have enabled the non-invasive mapping of activation in distinct layer-compartments in humans (Dumoulin et al., 2018). One such technique, Vascular Space Occupancy (VASO; Lu et al., 2003), is a promising cerebral blood volume (CBV)-sensitive signal. In recent years, this method (specifically slice-selective slab-inversion [SS-SI] VASO (Huber et al., 2014), henceforth referred to as VASO) has gained popularity for investigations of cortical depth-dependent signal changes at ultra-high fields (⁠$>$7T) due to its high specificity to the microvasculature and hence the underlying neural activity (e.g., Faes et al., 2023; Haenelt et al., 2023; Huber et al., 2017; Koiso et al., 2023; Persichetti et al., 2020; Pizzuti et al., 2023; Y. Yu et al., 2019; for other VASO variants at ultra-high field, see e.g., Chai et al., 2020; Hua et al., 2013).

The high microvascular weighting of VASO was established by using long stimulus durations (Beckett et al., 2020; de Oliveira et al., 2023; Jin & Kim, 2008), which were also beneficial for mitigating the challenging signal-to-noise ratio (SNR) and sampling efficiency (Huber et al., 2019). However, for neuroscience applications, the possibility of using short stimuli is often crucial (Huettel, 2012). For example, it allows accounting for carryover effects, enables researchers to study fleeting phenomena such as surprise, or correlates factors such as task performance with brain activity.

Next to challenges in SNR and temporal efficiency, a lack of knowledge regarding laminar CBV response-timings across vascular compartments in humans has hindered cognitive neuroscience applications with short stimuli using VASO. From the animal literature, it is known that micro-vessels (located most closely the neuronal tissue within gray matter) respond earliest upon neuronal activation but take long to fully dilate (Jin & Kim, 2008; Tian et al., 2010). This early response in capillaries is followed by a strong response in actively muscle-controlled feeding arteries with a slightly later onset but earlier saturation (Kennerley et al., 2012). In contrast to the microvascular response in gray matter, these larger macrovessels are located mostly in cerebrospinal fluid (CSF) above gray matter. As a result, responses in micro- and macrovessels occur on different timescales and, furthermore, differently influence responses in cortical layers. Therefore, the micro- versus macrovascular contribution to the VASO signal is expected to differ across stimulus durations and/or timepoint within the measured response (Jin & Kim, 2008).

Recently, we have shown that it is possible to overcome the difficulties in SNR and temporal efficiency and applied VASO with short (2 seconds) stimuli and an event-related paradigm (Dresbach et al., 2023; for an earlier example of stimuli shorter than 20 seconds, see Persichetti et al., 2020). Furthermore, we, indeed, found that stimuli with a duration of 2 seconds had a lower microvascular weighting compared to long stimuli. However, the relationship between (pial) macro- versus (intracortical) microvascular weighting at intermediate, or even shorter stimulus durations is currently unknown for the human VASO signal. Furthermore, responses were attributed to micro- and macrovessels based on cortical depth measurements, which only provides a rough estimate of vascular contribution.

While inferring the micro- or macrovascular weightings of the VASO signal from cortical depth measurements might provide initial insights, detailed information on the vascular architecture would enhance our understanding of signal origins even further. In animal models, this work has been pioneered in the rat somatosensory cortex by various groups. For example, X. Yu et al. (2016) established intracortical vessel maps by using multi-echo gradient-echo (ME-GRE) images and showed that blood oxygenation level dependent (BOLD) and CBV responses originate from venules and arterioles, respectively. On the other hand, Tian et al. (2010) used two-photon microscopy to measure the dilation of individual arterioles and capillaries and showed a backpropagation of arterial vessel dilation. They further compared the microscopy findings to BOLD data acquired with the same stimulation protocol and found good agreement. Although these investigations clearly demonstrate the benefit of combining vascular information with high-resolution functional imaging, comparable investigations in humans have been hindered so far by the difficulties of acquiring adequate vascular information in combination with multimodal functional data. Recently, Gulban et al. (2022) showed that acquiring high-quality ME-GRE images is feasible at the mesoscale in living humans. Specifically, the authors provided a scanning protocol that gives high-quality T2* maps acquired at 0.35 mm isotropic resolution and demonstrated that human vascular information can be resolved non-invasively down to the level of intracortical diving vessels (for an alternative approach with time-of-flight magnetic resonance angiography, see Bollmann et al., 2022). Combining this detailed vessel information with high-resolution BOLD and CBV fMRI allows us to investigate the neurovascular coupling in humans.

Therefore, the purpose of this study is twofold. First, we characterised the canonical haemodynamic response function of the layer-dependent VASO signal across stimulus durations common in neuroscience and at high spatiotemporal resolution in humans. Second, we showed a first step in the detailed investigation of submillimetre VASO and BOLD fMRI signals across stimulus durations in vascular compartments by using cutting-edge high-resolution anatomical information from ME-GRE data acquired at 0.35 mm isotropic resolution.

Five participants without any known neurological condition (3 female, 2 male, age range: 22–35 years, mean: 26.8 years) took part in the study after providing written informed consent. The study was approved by the Ethics Committee of the Medical Faculty of Leipzig University and the principles expressed in the Declaration of Helsinki were followed.

Scanning was performed in five 90-minute sessions per participant (25 sessions in total) using a 7T MRI scanner (Magnetom Terra, Siemens Healthineers, Erlangen, Germany), equipped with a 8Tx, 32Rx RF head coil (NOVA Medical, Wilmington, DE, USA) at the Max Planck Institute for Human Cognitive and Brain Sciences in Leipzig, Germany.

Functional scans were performed in 4 sessions with a 3D echo-planar imaging (EPI) sequence, with SS-SI VASO preparation (Huber et al., 2019; Stirnberg & Stöcker, 2021; VASO version 26dc5a59) and a nominal resolution of 0.9 mm isotropic (16 slices, time of inversion [TI]1 = 1047 ms, TI2 = 2462 ms, shot repetition time [shotTR] = 42 ms, echo time [TE] = 15 ms, partial Fourier factor = 6/8, variable flip angle scheme with reference flip angle = 33°, 3-fold generalised autocalibrating partially parallel acquisitions [GRAPPA] acceleration (Talagala et al., 2016) in z-direction with a “controlled aliasing in parallel imaging results in higher acceleration” [CAIPIRINHA] shift of 1, bandwidth = 1126 Hz/Px, field of view [FoV] = 133 mm, matrix size = 148 × 148; (Fig. 1A, D). B1 shimming was performed in Trueform mode (CP-mode) and the shimming volume included the circle of Willis, to ensure proper nulling of incoming blood. Slice position and orientation were chosen individually and aligned with the left calcarine sulcus (Fig. 1D, for the functional coverage in all participants, see Supplementary Fig. S1). To give participants enough rest between the long experimental sessions, we did not schedule multiple sessions in 1 day, but rather spaced the sessions across multiple weeks, if possible. Across all participants, the functional sessions were scheduled on average 8 days apart (range: 1–38 days; std. deviation: 9.48 days). To ensure similar slice positions between sessions, we made use of the head scout provided by the vendor.

Importantly, we implemented a “symmetric” VASO acquisition, in which an additional delay of 785 ms was inserted between the blood-nulled (further referred to as nulled) and BOLD EPI readouts (Fig. 1A). This delay was as long as the inversion delay and the volume TR, thus creating 4 equally spaced blocks (inversion delay, nulled EPI, additional delay, BOLD EPI) of 785 ms, each, within one acquisition cycle. Having 4 equally spaced blocks within the acquisition period had multiple benefits. First, it allowed us to systematically jitter the onset of our stimuli, thereby sampling the haemodynamic response to our stimuli with an effective sampling rate of 785 ms, uncoupling pair TR and effective temporal resolution. Second, including the delay before the BOLD EPI acquisition increased the temporal SNR (tSNR) of the BOLD data, as shown in previous experiments with a similar setup (Dresbach et al., 2023). Third, the longer T1 relaxation period between inversion pulses led to a better structural T1 contrast in the mean VASO EPI image despite a short volume TR (Dresbach et al., 2023). This was especially important in our case, because it would facilitate the registration of images across sessions.

In a fifth session, we acquired 4 3D ME-GRE images at 0.35 mm isotropic resolution (104 slices; TEs = [3.8, 7.6, 11.4, 15.2, 19, 22.8 ms]; TR = 30 ms, no partial Fourier, flip angle = 12°, GRAPPA factor for in-plane phase encoding = 2; bandwidth = 310 Hz/Px, FoV = 173 mm, matrix size = 496 × 496; volume acquisition time = 13 minutes; Fig. 6A). Acquiring 4 images at each echo had 2 purposes. First, we registered and averaged the images to increase SNR. Second, we rotated the phase-encoding orientation by 90° between acquisitions for blood motion artifact mitigation (Gulban et al., 2022; Larson et al., 1990). The slabs of these high-resolution images were centered on the functional coverage (see Figs. 1A and 6A).

Finally, whole-brain images (0.7 mm isotropic, 240 slices, TI1 = 900 ms, TI2 = 2750 ms, TR = 5000 ms, TE = 2.47 ms, flip angle 1 = 5°, flip angle 2 = 3°, bandwidth = 250 Hz/Px, GRAPPA acceleration factor = 3, FoV = 224 × 224 mm) were acquired in one or two of the sessions using MP2RAGE (Marques et al., 2010). For a full overview of the scanning sessions of all participants, see Supplementary Table S1. Scanning protocols are available on github under $<$https://github.com/sdres/neurovascularCouplingVASO/blob/main/protocols$>$⁠.

In a total of 21–25 runs spread over 4 sessions per participant, flickering checkerboards (contrast reversals at 8 Hz; Fig. 1B) were presented for 1, 2, 4, 12, and 24 seconds. Stimulus onset was systematically jittered 4 fold with respect to the TR, resulting in an effective stimulus sampling of 0.785 seconds (Fig. 1B), assuming stationarity of neurovascular responses. Each stimulus duration was presented once per jitter and run, resulting in at least 21 repetitions per participant and jitter. Based on a previous study (Dresbach et al., 2023), $∼$20 repetitions per stimulus constitute a good trade-off between efficiency and data stability. Inter-trial intervals (ITIs) were either 10, 12, 14, 20, and 24 seconds (short ITIs) or 20, 22, 24, 30, and 40 seconds (long ITIs, Fig. 1C). For an overview of which sessions contained runs with short or long ITIs, see Supplementary Table S1. For sessions with short ITIs, we generated a stimulation design per session and used it for all runs of the session. For the long ITI, we created one stimulation pattern per participant and used it for all runs with long ITIs. We argue that in the latter case, rest durations are sufficiently long to avoid carryover effects. On the other hand, using the same stimulation across more runs allows the averaging of nulled and not-nulled runs across sessions. This might help to reduce noise before dividing nulled by not-nulled time courses during the BOLD-correction. Initially, we had included short ITI sessions for higher efficiency but then opted for long ITIs due to higher signal integrity.

To keep participants’ attention, we implemented a cued button press task. To this end, a gray fixation dot (0.075 degree radius) surrounded by a yellow ring (0.15 degree radius) was presented in the center of the screen. In semi-random intervals, the yellow ring turned into a red target for 500 ms. Participants were instructed to press a button on a button box with their right index finger in a 2 second window after the target. Intervals between targets were chosen randomly between 40 and 80 seconds. The stimulation was controlled with scripts using Psychopy3 (v2022.1.3, Peirce, 2007), which can be found at $<$https://github.com/sdres/neurovascularCouplingVASO/blob/main/code/stimulation/stimulationParadigm.py$>$⁠.

Bias field correction was performed using N4BiasFieldCorrection as implemented in ANTs (Tustison et al., 2010; v2.3.5). Brain extraction was performed based on the second inversion image of the MP2RAGE sequence, using bet (v2.1) as implemented in the FMRIB Software Library (FSL; v6.0.5). The resulting brain mask was applied to the T1w UNI image, and the image was further cropped to the region of interest (ROI) to reduce data processing demands.

An initial tissue segmentation was performed on the cropped MP2RAGE UNI using FSL FAST (Zhang et al., 2001; v6.0.5). The initial segmentation was manually corrected in the ROI using ITK-SNAP (Yushkevich et al., 2006; v3.8.0) after upsampling to 0.175 mm isotropic using cubic interpolation. Layerification was performed at the upsampled resolution using the LN2_LAYERS program implemented in LayNii (Huber, Poser, Bandettini, et al., 2021; v2.2.1) following the equivolume principle. For the analysis of temporal data, we estimated 3 layers. To plot layer profiles, we estimated 11 layers.

All ROIs were defined manually using the macro-anatomical landmark of the calcarine sulcus, indicating primary visual cortex (V1). Specifically, we used the high-quality segmentation and the LN2_MULTILATERATE program in LayNii to grow discs with a cortical parametrisation (U,V coordinates), centred on the calcarine sulcus, for each participant. To validate that our ROIs were indeed located in V1, we made use of the fact that V1 is easily identifiable on inflated cortical surfaces and compared the definition to an atlas (Glasser et al., 2016). To define V1 for each participant individually, we followed a standard procedure in BrainVoyager (Goebel, 2012; Version 22.4.4.5188 for macOS) to obtain cortical surfaces. The pipeline is described in detail in the Supplementary Materials. Finally, V1 outlines were drawn corresponding to the definition in the Glasser atlas (Glasser et al., 2016) by author SD and quality controlled by author OFG. Manually drawn outlines of individual participants’ V1 are shown in Supplementary Figure S2. To check whether the ROIs fell within the definition of V1, we transformed the resulting surface representation of V1 to voxel space and visually checked the overlap of the V1 definition with our ROIs.

The VASO sequence provides two time series. The first consists of “nulled” images, in which the blood signal is as much minimised as possible—ideally 0, giving them a CBV-weighting. The second consists of T2*-weighted BOLD images. If not indicated otherwise, the two time series were treated individually. Data from individual sessions were motion-corrected in a multi-step procedure using ANTsPy (Avants et al., 2011; v0.3.3). First, all volumes of a run were aligned to the first volume of the run using 6 degrees of freedom and rigid body alignment. To guide the registration, we used a brain mask generated automatically using 3dAutomask in AFNI (Cox, 1996; v22.1.13). We then computed run-wise T1w images from the motion-corrected functional data by dividing 1 by the coefficient of variation of the concatenated nulled and BOLD time series on a voxel-by-voxel level. The resulting T1w image of each run was registered to that of the first run of the first session. Finally, the motion correction was performed again by concatenating the transformation matrices of each volume within a run and the transformation matrix generated by the between-run registration. This process was done to minimise consecutive interpolation steps which might introduce unnecessary blurring of the data. Then, the functional images were cropped to the caudal half to reduce computing and storage demands and we averaged all runs within each session to create a session-averaged run. After averaging, we temporally upsampled the data by a factor of 4 to match the stimulus onset grid, duplicated the first volume of the BOLD time series to match the nulled and BOLD data, and divided the nulled by the BOLD data to correct for BOLD-contamination. This resulted in the BOLD-corrected VASO data. Furthermore, we calculated tSNR maps using the LN_SKEW program in LayNii and, finally, we computed a T1w image as before, but now based on the averaged session data to increase SNR.

To estimate spatial stimulus responses, we ran a general linear model (GLM) analysis for the average run of each session using FSL (Woolrich et al., 2001). Here, we used an individual predictor for each stimulation duration, convolved with a standard gamma haemodynamic response function without temporal derivative (mean lag: 6 seconds, std. dev: 3 seconds), applied a high-pass filter (cutoff = 0.01 Hz) and no additional smoothing. A second-level analysis across sessions was run for each participant. No averaging of statistical maps was performed across participants.

To investigate spatio-temporal responses, we computed voxel-wise event-related averages for each participant. This was done in a multi-step procedure on a voxel-by-voxel level.

1. Responses to individual trials were extracted from the onset of stimulation until the end of the respective rest period for each temporally upsampled session-averaged run.

2. Responses to individual trials were averaged across sessions to create an event-related average for each stimulus duration separately.

3. The baseline was computed as the average of the first and final 30 seconds of each run across all session-averaged runs.

4. % signal-change was calculated for each timepoint of the event-related averages for each stimulus duration independently by dividing the signal at each timepoint by the voxel-mean calculated in step 3 and multiplying the result by 100.

5. For VASO, an increase in CBV is shown by negative signal changes. Therefore, we inverted the VASO signal changes, for easier comparison with the BOLD data.

Due to this averaging procedure, we observed an offset between the time points of the event-related average that were extracted from all runs (until the last volume of the sessions with short ITIs) and the time points that were only extracted from the runs with long ITIs (Supplementary Fig. S3). To account for this, we subtracted the difference in signal change between the last time point of the short ITI sessions and the first time point of the long ITI sessions from the time points only present in the long ITI sessions. This results in two time points with the same value for each participant and stimulus duration, which is not visible in the group data. Finally, we set the value of the first time point of a given event-related average to 0 and normalised the remaining time points to this. This procedure was done for each participant, layer, and stimulus duration independently before averaging across participants (For a graphical representation of the procedure, see Supplementary Fig. S3).

ME-GRE data were processed as described in Gulban et al. (2022). Key points of the analysis are the blood motion artifact mitigation and the correction of head motion between acquisitions. The blood motion artifact is a displacement of signal from arterial vessels due to the read phase encoding direction (Larson et al., 1990). This issue was mitigated by acquiring data with four different phase encoding directions and compositing new images with less influence of the artifact. To correct for head motion between acquisitions, we registered the averaged echoes of individual runs to the first run and applied the resulting transformation matrix to the individual echoes. In this process, we also upsampled the images to 0.175 mm isotropic using cubic interpolation and finally averaged images from four acquisitions to boost SNR. From the resulting images, we calculated T2* and S0 values. Based on the T2* data, we performed a manual vessel segmentation in ITK-SNAP restricted to our ROI. Specifically, we selected low-intensity voxels arranged in tube-like structures (Fig. 6B).

The statistical maps from the functional data analysis and the vessel segmentation based on the ME-GRE data were registered to the upsampled anatomical MP2RAGE data. To register the functional data, the upsampled UNI images of the MP2RAGE data were first registered to the T1w image in EPI space derived from all sessions. The registration was performed using the cross-correlation metric and the non-linear SyN algorithm of antsRegistration as implemented in the ANTS toolbox (Avants et al., 2011; v2.3.5) to account for the geometric distortions in the functional data. Furthermore, the brain mask generated during the motion correction was used to guide the registration. The resulting linear and non-linear transformation matrices were then used to inversely register the statistical maps and event-related averages to the anatomical MP2RAGE data at upsampled resolution (0.175 mm isotropic).

To register the ME-GRE data, the upsampled MP2RAGE inv-2 images were registered to the upsampled S0 images resulting from the ME-GRE data. Here, we first manually matched the images using ITK-SNAP. This initial alignment was then improved by registering the images using ITK-SNAP with the mutual information metric and a mask covering the posterior brain border generated during the ME-GRE processing pipeline. Finally, we used the inverse registration matrix to register the vessel segmentation to the anatomical data.

To keep participants’ attention during the functional sessions, they performed a cued button press task. Analysis of the attention task revealed that participants consistently performed well (Supplementary Fig. S4A). Only in 3 runs, the ratio of detected/total targets was $<$0.75. Given that targets could be presented alongside the task stimulus and therefore hard to detect, we did not exclude any data based on isolated undetected targets. We further examined how many consecutive targets were missed within runs, as missing multiple targets in a row could be indicative of participants falling asleep. We found that in 2 runs 2 targets, and in two other runs 3 consecutive targets were not detected. Given this low number and an average of 11 targets per run, we did not exclude data based on attention task performance.

Participant head motion is a common problem in fMRI, especially for high-resolution data. To assess head motion, we computed framewise displacements (FDs) from the rigid motion parameters given by the motion correction (Power et al., 2012). The distribution of volumes which were displaced with respect to the previous volume by more than the voxel size (0.9 mm) is shown for nulled and BOLD timecourses separately in Supplementary Figure S4B. Only a small number of volumes (92; $<$0.08 %) showed FDs greater than our voxel size (0.9 mm). 40 of those (in 16 individual runs) were in BOLD, 52 (in 12 individual runs) were in nulled time series, and FD never exceeded 1.88 mm (2 $×$ voxel size). We therefore did not exclude any data based on motion.

Figure 2A shows the high data quality in a representative participant (sub-06) after averaging data from an individual session. Both mean nulled and BOLD EPIs look crisp with minimal artifacts. Note the outstanding anatomical contrast in the mean nulled images. This helped in the correction of motion within runs and especially the registration between runs and sessions.

Figure 2B shows the tSNR across participants for VASO and BOLD within the ROIs (for an example ROI, see Fig. 3A). tSNR values were well above 10 for VASO and above 30 for BOLD. For tSNR values of individual participants and sessions within the ROI, see Supplementary Figure S5. Individual tSNR values are mostly uniform across sessions with minor differences between participants. Figure 2C shows an example of VASO and BOLD tSNR maps of one individual participant (sub-06) after averaging one session.

Taken together, we concluded that the functional data were of sufficient quality to proceed with further analyses.

Obtaining high-quality statistical maps is often difficult in layer-fMRI. Therefore, investigating voxel-wise results from the GLM analysis constitutes a rigorous assessment of data quality. For VASO, longer stimuli (⁠$>$12 seconds) evoked clear responses in visual cortices (Fig. 2D, top) showing the high data quality. Notably, the activation is centred in gray matter, as expected based on the high specificity of CBV measurements. However, we also observed sparse occurrences of highly active regions for VASO in CSF, likely reflecting a dilation of pial vessels (red arrows in Fig. 2D). While in response to 4 seconds stimulation, some activation can be found in gray matter for VASO, only weak to no visible patterns can be observed for stimuli with a duration of 1 or 2 seconds. For BOLD, all stimulus durations lead to clear responses whose amplitudes increased with stimulus duration (Fig. 2D, bottom). As expected for BOLD, strongest responses can be observed in CSF due to the draining vein effect. Taken together, the results from the statistical maps are in line with our expectations due to the limited number of trials and low tSNR of VASO compared to BOLD (Fig. 2B, C).

To quantify the signal changes across cortical depth and time, we performed a high-quality semi-manual segmentation of cortical gray matter (Fig. 3A). We then extracted event-related averages from 3 cortical layers in V1. Specifically, we computed a geodesic cylindrical ROI centered in the calcarine sulcus using LN2_MULTILATERATE as implemented in LayNii (green and layered shapes in Fig. 3A). Note that signals were averaged within the ROI, irrespective of whether all voxels were “activated” in all stimulus conditions or not. That was done in order not to bias the results with respect to certain cortical layers, BOLD, VASO, or stimulus durations.

The responses of individual participants are shown in Supplementary Figure S6 for VASO and Supplementary Figure S7 for BOLD. Most participants showed clear stimulus evoked responses for all stimulus durations and layer compartments for both VASO and BOLD. However, we did not observe reliable VASO activation in response to stimulation of 1 and 2 seconds for one participant (sub-08, see Supplementary Fig. S6). This participant was therefore excluded from further analyses. Group-level (n = 4), time-resolved VASO and BOLD responses for 3 cortical depths and all stimulus durations are shown in Figure 3B.

Most notably, VASO has slightly stronger responses in superficial compared to middle and deep layers for short stimuli (⁠$≤$4 seconds). On the other hand, longer stimulus durations evoked slightly stronger VASO responses in middle layers, which became more pronounced with increasing stimulus duration. For a more conventional representation of depth-dependent activity, we extracted laminar profiles from VASO and BOLD z-maps resulting from a GLM analysis with a predictor for each stimulus duration contrasted against the rest periods (Fig. 4). For VASO, the results show clear peaks in middle layers for all stimulus durations (Fig. 4A, left), whereas BOLD peaks are located in superficial layers (Fig. 4B, left). To show differences in profile shape irrespective of activation amplitude, we min–max normalised the responses across cortical depth (Fig. 4A, B, right). Also here, it can be seen that short stimulus durations show a slightly stronger response in superficial layers than long stimuli for VASO.

Note, that the clearer middle layer peaks for all stimulus durations in VASO based on the layer profiles (Fig. 4A) compared to the event-related averages (Fig. 3B) stems from the fact that the former are based on z-scores whereas the latter are based on raw % signal change. As superficial layers usually exhibit higher noise levels compared to middle or deep layers, responses close to the cortical surface are attenuated when the signal is plotted in terms of z-scores. For completeness, z-scored event-related averages are shown on the group level in Supplementary Figure S8 for VASO and BOLD (see Supplementary Figs. S9 and S10 for individual participants’ z-scored VASO and BOLD event-related averages, respectively). Here, middle layer activity for VASO is clearly the strongest for all stimulus durations, resembling the layer profiles. Furthermore, the clearer middle layer peaks in VASO based on the layer profiles (Fig. 4A) compared to the event-related averages (Fig. 3B) stem from the higher partial voluming when extracting signals from 3 (as done for the event-related averages) rather than 11 layers (as done for the layer profiles). For BOLD, response-strength increases from deep to superficial layers, irrespective of stimulus duration, or form of representation, as expected based on the draining vein effect.

We also find delayed VASO peaks in superficial layers compared to middle and deep layers for short stimuli (Fig. 3B). To quantify this delay, we defined the time to peak (TTP) as the time between stimulus onset and the time point of highest signal change for each stimulus duration and layer compartment separately. The results for VASO and BOLD are shown in Figure 5 (For data from individual participants, see Supplementary Fig. S13). For VASO, there are 3 noteworthy observations. Firstly, as described based on visual assessment above, the TTP was the highest in superficial layers for short stimuli. While this effect is small, it is generally consistent across participants, with 4/4, 3/4, and 3/4 participants showing increasing TTPs from deep to superficial layers for a stimulus duration of 1, 2, and 4 seconds, respectively. Secondly, the highest TTP is located in the middle layers for long stimuli. This effect is less consistent, but rather driven by individual participants. Thirdly, we find that for stimuli up to 12 seconds, the response peaks after stimulus cessation for all layers, while the response to 24 second stimulation drops off before the stimulation ends. This finding is highly consistent across participants, with only minor deviations. For BOLD, the response to stimuli with a duration of 24 seconds also peaks before stimulus cessation and shows the longest TTP in middle layers. While the shorter TTP with respect to stimulus cessation is robust, the longer TTP of middle compared to deep and superficial layers is driven by individual participants only.

In order to investigate the vascular origin of VASO and BOLD responses, we computed signal changes in vascular structures compared to gray matter. For this, we manually delineated vessels based on the ME-GRE data Figure 6B, which can be for example identified by voxels with very short apparent T2* values (due to strong dephasing inside the blood) and an arrangement in a tube-like fashion. Figure 6C shows the group level (n = 4) responses for gray matter and vessel-dominated voxels. Data from individual participants are shown in Supplementary Figures S11 and S12 for VASO and BOLD, respectively.

There are several noteworthy observations. Firstly, the vessel-dominated VASO signal change shows delayed peak times compared to signals originating predominantly from gray matter for all stimulus durations up to and including 12 seconds (Fig. 6C). This could explain the delayed peak in superficial layers of gray matter as observed in Figure 3B. Secondly, for short stimuli, the superficial vessel response is stronger than that of gray matter. On the other hand, gray matter VASO signals dominate for stimuli of 12 and 24 seconds. This is driven by a slow increase of VASO gray matter activity from short to long stimulus durations, while the response in vessel-dominated voxels is saturated after 2 seconds of stimulation and does not rise even for 24 second stimuli. For BOLD, vessel-dominated voxels show a markedly larger response than gray matter, irrespective of the stimulus duration, with responses that keep increasing for longer stimuli.

Figure 7 illustrates the ratio between peak-signal changes from voxels dominated by macrovessels and voxels dominated by gray matter for both VASO and BOLD, separately for each stimulus duration (for data from individual participants, see Supplementary Fig. S14). For VASO, the group-level data show again that only for long, and not short stimuli, the gray matter response is larger than the response in macroscopic vessels. While this trend is clearly visible for stimuli $>$4 seconds in all participants, the mean result of steadily decreasing macrovascular weighting is mostly driven by individual participants for 1 and 2 seconds of stimulation. Nevertheless, the overall pattern of vessel/gray matter signal $>$1, indicating higher macrovascular than microvascular weighting in VASO, for stimuli with a duration of 1, 2, and 4 seconds is supported by 4/4, 3/4 and 4/4 participants, respectively, while for 24 seconds of stimulation all participants show higher micro- compared to macrovascular weighting. For BOLD, the data show that the relationship between higher BOLD activity in vessels than gray matter stays very similar, irrespective of stimulus duration.

In the present study, we characterised laminar CBV responses to visual stimuli with a range of durations that are common in neuroscientific investigations. To do so, we used high spatiotemporal resolution (effective temporal sampling: 785 ms; spatial resolution: 0.9 mm isotropic) VASO fMRI, in combination with cutting-edge anatomical ME-GRE data acquired at 0.35 mm isotropic to glean vascular information. With these data, we showed that the microvascular weighting of the VASO signal depends on stimulus duration and could attribute this effect to the ratio between activity in gray matter and superficial, vessel-dominated voxels.

Our results have several implications for neuroscientific applications of layer-specific VASO, especially with respect to the stimulus duration. Firstly, we were able to obtain reliable responses to stimuli with a duration of 1 second in 4 out of 5 participants. After Persichetti et al. (2020) with a stimulus duration of 6 seconds and Dresbach et al. (2023) with a stimulus duration of 2 seconds, to our knowledge, this is the shortest stimulation duration with which VASO responses have been seen at submillimetre resolution in humans. For neuroscientific purposes, shorter stimulus durations are often desired, as they, for example, allow to investigate transitory phenomena like surprise or enable better modulation of stimulus detectability in cases where detection at threshold is desired (Huettel, 2012). Therefore, our results may facilitate future investigations of novel research questions using laminar-specific VASO.

However, the shorter stimulation time seems to come at the cost of spatial specificity to the underlying neural populations. Specifically, we further showed that the laminar weighting of VASO differed between stimulus durations. Shorter stimulation durations (⁠$≤$4 seconds) showed highest % signal changes in superficial layers, whereas during longer stimulation (⁠$≥$12 seconds) signal changes in middle layer activity were highest (Fig. 3B). Note that in contrast to raw % signal changes, the peak z-scores were in middle layers for all stimulation durations (Fig. 4A), however with a more shallow decrease towards the cortical surface for short stimuli. This can be explained by the higher noise level in superficial layers which downscales the response at the cortical surface when interpreting the data in terms of z-scores (for z-scored event-related averages, see Supplementary Figs. S8–S10). These findings are in line with our previous VASO results comparing block- to event-wise stimulation in humans (Dresbach et al., 2023) and animal work using MION contrast agent (e.g., Jin & Kim, 2008). However, in our previous study, we used either short (2 seconds), event-related or long (30 seconds), block-wise stimulation and could not determine the stimulation duration at which the relationship between micro- and macrovascular weighting changes is optimal. Based on the present study, stimulus durations of 12–24 seconds or longer might be necessary for the VASO signal to develop the maximal weighting towards the microvasculature. However, the gap between the 4 second stimulus and the 12 second stimulus studied here is still substantial. Therefore, future studies could investigate whether other durations within this range are sufficient to develop the desired spatial specificity.

Finally, we found that for stimuli up to and including 12 seconds, the VASO TTP was after stimulus cessation (Fig. 5). However, for 24 second stimuli, the response reached its highest point before the stimulus ended. Note that the TTP was longest in superficial layers for short stimulation durations, probably due to the largest contribution of larger vessels that react slightly slower than the microvasculature close to the neuronal tissue in gray matter. While this effect is small, it is generally consistent across participants (Supplementary Fig. S13). On the other hand, group-level TTP was longest in middle layers for longer stimuli, hypothetically due to the slow, passive dilation of the capillary bed. However, we found this effect to be highly variable across participants, which makes a solid interpretation difficult.

Taken together, several factors need to be considered when choosing the stimulus duration in future VASO experiments. Based on our experiments, longer stimuli (⁠$≥$24 seconds) are advisable for higher signal amplitudes, superior detection sensitivity, and good microvascular weighting. On the other hand, to capture transient brain activity dynamics, shorter stimuli are better suited. Finally, if one is interested in obtaining the highest possible magnitude of the response while being temporally efficient, stimuli with a duration between 12 and 24 seconds might be advisable.

In line with various studies using BOLD as measured with gradient-echo (Bause et al., 2020; Kay et al., 2019; Kemper et al., 2018; Moerel et al., 2018), spin-echo (Koopmans and Yacoub, 2019), or 3D-GRASE (Moerel et al., 2018), we find that BOLD signal changes are much larger in vessel-dominated voxels than in gray matter. Furthermore, our results show that the relationship between BOLD signal changes in vessel-dominated voxels and gray matter does not change substantially across stimulus durations (Fig. 7). As a result, the specificity of BOLD responses seems to be somewhat independent of the duration of stimulation.

For VASO, we also find signal changes in vessel-dominated voxels. Previously, CBV changes in upstream vessels have been shown in animals by several papers (e.g., Berwick et al., 2008; Huber, Poser, Kaas, et al., 2021; Kennerley et al., 2012). VASO responses in pial vessels of humans have also been highlighted in the past by us and other groups in numerous conference abstracts (Handwerker et al., 2016; Huber et al., 2015; Venzi et al., 2019), user manuals (Huber, 2018a, 2018b), and in a previous paper (Huber, Poser, Kaas, et al., 2021). However, they have never been systematically investigated. In the following we will discuss three main observations.

Firstly, we find consistently delayed onsets of vessel-dominated compared to gray matter VASO signal changes for all stimulus durations. This matches our previous results in humans (Dresbach et al., 2023) and results from Tian et al. (2010) who found earliest dilation of capillaries in middle gray matter, followed by a backpropagation towards the cortical surface and later activation of larger feeding vessels in rats, using two-photon microscopy.

Secondly, we find consistently delayed peaks of vessel-dominated compared to gray matter VASO signal changes for short stimuli only. In a previous study (Dresbach et al., 2023), we found similarly delayed VASO response peaks of superficial, compared to middle and deep layers. Potentially, the current finding indicates that this effect could have been (in part) driven by partial voluming with pial vessels. However, we cannot exclude a gradual process in line with the serial backpropagation proposed by Tian et al. (2010).

Note that the TTP for a stimulus of 2 seconds in Tian et al. (2010) is in the range of 2.5–3.5 seconds in superficial vessels, whereas we found TTPs for vessel-dominated voxels closer to 6 seconds. We believe the discrepancy in TTP can be explained by several aspects. Most importantly, Tian et al. (2010) investigated vessels immediately above the site of activation. On the other hand, we identified vessels close to, but also further away from the activated gray matter (see Fig. 6B for a 3D representation in one representative participant). Assuming that the backpropagation is at play, as proposed by Tian et al. (2010), this would lead to longer times-to-peak in our data when averaging across all identified vessel-dominated voxels. Additionally, Tian et al. (2010) conducted experiments with anaesthetised animals. During anesthesia, vascular autoregulation mechanisms tend to account for the drop in arterial blood pressure by pre-dilating early vascular compartments (during baseline of the task). This means that the task-induced vascular dilation might be differently distributed across fast-acting and slower vascular compartments compared to awake conditions. Furthermore, Tian et al. (2010) employed higher spatial resolution imaging methods and were able to identify individual vessels as arteries or veins. On the other hand, we only delineated vessel-dominated voxels irrespective of an arterial or venous contribution. Therefore, it is possible that we may have identified a substantial amount of veins, which is likely to delay the TTP. Finally, small rodents, like rats in Tian et al. (2010), have shorter arterial arrival times due to smaller brains.

Finally, we find that vessel-dominated activity saturates after 2–4 seconds of stimulation. This large initial VASO response in superficial vessels may result from the actively muscle-controlled dilation of superficial arteries, which reaches its maximum already after short stimulation times (Kennerley et al., 2012). On the other hand, we found that gray matter VASO activity continued to rise with stimulus duration, potentially through the slow and passive dilation of the capillary bed. As a result, the ratio of vessel-dominated over gray matter VASO activation is variable, with highest gray matter weighting for longer stimuli.

As discussed in the previous section, we have explored the combined use of high-resolution functional VASO data with ME-GRE data at 0.35 mm isotropic to investigate CBV and BOLD responses in superficial vessel-dominated voxels—so far without differentiation between arterial and venous contributions, as this was outside the scope of the current investigation. Importantly, previous studies in animals were able to differentiate (intracortical) arteries from veins and showed that CBV responses mainly originate from the arterial vasculature, whereas BOLD responses mainly originate from venous vasculature (Bolan et al., 2006; Chen et al., 2019, 2021; X. Yu et al., 2012, 2016). Therefore, future studies could investigate the contributions of superficial and intracortical arteries and veins to the CBV and BOLD signal acquired with VASO in humans.

Identifying signal changes in the vasculature is, however, hindered by the functional spatial resolution of 0.9 mm isotropic employed in the present study. Specifically, superficial cortical arteries and veins were found to have diameters up to 280 and 350 $μ$m, respectively (Duvernoy et al., 1981).^* Furthermore, intracortical vessels were found to have diameters up to 125 $μ$m and have been shown to be organised in vascular units with central veins (separated by approximately 1–1.5 mm), surrounded by rings of smaller arteries (Bolan et al., 2006; Duvernoy et al., 1981; Harel et al., 2010). Assuming that the arterial rings are roughly centred between two veins, this would yield an artery-vein distance of about 0.5–0.75 mm. As a result, both the superficial and intracortical vessels are significantly smaller and/or closer together than we can resolve with our voxel size of 0.9 mm isotropic. Therefore, partial voluming effects have to be taken into account for functional responses from superficial vessels investigated in the present study, and will limit the distinction between signal changes originating from arteries and veins in future investigations.

Importantly, the VASO signal is usually optimised to obtain functional responses in gray matter and the sequence protocol used here is exceptional in various ways; potentially leading to inflated signal changes above the cortical surface. For example, the use of very short TEs (15 ms) and incorporation of data from multiple sessions may lead to mis-corrections of BOLD contaminations in large superficial vessels due to BOLD-overcompensation and potential misregistration (Chaimow et al., 2019; Devi et al., 2022; Genois et al., 2021; Huber, 2014, 2018a). Furthermore, the strong (flip angle dependent) contrast between CSF and gray matter in our data may lead to apparent signal changes in superficial vessels due to volume redistributions upon activation (Donahue et al., 2006; Jin & Kim, 2010; Piechnik et al., 2009). Finally, the relatively short TR may lead to CBF-dependent outflow of not-inverted water into the post-capillary vasculature, resulting in lower gray matter signals (Wu et al., 2010) and incomplete nulling of draining veins (Huber, 2014 - Thesis chapter 5). Taken together, future studies with shorter experiments to avoid the integration of multiple sessions or different combinations of TEs, flip angles, and TRs will help with the interpretation of the presented results and with determining how far they generalise.

We characterised the spatiotemporal dynamics of CBV and BOLD responses across cortical depth and vessel-dominated voxels at 0.9 mm spatial, 0.785 seconds effective temporal resolution, and $>$7.5 hours of (f)MRI scanning per participant. We demonstrated that it is feasible to obtain reliable VASO responses to stimuli with a duration of 1 second and showed differences in vascular weighting of CBV signal changes depending on stimulus durations, which are fundamental to future neuroscientific applications of VASO. Furthermore, we found differences in signal dynamics across time, space, and stimulus duration for VASO and BOLD, which might inform laminar models of the haemodynamic response. Finally, we found strong VASO signal changes in superficial vessel-dominated, compared to gray matter voxels for short stimuli and stronger gray matter activity for long stimuli, which gives fundamental insights into mechanisms of the VASO contrast.

Analysis code is available on GitHub: $<$https://github.com/sdres/neurovascularCouplingVASO$>$⁠. Raw data are available on OpenNeuro: $<$https://openneuro.org/datasets/ds004943/versions/1.0.2$>$⁠.

According to the CRediT (Contributor Roles Taxonomy) system. Conceptualisation: S.D., R.H., O.F.G., D.I., A.P., R.T., N.W., R.G. Data curation: S.D., R.H., O.F.G., A.P., R.T., D.I. Formal analysis: S.D. Funding acquisition: R.G., N.W. Investigation: S.D., R.H., R.T., R.G., A.P., O.F.G. Methodology: S.D., R.H., O.F.G., D.I., R.G. Project administration: S.D., R.H., R.G. Resources: S.D., R.H., R.T., R.G., N.W. Software: S.D., R.H., O.F.G., R.G. Supervision: R.H., R.T., O.F.G., N.W., D.I., R.G. Validation: S.D. Visualisation: S.D., R.H., O.F.G., A.P., D.I., R.G. Writing—original draft: S.D. Writing—review and editing: S.D., R.H., R.T., O.F.G., A.P., D.I., N.W., R.G.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

We thank Domenica Klank and her team of medical technical assistants at the Max Planck Institute for Cognitive and Neuroscience in Leipzig for their help with participant handling, booking and general assistance while scanning. We thank Rüdiger Stirnberg and Tony Stöcker from DZNE in Bonn for sharing, optimising, and supporting us with the 3D-EPI VASO sequence used here. We thank Vojtěch Smekal, Till Steinbach, and Maite van der Miesen for helpful discussions on the manuscript. Special thanks goes to the remaining “Maastricht layer-seminar” members Lonike Faes, Lasse Knudsen, and Kenshu Koiso for countless discussions on topics related to layer-fMRI. Finally, we thank Kamil Uludag for the inspiring discussions on neurovascular coupling.

S.D. is supported by the “Robin Hood” fund of the Faculty of Psychology and Neuroscience and the department of Cognitive Neuroscience. R.G. is partly funded by the European Research Council Grant ERC-2010-AdG269853 and Human Brain Project Grant FP7-ICT-2013-FET-F/604102. O.F.G. is funded by Brain Innovation. N.W. has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 616905; from the European Union’s Horizon 2020 research and innovation programme under the grant agreement No 681094; from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project no. 347592254 (WE 5046/4-2); and from the Federal Ministry of Education and Research (BMBF) under support code 01ED2210. A.P. is funded by the EU-project H2020-860563 euSNN. R.H. was supported by the NIMH Intramural Research Program (#ZIAMH002783) and by NWO VENI project 016.Veni.198.032 during this project.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00263