Classic models of the hemodynamic response treat it as a linear system, an assumption that works well for long stimuli that are spaced far apart in time [1]. However, short stimuli and short interstimulus intervals (ISIs) elicit nonlinear responses: short stimuli cause proportionally larger fMRI responses [2,3], and short ISIs cause a reduction and broadening of the fMRI response [4,5,6]. This nonlinearity can vary substantially across cortical regions [7]. Improvements to high-resolution imaging techniques have led to an increased number of studies focusing on deep brain structures, and animal studies have shown that BOLD responses are different in subcortical structures as well [8,9]. We aimed to determine how hemodynamic nonlinearities vary in human thalamus and brainstem. We focused on the visual system, varying stimulus duration and ISI to compare the resulting nonlinearities in subcortical vs. cortical structures.Four subjects gave informed consent and were scanned on a 7T Siemens scanner with a custom-built 32-channel head coil array. Each session began with a 0.75 mm isotropic multi-echo MPRAGE [10]. Functional runs acquired 38 oblique slices, positioned to capture the superior colliculus (SC), lateral geniculate nucleus (LGN), and calcarine sulcus (primary visual cortex, V1). Functional scans consisted of a single-shot gradient echo SMS-EPI at 1.1 mm isotropic resolution (R=4 acceleration with FLEET-ACS [11], MultiBand factor=2, matrix=174x174, blipped CAIPI shift=FOV/2, TR=1.11 s, TE=26 ms, echo spacing=0.79 ms, flip angle=70°). Each functional run lasted 260 seconds. During stimulus presentation a simple fixation task was continuously administered to minimize eye movement. The first two runs presented a radial checkerboard flickering at 12 Hz for 16 seconds, followed by a blank gray screen for 16 seconds. These runs were averaged together and used to localize the visually driven voxels in the SC, LGN, and V1. SC and LGN ROIs were hand-drawn using the functional maps as guides; the V1 ROI was defined using the V1 labeling generated by Freesurfer, and thresholding the z-statistic map from the localizer runs above z=4. The following runs presented checkerboard stimuli lasting either 0.167, 0.5, 1, 2, or 4 seconds. Half the runs used a long ISI (17–19 s), and half the runs used a short ISI (2–3 s). The trial response to each stimulus type was deconvolved in MATLAB using a finite impulse response (FIR) basis set to model the mean timeseries in each ROI, and the impulse response was modeled in FSL using the linear optimal basis set (FLOBS). The first two functional localizer runs were used to identify visually driven voxels in SC, LGN, and V1 (Fig. 1). The trial response in each ROI was then identified across all remaining runs using FIR basis functions (Fig. 2a-c). Stimulus duration induced a nonlinear trial response in V1, as previously identified [1,2,3], with shorter duration stimuli eliciting larger amplitude responses than expected (2d,g). The nonlinearity was similar in LGN (2e,h), but substantially larger in SC (2f,i), in which the trial response to a 167 ms stimulus was almost as large as the response to a 1 s. stimulus. Using FSL's linear optimal basis set, we estimated the hemodynamic response in each region, and obtained physiologically plausible HRFs for the V1 and LGN ROIs (Fig. 3a,b). The hemodynamic impulse response in V1 was consistently narrower for brief stimuli, and was consistently broader in V1 than in LGN (Fig. 3c; FWHM=5.15 s in V1; 4.52 s in LGN). Comparing runs with short ISIs and long ISIs revealed an additional impact of ISI length on the FWHM of the hemodynamic impulse response, with narrower HRFs occurring in the short ISI condition (Fig. 4; Long ISI condition, LGN 4.56 s, V1=5.32 s; Short ISI condition, LGN=4.34 s, V1=5.23 s). We found that the nonlinearity of the fMRI response increases in brainstem components of the visual pathway, compared to cortical structures. In addition, hemodynamic impulse responses in LGN were narrower than in V1, indicating that the temporal dynamics of the response as well as the linearity of its magnitude may be altered. Further narrowing of the impulse response occurred in both structures in experiments using short ISIs, and in response to stimuli of short duration. Both neural and vascular components could contribute to these differences, as saccade-related neural activity can contribute to response nonlinearity, and future studies could combine electrophysiological recordings with functional imaging to dissociate these factors. Narrower hemodynamic response functions lead to increased fast temporal dynamics in the fMRI signal [12]; these results therefore suggest that fast temporal dynamics could be detected in subcortical structures when employing rapid event-related study designs.