Constructing Connectome Atlas by Graph Laplacian Learning

1.1. Related Works

Currently, most of the work on brain connectome atlas focuses on the construction of nodes in the brain network [21–23], instead of discovering the most representative node-to-node connectivities in the population. For example, a connectivity-based parcellation framework is proposed to better characterize functional connectivity by using fine-grained network nodes, thus resulting in a node-wise connectome atlas consisting of 210 cortical and 36 subcortical regions [21], Wig et al. [23] applied snowball sampling on the resting-state functional connectivity data to identify the centers of cortical areas, the subdivision of subcortical nuclei, and the cerebellum in a population. However, region-to-region brain connectivity in the population remains unexplored.

Among very few connectome atlas works that address the embedding of the most representative characteristics of brain connectivities in the population, the most popular approach is to average the individual connectivity values at each link separately by treating the brain network as a data matrix [24]. Albeit simple, such averaging-based solution might break down the network topology since the brain network is a structured data representation with complicated data geometry. It is worth noting that the similarity network fusion technique was used in [25] to propagate the connectivity information by random walk method from one subject to another until all networks become similar [26], Thus, the connectome atlas is the average over the diffused matrices, instead of the original network matrices. However, matrix diffusion suffers from the issue of network degeneration where each node has an equal probability of connecting with all other nodes after applying the excessive amount of matrix diffusions.

1.2. Our Contributions

To address the above challenges, we propose a first-time learning-based graph inference framework to unravel the connectome atlas in the population. To do so, we first introduce diffusion distance to measure the network connectivity, a time-dependent connectivity measurement, reflecting both local and global network geometry at different scales. Next, we employ a diffusion mapping technique [27] to embed the brain network into the Euclidean space, wherein each node of the network is represented as a vector in the Euclidean space, and the Euclidean distance between two vectors corresponds to the diffusion distance between two nodes they represent. After that, we propose a learning-based approach to discover the shared network that is prevalent in the population. Specifically, we regard the diffusion map vectors of each brain network as the graph signals [28] that reside on the latent common network. Thus, the graph-based inference is optimized to learn a graph shift operator such that the topology of the common network has the largest consensus with all aligned graph signals. Since the diffusion distance can integrate local connectivity to the global geometry of the whole network, our learning-based framework offers a novel solution for constructing multi-scale connectome atlas that can “think globally, fit locally” [28, 29],

We evaluate the accuracy and robustness of our learning-based connectome atlas construction method on both functional and structural network data in terms of various group comparisons including Autism Spectrum Disorders (ASD), Fronto-Temporal Dementia (FTD), and Alzheimer’s Disease (AD). Compared to the non-learning based counterparts (i.e., a network averaging method [24] and network diffusion method [25]), our method shows more reasonable atlas construction results.

2.1. Definition of Graph and Graph Laplacian

Suppose each brain network consists of N nodes V ={v_i|i = 1,…, N}. Affinity matrix W =[w_ij]_i,j=1,…,N is often used to encode brain network in the form of a matrix, where each element w_ij measures the strength of connectivity between nodes v_i and v_j. In the structural network, w_ij is related to the number of fibers traveling between two regions [4], On the other hand, in the functional network, w_ij measures the synchronization of functional activities in terms of statistical correlation between BOLD (blood oxygen level dependent) signals [30, 31]. For convenience, it is reasonable to assume that each adjacent matrix is symmetric (i.e., w_ij = w_ji), non-negative ( w_ij ≥ 0), and no self-connection (w_ij =0). To that end, each brain network can be represented by a graph structure $\mathbb{G}=(\mathit{V},\mathit{W})$ The corresponding Laplacian matrix L is also a N × N matrix, i.e., L = D – W where D is the degree matrix that contains the degrees of the vertices along the diagonal with the i^th diagonal element equals ${\Sigma }_{j=1}^N{w}_{ij}$.

Since the graph Laplacian L is symmetric, it has a complete set of orthonormal eigenvectors, where we use $\mathit{X}={\left[{\chi }_c\right]}_{c=1}^N$ denote the left eigenvector matrix of the Laplacian matrix L Each χ_c is a column vector corresponding the c^th eigenvalue of L. Each χ_c is associated to the eigenvalue ω_c satisfying Lχ_c = ω_cχ_c. Since we consider connected graphs here, we assume the eigenvalues of graph Laplacian matrix L are ordered as 0 = ω₁ < ω₂ ≤ ω₃ …≤ ω_N. We denote the entire graph spectrum by Ω(L) = {ω₁, ω₂,…ω_N}.

2.2. A Graph Fourier Transform and Notion of Frequency

For any signal f ∈ ℝ^N, the classic Fourier transform $\widehat{f}\left(\xi \right)={\int }_{\mathbb{R}}f\left(t\right)e^{-2\pi i\xi t}dt$ is the expansion of a function f in terms of complex exponentials, which are the eigenfunctions of the one-dimensional Laplacian operator −Δ(e^−2πiξt)=(2πξ)²e^2πiξt. In analogy, the graph Fourier transform $\widehat{f}$ of signal f on the vertices of $\mathbb{G}$ can be defined as the expansion of f in terms of the eigenvectors of the graph Laplacian L [28]: $\widehat{f}\left({\omega }_c\right)={\Sigma }_{i=1}^{N}f\left(i\right){\chi }_c^{\ast }\left(i\right)$, where ${\chi }_c^{\ast }$ is the column vector in the right eigenvector matrix of L. The inverse graph Fourier graph transform is given by $f\left(i\right)={\Sigma }_{c=1}^N\widehat{f}\left({\omega }_c\right){\chi }_c\left(i\right)$.

The eigenvalues {(2πξ)²} in the classic Fourier signal analysis carries a specific meaning of frequency. Specifically, the complex exponential eigenfunction associated with smaller ξ exhibits slower oscillating patterns, while higher ξ is often associated with complex exponential eigenfunction with much more rapid oscillations. In the setting of a graph, each element in χ₁ (associated with ω₁ = 0) is constant and equals to $\frac{1}{\sqrt{N}}$ at each graph vertex. As pointed in [28], the eigenvector χ_c with smaller eigenvalue ω_c (lower graph frequency) varies slowly across graph in term of that χ_c(i) and χ_c(j) are more likely having similar values if vertices V_i and V_j are connected in the graph. On the other hand, the eigenvectors associated with larger eigenvalues (higher graph frequency) oscillated more rapidly and more likely to have distinct values on the vertices connected in the graph even with high connectivity degree.

2.3. Diffusion Distance and Diffusion Mapping

2.4. Smooth Graph Signals

Graph signal processing [28] is an emerging research area for analyzing the structured data, where the signal values are defined on the vertex of a weighted graph. Taking the brain network as an example, we assume to have a mean cortical thickness for each node in the network. It is common to arrange the centralized cortical thickness degrees into a data array by following the order of brain parcellation, as shown in Fig. 3(a), where red and green denote for positive (higher than the mean) and negative (lower than the mean) values respectively. Note that the data array alone is unable to characterize the relationship of cortical thickness between two regions. In the setting of a graph, however, it is natural to regard the set of whole-brain cortical thickness as a graph signal since the interplay of cortical thickness between two brain regions can be interpreted in the context of the associated network connectivity.

Without a doubt, we can map the same data array to multiple graphs and potentially lead to different graph signals like the examples of two graph signals shown in Fig. 3(b)–(c). Although these two graph signals are both valid, the graph signal shown in Fig. 3(b) makes more sense since the signal is smooth in terms that two connected nodes have similar scalar values along with the connection defined by the underlying network. Given the network $\mathbb{G}=\left(\mathit{V},\mathit{W}\right)$, the smoothness of graph signal f can be measured in terms of a quadratic form of the graph Laplacian:

where f(i) and f(j) are the signal values associated to node v_i and v_j, respectively. The intuition behind Eq. (2) is that the graph signal f is considered to be smooth if the strongly connected nodes (with a large connectivity strength) generally have similar values. The smaller the quadric term in Eq. (3), the smoother the graph signal f on the network $\mathbb{G}$. As we will explain in the next section, we lavage this smooth graph model to find the latent common network from a set of individual brain networks.

2.5. Topological Distance Between Brain Networks

Since brain network is usually encoded in an adjacency matrix, many existing brain network distances are defined based on matrix norm that uses the sum of element-wise difference. For example, the L_∞ distance between two network χ¹ and χ² is defined as $D_{\infty }\left({\chi }^1,{\chi }^2\right)={\max }_{\forall i,j}\left|d_{i,j}^1-d_{i,j}^2\right|$, where $d_{i,j}^1$ and $d_{i,j}^2$ measure the correlation distance d_ij = 1 – ρ_ij where ρ_ij is the correlation between two regions. Such distance measurements have limited capability to quantify the topological differences such as connected components and modules in brain networks [32], Gromov-Hausdorff (GH) distance is proposed in [33, 34] to measure the topological difference two brain networks. Specifically, the shape of network is first defined using the graph filtration technique that iteratively build a nested subgraphs of the original network. As shown in the box of Fig. 4, we start the distance threshold ε from the smallest correlation distance and gradually increase ε until ε reaches the largest value of correlation distance. We connect two nodes x_i and x_j if their distance is less than current ε (i.e., d_ij < ε). If two nodes are already connected directly or indirectly via other intermediate nodes by smaller ε, we do not connect them. For example, we do not connect x₂ and x₅ in Fig. 4 at ε = 0.32 since we are already connected through other nodes at ε = 0.3. As ε increases, we can obtain a sequence of nested graphs which is called graph fdtration in algebraic topology [35]. After that, we construct a single linkage matrix S where each element quantifies the single linkage distance between nodes x_i and x_j in the network. Note, the single linkage distance s_ij is the minimum ε that makes nodes x_i and x_j are connected on the dengrogram (shown in the bottom of Fig. 4). For example, the single linkage distance between x₁ and x₆ in Fig. 4 is 0.3, although there is a direct link between x₁ and x₆ (the distance is 0.4 on the edge). GH distance between two networks X¹ and X² is then defined through the corresponding single linkage matrix L¹ and L² as $D_{GH}\left(x^1,x^2\right)={\max }_{\forall i,j}\left|l_{i,j}^1-l_{i,j}^2\right|$. It has been demonstrated in [33] that GH distance, characterizing the network topology, has superior performance over other graph theory based network similarity measures.

3.1. Overcome Irregular Data Geometry: from Network to Graph Signals

First, we address the challenge of irregular data structure by embedding each individual brain network ${\mathbb{G}}_m$ into its own native diffusion space. Specifically, we obtain a set of diffusion mappings {Ψ^t,m} by applying eigen decomposition on the adjacency matrix W_m and constructing the diffusion mapping ${\mathit{\Psi }}^{t,m}={\left[{\psi }_c^{t,m}\right]}_{c=1}^N$ at each scale t. As shown in Fig. 2I, a larger scale t renders more global connectivity characteristics in the entire network. In order to maintain the details of connectivity profiles in the common network, we propose to learn the common network $\overline{\mathbb{G}}$ at a relatively local scale (i.e., t = 1) based on the set of diffusion coordinates $\left\{{\mathit{\Psi }}^{t=1,m}|m=1,\mathrm{...},M\right\}$. After that, it is straightforward to construct a multi-scale connectome atlas based on the learned adjacency matrix $\overline{\mathit{W}}$. For simplicity, we omit the superscript t to the variables related to scale t in the following text. Note, there are alternative approaches that can be used here to transform a graph (network) into a collection of signals such as multi-dimensional scaling [36] used in [37]. The reason we use the diffusion mapping technique here is mainly because the diffusion distance is a function of scale t, which allows exploring the multi-scale connectome atlas in a local to global manner.

Fig. 5(a)–(c) demonstrate the procedure for embedding each network data to graph signals. For each individual network ${\mathbb{G}}_m$, its diffusion coordinates ${\mathit{\Psi }}^m={\left[{\psi }_c^m\right]}_{c=1,\mathrm{...},N}$ consist of N column eigenvectors, where each element in ${\psi }_c^m$ is associated with the corresponding node in the network. Instead of considering each ${\psi }_c^m$ as a column vector of diffusion coordinate, we further regard each ${\psi }_c^m$ as the graph signal residing on the latent common network $\overline{\mathbb{G}}$ (shown in Fig. 5(c)), where the element-to-element relationship in each ${\psi }_c^m$ is interpreted in the context of $\overline{\mathbb{G}}$ (shown in Fig. 5(d)). In addition, the graph signal ${\psi }_c^m$, which is structured data representation, allows us to quantify the representativeness of common network $\overline{\mathbb{G}}$ by inspecting the consistency between the diffusion coordinates {Ψ^m} that derive from individual networks and the common network $\overline{\mathbb{G}}$ that supports these diffusion coordinates.

3.2. Learn Graph Laplacian of Common Network: from Graph Signals Back to Network

Since the graph signals $\left\{{\psi }_c^m\right\}$ are defined in Euclidean space, the naïve approach to construct a common network consists of two steps: (1) calculate the mean graph signal y_c at each oscillation mode by arithmetic averaging, i.e., $y_c=\frac{1}{M}{\Sigma }_{m=1}^M{\psi }_c^m$, and (2) reconstruct the adjacency matrix $\overline{\mathit{W}}$ where each element $w_{ij}=1/{\Sigma }_{c=1}^N{\Vert y_c\left(i\right)-y_c\left(j\right)\Vert }^2$ (the inverse of diffusion distance in Eq. (2)) for i ≠ j. w_ii = 0 since we do not allow for self-connection.

However, such averaging-based approach has two limitations. (1) Since the graph signals associated with c^th largest eigenvalue might be associated with different eigenvalues across M individual networks, direct averaging over $\left\{{\psi }_c^m|m=1,\mathrm{...},M\right\}$ in each oscillation mode might introduce the spurious connectivity information that is not related to the connectome atlas. In analogy to the scenario that one cubic box undergoes various rotations and translations in the 3D Cartesian coordinates (the purple boxes in Fig. 5(b)), the success to unravel the accurate atlas shape (i.e., cubic box) is dispensable to the removal of the external variations (i.e., translations and rotations) which is not related to the intrinsic nature of the underlying shape. To that end, it is of necessity to align each graph signal ${\psi }_c^m$ in the common space, such that the graph signals at the same oscillation mode are associated with the same eigenvalue. (2) Simply averaging individual diffusion maps ignores the interplay between data and the latent common network. As shown in Fig. 3, a data array consisting of unordered scalar values can potentially reside on different graphs and then lead to entirely different graph signals. In this regard, the common network $\mathbb{G}$ should be consistent to all graph signals at all network frequencies. To tackle this issue, we propose a learning-based approach to learn the Laplacian matrix of common network $\overline{\mathbb{G}}$ from a set of graph signals $\left\{{\psi }_c^m||c=1,\mathrm{...},N,m=1,\mathrm{...},M\right\}$ which mainly consists of two alternative steps (shown in Fig. 5(c)–(e)).

3.3. Optimization

Although the energy function in Eq. 4 is non-convex, we opt to optimize the representative graph signals {y_c}, weights ${\pi }_c^m$ and the common graph Laplacian L in an alternative manner.

Description of dataset.

There are four datasets used in the following experiments in total. (1) ASD dataset. We use the resting-state fMRI (rs-fMRI) data of 45 ASD patients and 47 healthy controls (HC) from the NYU Langone Medical Center of Autism Brain Imaging Data Exchange (ABIDE) database (http://preprocessed-connectomes-proiect.org/abide/). (2) FTD dataset. We processed the rs-fMRI data of 90 FTD subjects and 101 age-matched HC subjects from the recent NIFD database, which is managed by the frontotemporal lobar degeneration neuroimaging initiative (https://cind.ucsf.edu/research/grants/frontotemporal-lobar-degeneration-neuroimaging-initiative-0). Furthermore, there are three major clinic subtypes in the FTD cohort: behavioral variant of FTD (BV), semantic variant of FTD (SV), and progressive nonfluent aphasia (PNFA). For the above functional neuroimaging data, we follow the partitions of AAF template [41] and construct the function network for each subject with 90 nodes ¹. The region-to-region connection is measured by Pearson’s correlation coefficient between the mean time courses. (3) ADNI dataset. In total, 51 subjects (31 HC and 20 AD) are selected from the ADNI database (http://adni.loniusc.edu/). For each scan, we parcellate the cortical surface into 148 regions [42] using FreeSurfer on T1-weighted MRI scans. Then we apply the probabilistic fiber tractography on DWI (difiusion weighted imaging) and T1-weighted image using FSF software library to obtain a 148×148 connectivity matrices of the structural networks, where each element reflects the number of fibers traveling between two brain regions. (4) HCP dataset. Resting stage fMRI data from 90 healthy young adults from Human Connectome Project (HCP) are used in the replicability test in Section 4.2. MRI scanning was done using a customized 3T Siemens Connectome Skyra with a standard 32-channel Siemens receive head coil and a body transmission coil. Resting state fMRI data were collected using a gradient-echo echo-planner imaging (EPI) with 2.0mm isotropic resolution (FOV=208×180mm, matrix=104×90, 72 slices, TR=720ms, TE=33.1ms, FA=52°, multi-band factor=8, 1200 frames, ~15min/run). Runs with left-right and right-left phase encoding were done to correct for EPI distortions. Here, we consider the left-right and right-left runs as test/retest data. For each fMRI data, we first partition whole brain into 268 regions using the spatially aligned T1 image. Then, the function connectivity network is constructed based on the mean BOLD signals of each parcellated region. The demographic data of these three datasets is summarized in Table 1.

Experiment setup.

As shown in Eq. 4, α, β, and ρ are three critical parameters required in the optimization. Due to the lack of ground truth, we determine the optimal set of parameters based on the simulated network data in Section 4.1. Then, we evaluate the replicability of connectome atlas on the ASD dataset in Section 4.2. The quality of connectome atlas is inspected visually and quantitatively in Section 4.3 for ASD, FTD, and AD population separately. Finally, we demonstrate the augmented statistical power of network analyses using our learning-based connectome atlas construction method. For comparison, the counterpart methods include (1) matrix averaging on original networks (called “matrix averaging”), and (2) graph-based averaging with matrix diffusion [25] (called “graph diffusion”). Note that none of the counterpart methods is learning-based.

4.1. Validation on Simulated Network Data

In this experiment, we first generate a set of networks from a predefine node distribution where the node locations and affiliations (clusters membership) are fixed. However, we randomly mislabel the node affiliations in each network instance. Since the mislabeling process is performed in a random manner, the average of these individual networks should respect the original node distribution, which becomes the criterion of measuring the accuracy of connectome atlas construction methods. As shown in Fig. 6, a set of simulated networks were generated from the predefined node distribution (Fig. 6(a)), with three separable clusters (displayed in red, blue, and purple). At each iteration, we selectively swap several data points (randomly, near the cluster boundaries) to another cluster. For the data points with correct cluster labels, we construct the connection based on the distance between two points. However, we specifically leave some points not sampled, resulting in the isolated nodes in the network. For the data points with the incorrect cluster labels, we only allow very few connections to the points with the same label. We show three typical simulated networks in Fig. 6(b).

We simulate networks with different proportions of node swapping percentage, each with 50 simulations. Given the swapping percentages, we apply our learning-based connectome atlas construction method to the simulated networks with different sets of parameters. Specifically, we apply a grid search strategy for α, β, and ρ separately, where the search range is set to [0.005,2.0] for each parameter, and the search step is set to 0.005. For each setting, we apply the classic k-mean clustering method to determine three clusters based on the common network. The criterion of parameter selection is the dice ratio between the ground truth and the clustering result based on common network $\overline{\mathit{W}}$. In this way, we fix using α = 0.5, β = 1.0, and ρ = 0.5 in the following experiments.

Furthermore, we apply matrix averaging (no parameter needed) and graph diffusion (using the recommended parameters in [25]) methods on the same simulated network with respect to different swap percentages. We show the dice ratio between the ground truth (Fig. 6(a)) and the clustering result based on the common networks in Table 2. Since the clustering result closely depends on the node-to-node connectivity encoded in the common network, it is reasonable to consider that our graph learning-based connectome atlas construction method achieves higher accuracy than the counterpart non-leaming-based methods.

4.2. Evaluation of Functional Connectome Atlas

4.3. Evaluation of Structural Connectome Atlas

Discussions.

As shown in Fig. 2 and Fig. 10, our connectome atlas construction offers the way to construct multi-scale atlas by adjusting the scales in diflusion mappings. It is worth noting that our proposed method constructs the connectome atlas at a given scale t separately. Although it is always favorable to jointly model the dependency of connectome atlas across scales and simultaneously construct the multi-scale connectome atlas, it is computationally prohibitive to do so in our current optimization framework as illustrated in Section 3.3. However, there might be several possible solutions to explore the new regime of multi-scale connectome atlas. (1) Unify multi-scale connectome atlases by fusing the connectivity information from local to global manner. First, we generate a set of connectome atlas with different setting of scale t. Then, we can apply the matrix fusion method in [50] to combine into single network which includes both local and global connectivity information. (2) Use supervised/unsupervised feature selection method to find the multi-scale connectome feature representation. Since most of the network-based studies is interested in finding network feature representations which can separate two or multiple cohorts, we can extend the current approach into the scenario of multi-scale by using feature selection technique such as sparse feature selection [51, 52] to determine the best feature representations from multi-scale connectome atlases.

Conclusions.

The construction of connectome atlas is much more challenging than image atlas, due to the irregular graph structure. To address this problem, we proposed a learning-based approach to discover the common network in the setting of graphs. Our proposed method offers a new window to investigate the population-wise whole-brain connectivity profiles in the human brain. Promising results have been obtained in various brain network studies such as ASD, FTD, and AD, which indicates the applicability of our proposed learning-based connectome atlas construction methods in neuroimaging applications.