Estimating average single-neuron visual receptive field sizes by fMRI.

The noninvasive estimation of neuronal receptive field (RF) properties in vivo allows a detailed understanding of brain organization as well as its plasticity by longitudinal following of potential changes. Visual RFs measured invasively by electrophysiology in animal models have traditionally provided a great extent of our current knowledge about the visual brain and its disorders. Voxel-based estimates of population RF (pRF) by functional magnetic resonance imaging (fMRI) in humans revolutionized the field and have been used extensively in numerous studies. However, current methods cannot estimate single-neuron RF sizes as they reflect large populations of neurons with individual RF scatter. Here, we introduce an approach to estimate RF size using spatial frequency selectivity to checkerboard patterns. This method allowed us to obtain noninvasive, average single-neuron RF estimates over a large portion of human early visual cortex. These estimates were significantly smaller compared with prior pRF methods. Furthermore, fMRI and electrophysiology experiments in nonhuman primates demonstrated an exceptionally good match, validating the approach.

An important contribution of functional magnetic resonance imaging (fMRI) in human neuroscience is the noninvasive in vivo estimation of the voxel-by-voxel organization of several cortical areas (1). Recent studies substantially advanced this field of research by using novel neurocomputational methods that can uncover neuronal properties previously only accessible by invasive electrophysiological techniques (2, 3). Estimating neuronal properties in vivo by fMRI is of great significance for understanding the functional organization of the cortex as well as cortical reorganization and plasticity in patients with diseases afflicting the brain (4).

A prime example of such methods is the estimation of population receptive fields (pRFs) in retinotopically organized visual areas (5–7). However, pRFs are only estimates of aggregate voxel-based averages of ten to hundreds of thousands of neurons within fMRI voxels and are a function of (i) the receptive field (RF) sizes of single units belonging to a voxel, (ii) the local scatter in the positions of RF centers, and (iii) the interactions [such as center-surround (8, 9)] between nearby connected units. Here, we present an approach that eliminates local scatter to estimate the average single-neuron RF sizes over a large portion of human early visual cortex. To achieve this, we exploited the spatial frequency (SF)-dependent fMRI responses of visual RFs modeled as Gabor functions and estimated single-unit RF (suRF) sizes from the properties (the SD of the Gaussian envelope) of the best-fitting Gabors. Thus, our estimates are not sensitive to RF scatter. Furthermore, we validated this approach by comparing noninvasive RF size estimates obtained using the same fMRI method in nonhuman primates with RF sizes obtained directly via intracranial electrophysiological recordings.

The relationship between RF scatter and the retinotopic organization of visual cortex has first been studied in primates by Hubel and Wiesel (10). These authors indicated that the random scatter of RF positions recorded from a vertical electrode penetration in the cortex is comparable to the RF size at a given eccentricity. However, following studies in cats that used better quantification indicated that random scatter may be smaller than previously believed (11). More recently, precise mapping of primate visual cortex with two-photon calcium imaging has investigated this further and found that random RF scatter is almost negligible (12). Instead, they observed that the positions of the RF centers of neurons within a small patch of cortex were smoothly changing depending on the local cortical magnification factor (10, 13, 14). This suggests that RF position is quasideterministic and can be potentially estimated directly by analytical formulations that have been developed to describe the global mapping from the visual field (retina) to the cortical space (15–17). Since typically the spatial sampling of fMRI voxels happens uniformly in cortical space (Fig. 1), the relationship of RFs within a voxel and the pRF can also be estimated if we know the RF sizes. To demonstrate that, we used the inverse of the simple model described by Schwartz (15), relating visual field to cortical location, to back-transform the coordinates of square pixels from cortical to visual space (Fig. 1 ). In addition, we used the linear relationship between single-neuron RF size and eccentricity () as reflected in the measurements of Hubel and Wiesel (10) to roughly estimate the expected RF sizes of neurons whose RF centers lie along the edges of the transformed fMRI-like pixels (Fig. 1, red circles). To this end, we assumed no additional random scatter given the recent findings by Nauhaus et al. (12). Dashed circles were then plotted encircling this group of RFs to represent each voxel’s expected pRF size. Note that expected pRF size strongly depends on the nonlinear inverse-Schwartz model transformation from cortical to visual space, clearly overestimating individual neuron RF size. Noteworthy, a similar relationship would also be present if we used other more recently suggested transformations that take into account more details of the known local magnification factor asymmetries (16, 17).

Fig. 1.

Visual field backprojection of fMRI voxels, RF modeling approach, and example of GLM fit of stimulus predictors. (A and B) Transformation of the fMRI voxels (squares) from the cortical space (A) to the visual field (B) by using the inverse-Schwartz model (15). In B, we plot estimated RF sizes along voxel borders at different eccentricities (red circles). Population RF size (pRF) depends not only on RF size but also on the nonlinear transformation from cortical to visual field. (C) We modeled suRFs as 2D-Gabor functions assuming that voxels contain an approximately homogenous representation of all orientations. Thus, voxel SF selectivity becomes independent of orientation and can be estimated according to the σ of the Gabor envelope. (D) Full-field checkerboard patterns were presented in blocks of 10-s ON and 20-s OFF. Blocks were pseudorandomly interleaved across six different SF conditions. (E) Raw BOLD-signal time course from a sample region of interest in area V1 of one subject (black) overlaid by the GLM fit of the six predictors (checker sizes) in different colors. (F) The GLM β-weights corresponding to each checker size predictor.

RFs of neurons in striate and early extrastriate visual cortex can be approximated by 2D Gabor functions (18, 19) whose SF selectivity can be analytically estimated (). We hypothesized that the amplitudes of the evoked blood oxygen level-dependent (BOLD) responses in voxels containing such neurons would thus depend on the SF content of the presented visual stimuli and that the recorded fMRI responses could be used to decode the average Gabor filter parameters of this population of cells including the SD σ of the Gabor envelope that is proportional to the RF size. To simplify this estimation, we took advantage of the fact that usual fMRI-voxel sizes are at least an order of magnitude larger than orientation columns, and assumed a homogeneous distribution of orientation selectivity within a voxel (Fig. 1 and ). Given this assumption, the Gabor SF selectivity becomes circularly symmetric and independent of orientation greatly reducing the Gabor parameters space to a single σ (approximately RF size) parameter (Fig. 1 and ). As stimuli, we chose to use binary white-noise checkerboard patterns and modulated their SF content by changing the checker size across visual stimulation blocks (Fig. 1). We then modeled the fMRI responses to these blocks using generalized linear modeling (GLM) to estimate the β-amplitudes for each condition and then used these β-values in voxel RF models to decode the average Gabor σ. We note that our Gabor-based modeling approach would also be valid if we used other stimuli with manipulated or selected SF content like sinewave gratings or filtered natural images. Here, we used checkerboard patterns since their frequency content could be estimated analytically, allowing easier computational modeling ().

To estimate the average suRF size in a voxel, we used the expected neuronal responses R for each checker size λ. These were calculated based on the integral of the product of the 2D SF content of the stimuli and the 2D SF sensitivity of the Gabor filter and given by the following:Furthermore, to fit the fMRI responses for each checker size , we added a compressive power-law static nonlinearity:where the exponent n < 1 acts as the compressive nonlinearity and γ is a single gain parameter across all conditions. We refer to this as the power-law compressive (PLC) model. This model has a single parameter σ that is proportional to RF size. In all RF estimations in this manuscript, we considered RF size to be equal to 2σ←←, as this is closer to the classical way of mapping RFs by detecting the first response to moving stimuli across the RF edges. demonstrates the effect of the compressive nonlinearity for different values of parameter σ.

Several studies have postulated that neural responses at a particular location are normalized by the responses of the surrounding neurons and that this normalization is a canonical neural computation (for a review, see ref. 20). In an attempt to also take into account responses of the surround we have used extended models that included additional gain and RF-size parameters for the surround. To this end, we either used a divisive normalization PLC model (dPLC) given by the following:or a subtractive normalization PLC model (sPLC) given by the following:where γ, σ, and γ, σ are gain and RF-size parameters for center and surround respectively. Since these models had more parameters than the simple PLC, we have set the exponent n = 0.325 based on previous research (21) to avoid overfitting. We note that, on average, this value was also very similar to the exponents fitted in the simple PLC model of Eq. . Simulations of the dPLC and sPLC for different values of σ and γ are presented in .

Results

Responses in human primary visual cortex were nicely modulated by our stimuli (Fig. 2). On average, V1 responded more strongly to blocks with smaller checkers (Fig. 2 ). More specifically, as expected from our hypothesis, voxels selected from anatomical and functionally defined areas close to the fovea were biased to the smallest checker sizes while gradually more peripheral voxels showed stronger responses to larger sizes (Fig. 2 ). Note that the amplitude of the response to each condition is reflected in the fitted GLM β-values that were then used for suRF modeling (Fig. 2 ).

Fig. 2.

Example of paradigm and V1 responses. (A) An example of the fMRI response in V1 of a subject is shown on top of a cartoon of the stimulus presentation with different conditions in colored blocks (10 s) and the interstimulus interval (20 s) in gray. (B and C) The event-related time courses for each condition are shown for left (B) and right V1 (C). (D and E) Similarly, the event-related time courses of selected single voxels in foveal (D) and peripheral (E) V1 are shown. (F and G) The GLM β-weights corresponding to each checker size predictor in D and E, respectively.

Examples of the PLC–suRF model fit for each subject H1–H4 are presented in Fig. 3 , respectively. On the top row, voxels with eccentricity closer to the fovea were selected based on the classical pRF model fit (ref. 5 and ); gradually higher eccentricities were selected in the middle and bottom rows. As expected from model simulations (), the shape of the model fit (solid lines) predicts small RF sizes for foveal voxels and gradually larger for more peripherally located ones. To better demonstrate the relationship between estimated voxel-based RF sizes and eccentricity, we performed linear regression for each subject (Fig. 3 , black lines). For comparison, linear regression was also performed for the pRF model across the same subjects (Fig. 3 , gray lines). The results demonstrated a significant linear relationship of suRF as well as pRF size with eccentricity for all individual subjects and across the population (). To test whether the suRF-estimated RF sizes were smaller than the pRF, as we hypothesized (Fig. 1), we performed analysis of covariance (ANCOVA) and second-level comparisons (Tukey–Kramer) of the intercepts [suRF, 0.69 ± 0.07 (95% CI); pRF, 0.94 ± 0.19] and slopes (suRF, 0.05 ± 0.01; pRF, 0.33 ± 0.02) of the two models across the population (H1–H4). We observed a significant difference of the intercepts (P = 0.01) as well as the slopes (P = 1.06 × 10−10), being smaller for the suRF model as expected.

Fig. 3.

Comparison of suRF with pRF. (A–D) Examples of suRF voxel fits for subjects H1–H4. Each column presents example voxel responses (β-weights per condition; dots) from gradually increasing eccentricity (Top to Bottom). Solid lines represent PLC–suRF model fits. (E–H) suRF size (black) and pRF size (gray) as function of eccentricity for subjects H1–H4. The eccentricity of the same voxel’s pRF was used for suRF that only estimates RF size. Lines represent linear regression fits. For parameter estimates and statistics, see .

Having established that the simple PLC model provides RF-size estimates that follow the linear relationship with eccentricity and these estimates are significantly smaller in comparison with the pRF model, we proceeded to fit the same data using the models dPLC and sPLC, which include RF-size estimates for the center (σ) as well as the surround (σ). Both models demonstrated exceptionally good fits as reflected in very high median coefficients of determination (CoDs) across all subjects, which were similar to the simple PLC model (Table 1). To select which of the three models was more supported by the data, we calculated Akaike’s information criterion (AIC), which takes into account the number of parameters included in each model (NPLC = 3, NdPLC = 5, NsPLC = 4) and penalizes models with higher number of parameters (22). Based on AIC, we also calculated the relative likelihoods of each model per voxel and plotted their distributions (). As the PLC model outnumbered the two normalization models in voxels with higher relative likelihoods and also demonstrated lower (i.e., better) mean AIC values, we chose to use this model for further reporting our results in the remaining sections of the manuscript. This is also supported by analysis of the CoD as a function of eccentricity for all voxels as subjects (). It is important to note, however, that there were some voxels, that is, in particular ones that demonstrated clear suppressive effects for intermediate checker-size conditions, for which the surround models clearly outperformed the PLC (center only) model. A comparison of the PLC and sPLC model for some selected example voxels is presented in Fig. 4. As expected and evident from these plots, PLC cannot capture the suppression effects (BOLD reduction at intermediate checker sizes), while sPLC performs very well at these voxels. However, it is important to note that the addition of a surround normalization in the model by and large did not affect the estimate of the size of the RF center, which remained similar to the estimate derived by the PLC model ().

Table 1.

Evaluation of different models in the human subjects H1–H4

Model	Nparams	Score	H1	H2	H3	H4
PLC	3	CoD	98.9%	98.8%	98.5%	99.1%
		AIC	−3.9078	−3.7378	−3.4641	−3.9815
dPLC	5	CoD	99.1%	99.1%	98.6%	99.1%
		AIC	−3.4779	−3.3309	−2.8225	−3.3747
sPLC	4	CoD	98.5%	98.0%	98.0%	98.1%
		AIC	−3.1736	−2.8642	−2.8409	−2.9591

For each model, the coefficient of determination (CoD) and Akaike's information criterion (AIC) are shown. The PLC model (bold) outperformed the other models based on the AIC and was selected for further analysis.

Fig. 4.

Comparison of PLC and sPLC. Selected example voxels for which the subtractive normalization model (sPLC) outperforms the no-normalization model (PLC) based on Akaike’s information criterion (AIC). Note that the responses as reflected in the normalized β-weights demonstrate suppression for intermediate checker sizes and that this could not be fit by the PLC model. H1–H4 reflect the different subjects.

To better understand the relationship between our proposed PLC–suRF model and electrophysiological measurements of RF sizes, we performed additional experiments in rhesus macaques. During anesthetized fMRI experiments, monkeys were presented with identical stimuli as human subjects. PRF and suRF models were then estimated with the same methodology as in humans (). As shown in the example voxels in Fig. 5 , the responses and suRF-model fits for monkeys were very similar to humans with lower eccentricity voxels demonstrating smaller RF sizes in comparison with voxels located at more peripheral locations. Linear regression analysis of the eccentricity vs. suRF and pRF sizes demonstrated identical results as those in humans (Fig. 5 and ) with all subjects showing significant linear relationships for both models (). Furthermore, we performed ANCOVA followed by the comparison (Tukey–Kramer test) of intercepts (suRF, 0.42 ± 0.05; pRF, 1.39 ± 0.15) and slopes (suRF, 0.06 ± 0.01; pRF, 0.19 ± 0.02). As in humans, we observed significant differences in both the intercepts (P = 9.56 × 10−10) and slopes (P = 9.56 × 10−10) of pRF vs. suRF sizes as a function of eccentricity. Importantly, the estimated suRF sizes approximated previously reported electrophysiological measurements of suRFs (10, 23, 24).

Fig. 5.

Validation of suRF with electrophysiology in rhesus macaques. (A–D) Example voxel responses (β-weights) and respective suRF fits (solid lines) for monkeys M1 and M2. (E and F) pRF (gray) and suRF (black) size as a function of eccentricity for monkeys M1 and M2. Each dot represents a voxel. Dashed lines: linear regression fits. (G) pRF and suRF fits across both monkeys M1 and M2 compared with electrophysiology obtained from two other animals (M3 and M4; red). (H) Example of electrophysiology RF estimate by reverse correlation (). ***P < 10−9; NS, P > 0.05. For parameter estimates and statistics, see .

To more directly investigate how suRF–fMRI estimates compare with single-neuron electrophysiological RF sizes, we performed electrophysiological RF measurements in two other monkeys (M3 and M4). To this end, we used the method of reverse correlation that has been extensively used in previous RF-mapping experiments in visual cortex (25–29). In addition, we also used a moving bar method that closely resembles the pRF mapping we used in fMRI (). An example of a recorded RF is presented in Fig. 5 (for additional examples, see ). Since the electrophysiology data from both monkeys and methods were consistent, we have collapsed them and performed linear regression for RF size vs. eccentricity like we did for the fMRI measurements, which are presented in the same figure (Fig. 5). Comparison of the suRF intercept and slope with electrophysiology (intercept, 0.16 ± 0.05; slope, 0.08 ± 0.02) showed no significant difference (intercepts, P = 0.47; slopes, P = 0.98; Fig. 5 and ).

To absolutely settle the correspondence between the suRF model and electrophysiology, we performed fMRI pRF and suRF experiments in monkey M3 that had MRI compatible implants (Fig. 6). To be able to coregister the physiology and MRI estimates, we have inserted an MRI-compatible guide (Fig. 6) on top of the grid in the chamber and filled it with an MRI contrast agent (MION). In this way, a reference frame was reconstructed, and we could use it to estimate the voxel corresponding to our recording electrode. In Fig. 6 , we show estimates of the RF size of a recording location based on all three methods (pRF, suRF, and electrophysiology). It is clear that this example demonstrates a close correspondence between the electrophysiology and suRF RF center estimates, and that they are both much smaller compared with the pRF (Fig. 6). This relationship between pRF, suRF, and electrophysiology was also present across the population (Fig. 6 ). Overall, our fMRI and electrophysiology results demonstrated that the suRF method can estimate reliably average single-neuron RF sizes in primary visual cortex.

Fig. 6.

Comparison of suRF with pRF and electrophysiology in the same animal. (A) Anatomical slice from the recording location with overlaid pRF eccentricity map. To localize the recording location, we used an MRI-compatible insert (B, Left) in the implanted recording chamber on top of the already present recording grid (B, Right). The chamber including the grid and holes of insert were filled with MRI contrast agent so we could reconstruct their position (B, Middle) and localize the position of the electrode recordings that happened in separate sessions outside the scanner (light blue arrow) to avoid artifact. (C) The pRF model fit and parameters are plotted on the fMRI signal time course from a voxel under the electrode tip (white circle). (D) suRF model fit for the same voxel. (E) Electrophysiology RF from the same location using the moving-bar method (). (F) Relative RF size for all three methods. (G) Comparison of the electrophysiology (cyan), suRF (green), and pRF (red) RF vs. eccentricity. (H) Similar to G, but only showing the low eccentricities corresponding to the electrophysiology data.

One interesting question is whether this method can generalize beyond V1, to higher visual areas. To this end, we used the PLC–suRF model to fit the V2 voxels of both human and monkey subjects. In general, the model demonstrated exceptionally good fits (Fig. 7) with very high median coefficients of determination (H1, 98.7%; H2, 97.6%; H3, 98.7%; H4, 98.7%; M1, 90.8%; M2, 97.2%). Similar to area V1, the results demonstrated a significant linear relationship of suRF as well as pRF size with eccentricity (Fig. 7 ) for all individual subjects and across populations (). We again tested whether the suRF-estimated RF sizes in V2 were smaller than the pRF by performing ANCOVA and second-level comparisons (Tukey–Kramer). In humans, the intercepts [suRF, 0.8 ± 0.09 (95% CI); pRF, 1.4 ± 0.15] and slopes (suRF, 0.10 ± 0.01; pRF, 0.33 ± 0.02) of the two models across the population (H1–H4) were significantly different both for the intercepts (P = 1.1 × 10−10) as well as the slopes (P = 1.06 × 10−10). Similar results were obtained in monkeys, with both intercepts [suRF, 0.48 ± 0.15 (95% CI); pRF, 2.07 ± 0.26; P = 9.56 × 10−10] and slopes (suRF, 0.11 ± 0.02; pRF, 0.25 ± 0.03; P = 9.61 × 10−10) found again to be significantly different. Furthermore, we compared the intercept and slope estimates in areas V1 and V2 using a t test. We found that, across the human population, the slope for V2 was significantly larger than V1 (T = 6.68; P = 3.03 × 10−10) and there was also a trend for a larger intercept in V2 (T = 1.95; P = 0.051). Similar results were obtained across the monkeys with significant differences in the slope (T = 5.81; P = 7.64 × 10−9). lists the results per human and monkey subject.

Fig. 7.

suRF model in area V2. (A) Examples of suRF voxel fits for subjects H1–H4 and M1 and M2. Each subplot presents V2 example voxel responses (β-weights per condition; red dots), and solid red lines represent PLC–suRF model fits. (B–G) suRF size in V2 (red) and suRF size in V1 (gray) as function of eccentricity for subjects H1–H4 and M1 and M2. The eccentricity of the same voxel’s pRF was used for suRF that only estimates RF size. Lines represent linear regression fits. For parameter estimates and statistics, see .

Last, to understand whether the suRF model is robust to potential inaccuracies in the model assumptions, we performed further modeling and modulated the model parameters. First, we wanted to know whether the compressive nonlinearity is necessary or whether a simplified linear filter model could perform similarly well. To this end, we have set the exponent n = 1 (i.e., no nonlinearity) and refitted the model. We found that this model performed very poorly with a median coefficient of determination 78.7% in comparison with the PLC model (98.8% for this subject). To directly observe the goodness of fit of the linear model, we then plotted the fits of the 16 best-fitting voxels (). Given these fits, we concluded that a compressive nonlinearity is necessary; otherwise, the model cannot follow the structure of the data. Second, we investigated how the results are affected by the choice of the bandwidth we used for the Gabor filters. To this end, we fitted the PLC model for a range of bandwidths within the physiological range [n = (0.15, 0.35); ]. We found that, although the slopes of the linear eccentricity vs. RF-size relationship changed slightly (increased with n), the effect was minimal in comparison with the differences between suRF and pRF (). Moreover, since an fMRI voxel would probably contain cells with more than a single bandwidth, we wanted to test how this could affect the suRF model. To do so, we created a multibandwidth model that simply uses the average response of arbitrary distributions of Gabor bandwidths entered by the user. For example, we used this model with a homogeneous distribution of bandwidths simulating 10,000 cells within the range [0.15, 0.35]. The result is shown in and demonstrated that this multiband model was strikingly similar to our original model using a bandwidth n = 0.25 that is approximately in the center of the physiological range.

Discussion

To date, only two studies have attempted to report measurements of RF sizes in human primary visual cortex either with invasive intracranial electrophysiology (potentially one to two cells) (30) or surface electrocorticography in patients (31). Moreover, RF-size estimation from two intracranial electrodes in areas V2/V3 was also recently performed in a human subject (32) and indicated consistent results to previous electrophysiology in macaques. Here, we developed a method to estimate suRF sizes across human early visual cortex (areas V1 and V2) using in vivo fMRI and validated our approach by dual fMRI and electrophysiology measurements in rhesus macaques.

To achieve this, we took advantage of the well-known property that visual cortical neurons can be approximated by Gabor functions, which act like SF filters. We combined this property with a very simple stimulus design using checkerboard patterns whose checker size was modulated across different presentation blocks, resulting in stimulus conditions with substantially different SF content. These stimuli contain a broad range of frequencies, but importantly the average SF content can be estimated analytically (). This allowed straightforward computational modeling that predicts BOLD responses based on the Gabor parameters, and specifically on the SD σ of the Gabor envelope that is proportional to RF size.

Previous work has demonstrated the utility of encoding models either for extracting voxel-based neuronal properties from population-recording data such as fMRI, or for decoding brain activity (5, 33–38). The pRF (5) is an example of such a model that uses the very simple binary presence or absence of a stimulus on a particular location and the corresponding responses to build the aggregate RF for each voxel, assuming it is a Gaussian function over visual field space. Another example is the Gabor wavelet pyramid model (34), which was used to decode natural images in the human brain while at the same time it returns the location, size, SF, and orientation preference of each voxel based on the weights of the different Gabor functions in the underlying battery. An important difference of all of these models to the suRF model we introduce here, however, is that they work in the spatial visual field domain and thus estimate voxel-based RFs that include the scatter of neuronal RFs within a voxel. The suRF model, on the contrary, sacrifices the estimation of the spatial positions of the RFs and works directly in SF domain to estimate the RF sizes. This is advantageous as it makes it the only model that is not sensitive to RF scatter (neither random nor organized) and, as we showed, results in exceptionally good estimates of V1 RF sizes closely matching invasive electrophysiological measurements. It is important to note that the difference in RF sizes estimated by pRF, or other spatial domain methods, compared with suRF, do not reflect a discrepancy between the methods but rather a difference in the quantities estimated: that is, the aggregate RF of all neurons together vs. the average single-neuron RF. In future studies, it will be interesting to study the relationship of these estimates and evaluate whether one can be predicted reliably by the other.

In this work, we also attempted to go beyond estimating a single RF size per voxel by fitting two slightly more complex models that included surround normalization, which has been proposed to be a canonical computation in the brain (20). To this end, we created models for both a divisive (dPLC) as well as subtractive (sPLC) normalization by a second pool of neurons with different average RF sizes (σ). In general, similarly to the simple PLC model, both of these models resulted in exceptionally good fits with very high median coefficient of determination (Table 1). Importantly, adding the additional parameters by and large did not affect the estimation of the size of the RF center (). When we examined voxels in which the normalization models performed substantially better in comparison with the PLC, we observed that such voxels typically demonstrated suppressive effects at intermediate checker size stimulus conditions. In particular, the sPLC model seemed to perform slightly better especially at voxels that showed strong suppressive effects eliciting very low or even negative activity in comparison with the blank baseline (Fig. 4). However, the proportion of voxels in which these models demonstrated higher likelihoods in comparison with the PLC model was quite low, and thus we further reported only the results of the PLC. Given that these surround models had more parameters than the PLC and we only measured six stimulus conditions to fit the parameters with, it is hard to provide strong evidence for those models. Future studies could tackle this problem by using a wider range of stimulus conditions.

An interesting result of our modeling approach was that the structure of the data supported using a compressive static nonlinearity (). In the simple PLC model, we allowed this parameter to be estimated and found that its value approximated previous V1 modeling studies that used the nonlinearity in a different context and specifically compressive spatial summation within the RF (21). Even if this may at first appear puzzling, we conjecture that the compressive nonlinearity we find likely reflects a homologous phenomenon seen from the perspective of the SF domain. Perhaps, partially activating the RF of neurons in space produces a similar effect as produced by the partial overlap between the SF content of the stimuli and the SF sensitivity of the Gabor functions. Another potential explanation of the nonlinearity is the estimation of neural activity R. A usual model of estimating neural activity is the complex energy model (39). Many studies have used the square root of the complex energy model to predict neural activity assuming linearity. Our model is using a similar approach in the frequency domain by using the spectral amplitude. The single PLC nonlinearity we use is applied after this stage and can thus reflect compression at the neural or hemodynamic level (i.e., conversion of neural response to BOLD signal) or even a combination of both. A recent study using an intelligent design in an achiasmic human and metaanalysis of a number of previous studies in humans as well as one in macaques claimed that compressive nonlinearity with approximate power of 0.5 can account for all results in early visual cortex (40). Since the power best fitting the results in our study was smaller (i.e., more compressive), it could indicate that our nonlinearity reflects a combination of neural as well as hemodynamic effects. However, identifying the exact relationship was beyond the scope of this work.

To understand whether the suRF modeling approach can generalize to other visual areas, we used our method in human and monkey area V2. We found that the model performed very well and returned RF sizes that demonstrated the linear relationship with eccentricity as expected. In addition, the slopes and sizes of RFs in V2 were larger than V1, further demonstrating the validity of the approach. Moreover, the estimates of RF sizes in macaque V2 are consistent with previous electrophysiology studies (24, 41). We should note, however, that no direct electrophysiological validation was performed in V2. We have not yet tried this approach in higher visual areas. We expect that, for higher visual areas, more complex models reflecting the local electrophysiological signatures would be necessary to retrieve accurate estimates of neuronal properties.

RFs of neurons in the visual cortex of experimental animal models have been studied extensively using a wide variety of invasive techniques. However, since invasive studies are not possible in healthy humans, much less is known about human single-neuron RFs. Our study provides comprehensive human suRF measurements over the whole primary visual cortex using a simple noninvasive fMRI technique that can be easily adapted to the general population. Visual RF sizes have been reported to reorganize after injury (42–49), and studies suggest that visual stimulation could induce further reorganization and plasticity (50, 51). Human studies using the pRF as a readout have also demonstrated changes; however, it is not yet known whether these are a result of single-neuron RF-size changes vs. changes in the spatial distribution of RFs within a voxel. The suRF methodology we propose here can be useful for further dissecting the mechanisms of neural plasticity in human subjects with lesions of the visual pathways.

Materials and Methods

Human Experiments.

Subjects.

Four healthy human subjects (H1–H4; 24–44; two females) with normal or corrected-to-normal visual acuity participated in the fMRI experiments. Before each session, subjects provided written informed consent. The local ethics committee of the University Hospital Tübingen approved the study.

MRI data acquisition and preprocessing.

MRI experiments were performed at the Max Planck Institute for Biological Cybernetics (Tübingen, Germany). Functional and anatomical images were acquired in a 3.0-T Tim Trio Scanner (Siemens) using a 12-channel coil. At least two T1-weighted anatomical volumes were acquired for each subject with a 3D magnetization-prepared rapid acquisition gradient echo (T1 MPRAGE scan) and averaged following alignment to increase signal-to-noise ratio [matrix size, 256 × 256; voxel size, 1 × 1 × 1 mm3; 176 partitions; flip angle, 9°; repetition time (TR), 1,900 ms; echo time (TE), 2.26 ms; inversion time, 900 ms]. BOLD image volumes were acquired using gradient echo sequences of 28 contiguous 3-mm-thick slices covering the entire brain (TR, 2,000 ms; TE, 40 ms; matrix size, 64 × 64; voxel size, 3 × 3 × 3 mm3; flip angle, 90°).

At least five functional scans were acquired for each subject, consisting of 195 image volumes, the first three of which were discarded. The functional images were corrected for motion between and within scans (52) and were aligned to the high-resolution anatomical volume using a mutual information method (53). The high-resolution anatomical data were used to segment the white/gray-matter boundary in itkGray software, and 3D cortical surface and flat mesh models were created and realigned with the functional data. The functional time series were spatially resampled in the volume space using nearest neighbor interpolation. This preserves the original signals but up-samples them in space leading to some voxels with the same time courses. Analysis was accelerated by analyzing a single voxel corresponding to the original fMRI resolution and assigning the results to the corresponding anatomical voxels. All subsequent analysis was performed in the segmented volume space restricted in the gray-matter voxels. The above preprocessing steps were performed in MATLAB using the mrVista software package that can be found at https://github.com/vistalab/vistasoft.

pRF mapping.

For retinotopic mapping, we used the pRF method (5). Briefly, the pRF model estimates the region of the visual field that effectively elicits a response in a small region of visual cortex (voxel). The implementation of the pRF model is a circularly symmetric Gaussian RF in visual space, whose center and radius are estimated by fitting actual BOLD signal responses to estimated responses elicited by convolving the model with the moving bar stimuli. We retained only those voxels in the visual area for which the pRF model explained more than 15% of the variance. This threshold was set after measuring the mean explained variance (EV) in a nonvisually responsive area and setting the threshold at 6 SDs above the mean. This method derived reliable and reproducible retinotopic and pRF size maps.

Stimuli.

Stimulus presentation.

Subjects were presented with visual stimuli in the scanner by using MRI-compatible digital goggles [VisuaStim; Resonance Technology Company; 30° horizontal and 22.5° vertical field of view (FOV); 800 × 600 resolution; minimum luminance, 0.3 cd/m2, and maximum luminance, 12.2 cd/m2].

pRF mapping.

Moving square-checkerboard bars (100% contrast) were presented within a circular aperture with a radius of 11.25° (full vertical extent of screen) around the fixation point. The bar width was 1.875° and traveled sequentially in eight different directions, moving by a step half of its size (0.9375°) every image volume acquisition (TR, 2 s). Stimuli were generated using Psychtoolbox-3 (54) and an open toolbox (VISTADISP) in MATLAB (The Mathworks). The subjects’ task was to fixate a small dot in the center of the screen (radius, 0.0375°; 2 pixels) and respond to the color change by pressing a button. The color was changing randomly with average frequency of one every 6.25 s. An infrared eye tracker was used to record eye movements inside the scanner (iView XTM; SensoMotoric Instruments).

suRF mapping.

Binary full-field checkerboards (black/white; 100% contrast) with a frame rate of 30 Hz (new random configuration was presented in each frame) were presented in a block design (10 s ON, 20 s OFF). For each block, the size (side) of the checkers was chosen to be one of six possible conditions [0.25°, 1°, 2°, 3.75°, 7.5°, 22.5°].

Monkey Experiments.

Four healthy adult rhesus monkeys (Macaca mulatta; M1–M4, 5–11 kg, one female) were used in MRI (n = 3; M1–M3) and electrophysiology experiments (n = 2; M3 and M4). The experimental and surgical procedures were performed with great care and in full compliance with the German Law for the Protection of Animals, the European Community guidelines for the care and use of laboratory animals (EUVS 86/609/EEC), and the recommendations of the Weatherall report for the use of nonhuman primates in research. The regional authorities approved our experimental protocols, and the institutional representatives for animal protection supervised all procedures. Animals were kept in large cages located adjacent to the training and experimental facilities. Space in these cages allows swinging and jumping, and enrichment equipment such as toys was changed frequently. Group housing was maintained to increase the quality of life by rich visual, olfactory, auditory, and social interaction and stimulation for play. Balanced nutrition and regular veterinary care and monitoring were provided. Recording chamber implantations were performed in two animals while the animals were under general anesthesia and aseptic conditions. To alleviate postsurgical pain, we administered analgesics for a week after the surgery (see below). MRI experiments were also performed under anesthesia (see below). Animals were not killed after the experiments.

Anesthesia.

Details on the anesthesia protocol have been given previously (55). Briefly, the animals were premedicated with glycopyrolate [0.01 mg/kg, intramuscular (i.m.)] and ketamine (15 mg/kg, i.m.), and then deep anesthesia was induced by fentanyl (3 μg/kg), thiopental (5 mg/kg), and succinyl chloride (3 mg/kg). Anesthesia was maintained with remifentanyl (0.5–2 μg⋅kg−1⋅min−1) under paralysis with mivacurium chloride (3–6 mg⋅kg−1⋅h−1) to ensure the suppression of eye movements. Heart rate and blood oxygen saturation were monitored continuously with a pulse oximeter. Body temperature was kept constant at 37–38 °C.

Surgical procedures.

Recording chambers were positioned over the operculum in area V1 according to stereotaxic coordinates. This was aided by high-resolution magnetic resonance anatomical imaging. The anatomical scan and recording chamber implantation were done while the animals were under general anesthesia. Details of the procedure can be found in our previous work (55).

fMRI experiments were performed in a 4.7- or 7.0-T vertical scanner (Bruker Biospec; Bruker Biospin). Multislice fMRI was performed using eight-segmented gradient-recalled echo-planar imaging (EPI). Volumes of 17 slices of 0.756 × 0.756 × 2 mm3 were collected, each with a FOV of 96 × 96 mm on a 128 × 128 matrix and 2-mm slice thickness, flip angle of 40° for 4.7 T and 47.6° for 7 T, a TE of 20 ms, and a TR of 750 ms per segment resulting in a volume acquisition time of 6 s. For anatomical measurements, we used a FLASH sequence with the same FOV, 96 × 96 mm2; matrix, 256 × 256; slice thickness, 2 mm; flip angle, 70°; TE, 10 ms; and TR, 2,000 ms. A high-resolution 3D modified driven equilibrium Fourier transform (MDEFT) anatomical scan with an isotropic resolution of 0.5 mm3 was acquired for coregistration with the FLASH and EPI images. For more details on the fMRI acquisition methods, see previously published papers (55, 56). Then, fMRI data were reconstructed and imported into a MATLAB-based toolbox (mrVista). The high-resolution 3D-MDEFT anatomical data were used to manually segment the white/gray-matter boundary in itkGray software, and 3D cortical surface and flat mesh models were created and realigned with the functional data by using mrVista. The functional images were corrected for motion in between and within scans (52).

For retinotopic mapping, we used the pRF method. The threshold of the EV was set to 15%, the same with our human experiments and in agreement with previous studies (46, 57).

MRI stimuli.

Visual stimuli were displayed using a custom-made fiber-optic projection system at a resolution of 800 × 600 pixels (30° × 22.5°) with a 60-Hz frame rate and mean luminance of ∼50 cd/m2. Stimuli were centered on the approximate location of the fovea by using a modified fundus camera (Zeiss; RC250). Animal eyes were fitted with appropriate contact lenses to ensure the stimulus remained in focus. At the beginning of each experiment, a polar pattern with radius 5° was presented monocularly to the left and right eyes in a block design (40 s ON–40 s OFF), and the results were analyzed on-line to select the eye with best alignment and the strongest responses. Further stimulus presentation was restricted to this eye.

Moving square-checkerboard bars (100% contrast) with width 1° were moving in 0.5° steps every volume acquisition on the full screen extent sequentially in four different directions (left–right, up–down, right–left, down–up). Outside the bar aperture, an isoluminant gray background was presented. Stimuli were generated using custom-made stimulus generation and presentation software (MRIstim). Each direction was presented twice for each scanning acquisition of 204 images. This was repeated two to three times.

The suRF mapping binary checkerboard stimuli were generated with the same program creating the stimuli for the human subjects using Psychtoolbox-3 (54) in MATLAB (The Mathworks). The only difference was that the block design ON–OFF periods were 24 and 36 s, respectively, and we used eight possible checker size conditions [0.25°, 0.5°, 1°, 2°, 4°, 6.5°, 13.0°, and 22.5°]. Since the animals were under anesthesia, no fixation spot was presented.

Electrophysiology data acquisition and preprocessing.

Electrophysiological recordings were acquired from area V1 in two monkeys (M3 and M4). Recording chambers were positioned stereotactically over the operculum with the aid of T1-weighted high-resolution 3D-MDEFT anatomical MRI images (0.5 mm isotropic). These images were used for targeting the region of interest, 3D-skull reconstruction, and design of the implants to accurately fit the contours of the skull. A five-axis CNC machine (Willemin-Macodel; W428) was used to build these form-specific implants that resulted in an excellent fit between the implants and the underlying skull surface. In one animal (M4), the implant was constructed from medical-grade titanium, which precluded MRI measurements. In the second animal (M3), the chamber was constructed from polyether ether ketone plastic (TECAPEEK; Ensinger), which is MRI compatible, and we could thus also acquire pRF and suRF measurements.

In M3, recordings were conducted primarily using custom-made tetrodes [see previous work for details (58); n = 42 experiments], and in some (n = 4), a custom-made multichannel linear probe with 10 platinum/iridium channels with interdistance ∼150 μm was used in parallel. After eliminating noisy channels with no or very little spiking activity and channels recorded from area V2, this resulted in n = 35 datasets from area V1. In M4, recordings were conducted using either single-channel tungsten electrodes (FHC; n = 8 experiments) or custom-made multichannel probes made with platinum/iridium wire (nine linearly arranged recording locations with interdistances of ∼150 μm; the eighth location was consisting of three nearby channels; n = 10 experiments). After eliminating noisy channels with no or very little spiking activity and channels recorded from areas V2 and V4, this resulted in n = 38 datasets that were classified as recordings from V1.

The electrodes were guided into the brain manually by custom-designed adjustable microdrives. Electrophysiological activity was sampled at 32 kHz, digitized (16 bits), and stored using the Digital Lynx data acquisition system (Neuralynx). Multiunit activity was defined as the events that exceeded a predefined threshold (25 V) of the filtered, digitized signal (digital bandpass, 600 Hz to 6 kHz). One animal (M4) was implanted with a scleral search coil (59, 60), and its eye movements were monitored on-line. In the second animal (M3), eye movements were monitored by an infrared eye tracker (IView XTM Hi-Speed Primate; SMI). The behavioral aspects of the experiment were controlled using the QNX real-time operating system (QNX Software Systems).

Electrophysiology stimuli.

Visual stimuli were displayed using a dedicated graphics workstation (TDZ 2000; Intergraph Systems) running an OpenGL-based stimulation program (STIM). Stimuli were presented on a LCD monitor positioned at 1-m distance from the animals’ eyes (width, 60 cm; height, 34 cm) with a resolution of 1,920 × 1,080 and a refresh rate set to 100 Hz. The monitor was gamma corrected with a mean luminance of 22.2 cd/m2.

Reverse correlation RF mapping.

To map the multiunit neuronal RFs, we used the reverse-correlation technique (25–29). Stimuli were small square dots (black or white) presented on a gray background while the monkey fixated a small red spot with 0.2° diameter at the center of the monitor. The dots were positioned inside a rectangular 5 × 5 grid with side dimensions of 0.5°–2.5° (i.e., 25 possible dot locations with dot sizes one-fifth the rectangle side) and location depending on a preliminary manual mapping of the RF location that was placed approximately in the center of the grid. Each dot was presented for 20 ms at pseudorandomized locations to have approximately equal number of presentations in each position inside the grid for each luminance (black and white).

Moving-bar RF mapping.

In one of the animals (M3), in parallel to the reverse-correlation mapping, we also performed RF mapping with moving bars; this was similar to the pRF mapping performed in humans with the bars moving in eight different directions across the monitor in steps of 45°. The bar width was 0.5° and was moving by steps of 0.25° every 10 ms.

Electrophysiology analysis.

Eye movement analysis.

First, we calculated the time series of eye velocities by differentiation of the position signals. Then, the horizontal and vertical angular velocities were independently thresholded at seven times their median-based SD to detect putative microsaccadic events. An event was classified as a microsaccade if the following additional criteria were satisfied: (i) it had a minimum duration of 8 ms, (ii) it had an amplitude between 1 and 60 min of a degree, and (iii) it had a maximum peak-velocity of 110° per second (61). These parameters provided accurate detection of the microsaccades. In addition, the extracted microsaccades satisfied the main-sequence criterion and showed high correlation of amplitude and velocities. Fixation locations were extracted as the mean positions between saccades.

Reverse-correlation RF mapping.

The multiunit RFs at each recording location were reconstructed by iterative construction of spike histograms for each stimulus position. Each bin had a time span of 10 ms, and in total we considered 12 bins (120 ms) post stimulus presentation. For each presentation of a dot (separately for black-and-white dots), the spikes following were binned according to the eye movement-corrected position of the dot. Eye movement correction improved substantially the estimation of RFs in M4 that had an eye coil. In contrast, in M3 (infrared eye tracker), this procedure did not demonstrate improvements or even worsened RF reconstruction, and thus it was not applied. RF sizes with or without the correction in M3 were very similar with some changes in RF position that could reflect drifts of the infrared eye tracker. From the visual field position histograms, we then created spatial maps for each time bin. Separate maps for the dark and bright dots as well as the average of the two were created. These maps were converted to z scores by subtracting the mean and dividing with the SD of the first three time-bin maps (i.e., <30 ms where no substantial response is expected or was observed). Then, the bin with maximal responses was identified and used for fitting a 2D-isometric Gaussian with its parameters reflecting RF position (x, y) and size (2σ). This procedure was performed for each map type (i.e., dark, bright, mean), and we selected the parameters of the model providing the best fit based on the coefficient of determination (R2). Examples of RF data and fits are presented in .

For the analysis of the moving-bar RF mapping, we have used a procedure similar to the pRF topographic mapping approach we proposed earlier for fMRI data (7). For the analysis of the electrophysiological data, an additional step was required before estimation of the RF. Specifically, a time delay reflecting the neuronal response latency was estimated so that the peaks of the responses when the bar was running in opposite directions were aligned. Also, in contrast to the fMRI analysis pipeline, no hemodynamic deconvolution was necessary. All other steps of analysis were identical to the fMRI.

Modeling and Analysis.

To estimate the suRF size, the amplitudes of the responses (βBOLD) to the block design of checkerboard stimuli with different checker sizes λc presented in each block were estimated using a GLM. Only voxels that demonstrated GLM EV greater than 15% were further considered. The β’s for all conditions were then used to fit models (Eqs. –) by using nonlinear least squares (lsqcurvefit function in MATLAB; range of suRF size between 0° and 10°). These models use a static nonlinearity to transform estimated neural responses to BOLD signal.

For the purpose of modeling the RFs, we assume that neurons in primary visual cortex (V1) can be modeled as linear combinations of even (g) and odd (g) 2D-Gabor functions, described with the following two equations, respectively:

with x and y as the Cartesian coordinates in the visual field, σ and σ defining the ellipsoid 2D Gaussian envelope, and ω, ω defining the center frequency and orientation. The SF sensitivities ω, ω of these Gabor filters can then be described by the respective Fourier transforms:

Furthermore, we assume that the distribution of orientation selective cells within an fMRI voxel is homogeneous; thus, the SF sensitivity of voxels can be described relative to a single center frequency ω0 independent of orientation as follows:

Thus, the SF sensitivity of the population of cells within a voxel can be estimated according to a linear combination of the Eqs. and . For simplicity, we take the average of the even and odd parts as follows: ().

Then, for estimating the neural responses R (Eq. ), we integrated the product of ΣG with an analytical form of the stimulus frequency content S ():

Data Deposition.

The data reported in this paper have been deposited in the MEGA database, https://mega.nz/#F!rjxGTYRb!6c_XFeslILmSlt_QMfulpQ (62).