To assess the effects of venous density on blood oxygen level-dependent (BOLD) dynamics in order to better understand variations in hemodynamic timing across the human visual cortex.

We performed MRI acquisitions in 6 healthy participants (5 males and 1 female aged 25.3 +/- 3.9 years) on a 3 Tesla Ingenia scanner (Philips, Netherlands). Participants underwent multi-band fMRI imaging covering the occipital cortex (gradient echo EPI TR/TE = 450/30 ms; FA = 55°; multi-band factor = 3; acceleration factor = 1.1; 64 x 64 x 18 voxels; 3.5 mm isotropic) and maintained fixation while a visual stimulus following a boxcar paradigm (20s active; 20s rest; 10 repetitions) was presented. During the active state, a drifting horizontal sinusoidal grating (7 visual degrees; 3 cycles/degree; 6 cycles/second lateral drift) was presented and its contrast was either 100% or sinusoidally modulated over 20 s (Fig. 1a and b). Each condition was performed in 2 separate imaging runs, for a total of 4 runs per subject. Anatomical T1-weighted imaging (shot interval = 3000 ms; TR/TE = 7.9/3.5 ms; 240 x 240 x 150 voxels; 1 mm isotropic) and susceptibility-weighted imaging with phase difference² (SWIp) using 3 echoes (TR/TE = 31/7.2 ms; ΔTE = 10 ms; 336 x 336 x 185 voxels; 0.7 mm isotropic) were performed in the same session.

Standard image preprocessing was performed using AFNI³ and non-local means denoising was performed using DIPY⁴ on all images. A mask of cortical grey matter in the occipital lobe was obtained using FreeSurfer⁵ and a vein mask was computed using an in-house segmentation tool based on vessel-enhancing diffusion^1,6,7 with a common threshold across all subjects. BOLD timecourses were resampled to a TR of 0.5 s and activation maps were computed for each imaging run. The hemodynamic delay (time by which to shift the contrast modulation function to maximize correlation with the BOLD timecourse) was estimated via cross-correlation in the 3-10 s range. Additionally, the single trial average BOLD timecourse was computed in each voxel for each imaging run. The vein masks, activation maps and delay maps were then registered to native T1-space. Finally, vascular density was computed in each fMRI voxel. Results were analysed separately in activated and deactivated regions of the occipital gray matter and compared using two-tailed paired t-tests.

As expected, both visual stimuli yielded activation in the posterior primary visual cortex (V1) and deactivation in anterior V1 (Fig. 1c). This was observed across all subjects and imaging runs. The BOLD response in the deactivation region peaked prior to that in the activation region for both stimulus types (Fig. 2). Although the observed BOLD temporal dynamics were stimulus-specific, the mean hemodynamic delay difference between activation and deactivation was similar.

Venous reconstruction yielded the typical vein pattern observed in humans¹ (Fig. 3a). Vascular density in the deactivation area was found to be statistically greater than that in the activation area (Fig. 3 b and c). This suggests that vascular density differences might drive BOLD latency differences between the two areas, independent of stimulus used.

Due to SWIp image resolution and susceptibility artifacts effectively limiting our ability to detect small venules, the vein density reported likely reflect the larger pial veins. Moreover, it has been suggested that the presence of large veins may cause apparent deactivation that is not neural in origin, particularly at low spatial resolution.⁸ This might explain the larger venous density observed in deactivation, compared to activation. Further, in accordance with previous works,^1,8 our results suggest that proper interpretation and modelling of BOLD signals require accounting for the venous vasculature.We report an association between higher venous density and shorter hemodynamic delay between activation and deactivation sites in the occipital cortex following visual stimulation. Results were robust across two similar visual stimuli that produced different BOLD dynamics due to differences in contrast modulation schemes. Accounting for the effects of venous density on BOLD timing could help establish region-specific hemodynamic response functions (HRFs) to better model local BOLD signals, notably in activation versus deactivation sites. Additionally, cross-correlation analysis combined with fast multi-band fMRI imaging could serve to compute correlation coefficients using the optimal hemodynamic delay, thereby mitigating the effect of BOLD variability across regions on the extent of uncovered activation⁹. Both approaches could help obtain more robust activation maps and further characterize neurovascular coupling.