Blood oxygenation level-dependent (BOLD) signal time courses in functional magnetic resonance imaging are estimated within the framework of general linear modeling by convolving an input function, that represents neural activity, with a canonical hemodynamic response function (HRF). Here we investigate the performance of different neural input functions and latency-optimized HRFs for modeling BOLD signals in response to vibrotactile somatosensory stimuli of variable durations (0.5, 1, 4, 7 s) in 14 young, healthy adults who were required to make button press responses at each stimulus cessation. Informed by electrophysiology and the behavioral task, three nested models with an increasing number of parameters were considered: a boxcar; boxcar and offset transient; and onset transient, boxcar and offset transient (TBT). The TBT model provided the best fit of the group-averaged BOLD time courses based on χ ² and F statistics. Only the TBT model was capable of fitting the bimodal shape of the BOLD response to the 7-s stimulus and the relative peak amplitudes for all stimulus lengths in key somatosensory and motor areas. This suggests that the TBT model provides a more comprehensive description of brain sensorimotor responses in this experiment than provided by the simple boxcar model. Work comparing the activation maps obtained with the TBT model with magnetoencephalography data is under way.