Second, the spatial laminar point spread function (SL-PSF) of the BOLD response presents a fundamental stumbling block for gaining laminar specificity in fMRI data. Lower layers always contribute signal to the upper layers, because the intracortical veins (ICV) are perpendicular to the surface, and the draining blood flows along the ICV into pial veins on the pial surface [52]. The interpolation-averaging method, wherein the fMRI volume is interpolated at certain cortical depths and the surface profiles are averaged, has been proposed for addressing this issue [18, 21, 22], but a more precise method to extract laminar signals is needed. As we briefly mentioned in the introduction, this is especially true for gradient echo EPI based fMRI which has a flatter PSF compared to spin-echo based EPI. Therefore, future studies may investigate whether spin echo EPI may be better for FC studies at the laminar level, even with the loss of sensitivity in spin echo compared to gradient echo. Recently, an extension of the Friston–Buxton hemodynamic model, which accounts for blood draining effects by coupling local hemodynamics across layers in dynamic