A convenient way to analyze blood-oxygen-level-dependent functional magnetic resonance imaging data consists of modeling the whole brain as a stationary, linear system characterized by its transfer function: the hemodynamic response function (HRF). HRF estimation, though of the greatest interest, is still under investigation, for the problem is ill-conditioned. In this paper, we recall the most general Bayesian model for HRF estimation and show how it can beneficially be translated in terms of Bayesian graphical models, leading to 1) a clear and efficient representation of all structural and functional relationships entailed by the model, and 2) a straightforward numerical scheme to approximate the joint posterior distribution, allowing for estimation of the HRF, as well as all other model parameters. We finally apply this novel technique on both simulations and real data.