We address the issue of jointly detecting brain activity and estimating underlying brain hemodynamics from functional MRI data. We adopt the so-called Joint Detection Estimation (JDE) framework that takes spatial dependencies between voxels into account. We recast the JDE into a missing data framework and derive a Variational Expectation-Maximization (VEM) algorithm for its inference. It follows a new algorithm that has interesting advantages over the previously used intensive simulation methods (Markov Chain Monte Carlo, MCMC): tests on artificial data show that the VEM-JDE is more robust to model mis-specification while additional tests on real data confirm that it achieves similar performance in much less computation time.