Highly accurate retinotopic maps of the physiological blind spot in human visual cortex

Abstract The physiological blind spot is a naturally occurring scotoma corresponding with the optic disc in the retina of each eye. Even during monocular viewing, observers are usually oblivious to the scotoma, in part because the visual system extrapolates information from the surrounding area. Unfortunately, studying this visual field region with neuroimaging has proven difficult, as it occupies only a small part of retinotopic cortex. Here, we used functional magnetic resonance imaging and a novel data‐driven method for mapping the retinotopic organization in and around the blind spot representation in V1. Our approach allowed for highly accurate reconstructions of the extent of an observer’s blind spot, and out‐performed conventional model‐based analyses. This method opens exciting opportunities to study the plasticity of receptive fields after visual field loss, and our data add to evidence suggesting that the neural circuitry responsible for impressions of perceptual completion across the physiological blind spot most likely involves regions of extrastriate cortex—beyond V1.


| INTRODUCTION
The visual cortex in the occipital lobe contains precise topographic maps, with neighboring locations of the observer's visual field encoded by adjacent populations of neurons. Damage at various stages of the primary visual pathway results in localized vision loss, such as retinal damage in age-related macular degeneration and glaucoma, as well as stroke or brain trauma that can impact structures at later stages of the human visual hierarchy. These regions of localized loss of visual sensitivity from damage, known as scotomas, are not entirely irreversible, as there is evidence for plasticity by which the brain reorganizes mappings between inputs and neural responses in regions surrounding the region of acquired blindness-either through prolonged experience or active perceptual learning (Barbot et al., 2021;Gilbert, 1993;Miller et al., 2015;Saionz et al., 2021).
Such plasticity may enable the visual system to compensate for some vision loss, by making better use of intact neural circuitry.
Functional magnetic resonance imaging (fMRI) theoretically enables researchers to study neural processing and map reorganization around the scotomas of human participants in a non-invasive fashion. Advances in neuroimaging technology and data analysis methods have greatly improved the knowledge we can glean about neural functioning in visual cortex. For over a decade, population receptive field (pRF) analysis has become an important method in the toolkit of visual neuroscientists (Dumoulin & Wandell, 2008;Wandell & Winawer, 2015). The pRF estimates the position, size, and shape of receptive fields of individual fMRI voxels, corresponding to the aggregate activity of all the neurons within the voxel. Several studies have used pRF analysis to map the extent of visual scotomas and compare them to behavioral measures (Barbot et al., 2021;Papanikolaou et al., 2014;Prabhakaran et al., 2021;Ritter et al., 2019;Silson et al., 2018). However, as a means for better understanding scotomas in patients this comes with several caveats, not least of which is that ground truth is unknown. Functional brain maps of a scotoma patient can only be compared to maps in healthy controls, or possibly to equivalent maps at an intact location of the visual field.
The physiological blind spot is a region of the visual field of each eye that corresponds with the eye's optic disc. There are no photoreceptors in this region of retina, which causes a naturally occurring scotoma. This seems to make it an ideal model for studying visual scotomas in healthy control participants; however, to date very few investigations of the neural representation of the human blind spot have been attempted. Post-mortem studies of ocular dominance patterns in human area V1 (Adams et al., 2007) have demonstrated that the physiological blind spot corresponds with a small monocular oval region of V1, subtending a cortical surface area of approximately 50 mm 2 . It is located at an eccentricity of approximately 15 , where cortical magnification is already much reduced compared to the fovea or parafovea . This means that even with high spatial resolution brain imaging, only a small number of voxels encode responses from this region of V1. Nonetheless, fMRI experiments have shown that the blind spot is a discontinuity in the visual field representation: with responses to spatially separate stimuli on either side of a blind spot exciting anatomically segregated responses in V1, despite observers perceiving such a configuration as single and continuous (Awater et al., 2005;de Hollander et al., 2020).
While some important preliminary observations have been made regarding representations of the physiological blind spot in V1, to date, no detailed retinotopic maps of this region exist. Here, we report the world's first detailed mappings of retinotopic organization, in and around representations of the physiological blind spot in V1 for nine individuals. Our novel data-driven method will provide researchers with the means for further investigations of this region. More generally, it enables investigations of retinotopic organization around regions of acquired localized blindness which, like the physiological blind spot, are often located in the periphery of human vision-and consequently can pose similar challenges to brain imaging, that we have now solved. Moreover, our data provide converging evidence, suggesting that neural circuits responsible for perceptual completions of form across physiological blind spots most likely have a substrate in extrastriate cortex-beyond V1.

| Participants
Seven observers (ages: 21-41, 3 Female, all right-handed) participated in experiments using a 3 Tesla scanner in the Centre for Advanced MRI at the University of Auckland, New Zealand. Two additional observers (ages: 30 and 37, 1 female, all right-handed) participated in experiments on a 7 Tesla scanner at the Centre for Advanced Imaging at the University of Queensland, Brisbane, Australia. All observers were recruited from staff and student populations at each site and gave written informed consent to take part. Experimental procedures were approved by the University of Auckland Human Research Participants Ethics Committee, as well as the ethics committee at the University of Queensland. As is common in retinotopic mapping studies, we present datasets from each participant as individual replications.
The purpose of including the two exemplar 7 Tesla datasets was to generalize the procedure to a different scanner setup with an improved spatial resolution. We did not make any quantitative comparisons between data from the two scanning sites as these were deliberately not matched in terms of experimental parameters.

| Procedure
All participants wore an ophthalmic eye patch over one eye (4 observers: left eye; 3 observers: right eye; both observers in 7 Tesla scans: left eye) before they were placed inside the scanner bore. After initializing the scan, we used a behavioral localizer (see details below) to determine the position and borders of the participant's physiological blind spot. The center location of the blind spot thus determined was then used to place the stimuli in the retinotopic mapping experiments. Participants underwent 4-6 runs, each lasting 4 min 10 s on 3 Tesla or 6 min 6.5 s on 7 Tesla, for mapping pRFs in the region in and around their blind spot. Following that, the patient bed was moved some ways out of the bore, but the participant remained on it. We moved the eye patch to cover the opposite eye and returned the participant into the bore. The participant then underwent an identical number of runs of pRF mapping, for the same visual field region but now exciting the control eye without a blind spot. Finally, we again moved out the patient bed to remove the eye patch and added the front visor to the 32-channel head coil (at Auckland site only; see scanning parameters). We then acquired a T1-weighted scan of the participant's brain anatomy. For the first participant at the Brisbane site, we found that the scanner setup restricted their field of view when the control (right) eye was patched. We therefore used the behavioral localizer to map out the visible portion of the stimulus in the control scan for this participant. In the second participant, we instead collected a binocular viewing control (i.e., no eye was patched) to ensure they could see the whole stimulus. All stimuli were generated and displayed in MATLAB using Psychtoolbox 3 (Brainard, 1997).

| Stimuli
At the Auckland site, stimuli were presented via an LCD screen (BOLD screen, Cambridge Research Systems, Rochester, U.K.) placed at the back of the scanner bore. Participants viewed the screen via a mirror mounted on the head coil at a viewing distance of 111 cm.
Thus, the whole screen subtended a visual angle of 35.5 horizontally and 19.9 vertically, although the top corners of the screen were obscured by the curvature of the bore.
At the Brisbane site, stimuli were presented via a LCD laser projector (Sony VPL-FHZ55) back-projected onto a semi-opaque screen placed at the back of the scanner bore. Participants viewed the screen similarly to that at the Auckland site, but with a total viewing distance of 120 cm. Here the whole screen subtended a visual angle of 22.6 horizontally and 12.7 vertically, also with the top corners of the screen obscured.

| Blind spot localizer
When the participant had been placed inside the scanner bore, but prior to the functional scans, we localized the physiological blind spot in each observer. The participant fixated a small (diameter: 0.17 of visual angle) black dot located 10.6 of visual angle from the screen center on the side ipsilateral to their patched eye. A radar grid pattern comprising radial lines and concentric circles Morgan & Schwarzkopf, 2019) that emanated from the fixation dot was also displayed to enhance fixation compliance. A larger red disc (diameter: 0.53 ) was then moved by the experimenter using the computer mouse. The participant responded verbally via the intercom whenever this red target vanished or reappeared. The experimenter marked this location with a left click on the mouse and the program then placed a small grey dot (diameter: 0.17 ) at that location. We moved the target in several directions to ascertain the boundaries of the blind spot, also moving it out of the approximate region to minimize response biases. Occasionally, the experimenter made multiple estimates of the boundary at the same location to improve accuracy.
Once the border of the blind spot had been estimated with sufficient detail, a right click on the mouse ended this phase. At this point, geometrically redundant points were removed from the estimated blind spot boundary. The program calculated the centroid location of these border dots and displayed a black dot (diameter: 0.17 ) at that location, as well as a line circle (diameter: 10.6 ) around it, marking the visual field region to be stimulated during retinotopic mapping. We debriefed the participant to report what they could see. An accurate estimate of the blind spot entails that none of the small dots would be visible, but the large line circle should be clearly visible (meaning it was sufficiently far from the edge of the blind spot). The procedure would be repeated if these criteria were not met although this did not occur.

| Retinotopic mapping
Retinotopic mapping stimuli were like those used in previous studies (Alvarez et al., 2015;Dumoulin & Wandell, 2008 (b) Analysis pipeline for reverse correlation pRF fit. Binary apertures indicating the stimulus location relative to the blind spot center at each time point of the experiment are convolved with a canonical hemodynamic response function. To determine the pRF for a given voxel, its actual measured time series is correlated with each pixel time series in convolved apertures. This results in a reverse correlation pRF profile. We then fit a two-dimensional Gaussian pRF model to this profile, to estimate the location (x, y) and size (σ) of the pRF changed length at each step. One sweep of the bar along a given direction lasted 25 s, with one discrete 0.4 step of the bar per second. Each sweep was along the axis perpendicular to the bar orientation. We used eight orientations/directions in sequence starting from horizontal moving upward (0 ), and then changing in steps of 45 until we reached the final orientation of 315 . Interspersed between the fourth and fifth sweep and after the eighth sweep, we presented a blank period during which only the fixation dot and radar grid were presented. Inside the bar, a checkerboard (side length: 0.95 ) was flashed with 6 Hz.
To ensure that participants maintained fixation, they performed a detection task. The duration of the scan was divided into brief 200 ms epochs. In each, there was a 10% probability that the small black fixation dot would change into a black grapheme. This could be any of the 26 letters of the English alphabet, or a digit from 0-9. Participants were instructed to press a button whenever they saw a number. The 200 ms epoch following a character/digit was always followed by the fixation dot to ensure such events were not too fast.

| Magnetic resonance imaging
At the Auckland site, we used a Siemens SKYRA 3 Tesla scanner with a 32-channel head coil where the front element had been removed to permit an unrestricted view of the screen. This setup results in 20 effective channels covering the back and the sides of the head.
Between 4-6 pRF mapping runs of 250 T2*-weighted image volumes were acquired for each eye. We used an accelerated multiband sequence with a TR of 1000 ms and 2.3 mm isotropic voxel resolution, with 36 transverse slices angled to be approximately parallel to the calcarine sulcus. The scan had a TE of 30 ms, flip angle of 62 , field of view 96x96, a multiband/slice acceleration factor of 3, an inplane/parallel imaging acceleration factor of 2, and rBW was 1680 Hz/Px. After acquiring the functional data, the front portion of the coil was put back on to ensure maximal signal-to-noise levels for collecting a structural scan (a T1-weighted anatomical magnetizationprepared rapid acquisition with gradient echo [MPRAGE] scan with a 1 mm isotropic voxel size and full brain coverage).

| Data preprocessing
Functional data were realigned and co-registered to the anatomical scan using default parameters in SPM12 (https://www.fil.ion.ucl.ac. uk/spm). We further reconstructed and inflated surface mesh models of the grey-white matter boundary Fischl et al., 1999;Fischl, 2012) using the automatic reconstruction algorithm in FreeSurfer (Version 7.1.1; https://surfer.nmr.mgh.harvard. edu). Functional data were then projected onto this cortical surface model: for each vertex in the surface mesh, we determined the voxel in the functional image that lay at the midpoint between the vertex on the grey-white matter boundary, and the same vertex on the pial surface boundary. We then applied linear detrending to time series for each vertex and run to remove slow drifts, and the time series were normalized to z-scores. We then averaged the runs for the blind spot and control eye, respectively, resulting in two runs of 250 volumes for each condition.
We restricted all further analyses approximately to the occipital lobe by choosing vertices in the inflated surface model whose ycoordinates in FreeSurfer space ≤ À 35. We calculated a noise ceiling as follows: first, we separately averaged odd and even runs from each condition, and then correlated the time series for each half for each vertex. This gives an estimate of test-retest reliability of the visual response; however, this is not an accurate estimate of the true reliability of the time series because it is only based on half the data.
Using the Spearman-Brown prophecy formula (Spearman, 1910), we can extrapolate the actual reliability, r max , that is, the maximum theoretically achievable correlation for each vertex.
where r denotes the correlation between time series for odd and even runs. However, since goodness-of-fit is expressed by the coefficient of determination, R 2 , to obtain the noise ceiling, we need to calculate the square of this measure: This noise ceiling is the maximally achievable R 2 for a given time series and directly related to whether vertices are visually responsive. Since we only mapped a small fraction of the peripheral visual field, we only expect responses in a small part of visual cortex. For expediency, we therefore further limited our analyses to only those vertices with a noise ceiling exceeding 0.15, as it simply does not make sense to conduct the pRF modelling on vertices without an appreciable visual response. This threshold is arbitrary but chosen as a good middle ground since it revealed consistent and contiguous clusters of activation in V1 without an excessive number of stray vertices outside the parts of cortex expected to respond to our stimuli. Using this preselection is particularly useful for the slow forward-modelling analysis approach as it speeds up the analysis by a factor of 3-5.

| Reverse correlation pRF analysis
In all following analyses, we treated the center of the blind spot as standard pRF mapping analyses would treat the center of gaze. That is, we define the blind spot center as the center of the mapped visual field. The circular visual field region surrounding the blind spot center had a radius of 5.3 . We created a sequence of stimulus apertures, which depict the location of the bar stimulus for each 1 s fMRI volume with a binary 100 Â 100 pixel matrix, where each pixel denotes whether a stimulus was present at that location, or not . We estimated the pRF location by determining the profile maximum. We also recorded the squared correlation R 2 RC between this pixel and the fMRI response.
Resulting pRF profiles can be asymmetric or have multiple peaks.
Given our use of a regular bar stimulus, profiles can also contain artifacts related to spatiotemporal correlations within the stimulus sequence. Consequently, we restricted further analyses to vertices where R 2 RC >0.1. The distribution of R 2 RC values ( Figure S1a) suggests that this threshold removes the overwhelming majority of vertices that do not contain distinct reverse correlation profiles.
We then fit a two-dimensional pRF model to these reverse correlation profiles, to estimate the location and extent of pRFs. We used a symmetric two-dimensional Gaussian profile (Figure 1b, bottom-right), defined by its x and y positions, its SD σ, and amplitude β. The goodness of pRF fit (R 2 pRF ) is determined by calculating the sum of squared pixel-wise residuals between the predicted pRF profile and the one estimated from reverse correlation. All further quantitative analyses of data are restricted to vertices where R 2 pRF >0.5. The distribution of these values ( Figure S1b) shows that this threshold includes the majority of vertices that passed the R 2 RC threshold in the first analysis step. Vertices with R 2 pRF below threshold contain only fuzzy, illdefined reverse correlation profiles than cannot be fit well with a 2D Gaussian model.

| Forward-model pRF analysis
We also conducted a standard pRF analysis using forward-modelling (Dumoulin & Wandell, 2008;Morgan & Schwarzkopf, 2019;Moutsiana et al., 2016;Schwarzkopf, 2018). This analysis fits a predicted time series of the fMRI response to the observed time series, based on known stimulus apertures and an assumed pRF model profile. In short, a neural response prediction is generated by overlaying the stimulus aperture onto a two-dimensional Gaussian pRF model. This predicted response time series is further convolved with the canonical hemodynamic response. The best fitting pRF model is determined by varying its visual field position (x, y) and size (σ) parameters, and calculating the correlation between the predicted and observed time series. This fitting uses a coarse-to-fine procedure: first, we generate several thousand predictions for a plausible range of permutations of the three pRF parameters and calculate an extensive grid search by finding the best-correlated combination. These parameters are then used to seed an optimization procedure (Lagarias et al., 1998;Nelder & Mead, 1965) to fine-tune the fit. Because correlation is scale-invariant, we also calculated a linear regression between this fit and the observed time series, to determine the amplitude (β 1 ) and baseline level (β 0 ) of the response. These pRF modelling procedures have been described in detail elsewhere (Morgan & Schwarzkopf, 2019;Moutsiana et al., 2016). However, here, we used a different biophysical model for predicting a pRF's neural response than in previous work by us and others. Specifically, we calculated the percentage of overlap between the pRF model profile and the stimulus aperture. This accounts for the fact that the response to a constant visual stimulus should vary with different pRF sizes. Note that this merely affects the modelled signal amplitude, not the pRF position or shape.
Data from the forward modelling were statistically thresholded as follows: we normalized the raw goodness-of-fit (R 2 ) by dividing it by the noise ceiling. This measure, nR 2 , can be intuitively described as the proportion of explainable variance explained by the pRF model. We selected only those vertices with nR 2 > 0.1. Note that this is not numerically comparable to the goodness-of-fit values used in the reverse correlation analysis. However, this threshold is extremely liberal as it includes the vast majority of vertices with any pRF fit ( Figure S1c).
Our main analysis used the stimulus apertures of the full bar stimuli as they would have appeared on the screen (and which were used for the reverse correlation analysis). Previous research has shown, however, that pRF estimates derived via forward-modelling are susceptible to potential biases when a scotoma is not explicitly modelled in the stimulus apertures (Binda et al., 2013). This could produce spurious estimates of pRFs inside the blind region. We therefore generated masked stimulus apertures, specific for each individual participant, based on the blind spot localizer. Within the polygon described by the blind spot border, the aperture was set to zero. This mask was smoothed with a Gaussian kernel of 0.4 , to make the edges of the scotoma more biologically plausible. We then repeated the forward-model pRF analysis using this masked aperture.

| Regions of interest
All quantitative analyses of the blind spot were restricted to a probabilistic estimate of V1, based on an anatomical atlas derived from postmortem brains (Hinds et al., 2008). We selected a continuous region of this atlas, predicted to be within V1 with ≥80%. In all participants, this revealed clusters of significant pRF fits (especially for the reverse correlation approach) at a location consistent with the expected representation of the blind spot in V1 (Adams et al., 2007;Awater et al., 2005). However, due to artifacts, a few stray vertices occasionally survived thresholding throughout the rest of V1. We therefore further defined a circular region of interest (ROI) around the main cluster. Using the noise ceiling map for the control condition, we manually selected the vertex approximately at the center of the cluster and then used a geodesic region growing procedure by taking 12 steps across the grey-white matter surface mesh away from that center. This ROI encompassed the significant cluster. All further analyses reported were conducted exclusively on vertices from this ROI.

| Blind spot localization
Prior to scanning, we used a behavioral localizer to delineate the borders and determine the center location of the blind spot in the visual field. Averaged across all nine participants, the blind spot center was 15.8 ± 1.5 (SD) displaced from the vertical meridian and 0.8 ± 0.6 below the horizontal meridian. This translates to an average eccentricity of 15.9 ± 1.5 . Its average width and height were 5.2 ± 0.7 and 6.1 ± 0.6 , respectively, with an average area of 24.1 ± 6.4 2 . Individual blind spot dimensions are listed in Table S1. Maps for distance from the blind spot center were less clear. Nevertheless, for the control eye, a concentric organization is visible (Figure 2b, top). In the blind spot condition, not much difference could be determined between vertices (Figure 2b, bottom). This is unsurprising, however, as the blind spot effectively removes central vertices from these clusters. Notwithstanding this, differences between blind spot and control maps reveals a more complete coverage of the mapped region in the control condition. Taken together, these maps demonstrate that our experiments activated the expected clusters of V1 vertices, and that despite their small size, these clusters contained retinotopic maps consistent with the corresponding visual field location.

| Retinotopic maps around the blind spot
To obtain more detailed maps and to generalize the applicability of this method to a different experimental setup, we carried out similar experiments in two participants at 7 Tesla using a finer voxel resolution (1.5 mm isotropic voxels). Indeed, results showed an even clearer pinwheel progression of radial coordinates (Figure 3a), and a distinct concentric gradient in pRF distances from the blind spot center in control scans (Figure 3b Next, we quantified the gradient of pRF size when moving from the foveal to the peripheral side of the blind spot. We selected all vertices with R 2 pRF >0.5 whose horizontal visual field position was within 6 of the blind spot center. Separately for each participant, and each condition, we fit a robust linear regression to estimate the slope of the gradient. To make data comparable across participants whose left and right hemispheres were scanned, we inverted the sign for hori-  (Hinds et al., 2008) significantly different (t[8] = 1.56, p = .15763). This demonstrates that in both conditions, pRF size increased from more central to peripheral locations, even when that gradient was not always obvious in all cortical maps obtained at 3 Tesla. Of note, pRF sizes from the 7 Tesla scans (red and brown curves) were smaller than for the other participants, although due to the small sample size, we refrain from drawing statistical conclusions about this observation. But tentatively, this difference in pRF size could be due to the smaller voxel sizes used at 7 Tesla, consistent with previous reports (Himmelberg et al., 2021).

| Visual field coverage of reverse correlation profiles
We next used the reverse correlation profiles of pRFs to reconstruct how much of the visual field was covered by these pRFs. Peak amplitudes of the profiles from all significant vertices were normalized to their maximum. This is standard practice in visual field coverage plots because they seek to show how the visual field is encoded by the pRF mosaic, irrespective of the responsivity of individual F I G U R E 4 Quantifying pRF size gradient across the visual field. pRF sizes from the control eye stimulation (a) or blind spot eye stimulation (b) were plotted against the horizontal pRF location relative to the blind spot center (negative and positive values, respectively, denote visual field locations on the foveal side and peripheral side of the blind spot). Each dot denotes parameters from one pRF. Solid lines indicate a robust linear regression fit. Colors indicate the 9 different participants. Red and brown data are from the participants scanned at 7 Tesla pRFs. However, weighting each pRF by its peak amplitude had only negligible effect on the appearance of these reconstructions (data not shown). Moreover, all values below 50% of the maximum were set to zero. This effectively restricts the pRF profile to the fullwidth-at-half maximum and was designed to minimize the influence of any spatiotemporal correlations in the stimulus sequence. Without this clipping step, the visual field reconstructions appeared somewhat smoother but they were qualitatively very similar to the main analysis (data not shown). We then averaged these restricted profiles and plotted visual field coverage. In the control condition, the whole mapped region of the visual field is covered by pRFs, albeit not homogeneously (Figure 5a). In particular, the density of pRFs was greater on the side facing the fovea than the periphery, consistent with the gradient of cortical magnification and receptive field density in V1 . We obtained similar maps for the blind spot condition (Figure 5b), but notably these receptive field profiles spared the locations of the blind spot scotoma (as determined by the behavioral localizer) in all participants.
Due to the inhomogeneity in coverage, the annulus regions around the blind spot contains some gaps. Nevertheless, this suggests the reverse correlation profiles precisely captured the extent of blindness.

| Visual field reconstruction of pRF fits
To visualize the position of pRFs directly, we next created scatter plots of pRF position and size (defined by a radius of 1σ) in visual space ( Figure 6). Consistent with the visual field coverage (Figure 5), this showed complete coverage of the mapped visual field region in the control condition (Figure 6a). In the blind spot condition, pRF centers spared the blind spot scotoma (Figure 6b), although many Gaussian pRF fits extend into the scotoma. As suggested by the reverse correlation profiles, the density of pRFs in the annulus region surrounding the blind spot varied. Nevertheless, the entire region outside the scotoma was covered by pRFs. To quantify the visual field coverage, we divided pRFs from the control condition into those whose centers fell inside the scotoma, and those that fell outside of it. Figure 7 plots the mean peak amplitudes for pRFs outside the scotoma against pRFs inside the scotoma for each participant. Data from the control condition (blue circles) clustered around the identity line with half of data points above and half below it. This shows that responses were comparable for pRFs inside and outside the scotoma. In contrast, data from the blind spot condition suggested weaker responses for pRFs inside the scotoma than outside it; only one participant had a data point above the identity line. It is, however, important to stress that using the control condition to define this region of interest incorporates not only the true pRF position but also any error inherent in mapping the blind spot's location and extent. Some pRFs estimated as falling inside the scotoma might truly be located outside of it, and vice versa. Hence, this analysis necessarily underestimates differences between these regions.
Please note that data from the first participant scanned at 7 Tesla are not included in this plot, as the incomplete field of view in the control condition (when the opposite control eye was patched) precluded us from carrying out this analysis. This incomplete field of view, however, afforded us with a serendipitous opportunity to further test the accuracy of this approach for mapping visual field scotomas. For this participant, we show the visual field coverage and pRF scatter plot separately in Figure 8. As with other participants, the mapping approach allowed precise localization of significantly activated pRFs sparing the blind spot scotoma (Figure 8c,d). Critically, in the control condition, these reconstructions also distinctly respected the boundary of the visible portion of the visual field (see green/black dots in Figure 8a,b). This demonstrates that the approach could also map the broader absence of stimulation in these scans. The method should generalize to other kinds of visual field loss than that due to the physiological blind spot.

| Comparison with forward-model pRF fits
We also analyzed our data using a more traditional pRF analysis based on forward-modelling Dumoulin & Wandell, 2008;Morgan & Schwarzkopf, 2019). Figure 9 shows scatter plots of pRFs analyzed using this approach, comparable to Figure 6.
As with our reverse correlation data, in the control condition, pRFs covered most of the mapped visual field region, albeit less completely ( Figure 9a). Moreover, in the blind spot condition, pRF centers exclusively spared the blind spot scotoma (Figure 9b). However, only a small number of vertices survived statistical thresholding, despite the F I G U R E 8 Reconstruction in visual space (a,c) and corresponding pRF model fits to reverse correlation profiles (b,d). The scatter plots denote the location and size of pRFs. All conventions as in Figures 4 and 5, respectively. Data is from a participant scanned at 7 Tesla who had an incomplete field of view (FOV) with their control eye. (a,b). Control eye stimulation. (c,d) Blind spot stimulation. Green/black dots denote the visible field of view (a,b) or the outline of the blind spot (c,d) as determined by the behavioral localizer fact that the threshold for this analysis was very liberal ( Figure S1c). This renders these maps generally sparse and incomplete. The standard forward-modelling approach was evidently less effective in estimating the extent of the scotoma. Reconstructions of pRF profiles in visual space were unexpectedly also blurry and ill-defined and, in the case of the blind spot condition, very sparse ( Figure S2).
To further quantify the reliability of pRF positions for the two analysis methods, we reasoned that for pRFs located outside the blind spot, the position estimates should be similar for the control and blind spot condition. We therefore selected the pRFs that survived statistical thresholding in both conditions, and then computed the mean Euclidian distance between the blind spot and control estimates for these pRFs, separately for the reverse correlation and the forwardmodel analysis. This analysis revealed significantly smaller distances in pRF positions (t[7] = À3.6, p = .0087) for the reverse correlation (mean = 1.66 ) than the forward-model (mean = 2.59 ).
Previous research suggested that when using standard pRF modelling techniques for mapping scotomas, it is important that the stimulus model incorporates the extent of the scotoma (Binda et al., 2013). Without this, the modelling approach can produce biased estimates of pRF position. Most notably, this should result in pRFs being falsely estimated to fall inside the scotoma, something we did not see in our data at all (Figure 9b). We nevertheless conducted a control analysis, using stimulus apertures for each participant in which the blind spot scotoma had been masked out. The results for these pRF maps were extremely similar to those obtained with the complete stimulus apertures (Figure 9c). This demonstrates that the poor quality of forward-model pRF estimates in the blind spot condition was not simply due to using an incorrect stimulus model.

| DISCUSSION
We used pRF analysis based on reverse correlation to map the retinotopic organization of the physiological blind spot in human observers.
We obtained highly accurate maps of the visual space surrounding the blind spot. In particular, the organization of radial (angular) pRF position relative to the blind spot center followed the pattern expected based on the known retinotopic organization of V1. Maps for the distance from the blind spot center, and especially pRF size, were less clear in participants scanned at 3 Tesla. This is likely due to the comparably coarse voxel resolution relative to the size of the blind spot.
Two participants scanned at 7 Tesla, with 1.5 mm isotropic voxel resolution, confirmed expected map gradients for these pRF parameters. We further projected pRFs back into visual space to plot the visual field coverage of these small retinotopic maps. This showed accurate delineation of brain activity to only those sub-regions of the visual field outside the scotoma. In contrast, a conventional forward-modelling approach for estimating pRFs was far less accurate. While pRFs estimated in that way also consistently spared the blind spot scotoma, they were sparse and afforded incomplete visual field coverage. Control analyses confirmed this was not trivially explained by failing to account for the scotoma in the pRF model (Binda et al., 2013).
These findings do not imply that the reverse correlation approach is categorically superior to forward-modelling in all scenarios. The latter approach has numerous advantages that make it the ideal choice for addressing some research questions. Unlike reverse correlation, the forward-modelling approach can fit pRFs located outside the stimulated part of the visual field. While such fits must be treated with caution, this can nevertheless be a great advantage. Forwardmodelling is probably also less susceptible to spatiotemporal correlations in the stimulus sequence (although see [Alvarez et al., 2015;Infanti & Schwarzkopf, 2020;Linhardt et al., 2021]). Moreover, forward-modelling is based on theoretical models for the shape and function of pRFs. This allows researchers to test specific models of sensory processing, such as divisive normalization (Aqil et al., 2021).
Forward-modelling is also better suited for sparse stimulus designs, such as might be used in one-dimensional tuning models, such as numerosity tuning (Harvey et al., 2013), or for studying somatosensory cortex (Puckett et al., 2020). However, in the context of measuring scotomas, especially small, peripheral ones like the physiological blind spot, our present results suggest the reverse correlation approach is the more optimal method.

| Mapping scotomas
The study of retinotopic maps and reorganization around regions of vision loss is hampered by several methodological limitations. Testing and validating the techniques for visual field mapping require a welldefined model of how the scotoma should behave. Some researchers have therefore used masks to block out pre-defined parts of the visual field to validate measurements (Binda et al., 2013;Hummer et al., 2017). These efforts are important, and have resulted in several notable insights-for example, demonstrating that estimates of preferred locations near scotomas can be biased, especially when the scotoma is not factored into the analysis (Binda et al., 2013). Yet masking the stimulus by replacing it with a blank region on the screen is not the same as an actual scotoma. This approach can therefore only inform about analytical issues, and not about any physiological peculiarities regarding the consequences of visual stimulation within and immediately surrounding the scotoma.
Other studies have investigated cortical representations of the foveal rod scotoma (Barton & Brewer, 2015;Baseler et al., 2002), which arises because macula contains predominantly cone photoreceptors. However, such experiments require scotopic stimulation (ideally after dark adaptation), rendering them difficult to apply broadly and limiting the generality of these results.
Another recent study used an elegant design to map artificial scotomas-a visual illusion in which a dynamically flickering background fills in a static region of the stimulus, effectively simulating temporary blindness within that region (Carvalho et al., 2021). While fascinating, this is probably not directly comparable to an actual scotoma, arising due to an absence of innervation/stimulation. Differences in terms of physiological processes might be considerable (Weil et al., 2007;Weil et al., 2008)-notwithstanding the conceptual interest of that study. In contrast, the physiological blind spot is a visual field region corresponding to the optic disc where there are no retinal photoreceptors. This makes it a perfect model for a small, localized scotoma, as might appear after retinal damage or lesions to primary visual cortex.
Our results demonstrate that even at relatively conventional voxel resolutions, it is possible to produce accurate maps of visual scotomas. While several studies investigated the potential of pRF mapping for estimating scotomas (Barbot et al., 2021;Papanikolaou et al., 2014;Prabhakaran et al., 2021;Ritter et al., 2019;Silson et al., 2018), our method can reconstruct even small, peripheral scotomas like the physiological blind spot with unprecedented precision.
This opens a promising avenue for studying map reorganization and plasticity after damage to the visual pathway (Barbot et al., 2021;Gilbert, 1993;Papanikolaou et al., 2014;Saionz et al., 2021). Previous research has found no evidence for large-scale changes in retinotopic maps after macular degeneration (Baseler et al., 2011). However, previous experiments might have been susceptible to limitations of previous analysis procedures (Binda et al., 2013), or plasticity effects might have been too subtle to be resolved by those experiments. The precision with which our reverse correlation maps could delineate the borders of the blind spot scotoma suggests it could be a more powerful tool for exploring small-scale changes.

| Future directions
Further refinements to our method could be considered. Using reverse correlation for mapping pRFs is not new. Rather, this technique has been used in a number of studies (van Es et al., 2018;de Haas et al., 2014;Lee et al., 2013;Ress et al., 2011). Some of these used analysis approaches (e.g. using ridge regression) that account for the considerable spatiotemporal correlation in conventional pRF stimulation paradigms (Alvarez et al., 2015;Infanti & Schwarzkopf, 2020).
We decided against this, as these modelling approaches drastically increase processing time. Despite spatiotemporal correlations in the stimulus sequence, our approach evidently produced precise maps of the blind spot scotoma. However, future studies might consider using such refinements, for example, when studying subtle changes in pRF parameters around the blind spot. Randomized stimulation sequences that break spatiotemporal correlations might further enhance the accuracy of these maps, although they are also associated with poorer signal-to-noise ratios (Binda et al., 2013;Infanti & Schwarzkopf, 2020;Ma et al., 2013).
It is also possible to fit more complex pRF models in the second step than the two-dimensional Gaussian we used here. For example, pRFs at the edge of the scotoma may be asymmetric or skewed. Such an analysis would again benefit from a stimulus design that minimizes spatiotemporal correlations. We also did not attempt fitting such asymmetric pRFs because we have no clear hypothesis for what shape these should have. This would increase the chance of overfitting noisy data. Many pRFs, especially those whose true centers fall outside the stimulated part of the visual field, will also produce artifactually asymmetric profiles.
Due to logistical limitations, we carried out our experiments at the two sites independently. This precluded any direct comparison of the 3 Tesla and 7 Tesla data in the same participants. In previous work, we showed that pRFs estimated through forward-modelling are very similar between 1.5 Tesla and 3 Tesla located in London, U.K., and Auckland, New Zealand, respectively (Morgan & Schwarzkopf, 2019). Similarly, another study showed similar results for ocular dominance columns in V1 using two different 7 Tesla scanners located in Amsterdam, Netherlands and Beijing, China. Based on these findings, we would expect measurements of the blind spot to be reliable across field strengths-but of course this hypothesis remains to be tested.

| Perceptual filling-in
Our method also opens opportunities for studying active perceptual processing. Under normal binocular viewing conditions, input from the fellow eye compensates for the region of blindness resulting from the physiological blind spot, and even with monocular viewing the visual system can fill in missing information by perceptually extrapolating patterns from the regions directly abutting the blind spot. Maps obtained using our method could be used to study the neural correlates of this perceptual completion, at a level of detail thus far reserved for invasive electrophysiological recordings in animal models (Komatsu et al., 2000;Komatsu, 2006;Matsumoto & Komatsu, 2005).
While we made no attempt to quantify the effect, our traversing bar stimulus presumably induced at least some perceptual completion, in that none of our participants reported that the moving bars used in mapping seemed to split in two as they passed through the blind spot.
It is therefore notable that no responses were observed in regions of V1 corresponding with the site of the scotoma. This is consistent with previous fMRI experiments (Awater et al., 2005), which similarly found no representation of the physiological blind spot in V1. The convergence of this evidence suggests that perceptual filling-in likely arises through circuits involving extrastriate brain regions-beyond V1. We make no attempt here to test this because a systematic investigation of the neural correlates of filling-in requires two things: a stimulus design that reliably and strongly induces filling-in, and a behavioral task to confirm that.