The major limitation of fMRI is that it is an indirect measure of neural activity, because it measures changes in blood oxygenation level that is modulated by the local vascular distribution (vessel size) and the activation-induced hemodynamic changes [21]. The BOLD fMRI signal can be modeled as the result of the convolution of a latent neural response and the hemodynamic response function (HRF). At the voxel level, the HRF varies across brain regions as well as across individuals [44, 45]. Some animal and human MRI studies at high field have shown that the response height and time-to-peak (TTP) of the HRF varies with cortical depth [18, 46–50]. It was shown that the deeper layers have faster and narrower hemodynamic response compared to the superficial layers. In addition, at the laminar level, gradient-echo BOLD signals have relatively poor laminar specificity, because they are more sensitive to larger vessels [51]. However, a recent investigation of the spatial point spread function (PSF) for the BOLD response showed that the laminar PSF of the gradient-echo BOLD signal had a “flat-tail” characteristic across layers, with
