We propose a new, modular, open-source, Python-based 3D+time fMRI data simulation software, \emph{SNAKE-fMRI}, which stands for \emph{S}imulator from \emph{N}eurovascular coupling to \emph{A}cquisition of \emph{K}-space data for \emph{E}xploration of fMRI acquisition techniques. Unlike existing tools, the goal here is to simulate the complete chain of fMRI data acquisition, from the spatio-temporal design of evoked brain responses to various multi-coil k-space data 3D sampling strategies, with the possibility of extending the forward acquisition model to various noise and artifact sources while remaining memory-efficient. By using this \emph{in silico} setup, we are thus able to provide realistic and reproducible ground truth for fMRI reconstruction methods in 3D accelerated acquisition settings and explore the influence of critical parameters, such as the acceleration factor and signal-to-noise ratio~(SNR), on downstream tasks of image reconstruction and statistical analysis of evoked brain activity. We present three scenarios of increasing complexity to showcase the flexibility, versatility, and fidelity of \emph{SNAKE-fMRI}: From a temporally-fixed full 3D Cartesian to various 3D non-Cartesian sampling patterns, we can compare --- with reproducibility guarantees --- how experimental paradigms, acquisition strategies and reconstruction methods contribute and interact together, affecting the downstream statistical analysis.