Coactivated Clique Based Multisource Overlapping Brain Subnetwork Extraction
Chendi Wang Rafeef Abugharbieh
q-bio.NC cs.SI


























Expected Precision of Europa Clipper Gravity Measurements
    [
    Received ...; accepted ...
=========================================================


















	
	
	


























epigraph
	

	-0.75in
	

	ifundefinedsinglespacing
    
	
	
	









arabic


Biomedical Signal and Image Computing Lab, UBC, Canada

chendiw@ece.ubc.ca, rafeef@ece.ubc.ca



Subnetwork extraction using community detection methods is commonly used to study the brain's modular structure. Recent studies indicated that certain brain regions are known to interact with multiple subnetworks. However, most existing methods are mainly for non-overlapping subnetwork extraction.
In this paper, we present an approach for overlapping brain subnetwork extraction using cliques, which we defined as co-activated node groups performing multiple tasks. We proposed a multisource subnetwork extraction approach based on the co-activated clique, which (1) uses task co-activation and task connectivity strength information for clique identification, (2) automatically detects cliques of different sizes having more neuroscientific justifications, and (3) shares the subnetwork membership, derived from a fusion of rest and task data, among the nodes within a clique for overlapping subnetwork extraction. 
On real data, compared to the commonly used overlapping community detection techniques, we showed that our approach improved subnetwork extraction in terms of group-level and subject-wise reproducibility. 
We also showed that our multisource approach identified subnetwork overlaps within brain regions that matched well with hubs defined using functional and anatomical information, which enables us to study the interactions between the subnetworks and how hubs play their role in information flow across different subnetworks. We further demonstrated that the assignments of interacting/individual nodes using our approach correspond with the posterior probability derived independently from our multimodal random walker based approach.

Keywords: Clique, Overlapping Brain Subnetwork Extraction, Multisource Fusion, Functional Connectivity, Hypergraph




§ INTRODUCTION

The mainstream of brain subnetwork extraction and standard definition of modularity focus on nonoverlapping definition. However, studies have shown evidences of the existence of overlapping brain subnetworks, hence the methods for nonoverlapping subnetwork extraction are limited by neglecting inclusive relationships <cit.>. There are emerging approaches for discovering overlapping modular network structure, which implies that single nodes may belong in more than one specific module. We here summarize some representative approaches used in brain subnetwork extraction application, and detailed information can be found in a review paper on general overlapping community detection <cit.> . 

The CPM is one of the earliest methods for overlapping community detection <cit.>. It is based on the assumption that communities tend to be comprised of overlapping sets of cliques, fully connected subgraphs. It identifies overlapping communities by searching connected cliques. First, all cliques of a fixed size k must be detected, and a clique adjacency matrix is constructed by taking each clique as a vertex in a new graph. Two cliques are considered connected if they share k-1 nodes. Communities are detected corresponding to the connected components of the clique adjacency matrix. Since a vertex can be in multiple cliques simultaneously, mapping the communities from the clique level back to the node level may result in nodes being assigned to multiple communities <cit.>.

The limitation of CPM is that it operates on binarized graph edges, thus cannot handle weighted graphs <cit.>.

A new definition of modularity has been proposed to discover the overlapping subnetwork based on unbiased cluster coefficients using resting state connectivity <cit.>. However, methods based on the modularity function Q suffer from degenerate partitions and resolution limit <cit.>. 
Another line of studies is to transform a network into its corresponding line graph, where the nodes represent the connections in the original network. Thus, the nonoverlapping community detection (modularity maximization used in <cit.> and agglomerative hierarchical clustering used in <cit.>) on the line graph will result in overlapping subnetworks in the original network. There exist inherent limitations in the nonoverlapping community detection used for the line graph (resolution limit for modularity maximization and local sub-optimum for hierarchical clustering).


Fuzzy community detection algorithms quantify the strength of association between all pairs of communities and nodes <cit.>. Fuzzy k-menas clustering <cit.> and fuzzy affinity propagation <cit.>

have been applied to detect overlapping brain subnetwork extraction. However, one has to use an ad hoc threshold for extracting interacting nodes or independent nodes from the membership vector.


Local expansion and optimization algorithms grow a natural community <cit.> or a partial community based on local benefit functions <cit.>. 
One example is CIS, which has been explored for brain subnetwork extraction <cit.>. Taking each node as a partial subnetwork, CIS expands the subnetwork by determining if any other nodes belong to this existing subnetwork using a local function to form a densely connected group of nodes. Its limitation is the sensitivity to a density factor that controls subnetwork size <cit.>. Another good example is the OSLOM <cit.>, which uses statistical significance of a subnetwork when tested against a global randomly generated null model during community expansion. OSLOM has been shown to outperform many state-of-the-art community detection techniques.

In a previous work from our lab, the RD concept from theoretical biology for modeling the evolution of interacting and self-replicating entities was used to identify subnetworks. Further, the RD formulation was extended to enable overlaps between subnetworks by incorporating a graph augmentation strategy <cit.>. This approach, SORD <cit.>, has demonstrated its superiority over many commonly used overlapping subnetwork extraction methods, including OSLOM.

Most of the algorithms aforementioned are based on one single source, such as resting state functional connectivity. CSORD, the multimodal version of SORD, is one of the few overlapping methods which considers multi-source information. CSORD is based on survival probabilities of different genders in evolution and graph augmentation <cit.>. However, its theoretical background for overlapping assumption based on graph augmentation has relatively indirect neuroscientific justifications. We here explore the direction of integrating multisource information for the overlapping subnetwork extraction by using the straightforward clique concept. 






§ CO-ACTIVATED CLIQUE BASED MULTISOURCE OVERLAPPING SUBNETWORK EXTRACTION

Recent study has indicated that repeatedly activated nodes in different tasks could be canonical network components in the pre-existing repertoires
of intrinsic subnetworks <cit.>, we argue that the clique concept closely resembles groups of nodes which are the canonical network components. Based on the basic observation that typical communities consist of several cliques that tend to share many of their nodes <cit.>, clique-based approach would be a straightforward way to find overlapping brain subnetworks. However, the existing clique-based subnetwork extraction approach CPM (kclique) <cit.> 



has three major limitations that it can only handle binary graphs, but not weighted graphs;  the size k of cliques is fixed, which needs to be adjusted for different types of networks; and it only uses uni-source information. In order to tackle the aforementioned limitations, we here propose a multisource subnetwork extraction approach based on co-activated clique, which (1) uses task co-activation and task connectivity strength information for clique identification, (2) automatically detects cliques with different sizes having more neuroscientific justifications, and (3) shares the subnetwork membership, derived from multisource hypergraph based approach we recently proposed <cit.>, among nodes within a clique for overlapping subnetwork extraction. The schematic illustration of our approach is shown in <ref>.




We first detect co-activated groups of brain nodes across different tasks based on an activation fingerprint idea, and then identify densely connected cliques based on task-induced weighted connectivity. Core cliques are further detected using clique properties we defined.



The nodes within a clique should belong to the same subnetworks due to the close relationship between nodes in a fully connected clique, we thus share the subnetwork membership of nodes within a clique to facilitate overlapping subnetwork assignment. The initial subnetwork membership for each node is derived from non-overlapping subnetwork extraction technique, which is based on the fusion of resting state connectivity and task information embedded with high order relations using hypergraph (see details and notations in <cit.>). 

The difference of our approach from the traditional uni-source kclique method is the utilization of both the task co-activation information

and the connectivity weights (only the binarized connectivity is used in kclique method). The co-activated cliques derived using our approach have flexible clique sizes, which has more neuroscientific justifications than the fixed k. Besides, we explore if our proposed clique node subnetwork membership sharing idea can generate more straightforward and biologically meaningful results than the existing multisource method CSORD.



 §.§ Clique Identification Based on Task Co-activation

We define cliques as co-activated groups of brain nodes that are densely connected in our approach. We first identify the co-activated groups of brain nodes (coarse cliques) using an activation fingerprint idea. Then we refine coarse cliques into cliques, within which nodes are densely connected to each other based on task-induced connectivity information. We denote the clique set as CS, and the coarse clique set as CCS.
Given T different tasks, one can construct a hypergraph with an N× T incidence matrix 𝐇, where N is the number of brain regions and T is the number of tasks (hyperedge e). h(v,e) = 1 when the brain region node v is activated in the task corresponding to hyperedge e. 
The task-induced connectivity matrix 𝐂^task is generated by removing all inter-block rest periods from all regions’ time courses and computing pairwise Pearson's correlations of time courses which were concatenated through block/event durations across all the tasks. 
The underlying assumption for our clique identification is that nodes in the same clique should be co-activated across tasks at times from t = 1… T, where t indicates the number of tasks, in which the nodes are co-activated. 




There are two steps involved in our clique identification, which (1) pre-selects sets of coarse cliques in all T layers, (2) and refines the coarse cliques into cliques.




The approach starts with a pre-selection of coarse cliques CCS, which might include loosely connected nodes that are co-activated. Take each row from the incidence matrix 𝐇 as a activation fingerprint vector f corresponding to the task activation pattern of a node. For example, if one node is activated in the 1st, 3rd and 6th out of the seven tasks, the corrsponding f = [1 0 1 0 0 1 0]. We next operate bit-wise and between the fingerprints from a node pair {i, j}, which gives us an output fingerprint vector of co-activation patterns 𝐒𝐅_ij:

    𝐒𝐅_ij = f_i   f_j,

where f_i and f_j are the activation fingerprint vectors of node pair i and j, and 𝐒𝐅 is the matrix containing the co-activation fingerprint vectors between the nodes in each node pair. We then define a matrix 𝐍𝐓 which counts the number of co-activated tasks between two nodes:

    𝐍𝐓_ij = ∑_t=1^T 𝐒𝐅_ij(t),

where 𝐒𝐅_ij is the co-activation fingerprint vector of length T. Next, we define the node set PS^=t which contains nodes that are co-activated together for t times as:

    PS^=t = ⋃_∀ i,j  s.t.  𝐍𝐓_ij=t{i, j},

and define the node set PS^>t which contains nodes that are co-activated together for greater than 
t times as:

    PS^>t = ⋃_∀ i,j  s.t.  𝐍𝐓_ij>t{i, j}.


Based on the definition above, we follow the four steps as below to identify the coarse cliques.

Step 1 We extract M_t pre-selected sets of co-activated coarse cliques from the nodes in PS^=t. We identify {CCS_1^=t, CCS_2^=t, …, CCS_M_t^=t} by ensuring all the node pairs within a certain set share the same co-activation fingerprint vector in 𝐒𝐅:






    CCS_m^=t = {p_1^m, p_2^m,…,p_N_m^m  |  ∃ p_i^m, p_j^m ∈ PS^=t,  s.t. 𝐒𝐅_p_i^m p_j^m = 𝐒𝐅_p_1^m p_2^m}.

where m=1,…,M_t.
The minimal rank of CCS_m^=t is 2, being only one node pair within a coarse clique. The nodes identified in a coarse clique are fully connected to each other defined by sharing the same co-activation pattern. 


Step 2 Similarly, we extract M_t extended sets of co-activation coarse cliques, 
{CCS_1^>t, CCS_2^>t, …, CCS_M_t^>t}, from the nodes in PS^>t, based on the co-activation patterns between nodes in CCS_m^=t and PS^>t: 






    CCS_m^>t = ⋃_∀ i∈ CCS_m^=t,  ∃ j ∈ PS^>t, s.t. 𝐒𝐅_ij𝐒𝐅_p_1^m p_2^m = 𝐒𝐅_p_1^m p_2^m{j}.

We do not consider the coarse clique set selection for the nodes which only exist in the node set PS^>t for the tth layer, since those will be selected in the pre-selected sets in t+1 layer.






Step 3 We then generate M_t coarse clique sets by merging the pre-selected and extended sets together as:

    CCS^t = {CCS_1^t, CCS_2^t,…,CCS_M_t^t}
    
    CCS_m^t = CCS_m^=t∪ CCS_m^>t, m = 1,…,M_t


Step 4 Extract the coarse clique set CCS across layers in the order from T to 1:

    CCS = ⋃_t=T,…,1 CCS^t.


The second part of clique identification is to refine the coarse cliques into cliques. When we extract CCS = {CCS^1, CCS^2,…,CCS^M}, there still exist loosely connected nodes in the coarse cliques, mostly from lower layers when t is small, especially when t=1. Hence, we subsequently extract cliques based on the strength information from task-induced connectivity matrix 𝐂^task and hypergraph properties. 
We formulate a coarse clique set, CCS^k = {p_1^k, p_2^k,…,p_M_k^k} where there are M_k nodes within, as an M_k× M_k simple graph with the weights between nodes being the task-induced connectivity pairwise edge strength. 
We next apply a local thresholding <cit.> on the M_k× M_k connectivity matrix 𝐂^task k to find out the most closely connected nodes to each node, and binarize the thresholded matrix to generate an adjacency matrix 𝐀^task k. We then transform 𝐀^task k into its hypergraph 𝐇^task k using 𝐀 = 𝐇𝐖𝐇^T - 𝐃_𝐯, where the locations with 1 in each hyperedge correspond to the nodes that comprise a fully connected subgraph, cliques CS^c. We extract N_c cliques:

    CS = {CS^1, CS^2,…,CS^N_c}.









 §.§ Clique Property Computation

We present three properties that can be derived to study the cliques for further network analysis.

(1) Co-activation times NT^c of a clique CS^c, the number of ones in the clique co-activation fingerprint:

    CSF^c = ⋀_∀ i∈ CS^c f_i,

then the co-activation times:

    NCOA^c = ∑_t=1^T CSF^c(t).








(2) Activation times in a clique:

    NA^c = 1/|CS^c|∑_∀ p∈ CS^c∑_t=1,…, T f_p(t).


(3) Clique overlap ratio - the times of a clique overlaps with other cliques divided by the size of a clique, the number of nodes within a clique.
We first define the set of cliques which node i belongs to as a label set: 

    LC_i = {c_1^i, c_2^i,…,c_N_i^i}, c_k^i ∈ 1,…,N_c,

where LC_i is an empty set when node i does not belong to any cliques. We then define the clique overlap ratio as:

    RCO^c = 1/|CS^c||⋃_∀ p∈ CS^cLC_p|.




 §.§ Core Clique Identification

Based on the clique properties, we further identify core cliques out of clique sets for the future overlapping subnetwork extraction. We argue that core cliques should have relatively high co-activation times, high activation times, and high clique overlap ratio. We then devise a core clique selection criterion based on the combination of the clique properties. We normalize all the property values into the range of [0, 1] by dividing individual values by the maximum across all the cliques. The criterion is set as below:

    ρ = _∀ i∈ CS{NCOA^i}/max_∀ i∈ CS{NCOA^i} + _∀ i∈ CS{NA^i}/max_∀ i∈ CS{NA^i} + _∀ i∈ CS{ROC^i}/max_∀ i∈ CS{ROC^i}.


For any clique c which satisfies the criterion:

    NCOA^c/max_∀ i∈ CS{NCOA^i} + NA^c/max_∀ i∈ CS{NCOA^i} + ROC^c/max_∀ i∈ CS{ROC^i} > ρ,

it is selected into the core clique set.




 §.§ Clique Based Overlapping Subnetwork Extraction

Based on the identified core cliques, we further deploy a subnetwork membership sharing technique to identify overlapping subnetworks.
The underlying rationale is that the nodes residing within the same clique behave very similarly to perform some basic functions in tasks, thus, they should be within the same subnetworks. 










In a brain graph with N nodes, let 𝐂^rest be an N× N resting state connectivity matrix, and we have already labeled the non-overlapping subnetwork membership for each node using 𝐂^rest. We have also defined the clique membership of a node i as LC_i in <ref>. We then share the subnetwork membership of the nodes within a clique to facilitate overlapping subnetwork assignment.

First, M_s subnetworks are extracted using non-overlapping community detection approach applied on 𝐂^rest. We define the subnetwork membership of a node i as:

    label(i) = s, i∈ 1,…,N, s∈ 1,…,M_s.

Next, we deploy a sharing scheme of the subnetwork membership label from label(i) of a node i, with the label set of the remaining nodes in the clique where node i belongs to:

    LS(i) = ⋃_∀ c∈ LC_i⋃_∀ p∈ CS^c label(p),

and 

    label(i) = label(i) ∪ LS(i).


We have also explored replacing the resting state connectivity matrix 𝐂^rest with the multisource connectivity matrix 𝐂^t-r defined in <cit.>. We argue that we should further integrate the activation information from task data with high order relation information presented by hypergraph and the rest data when identifying the non-overlapping subnetwork membership.






























§ EXPERIMENTS

We first compare our multisource clique based approach against the uni-source kclique method <cit.>, which is the closest straightforward way to identify overlapping subnetworks. Next we compare against SORD, which has been proven to outperform the state-of-the-art techniques such as OSLOM, and CSORD (the multisource version of SORD) <cit.> to see if our proposed approach have more direct biological intuition for the overlapping subnetwork extraction. We also examine the nodes within subnetwork overlaps derived by our approach by assessing the probability of a node belonging to subnetworks using our recently proposed multimodal RW approach <cit.>, to verify that our overlapping subnetwork assignments correspond with the posterior probability.


 §.§ Materials

We used the resting state fMRI and task fMRI scans of 77 unrelated healthy subjects from the HCP dataset <cit.>. Two sessions of resting state fMRI with 30 minutes for each session, and 7 sessions of task fMRI data were available for multisource integration. The seven tasks are working memory (total time: 10:02), gambling (6:24), motor (7:08), language (7:54), social cognition (6:54), relational processing (5:52) and emotion processing (4:32). 
Preprocessing already applied to the HCP fMRI data includes gradient distortion correction, motion correction, spatial normalization to MNI space with nonlinear registration based on a single spline interpolation, and intensity normalization <cit.>.
Additionally, we regressed out motion artifacts, mean white matter and cerebrospinal fluid confounds, and principal components of high variance voxels using compCor <cit.>. Next, we applied a bandpass filter with cutoff frequencies of 0.01 and 0.1 Hz for resting state fMRI data. For task fMRI data, we performed similar temporal processing, except a high-pass filter at 1/128 Hz was used. The data were further demeaned and normalized by the standard deviation.
We then used the HO atlas <cit.>, which has 112 ROIs, to define the brain region nodes. We chose the well-established HO atlas because it sampled from every major brain system, and consists of the highest number of subjects with both manual and automatic labelling technique compared to other commonly used anatomical atlases. Voxel time courses within ROIs were averaged to generate region time courses. The region time courses were demeaned, normalized by the standard deviation. Group level time courses were generated by concatenating the time courses across subjects. The Pearson's correlation values between the region time courses were taken as estimates of FC matrices. Negative elements in all connectivity matrices were set to zero due to the currently unclear interpretation of negative connectivity <cit.>.
For task activation, we applied the activation detection on the seven tasks available following the steps described in <cit.>. 

We further applied local thresholding <cit.> on 𝐂̅^task by setting graph density to be 0.1 to generate the hypergraph when we identified cliques from the coarse clique set. We selected a relatively strict threshold to only select those most closely connected nodes to form cliques. 0.1 has been chosen based on the cross-validation on inter-subject reproducibility within the range between 0.03 to 0.2 at the interval of 0.01.

The non-overlapping subnetworks were derived from the resting state connectivity matrix 𝐂^rest or multisource connectivity matrix 𝐂^t-r (generated using strategies from <cit.>) using Ncuts, when the number of subnetworks was set to 7, same as the abaialble number of tasks.



§ RESULTS

We compared the overlapping subnetwork extraction using our proposed MCSE approach with 𝐂^t-r, or 𝐂^rest against the uni-source kclique approach <cit.>, SORD <cit.>, which has been demonstrated to outperform state-of-the-art overlapping community detection methods including OSLOM when applied to brain subnetwork extraction, and CSORD <cit.>, the multisource extension of SORD. Two uni-source approaches extract overlapping subnetworks using resting state data. The parameters for kclique were set using the cross-validation on the clique size k based on inter-subject reproducibility from the suggested range [3,…,6] <cit.> and reasonable graph densities from 0.03 to 0.2 at the interval of 0.01. SORD and CSORD applied 100 bootstraps by sampling with replacement as suggested in <cit.>.
We also evaluated the probability of a node being assigned to a subnetwork using our recently proposed RW based approach <cit.> to examine the proposed clique-based overlapping subnetwork identification.
All statistical comparisons are based on the Wilcoxon signed rank test with significance declared at an α of 0.05 with Bonferroni correction.


 §.§ Comparison with Existing Overlapping Subnetwork Extraction Methods

We quantitatively evaluated the contrasted approaches based on test-retest reliability and inter-subject reproducibility, since ground truth subnetworks are unknown for the real data of human brain.


  §.§.§ Group-level Subnetwork Extraction Reproducibility

We first assessed the test-retest reliability based on group level subnetworks extracted separately from two sessions of rest and task data (each of the seven tasks includes two sessions of fMRI data) using DSC. 

The subnetworks extracted from the first session's data are taken as the “ground truth”, against which the subnetworks from the second session are compared. 

We found that our proposed MCSE outperforms all other contrasted approaches, by achieving a DSC between subnetworks extracted from two sessions of data at 0.8917 with 𝐂^t-r and 0.8865 with 𝐂^rest, against kclique at 0.7514, SORD at 0.8378, and CSORD at 0.8514, see <ref>. 









  §.§.§ Subject-wise Subnetwork Extraction Reproducibility

We assessed the inter-subject reproducibility by comparing the subnetwork extraction results using subject-wise data against the group level data, <ref>. 
The average DSC between subject-wise and group level subnetworks across 77 subjects based on five approaches are MCSE with 𝐂^t-r at 0.7024±0.0722, MCSE with 𝐂^rest at 0.6281±0.0583, kclique at 0.4967±0.0430, SORD at 0.5129±0.0774, and CSORD at 0.5952±0.0901, respectively. MCSE with both 𝐂^t-r and 𝐂^rest are found to achieve statistically higher inter-subject reproducibility than constrasted approaches based on the Wilcoxon signed rank test at p<10^-10 and p<0.005. respectively.
Further, MCSE with 𝐂^t-r outperforms 𝐂^rest at p<0.00001, which confirms the benefit of incorporating the task information embedded with higher order relations in assigning non-overlapping subnetwork membership.











 §.§ Biological Meaning - Analyzing Function Integration




We further examined the biological meaning of the overlapping subnetworks found using all five methods, our proposed MCSE with 𝐂^t-r, MCSE with 𝐂^rest, kclique, SORD and CSORD, <ref>. We first measured the overlapping ratio by dividing the number of nodes residing in the subnetwork overlaps by the total number of brain regions detected in subnetworks. The ratio of five methods are 0.3482, 0.3482, 0.4444, 0.4328 and 0.2885. Our proposed approach can generate the similar ratio of interacting nodes which reside within subnetwork overlaps to the existing overlapping methods. We note that CSORD generated relatively smaller number of interacting nodes, the possible reason is that the strict stability selection resulted in exclusion of some meaningful nodes, which were taken as false detected nodes arising from noise <cit.>.

By examining the locations of those interacting nodes, we found that our proposed MCSE with 𝐂^t-r approach identified subnetwork overlaps within pre- and postcentral gyri, medial superior frontal cortex, inferior frontal gyrus, superior parietal lobule, precuneous, lateral occipital cortex, occipital pole and frontal orbital cortex; which match well with functional hubs previously identified by graph-theoretical analysis based on the degree of the voxels <cit.>. Besides, brain regions of insula, putamen, thalamus, supramarginal gyrus have been found within subnetwork overlaps, which match well with the connector hubs identified using the centrality measures <cit.>. The results of using MCSE with 𝐂^rest is very similar to 𝐂^t-r, only that precuneous cortex was missed, and the temporal pole was misclassified into the subnetwork overlaps. This result confirms the benefit of integrating the information from both task and rest data. Both MCSE methods also identified 

lingual gyrus and fusiform cortex around as interacting nodes. Lingual gyrus was identified as a hub based on cortical thickness correlation <cit.> and the fusiform cortex within occipitotemporal cortex has been found to be intermediary “hub” linking visual and higher linguistic representations <cit.>. 















As for the traditional kclique approach, biologically meaningful subnetwork overlaps were found within inferior frontal gyrus, superior and middle temporal gyri <cit.>, supramarginal gyrus, insula <cit.>, inferior temporal gyrus <cit.>, and occipitotemporal cortex <cit.>. kclique failed to identify all the other aforementioned (connector) hubs which were found using our proposed methods. Instead, regions normally were not considered to reside in subnetwork overlaps were found, such as temporal fusiform cortex, central opercular cortex, and parietal operculum cortex. On the other hand, this kclique approach detected angular gyrus (functioning as a semantic hub) within subnetwork overlaps.




SORD was able to find subnetwork overlaps within inferior, superior and middle temporal gyri, superior parietal lobule, lateral occipital cortex, occipital pole and lingual gyrus that match well with functional hubs, but failed to find other hub regions identified by MCSE. Instead, SORD detected many regions as interacting nodes, which normally are not considered as hubs, such as intracalcarine cortex and cuneal cortex in the visual system, and regions in language related system, including central opercular cortex, parietal operculum cortex, planum polare, planum temporale, heschls gyrus, and supracalcarine cortex.



With relatively lower number of overlapping ratio, CSORD identified biological meaningful subnetwork overlaps within regions such as pre- and postcentral gyri, middle temporal gyrus, angular gyrus and lateral occipital cortex, while failed to find any other hubs. Similar to SORD, it included some regions in language related system to the subnetwork overlaps, such as central opercular cortex, parietal operculum cortex, and planum temporale. We did not discover the single subnetwork constituting the visual corticostriatal loop, striatothalamo-cortical loop, and cerebello-thalamo-cortical loop, which was found in <cit.>. The reason could be this connection was reflected in AC, instead of task functional connectivity. 





Collectively, our proposed MCSE approach is able to identify subnetwork overlaps which constitute more biologically meaningful brain regions, such as hubs, compared against contrasted methods.



 §.§ Comparison Between the Subnetwork Overlaps and the RW Posterior Probability


We also examined the overlapping subnetworks derived from our approach MCSE with 𝐂^t-r by assessing the probability of a node belonging to subnetworks using our own recently proposed multimodal RW approach <cit.>, to verify that our overlapping subnetwork assignments correspond with the posterior probability. The underlying rationale is that for an interacting node, which resides within the subnetwork overlaps, its probability of belonging to a subnetwork will be distributed across the subnetworks it resides in. On the other hand, an individual node, which does not reside within subnetwork overlaps, would have higher chances to possess a dominant probability of being assigned to a particular subnetwork. Hence the difference of probabilities of a node being assigned to the first two subnetworks with the first two highest probabilities indicates the possibility of a node residing within subnetwork overlaps. Interacting nodes tend to have a smaller value of difference of first two highest probabilities.

We here define the degree of overlapping confidence as the subtraction from one of the difference between the first two highest probabilities of a node being assigned to subnetworks. The nodes identified within the subnetwork overlaps (interacting nodes) and outside of the overlaps (individual nodes) are considered as two populations. For each population, the average overlapping confidence is defined as below in <ref>:

    overConf = 1/|S|∑_i ∈ S (1-(p_i^max-p_i^smax)),

where S is a set of nodes, either nodes residing within or outside the subnetwork overlaps, and p^max is the maximal probability of a node belonging to subnetworks, and p^smax is the second maximal probability. Thus, the interacting node population is expected to have higher overConf compared to individual nodes.

We first derived the probabilities of each node being assigned into all possible subnetworks using our recently proposed multimodal RW approach <cit.>, where two sources of connectivity matrices are 𝐂^rest and 𝐂̅^task, matching with how 𝐂^t-r was generated in our approach. The number of seeds within each subnetwork n_k was set to [2,…,9], where 9 is 75% of 12, the minimal number of nodes which are included in non-overlapping subnetwork extraction using 𝐂^t-r.

We found that the overlapping confidence of the interacting nodes with an average overConf of 0.6884 are statistically higher than the individual nodes with an average overConf of 0.6338 based on the Wilcoxon signed rank test at p=0.006,

see <ref>. This finding confirms that the overlapping subnetwork assignments based on our proposed MCSE match with the probability derived independently from our RW based approach.  

















§ DISCUSSION



 §.§ Benefits of Clique Identification Based on Task Co-activation

The traditional definition of clique is the fully connected subgraphs identified by the connections between brain regions mostly on resting state connectivity. In our approach, we present a novel way to identify cliques based on the similarity of activation patterns between nodes. We argue that the clique concept closely resemble the canonical network components that are recruited selectively and repeatedly in different task-induced activities <cit.>. Different from the traditional kclique method <cit.>, our clique-based approach is able to utilize both the task activation information and the task-induced connectivity strength rather than only the binarized connectivity information used in kclique method. Besides, the cliques derived using our approach have flexible clique size, which was determined automatically, having more neuro-scientific justifications than the fixed clique size in <cit.>. Moreover, we estimate the properties from cliques to indicate the importance of cliques, which gives us a better control over falsely including some fake cliques due to noise. We did find cliques within brain areas that well match with hubs, which indicates that our approach can identify subnetwork overlaps with more biological meaning than the traditional kclique method.


 §.§ Multisource Information Integration Improves the Overlapping Subnetwork Extraction

Compared to the widely used overlapping community detection methods, our approach integrates information from multiple sources. We used both task information (including task activation and connectivity strength) for clique identification and resting state connectivity information for subnetwork membership sharing. The results from reproducibility and biological meaning indicate that our multisource approach, especially MCSE with 𝐂^t-r, outperforms uni-source methods such as kclique and SORD, which has been proven to give better overlapping brain subnetwork extraction results compared to state of the art techniques such as OSLOM. We note that our multisource approach further outperformed the multisource version of SORD, CSORD. The reason could be that clique based idea and the sharing of the node subnetwork membership is more straightforward and have more direct biological intuition than relying on survival probabilities of different genders in evolution, which is used in CSORD.


 §.§ Overlapping Subnetwork Identification Corresponds with the RW Posterior Probability



We have identified subnetwork overlaps within brain regions that well match with hubs defined using functional, structural and anatomical information. The results enable us to study the interaction and integration between subnetworks and how interacting nodes (or important hubs) play their roles in the information flow across different subnetworks. We further demonstrated that the assignments of interacting/individual nodes using our proposed MCSE correspond with the posterior probability derived independently from our previously proposed RW based approach <cit.>. 

The finding of more distinguishable overlapping confidence between two populations of nodes when the number of seeds was set within a range of [6, 8] confirms the merit of using multiple seeds within a reasonable range (not including connector hubs) in the RW based approach.


 §.§ Other Considerations

We have also discovered that the uni-source traditional kclique approach has high computational complexity when the graph density increases, where there exist large number of fully connected subgraphs. The computational complexity of both SORD and CSORD increases when the bootstrap sampling increases <cit.>. However, the computation time of our proposed MCSE is quite reasonable and not sensitive to the graph densities.

In terms of the coverage of the brain area from the subnetwork extraction results, SORD and CSORD neglected some brain regions which are not selected as significant nodes by the stability selection. While these two approaches offered this extra feature, they sometimes falsely missed important nodes and failed to cover the whole brain for analysis. 



 §.§ Limitations and Future Work

In this work, we presented an approach to identify cliques based on task information and extract overlapping subnetworks using both task and rest data. However, the ideal multimodal framework would be able to integrate AC into the fusion for detecting overlapping subnetworks. The challenge is to discover the relationship between AC and task activation, which enables the clique identification to incorporate anatomical information. Our future work will focus on integrating AC into task-activation based clique identification, or into the multimodal subnetwork membership assignment, using the multimodal RW approach or the multislice approach <cit.>.


§ CONCLUSION

We proposed an approach for multisource overlapping brain subnetwork extraction using canonical network components, cliques, which we defined based on task co-activation. Based on the clique concept, we investigated overlapping subnetworks based on a label sharing scheme which incorporates the rest data information and task data embedded with higher order relations.
We have demonstrated that integrating multimodal/multisource information and using high order relations result in better subnetwork extraction in terms of the overlaps to well-established brain systems, test-retest repeatability, inter-subject reproducibility and biological meaning.


§ LIST OF ACRONYMS











































AALAutomated Anatomical Labeling
ACAnatomical Connectivity
ACCAnterior Cingulate Cortex
ADAlzheimer’s disease
BOLDBlood Oxygenated Level Dependent
BPBoundary Point
CBPconnectivity based parcellation
CgCingulate Cortex
CISConnected Iterative Scan
CNNConvolutional Neural Networks
COREGCO-training with REGularization
CPMClique Percolation Method
CSAConstant Solid Angle
CSORDCoupled Stable Overlapping Replicator Dynamics
DAEDeep Auto-Encoder
DALYsDisability-Adjusted Life Years
DBNDeep Belief Network 
DCphysical distance
DMNDefault Mode Network
dMRIDiffusion-weighted Magnetic Resonance Imaging
DNNsDeep Neural Networks
DOFdegree of freedom
DSCDice Similarity Coefficient
DTIDiffusion Tensor Imaging
ECEffective Connectivity
ECNExecutive Control Network
EEGElectroencephalography 
EPIEcho-Planar Imaging
FAFractional Anisotropy
FCFunctional Connectivity
fMRIFunctional Magnetic Resonance Imaging
GBDGlobal Burden of Disease
GLMGeneral Linear Model
GMMGaussian Mixture Model
GROUSEGrassmannian Rank-One Update Subspace Estimation 
GTGlobal Thresholding
HARDIHigh Angular Resolution Diffusion Imaging
HCPHuman Connectome Project
HOHarvard-Oxford
ICAIndependent Component Analysis
ICCIntra-Class Correlation
ICNsIntrinsic Connectivity Networks
IPInterior Point
IPLInferior Parietal Lobule
IQIntelligence Quotient
LECNLeft Executive Control Network
LMaFitLow-Rank Matrix Fitting
LTLocal Thresholding
MCMultimodal Connectivity 
MCNFMatrix Completion with Nonnegative Factorization
MCSEMultisource Clique-based Subnetwork Extraction
MEGMagnetoencephalography
MITKMedical Imaging Interaction Toolkit
mmRWmulti-modal Random Walker
MNIMontreal Neurological Institute
MRIMagnetic Resonance Imaging
MST-KNNminimal spanning tree and k-nearest neighbors
MVSCMultiView Spectral Clustering
NcutsNormalized cuts
NMDNeighborhood-information-embedded Multiple Density
NMFNon-negative Matrix Factorization
NMINormalized Mutual Information
NPNon-deterministic Polynomial-time
NRMSEnormalized root-mean-squared-error
ODFOrientation Distribution Function
OSLOMOrder Statistics Local Optimization Method 
PCAPrincipal Component Analysis
PCCPosterior Cingulate Cortex
PDParkinson’s disease
PDDPrincipal Diffusion Direction
PDFProbability Density Function
PETPositron Emission Tomography
RBMRestricted Boltzmann Machines
RDReplicator Dynamics
RMrandom parcellation
ROIregion of interest
RQResearch Questions
rs-fcMRIResting State Functional Connectivity based on MRI
RWRandom Walker
SAEStacked Auto-Encoder
sMRIStructural Magnetic Resonance Imaging
SNRsignal-to-noise ratio
SORDStable Overlapping Replicator Dynamics
SVRSupport Vector Regression
t-fcMRITask Functional Connectivity based on MRI
TRrepetition time
WHOWorld Health Organization
3Dthree-dimensional


















splncs



