Many cognitive neuroscience theories assume that changes in behavior arise from changes in the tuning properties of neurons (e.g., Dosher & Lu 1998, Ling, Liu, & Carrasco 2009). However, direct tests of these theories with electrophysiology are rarely feasible with humans. Non-invasive functional magnetic resonance imaging (fMRI) produces voxel tuning, but each voxel aggregates hundreds of thousands of neurons, and voxel tuning modulation is a complex mixture of the underlying neural responses. We developed a pair of statistical tools to address this problem, which we refer to as NeuroModulation Modeling (NMM). NMM advances fMRI analysis methods, inferring the response of neural subpopulations by leveraging modulations at the voxel-level to differentiate between different forms of neuromodulation. One tool uses hierarchical Bayesian modeling and model comparison while the other tool uses a non-parametric slope analysis. We tested the validity of NMM by applying it to fMRI data collected from participants viewing orientation stimuli at high- and low-contrast, which is known from electrophysiology to cause multiplicative scaling of neural tuning (e.g., Sclar & Freeman 1982). In seeming contradiction to ground truth, increasing contrast appeared to cause an additive shift in orientation tuning of voxel-level fMRI data. However, NMM indicated multiplicative gain rather than an additive shift, in line with single-cell electrophysiology. Beyond orientation, this approach could be applied to determine the form of neuromodulation in any fMRI experiment, provided that the experiment tests multiple points along a stimulus dimension to which neurons are tuned (e.g., direction of motion, isoluminant hue, pitch, etc.).