Despite the interest in simultaneous EEG-fMRI studies of epileptic spikes, the link between epileptic discharges and their corresponding hemodynamic responses is poorly understood. We applied two biophysical models in order to investigate the mechanisms underlying the neurovascular coupling in epilepsy: a metabolic hemodynamic model, and a neural mass model that simulates epileptic discharges. Analyzing the effect of epileptic neuronal activity on the BOLD response we focussed on the issues of linearity and on the origin of negative BOLD signals. In our BOLD simulation results both sub- and supra-linearity occur one after another. The size of these effects depends on the spike frequency, as well as on the amplitude of the excitatory part of the neural input. For the hemodynamic model used in this study, we found that the sign of the BOLD response is mainly determined by the area under the curve describing the excitatory neural activity. Therefore, a strong deactivation following the initial peak of the excitatory time course of an epileptic spike is necessary to obtain a negative BOLD