Topographical Estimation of Visual Population Receptive Fields by fMRI
Summary
It is important to obtain unbiased estimates of visual population receptive fields (pRFs) by functional magnetic resonance imaging. We use mild regularization constraints to estimate pRF topography without a-priori assumptions about pRF shape, allowing us to choose specific pRF models post-hoc. This is particularly advantageous in subjects with visual-pathway lesions.
Abstract
Visual cortex is retinotopically organized so that neighboring populations of cells map to neighboring parts of the visual field. Functional magnetic resonance imaging allows us to estimate voxel-based population receptive fields (pRF), i.e., the part of the visual field that activates the cells within each voxel. Prior, direct, pRF estimation methods1 suffer from certain limitations: 1) the pRF model is chosen a-priori and may not fully capture the actual pRF shape, and 2) pRF centers are prone to mislocalization near the border of the stimulus space. Here a new topographical pRF estimation method2 is proposed that largely circumvents these limitations. A linear model is used to predict the Blood Oxygen Level-Dependent (BOLD) signal by convolving the linear response of the pRF to the visual stimulus with the canonical hemodynamic response function. PRF topography is represented as a weight vector whose components represent the strength of the aggregate response of voxel neurons to stimuli presented at different visual field locations. The resulting linear equations can be solved for the pRF weight vector using ridge regression3, yielding the pRF topography. A pRF model that is matched to the estimated topography can then be chosen post-hoc, thereby improving the estimates of pRF parameters such as pRF-center location, pRF orientation, size, etc. Having the pRF topography available also allows the visual verification of pRF parameter estimates allowing the extraction of various pRF properties without having to make a-priori assumptions about the pRF structure. This approach promises to be particularly useful for investigating the pRF organization of patients with disorders of the visual system.
Introduction
Functional magnetic resonance imaging (fMRI) measures non-invasively the functional organization of visual cortex at a macroscopic scale (typically on the order of millimeters). Early fMRI retinotopy studies used a coherence measure between stimulus location and elicited BOLD responses4-7. These studies typically did not estimate population receptive field size. Later, Dumoulin and Wandell1 proposed a method to overcome such a limitation by explicitly modeling the pRF location and size, using a linear function of this model to predict the BOLD response. However, one limitation of this pioneering method is that the parametric pRF model has to be chosen a-priori, and may lead to erroneous pRF estimates if it turns out not to be appropriate.
To overcome limitations of the parametric pRF-model method, new methods have been developed recently. These methods directly predict the BOLD response to the stimulus by reconstructing the pRF topography. A method8 proposed by Greene and colleagues reconstructs the pRF topography by back-projecting the BOLD responses to the individual 1D stimulus spaces and building the pRF topography in the 2D stimulus space like a typical computer tomography technique. On the other hand, the method2 proposed by us directly estimates the 2D pRF topography by using linear regression and applying a regularization technique. In this method, the pRF topography is represented as a set of weights which is multiplied by the stimulus to estimate the neuronal population response of a given voxel. Then, the final Blood Oxygen Level-Dependent (BOLD) response evoked by the stimulus is estimated by convolving the neuronal population response and the canonical hemodynamic response function. In order to solve the under-constrained linear system, additionally, ridge regression regularization is used to enforce sparseness (see Figure 1 below). The regularization technique suppresses noise and artifacts and thus allows our method to estimate the pRF topography more robustly.
The topographical methods do not force the pRF shape to have a certain parametric shape, and therefore can uncover the actual pRF structure. An appropriate parametric model can then be chosen based on the pRF topography. For example, the pRF topography can be used to separate the pRF center and surround, and then the subsequent pRF center modeling can be more accurate by minimizing the influence of surround suppression as well as the influence of other potential artifacts arising in areas distant to the pRF center. We have recently performed a quantitative comparison between our method and several other methods that directly (i.e. before estimating the topography) fit isotropic Gaussian1, anisotropic Gaussian, and difference of isotropic Gaussians to the pRF9. It was found that the topography-based method outperformed these methods with respect to pRF center modeling by achieving higher explained variance of the BOLD signal time series.
Accurate estimation of pRF properties in various areas reveals how they cover the visual field and is important for investigating the functional organization of the visual cortex particularly as it relates to visual perception. Properties such as how pRF size changes with eccentricity1,10 and pRF center surround organization9 are well studied in the human literature. The proposed method for estimating the pRF topography results in more accurate pRF parameter modeling and is more likely to reveal unknown regularities, not easily modeled a-priori in the direct parametric models. This approach will be especially suitable for studying pRF organization in patients with visual pathway lesions, for whom pRF structure is not necessarily predictable a-priori. Below is described how to estimate the pRF topography and how to use the topography to model the pRF center.
Protocol
1. Data Acquisition
- Prepare a stimulus protocol that is effective in eliciting a reliable retinotopic visual response as previously described in Dumoulin and Wandell1 and Lee et al.2. However, other well established paradigms are also applicable depending on the specific experimental question to be addressed.
- Present bar stimuli drifting across the screen sequentially along 8 directions of space, in steps of 45 degrees. Ensure that the motion is in synchrony with scanner frame acquisition (TR ~2 sec) so that the bar moves a step once an fMRI frame starts and stays at the new location until the frame ends.
- To measure a correct baseline signal, add epochs without bar stimulation1.
- Define a field of view (10 to 15° radius) in visual angle over which the stimulus is presented. Present moving or flickering checkerboard patterns (checker size = 0.94 x 0.94 deg2, pattern update rate = 250 msec/pattern) within the bar to elicit strong visual responses.
- Input the following specific parameters: 8 evenly spaced directions of motion, bar width equal to 1.875 deg, and bars move by half the bar width per frame (2 sec). Additional details can be found in Lee et al.2.
- Generate a spot (~0.25°) in the screen center on which the subject’s eyes fixate during the experiment. Change color of the spot randomly in time.
- Scan the brain of a subject in an MRI scanner using a typical echo-planar-imaging (EPI) scan that has 192 frames duration (24 frames in each direction of motion). Repeat the scans 4-8 times to increase signal-to-noise ratio.
- Set parameters for the EPI sequence as follows: TR = 2 sec, TE = 40 msec, matrix size = 64 x 64, 28 slices, voxel size = 3 x 3 x 3 mm3, flip angle = 90°, Alternatively, apply sequences with a finer resolution (e.g., 2 x 2 x 2 mm3) or a short TR (e.g., 1-1.5 sec) covering only the visual cortex2.
- Track eye movements with an eyetracker system during functional scans to ensure fixation is maintained to within 1-1.5° of the fixation point.
NOTE: Here, a head-coordinate based eyetracker in a goggle system is used, but other suitable eyetracker systems can be used instead. - Instruct the subjects to fixate the spot on the screen center generated in step 1.3.2. To ensure the subjects are fixating, instruct them to report the color changes of the fixation spot.
- Obtain anatomical scans, at 1 x 1 x 1 mm3 resolution (e.g., T1-MPRAGE; TR = 1,900 msec, TE = 2.26 msec, TI = 900 msec, flip angle = 9°, 176 partitions).
NOTE: These anatomical scans will be used for segmentation as well as for aligning the functional images to the anatomy both within and across scans. For better alignment between functional (EPI) images and the anatomy, obtain also an inplane anatomy scan, with resolution identical to the EPI, using T1-weighted fast spoiled gradient echo (SPGR) sequence1.
2. Data Pre-processing
NOTE: Prior to estimating pRF properties, several typical fMRI data pre-processing steps are needed, such as head motion correction and alignment of functional volumes to the anatomical scan. In this article, all pre-processing, estimation, analysis and presentation of results obtained are performed using the open source MATLAB-based software toolbox VISTA LAB available on the VISTA software site. http://white.stanford.edu/newlm/index.php/Main_Page.
- Load the anatomical scan into MATLAB and prepare a volume anatomy using a function called createVolAnat.
- Segment Gray matter, White matter, and CSF using the function “ItkGray”.
- Prepare functional data by converting DICOM (i.e., raw MRI file format for Siemens) files into NIFTI (i.e., standard functional MRI file format) files, and load data into VISTA using a function called mrInit.
- Correct head-motion and align functional images to the anatomy loaded in step 2.1 using rxAlign based on an affine matrix transformation.
- Average functional motion-corrected scans for improving signal-to-noise ratio by clicking mrVISTA Analysis TimeSeries Average tSeries. Exclude from averaging scans during which eye movements deviates from fixation more than 1-1.5°. If signals from different runs have different dc-drifts, average functional scans after removing the dc-drifts.
- Calculate the mapping coordinates between functional scans and Gray matter and identify corresponding Gray-matter voxels in the functional scans by selecting the following menus: mrVISTA Window Open Gray 3-View Window. Assign BOLD signals in the Gray matter voxels by interpolation, choosing one of the options available in mrVISTA.
3. Estimation of pRF Topography and Parametric Modeling
- Download the code files through the following link: https://sites.google.com/site/leesangkyun/prf/codes.zip, extract the compressed file and place them in a preferred location of the local computer. Add the path of the folder in MATLAB.
- Set the stimulus parameters used in the experiment by selecting the following menus: mrVISTA Analysis Retinotopic Model Set Parameters. Specify the following parameters such as stimulus images, the stimulus size, the canonical hemodynamic function, the frame rate of the fMRI scanner.
- Prior to the pRF estimation, prepare the initial parameter sets (Figure 1B).
- Set the cross-validation sets in “tprf_set_params.m” from the code files. Divide timeseries into at least two subsets (one set for testing and the remaining sets for training) that are long enough for the bar to sweep the entire stimulus space. Alternatively, without averaging scans in step 2.4, validate scans by leaving out one scan for testing and using the remaining scans for training.
- Set a coarse parameter set (λ in Figure 1; λ = [10-2 10-1 1 101 102]) in “tprf_set_params.m”. Then, set a fine scale range ([0.1 0.3 0.5 0.7 0.9 1 3 5 7 9]) in “tprf_set_params.m”.
NOTE: The program uses the coarse set to select the λ resulting in the highest explained variance. Then, the program searches the space around the selected λ using the fine scale range, further refining the selection of λ that yields the highest explained variance. - Set a threshold (0.2) of the explained variance for visually responsive voxels in “tprf_set_params.m”.
NOTE: This threshold is used as the reference for selection of visually responsive voxels. Alternatively, make an ROI for a non-visually responsive region (e.g., by drawing a sphere with a radius of 1 cm in a non-visually responsive brain area), where the threshold can be automatically calculated. - Set a set of thresholds ([0.3, 0.5, 0.7]) for defining the pRF center region in the normalized topography in “tprf_set_params.m” (i.e., [0 to 1] or [-1 to 1] with epochs without bar stimulation in step 1.3.1).
NOTE: From the set of thresholds the program provided selects the “best” threshold, i.e. the threshold that defines a pRF central region for which the pRF center model explains the greatest signal variance. Alternatively, choose a different set of threshold values depending on the characteristics of the topography.
- Execute “tprf_runpRFest.m” calculate the pRF topography (Figure 1) and fit a 2D anisotropic Gaussian. After specifying all parameters described in this protocol, and running the code, obtain the final estimation results.
Figure 1: PRF estimation process. (A) Schematic illustration of the process followed for pRF topography estimation. h(t): hemodynamic response function, A(t): stimulus, m: pRF, Reg: L2-norm regularization. (B) Specific steps for pRF topography estimation and pRF center modeling. The set of parameters required for the estimation is listed in each step. A one-dimensional section of topography and its model are illustrated. Under “Model Fitting”, black and red curves represent the topography and its pRF center model with a center threshold of 0.5, respectively. The blue dashed line indicates a threshold for the pRF central region.
Representative Results
Accurate pRF modeling requires capturing pRF shapes correctly. Without knowing the pRF topography, the selection of circularly symmetric models used in prior studies1,9-11 is a reasonable choice. This is because, if the local retinotopic organization is homogeneous in all directions of visual field, a local population response could be represented as a circularly symmetric cumulative aggregate of neuronal responses. However, our observations demonstrate that this is not necessarily the case (Figure 2). Therefore, observation of the pRF topography can be critical for selecting an appropriate parametric function for a pRF model. This is an advantage of the pRF topography, and so the topography-based models outperform the direct-fit isotropic Gaussian models in pRF center modeling, resulting typically in higher explained variance (Figure 2; see Lee et al.2 for additional comparisons with other models). These examples demonstrate the advantage of estimating the pRF topography prior to fitting the model.
Figure 2: Examples of pRF topography estimation and fit of pRF center models. (A) A typical pRF topography. In the topography, red color indicates the most responsive area, which shows the pRF center lying on the middle right horizontal meridian. In the pRF topography, bar patterns across the pRF center structure with low weights are also sometimes observed. This relates to the fact that the area along the bar aperture passing through the pRF center is also stimulated simultaneously with the pRF center. They are easily eliminated in the thresholding step. (B) Comparison between a previous method (DIG; direct-fit isotropic Gaussian)1 and topography-based pRF center model (T-model). The corresponding percent of explained variance is shown above each model. T-models show higher explained variance in all examples, with more accurate pRF shape capture. See Lee et al.2 for more details and additional examples.
One important requirement is to ensure that the fMRI paradigm used provides good retinotopy data. Then the pRF topography method can be used to estimate retinotopic eccentricity and azimuth maps (Figure 3). These maps show similar basic retinotopic architecture as previous methods1,4-7, but they are more accurate because observation of the pRF topography allows us to better separate the pRF center from the surround and from potential noise or artifacts distant to the pRF center. This, among other things, results in better estimation of the retinotopic maps at high eccentricities (a detailed account of the observed differences can be found in Lee et al.2).
Figure 3: Retinotopic maps and pRF size. (A) Eccentricity and Polar angle maps in the left hemisphere of a subject. CS indicates the calcarine sulcus. In the right panel of Figure A, the black circle indicates a region-of-interest (ROI) from which the voxel whose pRF is illustrated in Figure 4 is taken. (B) Relationship between pRF size and eccentricity. The pRF size increases with eccentricity in visual areas V1-3. This plot is drawn from (A).
The topography-based model (T-model) method can be used to estimate various pRF properties such as pRF size, elongation, orientation, and surround suppression efficiently, without having to test many different parametric models. To aid visualization of such properties, a MATLAB function (tprf_plotpRF.m) is provided that plots the pRF topography, the corresponding pRF center model, and their fit to the raw BOLD signal (Figure 4). Note that in some cases, pRF properties may also be estimated directly from the topography, eliminating the need for pRF modeling.
Figure 4: Demonstration of the MATLAB toolbox developed by the authors. This plot shows the pRF topography and corresponding pRF model fit of a voxel selected by a user. The illustrated voxel was selected from the ROI shown in Figure 3A. raw: actual BOLD response, predt: prediction with the pRF topography, predm: prediction with the pRF center parametric model. Please click here to view a larger version of this figure.
Discussion
This article demonstrates how to estimate the topography of visual population receptive fields in human visual cortex and how to use it to select an appropriate parametric model for the receptive field. For a successful retinotopy, an appropriate stimulation protocol and an efficient analysis method should be selected, and the subject’s experimental parameters (motion and fixation) should be optimized. Bar stimuli moving sequentially across the visual field are an efficient stimulus paradigm for pRF estimation as it generates distinct BOLD responses from distinct stimulus locations. The provided method constructs the pRF topography. Since the problem of pRF estimation is generally under-determined, a mathematical tool called ridge regression3 is used to enforce the reasonable constraint of sparseness on the pRF weight solution. This regularization technique is very effective at estimating the pRF model when the number of observations (time points of the BOLD signal) is considerably smaller than the number of pixels covering the spatial dimension of the stimulus.
This method provides more robust estimation of the pRF center than previous methods. There are several reasons for this: 1) it first segments the pRF central region from the pRF topography and then fits an appropriate model, avoiding potential biases that may influence pRF model fits in direct models (i.e. surround suppression or noise artifacts far from the pRF center). 2) Having the ability to inspect the topography visually gives one the opportunity to validate the performance of the final model fit uncovering systematic errors, as well as 3) the possibility to detect features of the pRF structure that may otherwise go undetected. 4) By constraining the fitting area, this model is less likely to map the pRF inside the border of stimulus presentation incorrectly compared to direct fit models (see Figure 2B). Nonetheless, a user need be aware that the proposed method also has limitations for accurately capturing pRF shape near the stimulus border. This is due to the fact that near the border the bar stimuli activate partial receptive fields belonging to voxels whose pRF center would ordinarily be outside the stimulus presentation region. Any receptive field mapping method would be subject to this problem and show a relative peak at the border unless it can perfectly extrapolate from the part of the receptive field center that is mapped to the whole. Having said that, our method is more accurate than direct-fitting methods1,9, which tend to markedly overestimate the distance to the center of pRFs that lie near the stimulus presentation border (see Figures 5 and 6 of Lee et al.2 for more detail).
As discussed, to construct a robust pRF topography depends on the free regularization parameter, λ (Figure 1), which can be optimized separately of individual voxels, or as a common parameter across all voxels. The regularization parameter influences pRF topography by adjusting the extent of fitting (over-fitting or under-fitting) to the data. While a small λ leads to noisy pRF topographies (i.e., over-fitting) compared to the actual pRF, a large λ suppresses visual responses and thus result in more spread topographies than justified by the actual pRF size (i.e., under-fitting). Selection of the optimal lambda is crucial for successful pRF estimation. We estimated λ’s in different subsets of data and evaluated these estimates using a cross-validation strategy. This minimizes biases in pRF topography estimation. Potential residual biases are further reduced in the pRF center modeling step, where different topography thresholds are explored to select one that results in the highest explained variance (see Lee et al.2).
Finally, the topography approach proposed is computationally efficient. The estimation of pRF topographies over all voxels, including finding the optimal regularization parameter λ, takes only a few minutes in a PC environment. Identifying visually unresponsive voxels at this step excludes them from the more computationally demanding step of pRF-center modeling, further improving efficiency. Perhaps more importantly, investigators no longer need to test multiple different pRF models to find one that fits well, since they can be guided in choosing the appropriate model by the pRF topography.
The method demonstrated in this protocol measures population receptive field topography and uses it to guide population receptive field modeling. This approach reduces the bias present in direct population receptive field mapping methods, resulting in more robust and accurate pRF estimates. It also minimizes systematic errors and allows us to study the functional organization of the visual cortex with higher sensitivity. It is particularly applicable in the case of subjects with lesions of the visual pathways, in whom pRF structure may not be easy to anticipate a-priori.
Disclosures
The authors have nothing to disclose.
Acknowledgements
We thank the VISTA software group (Brian Wandell and associates, at Stanford).
S. S. was supported by McNair 2280403105,NEI R01-EY109272, and NEI R01-EY024019 and as HHMI Early Carrer Award. A. P. and G. K. was supported by the Max-Planck Society, G. K. was supported by the PLASTICISE project of the 7th Framework Programme of the European Commission, Contract no. HEATH-F2-2009-223524.
Materials
| Name of Equipment/Software | Company/Provider | Web address | |
| MRI scanner | Siemens/Philips/GE | ||
| MATLAB | The Mathworks, Inc. | http://www.mathworks.com | |
| VISTA software | VISTA software group | http://white.stanford.edu/newlm/index.php/Software | |
| PsychoToolbox | PsychoToolbox | http://psychtoolbox.org | |
| Eye Tracker(VisuaStimDigital) | Resonance Technology Inc | http://mrivideo.com/ |
References
- Dumoulin, S. O., Wandell, B. A. Population receptive field estimates in human visual cortex. Neuroimage. 39, 647-660 (2008).
- Lee, S., Papanikolaou, A., Logothetis, N. K., Smirnakis, S. M., Keliris, G. A. A new method for estimating population receptive field topography in visual cortex. Neuroimage. 81, 144-157 (2013).
- Hastie, T., Tibshirani, R., Friedman, J. H. . The elements of statistical learning : data mining, inference, and prediction. , (2009).
- Sereno, M. I., et al. Borders of multiple visual areas in humans revealed by functional magnetic resonance imaging. Science. 268, 889-893 (1995).
- Engel, S. A., Glover, G. H., Wandell, B. A. Retinotopic organization in human visual cortex and the spatial precision of functional MRI. Cereb Cortex. 7, 181-192 (1997).
- Engel, S. A., et al. fMRI of human visual cortex. Nature. 369, 525 (1994).
- DeYoe, E. A., et al. Mapping striate and extrastriate visual areas in human cerebral cortex. Proc Natl Acad Sci U S A. 93, 2382-2386 (1996).
- Greene, C. A., Dumoulin, S. O., Harvey, B. M., Ress, D. Measurement of population receptive fields in human early visual cortex using back-projection tomography. J Vis. , (2014).
- Zuiderbaan, W., Harvey, B. M., Dumoulin, S. O. Modeling center-surround configurations in population receptive fields using fMRI. J Vis. , (2012).
- Harvey, B. M., Dumoulin, S. O. The relationship between cortical magnification factor and population receptive field size in human visual cortex: constancies in cortical architecture. J Neurosci. 31, 13604-13612 (2011).
- Haak, K. V., Cornelissen, F. W., Morland, A. B. Population receptive field dynamics in human visual cortex. PLoS One. 7, e37686 (2012).
