V1 neurons are tuned to perceptual borders in natural scenes

To derive the RMS contrast, we computed the contrast of every natural image within each site/neuron’s RF or voxel’s pRF. The (p)RF was modeled as a circular symmetric Gaussian function, described by parameters for position (x_c, y_c) and size (σ), giving rise to a Gaussian weighting function w_i:

where x_c and y_c define the location of the center of the (p)RF in the visual field, σ determines the size of the (p)RF, and x_i and y_i define the location of the i-th pixel of the image. We computed each site/neuron/voxel’s contrast value to each natural image by calculating the RMS contrast (39, 40) of the part of the image inside their (p)RF. RMS contrast was defined as the SD of the luminance of the pixels relative to the mean. The RMS-contrast was weighted by the (p)RF Gaussian function to obtain the local contrast-energy value:

where N is the number of pixels in the stimulus window, L is the mean luminance of the pixels inside the spatial window, L_i is the luminance of the i-th pixel (2), and w_i comes from Eq. 1.

The magnitude of the complex-cell response was computed as

where E is the dot product of the stimulus (grayscale) pixel luminance with a 2D Gabor with coordinates x and y, F is spatial frequency, θ orientation, and Φ phase. For neurons in anesthetized monkeys, F and θ were computed using a separate set of stimuli containing Gabor patches with varying spatial frequency, orientation, and phase. For recording sites in awake monkeys and fMRI voxels, F and θ were free parameters of the model and were optimized for every recording site before we calculated the explained variance. We used Matlab’s lsqcurvefit function to compute the fits.

We extracted the activity of units of VGG-19’s layer conv3_1, which is state of the art in predicting V1 responses to natural images (6, 29), and followed the approach of ref. 6 with two modifications. We used a two-step mapping (41, 42), described by following the equation:

where $r_{VGG-19energy}$ is the predicted response of a V1 recording site/neuron, $f_{VGG-19}$ (input) is the output of VGG-19’s conv3_1 to our stimulus set (i.e., $\mathit{input}$), and $W_s$ and $W_d$ are two sets of weights defining spatial and feature selectivity, respectively. The spatial mask ($W_s$, set as the 2D Gaussian RF estimate) approximates the RF and a weighted sum of the nodes in the ANN ($W_d$) approximates the feature selectivity of the recorded sites/neurons (41). $W_s$ and $W_d$ were determined separately for every recording site/neuron, while training occurred at the same time for all of them. We trained the model to optimize $W_d$ to predict V1 responses to the training set (i.e., $\mathit{input}$ in Eq. 4). We cross-validated the model for each site/neuron and used the final model to extract VGG-19 energy for each trial. For recording sites in the awake monkeys, we cross-validated using random 50% splits of the trials. For neurons in the anesthetized monkeys, we used 217 images for training, and 27 images for testing. For the anesthetized monkeys, the training images were fixed and the test sets differed across neurons because the test images depended on the presence of object borders in the RF. To avoid overfitting, we fitted models (i.e., optimized $W_d$) using images within the 25th-75th percentiles of object border relevance (see below in the section on quantification of object borders) and used the other images during testing. We tested the object BoM using VGG-19 energy with the 27 held-out stimuli, only including neurons for which VGG-19 energy predicted the neural responses on the test images (N = 26, Pearson’s correlation r > 0.323, P < 0.05 uncorrected, t test).

For RMS contrast (and for the models of tuning), we computed the CRF of the MUA recording sites and single neurons in V1, and voxels of areas V1, V2, V3, hV4, LO-1/2, and V3-a/b (see below) by measuring their responses as a function of local contrast (or the model energy) inside their (p)RF. We chose energy bins such that every bin spanned 10% of the contrast (or energy) distribution and fitted the following equation (modified from ref. 43):

where R is the neural response, C is the contrast (or model energy) inside the pRF, Q represents the model energy value where R is at half of its maximum response, and q determines the slope (a, Q, and q are free parameters). The number of bins did not influence the outcome.

The BSD images have been annotated by 5 to 9 human observers who drew lines to segment object borders that are relevant for the scene’s interpretation (1, 19). We used these measurements to define the perceptually relevant borders of the scene. Every pixel i of the manually labeled images has a value for the confidence of observers that the pixel belongs to an object border, S_i, between 0 (not labeled by any observer) and 1 (labeled by all observers). The border relevance in the pRF is calculated as a weighted sum across pixels that fall in it:

Here, w_i are the weights of the (p)RF estimate (Eq. 1), and N is the total number of pixels in the (p)RF. For every (p)RF, we computed the distribution of border relevance across images, and we included the lowest quartile of the distribution as “nonborder patches” and the highest quartile as “object borders.” We then computed the CRFs within these classes (SI Appendix, Fig. S1). These thresholds were not critical because we replicated the results when we included all responses, i.e., with a median split (not shown).

We used a bootstrapping procedure to determine the significance of differences in CRFs between conditions. We sampled the images with replacement 1,000 times, fit the CRF for the two simulated conditions, and computed the mean difference in the area under the curve of the two conditions. We derived the P-value from this distribution of differences.

All procedures complied with the NIH Guide for Care and Use of Laboratory Animals and were approved by the institutional animal care and use committee of the Royal Netherlands Academy of Arts and Sciences. Two macaque monkeys (males, 7 and 13 y old) participated in the electrophysiological experiments. They were socially housed in stable pairs in a specialized primate facility with natural daylight, controlled humidity, and temperature. The home cage was a large floor-to-ceiling cage which allowed natural climbing and swinging behavior. The cage had a solid floor, covered with sawdust, and was enriched with toys and foraging items. Their diet consisted of monkey chow supplemented with fresh fruit. Their access to fluid was controlled, according to a carefully designed regime for fluid uptake. During weekdays the animals received water or diluted fruit juice in the experimental set-up upon correctly performed trials. We ensured that the animals drank sufficient fluid in the set-up and supplemented extra fluid after the recording session if they did not drink enough. On days of the weekend, they received at least 700 mL water in the home cage in a drinking bottle. The animals were regularly checked by veterinary staff and animal caretakers and their weight and general appearance were recorded daily in an electronic logbook during fluid-control periods.

We implanted both monkeys with a titanium headpost (Crist instruments) under aseptic conditions and general anesthesia as reported previously (44). The monkeys were trained to direct their gaze to a 0.5° diameter fixation dot and hold their eyes within a fixation window (1.1° diameter). They then underwent a second operation to implant 5 × 5 arrays of microelectrodes (Utah-probes, Blackrock Microsystems) over opercular V1. The interelectrode spacing of the arrays was 400 μm. We obtained good signals from 4 V1 arrays in each monkey (13).

We recorded neuronal activity of 192 recording sites in V1 (96 in Monkey M and 96 in Monkey B). We recorded the envelope of multiunit activity by digitizing the signal referenced to a subdural electrode at 24.4 kHz. The signal was band-pass filtered (2nd order Butterworth filter, 500 Hz to 5 KHz) to isolate high-frequency (spiking) activity. This signal was rectified (negative becomes positive) and low-pass filtered (corner frequency = 200 Hz) to produce the envelope of the high-frequency activity, which we refer to as MUA (45). The MUA signal reflects the population spiking of neurons within 100 to 150 μm of the electrode and the population responses are very similar to those obtained by pooling across single units (45–48).

To normalize MUA, we first subtracted the mean activity in the pretrial period in which the animal was fixating (200 to 0 ms relative to stimulus onset) and divided by the maximum smoothed (26 ms Gaussian kernel) peak response (0 to 150 ms after stimulus onset). Neuronal activity was normalized to the peak response elicited by stimulus onset. The data are therefore in normalized units, where, e.g., a value of 0.1 indicates 10% of the maximal MUA onset response. We only included recording sites on days with a sufficient signal-to-noise ratio (SNR_DAY). SNR_DAY was estimated by dividing the maximum of the initial peak response by the SD of the baseline activity across trials. When the SNR_DAY of a recording site was smaller than 2 on a particular day, we removed that session from the analysis of that recording site. Approximately, 80% of recording sites in monkey B and 77% in monkey M were included. To test for statistical differences between conditions and to compute the CRFs, MUA was averaged in a time window 25 to 75 ms after stimulus onset.

In the experiments with awake monkeys, stimuli were presented on a CRT monitor at a refresh rate of 60 Hz and resolution of 1,024 × 768 pixels viewed from a distance of 46 cm. The monitor had a width of 40 cm, yielding a field-of-view of 41.6 × 31.2°. All stimuli were generated in Matlab using the COGENT graphics toolbox (developed by John Romaya at the LON at the Wellcome Department of Imaging Neuroscience). The eye position was recorded using a digital camera (Thomas recordings, 250 Hz frame rate).

We mapped the RF of each MUA recording site in V1 of the awake monkeys using a drifting luminance-defined bar that moved in one of four directions. The response to each direction was fitted with a Gaussian function. The borders of the RF were then calculated as described previously (45). The signal-to-noise ratio (SNR_RF) of the response was taken as the peak of the Gaussian divided by the SD of the pretrial baseline response. We only included recording sites in the analyses with a reliable visual response (i.e., the responses to all four bar directions had an SNR_RF of at least 1). The median V1 RF size, taken as the square root of the area, was 1.8° (range 0.4° to 8.2°) and the median eccentricity of the RFs was 2.4° (range 0.6° to 12.9°) (49).

To compute the latency of the difference between the response to object borders and nonsegmented borders, a function was fitted to the time course of the difference response (12, 13, 24). The function was derived from the assumptions that the onset of the response difference has a Gaussian distribution and that a fraction of the response dissipates exponentially which yields the following equation:

where $G(\mathrm{t},\mu ,\sigma )$ is a cumulative Gaussian density with mean $\mu$ and SD$\sigma$, ${\alpha }^{-1}$ is the time constant of the dissipation, and $c$ and $d$ represent the contribution the nondissipating and dissipating components, respectively. The latency was defined as the point at which the fitted function reached 33% of its maximum.

Four BSD images were used in the electrophysiological experiments (11.6° radius visual angle; SI Appendix, Fig. S2). At the start of the trial, the screen was gray (26.8 cd m⁻²) with a red fixation point with a position that was randomly selected from uniformly spaced grid (with ~500 positions) covering the circular aperture of the image. The image appeared once the monkey had maintained fixation for 300 ms (fixation 1). After an additional 400 ms, the first fixation point disappeared and the next fixation points appeared, and these shifts of the fixation point occurred five times per trial. Reward was delivered after every correct fixation, with an extra amount at the end of the trial. Here, we mainly focused on the data during the first fixation to facilitate the comparison with the data from anesthetized monkeys and human fMRI. We collected a total of 11,783 correct trials for monkey M and 13,373 for monkey B.

We compared differences between the CRFs or model output between object borders and nonborder patches using a bootstrapping procedure (1,000 iterations), as described above.

Four participants (all male; ages 29 to 41 y) participated in the fMRI experiment. All participants had normal or corrected-to-normal visual acuity. We obtained informed written consent of the participants and the protocol was approved by the Human Ethics Committee of University Medical Center Utrecht.

The visual stimuli were generated in Matlab (Mathworks Inc.) using the PsychToolbox (39, 50) on a Macintosh MacBook Pro. The stimuli were backprojected on a display inside the MRI bore. The subject viewed the display through mirrors inside the scanner. The size of the display was 15.0 × 7.9 cm with a resolution of 1,024 × 538 pixels. The total distance from the subject’s eyes to the display was 41 cm. The stimuli were constrained to a circular area (radius, 5.5°) with the size of the vertical dimension of the screen. The area outside this circle was maintained at a constant mean luminance.

The MRI data were acquired with a Philips 7T scanner using a 32-channel head coil (18). We scanned the participants with a 2d-echo-planar-imaging sequence with 25 slices oriented perpendicular to the calcarine sulcus with no gap. The following parameters were used: repetition time (TR) = 1,500 ms, echo time (TE) = 25 ms, and a flip angle of 80°. The functional resolution was 2 × 2 × 2 mm and the field of view (FOV) was 190 × 190 × 50 mm. We used foam padding to minimize head movement. The functional images were corrected for head movement between and within the scans. For computation of the head movement between scans, the first functional volumes for each scan were aligned. Within-scan motion correction was then computed by aligning the frames of a scan to the first frame. The duration of the pRF mapping scans was 372 s (248 time frames), of which the first 12 s (8 time frames) were discarded to eliminate start-up magnetization transients. During the three sessions, we acquired 6 to 8 pRF mapping scans in total per subject. To increase the signal-to-noise ratio, we averaged across the repeated scans. During the three sessions in which we presented the natural images, we acquired 6 to 7 scans for each of the three stimulus sets. The duration of the scans with the natural images was 432 s (288 time frames). The first 12 s (8 time frames) were discarded to eliminate start-up magnetization transients. The images were presented in a block design. Each image was presented during a 9-s block. Within this block, the same image was shown 18 times for a duration of 300 ms followed by 200 ms mean luminance. The full-field stimuli were presented with three alternating different high-contrast patterns, to obtain a full high-contrast response that is not based upon one specific high-contrast pattern (SI Appendix, Fig. S3B). Specifically, the phase of the full-field pattern was randomized on different presentations in order to obtain a response that is not influenced by one specific dartboard pattern. The block in which the stimulus was presented was followed by a 12-s mean luminance presentation. Four longer blank periods of 33 s were also included during the scan.

The T1-weighted MRI images were acquired in a separate session using an 8-channel SENSE head coil. The following parameters were used: TR/TE/flip angle = 9.88/4.59/8. The scans were acquired at a resolution of 0.79 × 0.80 × 0.80 mm and were resampled to a resolution of 1mm3 isotropic. The functional MRI scans were aligned with the anatomical MRI using an automatic alignment technique (51). From the anatomical MRI, white matter was automatically segmented using the FMRIB’s Software Library (FSL) (52). After the automatic segmentation, it was hand-edited to minimize segmentation errors (53). The gray matter was grown from the white matter to form a 4 mm layer surrounding the white matter. A smoothed 3D cortical surface can be rendered by reconstruction of the cortical surface at the border of the white and gray matter (54).

We used bar apertures filled with natural images (18, 21) (SI Appendix, Fig. S3A) to train the pRF model of participants in the fMRI experiment. The width of the bar subtended 1/4th of the stimulus radius (1.375°). Four bar orientations (0°, 45°, 90°, and 135°) and two different step directions for each bar were used, giving a total of 8 bar directions within a given scan. The bar stepped across the stimulus aperture in 20 steps (with a distance of 0.55° and a duration of 1.5 s per bar position) so that each pass took 30 s. A period of 30 s mean luminance (0% contrast) was presented after every pass. In total, there were 4 blocks of mean luminance during each scan, presented at evenly spaced intervals. The participants performed a fixation dot task to make sure they fixated at the center of the display. A small fixation dot (0.11° radius) was presented in the middle of the stimulus. The fixation dot changed its color from red to green at random time intervals and subjects were instructed to respond to color changes using a button press.

The natural images came from the BSD (1, 19). The original resolution of the images was 321 × 481 pixels (both landscape and portrait). In the fMRI experiments (18), we selected a square region of 321 × 321 pixels from the images and upsampled it to a resolution of 516 × 516 pixels, which corresponds to a stimulus of 11 × 11° diameter of visual angle. The images were masked by a circle with a raised cosine faded edge (width of 0.9°), and the areas outside this circle were set to the mean luminance. The images were gamma-linearized and the mean contrast was set to 50%. We used 3 image sets in different scanning runs, each containing 15 different natural images (45 in total) and one full-field binarized bandpass-filtered noise stimulus. SI Appendix, Fig. S2 shows the image set. A fixation dot was presented at the center of the stimulus. We used the same fixation dot task as for the pRF mapping runs.

The pRF model was estimated for every cortical location from the measured fMRI signal that was elicited by the pRF mapping bar stimuli (SI Appendix, Fig. S3A) (18, 21). In short, the method estimates the pRF by combining the measured fMRI time series with the position time course of the visual stimulus. A prediction of the time series is made by calculating the overlap of the pRF and the stimulus energy (RMS contrast, see below) convolved with the hemodynamic response function (HRF). We estimated the parameters of the HRF that best describes the data of the whole acquired fMRI volume (55). The optimal parameters of the pRF model are chosen by minimizing the residual sum of squares between the predicted and the measured time series. We used the conventional pRF model, which consists of a circular symmetric Gaussian. This model has four parameters: position (x, y), size (σ), and amplitude (β). For further technical and implementation details, see ref. 21.

This approach works best in areas with small (p)RFs, such as V1. In higher areas, the probability that pRFs contain a border is higher, because they are larger. Some of the pRFs in higher visual areas contain, for example, a short piece of object border, whereas the remainder was without any border, weakening the contrast. In other words, the distribution of object borders has a lower variance in higher areas, and this reduces the contrast between pRFs with and without object borders.

We used the pRF method to estimate position parameters x and y of the pRF of every voxel. From these values, we derived the polar angle [atan(y₀/x₀)] and eccentricity [√(x₀² + y₀²)] values. We drew the borders between visual field maps on the basis of polar angle and eccentricity maps on the inflated cortical surface (56). We defined visual areas V1, V2, V3, hV4, LO-1/2, and V3-a/b as our regions of interest (ROIs) (57–60).

We measured fMRI responses to 45 natural images (SI Appendix, Fig. S2) and 3 full-field high-contrast stimuli (100% contrast; SI Appendix, Fig. S3B) (18). We first determined the voxel response amplitudes in %BOLD signal change elicited by each of these images. The voxel responses were calculated using a general linear model (GLM) (61, 62). To correct for differences in response amplitudes between voxels, we normalized the responses to the voxel’s response to the full-field (100% contrast) stimulus.

To determine the CRF in the fMRI experiment, we only used the voxels with an overall significant response (t-values > 4.0), a pRF eccentricity between 0.5 and 4° and for which the pRF model explained more than 40% of the variance. Based on previous work, for every area we used a threshold for the pRF sizes (21, 55, 63, 64). In V1, we included pRFs with a value of σ (which determines pRF size) between 0.25° and 0.8°, for V2 between 0.25° and 1.1°, for V3 between 0.25° and 1.75°, for hV4 between 0.45° and 3°, for V3-a/b between 0.45° and 3.75° and for and LO-1/2 between 0.9° and 5°.