Data acquisition. Data were acquired in the Center for Brain Imaging at NYU with a 3 T Siemens Prisma MRI scanner scanner using a 32-channel head-coil. Twenty functional series of 120 volumes were collected for the retinotopic mapping experiment and twelve functional series of 248 volumes were collected for the spatial working memory experiment. Each functional series experiment. Each functional run were acquired with 14 coronal slices and a gradient echo, echo planar sequence with a 128 square matrix, 192 mm field of view, and 2.0 mm slice thickness, leading to a voxel size of 1.5 x 1.5 x 2.0 mm (TR = 1.5 s, TE = 41 ms, flip angle = 66˚, bandwidth = 752 Hz/pixel). Hz/pixels). A partial Fourier factor of 7/8 was used to acquire an asymmetric fraction of k-space and GRAPPA parallel imaging with a factor of 2 was used to reduce acquisition time. The posterior edge of the acquisition volume was aligned in the mid-sagittal plane with the posterior edge of inferior colliculus. We also collected a high-resolution T₁-weighted MPRAGE (spin-echo, TOMMY PARAMS) for projecting results onto anatomy. In addition, for each scanning session we collected a single whole-brain-coverage functional image with the same spatial resolution as the partial-brain coverage (TR = 10.8 s) in order to align the partial-coverage functional images to the whole-brain anatomical images.

Data processing and analysis. Functional data were motion-corrected and co-registered with the anatomical images. For each voxel, the linear trend was removed and the time series converted to percent signal change. For the mapping experiment, we modeled each voxel in terms of a Gaussian receptive field using methods and tools previously described (\citealt*{DeSimone_2015}; \citealt*{DeSimone_2016}). The pRF model estimates provide a description of each voxel's BOLD response in terms of a retinotopic location and extent. To reconstruct the spatial representations from data measured during the spatial working memory task, we used a spatial inverted encoding model [(IEM), \citealt*{Sprague_2013}; \citealt*{Rahmati_2017})]. Each voxel's response was modeled as the weighted sum of 9 information channels representing polar angles equally spaced around the visual field spanning 365˚, with each channel defined as a one-dimensional cosine function. Subjects participated in 192 trials each for the delayed pro- and anti-saccade tasks. To increase the signal-to-noise ratio, we combined trials by computing a three-fold mean trial time-series, reducing the total number of trials to 64 composed of 8 exemplars at each of the 8 memory locations. To estimate the spatial IEM, we used a leave-one-out cross-validation procedure where we trained the model using 63 trials and reconstructed the response for the held-out 64th trial. We defined voxel activity as the mean BOLD response during the delay period after the visual sample was presented but before the initiation of the first saccade in the double-step saccade. Channel weights for the training data were estimated using a general linear model, where the response for each voxel was modeled as a set of regression coefficients given the sequence of stimuli. The voxel activity for the held-out trial was then projected onto these trained channel weights to produce new channel responses, which were then convolved by the basis functions to produce a reconstruction of the spatial response for the held-out trial. We aligned and combined reconstructions across trials for each subject. In addition, we used a bootstrap procedure to repeatedly train and reconstruct the spatial IEM with different arrangements of trials for computing the the three-fold mean time-series. This ensured that any effects were not simply due to bias in the sampling and combination of trials. To quantify spatial tuning, we used the representational fidelity metric (\citealt*{Sprague_2016}).

Retinotopic mapping. We found clear bilateral activation retinotopic activation in the SC. Figure 2A shows the detailed activation maps for a single subject. Activated voxels whose fMRI time-series correlated with the modeled time-series r² > 0.1 are shown. The zoomed activation maps are overlaid on T₁-weighted anatomical images with the inset of the SC highlighted with a red square. The rows show successive coronal slices moving anterior to posterior from the top to bottom rows. The columns show each of the pRF model parameters (polar angle, eccentricity, and pRF size). The topography of the pRF model estimates were congruent with the organization described previously in the macaque and human. We found a complete representation of the contralateral hemifield in each SC, with the a polar angle progression from the upper vertical meridian, the contralateral horizontal meridian, and the lower vertical meridian as one moves from the medial to the lateral extent of each SC.  We found the eccentricity representation organized along the anterior-to-posterior axis, with representations of the fovea at the anterior pole the left and right SC, with a steady progression toward peripheral representation moving towards the posterior pole. The representational axis of receptive field size conformed with eccentricity where smaller receptive fields are found at the fovea near the anterior pole and larger fields near  periphery near the posterior pole. We summarize these findings in Figure 2B, which shows the direct linear relationship between eccentricity and receptive field size—a hallmark of retinotopic organization throughout the visual system.