Connectopic mapping techniques do not reflect functional gradients in the brain

Functional gradients, in which response properties change gradually across a brain region, have been proposed as a key organising principle of the brain. Recent studies using both resting-state and natural viewing paradigms have indicated that these gradients may be reconstructed from functional connectivity patterns via "connectopic mapping" analyses. However, local connectivity patterns may be confounded by spatial autocorrelations artificially introduced during data analysis, for instance by spatial smoothing or interpolation between coordinate spaces. Here, we investigate whether such confounds can produce illusory connectopic gradients. We generated datasets comprising random white noise in subjects' functional volume spaces, then optionally applied spatial smoothing and/or interpolated the data to a different volume or surface space. Both smoothing and interpolation induced spatial autocorrelations sufficient for connectopic mapping to produce both volume- and surface-based local gradients in numerous brain regions. Furthermore, these gradients appeared highly similar to those obtained from real natural viewing data, although gradients generated from real and random data were statistically different in certain scenarios. We also reconstructed global gradients across the whole-brain - while these appeared less susceptible to artificial spatial autocorrelations, the ability to reproduce previously reported gradients was closely linked to specific features of the analysis pipeline. These results indicate that previously reported gradients identified by connectopic mapping techniques may be confounded by artificial spatial autocorrelations introduced during the analysis, and in some cases may reproduce poorly across different analysis pipelines. These findings imply that connectopic gradients need to be interpreted with caution.


Introduction
Functional gradients, marked by a topographic map of preferred stimulus or response features across brain regions, are proposed as a key organising principle of the brain ( Bednar and Wilson, 2016 ). For instance, gradients have been described in the form of retinotopic maps in visual cortex ( Wandell and Winawer, 2011 ), tonotopic maps in auditory cortex ( Formisano et al., 2003 ), and somatotopic maps in somatosensory cortex ( Sanchez-Panchuelo et al., 2010 ). Recently, "connectopic mapping " techniques have been proposed for recovering functional gradients from patterns of neural connectivity ( Haak and Beckmann, 2020 ;Huntenburg et al., 2018 ). These approaches regard voxels or vertices within a brain region as being distributed along a high-dimensional manifold within the connectivity space, such that functionally similar locations are represented close to one another on the manifold surface. The functional gradients underlying a brain region can thus be revealed by applying non-linear manifold learning techniques to extract the principal dimensions of this manifold. brain region, rather than necessarily reflecting a topography embedded within the connectivity itself. We were able to confirm the accuracy of the connectopic maps in V1 by comparing them to retinotopic maps measured through traditional visual field mapping techniques. Nevertheless, the ground-truth functional organisation of many other brain regions is unknown or poorly understood, so there needs a to be a degree of confidence that gradients reconstructed from connectopic mapping accurately reflect the genuine underlying topography of such regions.
It has been suggested that gradients could represent an inevitability of the connectopic mapping analysis, rather than accurately reflecting the genuine functional topographies. Ciantar and colleagues noted a confound introduced during data pre-processing which can give rise to illusory local functional connectivity patterns, which could in turn potentially bias connectopic mapping analyses ( Ciantar et al., 2022 ). Specifically, interpolation of the timeseries data from the volume to the cortical surface is liable to induce local spatial autocorrelations amongst neighbouring vertices. Because connectopic gradients reflect the local similarity in connectivity patterns between voxels or vertices, it may be that such gradients are biased by these artificially induced spatial autocorrelations. Importantly, such autocorrelations could be introduced by other processing steps too, such as spatial smoothing or interpolation to other volume spaces. Nevertheless, it remains unclear to what extent this issue practically influences connectopic gradients.
Here, we provide a test of how local spatial autocorrelations affect connectopic mapping. We obtained naturalistic imaging data from a publicly available dataset. We then generated random white noise timeseries matched to the real data in each subject's functional volume. These random datasets are entirely biologically implausible, with data being both spatially and temporally uncorrelated, and so provide a strong test for whether gradients can be produced by artifacts introduced during the analysis pipeline alone. For both real and random datasets, we then optionally applied spatial smoothing and/or interpolated the data to a different volume space or to the cortical surface. We then identified both local and coarse-scale gradients via connectopic mapping. If smoothing and/or interpolation induces spatial autocorrelations sufficient to confound connectopic mapping, we would expect to reconstruct gradients even from random data and these should appear similar to equivalent gradients derived from real data.

Dataset
Movie-watching and visual field mapping MRI data were obtained for 15 subjects from the publicly available StudyForrest project ( Hanke et al., 2016( Hanke et al., , 2014 ; https://www.studyforrest.org ). Briefly, functional data were acquired on a 3T Philips Achieva MRI scanner via an EPI sequence (TR = 2 s, TE = 30 ms, voxel resolution = 3 mm isotropic). The dataset comprises approximately 2 h of MRI data of each subject watching the "Forrest Gump " movie.

Pre-processing
The same pre-processing steps were applied to movie-watching and visual field mapping datasets. Some light pre-processing was already applied by the StudyForrest project: motion correction using FSL's MCFLIRT tool ( Jenkinson et al., 2002 ) and aligning each volume to a subject-specific EPI template. We then applied further pre-processing using FSL's FEAT v6.0 ( Smith et al., 2004 ;Jenkinson et al., 2012 ; https://fsl.fmrib.ox.ac.uk/fsl ): slice-timing correction using Fourierspace time-series phase-shifting, non-brain removal ( Smith, 2002 ), grand-mean intensity normalisation of the entire 4D dataset by a single multiplicative factor, and high-pass temporal filtering (Gaussianweighted least-squares straight line fitting with = 50 s). Spatial smoothing was not yet applied at this stage. The timeseries were then standardised by converting to units of percentage signal change, and then finally both the mean ventricular and white-matter timeseries and motion parameters were regressed out. We computed registrations from each subject's EPI space to both standard volume and surface spaces. First, we computed a volumebased registration between each subject's EPI and T1 anatomical images via boundary based registration ( Greve and Fischl, 2009 ), and then further to the MNI152 standard brain via FSL's FNIRT tool ( Andersson et al., 2010 ). Next, cortical surfaces were reconstructed from T1-and T2-weighted anatomical images using Freesurfer v6.0 ( Dale et al., 1999 ;https://surfer.nmr.mgh.harvard.edu ). Functional data were co-registered to each subject's native cortical surface via boundary based registration, and then further to the fsaverage surface via a surface-based registration ( Fischl et al., 1999a( Fischl et al., , 1999b. To synthesise random data, we generated random timeseries comprising Gaussian white-noise matched in mean and variance to the real movie-watching data in the functional EPI volume of each subject. These timeseries will approximate the signal amplitudes of the real data, but will lack any coherent spatial or temporal correlation structure. Both real and random datasets were then optionally transformed to the MNI volume or cortical surface using trilinear interpolation. Data were transformed to the fsaverage5 surface for whole-brain analyses, and the fsaverage6 surface for all other analyses. We then applied spatial smoothing in all three spaces (functional volume, MNI volume, fsaverage surface) at FWHM = 6 mm (twice the voxel resolution), using volume-based smoothing for volume spaces and surface-based smoothing for surface spaces. We also retained the unsmoothed datasets for comparison. Volume datasets were restricted to cortical and subcortical grey-matter masks, derived from individualised Freesurfer segmentations ( Fischl et al., 2002 ) for the functional volumes (back-transformed from the T1 images), or from the Harvard-Oxford atlas for the MNI volume. Thus, each subject had a total of 12 movie-watching datasets: a real and complementary random dataset, with a smoothed and unsmoothed version of each, and with a version of each in the functional volume, MNI volume, and on the fsaverage surface.

Regions of interest
We measured connectopic gradients in a number of different regions of interest: (1) primary visual cortex (V1), (2) hippocampus, (3) striatum, and (4) Schaefer parcels. We first generated regions of interest (ROIs) for primary visual cortex in each subject. Visual field mapping data were either retained in the native functional volume or transformed to the fsaverage surface. Both spatially smoothed (FWHM = 6 mm) and unsmoothed versions were generated in each space. In each case, a travelling wave analysis ( Engel et al., 1997 ) was conducted using the 3dRetinoPhase command in AFNI ( Cox, 1996 ;Saad et al., 2001 ). This produced phase maps of the eccentricity and polar angle tunings. The phase maps from the spatially smoothed surface analysis were used to define individualised V1 ROIs by tracing along the phase reversals. These surface ROIs were additionally back-projected to each subject's functional volume. These ROIs were then applied to the analysis of both smoothed and unsmoothed data in each subject's functional volume and surface spaces. The eccentricity and polar angle phase maps also provided a ground-truth estimate of the retinotopic maps within V1 to compare the connectopic maps against. Next, for analysis of the MNI volume data, we defined ROIs for the hippocampus and striatum from the Harvard-Oxford subcortical atlas. Finally, for further analysis of the surface data, we generated surface ROIs from the 100-area Schaefer parcellation ( Schaefer et al., 2018 ). This parcellation scheme partitions the cortical surface along boundaries in local functional connectivity while simultaneously maximising the functional similarity within each parcel, and has been shown to predict functional boundaries as well as or better than comparable functional, anatomical, and multimodal parcellations ( Zhi et al., 2022 ).

Connectopic mapping
We tested two separate pipelines for performing the connectopic mapping analysis, based on subtly different methods described in the literature ( Haak et al., 2018 ;. Both pipelines aim to reconstruct functional gradients from patterns of connectivity, but use differing approaches to achieve this.

Pipeline 1
Our first connectopic mapping pipeline was based on the methods of Haak et al. (2018) . A schematic illustration of this pipeline for ROI analyses is shown in Fig. 1 . For analyses of V1 data, timeseries were concatenated over odd and even scan runs separately to allow crossvalidated parameter selection for the manifold learning. For all other analyses, timeseries were concatenated over all scan runs. Timeseries were split between grayordinates within the ROI and those outside the ROI. For surface-based analyses, this entailed splitting the timeseries between cortical vertices inside and outside of the ROI. For volumebased analyses, the split was between voxels inside and outside of the ROI, with non-ROI voxels selected from a grey-matter masque. Because the number of non-ROI grayordinates exceeds the number of timepoints, we reduced the dimensionality via a lossless PCA -retaining all available components (one fewer than the number of timepoints) so that 100% of the variance remained explained. This amounts to rotating the samples within the feature space and dropping the un-used dimensions, and aids computational tractability for later processing stages.
Next, connectivity fingerprints were derived by correlating timeseries between the ROI grayordinates and non-ROI principal components, and correlations were converted to units of Fisher's z . We then measured the pairwise Euclidean distances between samples (ROI grayordinates) in the connectivity space, yielding a symmetrical distance matrix. These distances were then used to construct an unweighted and undirected nearest neighbour graph. This is represented as a symmetrical sparse affinity matrix, where two samples are assigned a value of one if either sample falls within the k nearest neighbours of the other, and a value of zero otherwise. This requires a parameter ( k ) to be selected for the neighbourhood size. For analysis of the V1 ROI, the parameter was selected via cross-validation. Specifically, a Bayesian optimisation algorithm (implemented using the scikit-optimize Python module; https://scikit-optimize.github.io ) selected the parameter that maximised the correlation between retinotopic and connectopic maps in one of the data splits (odd or even runs), and the selected parameter was then applied to the analysis of the other data split. As such, the parameter selection remains statistically independent from the measure of the prediction accuracy. For the analysis of other regions, the ground-truth functional topography is unknown, so it wasn't possible to optimise a prediction accuracy. Instead, the parameter was set to the log of the number of samples (rounded to the nearest integer) as this has previously been suggested as a reasonable default choice ( von Luxburg, 2007 ). The affinity Fig. 1. Schematic illustration of the first connectopic mapping pipeline for ROI analyses, based on the methods of Haak et al. (2018) . Timeseries are divided between grayordinates within versus outside of the ROI. Non-ROI timeseries are reduced in dimensionality via lossless principal component analysis. The ROI and compressed non-ROI timeseries are cross-correlated to derive the connectivity space. Pairwise distances are measured between samples in the connectivity space, which are then used to derive a nearest neighbour network graph (represented by a sparse affinity matrix). Spectral embedding is then used to derive the connectopic maps. To perform connectopic mapping at the group level, the distance matrices are averaged over subjects, and then the remaining pipeline proceeds as normal. matrix was then decomposed via spectral embedding to yield the connectopic maps. In all cases, we retained the first two components. Spectral embedding was implemented using the scikit-learn Python module ( Pedregosa et al., 2011 ;https://scikit-learn.org ). The analysis can either be performed at the individual subject level, or at the group-level by averaging the distance matrices over subjects then deriving an affinity matrix and applying spectral embedding as normal. We conducted the V1 ROI analyses at the individual-level, the Schaefer parcel analyses at the group-level, and the MNI volume ROI analyses (hippocampus, striatum) at both individual-and group-levels.
We also applied a modified version of this pipeline to connectopic mapping of whole-brain gradients. Connectivity fingerprints were estimated for all vertices within each hemisphere (excluding the medial wall) by correlating the timeseries in those vertices against principal components derived from those same vertices (via a lossless PCA). From here, the remaining analysis pipeline proceeded as described above. This yielded whole-brain connectopic maps for each hemisphere. We performed the whole-brain analyses at the group-level.

Pipeline 2
We also employed an alternative connectopic mapping pipeline based on the methods of  , primarily for the analysis of whole-brain gradients ( Fig. 2 ). Timeseries were correlated between all vertices in each hemisphere (excluding the medial wall), yielding a symmetrical connectivity fingerprints matrix. Correlations were left in units of Pearson's r because the matrix diagonal comprises perfect correlations which would invalidate a Fisher's z transformation. Next, the connectivity fingerprints were made sparse by thresholding the values at the 90th percentile for each sample (row) in the matrix. Where necessary, values where further thresholded at zero to ensure that all values were non-negative. Next, the thresholded connectivity matrix was converted into an affinity matrix using a similarity kernel: we employed kernels based on cosine, normalised angle, and Pearson's correlation similarity. Although the connectivity matrix is square in the whole-brain case, it is still necessary to derive a separate symmetric affinity matrix because the connectivity matrix is rendered asymmetric by the thresholding. This affinity matrix is dense, such that the similarity is measured from every vertex to every other vertex, and hence represents a weighted and fully-connected network graph. Spectral embedding requires the affinity matrix to be non-negative: for cosine and angular similarity this is already guaranteed by the non-negativity of the thresholded connectivity fingerprints, but for correlation similarity the affinity matrix was further thresholded at zero. Finally, the affinity matrix was decomposed via spectral embedding, and the first two components were retained. This can be performed at the individual-subject level, or at the group-level by averaging the connectivity fingerprints over subjects before thresholding and then proceeding with the remaining pipeline as normal. To obtain an unbiased average, the connectivity fingerprints were converted to units of Fisher's z before averaging, and then converted back to Pearson's r after averaging. The principal difference between the analysis pipelines is that the first pipeline uses dense connectivity fingerprints and a sparse affinity matrix, while the second pipeline uses sparse connectivity fingerprints and a dense affinity matrix.
A modified version of the second pipeline was also applied to surfaceand volume-based ROI analyses of V1. Here, connectivity fingerprints were derived by correlating timeseries between grayordinates within the V1 ROI and all other cortical grayordinates outside of the ROI. Beyond this point, the remaining pipeline stages proceeded as described above. We conducted the V1 ROI analyses at the individual-level, and the whole-brain analyses at the group-level.

Statistical analyses
For the V1 ROI analyses, we estimated retinotopic and connectopic maps in each subject. Prediction accuracies were measured by the absolute correlation between the eccentricity and polar angle retinotopic and connectopic maps; the absolute value was taken to account for the sign ambiguity in the connectopic maps. Correlations were converted to units of Fisher's z and averaged within each subject over hemispheres and the two cross-validation splits. These correlations were then entered into a four-way repeated-measures ANOVA with factors for the co-ordinate space (volume/surface), randomness (real/random), smoothness (smoothed/unsmoothed) and map (eccentricity / polar angle). Next, we measured the absolute correlations between the real and random connectopic maps. Values were again converted to Fisher's z and averaged within subjects. These correlations were then submitted to a three-way repeated-measures ANOVA with factors for the co-ordinate space, smoothness, and map.
Connectopic maps in the MNI volume ROIs (hippocampus and striatum) were estimated at both the group-and individual-level. In each case, we measured the absolute correlations between the real and ran-

Fig. 2.
Schematic illustration of the second connectopic mapping pipeline for whole-brain analyses, based on the methods of  . Whole-brain timeseries are cross-correlated to derive the connectivity space. The connectivity space is made sparse by thresholding at the 90th percentile for each sample (row) in the matrix. This is then converted to a dense (fully-connected) affinity matrix via a similarity kernel. Spectral embedding is then used to derive the connectopic maps. To perform connectopic mapping at the group level, the initial dense connectivity fingerprints are averaged over subjects, and then the remaining pipeline proceeds as normal.
dom connectopic maps, converted the values to Fisher's z , and averaged over hemispheres. Individual correlations were entered into a two-way repeated-measures ANOVA with factors for the smoothness and component (primary and secondary gradients). To assess the consistency of the maps between individuals, we also calculated the pairwise correlations of individual-level gradients between subjects, as well as the correlation of each subject's gradient to the group-level gradient. We additionally estimated noise ceilings from these values by averaging the between-subject correlations over pairwise combinations for each real connectopic map -this estimates the average similarity of each real gradient to itself between subjects. Correlations were converted to Fisher's z before averaging to obtain an unbiased mean, which was then converted back to Pearson's r for visualisation. We compared the Fisher transformed within-subject correlations (between real and random connectopic maps) against these noise ceilings via a series of one-sample t-tests, including Bayes factors calculated with the BayesFactor R package ( https://cran.r-project.org/package = BayesFactor ).
We also estimated group-level connectopic maps within each of the 100 Schaefer parcel surface ROIs. Absolute correlations were measured between real and random connectopic maps for each parcel, converted to Fisher's z , and entered into a two-way repeated-measures ANOVA with factors for the smoothness and component.
Finally, we estimated group-level connectopic maps over the wholebrain within each hemisphere. At the whole-brain level, we observed more variable degrees of similarity in real and random connectopic maps within and between components. We therefore measured the full 2 × 2 pairwise absolute correlations between real and random maps across both components within each hemisphere. We assessed the significance of these correlations using spin permutation tests (1000 permutations per hemisphere) implemented using the BrainSpace toolbox ( Vos de Wael et al., 2020 ; https://brainspace.readthedocs.io ). On each permutation, the real gradients were randomly rotated over the cortical sphere, and the absolute correlation was measured between the rotated real gradient and original random gradient (excluding invalid vertices within the medial wall). Repeating this process over all permutations produces an empirical null distribution of correlation values. Statistical significance was determined by identifying the proportion of correlations in the null distribution that were greater than the genuine correlation value.
In all ANOVAs, we report effect sizes in units of both partial and generalised eta-squared ( Bakeman, 2005 ;Olejnik and Algina, 2003 ). All statistical tests employed an alpha criterion of 0.05 for determining significance.

Connectivity searchlight
To assess how the pre-processing pipeline affects spatial autocorrelations, we used a searchlight analysis ( Ciantar et al., 2022 ;Kriegeskorte et al., 2006 ) to measure local functional connectivity. We employed spherical searchlights for volume spaces (5 mm / 2 voxel radius for functional volumes, 6 mm / 2 voxel radius for the MNI volume) and a disc searchlight (6 mm radius) for the surface space. Volume searchlights were restricted to a grey-matter mask. For each searchlight, we calculated the pairwise correlations in timeseries between all grayordinates within the searchlight. These correlations were converted to units of Fisher's z and averaged over grayordinates, and then this average was converted back to Pearson's r and assigned to the central grayordinate of the searchlight. This resulted in a map of the local functional connectivity over the whole brain. We performed searchlight analyses in each individual subject for all 12 versions of the dataset: real and random data, smoothed and unsmoothed, and in the functional volume, MNI volume, and surface spaces. Volume searchlights were implemented with the PyMVPA toolbox ( Hanke et al., 2009 ; http://www.pymvpa.org/ ), and surface searchlights were implemented using custom code leveraging the surfdist toolbox ( Margulies et al., 2016 ; https://github.com/NeuroanatomyAndConnectivity/surfdist ).

ROI connectopic mapping
We employed our first connectopic mapping pipeline to analyses of local regions of interest. We first examined gradients estimated by retinotopic and connectopic mapping in primary visual cortex ( Fig. 3 ). We employed both surface-and volume-based analyses. For surface analyses, the timeseries were interpolated to the cortical surface and retinotopic and connectopic mapping were performed on the surface. For volume analyses, retinotopic and connectopic mapping were performed within the native functional volume, and the resulting maps were interpolated to the surface for visualisation. We observed highly similar connectopic maps between the real and random data, indicating that spatial autocorrelations induced via spatial smoothing and/or interpolation to the surface were sufficient to produce illusory gradients. Retinotopic and connectopic maps are shown for a representative subject in Fig. 3 a. As previously reported ( Haak et al., 2018 ;Watson and Andrews, 2022 ), connectopic mapping of the real movie-watching data identified a posterior-anterior gradient corresponding to the eccentricity map, and a superior-inferior gradient corresponding to the polar angle map. This was true for both surface and volume analyses and using both smoothed and unsmoothed data. However, similar gradients were also obtained from connectopic mapping of the random data. The only exception to this was the case of the unsmoothed random data in the functional volume, where the connectopic mapping is operating directly on the raw random data with no further manipulation.
We next tested volume-based connectopic gradients in the hippocampus and striatum. We again observed very similar connectopic maps between real and random data, indicating that spatial smoothing and  interpolation to the MNI volume were sufficient to produce illusory gradients. Furthermore, these gradients were highly consistent between subjects. Group-level connectopic maps of the hippocampus are shown in Fig. 4 a, and individual-level gradients are shown for a representative subject in Supplementary Fig. 2a. Connectopic mapping of the real data identified a primary gradient running posterior-anterior along the long axis of the hippocampus, replicating previously reported gradients ( Borne et al., 2023 ;Prze ź dzik et al., 2019 ;Vos de Wael et al., 2018 ) in this region, and a secondary gradient that sub-divided the primary gradient. However, we also observed similar gradients using the random data. Unlike the V1 analyses in subjects' functional volumes, these gradients were now also apparent for unsmoothed random data -this indicates that interpolation to the MNI space was sufficient to induce the necessary spatial autocorrelations even without further spatial smoothing. Correlations between the real and random connectopic maps appeared high at both the group and individual levels ( Fig. 4 b). Entering the individual correlations into a twoway repeated-measures ANOVA revealed no significant main effect of smoothness (F(1,14) = 0.02, p = .883, 2 < 0.01, 2 < 0.01), but did indicate a significant main effect of the component (F(1,14) = 535.13, p < .001, 2 = 0.97, 2 = 0.70) due to higher correlations for the primary gradient, and a significant smoothness by component interaction (F(1,14) = 23.68, p < .001, 2 = 0.63, 2 = 0.21) due to a bigger difference between gradients for smoothed data. We additionally compared individual-level gradients between subjects ( Supplementary Fig.  3a) -these revealed very high spatial correlations both between subjects and to the group-level for both real and random gradients. We estimated noise ceilings by averaging over the between-subject correlations for the real gradients. The within-subject correlations between real and random gradients ( Fig. 4 b) approached the noise ceilings. The correlations fell significantly below the noise ceiling for the unsmoothed primary (t(14) = 4.76, p = .001, Hedges' g = 1.16, BF 10 = 109.00) and smoothed secondary gradients (t(14) = 4.11, p = .003, Hedges' g = 1.00, BF 10 = 36.54). However, the correlations for the unsmoothed secondary (t(14) = 0.95, p = 717, Hedges' g = 0.23, BF 10 = 0.39) and smoothed primary gradients (t(14) = 0.93, p = .717, Hedges' g = 0.23, BF 10 = 0.38) did not differ significantly from the noise ceilings, and Bayes factors indicated slight support for the null hypothesis. Thus, very similar hippocampal gradients were obtained from real and random data, and these were highly consistent over subjects.
Group-level connectopic maps of the striatum ( Fig. 4 c) indicated a primary gradient that curved around the striatum and a secondary gradient running superior-inferior, replicating previously reported gradients ( Marquand et al., 2017 ). However, these gradients were again observed for both the real and random data with and without spatial smoothing. Similar gradients were also found at the individual-level (Supplementary Fig. 2b). Correlations between real and random gradients ( Fig. 4 d) again appeared high at both the group and individual-level. A two-way repeated-measures ANOVA on the individual correlations indicated significant main effects of smoothness (F(1,14) = 53.15, p < .001, 2 = 0.79, 2 = 0.52) due to higher correlations for smoothed data, and of the component (F(1,14) = 400.06, p < .001, 2 = 0.97, 2 = 0.74) due to higher correlations for the primary gradient, while the smoothness by component interaction was not significant (F(1,14) = 3.58, p = .079, 2 = 0.20, 2 = 0.03). We again observed very high spatial correlations between subjects for both real and random gradients ( Supplementary  Fig. 3b). We also found that the within-subject correlations between real and random gradients approached the noise ceilings estimated from the between-subject correlations ( Fig. 4 b). Within-subject correlations fell significantly below the noise ceilings for the unsmoothed primary (t(14) = 9.09, p < .001, Hedges' g = 2.22, BF 10 = 5.20 × 10 4 ), unsmoothed secondary (t(14) = 6.31, p < .001, Hedges' g = 1.54, BF 10 = 1231.97), and smoothed primary gradients (t(14) = 10.55, p < .001, Hedges' g = 2.58, BF 10 = 2.76 × 10 5 ). However, correlations did not differ significantly from the noise ceiling for the smoothed secondary gradient, and the Bayes factor indicated slight support for the null hypothesis (t(14) = 0.95, p = .356, Hedges' g = 0.23, B 10 = 0.39). Therefore, again, striatal gradients appeared very similar between real and random data and highly consistent over subjects. Thus, both spatial smoothing and interpolation to a different volume space (the MNI brain) induced sufficient spatial autocorrelations to obtain local connectopic gradients in the hippocampus and striatum, and these replicated gradients previously reported in the literature. We next examined local surface-based gradients across the whole brain. We estimated connectopic maps for each region of the 100-area Schaefer surface parcellation ( Schaefer et al., 2018 ). Highly similar local connectopic maps were observed between real and random data throughout the brain, again indicating that spatial smoothing and/or interpolation to the surface can induce illusory gradients. Group-level connectopic maps are illustrated for the left hemisphere in Fig. 5 a, and for both hemispheres in Supplementary Fig. 4. Gradients typically ran smoothly and continuously over the extent of each region, with primary and secondary gradients usually appearing orthogonal to one another. These gradients were again obtained using both real and random smoothed and unsmoothed data. Correlations between real and random gradients appeared high in all parcels ( Fig. 5 b). Entering the correlations over parcels into a two-way repeated-measures ANOVA revealed a significant main effect of smoothing (F(1,99) = 34.77, p < .001, 2 = 0.26, 2 = 0.09) due to higher correlations for smoothed data, and of the component (F(1,99) = 101.65, p < .001, 2 = 0.51, 2 = 0.13) due to higher correlations for the primary gradient, while the smoothness by component interaction was not significant (F(1,99) = 0.62, p = .434, 2 < 0.01, 2 < 0.01). Thus, spatial smoothing and/or interpolation to the surface induced sufficient spatial autocorrelations to generate local connectopic gradients throughout the brain. Similar results were obtained using an alternative multi-modal brain parcellation ( Glasser et al., 2016 ; Supplementary Fig. 5).

Whole brain connectopic mapping
All of the gradients we have considered so far have been on a local scale within regions of interest. However, connectopic mapping has also been used to describe coarse-scale gradients spanning the whole brain . Such studies have typically reported a primary gradient that runs between unimodal and transmodal regions, and a secondary gradient that runs between primary visual and somatosensory regions. To test how such coarse-scale gradients are affected by spatial autocorrelations, we ran both our connectopic mapping pipelines for the whole-brain within each hemisphere. Similar to the region of interest analyses, the first connectopic pipeline (using a sparse affinity matrix) identified gradients from both real and random data. However, none of these gradients appeared biologically plausible and they did not replicate previously reported whole-brain gradients. The second pipeline (using a dense affinity matrix) successfully replicated previous gradients, and appeared less susceptible to spatial autocorrelations, but was sensitive to the choice of similarity kernel.
Group-level connectopic maps derived from the first connectopic mapping pipeline are shown for the left hemisphere in Fig. 6 a (data for both hemispheres are shown in Supplementary Fig. 6). Using this pipeline, we failed to replicate previously reported gradients even with the real data, and instead identified biologically implausible gradients that simply ran continuously across the brain along one axis (anteriorposterior and superior-inferior for the primary and secondary gradients respectively). Similar gradients were obtained between real and random spatially smoothed data. Gradients were still obtained for unsmoothed random data, but appeared less similar to their real data counterparts. We measured the pairwise absolute correlations between real and random gradients ( Fig. 6 b), and determined statistical significance via spin permutations ( Vos de Wael et al., 2020 ). With spatial smoothing, real and random gradients correlated strongly and significantly within each component. Without spatial smoothing, moderate correlations were obtained both within and between components, though spin permutations indicated none were significant. Fig. 6. Group-level connectopic mapping results for whole-brain surface analyses. (a,c) Primary and secondary gradients in the left hemisphere for the first and second pipelines respectively. (b,d) Absolute correlations between real and random gradients for the first and second pipelines respectively. Pairwise correlations were measured between both gradients, and significance was determined via spin permutation tests. Results for the second pipeline are illustrated for the cosine similarity kernel.
Our first connectopic mapping pipeline is based on the methods described by Haak et al. (2018) . It is possible that the failure to replicate previously reported whole-brain gradients reflects subtle methodological differences between this pipeline and the one employed by  . One key difference lies in the treatment of the connectivity fingerprints and generation of the affinity matrix ( Figs. 1 &  2 ). Haak et al. measure dense connectivity fingerprints (from all grayordinates to all principal components), then construct a sparse affinity matrix (measuring affinity only between neighbouring samples in the connectivity space). By contrast, Margulies et al. make the connectivity fingerprints sparse by thresholding them, then construct a dense / fully-connected affinity matrix (measuring the affinity between every sample and every other sample). To test the influence of these methodological differences on the gradients, we re-ran our whole-brain analyses using our second pipeline based on the methods of Margulies et al. We first constructed the affinity matrix using a cosine kernel, as per  . Group-level connectopic maps are shown for the left-hemisphere in Fig. 6 c and for both hemispheres in Supplementary Figure 7. Using either smoothed or unsmoothed real data, we were now able to replicate the unimodal-transmodal primary gradient and the visual-somatosensory secondary gradient. Using smoothed random data produced gradient-like patterns, but these did not correspond well to the gradients obtained with the real data and did not appear biologically plausible. Using unsmoothed random data did not produce convincing gradients -the maps instead converged on local brain regions and appeared largely uniform elsewhere. Correlations between real and random maps were generally low, and spin permutations did not indicate any of them to be significant. To ensure that the key factor in producing these gradients was the sparseness of the affinity matrix, we also ran a variant of the second pipeline in which the connectivity fingerprints were again made sparse by thresholding, but these were then used to generate a sparse affinity matrix as per the first pipeline. The resulting gradients again ran continuously along one axis and appeared similar to those obtained from the first pipeline ( Supplementary Fig. 8). Thus, the ability to reproduce previously reported whole-brain gradients critically depends on using a dense rather than sparse affinity matrix.
For the second pipeline, various similarity metrics can be used to construct the affinity matrix. In the above analyses, we used a cosine Fig. 7. Effect of varying the similarity kernel on the second connectopic mapping pipeline for group-level whole-brain surface analyses of the spatially smoothed real data. Plots illustrate primary and secondary gradients. Affinity matrices were derived using cosine, normalised angle, or Pearson's correlation similarity kernels. kernel as per  , however other metrics are available. To test how the choice of kernel affects the resulting gradients, we re-ran our analyses of the real and spatially smoothed data using cosine, normalised angle, and Pearson's correlation similarity metrics. Group-level gradients are illustrated in Fig. 7 . While gradients do appear similar in all cases, they also differ in some key features that alter the interpretation of the gradients. For instance, using a normalised angle metric, the primary gradient represented a divide between primary visual cortex and the rest of the brain, rather than a unimodal-transmodal distinction. The secondary gradient under normalised angle similarity more strongly emphasised somatosensory versus superior and anterior temporal regions, rather than visual versus somatosensory. Gradients from the correlation kernel appeared to fall between the cosine and normalised angle solutions. Thus, although we can reproduce previously reported whole-brain gradients from this dataset, they are sensitive to the choice of analysis pipeline. Specifically, they critically depend on using a dense rather than sparse affinity matrix, and are sensitive to the choice of kernel used to construct the affinity matrix.
Because the second pipeline appears less susceptible to spatial correlations at the whole-brain level, we also tested how this pipeline (using a cosine similarity kernel) performed at a local level by re-running our analyses of primary visual cortex (Supplementary Fig. 9). In general, the connectopic maps provided a poorer reconstruction of the retinotopic maps than those obtained with the first pipeline. At the same time, while the degree of similarity between real and random connectopic maps decreased relative to the first pipeline, the correlations were still moderate and appeared comparable to the retinotopic prediction accuracies. Thus, while the second pipeline appears to reduce the influence of spatial autocorrelations on local connectopic maps, it does not eliminate the problem, and meanwhile produces poorer reconstructions of the ground-truth retinotopic maps compared to the first pipeline.

The effect of smoothing and interpolation on local functional connectivity
Finally, to test how the analysis pipeline affects spatial autocorrelations and local functional connectivity, we performed a series of connectivity searchlight analyses ( Fig. 8 ). These indicated that spatial autocorrelations could be induced even in random data by applying spa-tial smoothing and/or interpolating between volume or surface spaces. Searchlight connectivity maps for a representative subject are shown in Fig. 8 a-c, and distributions of correlations over all subjects are illustrated in Fig. 8 d. As expected, local functional connectivity was near zero for the unsmoothed random data in the functional volume, confirming that the random data as generated are spatially uncorrelated. However, both spatial smoothing and interpolation to another co-ordinate space (either volume-or surface-based) increased spatial correlations in both real and random data. It is likely that these artificially induced spatial autocorrelations form the basis of the connectopic maps identified in the real and random data. As previously reported ( Ciantar et al., 2022 ), the distribution of local functional connectivity over the surface covaried with the surface curvature -correlations increased in regions with high curvature ( Fig. 8 e). We also correlated the average surface searchlight maps between the different datasets ( Fig. 8 f). Moderate-tostrong positive correlations were found between most real and random smoothed and unsmoothed datasets (with the exception of randomsmoothed against real-unsmoothed) indicating some commonality in the pattern of local connectivity between datasets. Thus, both spatial smoothing and interpolation to either the surface or a different volume space induce artificial spatial autocorrelations sufficient to produce illusory gradients with connectopic mapping.

Discussion
In this study we investigated how local spatial autocorrelations introduced during data analysis influence gradients produced by connectopic mapping. We synthesised spatially uncorrelated random data in each subject's native functional volume, then explored the effect of applying spatial smoothing and interpolation to other co-ordinate spaces. We found that artificial spatial autocorrelations could be induced by applying spatial smoothing and/or interpolating to the cortical surface or MNI volume space. These spatial autocorrelations were sufficient to derive local gradients from connectopic mapping of both volume-and surface-based brain regions. Connectopic maps derived from the random data appeared highly similar to those produced from the real data and replicated gradients previously reported in the literature. Coarserscale gradients over the whole-brain appeared less susceptible to artificial spatial autocorrelations but were sensitive to specific features of the analysis pipeline. We first investigated local gradients within both volume-and surface-based brain regions using our first connectopic mapping pipeline. Connectopic mapping of real movie-watching data reproduced previously reported volume-and surface-based gradients in primary visual cortex ( Haak et al., 2018 ;Watson and Andrews, 2022 ), the hippocampus ( Borne et al., 2023 ;Prze ź dzik et al., 2019 ;Vos de Wael et al., 2018 ), and the striatum ( Marquand et al., 2017 ). We also reported novel surface-based local gradients throughout the brain within regions of the Schaefer and Glasser parcellations ( Glasser et al., 2016 ;Schaefer et al., 2018 ). However, we were also able to obtain highly similar gradients from random data following either spatial smoothing or interpolation to other volume or surface spaces. This suggests the manifold structure within the connectivity space is biased by spatial autocorrelations (naturally present in the data and/or introduced during data processing), such that neighbouring samples within the connectivity space tend to also be anatomically proximal. Consequently, the connectopic mapping tends to yield gradients that flow continuously over the extent of the region, even when derived from random data. The only combination that did not yield gradients was applying connectopic mapping directly to the raw unsmoothed random data in the native functional volume without any further manipulation. However, the complete lack of any spatial correlations in the raw random data is biologically implausiblesome degree of spatial correlations will always be present in real data, even if data are retained in the native volume and no artificial spatial smoothing is applied. While a functional gradient does necessitate some degree of spatial smoothness in the data, the reverse is not true and spatially smooth data could still reflect non-topographic organisations. Thus, with real data it will be challenging to disentangle the effects of genuine functional topographies versus non-topographic spatial correlations on local connectopic gradients.
We note that many analyses in this work rely on similarity in the visual appearance of gradients, and it is possible that connectopic maps in real data may be distinguishable from smoothed/interpolated random data. Indeed, in primary visual cortex we observed significantly higher prediction accuracies of the retinotopic maps for real compared to random connectopic maps, providing support for real connectopic gradients being distinct from random gradients induced by smoothing and/or interpolation. This may indicate contributions of genuine functional signals from the real data, or could reflect the higher spatial autocorrelations in the real compared to random data -the connectivity searchlight analysis revealed that spatial correlations were generally higher in the real than the random data even after smoothing or interpolation (see Fig. 8 ). Furthermore, the correlations between the real and random connectopic maps were typically comparable to or higher than their correlations with the retinotopic maps, so the greatest similarity was observed between the real and random gradients themselves. As such, further work is needed to establish to what degree smoothing/interpolation explain connectopic mapping gradient findings. For instance, future studies could titrate the degree of spatial autocorrelations to more fully determine the extent of their influence on connectopic gradients. Nevertheless, our results indicate that standard data preprocessing steps, such as spatial smoothing and interpolation between co-ordinate spaces, are sufficient to induce illusory connectopic gradients.
The organisation of the gradients also appeared highly consistent between regions, with the primary gradient typically running continuously along a single axis from one side of the region to the other, and the secondary gradient usually either subdividing or running orthogonal to the primary gradient. This suggests that these local gradients may be an inevitability of the connectopic mapping technique, rather than necessarily reflecting genuine underlying functional topographies. The specific direction of the gradients could be influenced by interactions with the surface curvature ( Ciantar et al., 2022 ), or could be biased by the shape of the brain region -for instance being orientated relative to the long axis of the region. Such anatomical constraints might influence both connectopic mapping and functional brain organisation so that connectopic and functional gradients may coincide, as was found in primary visual cortex. However, this is not always guaranteed to be the case, and indeed connectopic mapping failed to reproduce the phase reversals in the retinotopic polar angle map (which delineate the boundaries between retinotopic regions) of an expanded early visual region including V1, V2, and V3 ( Supplementary Fig. 1). Even where functional and connectopic gradients do coincide, there remains the issue that the connectopic maps are likely to be only partially driven by the underlying functional topography. In a region where the functional topography is unknown or poorly understood, it would therefore be difficult to ascertain how accurately the connectopic mapping is able to reproduce such topographies. Previously, we noted that the source of connectivity signals had minimal impact on the gradients reconstructed from connectopic mapping ( Watson and Andrews, 2022 ), and interpreted this as indicating that connectopic mapping largely reflects the functional topography within a brain region, rather than necessarily a topography embedded within the connectivity itself. Our current results partially accord with this principle, but further imply that the topography within a region need not even have a functional origin. Instead, illusory local gradients may be driven by an artificial topography induced by spatial autocorrelations introduced during the analysis pipeline.
One key issue is how connectopic mapping might be expected to perform under the null hypothesis. Spatial autocorrelations within the data make it unlikely that a random white noise pattern would be produced even when derived from spatially correlated random data. However, it is not guaranteed that the technique would yield similar gradients from both the real and random data. The spatial extent of the autocorrelations will be less than the size of the brain regions, so spatially correlated random data might have produced a smooth but still disorganised pattern rather than a continuous gradient. Alternatively, gradients might have been produced but following a different axis or pattern between the real and random data. In practice, we did obtain continuous gradients in each region, presumably as the connectopic mapping chained the local spatial correlations together over the extent of the region, and these appeared highly similar between the real and random data. For the hippocampus and striatum regions, we observed that individual-level gradients were highly consistent both between subjects and between each subject and the group-level gradients ( Supplementary Fig. 3). We estimated noise ceilings by averaging the between-subject correlations for the real gradients, indicating how well each real gradient correlated with itself between subjects, and observed that the within-subject correlations between real and random gradients approached these noise ceilings. This not only indicates that the real and random hippocampal and striatal gradients were highly consistent over subjects, but also that the similarity between the real and random gradients was close to the ceiling level. No genuine gradients are present within the random data itself, so any gradients produced by the connectopic mapping must be illusory. The high similarity observed between real and random connectopic gradients therefore indicates that the spatial autocorrelations driving the random gradients likely also strongly influence the gradients obtained from real data.
We also tested how spatial autocorrelations influence coarser-scale surface-based gradients over the whole-brain. Using our first connectopic mapping pipeline, based on the methods of Haak et al. (2018) , we obtained whole-brain gradients that appeared highly consistent between spatially smoothed real and random data. Without spatial smoothing, gradients were again obtained for random data, though appeared less similar to their real data counterparts. However, in all cases, gradients simply ran smoothly and continuously along a single axis along the cortical surface. Such gradients appear biologically implausible, and failed to replicate the unimodal-transmodal and visualsomatosensory gradients previously reported by   Xia et al., 2022 ;Zhang et al., 2022 ), structural connectivity ( Park et al., 2021a( Park et al., , 2021b, and brain microstructure . However, all previous replications used a connectopic mapping pipeline based on the methods of  , which differs subtly from Haak et al.'s ( 2018 ) methods. Specifically, Haak et al. estimate a sparse affinity matrix in which the affinity is only measured between samples within local neighbourhoods of the connectivity space, while Margulies et al. estimate a dense affinity matrix in which the affinity is measured between all samples. We were able to reproduce the unimodal-transmodal and visual-somatosensory gradients when repeating our analyses with a second connectopic mapping pipeline based on Margulies et al.'s methods using a dense affinity matrix. Furthermore, these gradients were not produced by the random data, indicating that this pipeline is less susceptible to the spatial autocorrelations. Employing a sparse affinity matrix is intended to aid in identifying nonlinearities within the geometric structure of the manifold by restricting the network graph to connections between neighbouring samples lying along the manifold surface, and ignoring connections that may jump across gaps or pass through the manifold surface. However, if spatial autocorrelations in the data make it so that neighbouring samples in the connectivity space also tend to be neighbouring in the anatomical space, then a sparse affinity matrix may be more prone to yielding connectopic maps in which gradients simply run continuously over the extent of the brain region. By contrast, a dense affinity matrix, representing a fullyconnected and weighted network graph, provides a coarser-scale view of the manifold that may be less biased by local spatial correlations. However, such a representation may also make it more challenging to identify non-linearities in the manifold and will hence promote more linear solutions. Indeed, the second pipeline did perform more poorly at reconstructing retinotopic maps in primary visual cortex compared to the first pipeline ( Supplementary Fig. 9).
We further tested how the choice of similarity kernel affects the whole-brain gradients produced by the second connectopic mapping pipeline. While  constructed their affinity matrix using cosine similarity, other metrics such as normalised angle or correlation similarity are also available ( Vos de Wael et al., 2020 ). While the gradients appeared similar in all cases, they differed in a number of important ways that could alter the interpretation of the gradients. For instance, using the normalised angle kernel, the primary gradient shifted to run between visual cortices and the rest of the brain, rather than providing a unimodal-transmodal distinction. Similarly, the secondary gradient under the normalised angle kernel ran between somatosensory versus frontal and superior and anterior temporal cortices, rather than between visual and somatosensory regions. A correlation similarity kernel produced gradients appearing between the cosine and normalised angle kernel solutions. Thus, although the unimodal-transmodal and visual-somatosensory gradients appear reproducible across studies and datasets, they may be less robust to changes in the connectopic mapping pipeline, and indeed this has been described previously ( Bajada et al., 2020 ). Specifically, reproducing these specific gradients depends on using a dense affinity matrix representing a fully-connected network graph, and the gradients are sensitive to the choice of similarity kernel. Functional topographies are unknown or poorly understood in many brain regions, so it is important to be confident that the gradients produced by connectopic mapping accurately reflect the genuine underlying topographies. Given the issues discussed above, how confident can we be in connectopic gradients? There are multiple features of previously identified connectopic maps that would be difficult to explain away purely with spatial autocorrelations or specific features of the analysis pipeline. Previously, we reported better reconstructions of retinotopic gradients from connectivity measured during movie watching than at rest ( Watson and Andrews, 2022 ), demonstrating connectopic maps are influenced by the stimulus. Furthermore, individual variability in local gradients in both the hippocampus and striatum has been linked to human behaviour ( Borne et al., 2023 ;Marquand et al., 2017 ;Prze ź dzik et al., 2019 ). Meanwhile, the topology of whole-brain gradients corresponds well with the organisation of large-scale restingstate networks . Whole-brain gradients have also been shown to modulate with task demands ( Gao et al., 2022 ;Samara et al., 2023 ;Zhang et al., 2022 ), correlate with brain microstructure and structural connectivity ( Larivière et al., 2020 ), change with ageing ( Bethlehem et al., 2020 ;Park et al., 2021a ), and differ between neurotypical and neurodiverse populations Lee et al., 2023 ;Park et al., 2021b ;Xia et al., 2022 ). In some cases, these differences may interact with spatial autocorrelations -for instance different stimuli or tasks may evoke responses with differing signal to noise ratios which in turn alter the degree of spatial autocorrelation. Nevertheless, it seems likely that connectopic maps do at least partially reflect genuine underlying functional topographies and connectivity. However, it is also clear that connectopic maps can be strongly influenced by artificially induced spatial autocorrelations and specific features of the analysis pipeline. It may therefore be challenging to disentangle the influence of functional versus artificial signals on connectopic maps, thus the extent to which these gradients truly capture underlying functional topographies remains unclear. Connectopic mapping techniques may be able to reveal genuine functional gradients, but this may require further methods to account for spatial correlations within the data or at least a consideration of how such correlations may influence the resulting gradients. Comparing connectopic maps to behaviour or connectivity patterns with other brain regions may be informative ( Marquand et al., 2017 ;Prze ź dzik et al., 2019 ), though such measures remain correlational. Finally, while we discuss these issues in relation to connectopic mapping, it is possible that other analysis techniques concerned with the topographic organisation of the brain could be similarly affected.
In this study, we identified a potential flaw in connectopic mapping analyses. Local connectopic maps can be substantially biased by local spatial autocorrelations, for instance as introduced by spatial smoothing or interpolation between co-ordinate spaces, to the extent that illusory gradients can be produced even from random data. Gradients appear highly similar between real and random data, and replicate gradients previously reported in the literature. Coarser-scale whole-brain gradients can be less susceptible to local spatial autocorrelations, but nevertheless appear sensitive to specific features of the analysis pipeline. These results indicate that previously reported gradients from connectopic mapping studies may need to be interpreted with a degree of caution. Future connectopic mapping studies may aim to develop methods to address these issues, or provide a consideration of how gradients may be influenced them.

Declaration of Competing Interest
None.

Data availability
All analysis code is available on the Open Science Framework ( https://osf.io/4dpuw/ ). All MRI data were obtained from the already publicly available StudyForrest project ( https://www.studyforrest. org/ ).

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.neuroimage.2023.120228 .