Functional Parcellation of Human Brain Using Localized Topo-Connectivity Mapping

The analysis of connectivity between parcellated regions of cortex provides insights into the functional architecture of the brain at a systems level. However, the derivation of functional structures from voxel-wise analyses at finer scales remains a challenge. We propose a novel method, called localized topo-connectivity mapping with singular-value-decomposition-informed filtering (or filtered LTM), to identify and characterize voxel-wise functional structures in the human brain from resting-state fMRI data. Here we describe its mathematical formulation and provide a proof-of-concept using simulated data that allow an intuitive interpretation of the results of filtered LTM. The algorithm has also been applied to 7T fMRI data acquired as part of the Human Connectome Project to generate group-average LTM images. Generally, most of the functional structures revealed by LTM images agree in the boundaries with anatomical structures identified by T1-weighted images and fractional anisotropy maps derived from diffusion MRI. In addition, the LTM images also reveal subtle functional variations that are not apparent in the anatomical structures. To assess the performance of LTM images, the subcortical region and occipital white matter were separately parcellated. Statistical tests were performed to demonstrate that the synchronies of fMRI signals in LTM-derived functional parcels are significantly larger than those with geometric perturbations. Overall, the filtered LTM approach can serve as a tool to investigate the functional organization of the brain at the scale of individual voxels as measured in fMRI.

T O UNDERSTAND the principles and mechanisms of complex brain functions at a systems level, computational models seek to analyze and interpret interactions of components within neurobiological networks [1]. To date, most studies of brain networks have focused on the patterns and dynamics of functional connectivity between parcellated regions of cortex by measuring statistical associations in blood oxygenation level-dependent (BOLD) MRI signals that are related to neural activity [2]- [4]. The region-level connectivity is commonly represented by a correlation matrix or a topological graph that consists of nodes and edges, which provides an elegant framework for assessments of region-level connectivity profiles, including the 'action-at-a-distance' connections between non-adjacent elements. While numerous studies have demonstrated the effectiveness of the region-level connectivity analysis in network neuroscience [5]- [7], this approach is unable to capture the spatial organization of individual tissue elements below the level of the selected parcellations. An alternative strategy is to derive functional interactions within 3D brain images using prior knowledge of neural anatomy, assuming brain function and structure are intrinsically related, thus offering the possibility of anatomically constrained, voxellevel analysis of functional connectivity.
Voxel-level functional connectivity has the potential to identify intricate functional structures from groups of individual tissue elements (see [8], [9] for examples), but in practice it is very challenging due to the limited signal-noise-ratio of fMRI signals and the difficulties of capturing the spatial patterns of functional interactions among the very large number of voxels that could be involved. Recently, principles from the emerging field of graph signal processing (GSP) have been applied in novel voxel-level analyses [10]- [14]. Under the GSP framework, statistical associations in BOLD signals can be characterized by a topological graph where nodes and linking edges represent image voxels and functional connections between them; mathematical constructs based on graph Laplacians can then be derived, e.g., eigen vectors, or continuous gradients, which depict patterns of spatial continuity of the connectivity in the lattice space [11], [15]. The GSP framework has demonstrated initial successes in delineating functional boundaries within the human brain, which can be further used for parcellations of a functional atlas [12], [16]. However, artificial functional boundaries that often arise from specific geometries of the regions of interest may pose a significant challenge to its applications [11]. Some strategies have been proposed to remedy this inherent limitation, but with no proof of their effectiveness in general [10], [11]. Hence, developing new and more robust methods would be valuable, especially considering that additional information on connectivity at a finer scale is potentially obtainable.
In deriving functional structures based on voxel-level topographs, two main technical issues, critical to reliable boundary delineations between distinct functional structures, need to be considered. First, signal components shared by separable but adjacent functional structures may create cross-boundary connections, which tend to blur the functional boundaries delineated. Hence, appropriate signal filtering is required to enhance the functional boundaries so that the signal components shared by multiple functional structures in raw BOLD data are suppressed. Second, as emphasized in reference [17], fMRI-based brain parcellations mainly rely on two connectivity features: local connection and global similarity, which are inherently dichotomous to each other and thus capture complementary aspects of brain connections. On one hand, the local-connection-based parcellations [18]- [23] exploit abrupt changes of functional connectivity profiles between adjacent voxels, and previous work suggests that approaches of this kind are especially suitable for parcellating cortical areas due to the 2D nature of the brain surface. Moreover, the parcels derived using local connection features basically contain spatially contiguous voxels; this is consistent with neuroanatomic concepts that brain tissue is rich in dense local connections which tend to form individual functional parcels. On the other hand, the global-similarity-based parcellations [24]- [28] cluster voxels that bear similarities in fMRI time courses or global connectivity patterns. As such, the spatial continuity of neuroanatomy mentioned above is disrespected to ensure temporally coherent voxels are grouped into a single functional ensemble, which even has 'action-at-a-distance' functional connections. These two important features have been widely exploited for the purpose of brain functional parcellations or deriving other regional descriptors [29]- [31], but to date few studies have made full use of them by integrating both into a single framework. There was one preliminary attempt that combines both the features to delineate boundaries between functional areas on a 2D cortical surface [17], but the extension to 3D space is not trivial.
The motivation of this study is to establish a novel GSPbased framework that takes the above technical issues into consideration, with the goal to functionally parcellate the brain cortex, subcortex and white matter in 3D space simultaneously. First, we leverage singular value decomposition (SVD) to adaptively suppress signal components in fMRI time series that are shared by multiple functional structures. Second, we define functional boundaries by integrating information encoded in local connections and global similarities with the following two constraints. First, the connections of a vertex are confined to a small neighborhood to emphasize the local connections in the boundary delineation. Second, the connection of two vertices in the neighborhood is further modulated by the similarity of their connections with all other vertices involved in the parcellation. Finally, virtual contrasts are reconstructed for the graph processed with the above procedures to map the localized topo-connectivity, which is built upon a rationally designed contrast mechanism that is first proposed on the framework of GSP in this work. We refer to the entire process as localized topo-connectivity mapping with SVD-informed filtering, and filtered LTM for brevity.
In the following, we first provide a detailed mathematical formulation of filtered LTM and describe simulation experiments that illustrate its performance. Then, in vivo studies are reported using resting-state fMRI data acquired by a 7T scanner in the Human Connectome Project (HCP). Functional structures displayed by the resulting LTM images are compared with anatomical structures conveyed by T 1 -weighted (T1w) images and FA (fractional anisotropy) colored maps derived from diffusion MRI. Furthermore, to assess the ability of LTM images to characterize subtle functional structures, we provide two demonstrations of applications. First, LTM images are used to parcellate atlases of subcortical nuclei, which has been proven to be challenging in a recent study [10]. Second, functional structures are parcellated within the occipital white matter. For the two parcellations, rigorous statistical testing of hypotheses is performed to determine whether the synchronies of BOLD signals in LTM-informed parcellations are larger than expected due to chance, where random parcellations are used to estimate a distribution for homogeneity under a null hypothesis.

II. METHODS
The proposed LTM algorithm consists of two major steps in which filtering of noise and multiple-structure-shared signal components are followed by an LTM reconstruction ( Fig. 1(a)). In the first step, raw time-courses in each volume of interest (VOI) are processed by the singular value decomposition to generate a self-adaptive filter, which is exploited to suppress noise and global signal components in the raw time-courses. In the second step, a contrast metric is derived to capture and visualize spatial continuity patterns of localized topo-connectivity. The contrast mechanism is formulated using differential equations, which are transformed into matrix form to yield analytical solutions, or LTM reconstructions.

A. Filtering of Noise and Global Signal Components
The fMRI time-course data from a defined VOI are rearranged into a matrix (X) of dimension T × N, where T denotes the number of time frames and N the number of voxels within the VOI. SVD of matrix X is then given by where U (T × T) and P (N × N) are unitary matrices, and diagonal entries (σ i , σ i = i,i ) of the rectangular diagonal matrix (T × N) are the singular values of X in a descending order of magnitude. The columns of U (U i ) provide a complete set of orthonormal bases, whereby the time-courses can be The raw time-courses are finally projected on to the filter space F to suppress noises and multiple-structure-shared signal components. The second major step of the algorithm, LTM reconstruction, is described as follows. First, a similarity matrix of F-fingerprint S, defined as pair-wise Pearson correlations among N voxels in the VOI, is computed (see text). Then, a mask matrix M is derived from S with localizing and pruning processes to enhance boundary features in connectivity profile. Subsequently, M is applied on the initial connection matrix C that is calculated from the noise-filtered signals X U trun , yielding a boundary-enhanced matrix W. Finally, LTM images are reconstructed with a rationally designed contrast mechanism to capture continuousness patterns of the localized topoconnectivity encoded in the matrix W. (b), singular values plotted in the descending order of magnitude. These values are calculated from fMRI time-courses of the cortical voxels extracted from the dataset of an HCP participant. (c), projection coefficients of three typical time-courses on the orthogonal space U. The projection coefficients on the first orthogonal bases representing multiple-structure-shared signal components have significantly larger magnitudes than others. The projection coefficients with the indexes from 160 to 900 are very small, which indicate noise components. (d), projection coefficients of the three time-courses on the filter subspace F. Here, the dominance of the projection coefficients on the first orthogonal bases have been reduced, which is showcased by two data points indicated by the green arrows. Besides, the noise components have been removed.
represented as a linear combination of the bases, i.e., X U = U T X. This representation allows an exact reconstruction, that is, the measured signals can be recovered exactly as X = UX U . Note that each base vector in U is associated with a singular value in (U i → σ i ), and the base vectors in U are arranged in the same descending order of singular values. An important feature encoded by each base vector in U is that the base vectors corresponding to larger singular values represent these signal components shared by the voxels in the VOI with larger total weights, compared with the base vectors corresponding to smaller singular values. A typical set of singular values computed from the cortical region of human brain is plotted in descending order of magnitude in Fig. 1(b). Based on SVD-informed component analysis, a self-adaptive filter can be derived to suppress noise and signal components in raw time-courses that are shared by multiple functional structures. Fig. 1(c) shows three typical time-courses (each x i , a column of matrix X) projected on to the orthogonal bases (U), where projection coefficients on the first orthogonal bases have significantly larger magnitudes than others. It is well established that neural activities are spatially integrated into multi-regional functional networks, which could lead to blurring effects on the delineation of functional boundaries between two neighboring functional structures that share common signal components. In principle, these signal components project with relatively large coefficients on to the bases associated with larger singular values. In contrast, signal fluctuations caused by noise should have low spatial coherence, and therefore they tend to project on to the bases associated with smaller singular values. To filter the noise in the time courses, we remove the projection coefficients on the bases associated with small singular values, such that the raw signal matrix X is projected on to a subspace (U trun ) that consists of only the bases associated with large singular values. The filtering of noise can be written as where , and L is determined by taking the maximum k subject to the constraint where k is an integer less than T; and ξ is a selected threshold parameter related to noise level. Subsequently, multiplestructure-shared signal components with large singular values in the SVD are suppressed because these components would contribute the undesirable back-ground connectivity and weaken the conspicuity of functional boundaries. Thus, the matrix X is projected on the subspace U trun and then multiplied with the matrix the resulting projection coefficients associated with large singular values are rescaled to relatively small values. As shown in Fig. 1(d), the dominance of the projection coefficients on the first orthogonal bases have been reduced. Filtering of the noise and suppressing the multiple-structure-shared signal components simultaneously can be reformulated by where F = U trun −1/2 trun . We refer to X F as an F-fingerprint matrix. The columns of F (F i ) are a set of orthogonal bases, subject to F T i F j = 0, (i = j ) and Note that the variance of projection coefficients of X on the i -th base of U (U i ) is σ i , and the variance becomes √ σ i on the i -th base of F (F i ) and thus the sorting order of the variances is maintained. By this operation, the dominance of the projection coefficients on the first orthogonal bases are reduced, but the sorting order of variances is preserved so that the bases of F associated with large singular values still carry large weights in characterizing connectivity profiles. This avoids the potential that signal components that account for only a small portion of total variance in the connectivity characterization could become dominant.

B. LTM Reconstruction
To reveal the functional connectivity between voxels in 3D space using the GSP framework, a voxel-level topological graph is derived based on the use of the F-fingerprints, or the filtered time-courses. We denote the graph as G = {V, E, W}, which consists of a finite set of vertices V with |V| = N, a set of edges E, and a weighted connection matrix W. The entry W i, j represents the connection strength of edge e = (i, j ) that connects vertices i and j . On the vertices of the graph, a vector f ∈ R N is defined to represent the intensity of the LTM image, where the i -th component of the vector, f (i ), represents the value at the i -th vertex in V. Note that the fMRI time-courses on vertices are measured data while the 1D vector f is a set of auxiliary quantities/signals to be derived computationally (see details below), which is finally represented in 3D space as virtual contrasts to characterize spatial continuity patterns of localized topo-connectivity. An initial connection matrix C of a graph is computed based on the noise-filtered signals (X U trun ) associated with pairs of voxels in the VOI, which is derived as To enhance boundary features encoded in the denoised connectivity matrix C, sparsification is performed on C to yield the final weighted connection matrix W as follows. First, similarity between each voxel pair is calculated by Pearson correlation of their F-fingerprints such that pairwise similarities among N voxels can be represented by a symmetric matrix is S ∈ R N×N . The matrix S represents the forementioned global similarity of the modulated fMRI signals that measure the temporal coherence between voxels without the constraints of spatial distance. Note that a groupconsensus similarity matrix can be obtained by averaging over individual subjects aligned in a common space to achieve topoconnectivity mapping at the group level. A mask matrix M that contains boundary information in local connections and global similarities can then be derived from S by applying the following processes. First, to emphasize local connections in the boundary delineation, the connections of vertex i are confined in a neighborhood N i defined by the cubic lattices of size 5 × 5 × 5 centered at the vertex i. Second, the connections between two voxels are further pruned to integrate global similarities in the boundary delineation if their Euclidean distance defined by S i. − S j. is greater than ε, where ε is defined as the minimum value required for the graph to stay connected [32]. Formally, the localizing and pruning process can be formulated as where S i. and S j. are the i -th column and the j -th column of S that represent the vectors of similarity between vertex i and j and all the vertices in V, respectively; and S i. − S j. denotes the Euclidean distance between the two vectors. Finally, this mask matrix is applied on the connection matrix C to yield the weighted connection matrix W to be used for the LTM reconstruction, which is derived as At the beginning of LTM reconstruction, the auxiliary signals f are initialized to the state f 0 (i ) = c (i = 1, 2, . . . , N), hence initial values of LTM images are the same at all positions. Then, f evolves with time to yield LTM images, in which image contrasts are modulated by a rationally designed mechanism such that the resulting image contrasts convey the spatial continuity patterns of the localized topoconnectivity that is encoded in the matrix W. The rationally designed mechanism of the time-dependent evolution includes two principles: a spatially uniform decay and local connectivity-specific propagation. The former takes account of auxiliary signal decay at a constant decay rate r for all vertices. The latter requires that the auxiliary signal on the vertex i propagates at a localized-connectivity-dependent rate. Note that the voxels belonging to a connectivity community would tend to have similar values as propagations of signals within the voxels are coupled with each other. With these two principles, the time-dependent evolution of the signal on a single vertex i can be described by a differential equation, The time-dependent evolution of the signals on all the vertices of the graph can be formulated in a compact form as where J = W − r · I, and I is an identity matrix of N × N. The graph operator J is a real symmetric matrix, thus it has a complete set of orthonormal eigenvectors {u } =1,2,...,N and associated eigenvalues {λ } =1,2,...,N , satisfying Ju = λ u . Based on these eigenvectors, the graph signal f can be transformed into a spectral representation, N). (11) This spectral representation allows an exact reconstruction, i.e., the signal f can be recovered as Accordingly, the analytic solution of Eq. (10) with an initial state signal f 0 (t = 0) can be derived as Given an appropriate value of parameter r and the initial state f 0 (i ) = c (i = 1, 2, . . . N), the signalf (t) at time point t conveys the image value for characterizing the spatial continuity pattern of the localized topo-connectivity. Based on analysis of Eq. (9), it can be inferred that the image intensity of a vertex depends on the total connecting weight of this vertex. As a result, the values of isolated vertices and those between connectivity communities (e.g., functional networks) would have lower intensities than those within a connected vertex community. Hence, this feature would lead to contrast valleys along the underlying boundaries between two connectivity communities, which naturally allows the boundaries in reconstructed images to be visualized. In addition, the propagation of signal at a vertex is modulated by the total connecting weights with respect to its 5 × 5 × 5 neighbors as well as signal intensities of these neighbors, such that the propagations of signals are spatially coupled and thus the vertices within a community tend to have similar values. Furthermore, vertices across communities tend to have different values because different vertex communities tend to have specific connectivity levels. Therefore, the image values at the community boundaries will lead to rapid contrast changes across the communities, which offers an additional mechanism to distinguish different communities. Taken together, the two mechanisms are involved in delineations of functional boundaries, namely 'valleys' and 'cliffs', where 'valleys' represents low contrasts on the connectivity boundaries, and 'cliffs' represents rapid contrast changes across the connectivity communities.

C. Proof of Concept With Simulated Data
Simulated 4D data were generated to demonstrate the concept of characterizing spatial patterns of localized topo-connectivity using the proposed filtered LTM reconstruction. The data were simulated with spatial patterns of topo-connectivity specified as follows. Briefly, time-courses with 900 time frames were assigned to a 3D digital phantom of dimensions 81×81×5, along the first two of which artificial spatial patterns were synthesized. The synthesized 2D spatial patterns were defined by constructing ten areas such that time-courses within each area were identical and across areas were different. Subsequently, uncorrelated noise time-courses were added to the time-courses of the synthetic patterns to achieve a signal-to-noise ratio (SNR) of 20. The correlational relationships of the time-courses within each area and between the areas were then computed and represented by a correlation matrix. Finally, the proposed filtered LTM reconstruction was implemented on the simulated data, to validate its ability to recover the spatial patterns of synthesized topo-connectivity. Parameters of the simulation experiment were identical to those for the in vivo experiments, which are described in Section II in Supplementary material.

D. Filtered LTM Reconstruction of Human fMRI Data
1) Data and Preprocessing: Human fMRI data were sourced from the Human Connectome Project (HCP) database and included part of the HCP S1200 cohort which involved healthy young adults ranging between 22 and 37 years old.
HCP datasets used in this study included 153 participants collected on a Siemens Magnetom 7T MR scanner that passed quality control criteria. HCP datasets used in this study include three imaging modalities: Resting-state fMRI, T 1 w MRI and Diffusion MRI. Parameters of the three MRI scans are described in [33]. Minimally preprocessed volumetric fMRI data were downloaded from the HCP repository. The minimal preprocessing pipeline includes removing spatial distortions, realigning volumes to compensate for subject motion, registering the fMRI data to the structural T1w images, reducing the bias field, normalizing the 4D image to a global mean, and masking the data with a final brain mask. Further preprocessing before the filtered LTM reconstruction included three basic procedures as follows. First, spatial smoothing was performed with a Gaussian smoothing kernel of 4-mm full width at half maximum. Second, time-courses were bandpass filtered to include frequencies spanning from 0.01 to 0.1Hz. Finally, the means of the filtered time-courses were subtracted, and the data then normalized to unit variance. To investigate the influences of global signals, nuisance signals arising from cerebrospinal fluid and head motion on LTM reconstructions, we regressed out these signals and compared the results with our basic procedures in Section III in Supplementary material. In this work, the fMRI dataset was separated into two groups: the first 140 subjects for a primary dataset of LTM reconstruction and the remaining 13 subjects for an independent validation dataset.
2) LTM Reconstruction of Functional Volumes: The filtered LTM reconstruction on a group level was conducted on selected brain cortex and a composite region including white and subcortical gray matters separately, which are described in Section I in Supplementary material. The key parameters of the LTM reconstructions were determined according to Section II in Supplementary material. In this study, we compared functional structures derived from LTM images with anatomical structures conveyed in FA colored maps using the following rationale. First, compared with other modalities of MRI, the DTI-derived FA colored maps can characterize the finest structures in WM, which show clear anatomical boundaries between WM and GM, and between different fiber bundles/tracts. These anatomical boundaries are informative for the parcellations of subcortex and occipital white matter that were implemented in our studies. Second, previous studies [34]- [36] have shown the spatial organizations of WM functions are relevant to the anatomical structures at the tract level, so that anatomical structures conveyed in FA colored maps are informative for the validity of the LTM-derived parcellation in white matter. FA colored images used as a reference were collected from only one subject since accurate registration is challenging due to distortions stemming from eddy currents and inhomogeneities of the main magnetic field in DTI acquisitions [37], [38]. In addition, the group-consensus T 1 w images were obtained by averaging 140 T 1 w images from subjects in the primary dataset after registration to the MNI152 space, which also provide an anatomical reference.

3) Parcellation of Subcortical Nuclei and Occipital White
Matter: Based on the LTM-derived functional structures, parcellation was performed for subcortical nuclei and occipital A key marker of parcellation validity is the homogeneity of a parcel, which was defined as the average temporal correlation coefficient between all pairs of voxels in that parcel. (a), statistical testing with the null hypothesis. A violin plot is used to display inferential statistics. In the violin plot, red dots indicate homogeneity of LTM derived parcels. Bottom and top edges of the narrow box in violin plot indicate 25th and 75th percentiles of the distribution, respectively. The central mark indicates the median. The whiskers extend to the most extreme data points that are not considered outliers (1.5 × interquartile range). To determine whether the homogeneity of LTM-derived parcels was larger than geometrical perturbed parcels, 200 perturbed parcels are created to estimate a distribution of homogeneity as that were used for the null hypothesis in statistic testing. b-d illustrate the principle of geometrical perturbation in a 2D scenario. (b), ground truths of five functional parcels (R1-R5). Colored regions show the five functional parcels, among which some voxels that are not belong to any parcels are indicated by the blank. (c), three perturbed parcellations generated under geometrical constraints (P1-P3). It is obvious that a parcel that agree with the ground truth of R5 has the largest value of homogeneity compared with the perturbations. (d), two perturbed parcellations generated without geometrical constraints (P4&P5). P4 shows a parcel that cover multiple functional structures and potentially has a low homogeneity, which could lead to false positive errors in statistical test; and P5 shows a parcel that coincidently overlap with an unintended functional structure thus potentially has a high homogeneity, which could lead to false negative errors in statistical tests.
white matter separately. A key marker of parcellation validity is within-parcel homogeneity, which was defined as the average temporal correlation coefficient between all pairs of voxels in the parcel. A higher homogeneity value indicates greater similarity in time series within the region. To validate the LTM-informed parcellation, we compared its homogeneity value with that derived from geometric perturbations, as illustrated in Fig. 2. Colored regions in Fig. 2(b) show the ground truth of five functional parcels, in between which voxels were not assigned to any of the parcels and thus left blank; three dashed curves in Fig. 2(c) enclosed parcels derived from the ground truth parcel in the center (pink) with constrained deformations that will be described later. It is obvious that a parcel that agrees with the ground truth has a higher value of homogeneity than those with geometric deformations. Hence, our parcellations could be validated by using statistical tests that demonstrate the LTM-derived parcels have significantly higher homogeneities.
In this work, the geometric perturbations were achieved by the following steps: (1) dilating the LTM-derived parcels by three layers of voxels to generate a confining mask; (2) placing a seed voxel randomly in the LTM-derived parcels; (3) growing voxels randomly from the seed voxel in the confining mask, until the total number of voxels matched that of the LTM-derived parcels. More specifically, new voxels iteratively grew in immediate neighbors of existing voxels with equal probabilities. In each iteration, only one new voxel was recruited, and the voxel would be removed if outside the confining mask. Any holes within the grown volume would be filled using simple morphological operations. To obviate incorrect spatial correspondences between LTM-derived parcels and those from random perturbations, the distorted volume would be discarded if it overlapped with less than half of the LTMderived parcels. The constraints in the perturbation can be justified by two typical parcels derived from perturbations without the geometric or spatial constraints (see Fig. 2(d)). P4 is a parcel that covers multiple functional structures with comparable areas and thus has a low homogeneity, which could lead to false positive errors in statistical tests; and P5 is a parcel that coincidently overlaps with an unintended functional structure and so potentially has a high homogeneity, which could result in false negative errors in statistical tests.
Finally, a statistical test against the null hypothesis, i.e., the distribution of synchronies estimated by the 200 perturbed parcels, was performed using an independent validation dataset of 13 subjects. A schematic of the related inferential statistics is given in Fig. 2(a). Briefly, for each parcel, the overall homogeneity was averaged cross the 13 subjects. The fMRI timecourses in LTM-informed parcellations (alternative hypothesis) were deemed synchronous if the null hypothesis was rejected (P<0.05).

A. Experiments With Simulated Data
Proof-of-concept results with simulated data are demonstrated in Fig. 3, which shows the patterns of the simulated areas, the effects of SVD-informed filtering, and the reconstructed LTM images. In Fig. 3(a), the spatial patterns of the simulated phantom were defined by ten non-overlapping areas (denoted alphanumerically and in different colors), each of which was constructed with synthetic time-courses that encoded connectivity patterns. The time-courses across the ten color-encoded areas were correlated at various levels, which are displayed in Fig. 3(b) by a connectivity matrix representing Pearson correlations between the raw time-courses. Based on this correlation or connectivity matrix, the ten areas were separated into three groups, namely, {A1, A2, A3}, {B1, B2, B3} and {C1, C2, C3, D}, with specific properties including: each area within {A1, A2, A3} had a different intra-area connectivity profile; each area in {A1, A2, A3} and {B1, B2, B3} had negligible inter-area connectivity with all its adjacent areas; the time-courses bore common components in {C1, C2, C3, D}, and thus these areas had relatively large connectivity with each other, which is shown in the pink box of higher connectivity in the lower right corner in Fig. 3(b). The final LTM images reconstructed without and with SVD-informed filtering are shown in Fig. 3(e)&3(f), respectively. From both reconstructions, the ten areas with predefined patterns could be identified using the 'valleys or cliffs' for boundary delineation. Both reconstructed images show that contrasts are proportional to the connectivity level, which is particularly pronounced in the area group {A1, A2, A3}. The effects of SVD-informed filtering are demonstrated in Fig. 3(c), which shows the similarity matrix computed from the pairwise correlations of filtered time-courses, or Ffingerprints. It can be observed that, after SVD-informed filtering, the inter-area connectivity between the areas in {C1, C2, C3, D} have been reduced to relatively low levels. The threshold parameter ξ for the filtering was determined by taking relatively large singular values, i.e., the first ten singular values int this case, as seen in Fig. 3(d). In the LTM image reconstructed without filtering (Fig. 3(d)), there are obvious blurring effects on boundaries of the areas in the group {C1, C2, C3, D}, which is presumably induced by the common components these areas share. In Fig. 3(f), the blurring effects on the boundaries of the same areas are eliminated by the filtering. An additional nonnegligible phenomenon present in both reconstructed LTM images is that the average LTM contrast of each area in {C1, C2, C3, D} and {B1, B2, B3} is positively correlated with the number of voxels it contains. This phenomenon arises because the signal propagations of voxels belonging to a connectivity Fig. 4. Group-average LTM images reconstructed from 140 participants, group-average T1w images obtained by averaging over the 140 participants and FA maps from one representative participant. Functional structures revealed by LTM reconstruction are compared with anatomical structures conveyed by T 1 relaxation and diffusion anisotropy. In the LTM images, functional boundaries can be identified in the LTM images using the 'valleys or cliff' mechanisms. Furthermore, based on delineation of the functional boundaries, connected vertices can be parcellated to yield connectivity communities, or functional areas. Some typical connectivity communities are indicated by red arrows and concomitant Arabic numerals (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19), which are analyzed in the main text.
community are coupled in LTM reconstruction. Specifically, the voxels on boundaries with poor connections could reduce image intensities in the whole connectivity community, which is dependent on the proportion of voxels on boundaries in the connectivity community. Fig. 4 shows group-average LTM images reconstructed with the fMRI dataset collected from 140 participants in HCP, as well as group-average T 1 w images and FA colored maps from a representative participant. In the LTM images, the patterns of 'valleys or cliffs' show functional boundaries, so that functional structures can be further identified based on the boundary delineations. In Fig. 4, several selected functional structures are indicated by red arrows and numerical labels. For the cerebral cortex, these functional structures correspond with anatomical structures in grey matter as apparent in the T 1 w images. Two functional structures indicated by the labels 1 and 3 vhave similar boundaries with anatomical structures within the T1w images. Some typical functional structures of white matter are indicated by the labels {2, 4, 6, 9, 10, 12, 15}. It is noticeable that several white matter regions exhibit low intensities in LTM images, which implies that the voxels in these regions do not belong to any functional structures. In the cerebral subcortex, as indicated by the labels {5, 7, 13, 18, 19}, the boundaries of functional structures in LTM images grossly correspond with the boundaries of anatomical structures that can be identified in T 1 w images, but subtle structures within the LTM images (e.g., three subdivisions of the thalamus indicted by {7} can be appreciated in Fig. 6(a)) are not seen in the T 1 or FA colored maps. In the cerebellum, LTM images show functional structures labeled as {8, 11, 14, 16, 17}. Overall, functional structures derived from LTM images exhibit similar boundaries with anatomical structures depicted in the T 1 w images and the FA colored maps. For additional visualization, we display LTMbased connectivity communities on two types of brain surfaces in Fig. 5 (5(a) corresponding to ICBM152 brain surface, and 5(b) corresponding to ICBM152 smoothed brain surface). Clusters of voxels in red colors represent functional areas on the surfaces, in which associated BOLD signals are relatively more synchronous that those in blue colors.

B. Experiments With in Vivo Data
1) Parcellation of Subcortical Nuclei: LTM images of subcortical nuclei are shown in axial, coronal and sagittal sections in Fig. 6(a)-6(m), which are segmented from the whole brain with a coarse subcortical mask. Based on the apparent "valleys and cliffs", functional boundaries can be delineated in the LTM images (e.g., as indicated by red lines in Fig. 6(a)). Voxels with high intensities can be clustered into connectivity communities, whereas voxels with low intensities do not belong to any connectivity communities since BOLD signals in these voxels are relatively asynchronous. As indicated by red arrows in Fig. 6(d), 6(h) and 6(i), many voxels around the core of the hippocampus have low intensities, which implies a potential pitfall in previous studies that all voxels of hippocampus are engaged in regional connectivity analyses. In Figs. 6(h)&6(i), LTM images show a concentrically stratified structure of the hippocampus in red contours, which agrees with an earlier report that was based on tissue morphometry (redisplayed in Fig. 6(m)) [39]. Using the boundary delineations of 3D LTM images, parcellation of the subcortical region can be used to generate an atlas, as shown in Fig. 7(a). This atlas Fig. 6. Typical LTM images of subcortical nuclei shown from axial, coronal and sagittal perspectives (a-i) and anatomical structures of hippocampus (j). To display detailed structures, the LTM images of subcortical nucleus are segmented from the LTM images of whole brain with a coarse subcortex mask. Using the'valleys or cliff' mechanisms, functional boundaries can be delineated in the LTM images. Some typical boundaries are indicated by red lines in a. Note that the centrically stratified structures of hippocampus (indicated in d, h and i) agree well with anatomical structures of hippocampus (j) that is reproduced from [39]; use permitted under the Creative Commons Attribution License CC BY 4.0.
identifies well-known anatomic nuclei, namely, the nucleus accumbens, the globus pallidus, the caudate nucleus, the putamen, the amygdala, the hippocampus, and the thalamus. Note that, in the LTM images, the thalamus is subdivided into the anterior thalamus, the dorsal thalamus and the posterior thalamus. Validation of these parcellations was conducted on a separate dataset of 13 subjects. As shown in Fig. 7(b), statistical analysis demonstrates that the homogeneities of all the LTM-derived parcels of the nuclei are significantly larger than those of geometrically perturbed parcels (P<0.05: the nucleus accumbens, the caudate nucleus, the globus pallidus and the amygdala, the anterior thalamus, the dorsal thalamus, the posterior thalamus and the hippocampus; P<0.01: the putamen).
2) Parcellation of Occipital White Matter: Fig. 8 shows a functional parcellation of occipital white matter. Three typical LTM images are displayed in axial, coronal and sagittal sections in Fig. 8(a)-8(c), respectively, wherein clear boundaries can be observed between the gray matter and the white matter. With delineation of the boundaries, occipital white matter was then parcellated on the LTM images. The results are shown in Fig. 8(d). Fig. 8(e) shows a mask that was generated within the same region, from which geometrically perturbed parcels were derived to test whether these results were significantly different. The parcellation of occipital white matter was also validated using the second dataset of 13 subjects. In Fig. 8(f), the violin plot shows that the homogeneities of the LTMderived parcels are significantly larger than those of geometrically perturbed parcels.

IV. DISCUSSION
The recently emerging field of GSP provides novel, intuitive methodologies for structurally constrained signal processing, which can be built upon a core matrix operator, or graph Laplacian, and its derived concepts [40], [41]. Essentially, a graph Laplacian is a differential operator that characterizes . This atlas recapitulated well-known anatomic nuclei, namely, the nucleus accumbens, the globus pallidus, the caudate and the putamen, the amygdala, the hippocampus, the thalamus. Note that the thalamus is subdivided into the anterior thalamus, the dorsal thalamus and the posterior thalamus. (b), results of homogeneity assessment. In violin plot, red dots indicate homogeneity values of LTM informed parcellation. Bottom and top edges of the narrow box in violin plot indicate 25th and 75th percentiles of the distribution, respectively. The central mark indicates the median. The whiskers extend to the most extreme data points that are not considered outliers (1.5 × interquartile range). Statistical analysis demonstrates that the homogeneities of all the LTMderived parcels of the nuclei are significantly larger than those of geometrically perturbed parcels (P<0.05: the nucleus accumbens, the caudate nucleus, the globus pallidus and the amygdala, the anterior thalamus, the dorsal thalamus, the posterior thalamus and the hippocampu; P<0.01: the putamen).
global discontinuousness (or smoothness) of signals on a graph, where the measure of discontinuousness of signals on the vertices are dependent on connection weights of edges contained in the graph. This allows it to be exploited for capturing spatial patterns of connectivity (e.g., the connectivity gradient) [11], [12], [16] and structure guided spatial smoothing (e.g., diffusion-informed spatial smoothing of fMRI data in white matter) [14], [42]. Another typical application with clear physical pictures is the use of a Laplacian to describe the diffusion of virtual quantities/signals on discrete graphic spaces, where the gradient of the quantities -the 'driving force' of diffusion -can be described by this operator [43], and this concept has recently been used to define spatial scales of brain functional networks [44]. The distinct contribution of this work is the proposal of a novel operator based on a GSP framework, wherein virtual quantities/signals propagate on the graph at a rate modulated by localized connectivity that simultaneously contains local connectivity and global similarity. With the representation of the virtual signals in the anatomical space, this new operator can capture spatial patterns of connectivity in 3D spaces. Methodologically, our study debuts a new paradigm to construct graph operators with intuitive physical pictures for characterizing and visualizing the features that implicitly exist in the structure of the graph.
In our approach, two novel mathematical techniques are formulated to characterize the voxel-wise connectivity structures using fMRI data. First, an SVD-based filtering is used to suppress noise and signal components that are shared by multiple connectivity communities, thereby enhancing the boundary features. As demonstrated by the connectivity matrices in our simulation studies (Fig. 3(b)&3(c)), background connections between two voxel communities, that likely stem from common components shared by the two communities and tend to blur their connectivity boundaries, are suppressed to low levels by this SVD-informed filtering. Second, a new LTM reconstruction is proposed to visualize connectivity structures that are implicitly encoded in a localized topo-connectivity graph. The resulting LTM images reveal spatial connectivity patterns with the boundaries defined by local connectivity and global similarity. With proof-of-concept experiments on simulated that provide intuitive interpretations of the proposed approach, we conducted in vivo validations using 7T fMRI data sourced from HCP.
In our in vivo validation, the functional structures delineated by the boundaries in LTM images were compared with anatomical structures conveyed by T 1 w images and DTI FA colored maps. Overall, most of the functional structures revealed by LTM images show similar boundaries with the relevant structures in the two anatomical images, which agree with the neuroanatomic knowledge that the brain function is dependent on its structure [18], [20], [35]. In particular, the LTM images reveal some subtle structures that are not observed in the anatomical images, which demonstrates that LTM images can offer distinctive information in functional structures that are not implied in anatomical structures. Subtle connectivity structures delineated by group-average LTM images reconstructed from the fMRI dataset in HCP provide a basis for parcellation of a functional atlas of the human brain. The LTM contrast not only reveals functional structures, but also measures the synchrony of fMRI signals associated the voxels within a neighborhood. Voxels with low intensities in LTM images do not belong to any specific functional structures. As demonstrated in Fig. 4, a lot of voxels possess low intensities in the LTM images. However, the recently reported method based on functional connectivity gradients, which successfully parcellates subcortical structures, may not treat such voxels optimally [10]. For example, we considered parcellation of the hippocampus. As shown in Fig. 5, numerous voxels around the core of the hippocampus have low intensities, which agrees with findings that the hippocampus has a centric stratified structure [39], [45]. In this study, the voxels of low contrast in the LTM images are excluded to ensure that voxels contained in the parcellation atlas have synchronous BOLD signals.
Our first demonstration of LTM reconstruction is the parcellation of subcortical structures, which successfully recapitulates the well-known anatomic nuclei ( Fig. 7(a)). Notably, the statistical analysis we performed to determine whether the synchrony of fMRI signals in LTM-derived parcels are larger than geometrically perturbed parcels is a further development of the existing statistical approach [10], [21]. Compared with the approach of random parcellation under relatively loose constraints [10], [21], our statistical analysis have two important improvements: (1), voxels constituting geometrically perturbed parcels are confined in a mask derived from the dilation of the LTM-derived parcels with three layers of voxels and (2), most voxels in the perturbed parcels are required to be contained in the responding LTM-informed parcellation. The two constraints could reduce false negatives and false positives in statistical analysis. Hence, statistical analysis based on the parcellations with geometrically perturbations should be more valid, suggesting that our statistical inference on parcellations ( Fig. 7(b)) is more valid.
The second demonstration of LTM reconstruction is the parcellation of occipital white matter, which provides insights to the emerging study of BOLD effects in white matter [46], [47]. Clear functional boundaries in LTM images can be observed between gray matter and white matter, which suggests that BOLD signals of occipital white matter have independent origins. Furthermore, clear boundaries between gray matter and white matter can also be observed in the whole brain in LTM images. Based on this finding, we can exclude the possibility that BOLD signals in the white matter originate from effects of veins draining from the gray matter. The drainage effects from larger vessels in gray matter impair the spatial specificity of BOLD contrast in fMRI [48], [49], which have been proposed as potential sources of white matter fMRI signals [50].
LTM images of human brain demonstrate that functional structures also exist in various regions of white matter, which can provide two important principles for parcellation of white matter in a future study. First, voxels with low contrast intensities in LTM images should be removed from an atlas of white matter since BOLD signals associated with these voxels are not synchronous. Second, fine parcellations of white matter can be conducted using a combination of LTM images and diffusion tensor images. Improving the accuracy of atlas parcellation is essential to the region-level connectivity analysis of white matter since BOLD signals of white matter have low SNR and hence their detection is sensitive to small errors.
To improve the performance of the filtered LTM reconstruction in discerning subtle functional structures in the brain, some technical limitations in this work may be addressed. First, to enhance the boundary features of LTM reconstructed functional structures, the connection of two vertices with a small global similarity was pruned with an empirical threshold such that the graph stays connected with minimum edges. One potential strategy to address this limitation is to reconstruct LTM maps iteratively, in which a new reconstruction is implemented on a volume defined by a functional structure from the previous iteration, thus providing a better capability of fine parcellations of the functional structure with intricate organization. A second limitation of this work is that, in the pruning process, the global similarity of two voxels is measured using the Euclidean distance of the corresponding columns in the matrix S. There is an attempt using a non-Euclidean distance (η 2 coefficient) for this purpose, which exhibits an extraordinary ability to capture functional boundaries encoded in global similarity [10], [11] and may promise an LTM image with finer functional structures. In this study, large numbers of voxels could be involved in the 3D LTM reconstruction so that using the non-Euclidean was not feasible due to the consequential high computational cost. The non-Euclidean distance method, along with the iterative reconstruction scheme mentioned above, could be implemented in the future by designing efficient algorithms that serve these purposes.

V. CONCLUSION
The proposed LTM reconstruction provides a novel method to derive voxel-wise functional structures within the whole brain that includes the cortex, subcortex and white matter. Our statistical tests demonstrate that the functional parcels derived from such LTM reconstructions exhibit highest regional homogeneities compared with the parcels with small geometric perturbations. In the future, the LTM reconstruction can be used to parcellate functional structures in other brain regions of the brain, which could offer a basis for whole brain connectivity analyses at regional levels.