Topology-preserving smoothing of retinotopic maps.

Retinotopic mapping, i.e., the mapping between visual inputs on the retina and neuronal activations in cortical visual areas, is one of the central topics in visual neuroscience. For human observers, the mapping is obtained by analyzing functional magnetic resonance imaging (fMRI) signals of cortical responses to slowly moving visual stimuli on the retina. Although it is well known from neurophysiology that the mapping is topological (i.e., the topology of neighborhood connectivity is preserved) within each visual area, retinotopic maps derived from the state-of-the-art methods are often not topological because of the low signal-to-noise ratio and spatial resolution of fMRI. The violation of topological condition is most severe in cortical regions corresponding to the neighborhood of the fovea (e.g., < 1 degree eccentricity in the Human Connectome Project (HCP) dataset), significantly impeding accurate analysis of retinotopic maps. This study aims to directly model the topological condition and generate topology-preserving and smooth retinotopic maps. Specifically, we adopted the Beltrami coefficient, a metric of quasiconformal mapping, to define the topological condition, developed a mathematical model to quantify topological smoothing as a constrained optimization problem, and elaborated an efficient numerical method to solve the problem. The method was then applied to V1, V2, and V3 simultaneously in the HCP dataset. Experiments with both simulated and real retinotopy data demonstrated that the proposed method could generate topological and smooth retinotopic maps.

This is a PLOS Computational Biology Methods paper.

Introduction

Retinotopic mapping, i.e., the mapping between visual input on the retina to neuronal activations in visual cortical areas, is one of the central topics in visual neuroscience [1]. Although initially discovered in neurophysiology [2], Blood Oxygenation Level Dependent (BOLD) functional magnetic resonance imaging (fMRI) provides an excellent noninvasive tool to perform in vivo retinotopic mapping on human observers [3,4]. The procedure typically consists of (1) recording BOLD fMRI cortical activations generated by carefully designed slowly moving visual stimuli on the retina to uniquely encode the visual space [5,6], (2) co-registering the fMRI activations with structural MRI, (3) flattening the 3D structural MRI data to obtain visual cortical surface, and (4) decoding the retinal visual coordinates underlying the observed fMRI time series at each location on the cortical surface [3,7]. In the last two decades, the development of a variety of retinotopic mapping techniques has greatly extended our understanding of visual cortical organization in normal and abnormal populations [3,8-13].

Although the exact shape of visual inputs on the retina is not preserved on the cortical surface, neurophysiology studies [14-17] have shown that the mapping preserves local neighborhood geometric relationships. More specifically, neighboring points on the cortical surface should have neighboring retinal visual coordinates [18]. We adopted the term topology preserving (or topological) to refer to the preservation of neighborhood relationships of the observations in the output space [19]. In this work, the topological condition or topology-preserving are not related to the surface genus, handles or holes. The concept of topological condition is illustrated in Fig 1 for retinotopic mapping. Fig 1A and 1B show an example of topological mapping that preserves the structure of the source polygon (on the cortical surface) in the target (retina) domain (Fig 1B). Fig 1A and 1C show an example of non-topological mapping that breaks the source polygon structure after mapping: f is out of the polygon edge , resulting in what is called a flipped triangle Δfff (the orientations of the triangle are different in the source and target domains).

Fig 1

Illustration of topological and non-topological maps.

(A) A polygon patch on an inflated cortical surface (source domain). (B) A topological mapping of (A) in the target domain (retina space). It preserves the structure of the source polygon. (C) A non-topological mapping of (A). It breaks the structure of the source polygon: f is moved out of the edge , resulting in what is called a flipped triangle Δfff. (d) A typical retinotopic map of V1. Flipped triangles are labeled with red color.

The topological condition is a natural consequence of the hierarchical organization of the visual system: there cannot be duplicated neuronal representations of the same retinal location within each visual area (otherwise the visual area should be further divided into more subareas). Although neurophysiological studies on animals have found that the retinotopic maps in lower visual areas are topological [2], the “raw” retinotopic maps derived from BOLD fMRI recordings are usually not because of the low signal-to-noise ratio and spatial resolution of the technology [3]. Although the quality of fMRI data has improved during the last three decades, the signal-to-noise ratio is still relatively low even with ultra-high field MRI systems [20]. A typical raw retinotopic map of V1 from the state of the art dataset and method [20] still contains about 20% flipped triangles (Fig 1D; Table A in S5 Text). Violations of the topological condition not only make the results from fMRI inconsistent with those from neurophysiology but also make it difficult to perform quantitative analysis on retinotopic maps.

Typically, fMRI signals are analyzed on a voxel-by-voxel basis. The most prominent method is based on the population receptive field (pRF) model, which decodes perception parameters of the fMRI signals at each voxel. However, influenced by low SNR and low spatial resolution, the retinotopic maps from the analysis are typically not topological. In the past decade, significant efforts have been devoted to improving their accuracy and stability. For instance, the reliability of perception size was discussed [12,13], and methods that could decode perception parameters by model-free back projections with potentially more precise perception shape were proposed [21-23]. Although they have greatly improved the pRF model, none of the methods has explicitly investigated the topological condition.

In addition, a number of methods have been proposed to correct or reduce topological violations in fMRI-based retinotopic maps from the pRF analysis. Yet none of them has fully solved the problem. For instance, traditional spatial smoothing methods, including Gaussian smoothing [14,24], Laplacian smoothing [25], mesh-spectrum-based smoothing [26], cannot guarantee the topological condition because retinotopic maps are defined on two-dimensional surfaces, yet these methods only smooth the maps along one dimension (e.g., polar angle) at a time. Another promising approach registers the noisy “raw” retinotopic maps to an ideal template and assigns coordinates based on the registration function/morphing [8]. Although it can generate smooth and topological retinotopic maps, the method relies on diffeomorphic registration of a noisy parametrized surface to a predefined template, which is not easy to optimize, especially when the data are noisy. Furthermore, even if the registration issue can be solved, defining the right template remains a challenging problem. Other methods based on fitting predefined algebraic models (such as Schwartz’s complex “log” model [27] or Schira’s “Double-Sech” model [9]) can guarantee topological condition. However, the models introduce a lot of assumptions on the data [8,28]. The topological condition was considered in our previous work [29], but it was only applied to V1 with limited validation. And, because the manually drawn boundaries were not accurate, the results were not very precise near the boundary.

We propose a topological smoothing method to generate smooth and topological retinotopic maps for multiple visual regions without assuming algebraic models. A geometric concept, Beltrami coefficient (BC) [30], is adopted to formulate the optimization problem to define, quantify, and ensure the topological condition [31]. The Beltrami coefficient is a complex number that can quantify triangle-wise geometric conformality (i.e., the angle preserving property) and topological condition for surface-to-surface mapping. When we consider the topological condition in the mapping from the polygon in Fig 1A to 1B, the overall topological condition can be checked for each triangle attached to point P. Namely, one may check whether its mapping result, f, is moved out of any of the polygon edges, , and . If f is moved out of any of the edges, e.g., , the mapping is not topological. Without losing generality, when considering the topological condition for polygon edge , we can construct a coordinate system (Fig 2A and 2B). In Fig 2B, we illustrate three mapping results for P: Case (1) if f coincides with the green dot such that triangles ΔPPP and Δfff are similar (ΔPPP~Δfff), the magnitude of the Beltrami coefficient associated with triangle ΔPPP is zero (|μ| = 0). This is the case of conformal mapping. Case (2) if f is above the f(1) axis (f(2)>0), shown as the yellow dot, the magnitude of the Beltrami coefficient associated with triangle ΔPPP is greater than zero but less than one (0<|μ|<1). This is the case of quasiconformal mapping with respect to edge . And Case (3) if f is below the f(1) axis, shown as the red dot, the magnitude of the Beltrami coefficient associated with triangle ΔPPP is greater than one (|μ|>1) and the mapping is no longer topological. In summary, the complex-valued Beltrami coefficient encodes relative mapping information for each triangle and can be used to quantify the topological condition [32,33].

Fig 2

Beltrami coefficient and topological condition: (A) source triangle (triangle ΔPPP), (B) three target triangles Δfff from three different types of mappings, with different Beltrami coefficients.

To ensure the topological condition, we defined an objective function such that, after optimally adjusting the retinal visual coordinates of each point on the cortical surface, the maximum magnitude of the Beltrami coefficients on the new retinotopic map is less than one (i.e., topological). An efficient numerical method was then developed to solve the optimization problem.

Our proposed method can work on multiple visual areas simultaneously without assuming any prior model and allow us to draw more precise boundaries between visual areas. To validate the method, we conducted experiments on both synthetic and real retinotopic map data from the Human Connectome Project (HCP) [34] and compared the performance of our method with that of existing methods based on the quality of fits to the raw fMRI time series.

Methods

HCP retinotopy data

The Human connectome project (HCP) [34] provides a large publicly available retinotopy dataset collected on 7T MRI scanners. The data collection, conducted on 181 healthy young adults (22–35 years; 109 females, and 72 males) with normal or corrected-to-normal visual acuity, involved carefully designed retinotopy stimuli and resulted in a substantial amount of fMRI data (30 min, 1,800 time-points) acquired at very high spatial and temporal resolutions (1.6 mm isotropic voxels, 1-second temporal sampling). The dataset provides an exciting opportunity to evaluate our method.

Overview of the pipeline

The overall proposed processing pipeline is shown in Fig 3. In the HCP dataset, the visual stimuli consist of rotating wedges, expanding/shrinking rings, and moving bars [3,35]. The carriers of the stimuli are made of dynamic color textures that can better activate visual neurons (Fig 3A). We denote a point in the visual field by = (v(1),v(2))∈ℝ2, where v(1) is the eccentricity (distance to the fovea in degrees of visual angle), and v(2) is the polar angle relative to the positive horizontal line (Fig 3B). Both high-quality structural MRI and fMRI scans were acquired by the HCP group (Fig 3C and 3D). The processing pipeline consists of fMRI preprocessing [36], compressive population receptive field decoding [3], and topological smoothing. The raw retinotopic maps from pRF decoding (RM; Fig 3G) contain many topological violations (Fig 3H). Topological smoothing is the proposed method in this paper to fix topological violations and smooth the retinotopic maps (Fig 3I). In the following subsections, we briefly describe fMRI preprocessing and pRF decoding (although we used the publicly available decoded (raw) retinotopic maps from [20]), and then introduce the topological smoothing method. We list the key definitions and symbols in Table A in .

Fig 3

The processing pipelines.

(A) Visual stimuli. (B) Coordinate system of the visual field. (C) Structural (anatomical) MR image. (D) BOLD fMRI scans. (E) Preprocessed fMRI signals. (F) Reconstructed cortical surface with projected fMRI activations. (G) Retinotopic maps from the pRF model. (H) Level set of the raw retinotopic maps. (I) Level set of the smoothed retinotopy maps.

fMRI preprocessing

The goal of fMRI preprocessing is to detect, in a robust, sensitive, and valid way, the time series of brain activations of the voxels on the visual cortical surface that are associated with the visual stimuli. The HCP preprocessing pipeline [36] is standardized. First, a cortical surface is extracted from structural MRI using Freesurfer [37] and resampled (Fig 3F). Then, the raw fMRI data from each imaging session are co-registered across time to reduce the influence of head movements and other motion artifacts (Fig 3E). Finally, the co-registered fMRI data are projected onto the cortical surface. During the projection, spatial smoothing is applied along the cortical surface to improve the quality of the fMRI signal. The preprocessing eventually generates a resampled cortical surface as well as the fMRI activation time series on each point of the cortical surface.

pRF decoder

We briefly describe the population receptive field model (pRF) [3,38] to introduce some necessary notations. For readers familiar with pRF model, the visual coordinates we used here are related to the x-y coordinates by x = v(1) cos v(2),y = v(1) sin v(2). For each voxel P = (X,Y,Z)∈ℝ3 on the visual cortical surface, pRF is used to determine its receptive field, including its center location and size σ in the visual field. Namely, the collection of the mapping between P and (,σ) constitutes a raw retinotopic map. Assuming that the voxel’s population receptive field model is r(′;,σ) and the hemodynamic function is h(t), the predicted fMRI signal of the voxel can be written as: where β is the activation level, n is the exponent of the power function. We used the standard population receptive field model, i.e., , where x = v(1) cos v(2),y = v(1) sin v(2) (similarly for x′,y′). The center of the receptive field and population receptive field size σ can be estimated by minimizing the Least Square difference between the measured and predicted fMRI activation levels: where y(P) is the fMRI signal at voxel P. Iterations of this procedure across all the voxels on the visual cortical surface generate the raw retinotopic map. The quality of fit of the pRF model for voxel can be evaluated by .

Topological smoothing

Quasiconformal mapping and Beltrami coefficient

Most surface mappings can be modeled by quasi-conformal mapping, a generalization of conformal mapping. Specifically, the local deformation of a quasiconformal map can be characterized by its associated Beltrami coefficient [39]. Mathematically, f:ℂ→ℂ is quasiconformal if it satisfies the Beltrami equation, for some complex-valued Lebesgue measurable function μ satisfying ‖μ‖∞<1. The complex-valued is called the Beltrami coefficient of f, with u = u(1)+iu(2) and .

Intuitively speaking, if ‖μ‖∞<1, the map preserves the orientations of all triangles and therefore preserves the topological condition. However, if ‖μ‖∞>1 the topological condition is violated because the orientations of some triangles are flipped.

Surface flattening

Let f: v ↦ P be the retinotopic map between the retina and V1 cortical surface. To use the tool of the Beltrami coefficient, we need to flatten the cortical surface from 3D to a 2D unit disk. Because topological proprieties are transitive (that is if mappings f and c are both topological, the composed mapping f∘c is also topological) and mutual (if mappings f is topological, f−1 is also topological), we can reformulate the topological condition by mapping the 3-dimensional cortical surface to a 2D parametric space with a topological transformation [37,40]. Here, we performed a conformal mapping (angle preserving) of the visual cortical surface to a parametric space c: P ↦ , where . Although it is angle preserving and does not introduce any angle distortions, conformal mapping may introduce big metric distortions. To reduce metric distortions, we cut a geodesic disk patch that contains the region of interest from the cortical surface by: (1) picking a point on the cortical surface that roughly corresponds to the fovea as the center of the disk, (2) calculating the geodesic distances from all points on the cortical surface to the center of the disk [41], and (3) keeping the points of the cortical surface whose geodesic distances to the center of the disk are within a certain value. The geodesic disk patch (the region-of-interest—ROI—shown in gray color in Fig 4B) is mapped to a unit disk (Fig 4C). We refer the readers to [31,42] for conformal mapping and [43] for the computation of geodesic distance.

Fig 4

Illustration of (A) the visual field space, (B) the cortical surface space, and (C) the parametric space.

In this work, we focus on the 2D retinotopic mapping f = f−1∘c−1 instead of the 3D retinotopic mapping f.

Mathematical model

With the conformal parametrization in the last section, we use the Beltrami coefficient μ to model topological smoothing by finding new visual coordinates such that, where w is the weight function, ∇ = (∂/∂u(1),∂/∂u(2)) is the gradient operator, s is a positive scalar emphasizing the smoothness and is the Beltrami coefficient of written as (see ): In Eq 5, we have rewritten a point in the 2D parametric domain as a complex number u = u(1)+iu(2) and the corresponding visual coordinate as . Because we will only discuss the problem in 2D in the subsequent sections of the paper, we interchangeably regard each complex number as a position vector or vice versa, that is, u = u(1)+iu(2) is regarded as a complex number or as a position vector = (u(1),u(2)). We denote a position vector with a bold and a complex number with an italic u. The same notation applies to .

The Beltrami of can be used to quantify the Beltrami of f−1 because the Beltrami coefficient does not change with the conformal mapping, that is, if c is a conformal mapping (see for mathematical details). Intuitively, the Beltrami coefficient measures the angle distortion in a mapping (see ), incorporating a conformal mapping c does not introduce any angle distortion and therefore does not change the original Beltrami coefficient. We can directly estimate the Beltrami coefficient of the retinotopic mapping in the 2D parametric domain on the unit disk based on this property.

Model solver

To solve Eq (4) efficiently, we divide the problem into two subproblems, namely smoothing and topological projection. Specifically, during smoothing, we ignore the topological condition and find a smoothed map by, which can be solved by Laplacian smoothing [44,45] efficiently. For the second subproblem, we ignore the smoothness term and try to make the results topological.

Since we require the results to be both topological and smooth, we solve the subproblems iteratively. Namely, the output of topological projection is fed to smoothing, and so on so forth.

The key challenge is to solve Eq (7) efficiently. Because the topological condition is formulated with the Beltrami coefficient (), we need tools to compute the Beltrami coefficient and recover the function after processing the Beltrami. In the following, we describe the Beltrami coefficient computation when a function is given and the function recovery when the Beltrami coefficients are given in the discrete setting. Then we elaborate on the details about solving each subproblem.

Model discretization

We now convert the problem to the discrete domain because a cortical surface, S, is typically described by vertices V and triangular faces F: S = (F,V). The visual field is also discretized. Retinotopic mapping can be reformulated as finding the corresponding visual coordinate v in the visual field for every vertex P∈V. We denote W as w(v) to emphasize/deemphasize voxels with high and low quality of pRF fits.

is the coordinate to be solved, is the Beltrami coefficient at point u.

Discrete Beltrami coefficient computation

Given the analytical form of a mapping function f, we can directly compute the complex-valued Beltrami coefficient μ according to Eq (5). However, in the discrete case, the function is only defined on vertices, i.e., we only know the mapping function f that maps the vertices on the unit disk to the vertices on the visual field (Fig 5A), that is, v1 = f(u1), v2 = f(u2), and v3 = f(u3). To approximate the derivatives, f is interpolated linearly on each triangle, i.e., for any point u within the triangle, f(u) = B1()v1+B2()v2+B3()v3. The coefficients, B1, B2, and B3 are called the barycentric coefficients. Specifically, B1 (similarly for B2 and B3) is the ratio between the areas of triangles Δuu2u3 and Δu1u2u3: B1(u) = Area(Δuu2u3)/Area(Δu1u2u3). Now we can compute f’s partial derivatives to u(1) and u(2) and therefore the Beltrami coefficient μ according to Eq (5) (see for the explicit expression). Clearly μ is a constant within each triangle, as we approximate the mapping function linearly within each triangle.

Fig 5

Illustration of a mapping function and the divergence computation.

(A) An illustration of the mapping function in the discrete domain, and (B) the divergence approximation for a vertex.

Recovering mapping with Beltrami coefficients

We now explain how to recover function f = f(1)+if(2) with the given Beltrami coefficient function μ = ρ(u)+iτ(u), first introduced in [32]. According to the definition in Eq (5), we have:

Based on its real and imaginary components, Eq (9) is equivalent to: where , and . Applying on Eq (10A) and on Eq (10B), and then adding them together, we can simplify Eq (9) as: where is the gradient of f(1), and ∇ is the divergence operator (∇· = ∂G(1)/∂u(1)+∂G(2)/∂u(2) for vector = (G(1),G(2)) = A∇f(1)). By solving Eq (11) with Dirichlet boundary conditions (i.e., f(1) values on the boundary of the cortical patch), we can find a unique solution of f(1). Similarly, f(2) can be written as ∇·A∇f(2) = 0 and be solved uniquely.

In the discrete case, the mapping function is interpolated on each triangle. The gradient operator ∇f(1)(u) can be written out and is a constant vector on each triangle. The divergence ∇· is approximated on the dual of each vertex (a convex cell consisted of the circumcenters of the neighbor triangles). Specifically, consider a vertex with its neighbors (Fig 5B). We assume that ∇· is a constant on the dual of the center vertex, i.e., polygon D, enclosed by circumcenters of the neighbor triangles. Namely, we approximate the term with . According to the divergence theorem [46], the divergence ∫∇·dσ can be written as the integral on the boundary, i.e. , where is the dual cell unit normal. Since is a constant vector on each triangle (in the discrete case), the integration can be written as: where is the value of within triangle T, and = (−)⊥ denotes a vector perpendicular to − with the same norm as −. Substituting Eq (12) into Eq (11), we have a set of linear equations for each f(1)(u) and its neighbors. Finally, we can write ·A∇f(1) =0 in a matrix form (see ) and solve f(1) efficiently with the Dirichlet boundary condition. Similarly, f(2) can be solved efficiently.

Next, we describe how we can optimize Eq (7) by solving two subproblems iteratively.

Solving the first subproblem: Laplacian smoothing

The first subproblem can be solved by Laplacian smoothing of f(1) and f(2) separately. To get a smooth scalar from a noisy f(1), we solve the following Laplacian smoothing equation [44,45]: where the ∇ is defined in the parametric domain, i.e., ∇ = (∂/∂u(1),∂/∂u(2)). To find the solution, we set the partial derivatives of Eq (13) to zero. It leads to the following equation . It can be solved efficiently with linear algebra because ∇2 = ∇·∇ can be written in a matrix form (the details are provided in ), and σ is optimally chosen to minimize the generalized cross-validation score [47]. Similarly, f(2) can be smoothed to . Since the Beltrami coefficients quantify local deformations, Laplacian smoothing on visual coordinates will smooth the Beltrami coefficients. If most triangles are correctly orientated, the smoothing will also reduce the number of flipped triangles.

Solving the second subproblem: Topological projection

The second subproblem is to make f topological. It means that the associated Beltrami coefficients must satisfy the condition ‖μ‖<1 on the entire cortical surface. To ensure topological condition, for any mapping whose |μ| is greater than 1, we scale its magnitude by μ′ = μ/(|μ|+ϵ). The procedure is called “topological projection” [48,49] or “chopping” [50]. The goal is to find a topological mapping that is close to the original non-topological mapping. The new μ′ shares the same argument of μ, and|μ′| is less than 1 when we choose a small positive ϵ. Because the Beltrami coefficient uniquely encodes the mapping [31,51], the projection process changes the mapping when we solve equation Eq (11). The process of making the mapping topological can be illustrated with the example in Fig 2B, in which the target f is below the horizontal axis, and the magnitude of the associated Beltrami coefficient is greater than one (). If we shrink the Beltrami coefficient for Δfff and retain the Beltrami coefficients for all the other triangles, the new Beltrami coefficient is less than 1. When we solve the mapping function according to the new Beltrami coefficient (Eq (11)), f is moved and placed above the horizontal axis, making the triangle Δfff topological. Note that other target points, f, f, and f, also move to maintain the Beltrami coefficients of their associated triangles, e.g., when f moves up, f moves accordingly so that the Beltrami coefficient of triangle Δfff stays the same. In summary, given an initial non-topological mapping, we can project the Beltrami coefficient and solve Eq (11) to obtain a topological mapping that is close to it.

Removing phase jumping

We adopt the widely used polar coordinate system for the visual field, i.e., eccentricity r = v(1) and polar angle θ = v(2)∈[0,2π). The polar coordinate system has many advantages, but the polar angle θ is not continuous numerically near the x axis even when the mapping itself is. As shown in Fig 6A, points v1 and v2 share the same eccentricity and are both very close to the x axis, but their polar angles are very different because the polar angle jumps from 2π to 0 near the positive x axis. We call the phenomena “phase jumping”. It makes the coordinates around the x axis discontinuous.

Fig 6

Removing phase jumping and phase reversal.

(A, B) Phase jumping and the shift scheme used to remove it. (C) Phase reversal and the transformation used to remove it.

Phase jumping is not an intrinsic property of retinotopic mapping. The issue can be addressed with a transformed coordinate system: If we rotate the polar coordinate system, the phase jumping axis can be shifted. To avoid phase jumping, we shift the zero polar angle to the negative x axis for the left hemisphere (Fig 6B). Now phase jumping occurs near the negative x axis. Because retinotopic mapping of the left hemisphere is in Quadrants I and IV, the shift allows us to avoid phase jumping. For the right hemisphere, retinotopic mapping is in Quadrants II and III; we do not need any shift. By doing this, we can avoid unnecessary visual coordinate discontinuities.

Removing phase reversal

We apply another linear transformation on the polar angle to remove phase reversal across visual areas. Phase reversal is an important intrinsic property of retinotopic mapping. As shown in Fig 6C, if one walks from V3d to V3v along the yellow arc, the corresponding polar coordinate v(2) on the visual field first increases and then decreases. We assign a visual area with a positive (negative) visual field sign if the polar coordinate increases (decreases) when we move clockwise on an equal eccentricity curve, e.g., the yellow arc in Fig 6C. A phase reversal across the boundary of two cortical areas causes a change of visual field sign and therefore breaks the topological condition for the points near the boundary between visual areas. To make our proposed smoothing method applicable to multiple visual areas simultaneously, we apply another linear transformation to remove the visual field sign changes across visual area boundaries. The combination of two transformations, one for removing phase jumping and the other for eliminating visual field sign changes, is summarized in Table 1. We call the transformed polar angle “extended polar angle”. Because polar angles and visual area labels have a one-to-one correspondence with the extended polar angles, one can recover the original polar angles and visual area labels from the extended polar angles. With accurate extended polar angles, one can draw accurate interior boundaries between different visual areas (V1 and V2, V2 and V3, etc).

Table 1

Polar angle transformations are used to remove phase-jumping and changes of visual field sign for the left and right hemispheres in multiple visual areas.

Region Name	Transformation of Polar Angle θ
	Left Hemisphere	Right Hemisphere
V1v	T(θ) = θ+π	T(θ) = θ
V1d	T(θ) = θ−π	T(θ) = θ
V2v	T(θ) = −θ+2π	T(θ) = −θ+π
V2d	T(θ) = −θ−2π	T(θ) = −θ−π
V3v	T(θ) = θ+2π	T(θ) = θ+π
V3d	T(θ) = θ−2π	T(θ) = θ−π

In addition, these transformations make the smoothing more precise, especially near the boundaries of visual areas, because they remove the sharp polar angle changes across the boundaries.

Algorithm

We now summarize the overall topological smoothing algorithm in Alg. 1. The input data consist of the conformal parameterization of the cortical surface, and their associated visual field coordinates . The output is the updated retinotopic map coordinate . The algorithm is implemented in MATLAB 2019b [52,53].

Algorithm 1. Topological Smoother

Input Data: Retinotopic coordinates = {}, conformal parametrization = {} of the cortical surface, initial boundary retinotopic coordinates, , and boundary value change tolerance ϵ.

Output: Smoothed retinotopic coordinates

Correct the visual coordinates according to Table 1 to get adjusted visual coordinates ′

assign the initial retinotopic coordinates for each vertex on the unit disk

repeat

Compute Beltrami coefficient μ for the mapping using Eq (5).

Project non-topological mappings with μ′ = μ/(|μ|+ϵ), if |μ|>1.

Compute new mapping by LBS for μ′ with boundary value .

Apply Laplacian smoothing on and , respectively) to get .

Fit a linear function near the boundary and update the boundary value to ′

Check and restore each point in ′ if it deviates from the initial value by more than ϵ.

until max|μ|<1

return.

Results

Synthetic data

We first evaluated the performance of the algorithm on synthetic data generated with the following function,

We chose ln function (k = 0.5 in this paper) as it is a good approximation of the retinotopic mapping from the visual field to the flattened cortical surface [54]. Note that in the original study, Schwartz used the model to map from the visual field to the parametric coordinates of a flattened cortical surface. It is the inverse function that we studied in this paper. We followed the convention only in this subsection, i.e., generating a smoothed from a noisy u(v), instead of a smoothed from a noisy v(u). Note that, because the topological condition is reciprocal, the conclusion based on this approach still applies to the inverse function. More specifically, the visual field grid mesh is first created (Fig 7A) with the following configurations: eccentricity spans from 0° to 4.5° and polar angle ranges from −π/2 to π/2. We then sampled 12 points in each direction and formed the quadrilateral grid mesh (144 points in total). Based on the quadrilateral mesh, we created the triangular faces by connecting diagonal points for each quadrilateral. Then we mapped the triangular mesh to the parametric space (Fig 7B), according to Eq (14). Ideally, the mapping points should be located on the grid (Fig 7B). We then added a small amount of noise with a peak signal-noise ratio (PSNR) of 10 (Fig 7C), and a large amount of noise with a PSNR = 5 (Fig 7D).

Fig 7

Illustration of the synthetic data.

(A) The visual field domain. (B) Mapping without noise. (C) Mapping with a small amount of noise (PSNR = 10). (D) Mapping with a large amount of noise (PSNR = 5).

We applied several smoothing methods, including Average Smoothing, Median Smoothing, and Laplacian Smoothing [44], as well as our proposed Topological Smoothing, to the synthetic mapping data in both high and low signal-to-noise cases. All the existing methods are spatial smoothing procedures. Intuitively, Average Smoothing takes the average around the neighborhoods of each vertex. Median Smoothing takes the median value of the neighborhood. All existing smoothing procedures treat the u(1) and u(2) coordinates separately. In this section, we fixed the smooth parameter s = 0.001 (Eq (7)) and the boundary values were initially inferred from the Average Smoothing results and boundary change tolerance is ϵ = 0.01. We list the performance metric, i.e., the deviation from the ground truth of all smoothers in Tables 2 and 3. Value deviation is the average difference relative to the ground truth. Angle distortion (the angle spanned by the tangent vector of u(1)’s contour curves and the tangent vector of u(2)’s contour curves, at the same location; see for details) quantifies the local angle changes (in degree unit) after the mapping, and the number of flipped triangles measures violations of the topological condition. All the experiments were repeated 50 times to reduce fluctuation caused by random noise. The number of flipped triangles is the median of the 50 experiments. In addition, we ran a non-parametric permutation test [55] to conduct pairwise comparisons of the results from the proposed method and each of the other methods (including no smoothing). The null hypothesis is that the results from the proposed smoother are not significantly different from those of the other methods. We report the p value, i.e., the probability that the null hypothesis was true. In brief, if the p value is less than 0.01, the results from the proposed smoothing method are statistically different from those of the other methods.

Table 2

Comparison of different smoothing methods based on three metrics: value deviation, angle distortion, and number of flipped triangles with a small amount of noise (PSNR = 10).

Method	Value Distortion			Angle Distortion in deg.			Topology
	Mean	Std	p	Mean	Std	p	T^	p
No Smoothing	0.235	0.13	0.00001	45.001	25.25	0.00001	54	0.00000
Average	0.134	0.09	0.00001	30.903	22.52	0.00005	10	0.00000
Median	0.170	0.10	0.00001	40.124	24.04	0.00001	19	0.00000
Laplacian	0.558	0.48	0.00001	41.432	25.16	0.00001	18	0.00000
Proposed	0.143	0.09	------	18.313	16.41	-----	0	-----

Table 3

Comparison of different smoothing methods based on three metrics: value deviation, angle distortion, and number of flipped triangles with a large amount of noise (PSNR = 5).

Method	Value Distortion			Angle Distortion (degree)			Topology
	Mean	Std	p	Mean	Std	p	T^	p
No Smoothing	0.334	0.18	0.00001	47.067	25.87	0.00001	65	0.00001
Average	0.169	0.10	0.00001	36.099	24.09	0.00040	22	0.00001
Median	0.206	0.12	0.00001	43.314	24.76	0.00001	32	0.00001
Laplacian	0.604	0.46	0.00001	44.544	26.18	0.00001	25	0.00001
Proposed	0.169	0.11	------	23.226	20.70	------	0	-----

We found that (1) all smoothing methods reduced value deviation, angle distortion, and number of flipped triangles (Tables 2 and 3); (2) the proposed method achieved the best topological result and the smallest angle deviation (all p<0.00001); (3) when the noise was relatively small (Table 2), violations of the topological condition were largely corrected by most of the methods, and the average smoothing method achieved very good topological results; (4) when the noise was relatively high (Table 3), the average smoothing method lost its ability to correct topological violations; and (5) the proposed method always performed best in angle preserving, which was not a surprise as the magnitude of Beltrami coefficient |μ| is related to angle distortion of the mapping (see ), and we penalized big |μ|′s in the proposed method (Eq (7)) to reduce angle distortions. To summarize, the proposed smoother performed best to correct topological violations and reduce angle distortions. It also performed better or comparably with all the other methods in terms of visual coordinate distortions. Fig 8 illustrates the smoothing results for all the methods when the noise was relatively low (PSNR = 10).

Fig 8

Smoothing results for the synthetic data: (A)-(D) are the results with PSNR = 10 for (A) Average smoothing, (B) Median smoothing, (C) Laplacian smoothing, and (D) proposed topological smoothing; (E)-(H) are the results with the same order as (A)-(D) but with PSNR = 5.

Synthetic Data II

We performed another set of synthetic data experiment in which noise was directly added to the fMRI signals. We started with the ground-truth parameters: receptive field centers spanning the right visual space [0:0.08:6.4]°×[0,0.08: 4.8]°; receptive field size = 0.4°; gain = 1; and exponent of the power function n = 0.5. Then, we followed Kay’s pRF model (AnalyzePRF [38,56] with stimulus from the HCP group, 200×200 in spatial resolution with 1800 seconds) to generate the ideal fMRI signal. We further normalized the fMRI signal to unit variance, added a specific amount of white noise, decoded the pRF model for the fMRI signal (with noise), and computed the number of flipped triangles. We varied the noise level until the percentage of flipped triangle was about 20%. (Specifically, in our experiments, noise with a 2.5 standard deviation induced about 22% flipped triangles). With pRF decoding results at this noise level as input, we applied the topological smoothing method and generated new visual coordinates. Then we compared our work with the pRF model and several post-processing methods (Table 4).

Table 4

Comparison of different smoothing methods based on noisy fMRI decoding with about 22% of flipping triangles.

Method	Center Distortion			Angle Distortion (degree)			T^	R 2
	Mean	Std	p	Mean	Std	p
No Smoothing	0.407	0.520	0.352	48.55	30.20	0.00	32	13.28
Average	0.311	0.238	0.070	23.35	23.16	0.00	15	13.41
Median	0.298	0.176	0.010	22.40	16.57	0.00	11	13.41
Laplacian	0.429	0.509	0.149	42.82	29.48	0.00	23	13.09
Proposed	0.361	0.231	------	13.84	12.87	------	0	13.27

We found that (1) the proposed method can fully ensure the topological condition, with 0 flipped triangles, while other smoothing methods had at least 6% of flipped triangles; (2) the proposed method had minimal angle distortion; and (3) all smoothing methods improved the receptive field center estimates, with the median smoothing method achieving minimal center distortion. The last point is reasonable because the current work is mainly focused on the topological condition. Future fine-tuning of the pRF parameters may alleviate this problem. The new synthetic experimental results further demonstrated that the proposed method could remove topological violations and reduce angle distortions.

V1 retinotopic map

Although 7-Tesla MRI systems have dramatically improved the signal-to-noise ratio, the retinotopic maps in the HCP dataset are still quite noisy and non-topological. We first applied our proposed smoother on the V1 retinotopic maps of all observers in the HCP dataset.

Since there is no ground truth, we evaluated the performance of the method by the goodness of fit of the predicted fMRI signals with the smoothed visual coordinates. Specifically, we first decoded the visual coordinates , population perception field size σ, BOLD response level g for the 181 observers by the pRF method [3]. Then we applied several smoothing methods on their V1 retinotopic maps. Based on the smoothed coordinates, we computed fMRI time series predictions at every vertex of the cortical surface. To evaluate the impact of changing the receptive field centers, we set the response level g and the exponent of the power function n the same as the original pRF in computing fMRI time series predictions.

We measured the goodness of fit of the smoothing methods to the fMRI time series (with the code from AnalyzePRF [38,56]). The results are listed in Table 5, with the following metrics: (1) the average number of flipped triangles across subjects, (2) the average percentage of the non-topological area Ar, (3) the average Normalized-Root-Mean-Square-Error (NRMSE) of the fits, and (3) the average variance accounted for across all the vertices in V1.

Table 5

Comparison of different methods in fitting the fMRI time series.

Method	T^	Arn(%)	NRMSE × 10−3	R^2¯
Before Smoothing	73	20.7	0.07757	28.23
Average Smoothing	28	8.1	0.07779	27.86
Median Smoothing	35	9.7	0.07773	27.92
Laplacian Smoothing	64	18.0	0.07763	28.09
Proposed Smoothing	0	0.0	0.08139	20.35

We observe that: (1) all the smoothing methods reduced the number of flipped triangles, while only the proposed method eliminated topology violations completely; (2) the percentage of the non-topological area Ar was about 8~20% of the total area with all but the proposed smoothing method; (3) The proposed method decreased R2, mainly due to the constrain of the topological condition. The number of flipped triangles for individual observers by all methods can be found in the supplementary data (see Table A in ). The results also show that the proposed method was always better in preserving topology while the goodness of fit to the fMRI time series was reduced from the original result.

For a more intuitive comparison, Fig 9 shows the original retinotopic map of the left hemisphere of the first observer in the HCP dataset and the smoothed versions by different methods. The original retinotopic map (Fig 9A) is not topological because there are many flipped triangles. Fig 9B–9E show the results from the Average, Median, Laplacian, and proposed smoothing methods, respectively. Although all smoothing methods improved the smoothness of the maps, topological violation still occurred in the results of all but the proposed smoothing methods, especially near the fovea. Only our proposed method can generate topology-preserving results (Fig 9E). Fig 9F shows the smoothed retinotopic map with level sets on the inflated left hemisphere. To enable a clear comparison, we present enlarged level sets on the inflated mesh for each smoothing method in Fig 9H–9K. The green curves from the left to right represent 0.5° to 6° eccentricities with a 0.5° interval, and the blue curves, from the bottom to up, represent -90° to 90° polar angles with a 15° interval. We can also clearly see topological violations in Fig 9G–9J, especially in the fovea region. Only the proposed method can generate topological and smooth results (Fig 9K). The same experiment was conducted on the data of all observers and the results are available in . These results clearly show that the proposed method can generate smooth and topological retinotopic maps in V1.

Fig 9

Results on the first observer: (A) the raw retinotopic map on the 2D parametric domain; (B)-(E) results from the Average, Median, Laplacian, and proposed methods, respectively; (F) the whole inflated left mesh; (G) the raw map on the cortical surface; (H)-(K) results on the inflated surface, in the same order of (B)-(E). In (A)-(E) the x axis is the visual x coordinate, and the y axis is the visual y coordinate.

Retinotopic maps of the V1/V2/V3 complex

We took the initial visual area labels from the multimodal surface matching (MSM) results of the HCP group [57,58], removed phase-jumping and phase-reversal for the left and right hemispheres in multiple visual areas after polar angle transformations, and achieved topological and smooth results across them (Fig 10). Fig 10A–10E show the visual coordinates in the eccentricity-extended-polar coordinate space, in which the x axis is the eccentricity, and the y axis is the extended polar angle (Table 1), for the raw data, average smoothed data, median smoothed data, Laplacian smoothed data, and topological smoothed data, respectively. All smoothing methods reduced topology violations, compared to the raw data. However, only our proposed smoothing method generated topological mapping such that the contour curves (the green and blue curves in Fig 10K) were smooth and had no self-intersection. We performed the same comparison across all observers (Table A in S5 Text) and showed that our proposed topological smoothing was the best in topology preserving.

Fig 10

The retinotopic mapping of the V1/V2/V3 complex of the left hemisphere of the first observer: (A) the raw retinotopic map in the eccentricity-extended-polar coordinate space (with 232 flipped triangles), (B)-(E) the smoothing results of the Average (with 55 flipped triangles), Median (with 85 flipped triangles), Laplacian (with 151 flipped triangles), and the proposed smooth (no flipped triangles), respectively, (F) the entire left inflated mesh in medial view, (G) an enlarged graph of the raw retinotopic map on the inflated mesh, (H)-(K) smoothing results on the inflated surface of the four smoothing methods, in the same order of (B)-(E). In (A)-(E) the x axis is the eccentricity v(1), and the y axis is the extended polar angle v(2), in (H)-(K), the green and blue curves are levels sets, i.e., the contours of eccentricity v(1) and extended polar angle v(2), respectively.

Boundary delineation

In the eccentricity-extended-polar coordinate space, the boundaries between different visual areas in the V1/V2/V3 complex are defined by specific polar angles. Fig 11A shows the boundaries directly inferred from the raw data of the first observer. Fig 11B–11E show the boundaries generated from the smoothed data of the same observer, after average, median, Laplacian, and topological smoothing. The new boundaries largely coincided with manually labeled visual areas from the average data, represented by the colors of the cortical surface Fig 11). The green curves (from the left to the right) are eccentricity contours from 1° to 6° with a 1° interval, and the purple curves represent boundaries of visual areas (V3d/V2d/V1d/V1v/V2v/V3v, from up to down). As in the retinotopic maps, topological violations occurred in boundary delineation from the raw data and after smoothing with all but our proposed method. We present all observers’ area boundaries after topological smoothing in They were all smooth and continuous.

Fig 11

Visual area boundary delineation of the V1/V2/V3 complex of the first observer.

The colors of the surface represent manually labeled visual areas from the average data. The purple curves indicate boundaries between visual areas of V3d/V2d/V1d/V1v/V2v/V3v (from up to down), respectively. The green curves are eccentricity contours. (A) boundaries inferred from the raw data, and (B)-(E) boundaries after average, median, Laplacian, and topological smoothing.

Discussion

In this study, we modeled the topological condition and generated topological and smooth retinotopic maps. The Beltrami coefficient, a metric of quasiconformal mapping, was used to define the topological condition. We developed a mathematical model to quantify topological smoothing as a constrained optimization problem and elaborated an efficient numerical method to solve it. The method was applied to V1, V2, and V3 simultaneously, and was robust to inaccurate boundary labeling. Experiments with both simulated and real retinotopy data demonstrated that the proposed method could generate topological and smooth retinotopic maps. As a result, we were able to improve boundary delineation. To our best knowledge, conventional methods have not fully considered topological constraints for multiple regions in retinotopic smoothing. Our novel algorithm made the retinotopic maps from BOLD fMRI topological and consistent with results from neurophysiology. Our work improved the quality of retinotopic mapping and built a solid foundation for future retinotopy related research.

Quality assessment

In this study, we enforced the topological condition in retinotopic mapping. This enforcement could potentially introduce large visual coordinate changes to the raw data compared to other smoothing methods. It is important to monitor whether the smoothing results are acceptable or have relatively big deviations. Such quality assessment can provide the right level of confidence when drawing conclusions based on the smoothing results. In Fig 12, the visual coordinate changes of two typical observers are shown on the cortical and parametric surfaces respectively. Fig 12A shows the studied visual areas on the 3D cortical surface. Fig 12B and 12C shows the distortion in visual space. The colors indicate the amount of visual coordinate changes for the proposed smoothing process. The first observer’s results (Fig 12A) suggested that the smoothing was generally good, with on average an about 1.12° change of visual coordinates. In contrast, the other observer’s results (Fig 12C) suggested that some sub-regions have relatively big visual coordinate changes, with on average an about 1.75° change of visual coordinates. We also show the mean visual coordination distortion distribution of all 181 subjects in Fig 12D. Together with Table 5, the quality assessment results ensure our retinotopic mappings well respects the original fMRI signals.

Fig 12

Visual coordinate change after smoothing.

(A) The visual change is rendered on the inflated cortical surface (B) The visual change is rendered on the parametric surface for the first observer. Similarly, (C) is the visual change for the second observer. (D) is the histogram of average visual change for the HCP dataset (N = 181).

Out-most boundary condition

In this work, visual areas, V1/V2/V3, and the initial boundaries between them were determined by surface registration. Following the HCP protocol [36], MSM [58] was adopted to register each observer’s cortical surface to the retinotopy atlas obtained from the HCP [34]. After the registration, the boundaries for each observer were determined along equal polar angle curves, and our extended polar angles were initialized with these boundaries. The out-most eccentricity boundaries were manually determined based on the eccentricity signal range (in our case, approximately from 0.2° to 7°).

The initial boundaries might not be accurate since they were determined by surface registration. With extended polar angles, our algorithm was robust enough to fix the interior visual area boundaries. Our method requires that the visual coordinates of the out-most boundaries are specified when solving the LBS. Even so, if the visual coordinates of the out-most boundary are not properly given, we can still update the out-most boundaries and recover interior visual coordinates. We tested two approaches to set the out-most boundary values. One approach is to set a fixed boundary for all observers to simulate the case in which a good boundary configuration is given, and the other approach is to adjust the out-most boundary values by fitting a smooth spline of the interior visual values to simulate the case in which the boundary coordinates are given but not precise. Each approach has its own pros and cons. The first approach is numerically stable but ignores individual differences, while the second approach adjusts the boundary values based on retinotopic data, which has higher accuracy but may fail to set reasonable boundary values. We compared both approaches in the experiments. Fig 13A and 13C are the results from the first approach, and Fig 13B and 13D are the results from the second approach. We found a relatively large difference on the out-most boundaries (red boundaries in Fig 13A and 13B), while the interior regions were less affected. The mean difference between the results of the two approaches was 0.23°±0.15°, which is relatively small. In this paper, we took a combined approach: the boundary values were set within a small tolerance range (0.5°) of the initial average boundary values. Nevertheless, we found that, although the out-most boundaries may not be precise, the interior results were not much impacted by the boundary values. Therefore, we can use the pre-defined boundary configuration for all the subjects.

Fig 13

Comparison of two boundary approaches.

(A) The result of fixed boundary visual coordinates. (B) The result of changeable boundary visual coordinates. (C) Result of (A) on the inflated cortex. (D) Result of (B) on the inflate cortex.

It is worth noting that our method is applicable to the registration results from Freesurfer [37] or MSM [58] or Bayesian register [8]. Here we applied it to successfully process data from all 181 HCP observers. In our code repository, we have published all the necessary manuals and atlas for others to reproduce our work. We also provided a simple software tool to manually draw the out-most boundaries on visual cortical surfaces and generate corresponding visual coordinates.

Robustness with inaccurate interior boundaries

When smoothing the retinotopic maps of the V1/V2/V3 complex, we applied linear transformations to the polar angles based on the visual region labels derived from the anatomy. Here, we explain why our proposed smoothing method was robust even with inaccurate boundary specifications.

We used the first observer’s data as an example to illustrate that the proposed method was robust to inaccurate interior boundaries (Fig 14). In Fig 14A, we show the raw visual polar angle data on the parametric disk for the first observer, with a white curve at 0.3-degree eccentricity. Then we applied a linear transformation before smoothing, according to Table 1. For a vertex with an inaccurate boundary label, the polar angle would be transformed with a wrong offset. Since the transformation adds or subtracts multiples of π as the offset according to its visual region label in Table 1, a wrong label makes the transformed value dramatically different from its neighboring value. Therefore, an incorrectly transformed vertex looks like a spike and destroys its topological condition (red points in Fig 14B).

Fig 14

Smoothing is robust to the inaccurate interior boundary: (A) the raw visual polar angle v(2) is rendered on the parametric unit disk (i.e., radius is one), (B) v(2) as a function on the cropped region, (C) visual polar angle v(2) along the white arc by topological smoothing, (D) similar result as (C) but the curve is from Laplacian smoothing.

In the next step, we applied the topological smoothing (smoothed visual polar angle is the slope-shaped surface in Fig 14B). Since it enforces topology conditions, our topological smoothing method replaces the spike with neighboring values. Namely, the proposed method can fix the wrong offset and is therefore robust to inaccurate interior boundaries.

We further verified that the robustness was not due to Laplacian smoothing. In Fig 14C, we plot the raw data and topologically smoothed results of the white curve. The white curve is formed by points with eccentricity v(1) = 0.3° with tolerance 0.05°. We now consider the data near the white curve. If the eccentricity measurement is accurate, the 2D topological condition can be considered as a monotonicity requirement. In Fig 14C, the blue dots are the raw extended polar angle data, and the green curve is the topological smoothing result. Even though there were outliers (the red circled points, mainly due to the noisy eccentricity data), the main portion of the trend was monotonic. In Fig 14D, we plot the results with only Laplacian smoothing without topology constraints as the green curve, which is crooked in order to fit the data but not monotonic (topological). This means Laplacian smoothing can smooth but cannot fix topological violations. In contrast, our topological smoother ensures topological condition, which essentially makes the mapping “monotonic”. It is the topological constraint rather than Laplacian smoothing that makes the smoothing more robust to inaccurate interior boundaries.

The new topology-projection module in the proposed method makes it robust to inaccurate interior boundaries. The smoothed extended polar angles from topological smoothing can be used to infer visual area labels.

We computed visual coordinate changes following an alteration of the V1/V2 boundary. If the visual coordinates did not change too much when the boundary was altered, the results were robust to interior boundaries. We started from the first subject’s data (Fig 15A). We expanded V2v by moving points in V1v into V2v if they were adjacent to V2v. The expansion was repeated to simulate different boundary misplacements (three rounds in Fig 15B and five rounds in Fig 15C). Then we applied our smoothing method and compared the results with no perturbation. The boxplot in Fig 15D shows that the visual coordinate differences were relatively small (the average is less than 0.5°) with V2v boundary expansion. The results suggest that the proposed method is robust to interior boundary placement.

Fig 15

(A, B, C) Perturbations of the V1/V2 boundary, where (B) was obtained after three rounds of expansion and (C) was obtained after five rounds of expansion. (D) Visual coordinate changes with different degrees of perturbations.

The trade-off between smoothness and quality of fit

The smoothing weight s in Eq (7) controls the trade-off between smoothness and quality of fit under topological condition. Although other studies have attempted to find the optimal s for linear smoothers [59], the topology constraints made it challenging. We used the following way to set the proper s: find the largest weight while ensure that the average difference between the smoothed visual coordinates and the raw data is within a certain range (1° in this work). To empirically determine the s value, we computed the percentage of the non-topological area and the Normalized-Root-Mean-Square-Error (NRMSE) of the fit to the fMRI time series for the raw data, when s = 0, and when s was set properly for the first observer. For the raw data: NRMSE = 0.07944×10−3, 24.7% of the area is non-topological; when s = 0, NRMSE = 0.08546×10−3, 0% of the area is non-topological; and when s = 0.0121, NRMSE = 0.8551×10−3, 0% of the area is non-topological. The NRMSE values for s = 0 and s = 0.0121 were quite close. Therefore, we prefer the latter result in this work.

Beltrami coefficient vs. Jacobian

In this work, the topological condition is ensured by enforcing ‖μ‖∞<1. Another geometric concept, the determinant of the Jacobian matrix J can also be adopted to enforce the topological condition. For example, it has been adopted by previous work to monitor the local visual field signs of retinotopic maps [14], which indicate the trend of visual polar angle changes (increasing or decreasing) on an iso-level eccentricity curve. Constraining the determinant of the Jacobian to be positive can also ensure the topological condition and visual field sign within a neighborhood. Mathematically, the Beltrami coefficient and the determinant of the Jacobian matrix are related via the following equation: [60].

We prefer the Beltrami coefficient because it is more suitable for retinotopic mapping. First, with the Beltrami coefficient, we can define the topological constraint as a penalty term to smoothing energy (Eq 7) because the Beltrami coefficient quantifies angle distortions and retinotopic mapping is considerably angle preserving [54,61]. In contrast, J cannot be added as a term in smoothing energy because it quantifies area changes, and retinotopic mapping is not area-preserving. Numerically, one must separately treat the constraint and adjust visual coordinates to ensure that J is postive. Second, as indicated by Theorem 1 in S2 Text, the Beltrami coefficient map uniquely determines a diffeomorphic mapping between two 2D surfaces, making it possible to recover/solve the mapping if we know the Beltrami coefficients. Therefore, we can divide the smoothing process into two subproblems and solve each separately and efficiently. To our knowledge, there are no efficient and stable methods to utilize Jacobian J to reconstruct diffeomorphic mappings. In summary, although it may be possible to use J to quantify the topological condition, the Beltrami coefficient is more suitable for modeling retinotopic maps.

Smoothness and topological projection

The proposed topological smoother produces topological and smooth retinotopic maps. Although these two conditions are related, the topological condition does not ensure smoothness, nor does smoothness ensure the topological condition. In fact, the topological condition can be ensured without smoothing. Empirically, no one expects retinotopic maps to be topological but extremely irregular and rough.

Smoothing does to some degree improve the topological condition. Although we do not know the exact relationship between visual coordinate smoothing and the topological condition (monitored by Beltrami coefficient), the Beltrami coefficient quantifies the local deformations based on its definition. Smoothing can reduce the outliers in Beltrami coefficients. Because 80% of the triangles in the HCP dataset are orientation preserving (See Table A in S5 Text), smoothing the visual coordinates may reduce the norm of the Beltrami coefficients in average and improve the topological condition. If the data have more flipped triangles (e.g., ~50% triangles’ are flipped), smoothing may introduce more non-topological triangles and make the problem worse.

The topological condition cannot be ensured by reiterative smoothing. As shown in Fig 15D, repeated Laplacian smoothing did not generate fully topological results. The topological condition is mainly fulfilled by topological projection. However, after fixing the topological condition, a major concern is the results could be too irregular or deviate too much from the original pRF solutions. The proposed Topological smoother can produce results that are both topological and smooth. Roughly speaking, if the fMRI signal to noise ratio is very high (say SNR >20), the pRF solutions would be smooth and topological. For example, the pRF solution of the averaged fMRI signals from many subjects (https://osf.io/bw9ec/wiki/home/) has very few flipping triangles in V1. This is why the energy function in the topological smoother has a term that penalizes non-smoothness and therefore may improve the fit to the data and improve the topological condition. We believe the approach can balance the topological condition and smoothness.

Extension to higher visual regions

Although we focused on the V1/V2/V3 complex in this study, our method can be in principle extended to higher visual areas. However, there are some additional challenges. As we move to higher visual areas, the signal-to-noise ratio of the retinotopic map is further reduced. It is extremely challenging for all smoothing methods to balance the topological condition and goodness of fit. As shown in Fig 16, we applied our proposed smoothing method on the V3B retinotopic maps of the first four observers in the HCP dataset. Instead of smoothing V3B with V1-V3 together, we performed separate smoothing on V3B and concatenated the results with those on V1-V3. Although topological smoothing was feasible in higher visual regions, extra care must be taken because of the low signal-to-noise ratio of BOLD fMRI activation.

Fig 16

Smoothing extended to the retinotopic map of V3B (the blue color region) of the first four observers.

Benefits of the extended polar angles

Our previous work on topology correction [29] was only applicable to V1; it could not be used to improve boundary delineation. In this work, by introducing extended polar angles, we extended the framework to the V1/V2/V3 complex and achieved topology-preserving smoothing in multiple visual areas simultaneously. Although our previous work [29] might be potentially appliable to V2 and V3, accurate V1/V2 and V2/V3 boundaries were necessary. In the current work, given the out-most boundaries of the V1/V2/V3 complex, our work not only generated more precise smoothing within the interior of each region, but also automatically delineated the boundaries between them. The new development on boundary delineation may significantly improve the robustness and accuracy of retinotopic mapping. In addition, our prior work was mainly evaluated with visual coordinate the changes. How much impact it had on the goodness of fit of the fMRI time series was not evaluated. Here we added fits to the fMRI time series in performance evaluation, which may better justify our work and alleviate the concern that the proposed method may over-smooth the retinotopic maps. With extensive evaluation based on both synthetic data and real retinotopic map data from 181 HCP observers, the current study demonstrated that our proposed method may improve the accuracy and stability of retinotopic mapping, especially the interior boundary conditions.

Elaboration on the reduced R2 in

The topological smoother can completely remove topology violations in retinotopic maps. However, as shown in Table 4, it reduced the measure of explained variance, R2. In general, the pRF model should have a better R2, since R2 values are computed at each point on the cortical surface. There are no constraints, nor penalty for non-topological results in R2. However, neurophysiology studies strongly suggested that retinotopic mapping is topological for normal subjects. If the pRF results are not topological, the problem is very likely caused by the low spatial resolution and SNR of the fMRI data, an issue that cannot be addressed with current MRI technologies. In this work, we assume that the neurophysiology results are correct and attempted to fix topological violations in retinotopic maps by slightly updating the visual coordinates. We used pRF results as input data without 1) fine tuning pRF model parameters, such as the gain or the exponent of the power function [62], and 2) iterating with the pRF model to search for an optimal solution with balanced fitting quality and topological condition. Our ongoing work [62] tries to integrate the pRF model and topological conditions by defining a combined energy term. We expect such a balanced approach may improve topological condition without sacrificing data fitting quality.

Relation to registration

The Beltrami coefficient (BC) encodes 2D-to-2D mapping up to a normalization [60]. Manipulating the Beltrami coefficients of the registration function between a subject and a template can also achieve diffeomorphic registration and generate topological results [63]. The benefits of this approach was discussed in [8]. However, one main issue is that the template itself might be biased. One potentially promising solution is to combine registration and smoothing: while registration provides good boundaries, topological smoothing provides good interior smoothing [62].

Conclusion and future work

We adopted the Beltrami coefficients to quantify the topological condition in retinotopic mapping, formed a concise but fundamental model, and provided an efficient numerical solver to generate smooth retinotopic maps. To our knowledge, the proposed topological smoother is the first method that guarantees the topological condition in retinotopic mapping. Our work improved the quality of the retinotopic maps and built a solid foundation for future research based on retinotopy.

Notations and definitions.

(DOCX)

Click here for additional data file.

Beltrami coefficient.

(DOCX)

Click here for additional data file.

Discrete operators on a 2D triangular mesh.

(DOCX)

Click here for additional data file.

Angle Distortion.

(DOCX)

Click here for additional data file.

Supplementary data and code.

(DOCX)

Click here for additional data file.

27 Jan 2021

Dear Dr Wang,

Thank you very much for submitting your manuscript "Topology-Preserving Smoothing of Retinotopic Maps" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

Please particularly pay mind to points on clarity and validation. As Reviewer 3 points out there is a previous short conference paper, which should be cited and extensions relative to that work should be clarified.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Emma Claire Robinson

Associate Editor

PLOS Computational Biology

Wolfgang Einhäuser

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The paper presents a method for retinotopic maps that preserve topology. Beltrami coefficent is used to define the topology.

* Why use Beltrami coefficent, surface Jacobian can be computed directly and that could be a more direct measure of distortion. In fact surface Jacobian has been used by the authors in the past.

* More description of Beltrami coefficent should be added to the paper.

* How does the the Beltrami coefficent relate to the topology and distortion should be described in more detail possibly with some illustrations.

* How does the smoothing in fMRI data affect the Beltrami coefficent?

* I think the notation of first fundamental form in eq 3 might not be necessary. If the authors want to keep it, more description of them should be added.

Overall, this is a very interesting work. Addressing the above issues would improve the paper.

Reviewer #2: The manuscript by Tu, Ta, Lu & Wang introduces and develops a new method for analysing retinotopic maps. Essentially, the method enforces a topological mapping from retina/visual field space to cortical space, by using the Beltrami coefficient. Turning a measurement of many independent and noisy estimates into a consistent mapping bears many advantages and this has been previously recognized by Doughtery et al 2003 and Schira et al. 2007, using a technically different approach of morphing a map template on retinotopic data. The procedure suggested here has many advantages and some disadvantages over previous methods. The approach suggested seems very elegant and promising, the implementation points to great skills, and the method could indeed significantly improve our analysis of retinotopic maps. However, I am struggling with the current manuscript both in clarity but also structure. I suggest the authors make the purpose clear and structure the manuscript accordingly. Importantly, measures and analysis such as CMF or quality assessment need to be significantly improved or removed from the manuscript. At present they seem like an afterthought and that does neither these assessments nor the method justice.

1. The authors should do a better job explaining their core assumption that is that the mapping between retina/visual field and cortex is topological. I feel Fig 1 is ill suited and instead confusing in this context. Upon reading Fig 1 (for example F1b and F1c) I was trying to understand what topological concept the six nodes and eight edges in F1c violate, for example they could form a handle or a hole in the surface. However, that little set of triangles is a perfectly fine surface if rotated. F1d is really not adding anything either. Essentially the concept is that a sequence of locations on the retina, such as a,b,c are represented in this order on cortex and not for example as a,c,b. This is obviously hard to illustrate in a figure without retinal space. The same applies to F2, where I am not certain what it explains.

2. The project starts with the proposal that no assumptions are needed other than enforcing the “topological condition”, but unfortunately this is not exactly correct. In order to map across several visual areas, the accuracy of boundary labels is a major concern. An ill placed boundary can significantly impact the smoothing process, by introducing apparent breaches in topology that do not exist in the data. The section “Robustness with inaccurate interior boundaries” L664f tries to address this, but the section is entirely insufficient and has a number of flawed assumptions. I was also struggling with comprehension of this section.

3. Specifically, the authors imply that their method is robust against such errors, but it is entirely unclear what they mean by this. Will mislabelling a) simply be “smoothed over without larger consequences” or if b) the method somehow corrects these mistakes and instead finds the most correct boundary. I suspect it is a). If it is indeed b) the mechanism by which this is done needs to be clearly specified. In fact it is unclear how area labels are created in the very first place.

4. Assuming it is a) error will inevitably distort/smooth surrounding mapping, confounding attempts to quantify local CMF. It would be important to quantify these distortions, for example by systematically moving internal boundaries such as the V1/V2 border by medium to small amounts (15mm?) and investigate and quantify the effects.

5. Would it be possible to conceive a method that optimises the boundary placements? Maybe even use the Beltrami coefficient to find the optimal boundaries?

6. It remains unclear throughout the entire manuscript what Software (Language) the authors used for implementation. I had to go to the OSF site of the project to find matlab code. Searching the document matlab is cited as reference no 51, but [51] is actually not cited in the text.

7. The authors apply their method to datasets of 171 participants, but for what purpose? Only the results in Fig. 13 seem to rely on 171 datasets. It is a curious result but it breaks the consistency of the manuscript that otherwise seems to be a pure method paper.

8. The authors describe their technique as topological smoother. To me the two things topological correctness and smoothness are not necessary the same. The name seems to suggest the method smoothes until the topological condition is met. Does that mean datasets with large inconsistencies end up being more smoothed?

9. I was majorly struggling (failing) to follow large sections of the manuscript, in particular pages 11 to 16. I felt the authors are using a large number of symbols often inconstant and the purpose is not clear. I understand that it is hard to structure a logical problem that arose piece by piece when implementing a new method, but I feel the authors need to significantly improve this. I am unclear what the first subproblem exactly is, and what EXACTLY is smoothed. I understand what Laplacian smoothing is, but I am unclear what the mesh is that is smoothed.

10. I am unclear how it is the second subproblem is solved. Do I understand correct that nodes with negative coefficients are found, then their coefficient is changed by “shrinking” whatever that means. Shrinking sounds like an iterative process. L337 refers to Figure 2d which does not exist. I believe it should point to F2b red. Also the reference to fi seems wrong, fi denotes the node most to the right and about halfway up the positive axis. The red node does not seem to have a designation. How is that linked with the first subproblem?

11. L560F “angle distortion”. Angle distortion is typically called anisotropy, and yes, this has been “formally formalized”. It is an extremely common topic in the field of vision science with many authors discussing and analysing it, even Holmes and Lister (1916) were discussing it. The pattern anisotropy at eccentricities 0.5 and 1 degree displayed in Fig. 13 looks suspiciously like the one found by Schira et al. (2007, 2009 and 2010). The asymmetry for eccentricity ranges 1.5 and 2 degrees (yellow and purple curve) however has not been reported before. However, if the authors would like to present these results, they need to deliver a much more solid report on how it was obtained, they need to provide an analysis of the variance (i.e. error bars…) across subjects and make sure the results are not the consequence of a systematic error. Just glancing at these I would be very suspicious the anisotropy found at -pi/2 (ventral?) but not at pi/2 (dorsal?) could be a systematic error, as nothing like this has ever been reported. Extraordinary claims require extraordinary evidence.

12. L39 “The violation of topology condition is especially severe near the fovea, “ this remains a little fuzzy, not specified what ‘near the fovea’ means, like 1 degree eccentricity?

13. L62 “recent advances in Blood Oxygenation Level Dependent…” This statement feels a little bit out of time. These “recent” advances have allowed retinotopic mapping and estimation of magnification more than 25 years ago by several authors (Sereno et al, De Yoe et al, Engel et al.) “recent advances” is a phrase commonly used to highlight sub mm fMRI, and even in this field the phrase feels dated. This somewhat extends to the entire paragraph. In the last two decades many more visual areas than just V1, V2 and V3 have been discovered.

14. L120 suggests fi is “a mapping”. I find this confusing and suspect it’s wrong. I also do not understand how fi can be “moved out of any of the edges“ (L121). The way it looks to me, in F1c the edge fjfk (or the triangle fjfkfi) violates the surface. The authors really need to clarify their terminology here.

15. I suspect the authors may understand F1 as an illustration of a surface to surface mapping, where F1b) is the retina and F1c) is the cortical surface. Then for the node fi the topology of the mapping is violated. However, then F1a) is confusing and the entire Figure 1 does not consider the retina a part of the mapping. Then the triangle PiPjPk is in the retina/visual space and the tringle fifjfk in cortex?

16. L194 ”assuming that the voxel’s respective is (’;,)” what is a respective? I am not certain it is required for this manuscript to describe the equations of the pRF fitting procedure. A simple 4-5 line section with quoting the relevant manuscript describing the specific implementation used should suffice.

17. L291 “we now explain how to recovery…” should be “recover” I think.

18. L 595 “Quality assessment”. I am not certain what that paragraph is supposed to illustrate, other than apparently that for S1 the method changed the map only a little, while for S2 it changed a lot. What does that mean? What about the other 169 subjects? What it too much what is too little? Of course, quality control is an important step in any data processing/analysis and hence a worthwhile discussion topic. The authors should try and define some sort of criteria/logic. For example the authors could provide a quantitative measurement of the amount of smoothing and report how these values are distributed. This would allow to classify good vs bad datasets and provide a framework to judge this. One could argue this would go beyond the scope of the publication, but if included it must be of substance.

Reviewer #3: This is a method manuscript. It tackles a very specific problem in population receptive field (pRF) fitting algorithms. Basically, most of the pRF algorithms work at the voxel level, and per every voxel, the algorithm will provide the pRF center and size. By definition, those pRF centers need to be unique and retinotopically organized. Usually, these voxels are mapped to the surface, in a 2D surface mesh, so that every vertex corresponds to a time series. As these fMRI time series are very noisy, it can happen (in 20% of the cases according to the authors) that the pRF center values estimated from those time series break the topological condition. The authors propose a method to solve this issue, generating topology-preserving retinotopic maps.

pRF algorithms are very powerful but they are relatively modern (they were introduced in 2008). There are improvement opportunities, and it is important for the field that the authors identify a weakness of the method and propose a solution to it. The paper is clearly written, it is well organized and the figures are generally clear, which is a great pleasure as a reviewer.

Major comments

----------------------

- The authors published this method already in ISBI 2020 (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7406191/#!po=43.3333). I think that: (1) If the method has been significantly improved, the ISBI paper should be cited and here explain what those improvements consist of. (2) If there are no improvements on the methods but there are in the applications (CMF and angle distortion quantifications?, discussion, examples...), the ISBI paper should be cited for the methods part and add the new contribution here. I can see how much work went into this manuscript, but it is the job of the editors to decide if the contribution is big enough to be published in their journal, but as a reader I would expect the other paper to be acknowledged/explained in this one.

- Validation of the method:

--- Synthetic data: If I understand it correctly, you basically take a topology and map an invented set of pRF center parameters, with more or less noise. I don’t understand why is the topology cuadricular when the vertices in the cortex are organized in triangles. I guess they are “equipotential” (angle, eccen) lines and not the triangles of the surface?

------ How can you prove that those modifications are realistical? I think that it would make more sense simulating time series with known ground-truth values (using Lerma-Usabiaga 2020 Plos Comp. Biol., for example) and adding noise until 20% of the triangles are flipped. And then applying the correction methods. You would know the ground truth of the parameters to estimate and therefore quantify the correction of your method.

------ How can you know that typically 20% in an experimental dataset are flipped? Table S2?

--- Experimental data: I don’t think the R2 calculations are the adequate ones here… Are you basically saying that for the time series, the solver went to a local minima but that there was a better solution, and that your algorithm could find the better fit to the time series? Knowing how most of the algorithms work ( seed > grid search > fine search), this can happen, but maybe there are better ways to show this than the one you provide. And if they are not, the concept of but you think it is going on should be explained, I think.

------ How can we know that the fixed data is closer to the truth? We are modifying the estimated parameters according to one criterion: we don’t want to have more flipped triangles. Do we need to assume that the averaged improved R2s are saying that we are closer to the truth or we are only saying that now we don’t have flipped triangles only?

- Boundary definition: do we need to set the boundaries manually to initialize the method or this is only necessary for the tests in the paper? I think it would be very useful for the reader/final user that you would add a small “user guide”. I understand that the user continues with exactly the same pipeline he/she was using, and then this method is something that is applied on top of it. After reading the manuscript, the manual and automated, the old and the new steps are not perfectly clear.

Other comments

----------------------

- pRF Methods are not explained sufficiently: it seems from the text, as Dumoulin and Wandell 2008 is cited (Page 8, Line 168: Page 9, Line 192), that the authors use mrVista software to fit the time series, but in a second lecture, it seems that they use the same algorithm used in the HCP 7T publication, which I guess is Kendrick Kay’s AnalyzePRF?. I think that the software used should be explained specifically and maybe linking to publications/software repositories and versions, for clarity and reproducibility.

- There are no papers cited later than 2018, but there have been relevant publications in the field not acknowledged by the authors (for example Lage-Castellanos 2020, or the previously cited Lerma-Usabiaga 2020). In the same line, I missed discussing other approaches to pRF fitting, such as the paradigm free approaches (Merkel 2018, 2020).

- Page 2, Lines 44-45: that sentence is not correctly constructed.

- Ref 21 has an error it seems.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Anand Joshi

Reviewer #2: Yes: Mark M. Schira

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions, please see

19 May 2021

Submitted filename: ResponseLetter-Final.pdf

Click here for additional data file.

27 Jun 2021

Dear Dr Wang,

We are pleased to inform you that your manuscript 'Topology-Preserving Smoothing of Retinotopic Maps' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests. You may use this as an opportunity to make the minor change to figure 1, suggested by R2, if you wish.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Emma Claire Robinson

Associate Editor

PLOS Computational Biology

Wolfgang Einhäuser

Deputy Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #2: Up reading the review of “Topology-Preserving Smoothing of Retinotopic Maps” by Tu, Ta, Lu and Wang, I find that the authors made a genuine and successful effort to improve and clarify the manuscript. The new section on robustness to inaccurate internal and external boundaries is useful and interesting. I find the manuscript is now ready to publication. The authors have addressed most of my concerns, except perhaps for Figure 1. For Figure 1 I would suggest the authors replace the brain background in F1a by a retina or eye background. This would make it intuitively clear that it is the estimated projection between retina and cortex that violates the topology condition. One may add the cortex picture as background to F1b), but that may be to much.

Reviewer #4: The authors present a very intriguing case for preserving the topology of population receptive field maps after standard parameter estimation methods. This is a fairly novel concept, in that location parameters of population units are being optimized for the topological condition based on their neighboring units in 2-dimensional vertex space via Beltrami-optimization.

The authors did respond to the reviewers concerns in great detail. After examining the initial submission and the comments to the initial reviews I remain slightly concerned about potential impacts of the ultimate enforcement of the topological condition on the location maps of the pRFs and their relation to the ground truth: The authors are right in their initial assumption of the neurophysiological reality of the topological mapping within retinotopic areas. However, the adoption of this method to the HCP-dataset shows, that its ultimate enforcement on the population level comes with the cost of a reduced raw-data-fit. The Beltrami-optimization does only change the location parameters of the initial pRF-solution (a solution estimating location and size(sd)) in order to retain topology. The R^2 between the raw-timeseries and the Beltrami-optimized time-series has to be lower than the initial pRF-estimated timeseries (as shown by the authors). I see now the problem of inferring whether this reduction in explained variance is a consequence of the change in location parameters OR keeping the size parameter constant. Sort of an answer can be found in the ground truth analysis, in which the authors did not take RF size into account. Since R^2 does not change for those synthetic datasets, one can assume that ignoring RF size in the topological optimization process may distort the estimated maps relative to the ‚real‘ topography.

One important next step would be to test ground truth incorporating all parameters used for the adopted pRF model.

The current topology-retaining smoothing method is certainly of interest for the community but more extensive simulation studies are required to test its validity.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #2: Yes

Reviewer #4: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: Yes: Mark M. Schira

Reviewer #4: Yes: Christian Merkel

21 Jul 2021

PCOMPBIOL-D-20-02101R1

Topology-Preserving Smoothing of Retinotopic Maps

Dear Dr Wang,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Andrea Szabo

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol