Nonlinearity in evoked hemodynamic responses often presents in event-related fMRI studies. Volterra series, a higher-order extension of linear convolution, has been used in the literature to construct a nonlinear characterization of hemodynamic responses. Estimation of the Volterra kernel coefficients in these models is usually challenging due to the large number of parameters. We propose a new semi-parametric model based on Volterra series for the hemodynamic responses that greatly reduces the number of parameters and enables "information borrowing" among subjects. This model assumes that in the same brain region and under the same stimulus, the hemodynamic responses across subjects share a common but unknown functional shape that can differ in magnitude, latency and degree of interaction. We develop a computationally-efficient strategy based on splines to estimate the model parameters, and a hypothesis test on nonlinearity. The proposed method is compared with several existing methods via extensive simulations, and is applied to a real event-related fMRI study.