The hemodynamic response could be different between regions across subjects and it has been previously shown that this might impact the estimates of connectivity obtained between such regions [44, 54]. To assess the laminar variability of the HRF and recover the neural response, we performed blind deconvolution before functional connectivity analysis. To assess the effect of deconvolution, we compared the shape of region-specific HRFs across six layers (Fig. 5). Three parameters of region-specific HRFs we examined were response height, time-to-peak, and full-width at half-max (FWHM). The means and standard deviations of the three parameters were calculated separately for each layer across all subjects. As an illustration, we show the region-specific HRF results for left orbitofrontal cortex (Fig. 5a–d) and primary somatosensory cortex (Fig. 5e–h), which are two of the 68 parcellation regions.