1 Introduction

Connectomes can be described as a graph of organized regions and their connections that putatively have foundational roles in emerging functional and cognitive outcomes [1]. Hence, many analyses in cognition, learning, and brain diseases and disorders investigate the organization of the brain [2]. Graph theoretical approaches such as complex network analysis provide powerful tools to study structural and functional characteristics of the brain without losing its organizational features [3].

In traditional connectomes, when representing the brain as a network, the nodes of the network correspond to the brain regions, and the edges between the nodes correspond to connections between those regions. In this approach, networks encode only inter-regional connectivity. The regional information such as degree, strength, or centrality that may be crucial to the analysis are usually treated as confounding factors or covariates. This hinders interpretations regarding regional changes due to, for instance, an underlying pathology. However, graph theory facilitates a principled methodology to combine regional characteristics (node features) with interactions between regions (edge features), by means of annotating nodes of the network [4]. Hence, the first contribution of this study is to provide a rich brain network representation, an enriched connectome, that enables inclusion of such nodal features when modeling brain connectivity.

Such a rich representation of brain organization including nodal features requires a new set of tools such as a similarity measure between these networks (graphs) which is essential for classification, clustering, or regression tasks [5, 6]. As a second contribution, we propose a graph matching algorithm that provides a similarity measure between brain networks with nodal features. Among several approaches proposed in the literature to calculate graph similarity over brain data such as seeded graph matching [7] and graph embedding [8], graph edit distance (GED) is arguably the most effective method with promising results [9, 10]. However, high running time complexity of GED requires use of approximation techniques such as Hungarian algorithm in [11] and hinders a detailed analysis of edge features [12]. We approach the graph matching problem as an instance of the metric labeling problem [13] and provide an efficient approximation algorithm using the primal-dual scheme [14] by extending our previous study [15]. Our graph matching method achieves two goals simultaneously: finding a mapping between brain regions of different graphs and computing a similarity score. The enriched connectome along with the graph-based similarity measure facilitates its use in classification of samples and we demonstrate its effective application on a traumatic brain injury (TBI) dataset. Results show that our enriched connectome along with the proposed matching algorithm provides better classification between the groups than the traditional connectivity based connectome representation.

2 Materials and Method

2.1 Dataset

Participants: We use a traumatic brain injury dataset consisting 39 patients (12 female) with moderate-to-severe TBI examined at 3 months post injury and 30 healthy controls (8 female). Age of patients are in [18, 65] years with a mean of 35 years and standard deviation of 14.7 years, while the age of healthy controls are in [20, 56] years with a mean and standard deviation of 34.7 and 9.9 years, respectively. Duration of post-traumatic amnesia of patients, which can be considered as a measure of trauma severity, has a mean of 26.7 days with a standard deviation of 21.2 days.

Data Acquisition and Preprocessing: For each subject, DTI data was acquired on a Siemens 3T TrioTim scanner with a 8 channel head coil (single shot spin echo sequence, \(\text {TR/TE} = 6500/84\,\text {ms}, b = 1000\,\text {s/mm}^{2}, 30\) gradient directions). 86 region of interests from the Desikan atlas [16] were extracted to represent the nodes of the structural network. A mask was defined using voxels with an FA of at least 0.1 for each subject. Deterministic tractography was performed to generate and select 1 million streamlines, seeded randomly within the mask. Angle curvature threshold of 60 degrees, and a min and max length threshold of 5 mm and 400 mm were applied, resulting in an \(86\times 86\) adjacency matrix of weighted connectivity values, where each element represents the number of streamlines between regions.

2.2 Enriched Connectome

Given parcellation of the brain into 86 regions, we constructed a weighted undirected graph with 86 nodes corresponding to brain regions and weighted edges corresponding to the number of fibers connecting region pairs. We annotate each node with two set of features. First, we generated a 6 dimensional feature vector by calculating various graph theoretical measures for each node, namely degree, strength, betweenness centrality, local efficiency, participation coefficient, and local assortativity, using the Brain Connectivity Toolbox [17]. While calculating participation coefficient of nodes, we used association of structural regions with 7 functional systems as described in [18]. Second, we generated an 86 dimensional feature vector, representing the weighted connectivity of each node to the rest of the nodes in the graph, where we considered self edges to be nil. In summary, our graph representation, denoted enriched connectome hereby, incorporates graph theoretical measures of the connectome along with the weighted connectivity of the regions that are to be found in network representations. We normalized the values of each graph theory measure and the edge weights to [0, 1] in order to make them comparable across subjects.

2.3 Graph Matching Based Similarity Measure

We propose taking a graph matching approach to define a similarity measure between two enriched connectomes, while providing a mapping between their nodes. We note that since the brains are parcellated into a common atlas in our setup, mapping between the regions are known a priori. However, we expect to get several mismatching nodes between dissimilar enriched connectomes due to differences in the connectivity of the network, making the similarity of graphs and the ratio of mismatching nodes effective measures for identifying patients from controls.

To this end, we evaluate the graph matching as a special case of the metric labeling problem [13]. Translating the metric labeling into the domain of brain graphs, the problem reads as follows: Given two enriched connectomes A and B, find the optimal one-to-one mapping \(f:\mathcal {A}\rightarrow \mathcal {B}\) between their nodes while minimizing the following objective function:

$$\begin{aligned} \beta \displaystyle \sum _{a \in \mathcal {A}} c(a,f(a)) + (1-\beta )\sum _{a,a' \in \mathcal {A}} w(a,a')\cdot d(f(a),f(a')). \end{aligned}$$
(1)

The first summation term in (1) is regarded as the assignment cost with \(c:\mathcal {A}\times \mathcal {B}\rightarrow \mathbb {R}\) that determines the cost of mapping a brain region \(a \in \mathcal {A}\) to a region \(b \in \mathcal {B}\), which we define as \(||v1_a - v1_b||_2 + ||v2_a-v2_b||_2\) where v1 and v2 indicate the graph theoretical and edge weight based feature vectors described earlier, respectively. The second summation term stands for the separation cost, penalizing strongly connected brain regions \(a, a' \in \mathcal {A}\) in getting mapped to loosely connected regions \(b,b' \in \mathcal {B}\) with \(w:\mathcal {A}\times \mathcal {A}\rightarrow \mathbb {R}\) indicating edge weights in \(\mathcal {A}\) as a measure of connectivity strength and \(d:\mathcal {B}\times \mathcal {B}\rightarrow \mathbb {R}\) indicating the distance between nodes of \(\mathcal {B}\) as a measure of proximity between regions which is defined inversely proportional to the w in \(\mathcal {B}\). Thus, the first half of the cost function encourages mapping each node of \(\mathcal {A}\) to a node that resembles it most in \(\mathcal {B}\) while the second half discourages two strongly connected regions in \(\mathcal {A}\) getting mapped to two loosely connected regions in \(\mathcal {B}\). The variable \(\beta \) in (1) is a balancing term to adjust the contribution of the assignment and separation costs to the objective function which takes values in [0, 1]. Once optimized, summation of the costs in (1) defines a similarity score between the two graphs.

In their seminal paper, Kleinberg and Tardos presented the following quadratic optimization formulation for the metric labeling problem which they showed to be computationally intractable to solve [13]:

$$\begin{aligned} \begin{array}{lll} \min \quad &{}\displaystyle \sum _{\begin{array}{c} a \in \mathcal {A}\\ b \in \mathcal {B} \end{array}} c(a,b)\cdot x_{a,b} + \sum _{\begin{array}{c} a,a' \in \mathcal {A}\\ b,b' \in \mathcal {B} \end{array}} &{} w(a,a') \cdot d(b,b') \cdot x_{a,b} \cdot x_{a',b'} \\ \text {s.t.} \quad &{}\textstyle \sum _{b \in \mathcal {B}} x_{a,b} = 1, &{} \forall a \in \mathcal {A}\\ &{}\textstyle \sum _{a \in \mathcal {A}} x_{a,b} = 1, &{} \forall b \in \mathcal {B}\\ &{} x_{a,b} \in \{0,1\} , &{} \forall a \in \mathcal {A}, b \in \mathcal {B}\end{array} \end{aligned}$$
(2)

where \(x_{a,b}\) is an indicator variable taking value 1 if a is mapped to b and 0 otherwise. They also presented a linear programming formulation of the problem along with an approximation algorithm using hierarchically well-separated trees (HST), which was inefficient due to the computational time it takes to build the HSTs and to solve the linear program. Using another integer linear programming formulation of the problem along with a primal-dual approximation scheme [14], we recently presented an efficient approximation algorithm for the traditional metric labeling problem [15]. Here, we extend the latter study by altering the constraints of the metric labeling to account for the particular case of matching the enriched connectomes. Traditional metric labeling formulation allows many-to-one matching of the nodes between graphs, that is, several nodes of the first graph can be mapped to the same node in the second graph. In the setup of enriched connectomes where the brains are registered to a common atlas and parcellated into the same number of regions across subjects, it is known a priori that there should be a one-to-one mapping between the nodes of the graphs. Motivated by this observation, we impose additional constraints to the metric labeling formulation to enforce a one-to-one mapping between graphs. Our extended version of the metric labeling with the integer linear programming formulation is as follows:

$$\begin{aligned} \begin{array}{lll} \text {min} \quad &{}\displaystyle \sum _{\begin{array}{c} a \in \mathcal {A}\\ b \in \mathcal {B} \end{array}} c(a,b)\cdot x_{a,b} + \sum _{\begin{array}{c} a,a' \in \mathcal {A}\\ b,b' \in \mathcal {B} \end{array}} w(a,a') \cdot &{} d(b,b') \cdot x_{a,b,a',b'}\\ \text {s.t.} \quad &{}\textstyle \sum _{b \in \mathcal {A}}x_{a,b} = 1, &{} \forall a\in \mathcal {A}\\ &{}\textstyle \sum _{a \in \mathcal {B}}x_{b,a} = 1, &{} \forall b\in \mathcal {B}\\ &{}\textstyle \sum _{a' \in \mathcal {A}}x_{a,b,a',b'} = x_{a,b}, &{} \forall a \in \mathcal {A},b,b' \in \mathcal {B}\\ &{}\textstyle \sum _{b' \in \mathcal {B}}x_{a,b,a',b'} = x_{a,b}, &{} \forall a,a' \in \mathcal {A},b' \in \mathcal {B}\\ &{}x_{a,b,a',b'} = x_{a',b',a,b}, &{} \forall a \ne a' \in \mathcal {A},b \ne b' \in \mathcal {B}\\ &{} x_{a,b} \in \{0,1\},x_{a,b,a',b'} \in \{0,1\}, &{} \forall a,a' \in \mathcal {A},b,b' \in \mathcal {B}. \end{array} \end{aligned}$$
(3)

Note that, the formulation (3) replaces the quadratic term \(x_{a,b} \cdot x_{a',b'}\) in (2) with the indicator variable \(x_{a,b,a',b'}\), introducing \(O(n^4)\) new variables and \(O(n^3+n^4)\) additional constraints relative to the linear programming formulation. Despite the increase in the size of the problem, this formulation allows applying the primal-dual scheme to efficiently achieve approximate results.

In order to get a primal-dual approximation algorithm for solving the metric labeling in its extended version in (3), we first state the dual of the formulation as follows:

$$\begin{aligned} \begin{array}{lll} \text {max} \quad &{}\displaystyle \sum _{a \in \mathcal {A}} y_a + \sum _{b \in \mathcal {B}} y_b \\ \text {s.t.} \quad &{}y_a + y_b - \sum _{a' \in \mathcal {A}}y_{a,b,a'} - \sum _{b \in \mathcal {B}}y_{a,b,b'} \le c_{a,b}, &{} \forall a \in \mathcal {A},b \in \mathcal {B}\\ &{} \!\!\left. \begin{aligned} &{}y_{a,b,a'} + y_{a,b,b'} + y_{a,b,a',b'} \le w_{a,a'}\cdot d_{b,b'} \\ &{}y_{a',b,a} + y_{a,b,b'} - y_{a,b,a',b'} \le w_{a',a}\cdot d_{b,b'} \\ \end{aligned} \right\} , &{} \forall a,a' \in \mathcal {A},b,b' \in \mathcal {B}\\ &{} y_{p} , y_{a,b,a'} \, , \, y_{a,b,b'} \, , \, y_{a,b,a',b'}\, \text {unrestricted},&{} \forall a,a' \in \mathcal {A},b,b' \in \mathcal {B}\\ &{} y_{a} \ge 0, y_{b} \ge 0, &{} \forall a \in \mathcal {A},b \in \mathcal {B}\end{array} \end{aligned}$$
(4)

Since the variables of type \(y_{a,b,a',b'}\) appears as a summation and a subtraction in the second type of constraints of (4) which accounts for the balancing constraints in (3), strictly following the primal-dual method presented in [14] would require making assignments in tuples since it enforces dual feasibility throughout the algorithm, resulting in poor performance. As we previously suggested in [15], we relax the dual feasibility condition for the first type of the dual constraints that previously became tight and present an efficient primal-dual approximation algorithm for the problem in Algorithm 1.

figure a

The algorithm starts by initializing indicative variables \(x_{a,b}\), set of unassigned nodes \(\hat{\mathcal {A}}\) and \(\hat{\mathcal {B}}\), and an adjusted assignment cost function \(\phi \) where the value of \(\phi (a,b)\) is initially set to be the assignment cost of a to b (line 1). In each iteration of the loop in lines 2–7, the algorithm maps a node a to a node b that minimizes the adjusted assignment cost function \(\phi \) (lines 3–4). Before proceeding to the next iteration, assigned nodes a and b are removed from the sets \(\hat{\mathcal {A}}\) and \(\hat{\mathcal {B}}\) (line 5) and \(\phi \) function is updated for each of the unassigned nodes in the set \(\hat{\mathcal {A}}\) by an amount of separation cost with respect to the recently assigned nodes (line 6). Algorithm iterates until no unassigned node is left in \(\hat{\mathcal {A}}\).

We note that, \(\phi (a,b)\) is not updated for a node a once it gets assigned, rendering the summation \(\sum _{a,b}\phi (a,b)x_{a,b}\) to be equal to the similarity score between the two graphs since it is equal to the value of the objective function in (4) which in turn is equal to the value of the objective function in (3).

3 Results

Here, we demonstrate the utility of our brain network representation and similarity measure on a TBI dataset, where the goal is the binary classification of subjects into healthy controls and TBI patients. We used k-nearest neighbors (kNN) classifier.

3.1 TBI Classification

We used nested leave-one-out approach for cross validation, due to limited number of subjects. For each subject in the dataset, we used the remaining 68 subjects of the dataset as the training set. Using an inner leave-one-out cross validation with training set, we decided the balancing parameter \(\beta \) and the number of neighbors k to be used in the nearest neighbor search. Then, we tested each subject with the learned parameter tuple that achieved best classification accuracy.

For comparison purposes, we performed the experiment using two scenarios. First (baseline), we used only a traditional connectome where we represented edge weights in a vector form without a graph representation. Similarity between subjects is calculated using Euclidean distance between these vectorized edge weights (denoted \(L_2\)-dist). Second, we use enriched connectome with Algorithm 1 (denoted Graph-match). Note that, Graph-match allows regions of the first graph to get mapped to any one of the regions in the second graph while \(L_2\)-dist inherently assumes an identity matching between the nodes of two graphs. The comparison of two scenarios, i.e., traditional connectome with \(L_2\)-dist vs. enriched connectome with Graph-match, facilitates subsequent analysis and interpretation on regional matches between brains, possibly providing insights into TBI-induced regional differences.

Table 1. Classification results from leave-one-out cross validation for the two scenarios: traditional connectome with \(L_2\)-dist vs. enriched connectome with Algorithm 1.
Full size table

Classification performance is presented in Table 1 for the two scenarios, showing overall accuracy, specificity, and sensitivity. Comparing overall accuracy of the two scenarios suggests that our graph representation with the similarity measure captures more information to decide about the classification than the baseline. As suggested from the results, incorporating nodal features into the representation along with connectivity information improves the classification accuracy. In addition to this, relaxing node mappings between enriched connectomes in Graph-match makes it possible to capture subtle regional alterations, possibly associated with injury, which is reflected by the increased classification performance of Graph-match. We also note that, our approach achieves similar performance for classifying patients and controls as the sensitivity and specificity have similar values whereas traditional connectome with \(L_2\)-dist performs poorly for classifying patients. The comparison of ROC curves presented in Fig. 1(a) demonstrates the improved performance of our method over the baseline.

Nested leave-one-out cross validation scheme results in 69 different parameter tuples (\(\beta ,k\)) for our method and 69 k values for the baseline approach. In our experiments, we observed that parameter values were mostly consistent for our method across runs. Specifically, we observed that the inner loop of the experiment has chosen \(\beta =0.9\) without any exception and \(k=15\) with only five outliers out of 69 iterations for our method. This can be contrasted to \(k=6\) being chosen for the baseline approach along with 9 outliers, suggesting the robustness of our graph matching algorithm.

Table 2. Classification results of brain graphs with only graph theoretical features and with only edge weights as features, using Graph-match as the similarity measure.
Full size table

3.2 Effect of Feature Types

In order to demonstrate the contribution of graph theory measures and edge weights as node features, in Table 2, classification results for the brain networks with only graph theoretical features and only edge weights as features are contrasted to both feature types being combined in a single brain graph. We observe that combining both feature types improve the classification accuracy by 10% indicating that enriched connectome maintains more information by combining various features into a single structure relative to a network having either one of them as its only nodal feature. We also observe that using edge weights alone performs better than using graph theoretical measures alone, which can be attributed to larger number of features present in the former, providing a better feature set for classification. However, combination of the two providing an improvement over both of their individual classification accuracies indicate that the two sets of features represent unique aspects of the connectomics.

Fig. 1.
figure 1

(a) Comparison of ROC curves showing the classification performance of the baseline and the proposed method. (b) Z-score distribution of the matching accuracy for controls and patients with respect to the control population.

Full size image

3.3 Mapping Between Nodes of Graphs

We note that, graph matching provides a mapping between nodes of the two graphs in addition to the similarity score between them. One might expect the regions of a brain graph to match their counterparts in another brain graph (such as, precentral gyrus in one enriched connectome would be expected to match with the precentral gyrus of another enriched connectome) as the brain anatomy is similar across people, with occasional mismatches due to subject-specific differences in connectivity. Leveraging this observation, we define another similarity measure, denoted matching accuracy, as the ratio of regions that are accurately matched with their counterparts to total number of regions. Matching enriched connectome of every subject to the healthy controls, we hypothesize that the matching accuracy of healthy controls with respect to themselves should be higher than the matching accuracy of patients with respect to healthy control population, as structural alterations introduced by TBI is expected to cause mismatching regions. As shown in Fig. 1(b), we observe a statistically significant group difference between the patients and controls in their matching accuracy with respect to the healthy subjects. We also observe that the matching accuracy is lower and has a larger variance in patient population, which can be attributed to altered structural connectivity due to pathology.

4 Conclusions and Future Work

In this paper, we presented an enriched connectome that allows combining multiple features into a single structure. The nodes in our representation correspond to the brain regions that are annotated with graph theoretical measures and connectivity of nodes with other nodes as node features, while the edges correspond to the structural connectivity between regions. We also proposed an efficient graph matching algorithm providing two similarity measures over our new representation, one being a summary measure of overall graph similarity and the other quantifying the ratio of number of accurately matched regions to total number of regions. Using the first measure, we showed that proposed enriched representation provided a better classification than the traditional connectomes, demonstrating contribution of the nodal features to information about the samples. Using the second measure, we demonstrated a significantly lower matching accuracy across patients relative to controls, highlighting trauma induced structural alterations in brains of patients.

In this study, we utilized features obtained from a single modality, namely DTI. Our graph representation can easily be extended to combine multiple modalities (e.g., DTI and fMRI). Adding multiple modalities introduce not only new nodal features, but also new edge types that will provide even a richer representation of the brain organization. Although the data that we used in this study involves known correspondences between connectomes, our method can also be applied on connectomes with unknown correspondences, as in subject specific parcellations.