Improving the reliability and accuracy of population receptive field measures using a logarithmically warped stimulus

Abstract

The population receptive field (pRF) method, which measures the region in visual space that elicits a blood-oxygen-level-dependent (BOLD) signal in a voxel in retinotopic cortex, is a powerful tool for investigating the functional organization of human visual cortex with fMRI (Dumoulin & Wandell, 2008). However, recent work has shown that pRF estimates for early retinotopic visual areas can be biased and unreliable, especially for voxels representing the fovea. Here, we show that a log-bar stimulus that is logarithmically warped along the eccentricity dimension produces more reliable estimates of pRF size and location than the traditional moving bar stimulus. The log-bar stimulus was better able to identify pRFs near the foveal representation, and pRFs were smaller in size, consistent with simulation estimates of receptive field sizes in the fovea.

Introduction

The population receptive field (pRF) method is a powerful tool for investigating the functional organization of human sensory cortex with functional magnetic resonance imaging (fMRI) (Dumoulin & Wandell, 2008). When applied to visual cortex, the region of visual space that elicits a response in an fMRI voxel is estimated by measuring the blood-oxygen-level-dependent (BOLD) activity generated by slowly moving bars of flickering checkerboards and finding the location and size of the Gaussian in visual space that, after accounting for hemodynamic delay, best predicts the time course of each voxel. The pRF is thought to represent the collective response of the many thousands of neurons in that voxel. The pRF method has been used in various ways, such as examining spatial attention dynamics (Kay, Weiner, & Grill-Spector, 2015; Klein, Harvey, & Dumoulin, 2014; Sprague & Serences, 2013), changes due to neural plasticity (Baseler et al., 2011; Haak, Cornelissen, & Morland, 2012; Huber et al., 2019), and assessing clinical diagnoses (Barbot et al., 2021; Brewer & Barton, 2012; Schwarzkopf, Anderson, de Haas, White, & Rees, 2014). 

The reliability and validity of pRF estimates are still a matter of active research. Previous work has shown that pRF polar angle and eccentricity estimates are, in general, highly reliable (Benson et al., 2018; Senden, Reithler, Gijsen, & Goebel, 2014; van Dijk, de Haas, Moutsiana, & Schwarzkopf, 2016; Zeidman, Silson, Schwarzkopf, Baker, & Penny, 2018). However, pRF size estimates vary dramatically across different stimulus configurations (Alvarez, de Haas, Clark, Rees, & Schwarzkopf, 2015; Binda, Thomas, Boynton, & Fine, 2013; Infanti & Schwarzkopf, 2020; Linhardt et al., 2021) and model fitting procedures (Lage-Castellanos, Valente, Senden, & De Martino, 2020; Zeidman et al., 2018). PRF size estimates in foveal regions have particularly poor model fits and are substantially less reliable than in other regions of cortex (Linhardt et al., 2021; Schira, Wade, & Tyler, 2007). 

The unreliability of some pRF size estimates seems to be related to the large cortical magnification factor and small pRF sizes near the fovea (Alvarez et al., 2015; Linhardt et al., 2021). The cortical magnification factor (CMF) is the amount of visual cortex devoted to a fixed area of visual space (Daniel & Whitteridge, 1961). In all retinotopic visual areas, CMF decreases from the fovea to the periphery. Also, pRF sizes vary inversely with CMF, such that CMF ≈ k/σ_pRF, so that the smallest pRFs are found near the fovea where the CMF is largest (Harvey & Dumoulin, 2011). Roughly seven times as much cortex represents the central 1 degree of visual space as a 1° patch of cortex at 12° in the periphery. 

For a given fMRI voxel, the optimal speed and width of a moving bar stimulus depend on the size of the pRF in that location. If the bar moves too quickly, it will sweep through the pRF of a voxel too quickly and produce little or no response due to the slow hemodynamics (blurred over 3 to 5 seconds) and discrete temporal sampling (typically 1 to 2 seconds) of fMRI. However, using a very slowly moving bar will take a prohibitively long time to sample the entire visual field. Because cortical magnification and pRF sizes vary throughout the visual field, a bar of fixed width and speed will sweep through pRFs approximately seven times faster in the fovea than in the periphery. Thus, using a bar of constant width and speed cannot be optimal for the entire visual field. 

One method to address this issue has been to conduct separate runs using bars with varying speeds and widths (Benson et al., 2012; Linhardt et al., 2021). This method successfully estimates smaller pRFs near the fovea. However, these methods require long acquisition times because of the need to acquire data for multiple bar configurations. A second approach has been to design multifocal stimuli scaled to cortical magnification—resembling a black and white dartboard with segments pseudorandomly alternating between black and white ( Binda et al., 2013). Although the multifocal stimulus successfully produces small pRF estimates, the significantly reduced signal-to-noise ratio for multifocal stimuli greatly increases the acquisition time required to obtain reliable estimates (Binda et al., 2013; Vanni, Henriksson, & James, 2005). Finally, Alvarez et al. (2015) compared two CMF-scaled stimuli: (a) a combined wedge and ring and (b) a logarithmically scaled bar stimulus. The combined wedge and expanding/contracting ring produced pRFs similar in size to the fixed-width drifting bar stimulus. The logarithmic bar produced smaller pRFs, but again suffered from low signal-to-noise ratios. This may be because the entire bar width was scaled based on the most central location of the bar, so the bar region that extended into the periphery was optimized for more foveal regions. 

Here, we introduce a novel log-bar stimulus: a moving bar logarithmically warped along the eccentricity dimension. This results in a slowly traveling, thin bar near the fovea with a faster traveling wider bar out in the periphery. Because of the logarithmic nature of the mapping from visual space to cortex (Schwartz, 1980), the log-bar stimulus is designed to approximate a constantly moving fixed-sized bar on the cortical surface whose speed and size can be optimized across the entire visual field simultaneously. 

Finally, using a simple cortical model, we closely replicated fMRI results for both fixed-bar and log-bar stimuli, including the finding that the log-bar recovered pRF parameters more accurately and reliably than a fixed-sized bar stimulus, especially for more foveal pRFs. 

Methods

Participants

Data were acquired from 12 participants (six females, six males), ages 20 to 38 years (M = 28.67, SD = 5.66). All participants were paid and had normal or corrected-to-normal acuity. All gave informed consent in accordance with the human participant Institutional Review Board at the University of Washington, in adherence with the tenets of the Declaration of Helsinki. 

MRI acquisition

MRI data were acquired at the University of Washington Center for Human Neuroscience on a Siemens 3T Prisma MRI scanner with a 64-channel phase-array head coil (Siemens Medical Solutions, Malvern, PA). The protocol included three-dimensional (3D) volumetric navigator-guided T1-weighted (T1w, multiecho) and T2-weighted (T2w) images acquired at 0.8 × 0.8 × 0.8-mm isotropic resolution (field of view [FOV] = 166.4 × 240 × 256 mm³, matrix size = 208 × 300 × 320) (Tisdall et al., 2012). 

BOLD data were acquired using T2*-weighted 2 × 2 × 2-mm isotropic multiband gradient-echo (multiband acceleration = 4) echo-planar imaging sequences (repetition time [TR]/echo time [TE] = 1200/30 ms; flip angle = 64°; FOV = 212 × 212; matrix size = 106 × 106; 56 oblique axial slices). Each run collected 305 volumes (approximately 6.1 minutes), and we collected a total of 12 runs with six runs per stimulus condition across two sessions for each participant. Additionally, a pair of opposite phase-encoded echo-planar imaging (EPI) references were acquired for each session. 

Preprocessing

Data were processed using the fMRIPrep 21.0.1 pipeline (Esteban et al., 2019). The following details are based on the documentation for fMRIPrep. The T1w image was corrected for intensity non-uniformity with N4BiasFieldCorrection (ANTs 2.3.3) (Tustison et al., 2010) and used as the T1w reference throughout the workflow. The T1w reference was then skull-stripped with a Nipype implementation of the antsBrainExtraction.sh workflow. Brain tissue segmentation of cerebrospinal fluid (CSF), white matter (WM), and gray matter (GM) was performed on the brain extracted T1w using the FSL FAST Segmentation Tool (Zhang, Brady, & Smith, 2001). Brain surfaces were reconstructed using FreeSurfer recon-all (Dale, Fischl, & Sereno, 1999). A B0-nonuniformity map (or field map) was estimated based on two EPI references with FSL topup (Andersson, Skare, & Ashburner, 2003). A field map was estimated for each session, and the correction was applied to all functional data within the collected session. 

BOLD runs were slice-time corrected to 0.549 second (50% of the slice acquisition range 0–1.1 seconds) using AFNI 3dTshift (Cox & Hyde, 1997). The BOLD reference was then co-registered with the T1w reference using FreeSurfer bbregister, which implements boundary-based registration (Greve & Fischl, 2009). Co-registration was configured as a rigid-body transformation with six degrees of freedom. The BOLD time series was then resampled onto participant-specific surface space in a non-gridded resampling performed using FreeSurfer mri_vol2surf. 

Head-motion parameters with respect to the BOLD acquisitions were estimated and expanded by including temporal derivatives and quadratic terms for each parameter (Satterthwaite et al., 2013). Additionally, a set of physiological regressors was extracted to allow for component-based noise correction (CompCor) (Behzadi, Restom, Liau, & Liu, 2007). After high-pass filtering, the preprocessed BOLD time series (using a discrete cosine filter with 128-second cut-off), principal components were estimated for the anatomical variant (aCompCor). The aCompCor components were calculated separately within the WM and CSF masks. 

All BOLD data underwent temporal high-pass filtering (Fourier generalized linear model [GLM]; cut-off: three cycles) and were denoised with 24 head-motion parameters (rigid body, temporal derivative, quadratic terms) (Satterthwaite et al., 2013) and the first five aCompCor (WM + CSF) components (Behzadi et al., 2007). 

Stimuli

Stimuli were presented on a BOLDscreen 32 LCD for fMRI (Cambridge Research Systems, Rochester, UK) operating at a resolution of 1920 × 1080 at 60 Hz with a viewing distance of 133 cm. A Mac Mini Intel Core i7 (Apple, Cupertino, CA) running MATLAB R2022a (MathWorks, Natick, MA) controlled custom stimulus presentation code based on Psychtoolbox 3 (Brainard, 1997; Pelli, 1997). Behavioral responses were recorded using a button box (Current Designs, Haverford, PA). Both bar stimuli consisted of an 8-Hz contrast-reversing checkerboard texture. Bars were constrained to a central 16° diameter circular aperture. The display was uniform gray outside the aperture. 

The experiment consisted of two sessions containing 12 functional runs. In each run, two different types of bars were presented (fixed-bar and log-bar). Each run began with a 2-second blank period, followed by eight 45-second sweeps of the bar across the full aperture, each with a pseudorandomly chosen direction that was randomly seeded. Each run ended with a 4-second blank period, resulting in a total duration of 366 seconds per run. The exact stimuli used in the fixed-bar runs were transformed into log-bar stimuli. Figure 1 shows example frames from a single sweep of the fixed-bar and log-bar stimuli used in these experiments. See the OpenNeuro repository for MATLAB/Python code that demonstrates how to warp an image into log-eccentricity space. All participants participated in two 1-hour scanning sessions with up to a 2-week interval between sessions. Each session consisted of three runs of the fixed-bar stimulus and three runs of the log-bar stimulus. 

Figure 1.

 

Example frames from the fixed-bar and log-bar stimuli. Log-bar stimuli were produced by distorting the respective fixed-bar frame along the eccentricity axis, as described in Equation 1. Participants were instructed to fixate on the central cross and report with a button press when the cross changed color.

Fixed-bar stimuli

Fixed-bar stimuli were generated based on the drifting bar from Dumoulin and Wandell (2008). The fixed-bar stimulus was 2° in width and traveled at a constant speed of 0.4°/second. 

Log-bar stimuli

To create the log-bar stimuli, each frame of the fixed-bar stimulus was distorted along its eccentricity axis, expressed as:  

\begin{eqnarray}r^{\prime} = c \cdot \log \left( {1 + k \cdot r} \right), \quad \end{eqnarray}

(1)

where r is the eccentricity (in degrees of visual angle) of each pixel, k is a distortion factor, c is an overall scale constant, and r′ is the mapped eccentricity of that pixel for the log-bar stimulus. Note that this logarithmic warping is consistent with the mapping from visual space to cortical space in the primary visual cortex (Duncan & Boynton, 2003; Schwartz, 1980). Ideally, the log-bar stimulus should produce a wave of activity that moves at a constant speed and width when projected to the cortex. 

The distortion factor, k, determines the strength of the logarithmic warping. In the limit, as k goes to zero, Equation 1 goes toward the identity:  

\begin{equation} \mathop {\lim }\limits_{k \to 0} c \cdot {\rm{log}} \left( {1 + k \cdot r} \right) = r , \end{equation}

(2)

which leads to zero distortion. As k increases, more of the original stimulus is warped toward the fovea. For this experiment, our distortion factor was set at k = 5, which is roughly consistent with the mapping from visual to cortical space in human visual cortex (Duncan & Boynton, 2003). 

The scale factor, c, was set as \(c = \frac{{{r_{max}}}}{{\log ( {1 + k \cdot {r_{max}}} )}}\), where r_max = 8 (the maximum eccentricity of the stimulus). This choice of c ensures that a pixel at 8° in the fixed-bar stimulus is mapped to 8° in the log-bar stimulus. Figure 1 shows example frames from a single sweep of the fixed-bar and log-bar stimuli used in these experiments. As noted earlier, see the OpenNeuro repository for MATLAB/Python code that demonstrates how to warp an image into log-eccentricity space. 

Experimental task

A small cross (“+”) was presented at the center of the display and was present continuously throughout the experiment. The color of the cross switched from green to magenta on average 15 times per run, with a minimum inter-switch delay of 3.2 seconds. Participants were instructed to maintain fixation on the cross and to press a response button whenever the color of the cross changed. The purpose of this task was to encourage fixation for the entire scan duration. 

pRF analysis

Region of interest selection

Initial visual cortex region selection for pRF fitting was performed by aligning each participant's anatomical surface with Benson's atlas of the visual cortex using neuropythy (Benson & Winawer, 2018). The entire visual cortex mask was dilated to ensure that the regions of interest (ROIs) encompassed the edges of each visual area. This process created a visual cortex region of interest for each participant and hemisphere used in the pRF fitting process. 

Benson's atlas was used again after pRF fitting to fine-tune the visual area boundaries by including participants’ pRF information from all stimulus types and runs. Finally, manual examination and editing were performed to create final V1, V2, V3, V3a/b, hV4, lateral–occipital complex (LOC), and temporal–occipital complex (TOC) ROIs for each individual. 

Separately, an anatomical definition of the foveal confluence was drawn by dilating a seed placed at the occipital pole to be at least 350 mm² in surface area. The foveal confluence was then functionally constrained to vertices with pRF eccentricity estimates that were ≤1.5°. This procedure defined the foveal confluence ROIs for each individual. The foveal confluence intersected completely with vertices in V1 (M = 29.18%), V2 (M = 53.75%), and V3 (M = 17.07%). 

Stimulus preprocessing

Prior to model fitting, we binarized the stimulus movies, such that 1 second indicated the presence of the bar and 0 second represented the background. To reduce the computational burden, each frame was spatially linearly downsampled from 540 × 540 pixels to 108 × 108-pixel resolution and temporally linearly downsampled to 366 frames to match the temporal resolution of the fMRI data. 

Hemodynamic response function

We modeled each participant's individual hemodynamic response function (HRF) using the SPM gamma function (Friston et al., 1998). The HRF is characterized by six parameters expressed as:  

\begin{equation}\begin{array}{@{}l@{}} h\left( t \right) = \frac{{{\rm{\beta }}_1^{{{\rm{\alpha }}_1}}t_{\rm{*}}^{{{\rm{\alpha }}_1} - 1}{\rm{exp}} \left( {{-{\rm{\beta }}_1}{t_{\rm{*}}}} \right)}}{{{\rm{\Gamma }}\left( {{{\rm{\alpha }}_1}} \right)}} - \frac{{{\rm{\beta }}_2^{{{\rm{\alpha }}_2}}t_{\rm{*}}^{{{\rm{\alpha }}_2} - 1}\exp \left( {{-{\rm{\beta }}_2}{t_{\rm{*}}}} \right)}}{{c\cdot{\rm{\Gamma }}\left( {{{\rm{\alpha }}_2}} \right)}},\; {t_{\rm{*}}} = t - {\rm{\delta }}, \end{array}\end{equation}

(3)

where δ is the onset delay in seconds, α₁ is the time to peak response in seconds, α₂ is the time to undershoot response in seconds, β₁ is the response dispersion, β₂ is the undershoot dispersion, and c is the response-to-undershoot ratio. An HRF function was estimated for each hemisphere. This estimation resulted in 12 participants × 2 hemispheres = 24 sets of HRF parameters. There were no significant differences in participant HRF parameters. 

pRF modeling

We modeled each vertex as an isometric two-dimensional (2D) Gaussian with three parameters of interest. The model function was expressed as  

\begin{equation}\begin{array}{@{}l@{}} G\left( {x,y;{{\rm{\mu }}_x},{{\rm{\mu }}_y},{\rm{\sigma }}} \right) = \exp \left( { - \left( {\frac{{{{\left( {x - {{\rm{\mu }}_x}} \right)}^2} + {{\left( {y - {{\rm{\mu }}_y}} \right)}^2}}}{{2{{\rm{\sigma }}^2}}}} \right)} \right), \end{array} \end{equation}

(4)

where μ_x and μ_y define the center and σ defines the standard deviation of the pRF 2D Gaussian. 

The pRF method assumes that the fMRI time series for each vertex is the dot product of the stimulus time series convolved with each participant's HRF in time and the pRF of the vertex Gaussian model in space. For each vertex, the pRF model can be formally expressed as  

\begin{equation}\hat{R}\left( t \right) = {\rm{S}}\left( {x,y,t} \right) \cdot {\rm{G}}\left( {x,y;{{\rm{\mu }}_x},{{\rm{\mu }}_y},{\rm{\sigma }}} \right) \circledast h\left( t \right), \end{equation}

(5)

where \(\hat{R}( t )\) is the predicted fMRI time series, S(x,y,t) is the stimulus movie, G(x, y; μ_x, μ_y, σ) is a 2D Gaussian, and h(t) is the hemodynamic response function. 

Using custom MATLAB software, pRF estimates for each vertex were determined to be the values that maximized the correlation between the predicted and actual fMRI time course. The square of this correlation is the proportion of variance in the fMRI time series explained by the pRF model. The pRF estimation yields three parameters: pRF center location (μ_x, μ_y) and pRF size (σ). 

We conducted pRF model fitting in two stages. First, a coarse parameter grid-search was performed with seeds for pRF centers (μ_x, μ_y) linearly sampled from –8° to 8° in 20 steps and seeds for pRF size (σ) linearly spaced from 1 to 5 in 20 steps. The pRF size seeds were chosen to be larger than expected because Lage-Castellanos et al., 2020 reported steeper convergence gradients for pRF size when starting from larger values, irrespective of true pRF size. A predicted time course was generated for each seed and correlated with the actual fMRI time course. The seed parameters with the highest correlation were chosen as the initial starting parameters for the next stage of fitting. 

Second, we estimated each participant's HRF for each hemisphere by holding the pRF parameters from the grid-search fixed and fitting for the six HRF parameters (δ, α₁, α₂, β₁, β₂, c). Then, we held the HRF parameters constant and fit the pRF parameters. This process was repeated for three iterations to ensure that the parameters converged on a stable solution. To limit computation time, the HRF estimation process was carried out by selecting 15% of the vertices out of the subset of vertices with explained variance >20%. The median HRF parameters of these fitted voxels were used to estimate that individual's HRF. This HRF fitting procedure used all collected fMRI data for each participant, regardless of session and stimulus type, to prevent biasing the final pRF estimates toward any particular stimulus type or session. 

After participant HRFs were estimated, the final stage of pRF estimation was performed. The best-starting seed parameters were once again calculated from the coarse grid-search procedure, and then a nonlinear minimization routine (MATLAB fminsearch) was used to find the final pRF estimates. 

The pRF estimates were always fit separately for each stimulus type (fixed-bar and log-bar). We estimated these parameters for each session as well as collapsed across both sessions (session 1, session 2, and sessions 1 and 2). This procedure resulted in 12 participants × 2 hemispheres × 2 stimulus types × 3 session configurations = 144 sets of pRF estimates. 

After pRF fitting, the Cartesian coordinates of the pRF center (μ_x, μ_y) were transformed into polar coordinates (polar angle and eccentricity). Only vertices with estimated pRF centers less than 8° eccentricity (i.e., within the display), that had a pRF size greater than 0.05°, and whose model variance explained was >10% were retained for subsequent analyses. 

pRF model simulations

To gain a deeper understanding of why differences between fixed-bar and log-bar stimuli produced different receptive field sizes, we carried out a simulation experiment. First, we defined a unique set of 4800 pRF locations by sampling polar angle coordinates from 0° to 345° in 15° steps and logarithmically sampling eccentricity coordinates from 0.01° to 8° in 200 steps. pRF size was fixed to be linearly related to eccentricity value, such that σ = 0.15 * eccentricity + 0.1, derived from Himmelberg et al., (2021). Each pRF was defined as a Gaussian, as described in Equation 4. 

Simulated fMRI time courses were then generated for each simulated pRF, for both fixed-bar and log-bar stimuli, as the dot product of the stimulus aperture with the simulated Gaussian pRF followed by a convolution with a realistic hemodynamic response function. Finally, independent and identically distributed (IID) Gaussian noise was added to the simulated time courses with a standard deviation chosen to match the average variance explained of V1 (∼42%) from our fMRI data. These simulated time courses were generated 100 times to create a simulated dataset of similar size as our collected data. The simulated data were fit using the pRF fitting procedures described above. As for the real data, only simulated pRFs with estimated pRF centers less than 8° eccentricity, pRF sizes greater than 0.05°, and variance explained > 10% were retained for all subsequent analyses. 

Results

Behavioral task

Participants’ responses to the color detection task were analyzed across 12 functional runs. Performance was defined as  

\begin{equation} {Performance} = \frac{{{n_{hit}} - {n_{false{\rm{\ }}alarms}}}}{{{n_{total}}}} \cdot 100\end{equation}

(6)

where n_hit is the number of successful color change detections (defined as a button press within 1 second of a color change), n_{false alarms} indicates the number of additional responses without a color change, and n_total indicates the total number of color changes. Average behavioral performance was 96.39%, and the interquartile range of performance values across participants was 95.28% to 98.89%. 

pRF parameter maps

For both fixed-bar and log-bar stimuli, we generated pRF parameter estimates of polar angle, eccentricity, pRF size, and variance explained for all participants. Figure 2 shows parameter maps for polar angle, eccentricity, pRF size, variance explained, and difference in variance explained across stimulus types for an example participant. 

Figure 2.

 

Cortical maps of pRF estimates for an example participant. Solid black lines represent regions of interest: foveal confluence, V1, V2, V3, V3a/b, hV4, LOC, and TOC. pRF estimates were thresholded at >10% variance explained. The variance explained difference (last row) was calculated as fixed-bar – log-bar.

pRF estimate intersession reliability

Next, we compared the intersession reliability of pRF estimates derived from fixed-bar and log-bar stimuli using the procedure described in van Dijk et al. (2016). Briefly, Spearman's rank correlation coefficients were calculated between pRF parameter estimates from session 1 and session 2 of the same stimulus type, retaining vertices that had >10% variance explained in both the fixed-bar and log-bar conditions. For our regression analysis, vertex-by-vertex correlation coefficients (because polar angle is a circular statistic, we used the circular correlation coefficient) for each parameter were transformed into Fisher's z-scores via the arctanh function that converts the non-normal correlation (r) sampling distribution into a standardized statistic:  

\begin{equation*}z{\rm{\ }} = \frac{1}{2}{\rm{ln}}\left( {\frac{{1 + r}}{{1 - r}}} \right)\end{equation*}

 

Figure 3 shows intersession reliability for polar angle, eccentricity, and pRF size estimates for the fixed-bar and log-bar stimulus for the eight regions of interest. We compared the intersession reliability of parameter values for fixed-bar and log-bar stimuli using a mixed-effects ANOVA that treated participants as a random-effect variable and hemisphere and ROI as fixed-effect variables. There was a main effect of stimulus type (fixed-bar vs. log-bar) on polar angle, F(1, 335) = 37.22, p < 0.0001, η² = 0.10; eccentricity, F(1, 335) = 19.03, p < 0.0001, η² = 0.05; and pRF size, F(1, 335) = 109.59, p < 0.0001, η² = 0.25, such that polar angle, eccentricity, and pRF size estimates were more reliable for the log-bar stimulus than the fixed-bar stimulus. 

Figure 3.

 

Intersession reliability for pRF parameter estimates of polar angle, eccentricity, and pRF size across eight regions of interest. Red values represent the fixed-bar stimuli, and blue values represent the log-bar stimuli. Each dot represents a hemisphere, and the solid dots and error bars represent the mean and ±1 SEM across participants’ hemispheres.

Consistent with previous studies (Alvarez et al., 2015; Benson et al., 2018; van Dijk et al., 2016), there was a main effect of ROI on intersession reliability for estimates of polar angle, F(7, 335) = 105.70, p < 0.0001, η² = 0.69; eccentricity, F(7, 335) = 197.11, p < 0.0001, η² = 0.80; and pRF size, F(7, 335) = 47.79, p < 0.0001, η² = 0.50. A Tukey–Kramer post hoc analysis revealed that the foveal confluence ROI was consistently the least reliable for all parameters relative to all other ROIs. 

Finally, there were significant interaction effects on the intersession reliability for pRF size. There was a significant interaction effect between ROI and stimulus type for pRF size, F(7, 335) = 3.10, p = 0.004, η² = 0.06. A Tukey–Kramer analysis showed that this effect was driven by lower intersession reliability for the fixed-bar than the log-bar stimulus in all ROIs except the foveal confluence, where reliability was approximately equal for fixed-bar and log-bar stimuli. There was also a significant interaction effect between ROI and hemisphere for pRF size, F(7, 335) = 3.88, p = 0.0004, η² = 0.07. The post hoc analysis revealed that intersession reliability was greater in the right hemisphere of the foveal confluence and TOC than the left hemisphere counterparts. 

Overall, pRF estimates of polar angle and eccentricity showed very high intersession reliability across all ROIs regardless of whether the fixed-bar or the log-bar stimulus was used. In contrast, the intersession reliability for pRF size was better for the log-bar stimulus than the fixed-bar stimulus in all ROIs. In the foveal confluence, the log-bar and fixed-bar stimuli showed similar intersession reliability, although far more vertices survived the >10% variance criterion, as described in the next two sections. 

Frequency of successfully estimating pRFs

On average, the log-bar stimulus found more vertices that survived the >10% variance explained criterion than the fixed-bar stimulus across all ROIs, as shown in Figure 4. The superior performance of the log-bar stimulus seemed to be confined to the central 2° of visual eccentricity. Within the central 2° of eccentricity, the log-bar stimulus significantly outperformed the fixed-bar stimulus. The crossing of the curves in the fovea confluence ROI at 0.4° eccentricity was likely caused by a bias for the fixed-bar stimulus that shifted estimates of foveal pRFs toward more eccentric locations, as described below. 

Figure 4.

 

The number of vertices whose pRF fits survived the >10% variance explained criterion as a function of eccentricity. The solid lines represent the average values from the fixed-bar (red) and log-bar (blue) stimuli, and the shaded area represents the ±1 SEM across participants’ hemispheres.

pRF parameter estimates

The pRF parameters were estimated separately for each stimulus aperture, collapsing across sessions, to produce pRF estimates for each participant for each hemisphere. 

Variance explained

We used a repeated-measures ANOVA, with participants as a random-effect variable and hemisphere and ROI as fixed-effect variables, to examine the effects of stimulus type on variance explained. There was a significant main effect of ROI, F(7, 163) = 40.70, p < 0.0001, η² = 0.64, on variance explained by the pRF model fitting. Bonferroni-corrected post hoc paired t-tests for each ROI (α′ = 0.01/8 = 0.0013) revealed that the variance explained was significantly greater when using the log-bar stimuli within the foveal confluence, t(22) = −6.95, p < 0.0001, Cohen's d = 1.45. 

Eccentricity

Figure 5 shows pRF eccentricity estimates binned as a function of cortical distance from the foveal center at the occipital pole. We did not plot results in the LOC and TOC areas because they do not share the same foveal center as the other visual areas, so it is not obvious how to define the foveal center in these areas (Schira, Tyler, Breakspear, & Spehar, 2009). In general, there is a tendency for the fixed-bar stimulus to produce larger estimates of eccentricity, and this tendency is pronounced near the fovea: Within the foveal ROI, pRF eccentricity estimates for the fixed-bar stimulus are much larger than for the log-bar stimulus. Indeed, for the fixed-bar stimulus, there were few pRFs with eccentricity estimates < 0.4°, clearly suggesting systematic biases. We used a repeated-measures ANOVA, with participants as a random-effect variable and hemisphere and ROI as fixed-effect variables, to examine the effects of stimulus type on eccentricity estimates. We found that there was a significant main effect of ROI, F(7, 163) = 26.40, p < 0.0001, η² = 0.53, on eccentricity differences. Post hoc paired t-tests revealed that eccentricity values were significantly larger for the fixed-bar stimulus than the log-bar stimulus in all eight regions of interest (see Table 1). 

Figure 5.

 

pRF eccentricity estimates binned as a function of cortical distance from the foveal center. Solid lines represent the average for the fixed-bar (red) and log-bar (blue) stimuli, and the shaded area represents the ±1 SEM across participants’ hemispheres.

Table 1.

 

Summary of post hoc paired t-tests for eccentricity and pRF size differences. PRF estimate differences were calculated as fixed-bar – log-bar estimates. df, degrees of freedom.

 

pRF size

Figure 6 shows pRF size estimates binned across pRF eccentricity. The fixed-bar stimulus produced larger estimates of pRF size in the foveal confluence. In V1 and V2, there is a cross-over near 4° eccentricity, whereas in V3 the log-bar stimulus consistently produced smaller size estimates. We used a repeated-measures ANOVA, with participants as a random-effect variable and hemisphere and ROI as fixed-effect variables, to examine the effects of stimulus type on pRF size estimates. We found that there was a significant main effect of ROI, F(7, 163) = 25.13, p < 0.0001, η² = 0.52, on pRF size differences. Post hoc paired t-tests revealed that pRF size values were significantly larger for the fixed-bar stimulus than the log-bar stimulus in each of the eight ROIs (see Table 1). 

Figure 6.

 

pRF size estimates binned as a function of eccentricity. Solid lines represent the average across participants for the fixed-bar (red) and log-bar (blue) stimuli, and the shaded area represents the ±1 SEM across participants’ hemispheres.

pRF model simulations

Figure 7 compares the ability to recover accurate pRF estimates from fixed-bar and log-bar stimuli. As shown in Figure 7A, consistent with our data in Figure 5, the fixed-bar stimulus overestimated pRF sizes at eccentricities of less than 1.5°. Figure 7B plots these same data as a function of the original pRF sizes used to generate the simulations. Recovered pRF sizes from the log-bar stimuli closely match the simulation pRF sizes across all eccentricities. In contrast, estimates made using fixed-bar stimulus systematically overestimated pRF sizes when the simulated initial pRF values were <0.2°, which in our model represented regions of less than 1.5° eccentricity. As described in the introduction, the size of these biases is determined by the mismatch between the size of the pRF and the width and speed of the traveling bar. 

Figure 7.

 

Simulated pRF size estimates and comparisons to ground truth. (A) Simulated pRF size as a function of eccentricity. (B) Simulated pRF size versus ground-truth pRF size. The colors indicate pRF estimates from the fixed-bar stimulus (red), log-bar stimulus (blue), and ground truth (dashed gray).

Discussion

This study shows how pRF estimates can be improved by using a logarithmically distorted drifting-bar stimulus instead of the standard fixed-bar stimulus. We found that, across visual areas, the log-bar stimulus produced (a) more reliable pRF estimates of polar angle, eccentricity, and size; (b) more pRFs near the foveal representation; and (c) pRFs that were much smaller in size and closer to the fovea (more closely resembling the what is believed to be the underlying ground truth) than those produced with the traditional drifting-bar stimulus. 

The limitations of the current study arise from technical capabilities and design choices. First, we were unable to stimulate the visual field beyond 8° eccentricity due to visual stimulation constraints. Although we were able to obtain reliable pRF estimates in higher visual areas (i.e., V3a/b, hV4, LOC, and TOC), these pRFs were primarily found near the fovea, presumably due to the restricted size of the visual stimulus. Therefore, we were not able to fully consider the effects that the log-bar stimulus would have in more eccentric visual locations in higher visual areas with larger pRFs. 

A second limitation is that we chose to distort the log-bar stimulus with only a single factor (k = 5). This factor was chosen to approximately match CMF properties of human early visual cortex, but a more severe distortion factor might be appropriate for pRF mapping in higher visual areas. A third limitation is that we used a 100% contrast flickering checkerboard as the background texture for our stimuli. Flickering checkerboard stimuli are known to strongly drive neural activity but it may not be the optimal background texture. Visual cortical neurons respond strongly to 1/f spectrum noise (Isherwood, Schira, & Spehar, 2017), and the Human Connectome Project implemented their stimulus with 1/f noise with embedded objects (Benson et al., 2018). It is possible that a combination of an optimized background texture and optimized stimulus aperture, like our log-bar stimulus, could lead to better pRF estimates and model performance, especially in higher order visual areas, although it is not obvious that changing the background texture would affect the pRF parameters of location and size. 

A fourth limitation is that the spatial frequency of the checkerboard background differed between the fixed-bar and log-bar stimuli. The logarithmic distortion used to produce the log-bar stimuli decreased the check sizes toward the fovea. However, although the differing spatial frequency content between the two stimuli affected the ability to compare their effects on pRF maps directly, we argue that smaller checks near the fovea are a feature of the log-bar stimulus, as high spatial frequency is known to drive smaller receptive fields better (Henriksson, Nurminen, Hyvärinen, & Vanni, 2008; Welbourne, Morland, & Wade, 2018). 

Finally, we did not incorporate more complex visual computations in our pRF modeling procedure, such as compressive spatial (Kay, Winawer, Mezer, & Wandell, 2013) or temporal (Zhou, Benson, Kay, & Winawer, 2018) summation, surround suppression (Zuiderbaan, Harvey, & Dumoulin, 2012), or normalization (Aqil, Knapen, & Dumoulin, 2021). Incorporating more complex pRF models would be an interesting future direction. However, we do not have prior reason to believe that more complex pRF models would reduce the advantages of the log-bar stimulus, especially because the advantages of the log-bar stimulus are greatest near the fovea where the simpler pRF models appear to be the most appropriate (Kay et al., 2013). 

In summary, a logarithmically distorted drifting bar stimulus, locally optimized to match the receptive fields of the corresponding region of cortical space, produces more accurate and reliable pRF estimates than a standard drifting bar stimulus, especially for pRFs located near central vision. 

Acknowledgments

Supported by a grants from the National Eye Institute, National Institutes of Health (R01EY014645 to IF) and the University of Washington Center for Human Neuroscience. 

The data and code reported in this article are available in the OpenNeuro repository (https://doi.org/10.18112/openneuro.ds004698.v2.0.0). 

Commercial relationships: none. 

Corresponding author: Kelly Chang. 

Email: [email protected]. 

Address: Department of Psychology, University of Washington, 119 Guthrie Hall, Seattle, WA 98195, USA.