1. Neuroscience
Download icon

Population receptive fields in nonhuman primates from whole-brain fMRI and large-scale neurophysiology in visual cortex

  1. P Christiaan Klink  Is a corresponding author
  2. Xing Chen
  3. Wim Vanduffel
  4. Pieter R Roelfsema
  1. Netherlands Institute for Neuroscience, Royal Netherlands Academy of Arts and Sciences, Netherlands
  2. Psychiatry Department, Amsterdam UMC, Netherlands
  3. Laboratory for Neuro- and Psychophysiology, Department of Neurosciences, KU Leuven Medical School, Belgium
  4. Massachusetts General Hospital, Martinos Ctr. for Biomedical Imaging, United States
  5. Leuven Brain Institute, KU Leuven, Belgium
  6. Harvard Medical School, United States
  7. Department of Integrative Neurophysiology, Center for Neurogenomics and Cognitive Research, VU University, Netherlands
Research Article
  • Cited 1
  • Views 842
  • Annotations
Cite this article as: eLife 2021;10:e67304 doi: 10.7554/eLife.67304

Abstract

Population receptive field (pRF) modeling is a popular fMRI method to map the retinotopic organization of the human brain. While fMRI-based pRF maps are qualitatively similar to invasively recorded single-cell receptive fields in animals, it remains unclear what neuronal signal they represent. We addressed this question in awake nonhuman primates comparing whole-brain fMRI and large-scale neurophysiological recordings in areas V1 and V4 of the visual cortex. We examined the fits of several pRF models based on the fMRI blood-oxygen-level-dependent (BOLD) signal, multi-unit spiking activity (MUA), and local field potential (LFP) power in different frequency bands. We found that pRFs derived from BOLD-fMRI were most similar to MUA-pRFs in V1 and V4, while pRFs based on LFP gamma power also gave a good approximation. fMRI-based pRFs thus reliably reflect neuronal receptive field properties in the primate brain. In addition to our results in V1 and V4, the whole-brain fMRI measurements revealed retinotopic tuning in many other cortical and subcortical areas with a consistent increase in pRF size with increasing eccentricity, as well as a retinotopically specific deactivation of default mode network nodes similar to previous observations in humans.

Editor's evaluation

This study is a detailed, systematic comparison of visually evoked population receptive fields (pRFs) measured non-invasively with MRI and invasively with electrophysiology in the same primate species. The authors show that MRI pRFs provide a good estimate of receptive fields based on multi-unit spiking activity in early visual areas. These results make an important contribution to our understanding of human imaging data in research and in the clinic.

https://doi.org/10.7554/eLife.67304.sa0

Introduction

The concept of a receptive field (RF) is crucial for our understanding of the mechanisms underlying perception, cognition, and action. RFs (Hartline, 1938; Sherrington, 1906) typically describe stimulus locations that evoke or modulate neuronal responses, but they can be generalized to different stimulus features such as color or spatial frequency. RFs are usually measured by determining the neuronal firing rate elicited by visual stimuli (Hubel and Wiesel, 1998; Hubel and Wiesel, 1968; Hubel and Wiesel, 1959), but they can also be defined based on other neuronal signals such as subthreshold activity (Priebe, 2008), properties of the local field potential (LFP; Victor et al., 1994), or calcium levels that can, for instance, be measured with fluorescent calcium indicators (van Beest et al., 2021; Bonin et al., 2011).

Noninvasive methods lack the spatial resolution to measure the RF properties of single neurons, but they can characterize the RF properties of the aggregate neural signals being measured. The retinotopic organization of the human brain has now been characterized with functional magnetic resonance imaging for decades (Wandell et al., 2007; Wandell and Winawer, 2011). Early studies used phase-encoding with ‘rotating wedge’ and ‘expanding or contracting ring’ stimuli to identify RF position (Engel, 2012; Engel et al., 1994; Sereno et al., 1995), while later studies increasingly used the ‘population receptive field’ (pRF) method that estimates RF size in addition to position (Dumoulin and Wandell, 2008; Wandell et al., 2007; Wandell and Winawer, 2015; Wandell and Winawer, 2011). The method is popular and has been used to map a range of visual and cognitive functions (Binda et al., 2018; Ekman et al., 2020; Harvey et al., 2020; Harvey et al., 2015; He et al., 2019; Hughes et al., 2019; Mo et al., 2018; Poltoratski et al., 2019; Poltoratski and Tong, 2020; Puckett et al., 2020; Shao et al., 2013; Shen et al., 2020; Silson et al., 2018; Stoll et al., 2020; Thomas et al., 2015; Welbourne et al., 2018; Zuiderbaan et al., 2017), dysfunctions (Ahmadi et al., 2020; Alvarez et al., 2020; de Best et al., 2019; Dumoulin and Knapen, 2018; Green et al., 2019; Schwarzkopf et al., 2014), mechanisms of brain development (Dekker et al., 2019), cortical evolution (Keliris et al., 2019; Kolster et al., 2014; Zhu and Vanduffel, 2019), and information transfer across different brain areas (Haak et al., 2013).

The term ‘pRF’ highlights the analogy to neuronal RFs. It assumes that the blood-oxygen-level-dependent (BOLD) signal measured with fMRI reflects the aggregate response of a large population of neurons within a voxel. Indeed, pRFs of the human visual cortex are qualitatively similar to the RFs of single neurons or multi-unit activity (MUA) in animals (Dumoulin and Wandell, 2008). However, most of the previous comparisons were between species, and between studies that used different techniques to measure pRFs/RFs (Barlow et al., 1966; Hubel and Wiesel, 1968). Exceptions are pRF studies based on intracranial recordings in human patients (Harvey et al., 2013; Winawer et al., 2013), which used a limited number of surface electrodes to measure ECoG or intracranial EEG, but not spiking activity. The pRFs derived from the LFP exhibited similar properties to pRFs derived from BOLD signals, including similar spatial summation characteristics (Winawer et al., 2013). Another study compared BOLD-based pRFs in monkeys to the RF properties of single units in published work (Kolster et al., 2014). Furthermore, Keliris et al., 2019 found that single-unit RFs in one of their monkeys were smaller than pRFs measured with BOLD and proposed another method to estimate RF sizes. Their study included MUA but not the LFP. A systematic within-species comparison of pRFs derived from BOLD, MUA, and LFP has however never been carried out.

Here, we fill that gap with extensive pRF modeling based on BOLD, MUA, and LFP signals in macaque monkeys. The question that neuronal signal forms the basis of the fMRI-BOLD signal has far-reaching consequences for the interpretation of human neuroimaging results and is therefore a topic of ongoing debate and rigorous investigation (Arthurs and Boniface, 2002; Bartels et al., 2008; Boynton, 2011; Drew, 2019; Ekstrom, 2010; Goense and Logothetis, 2008; Logothetis, 2010; Logothetis, 2003; Logothetis et al., 2001; Logothetis and Wandell, 2004; Maier et al., 2008; Schölvinck et al., 2010; Sirotin and Das, 2009; Winawer et al., 2013; Winder et al., 2017). Some studies reported that properties of the BOLD signal resemble features of both neuronal spiking and the LFP (Mukamel et al., 2005; Nir et al., 2007; Rees et al., 2000), others that they resemble the LFP but not spiking (Bartolo et al., 2011; Maier et al., 2008; Niessing et al., 2005; Viswanathan and Freeman, 2007), and yet others that they resemble spiking rather than the LFP (Lima et al., 2014). Here, we examine the degree to which pRFs based on the BOLD signal resemble pRFs based on spiking activity and distinct frequency bands of the LFP (Buzsáki, 2006; Buzsáki and Draguhn, 2004; Einevoll et al., 2013; van Kerkoerle et al., 2014). In nonhuman primates, we measured BOLD-pRFs using whole-brain fMRI and determined neuronal pRFs with large-scale neurophysiological recordings in V1 and V4 (Figure 1). Besides showing the presence of retinotopic information throughout the brain based on the fMRI data, we could directly compare V1 and V4 pRFs obtained with fMRI and electrophysiology. This intraspecies comparison provides new insight into the neurophysiological basis of the BOLD-defined pRFs and offers a benchmark for visual field maps obtained with fMRI.

Figure 1 with 1 supplement see all
Experimental setup and study design.

(A) Monkeys maintained fixation on a red dot while bars with high-contrast-moving checkerboards moved across the screen in eight different directions behind a virtual aperture (dashed line, not visible in the real stimulus). (B) Two animals performed the task in the MRI scanner. Two other animals were each implanted with 16 Utah arrays (1024 electrodes/animal) in the left visual cortex. The approximate locations of 14 V1 arrays (red) and 2 V4 arrays (blue) for one animal are depicted on the NMT standard macaque brain. For more detailed array configurations, see Chen et al., 2020. (C) Four population receptive field (pRF) models were fit to the data, differing in their pRF parameters (location: x, y; size: σ) and spatial summation characteristics. The difference-of-Gaussians (DoG) pRFs are described by an excitatory center and an inhibitory surround (both 2D Gaussians with σ1 and σ2 as size parameters; dark and light gray circles in the left panel, respectively). Other models fitted single Gaussians that were either constrained to be positive (second panel: solid line) or allowed to be negative (unconstrained linear U-LIN, dashed line). The compressive spatial summation (CSS) model implemented nonlinear spatial summation across the receptive field (RF) (fourth panel: dashed line), while all other models implemented linear summation (solid line). (D) The pRF model fitting procedure. A model pRF is multiplied with an ‘aperture version’ of the bar stimulus to generate a predicted response. For fMRI data, this prediction was convolved with a monkey-specific hemodynamic response function (HRF). The difference between the recorded neural signal (blood-oxygen-level-dependent [BOLD], multi-unit activity [MUA], local field potential [LFP]) and the predicted response was minimized by adjusting the pRF model parameters. (E) Examples of data and model fits for a V1 voxel (top panel) and a V1 electrode (middle and bottom panels). Average activity (gray data points) depicts the BOLD signal (top), MUA (middle), and LFP power in the low gamma band (bottom). Black lines are the model fits for a P-LIN pRF model. Light and dark gray areas depict visual stimulation periods (bar sweeps as indicated by the icons above). In the white epochs, the animals viewed a uniform gray background to allow the BOLD signal to return to baseline (these epochs were not necessary in the electrophysiology recordings). Note that that in MRI trials the bar stimuli had a lower speed (1 TR or 2.5 s per stimulus location) than during electrophysiology (500 ms per stimulus location).

The original method of estimating pRFs from BOLD responses (Dumoulin and Wandell, 2008) uses a forward modeling approach to fit the location and size of a symmetrical two-dimensional Gaussian to the BOLD responses. This approach is often used to predict neuronal activity elicited by moving bar-shaped stimuli. The RF model minimizes the difference between measured and predicted responses by multiplying the pRF with the stimulus and convolving the result with a hemodynamic response function (HRF), which accounts for the time course of neurovascular coupling (Figure 1). Later refinements implemented a difference-of-Gaussians pRF profile (DoG) to account for center-surround interactions (Zuiderbaan et al., 2012; Figure 1C), with substantial improvements in early visual cortex. Another refinement is the introduction of a static nonlinearity that models nonlinear spatial summation across RFs (Britten and Heuer, 1999; Kay et al., 2013; Oleksiak et al., 2011; Winawer et al., 2013). In such a model, the best parameters indicate subadditive spatial summation in all visual areas (Kay et al., 2013). This means that if stimulus S1 elicits a response R1 and a nonoverlapping stimulus S2 elicits response R2, the response to the combined stimulus, S1 + S2, is smaller than the sum, R1 + R2. For this reason, the nonlinear spatial summation model has also been called the ‘compressive spatial summation’ (CSS) model. A third extension has been the modeling of negative pRFs. Standard approaches tend to only include voxels that show increases in the BOLD signal in response to a stimulus. The inclusion of ‘negative’ pRFs with decreased BOLD activity has revealed the retinotopic organization of a number of areas in the so-called default mode network (DMN) (Szinte and Knapen, 2020). In our analysis of pRFs based on BOLD, MUA, and LFPs, we explored several pRF models, allowing us to investigate the potential presence of nonlinear spatial summation and negative pRFs.

Results

Four macaque monkeys (Macaca mulatta) participated in this study. They were rewarded with fluid for maintaining their gaze inside a 1.5° window around a fixation point that was presented at the center of a frontoparallel screen. While they fixated, a 2° wide bar containing full-contrast moving checkerboards traversed the screen in eight different directions (Figure 1). Two animals performed this task in a 3T horizontal bore MRI scanner. Two other monkeys were each implanted with 1024 electrodes in the visual cortex (V1, V4). They performed the same task while neuronal activity (MUA and LFP) was recorded simultaneously from all electrodes. In the MRI setup, the stimulus covered 16° of the visual field, which was the maximum possible with the monitor located just outside of the scanner bore. The bar traveled across this screen in 20 steps of 2.5 s (1 TR). In the electrophysiology setup, the monitor was closer to the animal, allowing a visual field coverage of 28°. The stimulus bar moved across this aperture in 30 steps of 500 ms. For both the MRI and electrophysiology recordings, we only included data from epochs when the animals maintained fixation for >80% of the time.

After preprocessing (see Materials and methods), we independently fit four pRF models to the average BOLD time courses. These models were (1) a linear pRF model constrained to have positive responses (P-LIN) (Dumoulin and Wandell, 2008), (2) an unconstrained version of the linear pRF model that can also model negative responses (U-LIN), (3) a DoG pRF model (Zuiderbaan et al., 2012), and (4) a nonlinear CSS pRF model (Kay et al., 2013; Figure 1). We used the fitting method to determine pRF size, shape, and location. A cross-validated goodness of fit (R2) was determined by fitting the model to one half of the data and calculating R2 using the other half of the data. Cross-validation allows the comparison of fit quality between models with different numbers of parameters.

pRFs measured with BOLD-fMRI

All models provided good fits to the BOLD time courses in a range of cortical and subcortical areas known to be involved in visual processing. For both monkeys (M1 and M2), we found robust retinotopic information in occipital, temporal, and parietal cortex (Figure 2). pRFs in all these areas were in the contralateral visual field and retinotopic maps were consistent with previous reports (Figure 2B), some of which were more extensive (Arcaro et al., 2011; Arcaro and Livingstone, 2017; Brewer et al., 2002; Janssens et al., 2014; Kolster et al., 2014; Rima et al., 2020; Zhu and Vanduffel, 2019). Weaker and sparser retinotopic information was also observed in the frontal cortex, for example, around the arcuate sulcus (area 8, including the frontal eye fields [FEF]) and in the ventrolateral prefrontal cortex (VLPFC). Throughout this study, we will use a voxel inclusion criterion of R2 > 5% unless otherwise noted. While R2 was generally much higher in visual areas (Figure 2A, bottom panel), this relatively low threshold also reveals retinotopic information in more frontal areas and some subcortical regions (where the signal picked up by surface coils has a lower signal-to-noise ratio [SNR]). Figure 2C shows the number of voxels within a range of areas for which the models explained more than 5% of the variance (Figure 2—figure supplement 1 shows the proportion per region of interest [ROI]). The functional parcellation of visual areas based on field sign inversions around horizontal and vertical meridians lined up well with a probabilistic atlas, co-registered to the individual animal’s anatomy (D99, Reveley et al., 2017).

Figure 2 with 1 supplement see all
Population receptive field (pRF) model fits and retinotopic maps.

(A) R2 value map of the compressive spatial summation (CSS) pRF model projected on the surface rendering of the brains of two monkeys (M1, M2). The lower panel illustrates that the R2 value in the visual cortex is generally much higher than 5% (going up to ~90%). The color range in the upper panels also reveals the weaker retinotopic information elsewhere in the cortex. AS, arcuate sulcus; CS, central sulcus; IPS, intraparietal sulcus; LatS, lateral sulcus; LuS, lunate sulcus; STS, superior temporal sulcus. (B) Polar angle maps for both subjects from the CSS model, thresholded at R2 > 5%, displayed on the inflated cortical surfaces. Functional delineation of several visual areas is superimposed. FEF, frontal eye fields. (C) Number of (resampled) voxels with R2 > 5% per brain area and subject for each model. Note the logarithmic scale. See Figure 2—figure supplement 1 for proportions per region of interest (ROI), and Supplementary file 1 for a table with all ROI abbreviations.

Subcortically, we could segregate the lateral geniculate nucleus (LGN), pulvinar and some striatal regions from their surrounding areas on the basis of a higher R2. In both monkeys, the LGN of both hemispheres contained clear retinotopic maps (Figure 3A). Retinotopic information was also evident in the bilateral pulvinar of M1, but some of the pRFs in the pulvinar of M2 were large and crossed the vertical meridian, resulting in noisy polar angle maps (Figure 3B). Striatal retinotopy was less pronounced in M1 and even more variable in terms of polar angle maps (Figure 3—figure supplement 1).

Figure 3 with 1 supplement see all
Retinotopy in the thalamus.

Thalamic population receptive fields (pRFs) in M1 (A) and M2 (B). The lateral geniculate nucleus (LGN, top rows) contained retinotopic maps of the contralateral visual field in both monkeys (M1: 38/38; M2: 73/80 voxels with contralateral pRFs). Retinotopic information was also present in the pulvinar (PULV, bottom rows), but its organization was much less structured, especially in M2 (M1: 23/32; M2: 61/131 voxels with contralateral pRFs). Voxels were thresholded at R2 > 3% for these polar angle visualizations due to the generally poorer fits in subcortex compared to visual cortex. Results from the compressive spatial summation (CSS) model are masked by region of interest (ROI) and shown both in a ‘glass’ representation of the individual animals’ brains (left), and on selected sagittal and coronal slices (monkey-specific T1-weighted images, cross-hairs indicate slice positions). Dashed lines indicate the boundaries of the LGN and pulvinar.

We calculated cross-validated R2 values to compare the four pRF models: P-LIN, U-LIN, DoG, and CSS (Figure 1C). A comparison of model performance pooled over subjects and voxels confirmed that there were significant differences across models, with the CSS model outperforming the P-LIN model (Kruskal–Wallis test on all four models, H = 21.33, df = 3, p<0.0001; post-hoc Tukey’s HSD multiple comparisons of mean rank, R2CSS > R2P-LIN, p<0.0001). Indeed, the CSS model fit was better than that of P-LIN in all ROIs in M1 and in 38 out of 39 ROIs in M2 (Wilcoxon signed-rank, p<0.05; no difference in premotor area F7 in M2) (Figure 4). Since the CSS model also provided the best fits for the neurophysiological signals (described later), we report results from this model and only extend the analysis to other models where this is useful (e.g., in the case of negative gain values). The advantage of the CSS model over the P-LIN model generally increased in higher visual areas (Figure 4—figure supplement 1).

Figure 4 with 2 supplements see all
Comparison of the four population receptive field (pRF) models.

(A) Comparison across pRF models. R2 data are in bins of 1% × 1%, and color indicates the number of voxels per bin. The compressive spatial summation (CSS) model fits the data best, while U-LIN and difference-of-Gaussians (DoG) models give better fits for voxels that are poorly characterized by the P-LIN model (arrows). Dashed triangle, voxels for which R2 of the U-LIN/DoG models was at least 5% higher than that of the P-LIN model. (B) Gain values of the U-LIN pRF fits for all voxels (black) and the voxels in the gray triangle in (A.) The negative gain values indicate visual suppression of the blood-oxygen-level-dependent (BOLD) signal. Arrows, medians. (C) Clusters of voxels for which the U-LIN model fit better than the P-LIN model were located in the medial occipital lobe (left panel). The pRFs of these voxels had a high eccentricity according to the P-LIN model, beyond the stimulated visual field. (D) Clusters of voxels with negative pRFs in the U-LIN model (without positive pRFs in the P-LIN model) were present around the lateral sulcus (LatS), in the medial occipital parietal cortex (mOP) and in the lateral occipital parietal cortex (lOP).

The static nonlinearity parameter, or ‘pRF exponent’ that models the nonlinearity of spatial summation (Materials and methods: Equation 4), was in the range of 0.2–0.4 and significantly below 1 in all areas (Wilcoxon signed-rank, one-tailed, all ROIs with more than four voxels R2 > 5%, p<0.001). These values of the pRF exponent represent subadditive (compressive) spatial summation, in accordance with observations in the human visual cortex (Kay et al., 2013; Winawer et al., 2013). The pRF exponent in early visual cortex is comparable to previously reported values for human V1. The exponent value in higher visual areas was similar to that in early visual cortex of the two monkeys, and higher than previously observed in human extrastriate cortex, which suggests that spatial suppression is less pronounced in higher areas of the monkey visual cortex than in higher areas of the human visual cortex.

Both the U-LIN model and the DoG model also performed better than the standard P-LIN model (Kruskal–Wallis, Tukey’s HSD, both ps<0.0001). The DoG had slightly better fits across all pooled voxels than the U-LIN model (Kruskal–Wallis, Tukey’s HSD, p<0.0001). The advantage of the DoG model over the P-LIN model was most pronounced in V1 and decreased in higher visual areas.

Negative pRFs from suppressed BOLD responses

There was a subset of voxels with negative BOLD responses for which both the U-LIN and DoG models provided much better pRF fits than the P-LIN and CSS models (arrows in Figure 4A, Figure 4—figure supplement 1C). We inspected the voxels for which the R2 in the U-LIN/DOG models was at least 5% higher than in the P-LIN model (gray, dashed triangles in Figure 4A). The U-LIN model estimated a negative gain for these voxels (median gain = –0.31, Wilcoxon signed-rank, one-tailed, z = –43.9, p<0.0001) (Figure 4B) and the DoG model returned a high level of inhibition (median normalized suppressive amplitude = 1.14, interquartile range [IQR], 0.98–1.29).

There were two categories of voxels with negative responses (Figure 4—figure supplement 2). For the first category of voxels, the P-LIN model estimated pRFs outside the boundaries of the stimulated visual field. This result suggests that the negative response represents surround suppression that is particularly strong around the fovea (Sereno et al., 1995; Shmuel et al., 2006; Smith et al., 2004). The retinotopy of these voxels is consistent with this explanation. In V1, for instance, the voxels were on the medial side of the occipital pole, which represents the peripheral visual field (Figure 4C). The second category of voxels with negative pRFs appeared to be different. Here, P-LIN and CSS models could not fit any pRF, suggesting purely negative BOLD responses. These voxels were primarily located in the medial occipital parietal cortex (mOP), at the superior border of the superior temporal sulcus in the lateral occipital parietal cortex (lOP) and around the lateral sulcus (LatS), which includes parts of the insula, cingulate, parietal, and premotor cortices (Figure 4D). These areas have all previously been identified as being part of the monkey DMN (Mantini et al., 2011). This finding aligns with recent research in humans that revealed similar negatively tuned pRFs in corresponding nodes of the human DMN (Szinte and Knapen, 2020).

pRF size as a function of eccentricity

As expected, pRF sizes were larger at higher eccentricities. This relationship was evident in all areas with larger numbers of well-fit voxels. pRFs were also larger in higher areas, which exhibited a steeper slope of the eccentricity-size relationship (Figure 5, Figure 5—figure supplements 1 and 2). The differences between the slopes in V1 and V2 are smaller than expected based on previous electrophysiological studies, but this is not uncommon with fMRI (Kay et al., 2013). In one animal, we also unexpectedly observed retinotopy in a number of higher areas, such as the anterior cingulate cortex (Figure 5—figure supplements 1 and 2). This brain region has been studied predominantly in the context of decision-making (Amiez et al., 2006; Fouragnan et al., 2019), but it does have resting state correlations with V1 (Griffis et al., 2017). Our design lacked the power for a more detailed investigation of this retinotopic organization, but this result may inspire future work on brain-wide retinotopic tuning.

Figure 5 with 2 supplements see all
Population receptive field (pRF) size as a function of eccentricity according to the compressive spatial summation (CSS) model.

(A) Eccentricity-size relationship for early and mid-level visual areas. (B) Eccentricity-size relationship for areas in the temporal and parietal lobes. (C) Eccentricity-size relationship for subcortical areas. (D) Eccentricity-size relationship for frontal cortical areas. Data points are pRF sizes binned in 2-dva-eccentricity bins; error bars denote SEM. Lines are linear fits with a significant slope (p<0.01). Slope values are shown between square brackets in the legend. The dashed line that separates the white and gray areas indicates the extent of the visual stimulus used to estimate pRFs. Data points in the gray regions come for pRFs that fell partially outside the region with a visual stimulus. The small panels show the same data for individual areas and include confidence intervals of the linear fit (shaded red area). Numbers denote the number of voxels per animal. The displayed areas were selected based on the presence of at least 20 voxels in each animal, with R2 > 5% for (A, B) and R2 > 3% for (C, D) due to the generally lower fit quality in the subcortex and frontal lobe. Figure 5—figure supplements 1 and 2 plot all suprathreshold voxels.

Multi-unit spiking activity RFs

We determined the RFs of MUA recorded with chronically implanted electrode arrays (Utah arrays) in areas V1 and V4 in two additional monkeys (M3 and M4) that did not participate in the fMRI experiments. We used a 1024-channel cortical implant, consisting of a titanium pedestal that was connected to 16 Utah arrays (Rousche and Normann, 1998). Each Utah array had 8 × 8 shanks with a length of 1.5 mm. In both monkeys, 14 arrays were placed in V1 and 2 in V4 of the left hemisphere (Figure 6). The stimulus was similar to that used in the fMRI experiments with some small differences due to constraints of the two setups (in the electrophysiology setup, the stimulus covered a larger portion of the visual field because the screen was closer to the animal) and the intrinsic nature of the recorded signals (stimulus steps were faster in the electrophysiology experiments because the electrophysiology signals are much faster than the BOLD response). We fit the four pRF models to the MUA responses and to the LFP power in five distinct frequency bands. We compared the MUA pRFs to a more conventional MUA RF-mapping method (cRFs). For this cRF method, we selected channels with an SNR larger than 3 (i.e., visual responses that were more than three times larger than the standard deviation of the spontaneous activity) and derived the RF borders from the onset and offset of the neuronal activity elicited by a smoothly moving light bar (see Materials and methods). Whenever both methods were able to estimate a pRF and cRF (R2 > 25% for the pRF method, SNR > 3 and R2 > 25% for the cRF method), the estimated locations were highly similar (median distance between pRF and cRF center, V1: 0.34, IQR 0.18–0.49 dva; V4: 0.90, IQR 0.40–1.41 dva). Compared to the P-LIN pRF model, the moving bar method estimated smaller cRFs (median size difference pRFsz-cRFsz: 0.50, IQR: 0.08–0.92) (Figure 7). The CSS model, however, returned pRF size estimates that were very similar to the cRF sizes or even a little bit smaller (median size difference pRFsz-cRFsz: –0.13; IQR: –0.41–0.14), suggesting that nonlinear spatial summation might indeed be better at capturing the RF properties of a small population of neurons than linear summation.

Figure 6 with 1 supplement see all
Visual field coverage of population receptive fields (pRFs) with Utah arrays.

(A) In both monkeys (M3, M4), 14 Utah arrays were implanted on the left operculum that is partly V1. Different colors represent the center of multi-unit activity (MUA)-based pRFs for the individual arrays. Only electrodes with R2 > 50% in the compressive spatial summation (CSS) model are shown. (B) Same as in (A), but for V4 electrodes. Note the different scale in the lower-right panel, with the gray arc indicating the extent of the visual stimulus. See Figure 6—figure supplement 1 for pRF sizes.

Comparison of multi-unit activity (MUA) population receptive field (pRF) sizes with conventionally determined RF (cRF) sizes (moving bar stimulus).

Data points represent recording sites of individual animals (black: M3; blue: M4) and brain areas (closed circles: V1; open circles: V4). (A) pRF sizes estimated with the P-LIN model (X-axis in left panel) are larger than cRF sizes obtained with a thin moving luminance bar (Y-axis in left panel). The median difference between pRF and cRF sizes across all electrodes (pooled across animals and areas) was 0.50 (interquartile range [IQR]: 0.08–0.92) and the median ratio was 1.41 (IQR: 0.99–1.83), as shown in the top and bottom-right panels, respectively. (B) As in (A), but for pRF sizes estimated with the compressive spatial summation (CSS) model, which are slightly smaller than the cRFs (median difference: –0.13, IQR: –0.41–0.14; median ratio 0.90, IQR: 0.67–1.12).

The pRF models that were used in this study are all based on circular (symmetric) RFs. A recent study (Silson et al., 2018) suggested that pRFs in human early visual cortex might be elliptical rather than circular, although this suggestion goes against previous work (Greene et al., 2014; Merkel et al., 2018; Zeidman et al., 2018). A later study demonstrated that the elliptical fits were an artifact of the software that had been used in the analysis (Lerma-Usabiaga et al., 2021). The cRF method separately estimates the width and height of the RF, and can thus be used to calculate a simplified RF aspect ratio to investigate RF symmetry. While this measure differs from RF ellipticity (a 45° tilted ellipse has the same aspect ratio as a circle), it does provide some insight into the symmetry of the MUA RFs. We did observe a few cRFs with aspect ratios (σlargesmall) that were larger than 2 (M3: 18/753; M4: 10/527; together 2.2% of all cRFS), but the vast majority of cRFs in both animals had aspect ratios close to 1 (M3 median: 1.12, IQR: 1.04–1.20; M4 median 1.13, IQR: 1.03–1.22) indicating near-symmetric RFs.

We obtained excellent fits to the MUA for all pRF models (see Figure 1E for an example fit). These pRFs covered a large proportion of the lower-right visual field (Figure 6, Figure 6—figure supplement 1). As expected, the pRFs from electrodes of the same arrays (shown in the same color in Figure 6) were clustered in space, and their locations were in accordance with established retinotopy (e.g., Hubel and Wiesel, 1974). The average R2 (over all electrodes with R2 > 0) was 64% in V1 (M3: 54%; M4: 73%) and 53% in V4 (M3: 37%; M4: 68%), which is substantially higher than the average R2 of 11% in both V1 and V4 for the MRI data (all voxels with R2 > 0; V1 M1: 14%, M2: 8%; V4 M1: 14%, M2: 8%).

Cross-validated comparisons revealed significant differences between the four models (Kruskal–Wallis test on all four models: HV1 = 204, dfV1 = 3, pV1 < 0.0001; HV4 = 13.4, dfV4 = 3, pV4 < 0.01) (Figure 8; similar patterns were present in each individual animal). Post-hoc pairwise comparisons (Tukey’s HSD) revealed that the CSS and DoG models provided a better fit than the linear models, although for V4 the advantage of the CSS model over the P-LIN model was only significant when electrodes with a poor fit (R2 < 25%) were excluded (V1, all electrodes: CSS vs. P-LIN, p<0.001, DoG vs. P-LIN, p<0.001; V4, all electrodes: CSS vs. P-LIN, p=0.16, DoG vs. P-LIN, p<0.001; V1, electrodes with R2 > 25%: CSS vs. P-LIN, p<0.001, DoG vs. P-LIN, p<0.001; V4, electrodes with R2 > 25%: CSS vs. P-LIN, p<0.02, DoG vs. P-LIN, p<0.01) (Figure 8). The improved fit of the DoG model was caused by the suppressive surround (median normalized suppressive amplitude = 0.71, IQR 0.56–0.86).

Figure 8 with 3 supplements see all
Comparison of multi-unit activity (MUA)-based fit results from the four population receptive field (pRF) models.

Scatterplots compare R2 of pRF models. Each dot represents an electrode (black: V1; green: V4).

The pRF exponent of the CSS model was 0.38 ± 0.23 in V1 (mean ± SD; 1358 electrodes with R2 > 25%), which was significantly smaller than 1 (Wilcoxon signed-rank, one-tailed, z = –31.59, p<0.0001) and similar to the MRI-based values, which had a mean of 0.34 ± 0.19 (6341 voxels). Likewise, in V4, MUA pRF exponent values were significantly smaller than 1 (0.33 ± 0.19; z = 11.13, p<0.0001; n = 165), and comparable to MRI-based values in the same area (0.30 ± 0.16; n = 2324). The similarity in the values of the pRF exponent indicates that subadditive spatial summation is a prominent feature in these areas.

The size of the estimated MUA pRFs increased with eccentricity (Figure 9). However, the pRF sizes for two of the V1 arrays in both monkeys were approximately three times smaller than expected when compared to the data from the other arrays. These outlying arrays were located on the posterior-medial side of the surface of V1 in both monkeys where the gray matter is relatively thin (Figure 9A). We therefore suspect that the 1.5-mm-long shanks of the Utah arrays were pushed into the white matter, where they picked up activity of thalamic afferents (i.e., the geniculostriate pathway). We therefore excluded these electrodes from the analyses of pRF sizes. The remaining MUA pRF sizes and eccentricities were highly similar to RFs reported in previous electrophysiology studies at similar eccentricities (Gattass et al., 1987; Gattass et al., 1981; Van Essen et al., 1984; Victor et al., 1994). Next, we compared pRFs between the electrophysiological signals and the fMRI-BOLD signal.

V1 arrays with outlying population receptive field (pRF) sizes.

(A) Schematic representation of the location of the craniotomy made during surgery (dashed line) and the implanted electrode arrays (rectangles) depicted on the NMT standard brain. The color map indicates the thickness of the cortical gray matter of the NMT. (B) For both monkeys, the estimated pRF sizes of V1 electrode arrays in the posterior medial corner of the craniotomy (red and orange data points correspond to red and orange rectangles in panel A) were surprisingly small for their eccentricity compared to the size-eccentricity relationship seen in the other arrays (gray circles, linear fit with 95% CI as black line and gray area). Given the length of the electrodes (1.5 mm), the typical thickness of the striate cortex, and these small pRF sizes, it is likely that the outlying pRFs do not reflect the tuning of V1 neurons, but that of the geniculostriate pathway in the white matter. The pRF sizes were estimated by the compressive spatial summation (CSS) model (recording sites shown have R2 > 70%).

Local field potential pRFs

The LFP was split into five frequency bands: θ (4–8 Hz), α (8–16 Hz), β (16–30 Hz), γlow (30–60 Hz), and γhigh (60–120 Hz). We fit the four models to estimate pRFs for each frequency band. The results for γlow and γhigh resembled those for the MUA with high R2 values for a large proportion of the electrodes (especially in V1). The CSS model again outperformed the other models, and we did not observe negative pRFs in these V1 and V4 regions (Figure 8—figure supplements 1 and 2). Electrodes with good MUA-pRF fits usually also had good LFP-pRFs, but the opposite was not always true (Figure 8—figure supplement 3). The pRFs in lower frequency bands differed. Whereas θ generally yielded low R2 values, α and β yielded good fits for a substantial number of electrodes. Interestingly, there were two classes of electrodes in these frequency bands. For the first class, the power increased with visual stimulation and CSS model fits were best. In contrast, the second class of electrodes had negative pRFs, that is, the stimulus suppressed power (Figure 8—figure supplement 1). Given the few electrodes with good low-frequency LFP-pRFs in V4, we focused our analysis on the positive and negative responses on V1 (the split into positive/negative pRF was based on the parameters of the U-LIN model; Figure 10, Figure 10—figure supplement 1).

Figure 10 with 2 supplements see all
Characteristics of local field potential (LFP)-α population receptive fields (pRFs) in V1 split by positive and negative gain values.

(A) Distribution of gain values for LFP-α pRFs of V1 electrodes estimated with the U-LIN model. Electrodes with positive gain pRFs are classified as α+, electrodes with negative gain pRFs as α-. (B, C) pRF eccentricities (B) and size (C) for α- and α+ electrodes (yellow shades). Colored boxes indicate interquartile range (IQR), and the median is shown as a thick horizontal line. Wilcoxon rank-sum tests were used for comparisons between α- and α+ electrodes. (D) Distance between the centers of LFP-α and LFP-γ pRFs from the same electrode, divided by the sum of their respective sizes. Values smaller than 1 indicate overlapping receptive fields. Nonshaded boxes are comparisons with LFP-γlow pRFs, shaded boxes are comparisons with LFP-γhigh pRFs. (E) Difference in eccentricity between LFP-γ and LFP-α pRFs from the same electrodes (calculated as Eccγ-Eccα). Positive values indicate that LFP-α pRFs are closer to fixation than LFP-γ pRFs. Wilcoxon signed-rank, one-tailed test (ΔEcc > 0), for individual cases; Wilcoxon rank-sum test for comparisons between α- and α+ electrodes. See Figure 10—figure supplement 2 for a visualization of the shifts per recording site. (F) Eccentricity-size relationship for α- (left) and α+ (right) electrodes. Dots indicate individual electrodes. Figure 10—figure supplement 1 shows the same pattern of results for the LFP-β pRFs.

We observed a number of remarkable differences between the positive and negative pRFs. First, negative pRFs generally had lower eccentricities than positive pRFs (Figure 10B, Figure 10—figure supplement 1B). Second, negative pRFs were larger than positive pRFs (Figure 10C, Figure 10—figure supplement 1C). Third, there was a systematic difference in the distance from the γ-pRF of the same electrode. Specifically, we calculated a ‘separation index’ by dividing the distance between the centers of the low- and high-frequency LFP-pRFs by their summed size estimates (SI = Distance/(σlf + σhf)). A separation index of less than 1 indicates pRF overlap. For both α and β power, the negative pRFs were farther from the γ-pRFs than the positive pRFs (Figure 10D, Figure 10—figure supplement 1D) and generally shifted in the direction of fixation (Figure 10E, Figure 10—figure supplements 1E and 2). These results suggest that the positive pRFs represent visually driven activity, whereas the negative pRFs signal a form of suppression that is strongest at smaller eccentricities, close to the fixation point. One possible explanation is that the monkeys directed attention to the fixation point, which may have caused a ring of suppression closely surrounding it. Another possibility is that the negative pRFs might be a consequence of small eye movements around the fixation point. The size of α-pRFs did not depend on eccentricity (Figure 10F), whereas the size of β-pRFs increased with eccentricity, although this relation was very weak for positive pRFs (Figure 10—figure supplement 1F).

The pRFs derived from all electrophysiology signals on the same electrode had similar locations, with less than 1 dva between their centers on average (CSS model, Figure 11A). We next analyzed RF sizes, normalizing the size estimates to the MUA-cRF. All LFP-pRFs estimates were larger than MUA-pRFs, and lower frequency LFP components yielded larger pRFs than the higher frequencies (Figure 11B; we observed the same pattern present in each animal). The pRF exponent was well below 1 for all LFP components (V1 and V4 electrodes with R2 > 25%; Wilcoxon signed-rank, one-tailed <1, p<0.001) and smaller at lower frequencies, indicative of stronger CSS. We also compared the exponents to those of the MRI-pRFs. The exponent of the MRI was not significantly different from that of γlow in V1, and both γlow and γhigh in V4 (p>0.05, Figure 11C). Differences between the MRI exponent and those of the other LFP bands were significant (V1: Kruskal–Wallis, H = 505.49, df = 6, p<0.0001; V4: Kruskal–Wallis, H = 21.65, df = 3, p<0.001; Tukey’s HSD for multiple comparisons).

Comparison of population receptive field (pRF) location and size of different electrophysiological signals at the same electrode, estimated by the compressive spatial summation (CSS) model.

(A) Median distance between receptive field (RF) estimates. Electrodes were only included if R2 > 25% (multi-unit activity [MUA], local field potential [LFP]) or signal-to-noise ratio (SNR) > 3 (MUA-cRF). (B) RF size for the electrodes of (A), normalized to the MUA-cRF (dashed line). Horizontal lines indicate the median, colored rectangles depict the interquartile range (IQR). (C) pRF exponent from the CSS model (indicating nonlinearity of spatial summation). The horizontal dashed line indicates linear summation. The exponent was significantly lower than 1 for all signals. Gray arrows indicate signals for which the exponent did not significantly differ from that of the MRI pRFs. A similar pattern was present in the electrophysiological data of each animal.

Comparison of pRF eccentricity-size relationship between fMRI and electrophysiology signals

We next compared the eccentricity-size relationship of the electrophysiological signals to that of the BOLD-fMRI pRFs using linear mixed models (LMMs) to evaluate the pRF estimates of the CSS model, separately for V1 and V4 (electrodes with R2 > 50% and R2 > 5% for MRI). We included voxels for which the pRF was in the lower-right visual quadrant where we had electrode coverage. In V1, positive eccentricity-size relationships existed for MRI, MUA, β, γlow, and γhigh (Figure 12A). We compared these signals with a single LMM, revealing an interaction between signal type and eccentricity (F = 14.44, df = 4, p<0.0001), which indicated that the slopes differed. We further used pairwise LMMs to compare the slope of the MRI-pRFs to the electrophysiological signals with a significant eccentricity-size relationship. The fMRI eccentricity-size relationship was similar to that of MUA (F = 0.38, df = 1, p=0.54), whereas it was significantly different from all LFP signals (β: F = 7.97, df = 1, p<0.01; γlow, F = 20.02, df = 1, p<0.001; γhigh, F = 20.05, df = 1, p<0.001) (Figure 12C). We repeated this analysis in V4 (Figure 12B) where we had fewer electrodes. We obtained positive eccentricity-size slopes for MUA, γlow, and γhigh but did not obtain fits of sufficient quality for the θ, α, and β bands. In contrast to V1, there was no clear difference across these signals (F = 1.11, df = 3, p=0.24). Pairwise comparisons with the MRI size-eccentricity relationship did not reveal a difference between those for γlow (F = 0.36, df = 1, p=0.55), γhigh (F = 0.28, df = 1, p=0.60), or MUA (F = 2.40, df = 1, p=0.12). Because the visual field coverage in V4 differed between the two animals that were used for electrophysiology, we repeated the analysis in both individuals and observed similar results (not shown).

Figure 12 with 1 supplement see all
Eccentricity-size relationship for population receptive fields (pRFs) across signal types.

(A) The pRF size-eccentricity relation for V1 electrodes (compressive spatial summation [CSS] model, R2 > 50%). Dots are individual electrodes, colored lines represent the slope of the eccentricity-size relationship. The dashed black line represents the relationship for V1 MRI voxels with a pRF in the lower-right visual quadrant and an R2 > 5%. We only included signals with >25 electrodes meeting the R2 threshold. (B) Same as in (A), but now for the V4 electrodes and voxels. (C) Eccentricity-size slopes (left: V1; right: V4). The dashed line represents the slope for MRI-based pRFs with the 95% confidence interval depicted in gray shading. Colored rectangles indicate the 95% confidence intervals for the electrophysiological signals and the horizontal black line the slope estimate. The lower bound of the 95% confidence interval of the LFP-θ is not visible. Asterisks indicate significant difference with the size-eccentricity relation of the MRI-based pRFs.

To further investigate the robustness of these results, we repeated the analysis (1) with the inclusion of either all V1 and V4 voxels or a subset of the voxels in approximately the same anatomical location as the electrode arrays and (2) by varying the R2-based inclusion criteria for electrodes and voxels (Figure 12—figure supplement 1). MUA-pRFs in V1 and V4 were generally similar to the BOLD-pRFs, although in V4 γlow and γhigh also approximated the fMRI results in some of the comparisons.

Discussion

The current study conducted a systematic comparison of pRFs obtained with fMRI and electrophysiological recordings in V1 and V4 of awake behaving macaque monkeys. Within the same species, we fit several pRF models to seven different signal types (BOLD, MUA, and the power in five LFP frequency bands) to gain insight into the neuronal basis of MRI-BOLD-pRF measurements (Wandell and Winawer, 2015). Our results demonstrate retinotopic tuning in many brain regions and the presence of negative pRFs in areas of the DMN. We found that subadditive spatial summation is a prominent feature of many measures of brain activity. Furthermore, the results establish a clear relationship between BOLD-pRFs and electrophysiologically determined pRFs, with MUA-pRFs being most similar to BOLD-pRFs (as discussed below).

Cortical and subcortical retinotopic tuning of the BOLD signal

Retinotopic information was present in occipital, temporal, parietal, and frontal cortex as well as in a subcortical areas. We could differentiate the LGN and pulvinar from their surrounding areas based on their higher R2. Retinotopic maps in the LGN were roughly in line with previously published retinotopic organization of the macaque LGN (Erwin et al., 1999), but they comprised few voxels and the variability across animals prohibited detailed inferences. In the pulvinar, the organization of the retinotopic information was even more variable across animals. In humans, subcortical retinotopic maps have been observed in LGN, pulvinar, superior colliculus (SC), thalamic reticular nucleus (TRN), and substantia nigra (SN) (Cotton and Smith, 2007; DeSimone et al., 2015; Schneider et al., 2004). Here, we also observed voxels with pRFs located lateral and medial to the pulvinar, but more targeted investigations are necessary for detailed individual-subject segmentation of the thalamus (DeSimone et al., 2015; Tani et al., 2011).

We obtained good pRF fits in the occipital, temporal, and parietal cortical areas, regions for which previous studies demonstrated retinotopic organization with phase-encoded retinotopic mapping (Arcaro et al., 2011; Janssens et al., 2014; Kolster et al., 2014; Kolster et al., 2010; Kolster et al., 2009; Patel et al., 2010). Despite the fact that checkerboard stimuli may not be ideal for the activation of frontal areas, which are better driven by more complex stimuli (Janssens et al., 2014; Saygin and Sereno, 2008), we also observed retinotopy in several frontal areas, including the insula, cingulate cortex, FEFs (area 8), orbitofrontal cortex, ventromedial prefrontal cortex, and dorsolateral prefrontal cortex. These areas are thought to be involved in visual processing and visual attention. Retinotopic maps have previously been reported in the cerebellum in humans (van Es et al., 2019). We did not detect such maps, likely because of a lower SNR in the cerebellum due to lower field strength, a less optimal coil placement, and the sphinx position of the monkeys.

Negative pRFs

We observed negative pRFs in the MRI data and in the low-frequency LFP, but they were of a different nature. We identified two classes of negative pRFs in the MRI data. Negative pRFs in the visual cortex were usually accompanied by positive responses at peripheral visual field locations, close the boundaries of the visual display. Similar negative visual BOLD responses have been reported in human visual cortex (Smith et al., 2004) and are presumably caused by surround suppression (Allman et al., 1985; Cavanaugh et al., 2002; Hubel and Wiesel, 1962; Knierim and van Essen, 1992). We found a second class of negative pRFs around the LatS, in the mOP, and at the superior border of the superior temporal sulcus in the IOP, all areas previously implicated in the monkey’s DMN (Mantini et al., 2011). A similar retinotopy of negative responses was recently found in the human DMN (Szinte and Knapen, 2020). The retinotopy of nodes of the DMN in humans and monkeys implies a sensory-based organization of this part of the DMN, which could play a role in visual cognition (Arsenault et al., 2018).

Negative pRFs were also observed for a subset of the electrodes in the α and β range of the LFP, although pRF fits were generally of lower quality than those for the γ-LFP or MUA signals. It is conceivable that these negative pRFS reflect shifts of attention, which are known to modulate α and β power (Griffiths et al., 2019; Siegel et al., 2008; Womelsdorf and Fries, 2007; Worden et al., 2000). In accordance with this idea, negative α/β pRFs were larger than their positive counterparts and shifted towards the fixation point. Future studies are needed for a closer examination of this phenomenon.

The relation between BOLD-based pRFs and electrophysiologically determined pRFs

We compared BOLD-based pRFs to electrophysiological measures in terms of their spatial summation characteristics, size, and the relationship between eccentricity and size (Amano et al., 2009; Dumoulin and Wandell, 2008; Felleman and Van Essen, 1987; Gattass et al., 2005; Kay et al., 2013; Larsson and Heeger, 2006; Van Essen et al., 1984; Victor et al., 1994). Eye movements can have an effect on the pRF size estimates, and it is therefore important to note that (1) we only included data from recordings where the animals maintained a high fixation performance (IQR-span of the horizontal and vertical eye position, M1: 0.23, 0.34 dva; M2: 0.36, 0.49 dva; M3: 0.18, 0.38 dva; M4: 0.37, 0.60 dva), and (2) we averaged across multiple stimulus presentations to obtain robust response profiles. Hence, variations in eye position can have had only minor effects on the present results.

Compressive spatial summation

The CSS pRF model includes nonlinear spatial summation (Kay et al., 2013). Previous studies observed subadditive spatial summation, also known as spatial compression, throughout the human visual cortex in BOLD-based pRFs, with stronger compression in higher visual areas than in V1 (Kay et al., 2013; Winawer et al., 2013). In the present study, spatial compression was present in subcortical and cortical areas with good BOLD-pRF fits, indicating that it is a widespread, if not universal, characteristic of the primate visual system. The strength of compression did not differ much across areas, and it was similar to that in human V1. For the MUA data, the CSS model estimated pRF sizes that were very similar to RF size estimates derived from conventional methods that use moving luminance bars, whereas the P-LIN model systematically returned larger pRF estimates. This difference might indicate that spatial compression indeed better captures the neuronal RF properties, at least in V1 and V4.

In human iEEG recordings, spatial compression in the broadband iEEG signal was reported to be similar in strength to that of the BOLD signal (Winawer et al., 2013). We also observed CSS in both the MUA-pRFs and in all frequency bands of the LFP, and the pRF exponent tended to be smaller for lower-frequency components, indicative of stronger compression. The lower frequencies also exhibited stronger suppression than the BOLD signal. In contrast, the pRF exponent of MUA and γlow was similar to that of BOLD-fMRI in both V1 and V4, but in V1 the exponent for γhigh was larger. Hence, MUA and γlow were the two electrophysiological markers for which spatial compression resembled that of BOLD-fMRI most.

pRF size and the eccentricity-size relationship

As expected, pRFs were larger at higher eccentricities, with a larger slope in higher visual cortical areas. The pRF sizes and the eccentricity-size relationships that we found in visual cortex were in line with previous results from human and monkey neuroimaging studies (Amano et al., 2009; DeSimone et al., 2015; Dumoulin and Wandell, 2008; Kay et al., 2013; Keliris et al., 2019; Kolster et al., 2014; Kolster et al., 2010; Welbourne et al., 2018; Zhu and Vanduffel, 2019; Zuiderbaan et al., 2012). Furthermore, the results are in keeping with previous electrophysiological recordings of single units and MUA in early visual areas such as V1 (Gattass et al., 1987; Gattass et al., 1981; Van Essen et al., 1984; Victor et al., 1994), V2 (Burkhalter and Van Essen, 1986; Gattass et al., 1981; Rosa et al., 1988), V3 (Burkhalter and Van Essen, 1986; Felleman and Van Essen, 1987; Newsome et al., 1986; Rosa et al., 2000), and V4 (Gattass et al., 1988). Finally, they also match the results of human electrophysiology with subdural electrodes (Harvey et al., 2013; Yoshor et al., 2007).

The spatial scope of fMRI-BOLD, MUA, and LFP

The spatial scale of recorded neural signals depends both on the nature of the signal and the recording method. For extracellular electrophysiology recordings, the impedance and size of an electrode determine its sensitivity for single-unit spiking activity or MUA. This relationship is complex (Viswam et al., 2019; Ward et al., 2009), but high-impedance (~1 MΩ) electrodes are generally better suited for the detection of spiking activity, with smaller contact sites and higher impedance thought to sample from smaller neural populations and being more likely to pick up single-neuron activity. The spatial range of the MUA signal is ~140 µm (Buzsáki and Draguhn, 2004). Electrode impedance, size, and shape are less important for LFP recordings, at least within commonly used ranges (Nelson and Pouget, 2010; Viswam et al., 2019). The spatial extent of the LFP signal is a topic of ongoing debate (Kajikawa and Schroeder, 2011) with some authors estimating it as low as 120–250 µm in visual cortex (Katzner et al., 2009; Xing et al., 2009), while others suggest it may stretch up to several millimeters (Berens et al., 2008; Katzner et al., 2009; Kreiman et al., 2006), and may even be detected more than a centimeter away from the source (Kajikawa and Schroeder, 2011). The hypothesis that low-frequency components of the LFP have a larger spatial reach than high-frequency components is supported by computational modeling (Leski et al., 2013), although it has also been challenged (Dubey and Ray, 2016). Yet, the exact origin of the LFP and its relationship to spiking remain to be completely understood (Einevoll et al., 2013).

Given the larger spatial spread of the LFP compared to MUA, one would expect pRFs based on LFPs to be larger than those based on MUA. This is indeed what we found. The size of LFP-pRFs furthermore depended on the frequency component of the LFP, with lower frequencies yielding larger pRFs, suggesting a frequency-dependent spread of visual information in the LFP signal. In fMRI, there is a trade-off between spatial and temporal resolution. While higher magnetic field strengths and specialized acquisition methods are continuously increasing spatiotemporal resolution, studies with voxel sizes on the order of 1–2 mm isotropic and a repetition time of 2–3 s at 3T field strengths are still common. Due to the complex relationship between the hemodynamic BOLD signal and the underlying neural activity, it is difficult to predict its stimulus sensitivity from neuroimaging parameters such as voxel size. Nonetheless, a comparison of the eccentricity-size relationship between MRI and electrophysiology in V1 revealed that the BOLD-based pRFs resembled MUA in V1 and V4. In V4, however, the BOLD slope was also very similar to that of the LFP gamma power (Figure 12—figure supplement 1). One possible reason for this difference between areas is that V4 is smaller and more heterogeneous (Kolster et al., 2014; Zhu and Vanduffel, 2019) than V1. A V4 voxel; therefore, samples from a neuronal population with more heterogeneous spatial tuning than an equally sized V1 voxel may therefore reflect the activity of a larger population of neurons, which is better approximated by gamma power. Other factors that might play a role are potential differences in the quality of the recorded signal across areas or partial volume effects that are likely to be more prevalent in V4 than in V1.

Conclusions

Our comparison of fMRI with large-scale neurophysiological recordings in visual cortex revealed that pRFs derived from the BOLD signal resemble MUA RFs. Subadditive spatial summation is a general feature of many brain areas and occurs for BOLD, MUA, and LFP. We observed negative pRFs in the monkey DMN and as part of center-surround organization of pRFs in early visual areas. The spatial compression and the eccentricity-size relationship of BOLD resembles that of MUA, but also bears a resemblance to gamma power. We conclude that BOLD-pRFs accurately represent the spatial tuning of the underlying neuronal populations.

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Biological sample(Macaca mulatta)Rhesus macaque (Macaca mulatta), maleBiomedical Primate
Research Center,
the Netherlands
n/a-
OtherPhilips Ingenia 3.0T MR systemPhilipsn/aAt Spinoza Centre for Neuroimaging, Amsterdam, the Netherlands
Other8-channel phased array receive MR coil systemKU Leuvenn/aCustom-built
Other16-channel MR pre-amplifierMR Coils BVn/aCustom-built
OtherETL-200ISCANRRID:SCR_021044MR-compatible eye tracker
OtherE3X-NHOmronn/aFiber optic amplifiers
Other5-RLD-E1 Liquid Reward SystemCrist Instrument
Company, Inc
n/aJuice reward system
OtherBOLDscreen 32 LCD for fMRICambridge
Research
Systems
n/aMR-compatible display
OtherUtah array (electrodes)Blackrock
Microsystems
n/a-
Other128-channel CerePlex M head-stagesBlackrock
Microsystems
n/aData acquisition
Other128-channel CerePlex M head-stagesBlackrock
Microsystems
n/aData acquisition
Other128-channel Digital HubBlackrock
Microsystems
n/aData acquisition
Other128-channel Neural Signal Processor (NSP)Blackrock
Microsystems
n/aData acquisition
Software, algorithmBlackrock Central Software SuiteBlackrock
Microsystems
n/a-
OtherET-49CTomas
Recording
n/aEye tracker
Software, algorithmMATLABMathWorksRRID:SCR_001622-
OtherLISA clusterSURFsaran/aComputing cluster
Software, algorithmdcm2niixhttps://github.com/rordenlab/dcm2niixRRID:SCR_01409-
Software, algorithmNipypehttp://nipy.org/nipype/RRID:SCR_002502Used as the basis of the custom NHP-BIDS pipeline
Software, algorithmNHP-BIDSNetherlands
Institute for
Neuroscience
RRID:SCR_021813In-house developed, available via: https://github.com/VisionandCognition/NHP-BIDS
Software, algorithmFreeSurferhttp://surfer.nmr.mgh.harvard.edu/RRID:SCR_001847Used as the basis of the custom NHP-Freesurfer
Software, algorithmNHP-FreesurferNetherlands
Institute for
Neuroscience
RRID:SCR_021814In-house developed, available via: https://github.com/VisionandCognition/NHP-Freesurfer
Software, algorithmPycortexhttps://gallantlab.github.io/pycortex/n/aUsed as the basis of the customized NHP-Pycortex
Software, algorithmNHP-PycortexNetherlands
Institute for
Neuroscience
RRID:SCR_021815In-house developed, available via: https://github.com/VisionandCognition/NHP-pycortex
Software, algorithmanalyzePRFhttps://kendrickkay.net/analyzePRF/n/aToolbox was edited for this study and made available with the code and data
Software, algorithmJupyter Notebookhttps://jupyter.org/RRID:SCR_018315-
Software, algorithmFSLhttp://www.fmrib.ox.ac.uk/fsl/RRID:SCR_002823-
Software, algorithmTracker-MRI: Experiment control softwareNetherlands
Institute for
Neuroscience
RRID:SCR_021816In-house developed, available via: https://github.com/VisionandCognition/Tracker-MRI
OtherNMT v1.3NIH, AFNI
https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/nonhuman/macaque/template_nmtv1.html#nmt-v1-3
n/aMacaque Brain Template and Atlas

Subject details

Request a detailed protocol

Four male macaques (M. mulatta; 7–12 kg, 5–8 years old) participated in this study. Animal care and experimental procedures were in accordance with the ILAR’s Guide for the Care and Use of Laboratory Animals, the European legislation (Directive 2010/63/EU), and approved by the Institutional Animal Care and Use Committee of the Royal Netherlands Academy of Arts and Sciences and the Central Authority for Scientific Procedures on Animals (CCD) in the Netherlands (license numbers AVD8010020173789 and AVD8010020171046). The animals were socially housed in an enriched specialized primate facility with natural daylight, controlled temperature and humidity, and fed with standard primate chow, supplemented with raisins, fresh fruits, and vegetables. Their access to fluid was controlled, according to a carefully designed regime for fluid uptake. During weekdays, the animals received diluted fruit juice in the experimental setup. We ensured that the animals drank sufficient fluid in the setup and received extra fluid after experimental sessions if needed. On the weekends, animals received at least 700 ml of water in the home cage. The animals were regularly checked by veterinary staff and animal caretakers, and their weight and general appearance were recorded in an electronic logbook on a daily basis during fluid-control periods.

Surgical procedures

Request a detailed protocol

Two animals (M1 and M2) participated in the MRI experiments and were implanted with an MRI-compatible plastic (PEEK) head-post, fixed to the skull with ceramic bone screws and acrylic (Papageorgiou et al., 2014; Vanduffel et al., 2001). Anesthetics, analgesics, and monitoring procedures were similar to previous surgical procedures in our laboratory and are described in detail elsewhere (Klink et al., 2017; Poort et al., 2012; Supèr and Roelfsema, 2005). Two other animals (M3 and M4) participated in the electrophysiology experiments. They were implanted with a custom 3D-printed titanium head-post that was designed in-house based on a CT scan of the skull, aligned to a T1-weighted anatomical MRI scan of the brain (Chen et al., 2020; Chen et al., 2017). The titanium head-post was attached to the skull with titanium bone screws, and the skin was closed around the implant without the use of any acrylic. In a second surgery, each animal was additionally implanted with a total of 1024 electrodes in 16 Utah electrode arrays (Blackrock Microsystems) in their visual cortices (14 arrays in V1, 2 arrays in V4; Figure 1B). Each array contained an 8-by-8 grid of 64 iridium oxide electrodes with a length of 1.5 mm spaced at a distance of 400 μm from each other. Pre-implantation electrode impedances ranged from 6 to 12 kΩ. A custom-designed 1024-channel pedestal was attached to the skull with titanium bone screws, and the skin was closed around it. More details on the surgical procedures have been published elsewhere (Chen et al., 2020; Chen et al., 2017).

Visual stimuli and procedures

Request a detailed protocol

In the MRI experiment, animals were head-fixed, sitting in the sphinx position (Papageorgiou et al., 2014; Vanduffel et al., 2001), and viewing a 32″ screen (1920 × 1080 pixels, 100 Hz) (Cambridge Research Systems) at the end of the bore, 130 cm away. pRFs were measured using conventional moving bar stimuli that traversed the screen in eight different directions behind a large virtual circular aperture (Figure 1). The borders of this virtual aperture were invisible because both the foreground and background had the same gray level (22.3 cd/m2) (Dumoulin and Wandell, 2008). In the MRI experiments, the bar sweep spanned 16° (diameter) in 20 steps (Dumoulin and Wandell, 2008). The moving bars were 2° wide and contained a checkerboard pattern (100% contrast; 0.5° checkers; luminance of white checkers: 106.8 cd/m2; luminance of black checkers: 0.2 cd/m2) that moved parallel to the bar’s orientation. Each bar position was on the screen for 2.5 s (1 TR), making one full bar sweep last 50 s. Bar sweep series (all directions presented once) were preceded and followed by 37.5 s (15 TRs) of neutral gray background. Each horizontal or vertical bar sweep was followed by a neutral gray background period of 25 s. The order of the bar sweep directions was 270°, 315°, 180°, 225°, 90°, 135°, 0°, 45° on most runs, but for one animal we inverted the directions to 90°, 135°, 0°, 45°, 270°, 315°, 180°, 225° on some runs to compensate for the animal’s tendency to fall asleep near the end of runs. We included data from 8 scanning sessions for monkey M1 (34 runs, 268 bar sweeps) and 10 sessions for monkey M2 (66 runs, 406 bar sweeps). During stimulus presentation, the animals received fluid rewards (Crist Instruments, Hagerstown, MD) for maintaining fixation within a circular fixation window with a diameter of 2°, centered on a 0.15° red fixation dot, surrounded by a 0.75° (diameter) aperture of neutral gray background color. In the electrophysiology experiments, the stimulus and task were very similar, but bar sweeps now spanned 28°, which was possible because the animals were closer to the monitor (luminance values of this monitor were black: 0 cd/m2; white: 92.1 cd/m2, neutral gray: 14.8 cd/m2). Bars traveled along this path in 30 steps of 500 ms, and the neutral gray luminance intervals were reduced to 2.5 s due to the much faster neuronal responses (compared to the BOLD signal). The fixation window in the electrophysiology experiments was slightly smaller with a diameter of 1.5°. In the MRI experiment, eye position and pupil diameter were tracked with an MRI-compatible infrared eye-tracking system at 120 Hz (ISCAN ETL-200). Hand positions were also monitored using fiber optic amplifiers (Omron E3X-NH) and optic fibers. To reduce body movement-related imaging artifacts, the animals were trained to maintain their hands inside a response box by making reward delivery contingent on both eye and hand position. In the electrophysiology experiments, animals were head-fixed in a conventional vertical primate chair and viewed a 21″ CRT monitor (1024 × 768, 85 Hz) at a distance of 64 cm while their eye position and pupil diameters were tracked at 230 Hz using an infrared eye tracker (TREC ET-49B, Thomas Recording GmbH).

MRI acquisition

Request a detailed protocol

MRI was performed in a standard Philips Ingenia 3.0T horizontal bore full-body scanner (Spinoza Center for Neuroimaging, Amsterdam, the Netherlands). We used a custom-built eight-channel phased array receive coil system (Ekstrom et al., 2008; Kolster et al., 2009) (KU Leuven) and the scanner’s full-body transmit coil. Functional images were obtained using a gradient-echo T2* echo-planar sequence (44 horizontal slices, in-plane 72 × 68 matrix, TR = 2500 ms, TE = 20 ms, flip angle = 77.2°, 1.25 × 1.25 × 1.25 mm isotropic voxels, SENSE-factor of 2 in the AP direction, and phase-encoding in the AP direction).

fMRI preprocessing

Request a detailed protocol

All fMRI data were preprocessed with a custom-written Nipype pipeline that we have made available online (RRID:021813; https://github.com/visionandcognition/NHP-BIDS). In short, MRI scans were exported from the scanner as DICOM images and converted to NIFTI files with the dcm2niix tool (Li et al., 2016). The volumes were then reoriented to correct for the animal being in the sphinx position and resampled to 1 mm3 isotropic voxels. The resulting images were realigned using a nonrigid slice-by-slice registration algorithm based on AFNI tools (Cox, 1996) followed by an FSL-based motion correction procedure MCFLIRT (Jenkinson et al., 2002). Functional volumes were linearly aligned to the individual high-resolution anatomical volumes, which were in turn nonlinearly registered to the NMT standard space (Seidlitz et al., 2018). Preprocessed data were further processed with a combination of custom-written MATLAB (MathWorks, Natick, MA) and shell scripts (https://gin.g-node.org/ChrisKlink/NHP-PRF). BOLD time courses for each voxel were normalized to percentage signal change and averaged across runs (or parts of runs) for which fixation was maintained at 80% of the time or more. We averaged odd and even runs separately to allow for a cross-validation approach in the evaluation of the pRF model fits. Anatomical ROIs were defined based on a probabilistic atlas (Reveley et al., 2017; Seidlitz et al., 2018) and refined using individual retinotopic maps.

Post-fit comparisons across pRF models, HRFs, and ROIs were performed in MATLAB based on the volumetric results. For visualization of the fMRI data, volumetric results were also projected to the individual cortical surfaces. To create these surfaces, we averaged multiple anatomical scans (T1-weighted, 3D-FFE, TE = 6 ms, TR = 13 ms, TI = 900 ms, flip angle = 8°, 100 horizontal slices, in-plane 224 × 224 matrix, 0.6 × 0.6 × 0.6 mm isotropic voxels, and phase-encoding in the AP direction) and processed the result with customized tools based on FreeSurfer (Fischl, 2012) and Pycortex (Gao et al., 2015) that were adjusted to handle our NHP data. These tools and their documentation can be found at https://github.com/VisionandCognition/NHP-Freesurfer; (Klink, 2021; copy archived at swh:1:rev:8d4b89337b865fb194e196cf1b2af4967e14d607) and https://github.com/VisionandCognition/NHP-pycortex, respectively (Messinger et al., 2021).

Electrophysiology acquisition

Request a detailed protocol

Neuronal activity was acquired from 1024 channels simultaneously at a 30 kHz sampling rate. The 1024-channel pedestal was connected to eight 128-channel CerePlex M head-stages through an electronic interface board. Each head-stage processed signals from two 64-channel electrode arrays with a 0.3–7500 Hz analog filter at unity gain (i.e., no amplification). After analog-to-digital conversion, the signal from each head-stage was sent to a 128-channel Digital Hub (Blackrock Microsystems) where it was converted into an optical output signal and sent to a 128-channel Neural Signal Processor (NSP, Blackrock Microsystems) for storage and further processing. The eight NSPs were controlled with eight simultaneously running instances of the Blackrock Central Software Suite (Blackrock Microsystems) distributed over two computers (four instances each) (Chen et al., 2020).

Electrophysiology data preprocessing

Request a detailed protocol

The neuronal signal that was acquired using different software instances was temporally aligned using common TTL pulses sent by the stimulus computer. The data were then separated in (1) envelope MUA and (2) broadband LFP. MUA represents the spiking activity of a local population of neurons around the electrode (Cohen and Maunsell, 2009; Palmer et al., 2007; Supèr and Roelfsema, 2005). To extract MUA, we amplified the raw neuronal signal, band-pass filtered it between 500 Hz and 9 kHz, full-wave rectified it, and applied a low-pass filter of 200 Hz. The resulting time series were downsampled to 1 kHz. We subtracted the baseline MUA activity in a 1000 ms prestimulus time window. Baseline-corrected MUA responses were then averaged, first across runs and then within a 50–500 ms time window for each stimulus position. The broadband LFP signal was generated by low-pass filtering the raw signal at 150 Hz and downsampling it to 500 Hz. The LFP signal was further processed with a multi-taper method using the Chronux toolbox (Bokil et al., 2010). Power spectra were calculated in a 500 ms moving window (step size 50 ms), using a time bandwidth product of five and nine tapers. LFP power was averaged within five distinct frequency bands: 4–8 Hz (theta), 8–16 Hz (alpha), 16–30 Hz (beta), 30–60 Hz (low gamma), and 60–120 Hz (high gamma). Baseline power in a 1000 ms prestimulus period was subtracted and the power for each recording site was averaged across runs, within a 50–500 ms time window during each stimulus position.

pRF models and fitting procedure

Request a detailed protocol

We fit four pRF models to all the data (voxels and electrode channels) using a customized version of the analyzePRF toolbox (Kay et al., 2013) for MATLAB. In the fitting procedure, the stimuli were spatially downsampled to a resolution of 10 pixels per dva and converted to ‘effective stimuli,’ consisting of binary representations that encode stimulus position. Response predictions were calculated as the product of the effective stimulus and the pRF shape (Equation 1). This method also yields a prediction for pRFs that are only partially stimulated, and the best-fitting pRF can have a center location outside of the stimulus aperture. Because pRF estimates are generally more reliable for pRF that are strongly driven by a visual stimulus, we include indications of the directly stimulated visual field in our figures. The four pRF models differed in the pRF shape. Linear models (Equations 13) describe a single isotropic 2D Gaussian-shaped pRF and assume linear spatial summation across the visual field (Dumoulin and Wandell, 2008). We implemented two linear model versions. For the first model, responses were constrained to be positively related to the visual stimuli (P-LIN). A second version lacked this constraint and also allowed negative responses, that is, stimulus-driven activity reductions (U-LIN). Negative BOLD responses have been demonstrated in some brain areas (Shmuel et al., 2006; Szinte and Knapen, 2020). The nonlinear spatial summation model (Kay et al., 2013) expands the linear model by capturing nonlinear summation of signals across the visual field. It has previously been shown that the value of the exponent is generally smaller than 1 in human visual cortex, indicating subadditive spatial summation or CSS (Kay et al., 2013). This model is therefore generally referred to as the CSS model. Because the pRF size and static nonlinearity interact in the nonlinear model, the pRF size is defined as the standard deviation of the predicted Gaussian response profile to a point stimulus for all models (Equation 3; Kay et al., 2013). The mathematical descriptions of the linear and nonlinear pRF models are in Equations 13, where Resppred indicates the predicted response, g is a gain factor to scale the response, S(x,y) is the effective stimulus, G(x,y) is the Gaussian pRF profile, and n is the exponent that determines the static spatial nonlinearity. For the P-LIN model, the gain g was constrained to positive values, while for the U-LIN model, gain values could be negative as well. Negative gain factors imply stimulus-induced reductions of activity. In both linear models, the exponent n was fixed to be 1. In the definition of the Gaussian, (x0,y0) defines the center and σ the standard deviation of the pRF.

Linear and nonlinear spatial summation pRF models:

(1) Resppred=g[x,yS(x,y)G(x,y)]n,P-LIN,U-LIN:n=1P-LIN,CSS:g>0
(2) Gx,y=e-x-x02+y-y022σ2
(3) pRFsize=σn

The DoG model uses a standard linear Gaussian (G1 in Equation 4) to describe the excitatory center of a pRF and subtracts a second Gaussian profile (G2) to model an inhibitory surround component (Zuiderbaan et al., 2012). This second Gaussian is by definition broader than the one describing the center. The size σ2 and amplitude a of the surround Gaussian are additional parameters (Equations 46).

DoG pRF model:

(4) Resppred=gx,ySx,yG1x,y-aG2x,yn
(5) G1(x,y)=e(xx0)2+(yy0)22(σ1)2,G2(x,y)=e(xx0)2+(yy0)22(σ2)2;σ2>σ1

The models assume a near-instantaneous link between the stimulus and the response dynamics, which holds for the electrophysiological signals. For these signals, the stimulus changes position every 500 ms, much slower than the onset or decay of activity. The BOLD response is much slower than the speed with which the stimulus traverses the screen in the fMRI experiments (2500 ms per position). We therefore convolved the predicted response with an HRF at a resolution of 1.25 s per sample (twice the acquisition rate of TR = 2.5 s). We used both a canonical human HRF and a standard monkey HRF that we derived from separate scanning sessions (Figure 1—figure supplement 1). In short, we presented the animals with brief (0.1 s) full-contrast and full-screen checkerboard stimuli. We then used FMRIB’s Linear Optimal Basis Sets (FLOBS) (Woolrich et al., 2004) toolkit from the FSL software package to estimate the relative contributions of a set of basis functions for those voxels in the primary visual cortex that were activated by the stimulus. We then calculated a single weighted average HRF function based on these basis functions and used it as the standard monkey HRF. The monkey HRF was narrower than the canonical human HRF. It had a faster time-to-peak (4.2 s vs. 4.8 s) and peak-to-fall time (6.2 s vs. 12.6 s). The fits with the monkey HRF were slightly better, especially in the lower visual areas (Figure 1—figure supplement 1). However, size and location estimates were highly similar for the two HRFs for all models, and we only report results from model fits based on the monkey HRF.

Model fitting was performed on a cluster computer (LISA, SURFsara) using nonlinear optimization (MATLAB Optimization Toolbox). The accuracy of the different model fits was quantified as the cross-validated percentage of variance (R2) explained of the BOLD response (Equation 6: DATA) by the model prediction (Equation 6: MODEL). For cross-validation, we divided the data into two nonoverlapping sets (odd and even runs) and tested the prediction of a model that was fit one data set against the time series of the other data set and vice versa. This yielded two R2 values per voxel or electrode that were averaged. The cross-validated comparison of model performance is valid for models with different numbers of parameters and prevents overfitting. Fit results are available as voxel-based maps warped to the NMT template space on Neurovault.org (Fox et al., 2021; Gorgolewski et al., 2015) at https://identifiers.org/neurovault.collection:8082.

Model accuracy:

(6) pRFcenter-size=σ1,pRFsurr-size=σ2

Comparison of pRF (and cRF) estimates

Request a detailed protocol

After fitting the pRF models to all voxels and recording sites, we compared the pRF estimates both within and across recording modalities. We pooled voxels and recording sites across subjects and used nonparametric statistical tests (Wilcoxon signed-rank, Wilcoxon rank-sum, or Kruskal–Wallis) with post-hoc Tukey’s HSD tests to correct for multiple comparisons. For the fMRI data, we compared R2 across models and HRFs. We constructed retinotopic maps using the best pRF model and HRF and investigated the relationship between pRF eccentricity and size for a subset of ROIs with good model fits (R2 > 5%). For the electrophysiological data, we compared model accuracy and pRF estimates across MUA and LFP components to unravel to which extent retinotopic information is available in the different neuronal signals. We also compared the pRF estimates to a more conventional RF mapping technique for MUA based on responses to thin, moving bar stimuli (cRF). For recording sites with an SNR larger than 3 (i.e., visual responses that were more than three times larger than the standard deviation of the spontaneous activity), we fitted a Gaussian to the averaged MUA traces and determined the onset and offset of the visual response as the mean of this Gaussian plus or minus its standard deviation (SD). Horizontal and vertical RF boundaries were then derived from the onset and offset times for stimuli moving in opposite directions (Supèr and Roelfsema, 2005). RF centers were defined as the midpoint between the horizontal and vertical borders. For comparison with the pRFs, we calculated RF sizes as half the diagonal of the rectangular area between the horizontal and vertical cRF borders (cRFsz=width2+height2/2). This measure approximates the RF radius based on the SD of a Gaussian response profile and can be directly compared to the sigma estimated by pRF models (Figure 7). It is smaller than cRF sizes typically reported in electrophysiological studies because neurons are activated by stimuli farther than 1 SD from their RF center (our lab usually defines cRF diameter as 3.3 * SD).

To compare pRFs based on fMRI-BOLD and electrophysiology, we combined data from individual animals to create one pool of BOLD-based voxel pRFs and six pools of electrophysiology-based electrode pRFs (MUA, and the different frequency bands of the LFP) for V1 and V4 data. We compared the relationship between RF eccentricity and size with a set of LMMs. We first tested for a correlation between eccentricity and size for the different signal types. Signal types with a significant correlation were subsequently tested together in a single LMM to determine whether the eccentricity-size relationship differed between signal types (interaction SIGNAL × ECC). Finally, we compared electrophysiological signals with the MRI results. For this analysis, we only selected V1 and V4 voxels with pRFs falling within the eccentricity range of the electrode arrays and voxels and electrodes with a fit accuracy above a predetermined threshold (fMRI threshold: R2 > 5%, electrophysiology threshold: R2 > 50%). To test the robustness of the results for this analysis, we repeated it with different combinations of criteria for initial voxel inclusion (i.e., including only voxels that roughly correspond to the location of the electrode arrays, or including all V1 and V4 voxels) and fit quality (fMRI threshold: R2 > 5% or 10%, electrophysiology threshold: R2 > 25% or 50%; Figure 12—figure supplement 1).

Data availability

All data and code are available on GIN: https://doi.gin.g-node.org/10.12751/g-node.2j01af. Unthresholded fMRI model fitting results are available on Neurovault: https://identifiers.org/neurovault.collection:8082.

The following data sets were generated
    1. Klink PC
    2. Chen X
    3. Vanduffel W
    4. Roelfsema PR
    (2021) NeuroVault
    ID 8082. Visual population receptive fields - macaque fMRI.
    1. Klink PC
    2. Chen X
    3. Vanduffel W
    4. Roelfsema PR
    (2021) GIN
    Population receptive fields in non-human primates from whole-brain fMRI and large-scale neurophysiology in visual cortex.
    https://doi.org/10.12751/g-node.2j01af

References

    1. Brewer AA
    2. Press WA
    3. Logothetis NK
    4. Wandell BA
    (2002)
    Visual Areas in Macaque Cortex Measured Using Functional Magnetic Resonance Imaging
    The Journal of Neuroscience 22:10416–10426.
    1. Burkhalter A
    2. Van Essen DC
    (1986)
    Processing of color, form and disparity information in visual areas VP and V2 of ventral extrastriate cortex in the macaque monkey
    The Journal of Neuroscience 6:2327–2351.

Decision letter

  1. Kristine Krug
    Reviewing Editor; University of Oxford, United Kingdom
  2. Timothy E Behrens
    Senior Editor; University of Oxford, United Kingdom
  3. Holly Bridge
    Reviewer; University of Oxford, United Kingdom

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Decision letter after peer review:

Thank you for submitting your article "A comparison of population receptive fields from fMRI and large-scale neurophysiology recordings in non-human primates" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Timothy Behrens as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Holly Bridge (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) Framing of the study:

– Title and Abstract should differentiate more clearly that the direct comparison between neuronal activity and MRI pRFs was done only for V1 and V4 only. The wider MRI dataset was used for assessment of pRF models. The abstract felt somewhat misleading because it suggests there is neurophysiological data from areas other than V1 and V4. It should be made clear from the outset that neurophysiological data are only available in these areas. The statement 'We found that pRFs derived from BOLD-fMRI were most similar to MUA-pRFs in areas V1 and V4…..' makes it sounds like they were compared in other areas, but were most similar in those ones.

– In terms of improvement in the manuscript, we had difficulty determining the precise hypothesis or aim of the manuscript because there seem to be 3 different ideas (1) the relationship of pRF with eccentricity across many cortical and subcortical areas (2) 4 different models for fitting the pRF model and (3) the comparison of BOLD pRF with MUA and LFP in areas V1 and V4. It would therefore be helpful to understand better from the manuscript what the authors intend to be the major question.

2) In some cases, pRFs from higher areas in the frontal lobe etc appear to be coming from just one animal and the number of voxels meeting the R2 criterion is very low. We would suggest that analyses are limited to regions where pRF can be fit in both animals (and with a threshold for the minimum number of voxels). The entry criteria used should be reported in the Results.

3) Figure 3. LGN and Pulvinar:

a) The legend states that both areas have retinotopic maps of the contralateral visual field. We would definitely agree with this assertion for the LGN, but this does not appear to be the case for the pulvinar. Is this quantified somewhere?

b) Figures 3A and 3B suggest that M1 and M2 had inverted retinotopic maps in the LGN in one hemisphere between the two animals. M1 shows upper visual field (yellow) dorsally and M2 shows upper visual field ventrally (yellow; sagittal view). The LGN retinotopic maps do not look consistent between the two monkeys, which suggests that they might not as reliable as stated in the results.

4) Figure 5 really needs some data points rather than just the fits. We realise all the voxels are shown in the supplementary materials, but showing binned data with errorbars would be very helpful. It would also be more transparent, given the considerable differences in voxel numbers (mentioned above), to have the saturation of the data and fits reflect the voxel counts. That would also remove the issue with the lines intersecting the y-axis at or below zero as this is difficult to explain.

5) For the neurophysiology and MRI comparisons, we wonder whether the authors could provide a more direct comparison of the data considering only the pRF data from the regions contributing to the neurophysiological recordings.

6) The discussion of the pulvinar seems out of place in this manuscript as it does not seem to be the point of the paper and there is little analysis of the data from this area (noted above that representation does not appear to be consistently contralateral).

7) Some critical details are missing in the methods with regards to the assessed visual field and the stimulus used:

a) In the results and in methods, the covered visual field is described as 16 degrees (MRI) and 28 degrees (neurophysiology). Since this critically limits the RF size-eccentricity relationship that can be assessed, the authors should clearly state whether this is the diameter or radius of the visual field and whether (and how) horizontal and vertical dimensions of the available visual field differed. We assume that the stimulated visual field is only up to 8 degrees eccentricity for the MRI data – this needs to be clearly stated and in every data figure with actual data points and fits/extrapolations distinguished.

b) If partially stimulated receptive fields are extrapolated, then these need also to be indicated where applicable in the results, their generation explained in the methods, and the implications discussed, as the precise shape of these depends on the assumptions made.

c) Also, the results (line 130) refer to a "large circular aperture" behind which the stimulus was shown and the results refer here to Figure 1. The stimulus in Figure 1 shows no such aperture. Details of this (size, position, contrast) should be provided.

d) Mean luminance of the monitors used and the stimulus contrast should be provided for both set ups in the Methods.

8) Following on from 7), given the methods descriptions, it is not clear how some foveal or more eccentric RFs were assessed.

a) The results state clearly that the animals fixated the screen centrally. If the visual field refers to the diameter, this means that RFs more eccentric than 8 degrees (MRI) and 14 degrees (neurophysiology) could not be assessed.

This potentially affects the results in Figures 5 and 6B (M4). In Figure 5, the authors should distinguish between the range based on available data and extrapolation for each area.

b) Conversely, according to the methods around the fixation dot there was mask of mean background colour (aperture of 0.75 degrees – we assume radius?). But in Figure 6A, B and C, there are RF centres within 0.5 degrees of the fixation point. The authors should add an explanation to the methods how they assessed the RFs in these position, if they did not fully stimulate this part of the visual field.

9) Figure 7 suggests a systematic overestimation of pRF size over cRF? It would be good to compare both measures quantitatively to the neurophysiological literature, so the readers can understand what "the different definitions of RF size" refer to (p11, l246-247; p13, l281-283). The comparison as is in results and discussion remains vague as to RF definition and systematic differences across neurophysiology – MRI literature in estimated RF sizes.

10) Could the authors provide measures about the mean, range and variance in eye position for the MRI and neurophysiology experiments, as this would have potentially affected the size estimates of the smaller RFs. It would be good to discuss this issue.

11) Figure 6B suggests that the neurophysiology data for V4 samples two separate parts of the visual field in the two monkeys.

(a) Could the authors state in the paper whether there is in any of their analyses a significant difference in the data from the two monkeys and whether this could have affected the comparison with the MRI?

(b) Figure 6B suggests that in M4, by far the more eccentric visual field positions were sampled. However, Figure 7B suggests that most of the measured RF sizes are the same or smaller than in M3. Is this correct? And how does this relate to the literature on V4 RF sizes?

Reviewer #1:

This study provides a systematic assessment of visual population receptive fields (pRF) obtained with functional magnetic resonance imaging (MRI) in primates.

Critically, pRFs from early visual areas V1 and V4 are compared to systematic neurophysiological measurements obtained in individuals from the same primate species (but not the same individual animals) across a range of visual field positions. This allows at a direct comparison to "ground truth".

Furthermore, the paper provides an assessment of four different models for generating receptive fields from imaging and neurophysiological data, including linear and non-linear models, and investigates multi-unit (MUA) and local field potential (LFP) data. The authors show the greatest correspondence of pRFs between MUA and fMRI data, but the reasons for this could be more fully explored.

The paper could situate the results more clearly in the wealth of data we have about neurophysiological measures of receptive fields and retinotopy for the macaque visual system. The presented results could differentiate more clearly between data for which MRI and neurophysiological data are available (central visual field representations for areas V1 and V4), data for which extrapolation partially stimulated RFs was included and results that depend on MRI data only. Also, outside the early visual cortical areas the MRI retinotopic maps appear potentially less reliable and with this the pRF measurements.

These results are are equally important for research and clinical assessments of the neuronal basis of visual functional and perception.

Reviewer #2:

This paper by Klink and colleagues investigates the relationship between visual receptive fields defined with non-invasive magnetic resonance imaging (MRI) and neurophysiology in the non-human primate. Population receptive field mapping (pRF) is a non-invasive imaging technique designed to estimate the area of the visual field to which particular visual brain areas respond. Since these pRFs are derived from MRI data, they have provided a window into human visual brain function. In contrast, using invasive neurophysiology, it is possible to measure the receptive field from a group of nearby neurons (multi-unit activity; MUA) or the local field potential (LFP).

Data from four animals are presented, with pRF MRI data acquired from two animals and neurophysiological measurement of both MUA and LFP in the other two with some comparison between the measurements.

Strengths of the study:

1. The comprehensive analysis of the pRF mapping, which compares the results of four different models. The data are of very high quality due to the relatively high resolution of the imaging and the large quantity of data acquired. This allowed the identification of many areas in the occipital lobe, along with parietal, temporal and even some frontal regions. It was also possible to map pRFs in both LGN and pulvinar.

2. The finding that negative BOLD signals were due to two different sources – (i) being outside of the region that was visually stimulated and (ii) regions that are part of the default mode network – is an elegant finding and is consistent with previously published human data.

3. Recording both MUA and LFP power at different frequency bands provides a number of different measures of receptive field size that can be compared to the pRF data, and indeed the correspondence seems high.

4. Number of recording sites in V1 is good giving a good range of eccentricities over which to measure receptive fields.

Weaknesses of the study:

1. The authors are not clear about the main message of the paper. The title suggests that it is the comparison of the pRF and neurophysiological findings, but the neurophysiological recordings are (understandably) predominantly in V1, with some also in V4. This means that the majority of the pRF data cannot be compared to the neurophysiological recordings. Additionally, that the data were acquired in different animals is unfortunate as it would have been particularly valuable to measure the neurophysiological data from cortical tissue where pRFs have been measured non-invasively.

2. Four different methods for measuring pRFs were used, but the authors do not draw strong conclusions about which model may be the most successful to use in the discussion. Having this information would be very valuable for other researchers.

3. The discussion could be more focussed and does not really provide the reader with an answer to the question of how well the pRF compares to the neurophysiological data. Much of the data on pRF size and eccentricity has been published many times in the human and there is already some data in the non-human primate, so it is important to emphasise the novel findings of the current data. I would expect some more in-depth explanation of the MUA and LFP relationship with pRF and why MUA appears better.

4. Since data are from 2 animals, it would be helpful to see which animal each datapoint comes from to ensure that there are no major biases due to the contribution of each animal.

While the conclusions are supported to a certain extent, the key point that is critical to highlight is that the neurophysiological data are predominantly from V1. It is too general to state that MRI-based pRF measurements reliably reflect receptive field properties of the primate brain.

Reviewer #3:

The authors present fMRI data from 2 macaques and neural electrophysiology data from 2 different macaques implanted with 16*64 electrodes covering the lower right visual field quadrant of V1, and V4. Stimuli consisted of a traveling bar that systematically stimulated the retina while awake animals fixated centrally. pRFs fit to BOLD signals evoked by these stimuli, were compared with pRFs fit to time courses from multi unit activity (MUA) as well LPF theta, α, β, low γ and high γ frequency power.

This is an impressive study that offers a first systematic and very thorough comparison of pRFs estimated from FMRI versus from LFP as well as MUA electrophysiology measures. It thereby contributes to validation and interpretation of an increasing body of human pRF mapping research. The analyses are thorough and conclusions drawn are well supported by the data in my view (though note that I am not an expert on electrophysiological recording methods). A particular strength of the approach taken, is that the authors have fit 4 different pRF models to all these data using a robust cross-validation approach, which allows them to capture and compare various types of deactivation and spatial summation dynamics, and make inferences about parallels and differences in spatial information processing across brain regions and neural activity signatures. Most crucially, results convincingly reveal accuracy of fMRI pRFs with electrophysiology-recordings of MUA and γ LFP, in terms of similar relationships between pRF size, eccentricity, and the best-fitting pRF model parameters (a 'compressive spatial summation' model). Another strength includes the investigation of wide-spread retinotopy with fMRI, confirming recent findings of retinotopic selectivity in brain regions not classically considered visual or even sensory, such as the default mode network, in humans.

In terms of weaknesses, I see no major issues, although there are some aspects of the work that could be clarified. First, it would be helpful if the authors could elaborate on how pRF size estimated from MUA and LPF measures depend on the size or sensitivity of the electrodes to various neural signals, and what this means for expected comparability of pRF parameter estimates across these measures and to MRI. Second, there were a few points I felt lacked some detail needed to understand the results: one intriguing finding, is that for some electrodes in V1, pRFs can be derived from α and β frequencies in LFP signal that differ from those in the γ frequency. Interestingly, there are negative or 'deactivating' pRFs measurable in α and β LFP signals, which are substantially shifted in position with respect to 'activating' pRFs derived from the γ frequency LFP at the same electrode. The authors speculate this reflects foveally-directed shifts due to attention or eye-movement, each of likely relevance to studies mapping near-foveal regions. However, I found the relevant plots collapsing across electrodes (Figure 10) quite hard to interpret, and wondered whether differences in pRF fit across the two data types (Supplementary Figure 7) may contribute to these shifts. Similarly, in another very interesting analysis the authors show that pRF sizes fit to the MUA response evoked by the mapping stimulus, correlates well with RF size estimates obtained with the 'classic' approach of recording the MUA response to a thin moving light bar. Data reveal a clear linear relationship in RF size estimates across measures, but an analysis of expected vs. observed differences in size and shape (i.e., the true identity line) that would be useful for assessing methodological biases is not provided.

So, barring some minor points of clarification, I believe this work offers a wealth of data and analyses that provide important validation of commonly used pRF mapping methods, and will without doubt spark many further investigations into population receptive field dynamics across the primate and human brain. This will be hugely facilitated by the free availability of data and code for this study, made accessible via a link in the paper.

https://doi.org/10.7554/eLife.67304.sa1

Author response

Essential revisions:

1) Framing of the study:

– Title and Abstract should differentiate more clearly that the direct comparison between neuronal activity and MRI pRFs was done only for V1 and V4 only. The wider MRI dataset was used for assessment of pRF models. The abstract felt somewhat misleading because it suggests there is neurophysiological data from areas other than V1 and V4. It should be made clear from the outset that neurophysiological data are only available in these areas. The statement 'We found that pRFs derived from BOLD-fMRI were most similar to MUA-pRFs in areas V1 and V4…..' makes it sounds like they were compared in other areas, but were most similar in those ones.

– In terms of improvement in the manuscript, we had difficulty determining the precise hypothesis or aim of the manuscript because there seem to be 3 different ideas (1) the relationship of pRF with eccentricity across many cortical and subcortical areas (2) 4 different models for fitting the pRF model and (3) the comparison of BOLD pRF with MUA and LFP in areas V1 and V4. It would therefore be helpful to understand better from the manuscript what the authors intend to be the major question.

Thank you for pointing this out. Our main objective was to compare BOLD-pRFs and neurophysiology pRFs. However, because the data-set is rich, we also explored it by fitting different pRF models and looking at brain areas beyond V1 and V4 in the fMRI results. While we feel that all these aspects of the data yield interesting insights, we agree that this broadening of the scope may have obscured the main question a bit.

We changed the title of the revised manuscript to ‘Population receptive fields in non-human primates from whole-brain fMRI and large-scale neurophysiology in visual cortex’ and more clearly differentiate the main question of the BOLD-ephys comparison from the complementary findings in the revised abstract. The title is succinct as prescribed by eLife, while also explicitly mentioning the scope of neurophysiological and neuroimaging recordings. The abstract now also clarifies that direct comparisons were made only for visual areas V1 and V4. It now also makes clear that this comparison was our main question, but that additional insights were derived from the broader scope of the fMRI data. This distinction is now also properly specified in the introduction.

2) In some cases, pRFs from higher areas in the frontal lobe etc appear to be coming from just one animal and the number of voxels meeting the R2 criterion is very low. We would suggest that analyses are limited to regions where pRF can be fit in both animals (and with a threshold for the minimum number of voxels). The entry criteria used should be reported in the Results.

This assertion is correct. Selecting data with a stricter criterion primarily affects Figure 5, which we have now changed to address this and other reviewer suggestions. As suggested by the reviewers, we now only show areas for which at least 20 voxels reach the R2 threshold in both animals and explicitly mention how many voxels meet this threshold in each individual. This inclusion criterion is also explicitly mentioned in the caption of the new Figure 5. For more details see point 4.

3) Figure 3. LGN and Pulvinar:

a) The legend states that both areas have retinotopic maps of the contralateral visual field. We would definitely agree with this assertion for the LGN, but this does not appear to be the case for the pulvinar. Is this quantified somewhere?

We thank the reviewers for these observations and agree with them that the pulvinar results are not very consistent across animals. We have changed our phrasing in the legend to reflect this better. We also added a quantification of the proportion of voxels with significant contralateral pRFs (ncontra/ntotal). The new text in the legend reads:

“The lateral geniculate nucleus (LGN, top rows) contained retinotopic maps of the contralateral visual field in both monkeys (M1: 38/38; M2: 73/80 voxels with contralateral pRFs). Retinotopic information was also present in the pulvinar (PULV, bottom rows), but its organization was much less structured, especially in M2 (M1: 23/32; M2: 61/131 voxels with contralateral pRFs).”

b) Figures 3A and 3B suggest that M1 and M2 had inverted retinotopic maps in the LGN in one hemisphere between the two animals. M1 shows upper visual field (yellow) dorsally and M2 shows upper visual field ventrally (yellow; sagittal view). The LGN retinotopic maps do not look consistent between the two monkeys, which suggests that they might not as reliable as stated in the results.

We did not perform a detailed analysis of the LGN and pulvinar maps, as we considered the resolution of our data to be unsuitable for such inferences. When we compared our results to previously published retinotopic maps of the LGN (e.g. Erwin et al., 1999; Schneider et al., 2004; DeSimone et al., 2015), we confirm the representation of the upper visual field in more anterior-ventral LGN, and of the lower visual field in posterior-dorsal LGN. This is less clear in M2 where fewer voxels could be fit, and while a similar pattern might be present in some sagittal slices, others appear to show the opposite pattern. Whereas pRF model fits revealed retinotopic information for the pulvinar as well, these did not really constitute a clear ‘map’. We added a brief discussion about the reliability of the subcortical results to the Discussion (p19, ln 387-390):

“Retinotopic maps in the LGN were roughly in line with previously published retinotopic organization of the macaque LGN (Erwin et al., 1999), but they comprised few voxels and the variability across animals prohibited detailed inferences. In the pulvinar, the organization of the retinotopic information was even more variable across animals.”

4) Figure 5 really needs some data points rather than just the fits. We realise all the voxels are shown in the supplementary materials, but showing binned data with errorbars would be very helpful. It would also be more transparent, given the considerable differences in voxel numbers (mentioned above), to have the saturation of the data and fits reflect the voxel counts. That would also remove the issue with the lines intersecting the y-axis at or below zero as this is difficult to explain.

We have re-plotted the results in a revised version of Figure 5 to address both this point and point 2 above. Visualizing the number of voxels with saturation in addition to indicating the visual area with different colors did not work very well, but we found a way to represent the requested information. In revised Figure 5, we (1) only show areas for which both animals contribute at least 20 voxels, (2) mention the voxel-counts per area per animal, (3) plot binned data with error bars, (4) show linear fits for multiple areas in different colors in the same panel without confidence intervals, (5) show linear fits with confidence intervals for separate areas in smaller panels on the side, and (6) indicate the extent of the visual stimulus to clarify which pRFs were only partially stimulated (those in the gray areas). We also added these binned data-points to the corresponding Figure 5—figure supplement 1, which we expanded with a version of Figure 5—figure supplement 2 in which we lowered the R2 threshold to 3%. Since we make all data and scripts available online, interested readers will be able to investigate the effect of different voxel selection criteria on the results.

5) For the neurophysiology and MRI comparisons, we wonder whether the authors could provide a more direct comparison of the data considering only the pRF data from the regions contributing to the neurophysiological recordings.

Ideally, one would like to compare the BOLD-pRFs and ephys-pRFs from identical brain locations. However, because the recordings come from different individuals, the comparison depends on the precision with which we can estimate the correspondence of voxels in M1 and M2 with the location of electrodes in M3 and M4. We chose to approach this comparison from a functional angle: as the ephys-pRFs were located in the lower right visual field, we only included V1 and V4 voxels for which the BOLD-pRFs were also located in the lower right visual field. We have now expanded the analysis with the subset of V1 and V4 voxels in the left hemisphere that roughly corresponded to the location of the electrodes implanted in M3 and M4. Data from these three different voxel inclusion criteria were furthermore analyzed with four different combinations of R2-based thresholds for the MRI and electrophysiology results. As a consequence, the analysis was run on 12 different combinations of inclusion criteria, yielding similar results. This approach is described on p18 ln 369-373 and an overview of the results is presented in the revised Figure 12—figure supplement 1:

“To further investigate the robustness of these results, we repeated the analysis, (1) with the inclusion of either all V1 and V4 voxels or a subset of the voxels in approximately the same anatomical location as the electrode arrays, and (2) by varying the R2-based inclusion criteria for electrodes and voxels (Figure 12—figure supplement 1). MUA-pRFs in V1 and V4 were generally similar to the BOLD-pRFs, although in V4 γlow and γhigh also approximated the fMRI results in some of the comparisons.”

6) The discussion of the pulvinar seems out of place in this manuscript as it does not seem to be the point of the paper and there is little analysis of the data from this area (noted above that representation does not appear to be consistently contralateral).

We agree that the discussion of the (noisy) pulvinar could distract from the manuscript’s main points and have removed it. In addition, we extended the discussion of the subcortical results (as mentioned above in A3b).

7) Some critical details are missing in the methods with regards to the assessed visual field and the stimulus used:

a) In the results and in methods, the covered visual field is described as 16 degrees (MRI) and 28 degrees (neurophysiology). Since this critically limits the RF size-eccentricity relationship that can be assessed, the authors should clearly state whether this is the diameter or radius of the visual field and whether (and how) horizontal and vertical dimensions of the available visual field differed. We assume that the stimulated visual field is only up to 8 degrees eccentricity for the MRI data – this needs to be clearly stated and in every data figure with actual data points and fits/extrapolations distinguished.

Yes, the assumption is correct and we have clarified the methods description accordingly (p25 ln 542-546). We have also added an indication of the extent of the visual stimulus in Figures 5,6 and Figure 5—figure supplements 1,2.

b) If partially stimulated receptive fields are extrapolated, then these need also to be indicated where applicable in the results, their generation explained in the methods, and the implications discussed, as the precise shape of these depends on the assumptions made.

The fitting procedure is indeed capable of estimating partially stimulated pRFs and all reported pRF estimates are fully data-driven and not extrapolated. Especially for larger pRFs just outside the visually stimulated region this works well. We have added an explanation of how the method deals with partially stimulated pRFs to the description of the methods (p28 ln 632-635).

“This method also yields a prediction for pRFs that are only partially stimulated, and the best fitting pRF can have a center location outside of the stimulus aperture. Because pRF estimates are generally more reliable for pRF that are strongly driven by a visual stimulus, we include indications of the directly stimulated visual field in our figures.”

c) Also, the results (line 130) refer to a "large circular aperture" behind which the stimulus was shown and the results refer here to Figure 1. The stimulus in Figure 1 shows no such aperture. Details of this (size, position, contrast) should be provided.

We have improved the description. The foreground and background of this aperture had the same color, which makes the aperture ‘virtual’. It clipped the bar stimuli, but its edges were not visible on the screen. We have now added a dashed line to the revised Figure 1A to indicate the virtual aperture, removed the mentioning of the aperture from the main text of the Results section to avoid confusion, and described the stimulus better in the Methods section:

“Population receptive fields were measured using conventional moving bar stimuli that traversed the screen in eight different directions behind a large virtual circular aperture (Figure 1). The borders of this virtual aperture were invisible, because both the foreground and background had the same gray level (22.3 cd/m2) (Dumoulin and Wandell, 2008).”

d) Mean luminance of the monitors used and the stimulus contrast should be provided for both set ups in the Methods.

We have added this information to the Methods.

8) Following on from 7), given the methods descriptions, it is not clear how some foveal or more eccentric RFs were assessed.

a) The results state clearly that the animals fixated the screen centrally. If the visual field refers to the diameter, this means that RFs more eccentric than 8 degrees (MRI) and 14 degrees (neurophysiology) could not be assessed.

This potentially affects the results in Figures 5 and 6B (M4). In Figure 5, the authors should distinguish between the range based on available data and extrapolation for each area.

b) Conversely, according to the methods around the fixation dot there was mask of mean background colour (aperture of 0.75 degrees – we assume radius?). But in Figure 6A, B and C, there are RF centres within 0.5 degrees of the fixation point. The authors should add an explanation to the methods how they assessed the RFs in these position, if they did not fully stimulate this part of the visual field.

We have improved the methods description. The pRF fitting method is a forward modeling approach that predicts the responses to a visual stimulus for a hypothesized pRF (location and size) and minimizes the difference between the observed and predicted activity by changing the parameters of the hypothesized pRF. This means that (1) all estimated pRFs are data-driven, and (2) pRFs can be estimated even if their center is not stimulated. This holds for both the foveal (0.75 deg aperture was defined in diameter, which has been added to the Methods) and eccentric pRFs. Moreover, while we take the σ of the estimated pRF as a description of its size, the 2D Gaussians are larger than one σ, and a response occurs if a visual stimulus is farther from the center of the pRF. Figure 5 has been revised in light of other comments (see A4). In the new figure (as well as in other figures that deal with eccentricity), we added an indication of the extent of the visual stimulus to all figures that deal with eccentricity. We also added a few lines about cRF and pRF size definitions to the Methods section that will hopefully clarify this issue as well (p31 ln 718-722):

“This measure approximates the RF radius based on the SD of a Gaussian response profile, and can be directly compared to the σ estimated by pRF-models (Figure 7). It is smaller than cRF sizes typically reported in electrophysiological studies, because neurons are activated by stimuli farther than one SD from their RF center (our lab usually defines cRF diameter as 3.3*SD).”

9) Figure 7 suggests a systematic overestimation of pRF size over cRF? It would be good to compare both measures quantitatively to the neurophysiological literature, so the readers can understand what "the different definitions of RF size" refer to (p11, l246-247; p13, l281-283). The comparison as is in results and discussion remains vague as to RF definition and systematic differences across neurophysiology – MRI literature in estimated RF sizes.

This comment inspired us to recalculate the cRFs with a measure that is more comparable to the pRFs. In the previously submitted version of the manuscript, we took the position of a horizontally or vertically moving bar stimulus at visual field positions that included 95% of the area under a Gaussian fit to the MUA traces as the borders of the receptive field. Now we use the standard deviation of this Gaussian as a more direct analog to the pRF-size which is also quantified as 1 standard deviation. Interestingly, the cRFs calculated this way are smaller than pRFs estimated with the traditional P-LIN model (median difference 0.50 dva). This difference largely disappears when pRFs are calculated with the CSS model. Now, pRFs are slightly smaller than cRFs (median difference -0.13 dva). We have revised Figure 7 to show this direct comparison and changed the text where appropriate (p11 ln 247-251; p20 ln 439-443; p31 ln 714-722).

[p11 ln 247-251]”For this cRF method, we selected channels with a signal-to-noise ratio (SNR) larger than three (i.e., visual responses that were more than three times larger than the standard deviation of the spontaneous activity) and derived the RF borders from the onset and offset of the neuronal activity elicited by a smoothly moving light bar (see Materials and methods).”

[p20 ln 439-443]”For the MUA data, the CSS model estimated pRF sizes that were very similar to RF size estimates derived from conventional methods that use moving luminance bars, whereas the P-LIN model systematically returned larger pRF estimates. This difference might indicate that spatial compression indeed better captures the neuronal RF properties, at least in V1 and V4.”

[p31 ln 714-722]”Horizontal and vertical receptive field boundaries were then derived from the onset and offset times for stimuli moving in opposite directions (Supèr and Roelfsema, 2005). Receptive field centers were defined as the midpoint between the horizontal and vertical borders. For comparison with the pRFs, we calculated RF sizes as half the diagonal of the rectangular area between the horizontal and vertical cRF borders (cRFsz=(width2+height2)2). This measure approximates the RF radius based on the SD of a Gaussian response profile, and can be directly compared to the σ estimated by pRF-models (Figure 7). It is smaller than cRF sizes typically reported in electrophysiological studies, because neurons are activated by stimuli farther than one SD from their RF center (our lab usually defines cRF diameter as 3.3*SD).”

Various methods exist for the determination of RFs from spiking data, making comparisons across different labs somewhat difficult. The method we describe here is a procedure in use in our lab. To facilitate the comparison with the σ of pRFs we here use the σ value of a fitted Gaussian as measure for cRF size. However, we and other labs typically define the RF border at a multiple of this value (e.g., 1.65*σ) since neurons are activated by stimuli well beyond one σ from the center of their RF. In many of the papers we cite, RFs were manually marked based on subjective criteria derived from listening to neural activity (Gattass et al., 1987, 1981; Van Essen et al., 1994; Burkhalter and Van Essen, 1986; Rosa et al., 1988; Felleman and Van Essen, 1987; Newsome et al., 1986; Rosa et al., 2000; Gattass et al., 1988). Conversely, it is common in pRF studies to report the σ of the fitted Gaussian as measure of pRF size, but some studies also report pRF sizes as FWHM.

10) Could the authors provide measures about the mean, range and variance in eye position for the MRI and neurophysiology experiments, as this would have potentially affected the size estimates of the smaller RFs. It would be good to discuss this issue.

We used eye-trackers to monitor eye position and fixation performance determined reward delivery. We excluded epochs with insufficient fixation performance (<80%) from the subsequent model-fits and analysis. For all the recordings that were included in our pRF modeling pipeline, the monkeys had their gaze predominantly inside a fixation window that was slightly smaller in the electrophysiology recordings (0.75° radius) compared to the MRI experiments (1° radius) for practical reasons (vibrations from the scanner add a bit of noise to the eye position signal). We now determined the median eye position and the span of the IQR (in X and Y direction separately) for each animal across all sessions in both the MRI and electrophysiology experiments. The median eye positions for all animals was within 0.2° from the center of fixation, while the maximum IQR we observed was 0.6° for M4 in the electrophysiology recordings. Importantly, the IQR for eye position were similar for MRI and electrophysiology (MRI M1: IQRx = 0.2°, IQRy = 0.3°; M2: IQRx = 0.4°, IQRy = 0.5°; Electrophysiology M3: IQRx = 0.2°, IQRy = 0.4°; M4: IQRx = 0.4°, IQRy = 0.6°).

Brief and occasional eye-movements outside the fixation window and small residual eye-movements within the fixation window are unlikely to substantially affect the slow BOLD signal and BOLD-pRF estimates. Furthermore, they are likely to average out over stimulus set repetitions. Small residual eye-movements could cause a slight over-estimation of pRF-size. Moreover, the comparison of pRF size estimates across the various electrophysiology signals are unlikely to be influenced by eye movements, because the signals on which they are based were recorded at the same time (i.e., during the same eye-position variations).

Combining all of the above considerations, we consider it unlikely that eye movements play a significant role in our findings. In the discussion of the comparison of pRFs across signal types, we now briefly mention the potential effect of eye movements (p20, ln 427-432):

“Eye-movements can have an effect on the pRF size estimates and it is therefore important to note that (1) we only included data from recordings where the animals maintained a high fixation performance (IQR-span of the horizontal and vertical eye position, M1: 0.23, 0.34 dva; M2: 0.36, 0.49 dva; M3: 0.18, 0.38 dva; M4: 0.37, 0.60 dva), and (2) we averaged across multiple stimulus presentations to obtain robust response profiles. Hence, variations in eye position can have had only minor effects on the present results.”

11) Figure 6B suggests that the neurophysiology data for V4 samples two separate parts of the visual field in the two monkeys.

(a) Could the authors state in the paper whether there is in any of their analyses a significant difference in the data from the two monkeys and whether this could have affected the comparison with the MRI?

We ran the underlying analysis separately for M3 and M4 and found the same qualitative pattern of results in both animals. This is now also mentioned in the paper (p11 ln279, p16 ln345). In the comparison with fMRI, V4 differs from V1 since the eccentricity-size relationship of the γ LFP signals are sometimes not significantly different from the MRI results. We also ran the comparison with fMRI separately including either only the electrophysiology pRFs from M3 or M4. The results of these analyses were consistent with the original results based on the data from both animals (p18 ln 367-368).

“Because the visual field coverage in V4 differed between the two animals that were used for electrophysiology, we repeated the analysis in both individuals and observed similar results (not shown).”

(b) Figure 6B suggests that in M4, by far the more eccentric visual field positions were sampled. However, Figure 7B suggests that most of the measured RF sizes are the same or smaller than in M3. Is this correct? And how does this relate to the literature on V4 RF sizes?

We thank the reviewers for catching this inconsistency. The labels for M3 and M4 were flipped in the legend of Figure 7. We corrected this in the revised version of the figure that now makes a more direct comparison between cRF and pRF (see also point 9).

https://doi.org/10.7554/eLife.67304.sa2

Article and author information

Author details

  1. P Christiaan Klink

    1. Netherlands Institute for Neuroscience, Royal Netherlands Academy of Arts and Sciences, Amsterdam, Netherlands
    2. Psychiatry Department, Amsterdam UMC, Amsterdam, Netherlands
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Visualization, Writing – original draft, Writing – review and editing
    For correspondence
    c.klink@nin.knaw.nl
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6784-7842
  2. Xing Chen

    Netherlands Institute for Neuroscience, Royal Netherlands Academy of Arts and Sciences, Amsterdam, Netherlands
    Contribution
    Formal analysis, Funding acquisition, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3589-1750
  3. Wim Vanduffel

    1. Laboratory for Neuro- and Psychophysiology, Department of Neurosciences, KU Leuven Medical School, Leuven, Belgium
    2. Massachusetts General Hospital, Martinos Ctr. for Biomedical Imaging, Charlestown, United States
    3. Leuven Brain Institute, KU Leuven, Leuven, Belgium
    4. Harvard Medical School, Boston, United States
    Contribution
    Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
  4. Pieter R Roelfsema

    1. Netherlands Institute for Neuroscience, Royal Netherlands Academy of Arts and Sciences, Amsterdam, Netherlands
    2. Psychiatry Department, Amsterdam UMC, Amsterdam, Netherlands
    3. Department of Integrative Neurophysiology, Center for Neurogenomics and Cognitive Research, VU University, Amsterdam, Netherlands
    Contribution
    Conceptualization, Funding acquisition, Resources, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1625-0034

Funding

Nederlandse Organisatie voor Wetenschappelijk Onderzoek (VENI 451.13.023)

  • P Christiaan Klink

Nederlandse Organisatie voor Wetenschappelijk Onderzoek (STW-Perspectief P15-42 "NESTOR")

  • Xing Chen
  • Pieter Roelfsema

FP7 Ideas: European Research Council (ERC 339490 "Cortic_al_gorithms")

  • Pieter R Roelfsema

Human Brain Project (Agreements 720270 and 785907 "Human Brain Project SGA1 and SGA2")

  • Pieter R Roelfsema
  • Wim Vanduffel

Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Crossover Program 17619 "INTENSE")

  • Pieter R Roelfsema
  • P Christiaan Klink

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Jonathan Williford for his contributions to the fMRI preprocessing pipeline; Pieter Buur, Wietske van der Zwaag, Diederick Stoffers, and the Laboratory of Neuro- and Psychophysiology of KU Leuven for technical assistance in setting up the nonhuman primate MR infrastructure; Kor Brandsma, Anneke Ditewig, and Lex Beekman for animal care and biotechnical assistance; Feng Wang for help with electrophysiology data collection; Chris van der Togt for help with data management; and Tomas Knapen and Serge Dumoulin for fruitful discussion and comments on an earlier version of the manuscript. This work was supported by NWO (Crossover Program 17619 'INTENSE'; STW-Perspectief P15-42 'NESTOR'; VENI 451.13.023), the European Union FP7 (ERC 339490 'Cortic_al_gorithms'), the Human Brain Project (agreements 720270 and 785907, 'Human Brain Project SGA1 and SGA2'), and the Friends Foundation of the Netherlands Institute for Neuroscience.

Ethics

Animal care and experimental procedures were in accordance with the ILAR's Guide for the Care and Use of Laboratory Animals, the European legislation (Directive 2010/63/EU) and approved by the institutional animal care and use committee of the Royal Netherlands Academy of Arts and Sciences and the Central Authority for Scientific Procedures on Animals (CCD) in the Netherlands (License numbers AVD8010020173789 and AVD8010020171046).

Senior Editor

  1. Timothy E Behrens, University of Oxford, United Kingdom

Reviewing Editor

  1. Kristine Krug, University of Oxford, United Kingdom

Reviewer

  1. Holly Bridge, University of Oxford, United Kingdom

Publication history

  1. Preprint posted: September 8, 2020 (view preprint)
  2. Received: February 7, 2021
  3. Accepted: October 24, 2021
  4. Accepted Manuscript published: November 3, 2021 (version 1)
  5. Accepted Manuscript updated: November 4, 2021 (version 2)
  6. Version of Record published: December 3, 2021 (version 3)

Copyright

© 2021, Klink et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 842
    Page views
  • 164
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Neuroscience
    Blaise Robert et al.
    Research Article Updated

    Basal forebrain cholinergic neurons (BFCNs) project throughout the cortex to regulate arousal, stimulus salience, plasticity, and learning. Although often treated as a monolithic structure, the basal forebrain features distinct connectivity along its rostrocaudal axis that could impart regional differences in BFCN processing. Here, we performed simultaneous bulk calcium imaging from rostral and caudal BFCNs over a 1-month period of variable reinforcement learning in mice. BFCNs in both regions showed equivalently weak responses to unconditioned visual stimuli and anticipated rewards. Rostral BFCNs in the horizontal limb of the diagonal band were more responsive to reward omission, more accurately classified behavioral outcomes, and more closely tracked fluctuations in pupil-indexed global brain state. Caudal tail BFCNs in globus pallidus and substantia innominata were more responsive to unconditioned auditory stimuli, orofacial movements, aversive reinforcement, and showed robust associative plasticity for punishment-predicting cues. These results identify a functional topography that diversifies cholinergic modulatory signals broadcast to downstream brain regions.

    1. Neuroscience
    Rawan AlSubaie et al.
    Research Article Updated

    Projections from the basal amygdala (BA) to the ventral hippocampus (vH) are proposed to provide information about the rewarding or threatening nature of learned associations to support appropriate goal-directed and anxiety-like behaviour. Such behaviour occurs via the differential activity of multiple, parallel populations of pyramidal neurons in vH that project to distinct downstream targets, but the nature of BA input and how it connects with these populations is unclear. Using channelrhodopsin-2-assisted circuit mapping in mice, we show that BA input to vH consists of both excitatory and inhibitory projections. Excitatory input specifically targets BA- and nucleus accumbens-projecting vH neurons and avoids prefrontal cortex-projecting vH neurons, while inhibitory input preferentially targets BA-projecting neurons. Through this specific connectivity, BA inhibitory projections gate place-value associations by controlling the activity of nucleus accumbens-projecting vH neurons. Our results define a parallel excitatory and inhibitory projection from BA to vH that can support goal-directed behaviour.