This paper presents an approach to characterizing evoked hemodynamic responses in fMRI based on nonlinear system identification, in particular the use of Volterra series. The approach employed enables one to estimate Volterra kernels that describe the relationship between stimulus presentation and the hemodynamic responses that ensue. Volterra series are essentially high-order extensions of linear convolution or "smoothing." These kernels, therefore, represent a nonlinear characterization of the hemodynamic response function that can model the responses to stimuli in different contexts (in this work, different rates of word presentation) and interactions among stimuli. The nonlinear components of the responses were shown to be statistically significant, and the kernel estimates were validated using an independent event-related fMRI experiment. One important manifestation of these nonlinear effects is a modulation of stimulus-specific responses by preceding stimuli that are proximate in time. This means that responses at high-stimulus presentation rates saturate and, in some instances, show an inverted U behavior. This behavior appears to be specific to BOLD effects (as distinct from evoked changes in cerebral blood flow) and may represent a hemodynamic "refractoriness." The aim of this paper is to describe the theory and techniques upon which these conclusions were based and to discuss the implications for experimental design and analysis.