A linear HRF largely simplifies the analysis and interpretation of BOLD fMRI data. It serves as the central assumption of the classic GLM analysis [110, 111] and the rapid event-related fMRI design [109]. Note that, in such fMRI analyses and designs, the linear system is specified to the relationship between BOLD responses and external stimuli (or tasks), while implicitly assuming linear, yet unknown, neural responses to the stimuli. However, a nonlinear stimulus-to-BOLD relationship has been found in a number of studies [117–123]. The observed nonlinearity may originate from neural and/or vascular sources (e.g. neuronal and vascular refractory effects), since the BOLD response is the ultimate consequence of cascaded processes involving the neural response, the neurovascular coupling and the vascular response [119, 123].