Laminar fMRI using the BOLD contrast enables the non-invasive investigation of mesoscopic functional circuits in the human brain. However, the laminar neuronal activity is spatiotemporally biased in the observed cortical depth profiles of the BOLD signal. In this study, we propose a generative fMRI signal model, comprehensively covering the relationship between cortical depth-dependent changes in excitatory and inhibitory neuronal activity with the sampling of the BOLD signal with finite voxels. The generative model allowed us to investigate pertinent questions regarding the accuracy of the laminar BOLD signal relative to the neuronal activity, and we found that: a) condition differences in laminar BOLD signals may be more reflective of neuronal activity than single condition BOLD signal depth profiles; b) angular dependence of the BOLD signal induces significant signal variability, which can mask underlying activity profiles; c) even if only three neuronal depths are of interest, more BOLD signal depths should be considered in the analysis. In addition, we recommend that the laminar BOLD data should be displayed using the centroid method to appreciate its spatial distribution in the original resolution. Finally, we showed that Bayesian model inversion of the generative model can improve sensitivity and specificity of assessing depth-dependent neuronal changes both for steady-state and dynamically.