Background: Functional magnetic resonance imaging (fMRI) is commonly used to infer hemodynamic changes in the brain after increased neural activity, measuring the blood oxygen level-dependent (BOLD) signal. An important challenge in the analyses of fMRI data is to develop methods that can accurately deconvolve the BOLD signal to extract the driving neural activity and the underlying cerebrovascular effects. New method: A biophysically based method is developed, which combines an extensively verified physiological hemodynamic model with a Wiener filter, to deconvolve the BOLD signal. Results: The method is able to simultaneously obtain spatiotemporal images of underlying neurovascular signals, including neural activity, cerebral blood flow, cerebral blood volume, and deoxygenated hemoglobin concentration. The method is tested on simulated data and applied to various experimental data to demonstrate its stability, accuracy, and utility. Comparison with existing methods: The resulting profiles of the deconvolved signals are consistent with measurements reported in the literature, obtained via multiple neuroimaging modalities. Conclusions: The method provides new testable predictions of the spatiotemporal relations of the deconvolved signals for future studies. This demonstrates the ability of the method to quantify and analyze the neurovascular mechanisms that underlie fMRI, thereby expanding its potential uses.