Despite the interest in simultaneous EEG-fMRI studies of epileptic spikes, the link between epileptic discharges and their corresponding hemodynamic responses is poorly understood. In this context, biophysical models are promising tools for investigating the mechanisms underlying observed signals. Here, we apply a metabolic-hemodynamic model to simulated epileptic discharges, in part generated by a neural mass model. We analyze the effect of features specific to epileptic neuronal activity on the blood oxygen level dependent (BOLD) response, focusing on the issues of linearity in neurovascular coupling and on the origin of negative BOLD signals. We found both sub- and supra-linearity in simulated BOLD signals, depending on whether one observes the early or the late part of the BOLD response. The size of these non-linear effects is determined by the spike frequency, as well as by the amplitude of the excitatory activity. Our results additionally indicate a minor deviation from linearity at the neuronal level. According to a phase space analysis, the possibility to obtain a negative BOLD response to an epileptic spike depends on the existence of a long and strong excitatory undershoot. Moreover, we strongly suggest that a combined EEG-fMRI modeling approach should include spatial assumptions. The present study is a step towards an increased understanding of the link between epileptic spikes and their BOLD responses, aiming to improve the interpretation of simultaneous EEG-fMRI recordings in epilepsy.