Interpretation of fMRI data depends on our ability to understand or model the shape of the hemodynamic response (HR) to a neural event. Although the HR has been studied almost since the beginning of fMRI, we are still far from having robust methods to account for the full range of known HR variation in typical fMRI analyses. This paper reviews how the authors and others contributed to our understanding of HR variation. We present an overview of studies that describe HR variation across voxels, healthy volunteers, populations, and dietary or pharmaceutical modulations. We also describe efforts to minimize the effects of HR variation in intrasubject, group, population, and connectivity analyses and the limits of these methods.