Measurement of population receptive fields in human early visual cortex using back-projection tomography

Abstract

Abstract

Abstract:

Abstract  Properties of human visual population receptive fields (pRFs) are currently estimated by performing measurements of visual stimulation using functional magnetic resonance imaging (fMRI), and then fitting the results using a predefined model shape for the pRF. Various models exist and different models may be appropriate under different circumstances, but the validity of the models has never been verified, suggesting the need for a model-free approach. Here, we demonstrate that pRFs can be directly reconstructed using a back-projection-tomography approach that requires no a priori model. The back-projection method involves sweeping thin contrast-defined bars across the visual field whose orientation and direction is rotated through 0°–180° in discrete increments. The measured fMRI time series within a cortical location can be approximated as a projection of the pRF along the long axis of the bar. The signals produced by a set of bar sweeps encircling the visual field form a sinogram. pRFs were reconstructed from these sinograms with a novel scheme that corrects for the blur introduced by the hemodynamic response and the stimulus-bar width. pRF positions agree well with the conventional model-based approach. Notably, a subset of the reconstructed pRFs shows significant asymmetry for both their excitatory and suppressive regions. Reconstructing pRFs using the tomographic approach is a fast, reliable, and accurate way to noninvasively estimate human pRF parameters and visual-field maps without the need for any a priori shape assumption.

Introduction

Discovering how the primate visual system mediates vision requires an understanding of how the visual field is encoded and represented in the brain. In early visual cortex, the visual field is topographically represented many times. Each visual-field map is, in turn, composed of an ensemble of receptive fields (RFs). Each neuron has a receptive field, which is a discrete area within which a change in a visual stimulus alters that neuron's activity. RFs for individual neurons can have complex shapes and configurations that include both excitatory and suppressive regions. In nonhuman primates, these measurements have largely been performed using extracellular electrodes (Hubel & Wiesel, 1974; Van Essen, Newsome, & Maunsell, 1984) or intracellular electrodes (Gilbert & Wiesel, 1985; Long & Lee, 2012). Within early visual areas, RFs are smaller near foveal regions and increase in size with eccentricity. Moreover, as one ascends the visual hierarchy from primary visual cortex, the RFs increase in size and complexity (Hubel & Wiesel, 1977; Van Essen & Maunsell, 1983; B. Wandell, 1995). 

The retinotopic organization of visual cortex causes neurons with similar RF locations and sizes to be spatially grouped. Therefore, RFs can also be discerned for large populations of neurons (Victor, Purpura, Katz, & Mao, 1994). In nonhuman primates, population RF (pRF) measurements have been performed using intrinsic optical imaging (Blasdel, 1992; Gibson, Hebden, & Arridge, 2005; Ts'o, Frostig, Lieke, & Grinvald, 1990) and voltage-sensitive dyes (Blasdel & Salama, 1986; Chemla & Chavane, 2010; Shoham et al., 1999). The pRFs observed in these experiments can also have asymmetric shapes and regions of excitation and suppression (Motter, 2009). 

In humans, pRF properties can be estimated using invasive subdural electrocorticography (Harvey et al., 2012; Yoshor, Bosking, Ghose, & Maunsell, 2007), or noninvasively using fMRI (Amano, Wandell, & Dumoulin, 2009; Dumoulin & Wandell, 2008; Harvey & Dumoulin, 2011; Kay, Naselaris, Prenger, & Gallant, 2008; Zuiderbaan, Harvey, & Dumoulin, 2012). All noninvasive human approaches partially capture the pRF shape using a predefined parametric model. These models vary in complexity; responses from a single cortical location can be characterized using a 2-D Gaussian (Amano et al., 2009; Dumoulin & Wandell, 2008), difference of Gaussians (Zuiderbaan, Harvey, & Dumoulin, 2012), or Gabor wavelet pyramids (Kay et al., 2008). The ability of the models to accurately describe the pRFs also varies systematically across the visual hierarchy in healthy subjects. In clinical populations, pRF properties may differ even more, requiring development of novel pRF models—for example, bilateral pRF models mirrored across the vertical meridian (Hoffmann et al., 2012). Consequently, no single pRF model is sufficient. 

Here we demonstrate that pRFs can be directly reconstructed from fMRI data using back projection, a computed-tomography (CT) technique that requires no a priori shape assumption. Consequently, our tomographic approach has the power to guide and validate model-based shape assumptions and to explore properties of asymmetric pRFs. pRF models can be fitted to the same data, allowing comparisons between model-based and model-free estimates of pRF parameters. 

In back-projection CT, a series of 1-D projections is formed from an unknown 2-D object. A set of these projections that encircle the object is called a Radon transform or sinogram. The method assumes that the unknown object is entirely contained within a particular region of visual space, which we will term the field of view (FOV). The sinogram can then be inverted (reconstructed) using a linear technique called weighted back projection to get a direct estimate of the detailed shape of the original 2-D object (Kak & Slaney, 1988; Smith & Webb, 2010). 

Our experiments involved sweeping thin bars of contrast slowly across the FOV at 12 different orientations spanning 0°–165°. The stimulus bar must be thin and move slowly to acquire high-quality data and useful visual-field resolution. By including corrections for hemodynamic response and bar-thickness effects, we interpret the fMRI time series obtained during these sweeps as a set of 1-D projections of the pRF, a sinogram that can then be reconstructed and analyzed. Our new method of measuring pRFs is computationally fast and straightforward. A similar CT approach was recently used to reconstruct RFs using spikes recorded from V1 neurons in monkey V1 (Pipa, Chen, Neuenschwander, Lima, & Brown, 2012). Very recently, an alternate method for tomographic reconstruction of pRFs was presented (Lee, Papanikolaou, Logothetis, Smirnakis, & Keliris, 2013). This method used an algebraic reconstruction scheme rather than the back-projection methods utilized in our approach. 

Our polar-angle and eccentricity maps agree well with the visual-field maps produced by a conventional model-based approach (Dumoulin & Wandell, 2008) from the same data. Tomographic field mapping of the early visual areas is reliable. Within-session positional variability was found to be small, indicating that very small portions of the visual field (≪1°) can be reliably discerned for experimental manipulations. As expected, tomographic calculation of pRFs shows that sizes increase with eccentricity and at successive stages of the visual hierarchy. In comparing pRF size, the methods diverge because the tomographic approach cannot resolve pRFs below our full-width-at-half-maximum (FWHM) resolution limit of ∼1.5°. Nevertheless, larger pRFs can be resolved throughout the visual field. A substantial subset of these pRFs exhibits asymmetry in both their excitatory and suppressive regions. Thus, the back-projection approach shows great promise as a means to noninvasively and reliably estimate the shapes of pRFs in human visual cortex. 

Materials and methods

Subjects

Measurements were obtained from four subjects (males, ages 26–53), who were all trained psychophysical observers with corrected-to-normal visual acuity. Informed consent was obtained from all subjects based on a protocol approved by the University of Texas, Austin, Institutional Review Board. 

Stimulus

Stimuli were generated using PsychToolbox 3 (Brainard, 1997) running in the MATLAB programming environment on a Macintosh Pro computer. The stimulus was projected using a 60-Hz LCD projector onto a rear-projection screen mounted inside the scanner bore 49 cm from the subject's eyes. 

Stimulation was a thin bar (1° wide) containing a moving black-and-white checkerboard pattern whose columns moved up and down by one plaid panel at a rate of 2 Hz (100% contrast, Figure 1A). The bar swept slowly (0.25°/s) across a 12° FOV. Motion of the bar was perpendicular to its long axis. The motion direction and bar orientation of each sweep rotated through angles 0°–165° in 15° increments for a total of 12 bar orientations. Six bar-sweep orientations were measured per run, so that on odd runs the angles 0°, 30°, …, 150° were collected, and on even runs the angles 15°, 45°, …, 165° were traversed. 

Figure 1

 

Example pRF reconstruction for a voxel in area V2d. (A) Stimulus was a thin, flickering-checkerboard bar that swept across the FOV at 12 angles uniformly spanning 0–165°. (B) Initial representation of the projections, the raw sinogram. (C) A corrected sinogram can be obtained from the raw sinogram by using a Wiener filter to remove bar-width and HRF effects. (D) Estimated pRF reconstructed from the sinogram by weighted back projection. A contour was interpolated from the pRF at half maximum (dashed line). (E) A still from the movie illustrating the 12 bar sweeps from two runs and how the corresponding backprojection reconstruction is built-up from the sequence of projections (see Supplemental data for full movie clip).

An “off time,” consisting of mean-luminance periods (10 s), was inserted after each sweep. The choice of off-time duration between sweeps is a trade-off between data-collection efficiency and some limited reconstruction artifacts for peripheral pRFs located near the boundary of the visual FOV. On all sweeps, the hemodynamic response evoked by pRFs near the center of the FOV will largely return to baseline by the end of the sweep, and little or no off time will be required. However, for some sweep directions, peripheral pRFs will evoke a hemodynamic response that must be allowed to subside to avoid artifacts. To preserve the quality of peripheral pRFs, we provided a 10-s off time, which is sufficient for the majority of the hyperoxic peak of the response to decay before the next sweep is initiated. 

Each session included 12 tomography runs, to create six repeats of each sinogram measurement. 

During each session, the hemodynamic response function (HRF) was also measured using a brief 1.7-s full-field flickering dot stimulus with a 28-s interstimulus interval, 18 events per run, and one or two runs per session. Flickering dots were used to avoid the strong afterimages introduced by brief exposures to typical high-contrast plaid stimuli. 

To control the effects of visual attention and encourage stable fixation during fMRI and HRF runs, subjects performed a fast-paced color-detection task upon the fixation point (0.2° diameter). The dot color changed randomly every 300 ms among 10 colors uniformly distributed along RGB coordinates. Subjects were instructed to push a button each time the fixation point changed to the particular target color. Subjects correctly detected ∼96% of the target appearances. 

MRI protocols

fMRI images were obtained using a GE 3T (Signa Excite HD) scanner with a custom-designed seven-channel surface coil array (ScanMed, Omaha, NE) constructed within a flexible coil former so that it closely and comfortably fit the back of each subject's head. This array produces images with a signal-to-noise ratio (SNR) advantage of 2–4 over the GE eight-channel product coil, depending upon depth. Head movement was minimized with padding. The prescription consisted of 16 2.2-mm-thick slices (140-mm field of view) oriented orthogonal to the calcarine sulcus, covering the majority of early visual cortex. We used a two-shot outward-spiral acquisition (Glover, 1999b; Glover & Lai, 1998) to obtain an in-plane pixel size of 2.2 mm, with echo time (TE) = 30 ms. Acquisition bandwidth was limited to 41.67 kHz to improve SNR and reduce peak gradient current that causes unwanted heating on our scanner. We chose repetition time (TR) = 1 s, so with two shots a volume was acquired every 2 s. 

A set of T₁-weighted structural images was obtained on the same prescription at the end of the session using a three-dimensional RF-spoiled GRASS (SPGR) sequence. These images were used to align the functional data to a structural 3-D reference volume, which was acquired for each subject in a separate session. The structural reference volume was T₁ weighted with good gray-white contrast and was acquired using a 3-D, inversion-prepared SPGR sequence (minimum TE and TR, inversion time = 450 ms, 15° flip angle, isometric voxel size of 0.7 mm, two excitations, duration ∼28 min). These volumes were segmented using ITKSnap and FreeSurfer to delineate the white and gray matter (Fischl, 2012; Yushkevich et al., 2006), and surface models were constructed at the gray-white interface. The gray matter of posterior occipital lobe was also flattened (B. A. Wandell, Chial, & Backus, 2000), to facilitate visualization. 

Image analysis

Analysis of the fMRI data was done using the mrVista software package (http://white.stanford.edu/software) as well as additional tools developed on the mrVista framework in our lab. The first 8 s and 28 s of data were discarded for the moving-bar and HRF runs, respectively, to reduce onset transient effects. We then estimated in-scan motion using a robust scheme (Nestares & Heeger, 2000). Between-run motion was corrected using the same intensity-based scheme, this time applied to the temporal average intensity of the entire scan. The last run of the session was used as the reference. After motion correction, the repeated runs of each condition (odd angles, even angles, and HRF) recorded during each session were averaged together to improve SNR. The intensity of the averaged data was spatially normalized to reduce the effects of coil inhomogeneity. The normalization used a homomorphic method, that is, dividing by a low-pass-filtered version of the temporally averaged volume image intensities with an additive robust correction for estimated noise. A form of high-pass filtering was used to remove the low-frequency drift common in fMRI data. 

A distance map was calculated for gray-matter tissue voxels and the vertices of the surface constructed at the gray–white interface. We used these distances to measure laminar position in the reference volume. Thus, the functional data at each volume voxel were now associated with a depth coordinate (Khan et al., 2011; Ress, Glover, Liu, & Wandell, 2007). For each point on the cortical surface, we could use the depth associations to average the time series over the entire thickness of the gray matter. This procedure was utilized in the tomographic data to reduce the computational intensity of bootstrapping procedures (see later). 

Model-based analysis

The mrVista package includes software for model-based analysis of any collection of fMRI time-series data. The software accepts a description of the stimulus aperture and a description of the HRF. It then assumes a simple parametric model for the pRF, which we chose to be a circular Gaussian form characterized by three parameters: center position (x₀, y₀) and a size scale factor σ which can also be expressed as an FWHM diameter d_FWHM = 2σ Display Formula . The model-based method then assumes that the fMRI response can be predicted by multiplying the stimulus aperture by the model pRF and convolving with the HRF. A combination of gridded search and nonlinear optimization is then used to find the best parameters to match the fMRI time-series data obtained from each voxel (Dumoulin & Wandell, 2008). 

Tomographic back-projection analysis

Theory

Assume that the neuronal population response n(t) to a time-varying two-dimensional visual stimulus s(x, y, t) is given by  

 

where p(x, y) describes the population receptive field (pRF), the spatial-response pattern of the neuronal population. Assume further that the stimulus is constant along the y-axis, so that s(x, y, t) = s(x, t), and define the projection of the pRF as  

 

then  

 

Now let us define the stimulus as a bar of contrast with its long axis along y and width w, moving with speed v:  

 

where R is the radius of the stimulus aperture. At t = 0 and a sweep angle of 0°, the bar starts from the left edge of the aperture and moves to the right. Define a spatial variable x₁ = vt − R and, noting that Rect is an even function,  

 

 

 

where the ∗ operator denotes convolution and the time dependence is now implicit in x₁. If we further assume that the fMRI response f(t) is given by the convolution of the neuronal response with a hemodynamic impulse response function (HRF) h(t), then  

 

 

 

where b(t) = h(t) ∗ Rect(x₁/w) is a blur kernel that includes the effects of both the finite bar width and the HRF. The same relationship holds for any sweep angle θ, with the pRF projection more generally defined as  

 

where the (x′, y′) coordinates correspond to the original coordinates rotated by angle θ. 

We can mitigate the blurring effects from b(t) using a Wiener filter in the Fourier-transform domain:  

 

where k_w is the noise variance, B(ω) is the Fourier transform of time-domain variable b(t), and similarly for other quantities; and ω is temporal frequency. Applying the Wiener filter to the blurred projection data gives an estimate of the projected pRF in the frequency domain:  

 

Thus, for each raw projection obtained from the fMRI data, we can use the Wiener filter in the frequency domain to mitigate the HRF and bar-width effects and then take the inverse transform to estimate the true projection of the pRF. By forming a discrete series of projections around the pRF on the semiopen interval [0°, 180°), we then obtain an estimate of the Radon transform of the pRF, a sinogram that can be inverted using various means, most commonly weighted back projection with a Ram-Lak filter (Kak & Slaney, 1988). 

In Equation 9, the point-spread function (PSF) is given by the first term:  

 

This result in the transform domain can be inverted and back-projected to yield an estimate of the PSF in the spatial domain. Thus, for a particular set of bar-width, HRF, and noise-level parameters, one can quantify the spatial resolution of the system as the full width at half maximum of the PSF. 

Reconstruction

In our experiment, the unknown 2-D object is the pRF and the set of 1-D projections are derived from the fMRI time series collected at 12 different sweep angles. The sweeping bar of contrast (Figure 1A) acts as a line integral, averaging the pRF along its long axis; this is assumed to be the neural population response. We then assume that the measured fMRI time series is a convolution of the stimulus-evoked neural response with the HRF, which is measured in each session. The time series was upsampled to a standard 32 samples per sweep to create a “raw” sinogram (Figure 1B). As discussed previously, the raw sinogram was passed through a Wiener filter that incorporated HRF waveforms measured during each scanning session to remove the blur caused by the HRF and the finite bar width, as well as the usual ramp weighting necessary for back projection (Kak & Slaney, 1988). For illustration, Figure 1C shows the corrected sinogram without the ramp filtering; this is the usual concept of a Radon transform of the unknown pRF. pRFs were then reconstructed from the corrected weighted sinograms using back projection (Figure 1D). Contours were created around the pRF half maxima, and an ellipse was fitted to each contour (Fitzgibbon, Pilu, & Fisher, 1999). Animation shows the pattern of bar sweeps and how the corresponding back-projection reconstruction is constructed from each succeeding projection sweep (Figure 1E). 

Resolution and noise

Our visual-field resolution is determined by the point-spread function (PSF) of the reconstruction (see the Theory subsection previously), which can be quantified by its FWHM diameter. The resolution is affected by the step size of the bar across the visual field, the speed of the bar's movement, and the noise variance k_w chosen for the Wiener filter. Reconstruction noise can be quantified by examining the fraction of variance explained by the reconstruction in the raw sinogram, the set of fMRI time series. To perform this comparison, each reconstruction is multiplied by the stimulus aperture then convolved with the HRF to create a “reconstruction” time series. This is then correlated with the original raw sinogram. We report the mean variance explained for all voxels above our data-quality threshold (see later). 

Visual-field parameters

After ellipses were fitted to the contours, the position of the pRF in Cartesian coordinates was taken as the center of the ellipse. We estimated the size of the pRF as the equivalent circular diameter d = 2 Display Formula , where a and b are the major and minor radii of the ellipse, respectively. We used the elliptical fit to quantify the asymmetry of each contour, defining an aspect ratio AR = a/b. 

Sinogram data quality

The quality of the fMRI data varied substantially from voxel to voxel, so we used correlation upon the mean as a relative metric of tomography data quality. For each raw-sinogram time-series data set, we calculated the correlation coefficient upon the mean sinogram time series. In other words, because a pair of runs generates each sinogram, we calculated the correlation of the mean time series across all adjacent pairs of runs with time series from each individual pair. Because the mean is obtained from only a small number of repeats, typically n = 6, the correlation coefficients obtained in this fashion have a minimum value 1/ Display Formula , but they nonetheless serve to provide a relative quality metric so that we could choose only the higher quality voxels. Accordingly, comparisons were performed using the top 50% of the data based upon the sinogram quality metric for the tomography data and the fractional variance explained for the model-based data. 

Positional reliability

To estimate the reliability of the center position for the tomographic data, each set of six sinograms collected in a session was individually reconstructed into a set of pRFs, one for each voxel. Thus, for each fMRI voxel we obtained a set of six reconstructed pRF images. Next, we randomly resampled these pRF images with replacement from the set of six and then averaged across the resampled set. Each resampled average pRF was then analyzed as usual to compute x and y field-map positions, yielding an estimated sampling distribution for position data and their 68% confidence intervals, p_x and p_y. We chose 68% confidence intervals because they are analogous to the standard error of the mean. We then formed an overall positional variability p = Display Formula . Because this resampling was applied to an entire data set consisting of more than 10,000 voxels, we used only 500 resampled examples, but this is satisfactory to accurately estimate 68% confidence intervals. 

Size versus eccentricity

Because the model-based analysis utilizes a mixture of subsampling and nonlinear methods, a scheme to estimate error distributions by bootstrapping was not attempted. Instead, we used a spatial estimate of the sampling error. Plots of size versus eccentricity were generated by gridding the data into uniform 1° bins. Confidence intervals were estimated by bootstrapping the regridded size data (resampling 10,000 times with replacement). The 68% confidence intervals from the bootstrapped data were then corrected for multiple comparisons based on the oversampling of the reference-volume voxel size (0.7 mm) relative to the acquired fMRI voxel size (2.2 mm), and a further estimate of the spatial correlation of fMRI noise of 1.5, which was estimated by assuming a 1.4-mm-diameter point-spread function for the blood oxygen-level dependent (BOLD) response (Kruger & Glover, 2001; Thompson, Peterson, & Freeman, 2005). This procedure was utilized for both tomographic and model-based size measurements to permit comparison. 

Contour reliability

To estimate the spatial reliability of the half-maximum contours for the tomographic data, we used another bootstrapping approach. For each randomly resampled pRF average, a contour was generated at half maximum. This procedure was repeated 10,000 times to create a large ensemble of contour paths, and we then used a gridding procedure to create a contour-density map. Specifically, we counted the number of contours that passed through each pixel in the map. We then normalized the density map by dividing it into 36 10° angular bins encircling its centroid. Within each angular bin, we found the maximum value and divided the values of the density map within that wedge by the maximum. The maximum of this normalized map followed a path that was very similar to the mean contour obtained from the ensemble. Contours generated from the normalized density map could then be used to estimate spatial confidence intervals; the 5% contours provided 95% confidence intervals. 

Reconstruction cross validation

The sinogram-by-sinogram bootstrapping schemes described previously provide a direct means to assess the reliability of each reconstruction, but we performed a cross-validation procedure as well. Specifically, we split the data from a particular session into half sets each containing three sinograms, and reconstructed the mean of each trio. We then compared the variance explained by each half-set reconstruction against the other. The procedure was repeated for all possible half-set splits, a 10-fold cross validation. This validation was performed among the reconstructions to directly measure their reliability. We also converted the reconstructions into time series to compare the variance explained by the time series. The results measure the reliability of the reconstruction method for smaller data sets consisting of only three repeats. 

Simulation

We simulated the tomographic process in detail to quantify its performance as various parameters were adjusted. The process begins by taking the sweeping-bar stimulus aperture as a time series of images. At each time point in the simulation, we multiplied this aperture by any particular choice of pRF. The result was a simulated neural response. Next, the stimulated neural response was convolved with a model HRF, which we chose as a popular difference-of-gamma functional form (Glover, 1999a) to yield a simulated time series for all sweeps. These data were then reconstructed by the tomographic methods described above. To illustrate the simulation results, we chose four example shapes: a circular Gaussian, an ellipsoidal Gaussian with 2:1 asymmetry, a difference of Gaussian, and an on-center Gabor patch with an envelope of ∼2° FWHM having 1.3:1 asymmetry. 

We also performed simulated reconstructions of a noisy model pRF. Normally distributed noise was added to the simulated hemodynamic response to produce a variance explained of 0.78, similar to that estimated for our fMRI data. 

Results

Simulations

Simulation results confirm that reasonable choices of pRF shapes can be reconstructed with good fidelity by our methods (Figure 2). The first column shows the source pRFs. The second column shows the reconstruction if only six uniformly spaced sweep angles are used. All of the six-projection reconstructions show a sequence of banding artifacts that are typical for undersampled CT. These artifacts are particularly pronounced for the Gabor pRF, which has the largest amount of orientation-dependent structure. High-resolution (0.9° FWHM) reconstructions obtained from 12-sweep data (third column) show much less artifacts, although some banding is still visible in the far periphery for the Gabor pRF. Finally, 12-sweep data reconstructed at 1.5° resolution show minimal artifacts, but some amount of blurring is now evident. 

Figure 2

 

Simulated pRF tomographic reconstructions of different model shapes. From top to bottom: circular and elliptical Gaussians, difference of Gaussians, and a Gabor. Left column shows model shape, second column shows reconstruction for six sweep angles, and third column for 12 sweep angles, both uniformly distributed on [0, 180°) with a spatial resolution of 0.9° FWHM. Fourth column shows 12-sweep reconstruction with 1.5° resolution.

Reconstruction of a noisy pRF (Figure 3C) is precisely equal to the sum of reconstructions of the noise-free pRF (Figure 3A) with a reconstruction of the noiAe (Figure 3B), confirming the linear treatment of noise by back-projection reconstruction. Correlation of the noisy reconstruction with the noise-free reconstruction yields a value of 0.74, similar to the actual value in the simulated input data. We then converted these reconstructions to simulated time series. Correlation between the noisy and noise-free time series is much higher, 0.92 (Figure 3D). Clearly, the projection process greatly reduces the amount of noise that is evident in the reconstruction. 

Figure 3

 

Response of back-projection reconstruction to additive noise. (A) Reconstruction of a noise-free model pRF. (B) Reconstruction of noise. (C) Reconstruction of the sum of signal and noise. The image in (C) is precisely equal to the sum of (A) and (B). (D) Time series obtained from the noisy (blue) and noise-free (black) reconstructions are very similar (R² = 0.92), illustrating the noise filtering produced by the projection process.

We examined the effects of HRF mismatch by reconstructing simulated data generated from a small Gaussian pRF using an HRF that was temporally shifted by ±2 s from the original HRF (Figure 4). Both mismatches produce artificial asymmetry, blur, and regions of apparent suppression. Because the projections only traverse the range of 0°–165°, the distortions have a near-vertical symmetry axis. Center position is also shifted upward for the advanced HRF and downward for the delayed HRF. It is interesting to note that the delayed HRF causes stronger artifacts than the advanced HRF. 

Figure 4

 

Effect of HRF time to peak on reconstruction of a 1° FWHM circular Gaussian model pRF. “Standard” HRF is a canonical difference-of-gamma HRF with 5-s time to peak. Advanced and delayed HRFs are shifted by ±2 s. Spatial resolution is 0.9° FWHM.

Experiment

We first examined the trade-offs between noise and visual-field resolution by examining spatial resolution and time-series variance explained as a function of the Wiener-filter noise-variance parameter (Figure 5). A broad maximum in the variance explained is generally observed for k_w in the range of 0.01–0.1. The variance explained varies from subject to subject in the range of 0.72–0.88. If k_w is chosen too small, the deconvolution becomes ill conditioned and adds noise to the reconstruction. In this case, the corrected sinograms of many voxels become too noisy to adequately reconstruct the pRF; when contours are constructed around the pRF maxima for such voxels, they tend to show sizes near the resolution limit for that value of k_w. Larger values of k_w increase the diameter of the PSF, reducing the visual-field resolution and decreasing noise but preventing the resolution of the smaller pRFs in early visual cortex. The blur also reduces the variance explained by the reconstruction. Based upon examination of the variance explained and the distribution of pRF sizes, we chose k_w = 0.03 because it provided a satisfactory trade-off between visual-field resolution and noise. The resolution varies somewhat from subject to subject because of variations in their HRFs. The mean resolution for the data across all four subjects was 1.6°. 

Figure 5

 

Resolution and noise in back-projection pRF reconstruction as a function of Wiener-filter noise variance. Solid line (right y axis) shows visual-field resolution; dashed line (left y axis) shows time-series variance explained by the reconstruction (mean across voxels). Subject 3 (blue) had the strongest data and greatest variance explained; Subject 4 (red) had the weakest.

Cross-validation results for the reconstructions measures the reliability of reconstructions based on the means of half of the six sinograms obtained in each sessions. The correlations (variance explained) between reconstructions were as follows: Subject 1, 0.46; Subject 2, 0.47; Subject 3, 0.61; Subject 4, 0.43; mean across subjects, 0.49. The correlations between the time--series data were notably higher: Subject 1, 0.81; Subject 2, 0.82; Subject 3, 0.92; Subject 4, 0.69; mean across subjects, 0.81. 

The eccentricity data (Figure 6A) span 0°–6°, and show the expected posterior-to-anterior progression in both methods. In the polar-angle data (Figure 6B), phase reversals are clearly evident at the vertical and horizontal meridians and delineate early visual-area boundaries on the surface rendering and the flat view of the occipital cortex. The maps derived by the two different methods are quite similar at eccentricities ≤6°, but the model-based method gives larger eccentricity values in the periphery. 

Figure 6

 

(A) Eccentricity and (B) polar angle mapped on an inflated rendering of the left cortical (gray-white) surface of Subject 2 (left two columns) and on a flattened section of visual cortex from the right hemisphere of Subject 1 (right columns). Tomographic data are shown on the inner columns and model-based data on the outer columns. Gray scale on the flat maps shows curvature—dark regions are sulci and light regions are gyri.

To provide a quantitative comparison between the methods, we plotted the model-based positional data against the tomography positional data for V1 and V3 (Figure 7). The x and y positions estimated by the two methods agree well (R² = 0.97; median RMS error = 0.34°) for Subject 1, with similar results for the other subjects (Table 1). Comparison data for polar angle are subject to artificial phase-wrapping errors and are therefore not shown. However, the eccentricity data show that the tomographic field maps systematically give smaller values toward the edge of the FOV, yielding somewhat poorer agreement overall (R² = 0.85; median RMS error = 0.41°). 

Figure 7

 

Quantitative comparison of pRF center position between the model-based and tomography methods for Subject 2. In V1 and V3, there is excellent agreement for both x and y positions except at larger eccentricities, where clipping effects occur. Unity slope line is marked in gray.

Table 1

 

Comparison of tomography and model-based pRF center positions (x, y) and eccentricity ϵ. Average R² and RMS values for all subjects for (x, y) are 0.94 and 0.27°, respectively, and for ϵ, 0.79 and 0.32°.

	V1		V2		V3
	R2	RMS	R2	RMS	R2	RMS
Subject 1
x	0.97	0.09	0.97	0.11	0.96	0.13
y	0.97	0.64	0.98	0.56	0.96	0.55
ϵ	0.85	0.41	0.64	0.49	0.82	0.47
Subject 2
x	0.99	0.15	0.98	0.23	0.91	0.37
y	0.95	0.42	0.97	0.40	0.86	0.49
ϵ	0.87	0.30	0.82	0.26	0.63	0.48
Subject 3
x	0.98	0.12	0.97	0.12	0.93	0.17
y	0.96	0.22	0.98	0.29	0.89	0.38
ϵ	0.91	0.22	0.88	0.25	0.70	0.37
Subject 4
x	0.96	0.08	0.95	0.08	0.89	0.14
y	0.92	0.19	0.94	0.21	0.65	0.38
ϵ	0.88	0.30	0.82	0.26	0.63	0.48

For the tomography reconstructions, the within-session positional variability was very small, far less than 1°, at eccentricities below 6° (Figure 8). Median reliability was 0.13° for Subject 3. However, larger positional variability was evident toward the edges of the FOV; at eccentricities >6° in all visual areas, the variability was >3°. At larger eccentricities, dorsal visual areas exhibited larger variability than ventral areas. Some regions of large variability were also evident at all eccentricities near the vertical meridian in V3d and V3v. Similar results were obtained in the other three subjects. 

Figure 8

 

Within-session positional reliability in V1–V3. Positional errors are small throughout the early visual areas, median error = 0.13°.

The tomographic method often estimates larger pRF sizes (Figure 9A) compared with the model-based method (Figure 9B), particularly at smaller eccentricities (Figure 9C through F). In area V1, the tomographic pRF sizes are close to the estimated resolution limits (marked with arrows) across the full range of eccentricities (Figure 9C, E). Nonetheless, the tomographic data show small but significant (p < 0.005) trends for pRF size increasing with eccentricity on the eccentricity range of 0°–6°. The trends are substantially larger in area V3, but again the sizes are limited by the resolution at small eccentricities. At large eccentricities, the rising size trend reverses; this is likely a clipping effect associated with the finite FOV analyzed by the tomographic reconstruction (see Discussion). In contrast, the model-based approach shows much larger trends for increasing size with eccentricity, particularly in V1. For Subject 1 (Figure 9D), the size estimates from the two methods do not agree except at the larger eccentricities in V3. However, for Subject 2 there is fairly good agreement throughout the eccentricity range in area V3 (Figure 9F). In general, the size estimates agree when both tomography measurements and model predictions are above the resolution limit for the tomography data. The size trends for Subject 1 were similar to those for Subject 3, but in extrastriate areas, the tomographic approach gave smaller sizes, close to the resolution limit (Figure A1). 

Figure 9

 

PRF size comparisons. Size estimates for (A) tomography and (B) model-based analysis on flattened cortical surface for Subject 3, left hemisphere. The size overlay displays only the top 50% of data based upon the sinogram quality metric. (C–D) Size versus eccentricity for tomographic (black) and model-based (gray) size data averaged in polar angle in V1 and V3 for Subject 1. (E–F) Same data for Subject 2. The small arrows on the y-axis show tomographic resolution limits.

We want to use the tomographic approach to discern not only the position of the pRFs but their shape as well. Because the method has limited spatial resolution and is sensitive to clipping effects at the edge of the FOV, we must confine analysis to a particular subset of the reconstructions that can be regarded as well resolved. Specifically, we chose reconstructions with equivalent diameters d > 1.5d_psf to ensure our ability to resolve their shape, as well as visual-field eccentricities below (FOV/2 − d) to minimize clipping effects. Typically, about 30%–50% of the pRFs in areas V1–V3 met these constraints. Such pRFs typically were concentrated in a band of eccentricities of 1°–3.5°. 

We frequently observe asymmetric pRFs in this subset of well-resolved reconstructions (Figure 10). Example pRFs from three subjects (Figure 10A through C) show an asymmetrical excitatory region with a nearby suppressive surround. The contours at half maximum in these data are quite asymmetric (aspect ratios: Figure 7A, 2.04; 7B, 1.82; 7C, 1.72; 7D, 1.75), and 95% confidence intervals on these contours demonstrate that the asymmetry is reliable across runs. Often, the excitatory region and suppressive surround have parallel orientations (Figure 10A, C). Occasionally, the asymmetrical pRFs have different orientations for the excitatory and suppressive regions (Figure 7B). The magnitude of the suppression is often a substantial fraction, as much as 50%, of the excitation, and it is statistically significant across runs (Figure 7A through D, p < 0.0001). These patterns of excitation and suppression generally repeat well across sessions in the same subject (Figure 7C, D). Asymmetry is fairly common within this subset of our data. For areas V1–V3 combined, a substantial fraction of the measured pRFs exhibit asymmetry (Figure 7E, F), with 45% having an aspect ratio above 1.5 and 11% having an aspect ratio above 2. Asymmetry shows little variation across visual areas or subjects. 

Figure 10

 

Asymmetrical pRFs in three subjects. Solid contours are half maximum (blue) and half minimum (yellow) of pRF; dashed lines mark 95% confidence intervals. (A) Subject 1: excitatory region and suppressive surround with parallel orientation. (B) Subject 2: orthogonally oriented excitatory region and suppressive surround. (C–D) Subject 3: parallel pattern of excitation and suppression repeats across sessions. (E) Histogram of pRF aspect ratios (AR) from same session as (C). (F) Histogram of pRF ARs for all subjects; 11% of the ARs are greater 2.

Discussion

The tomographic approach operates by sweeping bars of contrast across the visual field at 12 angles that encircle the FOV to generate a set of 1-D projections of the pRFs, then performing reconstruction using weighted back projection after correction for blur caused by the HRF and stimulus-bar width. Mapping visual cortex using our method is a fast and reliable way to directly measure population receptive-field positions, but the size measurements are subject to well-defined resolution limits. Nonetheless, the tomographic approach provides model-free estimates of the pRF and thus provides a means to test pRF shape assumptions. 

Previous methods estimating pRF properties from fMRI time series relied on fitting an assumed model form to estimate visual-field parameters and maps (Dumoulin & Wandell, 2008). The model-based method is a more precise approach for measuring visual-field maps compared with phase-encoded methods. The good agreement between the tomographic and model-based estimates of the center position of the pRF confirms that the tomographic approach is another effective way to map population receptive-field locations and visual-area boundaries. Furthermore, the mapping process is computationally somewhat faster than the model-based method while yielding potentially richer information about the shapes of the pRFs. 

The greater speed and linear character of the tomographic reconstruction facilitates bootstrapped estimates of the reliability of visual-field position and other parameters. Using these methods, we established that the tomographic approach yields very reliable visual-field map estimates (Figure 5). Variability emerges at large eccentricities, near the lower vertical meridians, and in dorsal visual areas. The observed variability near the meridians might occur for a variety of reasons, such as lateralization of receptive fields and eye movements. pRFs at larger eccentricities in the lower visual field may be more strongly affected by clipping effects because the sweeps tend to terminate in the lower visual field. Interpretation of these pRFs is much more strongly dependent on the no-stimulation period that follows the end of the sweep, and this may explain the reduced positional reliability in the dorsal visual areas. Reversing the sweep direction in one set of the sweeps so that the visual field is more uniformly sampled on a temporal basis could circumvent this dorsal–ventral asymmetry in reliability. Randomizing the ordering of the sweep angles might also be useful to avoid artifacts caused by long-term slow hemodynamic responses. 

Lee and colleagues (2013) demonstrated an alternative scheme for pRF tomography that relied upon a variety of algebraic reconstruction. Like our own, their approach did a good job of determining visual-field maps. However, their reconstructions generally exhibited substantial banding artifacts, which may have partly been the consequence of their use of a lesser number of projections (eight as compared to 12), similar in character to those seen in our own simulation results. Moreover, their work did little to characterize the visual-field resolution of the tomographic technique as we have done here. Algebraic reconstruction techniques can offer smoother and more consistent reconstructions, but generally they are computationally more intensive and offer lower spatial resolution than back projection (Kak & Slaney, 1988; Smith & Webb, 2010). 

We have shown that our method can explain a large fraction of the variance, ranging from 77% to 88%, of the time-series data. Lee et al. (2013) noted that their reconstructions could explain a larger fraction of the time-series variance than model-based methods. They claimed that the greater correlations were a notable feature of tomographic reconstruction. We take issue with this claim. Our cross-validation analysis shows, in particular, that the reliability of the reconstructions is much lower than that of the time series obtained from the reconstructions. The time series are, by construction, projections of the reconstruction, and the projection process greatly filters the data, as was illustrated by our simulation results. Accordingly, it is critical to examine the reliability of the reconstructions themselves, as we have done here using a nonparametric statistical approach. 

Tomographic reconstruction, in general, is much more sensitive to data quality, because the estimation process extracts a full image from the measurements. Noisy input data thus directly result in a noisy estimated pRF. Essentially, tomography forces the reconstruction to extract a much higher dimensional estimate of the pRF than is attempted in the model-based approach using only three to five parameters. In our approach, the use of Wiener filtering provides an adjustable means to reject some of the noise, but such filtering assumes that the noise power spectrum is uniform in the temporal-frequency domain, which is not an accurate description of fMRI time-series noise (Glover, Li, & Ress, 2000; Kruger & Glover, 2001). Similarly, in the algebraic scheme presented by Lee et al. (2013), the regularization of the ill-conditioned matrix inversion provides some measure of noise rejection. However, the statistical assumptions inherent in these regularizers are also not matched to the complex structure of fMRI noise. 

Both model-based and tomographic methods show evidence for clipping effects at the edge of the FOV, particularly in area V3, but they are stronger for the tomographic estimates. The tomographic approach is directly sensitive to clipping of pRFs at the edges of the FOV, because the reconstruction only operates within the FOV; portions of the pRF that lie beyond the FOV are ignored. Thus, large pRFs at the edges of the FOV appear to have much smaller sizes, and the eccentricity of their center position is also underestimated. The model-based approach has the freedom to test possible pRF locations outside the FOV, resulting in better estimates. However, even with model-based methods, estimates at the edge of the FOV are less reliable, as less information is available to the model than for pRFs that are entirely within the visual field. 

The tomographic approach diverges from the model-based approach in measuring sizes of the pRFs. Above the tomographic resolution limit, pRF sizes measured by the two methods show closer agreement; this is typically more apparent at intermediate eccentricities and in higher order visual areas such as V3. At small eccentricities and in area V1, pRF sizes are typically below the resolution limit of the tomographic method, given the available data quality. To further investigate the hypothesis that resolution limits explain the disagreement between size estimates for tomography and the model-based approach, we performed two studies (see Appendix). First, a deblurring calculation to approximately remove the PSF effects demonstrated good quantitative agreement between the tomographic and model-based results. Second, a size study using simulated data matched the character of the observed data and confirmed our methods to estimate spatial resolution in the reconstruction. 

Model-based size estimates suffer in a similar fashion from a reduction in information related to the speed of the bar sweep. If the bar moves faster—that is, step size is increased—across the visual field, less information (smaller signal change) is available to resolve small differences in pRF size. However, some information is always available to differentiate pRF size estimates, so this is not an absolute limit: The reduction in signal change can be overcome with an increased contrast-to-noise ratio. Likewise, model-based size estimates are limited by the stimulus FOV: pRFs much larger than the FOV are less reliable. Again, there is always some information available in the fMRI signal, and this reduction can be overcome with an improved contrast-to-noise ratio. These issues highlight the differences between a tomographic reconstruction approach and a model-based fitting scheme. 

The limits of tomography regarding visual-field resolution and clipping must guide its application in studies of spatial vision. To examine the shape of the excitatory center of each pRF, one must choose those receptive fields that are large enough to be resolved; blurring will cause all small pRFs to appear very nearly circular. Tomography is also likely to be a useful tool for studying the pRF suppressive surround that has been characterized using model-based methods (Zuiderbaan et al., 2012). For all such shape studies, it is necessary to avoid or compensate for data at the edges of the FOV because of clipping effects. When these constraints are applied, fortunately, a substantial fraction of the reconstructed pRFs is still available for study. The majority of pRFs are fairly circular, confirming the efficacy of symmetric models that are often used. 

However, asymmetry is clearly present within the well-resolved subset of our reconstructed data. The observed asymmetry is too large to be the consequence of an HRF mismatch. Asymmetric pRFs tend to exhibit an elongate excitatory subregion with an adjacent, often parallel suppressive subregion. Thus, there appears to be a population of pRFs that is more Gabor-like than Gaussian in character. Such asymmetry may help to explain many previous observations, such as the ability to elucidate visual-stimulus orientations from fMRI data (Swisher et al., 2010), a retinotopic organization to orientation selectivity (Freeman, Brouwer, Heeger, & Merriam, 2011), and a preference for radial pattern orientation in the periphery (Sasaki et al., 2006). 

The tomographic pRF-analysis approach is subject to many possible confounds that require further study. Most importantly, the assumption of spatiotemporal linearity is poorly justified. Strict linearity may not be required. For example, neuronal spiking is clearly nonlinear, but receptive fields can be resolved tomographically with appropriate processing (Pipa et al., 2012). Recent evidence for BOLD imaging in visual cortex suggests a compressive nonlinearity for BOLD spatial summation in early visual cortex, and this nonlinearity varies between visual areas (Kay, Winawer, Mezer, & Wandell, 2013). This form of nonlinearity will cause distortions in pRF size estimates. Detailed size studies using tomographic (or model-based) methods will likely need to compensate by a nonlinear rescaling of the time-series data. This could be accomplished within the linear-systems framework by using a nonlinear kernel to describe the convolution with the stimulus aperture. For example, one could utilize the sum of a decaying exponential with the current uniform Rect function. 

In the present work, we used a measured HRF, generated by brief full-field stimulation, as the critical impulse response needed for reconstruction. Because fMRI responses are only partially linear, a modified HRF may be needed to more accurately reconstruct the projections obtained by each bar sweep. Sweep reversal may be a useful tool to adjust the HRF; sweeps obtained from opposing sides of the FOV should produce identical projections, if the HRF is appropriate. Many other hemodynamic issues could also confound this approach. For example, the linearity of the BOLD response has generally only been investigated for stationary stimuli. The use of moving stimuli intentionally generates slow waves of neural activity across visual cortex, and in fMRI we measure the corresponding slow waves of hemodynamic activity. These hemodynamic waves could be systematically modulated by perfusion variations, and tomographic unfolds will be particularly affected by such nuisance modulations. Accordingly, it will be important to perform experiments that confirm pRF structure. For example, if a subject shows substantial asymmetry in pRFs that sample a portion of their visual field, it would then be useful to test if their visually guided behavior is affected by the asymmetry, using psychophysical methods. 

Conclusion

We have demonstrated an alternative method to directly quantify the locations and shapes of population receptive fields in human visual cortex. The method is fast and does not rely on an a priori shape assumption for the pRF. The tomographic approach provides accurate measurements of visual-field positions, although the visual-field resolution of the tomographic method is limited. Spatial resolutions of ∼1.5° FWHM are feasible across a 12° FOV, and higher resolution may be possible across smaller portions of the visual field or by using still slower bar sweeps. With the model-free method, assumptions about the shape of the pRF can be tested. Furthermore, the new method opens up many new possibilities for experiments that attempt to quantify not only the size of pRFs but details of their structure, such as asymmetries in their excitatory centers and the shapes of their suppressive subregions. 

Supplementary Materials

mov S1 

Acknowledgments

Support was provided by National Science Foundation grant BCS 1063774. We thank the Wandell lab for their assistance in understanding and applying the model-based pRF approach, and Gary Glover for the spiral pulse sequence. 

Commercial relationships: none. 

Corresponding author: David Ress. 

Email: [email protected]. 

Address: Department of Neuroscience, Baylor College of Medicine, Houston, TX. 

Appendix: Size and spatial resolution

For completeness, we first plot size versus eccentricity for Subjects 3 and 4 (Figure A1). In general, these subjects yielded somewhat noisier data than the first two subjects, and Subject 4 noted that he had difficulty maintaining attention on the task throughout the session. Size data from Subject 3 shows a trend similar to that of Subject 1, with virtually all V1 pRFs below the tomographic resolution limit, while V3 pRFs show good agreement between methods at intermediate ranges where the pRFs are both large enough to be resolved, and not much affected by clipping issues at the edge of the FOV. 

Figure A1

 

PRF size comparison. (A, B) Size vs. eccentricity for tomography (black) and model-based (gray) in V1 and V3 for Subject 3. (C, D) Same data for Subject 4. Tomographic resolution limits are delineated by the small arrow on the y-axis.

To further investigate the hypothesis that resolution limits explain the disagreement between size estimates for tomography and the model-based approach, we performed the following procedure to approximately compensate for the blurring that remains after the reconstruction. We assume that the observed size estimates are given by d_tomo ≈ Display Formula , where dpRF is the actual size of the pRF and d_FWHM is the FWHM resolution. We can then roughly deblur the data by solving this approximation for d_pRF ≈ Display Formula . When this procedure is applied to the data (discarding a very small fraction of the data where d_tomo ≤ d_FWHM), we obtain good agreement between the tomographic and model-based size estimates for the data from Subject 2, confirming the resolution limits of the tomographic method (Figure A2). Similar results were obtained for Subjects 1 and 3, but not for Subject 4. 

Figure A2

 

Comparison of size versus eccentricity for V1 and V3 in Subject 2, for deblurred tomographic (black lines) and model-based (gray lines) analysis. Tomographic resolution limits are delineated by the small arrow on the y-axis.