The inferred retinotopic maps, predicted by Bayesian inference, provide high-quality descriptions of the full retinotopic topology for each subject’s V1-V3 regions. These maps can be produced even in the absence of observed retinotopic measurements (i.e., using the prior alone), or by combination with retinotopic data. In this latter case, the inferred maps have all the advantages of the retinotopic prior (topologically smooth maps, predictions beyond the stimulus aperture, etc.), while also accounting for idiosyncrasies in individual maps. Three examples of maps that demonstrate this advantage are shown in Figure 3. The first two columns show maps in which, relative to the validation maps, the predictions made from data alone have highly curved iso-eccentricity contours. These contours reflect noise rather than true curvature in the iso-eccentricity contours, as shown by the validation data. For these two columns, the predictions from the prior alone have iso-eccentric contours that are too smooth (as compared to the validation data). The correct lines appear to lie between the training data and the prior. Hence when data is combined with the prior (Figure 3, third column) the iso-eccentric contours resemble those of the validation dataset. The third row of Figure 3 shows an instance in which, even lacking a coherent polar angle reversal to define the ventral V1/V2 boundary near the fovea in the predictions made from data alone, combination of the data with prior more accurately predicts that boundary in the validation dataset than the prior alone.

Figure 3

Download asset Open asset

In constructing the Bayesian-inferred retinotopic map for a single hemisphere, we perform two deformations, detailed in Figure 4 (see also Materials and methods): (1) we first deform that hemisphere’s inflated surface to register it to an average anatomical atlas using FreeSurfer (leftmost arrow in Figure 4C); and (2) we then further deform the surface to register it to the retinotopic prior, based on the hemisphere’s retinotopic measurements (second arrow in Figure 4C). These steps together account for the individual differences in the organization of the subjects’ retinotopic maps. In step 1, we account for structural differences across subjects—the deformation that occurs for this registration is unique for each subject. This is where prior work ended (Benson et al., 2012; Benson et al., 2014). In step 2, we account for the differences in the relationship between structure and function across subjects. Although it is possible that deformations in step two partly compensate for imperfections in the first two registrations, we propose that, to a first approximation, the deformations applied in the last step indicate meaningful individual differences in the structure-function relationship.

Figure 4

Download asset Open asset

The Bayesian model fitting allows us to parcellate two sources of variation between individuals: differences in surface anatomy (sulcal patterns) and differences in structure-to-function mapping, that is how retinotopic features map onto the surface anatomy. These two sources of variation map approximately to the two deformations in our atlas fitting: differences in surface anatomy are reflected in the deformation used for the surface alignment, and differences in structure-to-function mapping are reflected in the deformation for retinotopic alignment (the Bayesian fit). Note that in our implementation, both alignments are achieved by warping the individual subject’s vertices, the former to minimize error in surface curvature, the latter to minimize errors in retinotopic measures. (After the process is complete, the alignments can be reversed, thereby bringing the retinotopic predictions back into the native anatomical space.) Because both the anatomical alignment and retinotopic alignment are computed as changes in the position of surface vertices, it is straightforward to compare the two processes.

There are some subjects for whom there are large differences between the retinotopic atlas defined by the prior and the atlas defined by the Bayesian fit (Figure 5A and B). In this example subject, the iso-eccentricity lines in the Bayesian atlas are substantially more compressed along the posterior-anterior axis compared to the anatomical atlas, and the iso-angle lines in V2/V3 are more dorsal compared to the anatomical atlas. Where there are discrepancies, the Bayesian inferred map is more accurate. For example, the polar angles and the eccentricities in the validation data are approximately constant where the Bayesian map predicts iso-angle and iso-eccentricity lines, but not where the prior map predicts them (Figure 5B). This indicates that even after anatomical alignment, the retinotopy in this subject differs systematically from the prior. For some other subjects, the two atlases are in closer agreement such that the prior alone is a good fit to the retinotopic data (Figure 5C).

Figure 5

Download asset Open asset

To quantify the two types of deformations, we calculated the mean 3 × 3 distance matrix between (1) a vertex’s native position, (2) its position in the anatomical alignment (fsaverage position), and (3) its ‘retinotopic position’ after alignment to the retinotopic prior. All vertices in the V1-V3 region within 12° of eccentricity, as predicted by the Bayesian inference on the full dataset (all retinotopy scans combined), were used. We then performed 2D metric embedding to determine the mean deformation steps and the mean angle between them (Figure 6A). Overall, the mean deformation distance across vertices is 3.3° ± 0.6° (µ ± σ across subjects) of the cortical sphere for the anatomical alignment and 3.4° ± 0.4° for the retinotopic alignment. The mean angle between these deformations is 83.0° ± 9.6° (note that this last measurement is in terms of degrees of rotational angle rather than degrees of the cortical sphere). The anatomical alignment corresponds to variation in the surface topology and is accounted for in anatomical atlases (Benson et al., 2012; Benson et al., 2014); the retinotopic alignment corresponds to variation in the structure-function mapping and is accounted for in the Bayesian model. An additional summary measurement of these deformations, the root-mean-square deviation (RMSD) distances across vertices near the occipital pole in a particular retinotopy dataset, provides a metric of the total warping applied during each step of the alignment process for each subject. A summary of this measurement, as well as various other summary statistics is provided in Table 1.

Figure 6

Download asset Open asset

The retinotopic deformation distances were not significantly different than the the anatomical deformation distances, which were still substantial; this is true whether one calculates the mean deformation distance over the entire patch of cortex immediately around the occipital pole (shown in the maps in Figure 3) or over only the vertices that are predicted to be in V1-V3; Figure 6 is calculated over the latter of these two ROIs. Note that if the retinotopic deformation distances had been much larger, the prior anatomical atlases would have been less accurate. Had they been close to 0, then the anatomical atlas alone would have been as accurate as the Bayesian model.

We interpret the warping performed to align retinotopic data to the anatomical prior (the retinotopic alignment) as evidence of individual differences in the way in which retinotopy maps to sulcal topology. An alternate possibility is that the retinotopic alignment corrects for incomplete or incorrect warping performed by FreeSurfer during alignment of each subject’s sulcal topology to that of the fsaverage hemispheres (the anatomical alignment). We have found FreeSurfer to be a well-vetted and reliable tool for functional alignment; however, no anatomical alignment process is optimal, and more improvements to alignment algorithms, such as those being pursued using the HCP database (Glasser et al., 2016), may reduce the length of our retinotopic alignment step. (Note, additionally, that though we use the fsaverage here, our tools are compatible with other possible alignments.) It is thus possible that we have overestimated the amount of functional variance remaining across subjects after anatomical alignment. Two observations suggest that the functional variance due to imperfect anatomical alignment is small, however. First, the angle between the alignments is roughly orthogonal, meaning that there is very little movement along the axis of the first (structural) alignment during the second (retinotopic) alignment. Had the first alignment been in the correct direction but too conservative, then we would have expected the retinotopic alignment to be in the same (or similar) direction, rather than orthogonal. Second, if the retinotopic alignment served to correct an incomplete or incorrect anatomical warping, then anatomical metrics such as curvature would become more uniform across subjects after retinotopic alignment compared to after only anatomical alignment. Figure 6B demonstrates that this explanation is unlikely by showing the distribution across surface vertices of the standard deviation of curvature across subjects. When curvature is compared across subjects without retinotopic or anatomical alignment (‘Native’ alignment in Figure 6B) the standard deviation is quite high. When subjects are compared after anatomical alignment, the standard deviations are much lower. After further alignment to the anatomical prior of retinotopy (Retinotopic alignment), the standard deviation of curvature across subjects is between these two extrema. This suggests that the retinotopic alignment is sacrificing some amount of structural uniformity across subjects in order to accommodate the individual differences in subjects’ structure-to-function mapping, and is consistent with our interpretation that there are substantial individual differences in the mapping between retinotopy and surface topology.

The large individual differences that remain, even after structural co-registration (Figure 6C), point to the importance of using at least some individual subject functional data when inferring the maps, rather than assuming the atlas (prior) is correct. The specific nature of these deformations, and whether, for example, they fall into a few basic patterns, is an important question about the natural variation of individual brains. Our new method, combined with large datasets such as the HCP retinotopy data set (Benson et al., 2018) and new alignment tools such as MSMAll (Robinson et al., 2014), could be used to address this question.

To quantify the accuracy of our Bayesian-inferred retinotopic maps, and to compare the accuracy against other predictions, we used a cross-validation scheme such that predictions from data alone, the prior alone, or via Bayesian inference were compared against independent validation datasets (Supplementary file 1). The validation datasets were derived from 6 of the 12 scans; the predictions from data alone and from Bayesian inference were derived from training datasets, which were comprised of various combinations of 1–6 independent scans (between 3.2 and 19.2 min of data). The predictions from the prior alone did not use training data.

To compare methods of predicting retinotopic maps (Figure 1), vertices of interest were identified using the maps inferred from the validation dataset. All vertices from the inner 12° of eccentricity of these maps were compared to their counterparts in the predicted maps. Note that the inferred maps from the validation dataset were used only to identify the vertices included in the comparison; the retinotopic coordinates from the validation datasets themselves were taken as the ‘gold-standard’ measurements. In computing prediction accuracy, we weighted the vertices by the fraction of variance explained for each vertex’s pRF solution in the validation dataset. For the three types of training datasets (prior alone, data alone, Bayesian inference), we assume that each vertex makes a prediction regardless of the variance explained. To prevent errors at high eccentricity from dominating the error metric, we calculated the scaled error for a vertex to be the angular distance in the visual field between its retinotopic coordinates in the predicted map and the validation dataset, divided by the eccentricity in the validation dataset. Figure 7 shows the scaled mean squared error (MSE) for the various predicted maps in terms of the amount of time spent collecting retinotopic data for the map.

Figure 7

Download asset Open asset

The maps predicted via Bayesian inference were highly accurate irrespective of the amount of data used to inform the fits (Figure 7). Between those inferred maps informed by 3.2 min of scan time (one scan) and 19.2 min (six scans), the scaled MSE of the prediction remains in the range of ~0.4–0.5. These scaled errors are larger near the fovea because the denominator used for scaling the error metric (i.e., the eccentricity) could be very small; when the range is limited to 3 to 12 deg, the MSE is much lower,~0.20–0.26. Expressed separately in units of polar angle and eccentricity, the mean absolute polar angle error from a Bayesian map derived from a single 3.2 min scan is 25° ± 11° and the mean absolute eccentricity error is 0.76° ± 0.34° (µ ± σ across 96 datasets). For the prior alone, these errors are substantially higher: 34° ± 12° for polar angle and 1.3° ± 0.17° for eccentricity. Note that these errors are approximately 3 × higher than those reported for previous versions of the anatomical prior (11° for polar angle and 0.37° for eccentricity) (Benson et al., 2014); however, these discrepancies are due to differences in the metric used, the amount of data collected, the thresholding applied, and the use of smoothing. Some of these factors we cannot reproduce exactly (amount of data) or have deliberately abandoned (smoothing), but by using the same metric (median absolute error across all vertices) and similar thresholding (1.25°<predicted eccentricity<8.75°), we obtain errors very close to those previously reported: 5.9° of polar angle and 0.46° of eccentricity. In contrast to the inferred maps for which the accuracy is largely independent of scan time, the accuracy of the predictions from data alone was highly influenced by scan time. The scaled MSE of the maps predicted from the training datasets alone for the same range of scan times ranged from ~2.2 (3.2 minutes of training data) to ~0.3 (19.2 minutes of training data). With more than ~11 min of scan time, the predictions made from the training datasets alone have a slightly lower scaled MSE than those made from Bayesian inference; although, notably, the improvements in scaled MSE for more than 11 min of scan time are small.

The prediction accuracy from the prior alone (0 min. scan time) was generally intermediate in accuracy between the predictions from the Bayesian model and the data alone. Predictions derived from the anatomically-defined atlas alone are more accurate than predictions from 3.2 min of scanning and only slightly less accurate than predictions from training datasets derived from 6.4 min of scanning; this is in agreement with previous analyses of prediction error versus measurement error in retinotopic mapping experiments (Benson et al., 2014). The predictions using the prior alone were universally less accurate than the Bayesian predictions (Figure 7, cyan)—for all subjects and all training datasets, the combination of prior with data improved prediction accuracy compared to the anatomical prior alone. These data demonstrate that the application of our new method to a small amount of retinotopic mapping data yields a higher quality retinotopic map than can be derived from other sources alone, with the possible exception of data derived from a long retinotopic mapping session.

The fact that the increased prediction accuracy from the Bayesian maps is almost independent of the amount of scan time used for the observations suggests that much of the individual variability is captured by a low dimensional warping from the template, which can be inferred from a modest amount of data. This hypothesis is further supported by visual inspection of different datasets, such as in Figure 8A. Although the amount of noise in the maps clearly varies between the validation dataset (19.2 min. scan time, left column), training dataset 1 (3.2 min., second column), and training dataset 10 (6.4 min., third column), the signal is clear enough that a human expert would likely draw similar boundary lines for each map; our method does as well. Importantly, some warpings are not permitted by the fitting algorithm. For example, the topology of the template is a hard constraint, such that vertices cannot pass through one another. This puts an upper limit on the accuracy of the template: the best predictors of left out data might require a change in topology, which is not permitted. We discuss the significance of these issues in the section subsequent section, 'What is ground truth’.

Figure 8

Download asset Open asset

Another significant advantage of the method is that it eliminates the need for human intervention in the process of delineating retinotopic maps and visual areas. In most studies that require retinotopic mapping data, one or more experimenters hand-label the visual area boundaries. While human raters are better able to understand atypical retinotopic boundaries than our method, they are nonetheless subject to inter-rater disagreement and human error. Furthermore, although expert human raters have a much more nuanced prior about retinotopic map organization than our method, and thus may sometimes draw boundaries better than our method, our method at least makes its prior explicit and quantifiable, and, thus, comparable and replicable across studies.

The motivation for a Bayesian approach to retinotopic mapping can be found most clearly in the measured retinotopic maps themselves. In all of our measured retinotopic maps, there are numerous systematic imperfections (Figure 8A), and the literature contains many reports of similar errors (Winawer et al., 2010; Press et al., 2001; Gardner, 2010; Boubela et al., 2015). These imperfections can arise from a variety of sources, including partial voluming (Huettel et al., 2014), negative BOLD (Shmuel et al., 2002), and large, draining veins. Imperfections due to blood vessel artifacts can have effects over large distances (Winawer et al., 2010), and most perniciously, they may lead to large and reliable responses that nonetheless differ from the local neuronal activity in the voxel (Boubela et al., 2015). Such artifacts can be difficult to track down and are often not eliminated by typical methods of cleaning up maps such as smoothing, thresholding, or simply collecting larger datasets.

Although with large datasets (>19 min of scan time), the prediction accuracy for the validation dataset is highest using the data alone rather than the data and prior, we believe that even in these cases the combination of data + prior is probably closest to ground truth. We defined accuracy operationally as the difference from the validation set, as this provides a single set of independent measures that can be used to assess the accuracy of all three types of models. The validation dataset is defined by at least as many scans as any of the training datasets, and hence is our best measurement. However, the validation dataset is not ground truth, as it is subject to errors from systematic and random measurement noise. The inferred maps, unlike the maps from the validation and training datasets alone, produce a topologically smooth transformation into the visual field: they represent the complete visual field with no holes, redundancies, or discontinuities. Hence, an enclosed region in the visual field will project to an enclosed region in the inferred map on the cortical surface. This is not the case for the maps predicted by the data alone without the use of an atlas (Figure 8B). We consider this difference to be an advantage of the inferred maps, since it is generally accepted that the cortical surface of V1 is a topological map of the visual field.

In short, since we do not correct for all of these potential sources of systematic error, we consider our estimates of error from the Bayesian-inferred maps to be conservative, and the estimate from the data-to-data predictions to be liberal.

One important advantage of using the method of Bayesian inference outlined in this paper is that it provides predictions beyond the extent of the stimulus aperture in the retinotopy experiment. These peripheral predictions extend to 90° of eccentricity, even though the data used to derive the prior was based on stimuli that only extended to ~8° of eccentricity. Hence, it is important to ask whether the model makes accurate predictions in the periphery. We demonstrate this in two ways. First, when our registration algorithm is run using only a subset of the eccentricity range (e.g., only data within the first 3° or the first 6° of eccentricity), the predicted maps remain accurate to 12° of eccentricity (Figure 9A). Second, we compared wide-field retinotopy data, collected out to 48° of eccentricity from subject S1201 to the Bayesian-inferred map predictions made using our data with a 12°-aperture (Figure 9B). We find that in both cases, our method is highly accurate despite lacking training data for peripheral measurements (though note that in the latter case, the extrapolation was only tested on one subject; in principle, subjects with poorer data quality or unstable fixation could result in less accurate extrapolation). Because the extrapolations to untested eccentricities are generally accurate, we conclude that even if prediction accuracy within the measured regions were similar for the Bayesian model and the data-to-data predictions, the Bayesian model is advantageous because it includes predictions for regions of the visual field beyond training data.

Figure 9

Download asset Open asset

Another aspect in which our work here extends previous methods is the addition of the pRF size to the retinotopic quantities predicted by the model in the inferred maps—previous work predicted only the pRF centers (Benson et al., 2012; Benson et al., 2014). Here, we predict the pRF sizes for the vertices based on the eccentricity inferred from the Bayesian map and the assumed linear relationship between eccentricity and pRF size. The inferred pRF size of a vertex is the best linear fit to the measured pRF size versus the vertex’s inferred eccentricity (Figure 10A). While an approximately linear relationship is reasonable given the literature, the absolute scale is likely dependent on variety of measurement factors such as voxel size, stimulus spatial frequency, and subject fixation (Alvarez et al., 2015). Hence, we do not attempt to infer the slope or intercept based on prior measurements.

Figure 10

Download asset Open asset

Another metric inversely related to pRF size is the cortical magnification, usually measured in terms of mm² of cortex representing one degree² of the visual field. We summarize these measurements in Figure 10B and C. Our measurements of cortical magnification are broadly in agreement with previous work by Horton and Hoyt, 1991, shown by the dotted black line panels B-E of Figure 10. The cortical magnification of the inferred maps is quite similar to that of the observed retinotopic maps. In both cases, V1 has slightly lower cortical magnification than V2 and V3 near the fovea, but higher magnification in the periphery. This difference is slightly exaggerated in the inferred maps relative to the observed maps; although this difference is slight and is in agreement with previous examinations of cortical magnification (Schira et al., 2009); however, note that in our maps, the crossover between V1 cortical magnification and V2/V3 cortical magnification occurs at a higher eccentricity (~3°) than previously reported (~0.7–1°). This is emphasized in Figure 10D and E, which shows the curves from Figure 10C and D in terms of their difference from the prediction of Horton and Hoyt, 1991 (the black dashed line in panels B-E). Note that in the inferred maps, although the cortical magnification in V1 is lower than V2 below 3° of eccentricity, the difference between them is small between 1.5° and 3°.

Previous research on the retinotopic organization of visual cortex used a model of V1, V2, and V3 retinotopy described by Schira et al. (2010) to produce a template of retinotopy that included only those visual areas. This ‘banded double-sech’ model accurately describes the anisotropic magnification of the visual field on the cortical surface, particularly near the fovea. However, we have observed, particularly in individual data, that retinotopic data from outside the V1-V3 region described by the Schira model has a large impact on the quality of the inferred map. Accordingly, in creating our retinotopic prior, we constructed a new model that includes nine additional visually active regions: hV4, VO1, VO2, V3a, V3b, LO1, LO2, TO1, and TO2. This model employed a new method of constructing 2D models of retinotopy that was specifically designed to accommodate the distortions caused by anatomical alignment, inflation, and flattening of the cortical surface. The new method is much simpler to extend to many more visual field maps, as it does not rely on an analytic description of the flattened (2D) retinotopic maps, which are only available for V1-V3. Rather, it requires as input a cortical map image on which estimates of the visual area boundaries have been drawn manually and labeled as either a foveal boundary, a peripheral boundary, an upper vertical meridian, a lower vertical meridian, or a horizontal meridian. A minimization technique is then used to fill in the retinotopic coordinates between the drawn boundaries (see Models of Retinotopy in Materials and methods). The new retinotopic prior, including all new areas can be seen in Supplementary file 3. Although we consider the addition of 9 retinotopic areas to be an important development, we consider these areas preliminary and do not analyze them in detail here. One reason for this is that the organizations of many of these areas remain under dispute. Additionally, the responses to our stimuli in these areas is of a considerably lower quality than in V1-V3; thus even were we to analyze the accuracy of the predictions in these ares, our validation dataset would be a particularly poor standard. We do, however, include these areas in the predictions made by the Bayesian inference method so that interested researchers may analyze or use them. These areas are included in the data provided with this paper, and predictions of these areas are included when using the tools provided in the Data Availability Statement.

Our method has the advantage of allowing the retinotopic atlas to act as a prior constraint on new observed data. This is a Bayesian model in the general sense of combining a prior belief with a measurement in order to make an inference. The computation can be formulated in an explicit Bayesian framework. We define a hypothesis H to be a particular warping of the cortical surface, and we define the evidence E to be a particular set of retinotopic measurements. We then convert the cost functions from Table 2 into probabilities by assuming an exponential relationship. Hence, the prior probability of H is defined in terms of the deviation from the retinotopic prior: ${\displaystyle P(H)=\mathrm{e}\mathrm{x}\mathrm{p}(-(F_e(x)+F_{\theta }(x)+F_p(x)))}$, and the likelihood of the evidence under a given hypothesis, ${\displaystyle P(E|H)}$, is defined in terms of the fit between the retinotopic model and the retinotopic measurements: ${\displaystyle P(E|H)=\mathrm{e}\mathrm{x}\mathrm{p}(-(E_e(x)+F_{\theta }(x)+F_p(x)))}$. During registration, we seek the hypothesis ${\displaystyle H}$ that maximizes the posterior probability ${\displaystyle P(H|E)=P(E|H)P(H)/P(E)}$. Because ${\displaystyle P(E)}$ is a constant, we can ignore it and instead maximize the function given in Equation 1, which is equivalent to minimizing ${\displaystyle F(x)}$. This operation is performed during registration. Thus, to derive our cost function from Bayes’ rule, we write:

The explicit Bayesian formulation above clarifies several features of our model. First, the prior probability distributions assumed for vertex lengths are the same for all vertices (rows 1, 5, and six in Table 2; P(H) in Equation 1). If we had ground truth maps for a large population, we could, in principle, derive edge-specific probability distributions for Equation 1, and convert these to edge-edge-specific cost functions (Table 1) for the minimization process. We can get a sense of how these distributions might differ across occipital cortex by visualizing the warp fields from our data set (Supplementary file 1). These fields show that our registration process causes some vertices to move much more than others, at least in our small subject pool (n = 8). These warp fields are not sufficient to derive edge-specific priors because the number of subjects is small, and because we do not know the end-points of the registration reflect the ground truth maps. The use of a large dataset, such as the 181 HCP subjects (Benson et al., 2018) might be helpful in future work to derive edge-specific priors. A further challenge to incorporating realistic priors would be to capture the dependencies across edges in the prior distribution (the joint probability distribution, which would be a function of thousands of variables, one per edge, imposing an enormous computational burden).

A second feature of our method made explicit by the Bayesian formulation is that the prior probability distributions are 0 for solutions that violate the atlas topology. This assumption is implemented implicitly in the cost function, which rises to infinity as the length of an edge approaches 0 or the angle between edges approaches 0. This aspect of the cost function prevents vertices or edges from crossing, thus preserving topology. Because we assume that the cost function is the negative logarithm of the prior probability distribution, the infinite cost indicates an assumed probability of 0. If ground truth data contradicted this assumption (that is, if there were ground truth maps which violated the topology of the model), the prior probability distributions and corresponding cost functions could be changed accordingly.

A third feature of the method is that the likelihood functions depend on the data quality. In Table 1, line 4, the weighting of each vertex (w) is proportional to the variance explained by the pRF model. PRF solutions with high variance explained lead to a higher cost when the atlas vertex is far from the corresponding data point. This part of the cost function shows up in the likelihood, F_φ(x), in the Bayesian formulation (Equation 1). The interpretation is that there is a low likelihood of observing a high-variance-explained pRF solution in a location far from the template solution. A more realistic likelihood calculation (but one that is beyond the scope of our current knowledge and computational resources) would require a noise model that allowed one to compute how likely a pattern of pRF solutions was given a hypothesized map.

An important question in human neuroscience is the degree to which different brains, when brought into anatomical registration, share the same functional mapping. There is no single, agreed-upon method to register the brains of different individuals, but a general finding is that cortical function shows better inter-observer agreement when the brains are aligned based on sulcal topology (surface registration) rather than volume registration (Wang et al., 2015; Van Essen et al., 1998). Here, using surface registration, we find that substantial individual differences in functional mapping remains, for example as evidenced by the amount of additional warping needed to align individual brains to retinotopic measurements (Figure 6). These results are consistent with studies showing differences in structural and functional alignment of primate and human area MT (Large et al., 2016). They are also consistent with studies combining structural and functional alignment in the absence of a model or template of the underlying function (Frost and Goebel, 2013; Haxby et al., 2011). Such studies show that the extra warping driven by functional alignment leads to better predictions of functional responses in cross-validated data.

The anatomical atlas, described previously (Benson et al., 2014), was adapted into the first step of our method (Figure 4C) and is equivalent to the prediction using the retinotopic prior alone that we present here (Supplementary file 2B). Consistent with previous results, we find that the prior alone produces reasonably good predictions for the retinotopic maps of most subjects (Figure 7). Additionally, maps predicted using the prior alone contain many of the advantages described here such as topological smoothness, complete coverage of the visual field, and prediction of peripheral retinotopy. However, the idiosyncrasies of individual subjects’ retinotopic maps are often not well predicted by the retinotopic prior alone (Figures 3, 5A and B). By combining both the retinotopic prior and a small amount of measured data, we are able to produce higher-quality predictions that not only share these advantages but also improve the prediction accuracy of the maps beyond that of the measurement error (Figure 7).

We have shown in this paper that the application of Bayesian inference to retinotopic mapping data can yield substantial rewards for basic researchers interested in quantitative retinotopic analysis. However, a number of questions about the scope and limitations of the method remain. For one, it is unclear how our method, developed using data from a small subset of the population (eight subjects), would cope with subjects whose retinotopic organizations are much different than those that are typically assumed in vision science. Such edge-case subjects could include members of clinical populations, such as individuals lacking an optic chiasm (Hoffmann et al., 2012; Bao et al., 2015; Olman et al., 2018), or healthy individuals whose retinotopic boundaries are merely unusual (Van Essen and Glasser, 2018). We consider how our Bayesian models perform for edge-cases by fitting the models to two subjects from the Human Connectome Project whose retinotopic maps were recently noted for their peculiarity by Van Essen and Glasser (2018). The retinotopic maps for these subjects as well as the inferred iso-angular and iso-eccentricity lines in V1-V3 are shown in Supplementary file 5. Both of these subjects have atypical polar angle organization in their left hemisphere dorsal V2 and V3 maps. Regarding the limitations of the method as applied to subjects such as these, two things are clear: first, the method is unable to reproduce the precise topology of the subjects’ unusual dorsal maps, and, second, it is nonetheless capturing most of the maps accurately. In particular, the inferred eccentricity maps are highly accurate despite the mismatched polar angle maps. The polar angle maps cannot be accurately captured by the Bayesian model, because the Bayesian model assumes a prior probability of zero to any solutions that differ topologically from the template, such as these.

In the case of more extreme departures from typical retinotopic organization, such as the overlapping maps observed in achiasmic patients, we cannot be certain that our method would yield coherent inferences. However, we note that this problem is not unique to our method; in fact, generic voxel-wise pRF models also fail for subjects with highly unusual retinotopic organizations such as achiasma. Typically, these model failures are identified, and updates to the model are proposed to account for the relevant conditions. In the case of achiasmic patients, the usual assumption that a pRF can be described as a single Gaussian fails, and new models were developed with two spatially displaced Gaussians (Hoffmann et al., 2012). In the case of our method, one might use an alternate formulation of our retinotopic prior in order to better model the maps of the subjects shown in Supplementary file 5 or subjects from a particular clinical population.

A separate but equally critical question about the method we present is whether it encapsulates any systematic biases about retinotopic organization. Given the field’s imperfect knowledge about precise retinotopic organization across individuals, this can be a difficult question to answer; however we note a number of features and assumptions along these lines. With respect to the retinotopic prior, one critical assumptions that was employed during its creation involves the structure of the polar angle reversals at map boundaries (e.g. the V1/V2 boundary or the V3/hV4 boundary). Our prior assumes that the polar angle at these boundaries lies on the vertical meridians (or horizontal meridians in the case of the V2/V3 boundaries). In the group-average retinotopic maps from the 181 HCP subjects (Supplementary file 3), however, it is clear that many polar angle reversals occur several degrees away from the vertical (or horizontal) meridian, differing from the prior. In these cases, our Bayesian maps differ systematically from the data. Because the Bayesian computation allows the vertex positions to change but does not allow the retinotopic quantities to change, all retinotopic locations contained in the prior are assigned to some cortical location in the Bayesian map, differing from the observed validation data, which is often missing representations of the vertical meridians. In addition, ‘notches’ of missing representation of the lower vertical meridian can be seen in several polar angle maps (for example, inSupplementary file 3A, black arrows). In both cases—maps that do not quite reach the vertical meridian and maps that have large notches along the boundaries—it is not yet known whether these properties reflect unusual features of the underlying neuronal maps or limits of the fMRI acquisition or analysis.

The method described in this paper has been made publicly available in a number of ways in order to enable easy use by other researchers. The Bayesian inference method itself is implemented as part of a free software library called Neuropythy; we have publicly archived the version of this library used in the preparation of this manuscript at DOI:10.5281/zenodo.1312983. Additionally, we have created a universally executable Docker image and have publicly archived it at DOI:10.5281/zenodo.1313859. Detailed instructions for the use of both the library and the universally executable Docker image are available at the Open Science Foundation repository associated with this paper (https://osf.io/knb5g/). In brief, the method may be run with only a few inputs: a FreeSurfer directory for the subject (which can be generated from an anatomical image), and a set of files containing the measured retinotopic parameters for the subject’s cortical surface. The outputs produced are a similar set of files describing the inferred retinotopic parameters of the subject’s cortical surface. Runtime on a contemporary desktop computer is less than an hour per subject. Detailed instructions on how to use the tools documented in this paper are included in the Open Science Foundation website mentioned above.

All source code and notebooks as well as all anonymized data employed in this Methods section and the preparation of this manuscript have been made publicly available at the Open Science Foundation: https://osf.io/knb5g/. Version 0.6.0 and later of the Neuropythy library can automatic download these data and interpret them into Python data structures.

This study was approved by the New York University Institutional Review Board, and all subjects provided written consent. A total of eight subjects (4 female, mean age 31, range 26–46) participated in the experiment. All scan protocols are described below.

All MRI data were collected at the New York University Center for Brain Imaging using a 3T Siemens Prisma scanner. Data were acquired with a 64-channel phased array receive coil. High resolution whole-brain anatomical T1-weighted images (1 mm³ isotropic voxels) were acquired from each subject for registration and segmentation using a 3D rapid gradient echo sequence (MPRAGE). BOLD fMRI data were collected using a T2*-sensitive echo planar imaging pulse sequence (1 s TR; 30 ms echo time; 75° flip angle; 2.0 × 2.0 × 2.0 mm³ isotropic voxels, multiband acceleration 6). Two additional scans were collected with reversed phase-encoded blips, resulting in spatial distortions in opposite directions. These scans were used to estimate and correct for spatial distortions in the EPI runs using a method similar to (Andersson et al., 2003) as implemented in FSL (Smith et al., 2004).

Anatomical images were processed using the FreeSurfer image analysis suite, which is freely available online (http://surfer.nmr.mgh.harvard.edu/) (Dale et al., 1999; Fischl et al., 1999a; Fischl et al., 1999b; Fischl and Dale, 2000). Subject brains were inflated and aligned to FreeSurfer’s anatomical fsaverage atlas.

Each subject participated in 12 retinotopic mapping scans using the same stimulus employed in the Human Connectome Project (Benson et al., 2018). Briefly, bar apertures on a uniform gray background swept gradually across the visual field at four evenly-spaced orientations while the subject maintained fixation. Bar apertures contained a grayscale pink noise background with randomly placed objects, faces, words, and scenes. All stimuli were presented within a circular aperture extending to 12.4° of eccentricity. The bars were a constant width (1.5°) at all eccentricities. Subjects performed a task in which they were required to attend to the fixation dot and indicate when its color changed.

The 12 scans were split into several subsets and analyzed as independent datasets. Six of the scans (two of each bar width) were allocated to the subject’s validation dataset, while the remaining three scans were used for 21 training datasets: 6 datasets with one scan each, 5 datasets with two scans each, 4 datasets with three scans each, 3 datasets with four scans each, 2 datasets with five scans each, and one dataset with all six training scans. Additionally, all 12 scans were included in a full dataset which was used for all analyses not related to the accuracy of the Bayesian inference method.

Additionally, one previously published retinotopy dataset for with a wide field of view (48° of eccentricity) was re-analyzed (Wandell and Winawer, 2011) (their Figure 3). The subject for this dataset was also included in the newly acquired data, S1201. The wide-field-of-view dataset was used as a further validation set for models derived from the newly acquired data for S1201, as it enabled us to test the accuracy of model predictions in the far periphery from models derived from data with limited eccentricity.

Spatial distortions due to inhomogeneities in the magnetic field were corrected using in-house software from NYU’s Center for Brain Imaging (http://cbi.nyu.edu/software). The data were then motion-corrected by co-registering all volumes of all scans to the first volume of the first scan in the session. The fMRI slices were co-registered to the whole brain T1-weighted anatomy, and the time series resampled via trilinear interpolation to the 1 mm³ voxels within the cortical ribbon (gray matter). Finally, the time series were averaged for each voxel across all scans with the same stimulus within a given dataset.

Retinotopic maps were produced by solving linear, circularly symmetric population receptive field (pRF) models for each voxel using Vistasoft, as described previously (Dumoulin and Wandell, 2008). pRF models were solved using a two-stage coarse-to-fine approach on the time series in the 1 mm³ gray matter voxels. The first stage of the model fit was a grid fit, solved on time series that were temporally decimated (2×), spatially blurred on the cortical surface using a discrete heat kernel (approximately equal to a Gaussian kernel of 5 mm width at half height), and subsampled by a factor of 2. The decimation and blurring helps to find an approximate solution that is robust to local minima. The parameters obtained from the grid fit were interpolated to all gray matter voxels and used as seeds for the subsequent nonlinear optimization. Finally, the pRF parameters were projected from the volume to the cortical surface vertices for white, mid-gray, and pial surfaces using nearest-neighbor interpolation; values were then averaged across the three layers using a weighted-mean in which the fraction of BOLD signal variance explained by the pRF model was used as a weight. All vertices with a pRF variance explained fraction less than 0.1 were ignored.

To generate our initial models of the retinotopic maps, we begin by hand-drawing boundaries for 12 retinotopic maps. These boundaries need only be drawn once for a single group-average retinotopic map. The boundaries are projected onto the cortical surface, and the retinotopic coordinates for each vertex on the surface are deduced via a minimization procedure. This minimization is motivated by two principles: (1) the retinotopic fields (polar angle and eccentricity) should be as orthogonal to each other as possible and (2) the retinotopic fields should be as smooth as possible. To this end, the minimization routine simultaneously maximizes both the smoothness of the retinotopic fields between vertices connected by edges in the mesh as well as the overall orthogonality between the polar angle field and the eccentricity field. The hand-drawn boundary values are held constant during the minimization. Equation 2 gives the function ${\displaystyle f}$ that is minimized, where ${\displaystyle \theta }$ and ${\displaystyle \varrho }$ represent vectors of the polar angle and eccentricity values, respectively; ${\displaystyle \mathrm{E}}$ represents the set of edges between vertices in the mesh; and ${\displaystyle \mathrm{x}}$ represents the matrix of vertex coordinates (i.e., ${\displaystyle {\mathrm{x}}_u}$ represents the coordinate vector of the vertex u). During minimization the values of the polar angle and eccentricity vectors are scaled such that both fields ranged from −1 to 1 (e.g., polar angle boundary values of 0° and 180° were assigned values of −1 and 1, respectively) so that the fields could be evaluated easily for orthogonality; after the minimization, the polar angle was linearly rescaled back to the range 0°−180° while eccentricity was rescaled so as to have an exponential distribution that best fit the group-average. We employed this model generation routine using boundaries drawn over the group-average data (see Anatomically-defined Atlas of Retinotopy: Group Average, below) as well as the Wang et al. (2015) atlas and the Hinds et al. (2008) V1 boundary as rough guidelines.

A full description of the model, including how it can be projected onto an fsaverage spherical surface or an individual subject’s fsaverage-aligned spherical surface, is provided in the open source Neuropythy library (https://github.com/noahbenson/neuropythy).

Construction of the anatomically-defined atlas of retinotopy is summarized in Supplementary file 2. Previous work employed a mass-spring-damper system combined with a nonlinear gradient-descent minimization in order to register group-average retinotopic data, averaged on FreeSurfer’s fsaverage_sym hemisphere (Greve et al., 2013), with a model of V1, V2, and V3 retinotopy (Schira et al., 2010). In this paper, we modify this technique slightly to bring it more in line with previous established methods such as those used by FreeSurfer for surface-based anatomical registration (Dale et al., 1999; Fischl et al., 1999a). In brief, retinotopy is measured in a group of subjects via fMRI; the subjects’ cortical meshes are aligned to the fsaverage surface via FreeSurfer’s surface registration; the retinotopic coordinates are then averaged across subjects at each vertex on a single atlas of the cortical surface; a 2D atlas of retinotopy is then placed on this cortical surface; and finally, the cortical surface is warped to match the retinotopic atlas as best as possible given constraints on the warping. Each of these steps is described in more detail below.

Group-average Data. Group-average retinotopic maps (Supplementary file 2B) were obtained from 181 subjects whose data were published and made freely available as part of the Human Connectome Project (Benson et al., 2018). The resulting group-average retinotopic maps are shown in Supplementary file 2B.

Cortical Map Projection. The cortical surfaces of the fsaverage left and right hemispheres, on which the group-average data were constructed, were inflated both to a smooth hemisphere (FreeSurfer’s ‘inflated’ surface) as well as to a sphere (FreeSurfer’s ‘sphere’ surface); the vertices on the spherical surfaces were then flattened to 2D maps using an orthographic map projection. Precise parameters of this projection and the source code used to generate it are included in the Data Availability Statement. We refer to the 2D vertex coordinates in this resulting map as the ‘initial vertex coordinates’ because they precede the warping of the vertex coordinates that occurs during registration.

Registration. The initial vertex coordinates of the map projections described above were warped in order to bring the polar angle and eccentricity measurements of the vertices into alignment with the 2D model’s predictions of retinotopy while maintaining topological constraints: that is preventing triangles in the triangle mesh representing the 2D cortical map from inverting and penalizing excessive stretching or compression of the map. This process was achieved by minimizing a potential function defined in terms of the edges of the triangle mesh, the angles of the triangle mesh, and the positions of the vertices with polar angle and eccentricity measurements above the weight threshold (see Group-Average Data, above). Equation 3 gives this potential function, ${\displaystyle F(x)}$, which is further broken down into four components detailed in Table 2. Fundamentally, the potential function F is a sum of two kinds of penalties: penalties for deviations from the reference mesh and penalties for mismatches between the vertices with retinotopic coordinates and their positions in the retinotopic model. In the case of the former, the reference mesh is gi ${\displaystyle x_0,E,\mathbf{\Theta },P,\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{\Phi }}$ and the potential of the deviations are defined by ${\displaystyle f_e}$, ${\displaystyle f_{\theta }}$, and ${\displaystyle f_p}$. The latter is described by ${\displaystyle f_{\phi }}$. In these functions, ${\displaystyle \mathbf{x}}$ represents the n × 2 matrix of the 2D-coordinates of each vertex while ${\displaystyle {\times }_0}$ represents the same coordinates in the reference mesh; E represents the set of undirected edges (represented as ${\displaystyle (u,v)}$ pairs such that ${\displaystyle (u,v)}$ and ${\displaystyle (v,u)}$ are not both in E) in the reference mesh; ${\displaystyle \mathbf{\Theta }}$ represents the set of angle triples (a, b, c) such that the angle is between edge ${\displaystyle (a,b)}$ and edge ${\displaystyle (a,c)}$; P is the set of vertices that lie on the perimeter of the 2D map projection; r_x(u,v) is the Euclidean distance between vectors ${\displaystyle (\mathbf{x}{)}_u\mathrm{a}\mathrm{n}\mathrm{d}(\mathbf{x}{)}_v}$; and ${\displaystyle {\alpha }_x(a,b,c)}$ is the counter-clockwise angle between vectors ${\displaystyle ((x)b-(x)a)\mathrm{a}\mathrm{n}\mathrm{d}((x)c-(x)a);\mathbf{\Phi }}$ represents the set of anchors defined by the retinotopic model in which each anchor is represented by a tuple ${\displaystyle (u,\mathbf{y},\sigma ,w)}$ where w is the weight of the anchor, ${\displaystyle u}$ is the vertex drawn to the anchor, ${\displaystyle \sigma }$ is the standard deviation of the anchor’s Gaussian potential well, and ${\displaystyle y}$y is the 2D point to which the anchor is attached; the constants ${\displaystyle q_0\mathrm{a}\mathrm{n}\mathrm{d}q1}$ are the minimum and maximum allowable edge lengths, respectively.

The term of the potential function devoted to the retinotopic model is given in F_φ (Equation 3; Table 2). This potential term is a set of inverted-Gaussian potential wells called anchors. Each anchor represents the attraction of a single vertex u, with measured polar angle ${\displaystyle \theta }$, eccentricity ${\displaystyle \varrho }$, and weight w, to a 2D point y, at which the retinotopic model predicts a polar angle value of ${\displaystyle \theta }$ and an eccentricity value of ${\displaystyle \varrho }$. Note that each visual area represents every point ${\displaystyle (\theta ,\varrho )}$ in the visual field, there are multiple anchors per vertex with retinotopic data. In fact, the retinotopic model used in this paper defines nine maps in addition to the V1-V3 maps (see Model of Retinotopy, above), bringing the total number of anchors per retinotopic vertex to 12. The additional areas are intended partly to prevent vertices immediately outside of V1-V3 from being drawn incorrectly into the V1-V3 section of the model and are not analyzed in detail in this paper. Each anchor additionally defines a parameter σ; this value is the width (standard deviation) of the anchor’s Gaussian potential well; σ is defined as the minimum distance from the given anchor to any other anchor to which u is also attracted; this value was given a maximum value of 20ε where ε is the mean edge-length in the projected map.

The potential function was minimized using a gradient descent algorithm sensitive to the singularities in the terms G_e, and G_θ (Table 2); whenever the singularity is accidentally crossed, the minimizer backtracks and chooses a smaller step-size. This approach prevents the inversion (from counter-clockwise ordering to clockwise ordering) of any triangle in the mesh, as such an inversion would require the minimization trajectory to pass through a singularity at the point where α = 0 or α = π. The source code used to minimize the potential function as well as specifications of the gradients of each term is provided in the open-source library included with the Neuropythy and Neurotica libraries (https://github.com/noahbenson/nben).

Minimization was run for at least 2500 steps in which the step-size was constrained such that the displacement of each vertex in each step was at most 1/50^th of the average edge-length in the map projection. A small amount of exponentially distributed random noise was added to the gradient at each step with the constraint that the gradient direction at each vertex be conserved; this noise did not affect the minimum obtained by the search but did speedup convergence significantly (see associated libraries for further details). Convergence was generally observed within 1000–2000 steps. The set of vertex coordinates that resulted from this minimization brings the retinotopic measurements associated with the vertices in V1, V2, and V3 referred to as the registered vertex coordinates.

Prediction. The registered vertex coordinates, once obtained, give the alignment of the subject’s cortical surface to the model of retinotopy; accordingly, a prediction of any vertex’s associated pRF and visual area label can be derived by comparing the the vertex’s registered coordinates with the model. Every vertex whose registered coordinates fall within the model’s V1 boundary, for example, is labeled as part of V1. Because only the vertex coordinates, and not the vertex identities, are changed during the registration process, there is no need to invert the registration: visual area label, polar angle, and eccentricity values assigned to each vertex apply as readily to the vertices whether they are visualized in the registered vertex coordinates or in the coordinates that define the subject’s white-matter surface, for example. The retinotopic map predictions for the group-average data is shown in Supplementary file 2D (left column).

Because the group-average retinotopic data were used in the registration, the predicted map that results provides a reasonable estimate of any subject’s expected retinotopic map, as shown previously (Benson et al., 2012; Benson et al., 2014); although the predicted map does not account for further individual differences in the structure to function relationship, as we show in this paper. Additionally, because the predicted map from the group-averaged data is defined on the fsaverage subject’s cortical surfaces, a retinotopic map prediction for any new subject, for whom retinotopic mapping measurements may not be available, can be easily obtained: one can use FreeSurfer to align the new subject’s cortical surface with the fsaverage subject’s surface (anatomical structure alignment) then to project the retinotopic maps from the fsaverage subject to the new subject based on the anatomical similarity between them. Because of this, we refer to this group-average retinotopic prediction as the anatomically-defined atlas of retinotopy. This atlas is used as the prior for the Bayesian model fit, described below. The atlas is similar but not identical to one presented previously (Benson et al., 2014).

The anatomically-defined atlas of retinotopy, while providing a good prediction for most subjects’ individual retinotopic maps, nonetheless does not account for individual differences in the mapping between anatomical location and retinotopic coordinates. Accordingly, predicted retinotopic maps for individual subjects were refined starting from the anatomically-defined atlas of retinotopy using a similar method as was used to generate the atlas originally; this process is detailed in Figure 4C. For each subject, their cortical surface was aligned based on anatomical structure to the fsaverage subject’s cortical surface using FreeSurfer, then their retinotopic data were projected to a map using the identical map projection described above in the section on the anatomically-defined atlas of retinotopy. Note that, in this case, the anatomical alignment to the fsaverage subject serves to make the map projections as similar as possible between individual subjects and the anatomically-defined atlas of retinotopy. If we were not interested in incorporating information obtained from the anatomically-defined atlas of retinotopy (which represents a prior belief of retinotopic organization based on group-average data), this step would not be necessary.

The individual subject’s projected map is then arranged according to the registered vertex coordinates from the anatomically-defined atlas of retinotopy; this step reflects the prior belief that the group-average registration to the retinotopy model is generally accurate for an individual subject when that subject’s anatomical structure has been aligned to the fsaverage subject’s. Critically, none of the steps taken so far in processing the individual subject’s data relies on any measurements of retinotopy that might be associated with that subject. Rather, these steps have relied only on anatomical structure. If, for a subject, no retinotopic measurements are made, then there is no data with which to modify this prior belief; accordingly, the prediction of retinotopy for that subject would be identical to the prediction of retinotopy contained in the anatomically-defined atlas. In other words, without observation, the prior remains the prediction.

The next step registers the individually measured retinotopy data to the anatomically-defined atlas. Before registration, the individual subject’s data is resampled onto a uniform triangular mesh, and each vertex whose retinotopic measurements are above threshold are given a weight, w, based on the variance explained, ${\displaystyle \omega }$, of its pRF model solution. The mesh is resampled to the same uniform triangle mesh used as the initial vertex coordinates in the registration of the anatomically-defined atlas of retinotopy in order to speedup registration. Triangles that are tightly pinched (i.e., triangles with internal angles near 0 or π) can drastically slow the registration progress by forcing the minimizer to frequently backtrack steps; resampling makes such behavior much less likely during the initial minimization. Aside from the weight, other parameters tracked by the potential field, including anchors parameters used by the function F_φ(x), are obtained identically as with the anatomically-defined atlas of retinotopy. These anchors inherit the weight of the vertex to which they apply, but are reduced when the field sign of the triangles adjacent to the vertex does not match the field sign of the visual area to which it is tied by the anchor or when the pRF size predicted by the model does not match that of the vertex’s measured pRF. Details regarding the weights on anchors are provided in the neuropythy library.

For each training dataset of each subject, minimization was run for 2500 steps using the same protocol that was used with the anatomically-defined atlas of retinotopy. Retinotopic map prediction, based on the positions of the registered vertex coordinates in the retinotopy model, were also computed identically to those in the anatomically-defined atlas. Identical minimization and prediction methods were run for each test dataset as well, but these results were not used to measure the accuracy or effectiveness of the prediction methods.

Cortical magnification was calculated using both the observed retinotopic maps and the inferred maps that were produced by combining each subject’s full retinotopy dataset with the retinotopic prior. This combination of data should, in theory, produce the highest-quality retinotopic map predictions of which we are capable (see Results and Discussion). Cortical magnification was calculated by first projecting all vertices in a single visual area (such as V1) into the visual field based on their pRF centers. The cortical magnification of a particular polar angle and eccentricity is then the total white vertex surface-area (as calculated by FreeSurfer) of all pRF centers within a disk of some radius α, divided by the area of the disk (πα²). For an eccentricity ρ, we used a radius α = ρ/3.

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

The research was supported by NIH grants R00-EY022116, R01-MH111417, R01- EY027964, and R01-EY027401. The authors declare no competing financial interests.

All figures and plots were produced using the Neurotica library (https://github.com/noahbenson/neurotica) for Mathematica 11 (Wolfram Research, 2018). Registrations and interpolations were performed with the assistance of the nibabel (DOI:10.5281/zenodo.1287921), numpy (Oliphant, 2006), and scipy (Jones, 2001) libraries for Python.

Human subjects: This study was conducted with the approval of the New York University Institutional Review Board (IRB-FY2016-363) and in accordance with the Declaration of Helsinki. Informed consent was obtained for all subjects.

1. Received: July 18, 2018
2. Accepted: November 29, 2018
3. Accepted Manuscript published: December 6, 2018 (version 1)
4. Version of Record published: January 21, 2019 (version 2)

© 2018, Benson and Winawer

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