Connectopic mapping with resting-state fMRI

Brain regions are often topographically connected: nearby locations in one brain area connect with nearby locations in another area. Mapping these connection topographies, or 'connectopies' in short, is crucial for understanding how information is processed in the brain. Here, we introduce principled, fully data-driven methods for mapping connectopies using functional magnetic resonance imaging (fMRI) data acquired at rest by combining spectral embedding of voxel-wise connectivity 'fingerprints' with a novel approach to spatial statistical inference. We applied the methods in human primary motor cortex and show that it can trace biologically plausible, overlapping connectopies in individual subjects with exceptionally high reproducibility. As a generic mechanism to perform inference over connectopies, the new spatial statistics approach enables rigorous statistical testing of hypotheses regarding the fine-grained spatial profile of functional connectivity and whether that profile is different between subjects or between experimental conditions. The combined framework offers a fundamental alternative to existing approaches to investigating functional connectivity in the brain, from voxel- or seed-pair wise characterizations of functional association, towards a full multivariate characterization of spatial connectopies.


Abstract
Brain regions are often topographically connected: nearby locations in one brain area connect with nearby locations in another area. Mapping these connection topographies, or 'connectopies' in short, is crucial for understanding how information is processed in the brain. Here, we introduce principled, fully data-driven methods for mapping connectopies using functional magnetic resonance imaging (fMRI) data acquired at rest by combining spectral embedding of voxel-wise connectivity 'fingerprints' with a novel approach to spatial statistical inference. We applied the methods in human primary motor cortex and show that it can trace biologically plausible, overlapping connectopies in individual subjects with exceptionally high reproducibility. As a generic mechanism to perform inference over connectopies, the new spatial statistics approach enables rigorous statistical testing of hypotheses regarding the fine-grained spatial profile of functional connectivity and whether that profile is different between subjects or between experimental conditions. The combined framework offers a fundamental alternative to existing approaches to investigating functional connectivity in the brain, from voxel-or seed-pair wise characterizations of functional association, towards a full multivariate characterization of spatial connectopies.

Introduction
Functional magnetic resonance imaging (fMRI) connectivity patterns are widely derived in an attempt to characterize the brain into patches of homogeneous connectivity (Smith et al., 2013b;Power et al., 2014), thereby providing a basis for cortical cartography, brain parcellation and subsequent network-based connectomics. However, it is well acknowledged that such a strict parcellation model is generally overly simplistic since brain organisation is known to be not piecewise constant in patches. Instead, information processing in the brain is often topographically organised: nearby locations within one brain region typically connect with nearby locations elsewhere in the brain.
For instance, the primary cortical sensory-motor areas contain somatotopic, retinotopic and tonotopic maps, and these maps are generally continued through to higher-level sensory-motor cortex via topography preserving corticocortical connections (Jbabdi et al., 2013;Kaas, 1997). Such connection topographies, which we will refer to as 'connectopies' in short, are also found in the parietal, entorhinal, parahippocampal, and prefrontal cortices, as well as the basal ganglia, cerebellum and corpus callosum (Jbabdi et al., 2013;Thivierge and Marcus, 2007).
Connectopies are thought to be ubiquitous because they follow the principles of good biological design, allowing for context-dependent computations using metabolically efficient short-range neural circuitry. Multiple, overlapping connectopies can co-exist within the same brain area, enabling complex computations, such as coincidence detection, with relatively simple spatial rules (Jbabdi et al., 2013;Kaas, 1997). Thus, connectopies are widely thought to be crucial to brain function. While connectivity-based parcellation continues to provide important information about the brain's functional anatomy, therefore, it appears that a more complete understanding can be achieved by investigations into the more fine-grained, topographic organization of functional connectivity within brain areas.
Only few attempts have been made to characterize the topographic organization of connectivity within a cortical region (Anderson et al., 2010;Buckner et al., 2011;Cauda et al., 2011;Gravel et al., 2014;Haak et al., 2013;Heinzle et al., 2011;Raemaekers et al., 2014;Taren et al., 2011;van den Heuvel and Hulshoff Pol, 2010). In those studies, the typical approach was one of gradually changing the position of a seed, thereby characterising gradual changes in ensuing functional connectivity. While this approach has successfully revealed known connectopies, for instance in sensorimotor cortex (van den Heuvel and Hulshoff Pol, 2010), it is less applicable to cortical patches with unknown topographic organization: if all goes well, the topography emerges after carefully looking at the different seed connectivity patterns linked to the seed locations, but there are many ways of traversing a region in the brain, and an exhaustive test of all possible trajectories is prohibitive. In addition, many areas are known to contain multiple modes of topographic organisation within the same cortical patch (Jbabdi et al., 2013). This presents a formidable challenge to basic seed-based connectivity approaches: if multiple overlapping connectopies simultaneously exist within the area of interest, the 'moving-seed' approach will erroneously uncover their single superposition instead of revealing the true multiplicity of modes of organisation ( Fig. 1). A further disadvantage is that smaller seeds produce noisier results, so the resolution at which connectopies can be reliably detected is limited.
Here, we draw upon the idea that all of these limitations can be overcome by reformulating the problem in terms of finding the intrinsic degrees of freedom of the high-dimensional connectivity dataset. Thus, the problem of mapping connectopies is recast as a manifold learning problem (Lee et al., 2007). Rather than moving a seed around while monitoring for gradual changes in connectivity, therefore, we instead compute the pair-wise similarities among the functional connectivity 'fingerprints' of all voxels within a pre-specified region of interest (ROI) and then employ manifold learning in order to find a limited set of overlapping connectopies. Although this strategy is broadly similar to procedures employed in connectivity-based parcellation, its aim is fundamentally different. Connectivity-based parcellation is employed to delineate brain areas of homogeneous function, while connectopic mapping aims to characterise the internal topographic organisation of connectivity within such brain regions-rather than segregating brain areas into piece-wise constant parcels, we characterise them into (non-constant) connection topographies. This novel view necessitates alternative approaches to statistical inference: we no longer assume piece-wise constant connectivity, but characterise the topographic organization of connectivity in terms of a multivariate estimate (map).
Thus, we propose to employ a spatial statistics approach known as trend surface analysis to perform inference over the topographic organisation of connectivity.
In what follows, we demonstrate (1) that the proposed data-driven approach indeed represents a useful strategy for connectopic mapping, producing biologically valid, overlapping connectopies based on resting-state fMRI with exceptionally high fidelity in single subjects, and (2) that the ensuing connectopies can be adequately characterised into a small number of parameters using trend surface analysis, thereby enabling statistical inference over connectopies. The combined framework offers a fundamental alternative to existing approaches to investigating functional connectivity in the brain-from voxel-or seed-pair wise characterizations of functional association and connectivity-based parcellation, towards a full multivariate characterization of spatial topographies that can be compared across subjects and experimental conditions using rigorous statistical hypothesis testing. Illustrated is an extreme case where one connectopy is orthogonal to a second connectopy (similar graytones indicate similar connectivity patterns). Unaware of the fact that these two overlapping connectopies underlie the measurements, moving a seed along the region results in the erroneous inference that the change of connectivity occurs along the diagonal. If one were to use these measurements as a basis for cortical parcellation, for instance into two sub-regions, they would lead to highly reproducible, yet erroneous subdivisions. That is, the regions would be split along the diagonal perpendicular to the diagonal direction of connectivity change (white dashed line).

Methods
The proposed framework consists of three fundamental elements: connectivity fingerprinting, manifold learning, and spatial statistics for inference over connectopies (Fig. 2). We demonstrate and evaluated this approach by mapping the well-known somatotopic organisation of human primary motor cortex (M1) using the resting-state fMRI data of the first 60 subjects of the WU-Minn Human Connectome Project (HCP).

Connectivity fingerprinting
Voxel-wise connectivity fingerprints were derived according to the following steps. is computed as the Pearson correlation between the voxel-wise time-series and the SVD-transformed data. This yields matrix C, whose rows convey correlation maps; one map for each voxel within the ROI.

Manifold learning
We elected to employ non-linear manifold learning using the Laplacian Eigenmaps (LE) algorithm (Belkin and Niyogi, 2003), which performed well in a precursor of our framework (Navarro-Schröder et al., 2015). Known for its computational simplicity, the LE algorithm effectively represents the initial data transformation step of spectral clustering (von Luxburg, 2007), and has previously been shown to be useful for tracing changes in probabilistic white-matter tractography connectivity (Cerliani et al., 2012;Johansen-Berg et al., 2004) and connectivity-based parcellation (CBP) (Craddock et al., 2012). Our implementation of the LE algorithm involves the following steps.
First, the between-voxel similarity of the correlation maps is computed, yielding a matrix S that characterizes the within-ROI similarity of connectivity. To compute S, we used the η 2 coefficient (Cohen et al., 2008): where µ j = (C α,j -C β,j ) / 2 and is the mean of µ across all p SVD-components. The η 2 coefficient represents the fraction of the variance in one connectivity profile that is accounted for by the variance in another, and ranges between 0 (entirely dissimilar) to 1 (entirely similar). In principle, other similarity measures such as the Pearson correlation or simply CC T could also be used. However, for these measures, negative values should first be converted to positive values so that S can be transformed into a graph with vertices that carry only non-negative weights (von Luxburg, 2007). Next, S is transformed into a connected graph represented by matrix W: where ε is defined as the minimum value required for the graph to be connected (Cerliani et al., 2012). Note that this is done so as to meet the connected graph assumption of the next step of the analysis-or the following steps will need to be performed for each connected component separately. Ensuring a single connected component at this stage of the analysis fits with the intuition that the ensuing connectopies should cover the entire ROI instead of multiple connectopies that each might cover only a restricted portion of the ROI, depending on parameter ε. Given matrix W, . , ! of the graph Laplacian are then found by solving the generalized eigenvalue problem = . As such, the eigenvectors ! , . . , ! associated with the . That is, the eigenvectors convey mappings wherein voxels with similar connectivity profiles stay as close together as possible-they convey connectopies, wherein similar values represent similar connectivity patterns.
For comparison, connectopies were also derived using the Isomap algorithm (Tenenbaum et al., 2000). To this end, S was transformed into graph G, wherein connections exist between k nearest neighbours and where k was chosen as the minimum value required for the graph to be connected. Next, Dijkstra's algorithm (Dijkstra, 1959) was employed to compute the shortest path between each pair of nodes. Finally, the ensuing distance matrix D was subjected to multi-dimensional scaling. Specifically, a matrix Q was set-up by first computing = ∘ , and then computing = − 1 2 , where = − !! ′, with n being the number voxels inside the ROI. From Q, the eigenvectors ! , . . , ! associated with the m largest eigenvalues were extracted, which were then transformed into the final embedding = / , where Ν is the matrix of m eigenvectors, is the diagonal matrix of the m eigenvalues of Q, and the columns of convey the connectopies according to Isomap.

Spatial inference
Connectopic mapping characterizes the topographic organization of connectivity in terms of a multivariate estimate (map). To enable statistical hypothesis testing on these maps, therefore, we propose a spatial statistics approach. The approach involves finding a parsimonious representation of the connectopy, governed by coefficients that can be tested either parametrically or non-parametrically or employed as features in other analyses. Finding this representation is achieved by estimating the parameters of a spatial model that describes the connectopies in terms of polynomial basis functions for the spatial trend and a Gaussian process that models more detailed spatial variation.
Thus, we approximate each connectopic map using a spatial model where the value of the connectopy at spatial location is given by: Here, ! is a connectopy, ( ) is a spatial basis function with coefficients that collectively describe the lowfrequency spatial trend; , ∼ ( , , ′ ; ) is a zero-mean Gaussian process that models more detailed spatial variation and ( ) ∼ (0, ! ! ) are spatially uncorrelated residuals. A model of this form is known as a 'trend surface model' in the spatial statistics literature (Gelfand, 2010). For the specification of the covariance function for the Gaussian process, , ′ ; , we use a Matérn covariance function of the form: where ! , and ℓ are respectively scaling, smoothness and length scale parameters and ! is a modified Bessel function (Rasmussen and Williams, 2006). This covariance function is preferred in spatial statistics over the 'squared exponential' covariance function more common in machine learning because the latter is regarded as too smooth for spatial applications (Gelfand, 2010;Wackernagel, 2003). We follow the convention in spatial statistics and fix the smoothness parameter to = 5 2, which shows good performance in spatial applications (Wackernagel, 2003) and performed well in preliminary tests. We estimated the remaining parameters ( , ! ! , ! ! and ℓ) using nonlinear conjugate gradient optimization. See (Rasmussen and Williams, 2006) for more details.

Evaluation data and ROI definition
We evaluated our data-driven connectopic mapping approach using a data-set that comprised the first 60 subjects of the WU-Minn Human Connectome Project (Van Essen et al., 2013), with two sessions of two 15-minute multi-band accelerated (TR = 0.72s) resting-state fMRI scans per individual. This 3T whole-brain dataset, with an isotropic spatial resolution of 2mm, is publicly available and has been pre-processed as detailed in (Smith et al., 2013a).
Briefly, pre-processing steps included corrections for spatial distortions and head motion, registration to the T1w structural image, resampling to 2mm MNI space, global intensity normalization, high-pass filtering with a cut-off at 2000s, and the FIX artefact removal procedure (Griffanti et al., 2014;Salimi-Khorshidi et al., 2014). For this work, we additionally smoothed the images using a 6mm FWHM Gaussian kernel and removed the mean ventricular and white-matter signal from the time-series data. After converting the time-series of each voxel to percent signal change by dividing by and subtracting its mean amplitude over time, the data from the two scans in each session were concatenated, resulting in two 30-minute functional scans per subject (one for each session day).
Primary motor cortex (M1) was defined based on anatomical criteria using the Freesurfer toolbox. Specifically, the T1w MNI template image (1mm isotropic resolution; as provided by FSL) was subjected to Freesurfer's cortical reconstruction procedure ('recon-all') to create a number of MRI volumes wherein voxels are assigned a neuroanatomical label (e.g., left precentral gyrus). The relevant volumes were subsequently converted to the NiFTI file format, resampled to 2mm MNI space and binarized.

Connectopic mapping at the group-level
To evaluate the proposed analysis framework, we applied it to the resting-state fMRI time-series from the human motor strip (M1), a brain region with a well-known topographic (i.e., somatotopic) organization, well-established topographic connectivity with the motor strip in the opposing hemisphere (van den Heuvel and Hulshoff Pol, 2010) and cerebellum (Buckner et al., 2011), and whose borders are well-defined by anatomical criteria (i.e., the precentral gyrus). The dataset comprised the first 60 subjects from the WU-Minn Human Connectome Project (HCP) (Van Essen et al., 2013). We first describe the results at the group-level, which were obtained by computing the pair-wise similarities among the voxel-wise, whole-brain, gray-matter connectivity fingerprints in the anatomically defined ROIs (separate ROIs for each hemisphere) in each of the 60 subjects and averaging these values across them. Note that this is a valid approach for pooling data across subjects because contrary to the heterogeneous input time-series data, the similarities among the connectivity fingerprints should be broadly similar across subjects. Feeding the averaged similarity scores to the LE algorithm revealed many overlapping connectopies, sorted in the order of how well how well they preserved the topology of the original, high-dimensional, functional connectivity fingerprints. Of these, we found that those with the greatest and second-greatest topology preservation exhibited a spatial profile directly related to known functional and structural features of the human motor strip.
The connectopy that exhibited the greatest topology preservation showed a clear correspondence with M1's somatotopic map (Fig. 3A), which follows the superior-medial to inferior-lateral axis of the precentral gyrus.
Notably, this somatotopic organisation of connectivity could not be uncovered by simply computing the correlation between the voxels' rfMRI time-series and the first principal component time-series of the ROI (Fig. 4). For further validation, we also assessed the connectopy's underlying connectivity patterns by colour-coding the voxels outside M1 according to the voxels inside M1 that they correlate the most with (Jbabdi et al., 2013). This analysis confirmed that the connectivity patterns underlying the first connectopy can be characterized by mirror-symmetric, interhemispheric topographic connectivity with M1's contralateral counterpart in the opposite cerebral hemisphere (Fig.   3B) as well as its topographically organised connectivity with anterior cerebellum (Fig. 5). The connectopy with the second-greatest amount of topology preservation also represents a biologically plausible result (Fig. 6), as it reflects the distance from the superior longitudinal fascicle (subdivision II), a large white-matter fibre tract connecting the caudal middle frontal and precentral gyri with the parietal lobe (Makris et al., 2005;Thiebaut de Schotten et al., 2011). Thus, the group results demonstrate that the proposed approach is capable of simultaneously tracing biologically valid, overlapping connectopies from resting-state fMRI data.  Computing individual voxel-wise correlation with a region's single principal component represents a univariate approach to characterising regional functional connectivity change; see also other seed-based methods such as the one proposed by Cohen et al. (2008). Univariate approaches do not incorporate the entire set of overlapping topographies of connectivity similarity but consider only single voxel-wise estimates. As such, they cannot disentangle overlapping modes of connectivity change that may coexist within the same area of interest, leading to biologically invalid solutions. The connectopic mapping approach described here aims at identifying multiple modes of functional connectivity change in a fully multivariate fashion and thereby effectively models secondary modes of overlapping organisation that-when resorting to univariate techniques-end up confounding the estimates and thereby reducing reproducibility.

Fig. 5. The dominant connectopy in M1 reflects somatotopic organisation in both cerebral and cerebellar cortex. (A)
Group-level fMRI task-activation maps (5 < z < 10) for the HCP motor mapping task; see  for details. (B) Projection of M1's dominant connectopy onto the portions of cerebral and cerebellar cortex that were activated during the motor mapping experiment. Cortical voxels are colour-coded according to the M1 voxels that they correlate the most with. Although the left-right asymmetry in task-activation cannot be resolved, the rfMRIbased colour-gradient follows the motor mapping results across the entire motor network, including cerebellum.

Connectopic mapping in single subjects
We next assessed our approach' capability of mapping connectopies at the single-subject level. To this end, we applied the LE algorithm directly to the individual connectivity similarity scores, separately for each of two independent scan sessions, and computed the normalized root mean squared error (nRMSE) between the ensuing connectopies. The dominant, individual connectopies showed considerable resemblance to the group result, to each other, and across the two independent sessions within the same subject (Fig. 7). Indeed, Table 1 shows that across all 60 subjects, the dominant, individual connectopies were highly similar across sessions and also showed substantial correspondence across subjects. A similar pattern of results was found for the second-dominant connectopy.
Moreover, shortening the input time-series and/or sampling every third time-point suggests that the high crosssession stability of the subject-level results also applies to scans as short as ~7.5 minutes with a TR of ~2s (Fig. 8).

Fig. 7.
Subject-level connectopies. The dominant connection-topography in M1 (y 1 ) is highly reproducible across sessions within the same subject (compare top with bottom rows). The dominant connectopy also exhibits considerable reproducibility across subjects (compare left with right panels), but to a lesser extent than across sessions within the same subject, suggesting subject-specificity.

Single-subject specificity of connectopic organization
Notably, the subject-level difference scores presented in Table 1 were roughly twice as low across sessions within the same subject as they were between subjects. This suggests that the estimated connectopies capture subjectspecific features of connectopic organization. To test this idea, we conducted a mate-based retrieval experiment. In the context of the evaluation dataset considered here-where each subject was scanned twice-this involved aiming at retrieving the matching resting-state fMRI run for each subject. We employed a stringent (exact) matching criterion where we considered a match successful if the connectopy based on any given fMRI run achieved maximal correlation with the connectopy based on the second run from the same subject. We then computed matching accuracy as the sum of correct matches divided by the total number of fMRI runs across all subjects.
The mate-based retrieval experiment was conducted using the spatial model approximations of the connectopies. To obtain the spatial model approximations, we estimated an independent spatial model of the dominant connectopies for each session of each subject using polynomial spatial basis functions of degrees 1 through 4 and selected the optimal polynomial degree using the Bayesian information criterion (Schwarz, 1978). The optimal polynomial degree yielded fits that were nearly exact (nRMSE < 10 -3 ) and was either two or three for all subjects, with a much smaller difference between these degrees relative to other values (e.g. first or fourth degree). To enable comparisons across subjects, we used a cubic polynomial to model the surface trend for all subjects.
We found that the spatial models allowed for correctly recognizing 62% and 53% of fMRI runs with the matching run from each subject (for the left and right ROIs, respectively), which greatly exceeds the chance level of 1/60 = 1.67% (p < 10 -16 ; binomial test; note that this is a much more difficult problem than standard binary classification, because chance level is 1/N, where N is the number of subjects). To establish that this identification rate ensued from differences in connectopic organisation rather than individual differences in anatomy, we repeated the matematching experiment asking whether a subject's spatial pattern of the mean BOLD signal over time (derived from the time-series without conversion to percent signal change) could be used to retrieve that subject's connectopy.
That is, even though the functional data have been normalised to MNI space, anatomical idiosyncrasies could have

Comparison with Isomap
Non-linear manifold learning algorithms fall broadly in two categories: local and global approaches (de Silva and Tenenbaum, 2003). The LE algorithm represents a local approach because it maps nearby points in the highdimensional dataset to nearby points in the low-dimensional embedding. Isomap represents a global approach because it maps nearby points in the high-dimensional dataset to nearby points in the low-dimensional embedding as well as faraway points to faraway points (de Silva and Tenenbaum, 2003). The principle advantages of local approaches such as LE are computational efficiency and broader applicability, while global approaches such as Isomap tend to yield more faithful representations of the data's global geometry. However, global approaches such as Isomap are relatively prone to "short-circuiting" when the data are noisy, which can have major negative consequences for the ensuing embedding (de Silva and Tenenbaum, 2003;Belkin and Niyogi, 2003). It is for these reasons that we focused on a local approach. However, it is still worth investigating whether the results presented so far crucially depend on this choice. We therefore repeated our analyses using an alternative 'global-approach' implementation of our framework utilising Isomap rather than LE. Comparisons between the dominant and seconddominant connectopies obtained with LE and Isomap indicated highly similar results at the group-level (nRMSE = 0.078 and nRMSE = 0.087 for the dominant connectopy in the left and right precentral gyri, respectively; nRMSE = 0.280 and nRMSE = 0.316 for the second-dominant connectopy). Similar performances could also be observed at the individual level (see Table 2), though Isomap tends to produce maps that are more homogeneous across subjects.
Finally, figure 9 shows how the results obtained with Isomap generalise to rfMRI time-series of shorter duration and a lower sampling frequency, indicating that Isomap exhibits slightly better reproducibility for ~30min of data with sub-second sampling, but poorer reproducibility for more conventional scan durations (5-20min) and/or TR.

Discussion
We have demonstrated that it is possible to map biologically meaningful connectopies with resting-state fMRI in a principled, fully data-driven manner. This innovation builds on previous characterizations of probabilistic white matter tractography change using diffusion imaging (Cerliani et al., 2012;Johansen-Berg et al., 2004) as well as the established utility of spectral techniques in connectivity-based parcellation (CBP) (Craddock et al., 2012). It is important to note, however, that connectopic mapping cannot be perceived as an instantiation of CBP. The crucial distinction is that connectopic mapping exposes the fine-grained topographic organization of a brain region's connectivity, whereas CBP aims to yield strictly segregated cortical patches of piece-wise constant connectivity. We demonstrate that this shift from a piecewise constant model of connectivity towards a map-wise characterisation of change in connectivity results in biologically interpretable modes of organisation. This fundamental shift necessitates novel inference procedures and we have introduced a trend-surface modelling approach that can accurately condense the high-dimensional connectopy imaging phenotype into a low number of trend coefficients.
Previous work on connectivity-based cortical cartography has focused almost exclusively either on parcellation or used a multitude of separate seeds. The present and previous characterizations of topographically organized connectivity show that CBP discards potentially meaningful information about the functional organization of the area under investigation-information that may in actual fact refute the assumption of homogenous connectivity, and thus the validity of the ensuing parcellation (Smith, 2012). Moreover, CBP approaches generally generate parcellations using an aggregated, multi-dimensional representation of the connectivity data, but this can lead to biologically invalid solutions when multiple overlapping connectopies simultaneously exist within the same area (see Fig. 1). Thus, while connectopic mapping and CBP can both be based on connectivity fingerprinting and manifold learning, only connectopic mapping provides a detailed, biologically plausible picture of cortical functional organization.
Here, our evaluations are largely based on the LE algorithm for manifold learning. In principle, however, one might also rely on one of various alternative techniques such as kernel principal component analysis (KPCA) (Scholkopf et al., 1998), isometric feature mapping (Isomap) (Tenenbaum et al., 2000), locally linear embedding (LLE) (Roweis and Saul, 2000), or the more recently introduced structural LE approach (Lewandowski et al., 2014). Each of these techniques has its own benefits and limitations and it is difficult to predict which is fit for the dataset at hand. Given the context of resting-state fMRI data, we chose to focus on the traditional LE algorithm because it has previously been shown to be effective in the context of tracing changes in white-matter tractography connectivity (Cerliani et al., 2012), its computational simplicity, and because of its close connection to spectral clustering (Belkin and Niyogi, 2003), which is widely applied to characterize resting-state fMRI connectivity patterns. Indeed, the biological plausibility of the present results demonstrate that the effectiveness of the LE algorithm also applies to capturing the fine-grained topographic structure of functional connectivity using fMRI data acquired at rest.
Furthermore, a comparison with the connectopies ensued by Isomap (Tenenbaum et al., 2000) suggests that the LE algorithm or other local non-linear manifold learning approaches should be preferred when applying the proposed connectopic mapping methods to datasets obtained using more conventional acquisition protocols, which is likely due to the fact that LE is less prone to "short-circuiting" when the data are noisier.
To validate the biological plausibility of the ensuing connectopies across information processing hierarchies, we projected the connectopies onto the rest of cortex by colour-coding voxels outside the ROI according to the voxels inside the ROI that they correlate the most with (Jbabdi et al., 2013). This simple analysis revealed that the primary motor cortex in one hemisphere of the brain is topographically connected with the opposing hemisphere as well as anterior cerebellum, replicating previous work (Buckner et al., 2011;van den Heuvel and Hulshoff Pol, 2010). This result not only demonstrates the biological plausibility of the estimated connectopies. It also indicates that the procedure of first performing connectopic mapping in one brain region and then projecting that map onto the rest of the brain is an effective generic approach to discover new topographically connected information processing networks. The discovery of such brain networks could be important to better understand the neural underpinnings of various perceptual and cognitive functions, e.g. by mapping such hierarchical organisation in areas beyond simple sensory cortices. In this context, one might also consider augmenting the connectopic mapping framework with more sophisticated approaches to mapping topographic organisation across information processing hierarchies, such as connective field modeling (Haak et al., 2013).
The results further demonstrate that the proposed analysis framework can map the connectopic organization of a brain region in single subjects with an exceptionally high cross-session reproducibility rate. Thus, it is possible to meaningfully compare connectopic maps between subjects and across conditions. In addition, the high cross-session reproducibility rate suggests that-at least within the context of motor cortex-the proposed connectopic mapping framework is relatively robust to inaccurate ROI definitions. That is, to generate group results and to allow for cross-subject comparisons, the subject-native functional data were normalized to standard MNI space and a single ROI was used for all subjects. This ROI approximates the locality of an individual's true motor strip, and so the fact that we were able to correctly trace the connectopic organization of this region in the vast majority of subjects with high cross-session reproducibility suggests that it does not require a very precise delineation of the area under investigation. However, this may be different for areas that are more variable in their shape and location. In principle, the effect of inaccurate ROI definitions will generally be limited if the ROI extends only a little beyond the area of interest, but if the ROI includes extended portions of multiple areas with distinct connectivity patterns, the dominant connectopy will reflect this distinction (Eickhoff et al., 2015).
To allow for statistical inference over the ensuing connectopies, we propose a spatial statistical approach known as trend surface analysis (Gelfand, 2010). Trend surface analysis is a standard technique in the geosciences that is used to test hypotheses about the spatial variation of, for instance, physical geography, rainfall, temperature, political climate and so on. Here, we adopted trend surface analysis as a generic approach to parameterize the connectopic mapping results, thereby opening up the possibility to test hypotheses about the spatial variation of functional connectivity and whether that spatial variation is different between subjects and experimental conditions. For instance, having access to a low-parametric characterisation of connectopic organisation offers the opportunity to explicitly formulate and test anatomically relevant hypotheses beyond what can be achieved using traditional voxelor cluster-wise testing procedures such as testing if the gradient of connectivity follows a certain anatomical orientation. Indeed, this can also be used to explicitly test the implicit assumption of CBP approaches that brain areas exhibit piece-wise constant organisation. It also affords an economical description of the spatial variation of functional connectivity for purposes such as classification. These features are likely to be of great interest to the connectivity-based cortical cartography community, as methods to perform statistical inference over the spatial layout of the brain's functional anatomy have thus far been markedly lacking.
Previous work on probabilistic tractography data acquired with diffusion imaging has also applied manifold learning in an attempt to find anatomical gradients of connectivity (Cerliani et al., 2012;Johansen-Berg et al., 2004).
However, white-matter tractography is limited in the ability to trace the precise site-to-site connections required for mapping connectopies (Jbabdi et al., 2013;Jbabdi and Johansen-Berg, 2011). Resting-state fMRI measures connectivity directly at each voxel and these measurements are one of function. The ability to map and perform inference over resting-state fMRI connectopies thus opens up a wide range of novel research opportunities. For instance, our framework can be readily employed to test for the presence of topographic maps in association cortex, which was not possible in the past because association cortex represents information that cannot be easily targeted with stimulus-driven or task-based experiments. Testing for resting-state connectopies instead could offer a first time glance into these regions' functional organizational structure and drive future studies for better understanding the fundamental nature of the regional computations. Likewise, connectopic mapping provides a translational avenue to patients with impairments that preclude stimulus-driven or task-based experiments. Because topographic maps are widely thought to be crucial to healthy brain function, connectopic mapping also holds great promise for developing more sensitive markers of disease, advancing both cognitive and clinical imaging neuroscience.