Resting-State Brain Network Analysis Methods and Applications

Abstract

Resting-state fMRI has been widely applied in clinical research. Brain networks constructed by functional connectivity can reveal alterations related to disease and treatment. One of the major concerns of brain network application under clinical situations is how to analyze groups of data to find the potential biomarkers that can aid in diagnosis. In this paper, we briefly review common methods to construct brain networks from resting-state fMRI data, including different ways of the node definition and edge calculation. We focus on using a brain atlas to define nodes and estimate edges by static and dynamic functional connectivity. The directed connectivity method is also mentioned. We then discuss the challenges and pitfalls when analyzing groups of brain networks, including functional connectivity alterations, graph theory attributes analysis, and network-based statistics. Finally, we review the clinical application of resting-state fMRI in neurorehabilitation of spinal cord injury patients and stroke patients, the research on the mechanism and early diagnosis of neurodegenerative diseases, such as multiple system atrophy, as well as the research on brain functional network alteration of glioma patients.

1. Introduction

Magnetic resonance imaging (MRI) is a multimodal technique that can noninvasively reflect the structure and function of the human brain. Structural MRI (sMRI), including longitudinal (spin-lattice) relaxation time T1-weighted and transverse (spin-spin) relaxation time T2-weighted imaging, has been applied to investigate the structural features of the brain. Based on the different relaxation times of different tissue, T1-weighted and T2-weighted images can be used to reflect the volume of grey matter, white matter, as well as lesions caused by infarction or hemorrhage. Diffusion MRI (dMRI), such as diffusion-weighted imaging (DWI) and diffusion tensor imaging (DTI), can be used to measure water diffusion along different directions and tract neural fiber counts and orientation. Functional MRI (fMRI) reflects neural activity during a period of time by measuring the relative amount of deoxygenated hemoglobin and oxygenated hemoglobin in the blood flow, which is also called the blood-oxygen-level-dependent (BOLD) signal. The fMRI is becoming popular in clinical situations to investigate the functional alterations following disease or treatment.

The fMRI experiment can be categorized into task fMRI and resting-state fMRI (rs-fMRI). For task fMRI, subjects need to perform a specific task, such as finger tapping or receive external stimulation like heat or sound during the scanning session. Resting-state fMRI, on the other hand, is collected when the subject lies still in the scanner, without doing any movement or thinking anything particular, and keeping awake at all time. Researchers focus on the spontaneous neural activity reflected by the BOLD signal under resting conditions. The correlation of signals related to spatially distinct regions is commonly defined as functional connectivity (FC) [1].

In the recent two decades, several methods have been developed to analyze functional connectivity in the resting state, including seed-based analysis, independent component analysis (ICA) [2], and resting-state network (RSN) method [3, 4]. The network method characterizes brain spontaneous activity as a graph, where nodes are defined as brain regions and edges are represented as connectivity between regions. There are different ways to calculate the connectivity, including static and dynamic functional connectivity and directed connectivity. Furthermore, features proposed in network science can be adopted to characterize the brain network topology, such as graph theory attributes [5].

Resting-state fMRI has been applied to clinical research and applications [6, 7]. In clinical situations, a common research paradigm is performing group comparison and searching for inter-group significant different features. Researchers are interested in whether a group of patients is significantly different from a group of healthy controls, or whether the same group of patients shows significant recovery after treatment. The identified significant different features may be the potential biomarker to aid in diagnosis as well as treatment. More importantly, the location of the significant different feature is of great interest, since each brain region has its unique function. As a result, this requires comparing groups of brain networks and other extracted network features. In clinical research, there are two key techniques of brain network analysis, the method of network construction and significant difference analysis of groups of brain networks.

In the following sections, we first describe how to construct brain networks from resting-state fMRI data, including different node definitions and a range of connectivity measurements. Then, we present common group analysis methods of brain networks. The clinical application of brain network analysis is also reported. We also propose several future directions in brain network research and end the chapter with a conclusion.

2. Constructing brain networks

Unlike structural and diffusion MRI, the fMRI scanning captures the BOLD signal in a period of time that typically lasts for several minutes. The collected data are a time series, and the “sampling period” is called repetition time (TR). That is, whole-brain data are collected every TR seconds. Before constructing brain networks, the data need to be preprocessed to clean out non-neural artifacts, including physiological signals like breath and heartbeat, head movements, and scanner noise. Then the nodes of the network are defined and connectivity between each pair of nodes is calculated. The whole data processing pipeline is shown in Figure 1.

Figure 1.

Resting-state fMRI data process pipeline.

2.1 fMRI preprocessing

The preprocessing of fMRI data is necessary since there are non-neural noises in the signal. There are openly available toolboxes to carry out preprocessing, such as Statistical Parametric Mapping (SPM), FMRIB Software Library (FSL), and Data Processing Assistant for Resting-State fMRI (DPARSF)[8]. Common preprocessing procedures begin by removing the first 10 time points to let the subject be familiar with the scanning environment. Since the scanning of fMRI data within a repetition period (2s) is done in a slice-by-slice manner, the exact collection time of the first slice and the last slice has a time difference. To correct this difference, a procedure called slice timing correction needs to be performed. Then the head motion is corrected so that each voxel corresponds to the same brain location in the scanning series.

For group analysis, the data of different subjects need to be co-registered or normalized to the Montreal Neurological Institute (MNI) standard space. The data then undergoes smoothing using a Gaussian filter with a specified full-width-half-maximum (FWHM) value. After that, the linear trend in the signal is removed and nuisance covariates, such as white matter, cerebral spinal fluid (CSF), and global signal, are regressed out. At last, the data are filtered to keep signals within 0.01-0.08 Hz, since signals within this frequency range are reported to reflect spontaneous neural activities.

Although numerous preprocessing steps have been developed, there is still no consensus on the standard fMRI data preprocessing pipeline. The controversy is centered on the nuisance covariates regression, especially global signal regression (GSR)[9] and white matter signal regression [10]. Other researchers tried to optimize the preprocessing across multiple outcome measures [11], for low-frequency fluctuation analysis [12] and specific patients, such as stroke patients [13]. We have also investigated how the choices of preprocessing parameters and steps influence statistical analysis results [14]. The preprocessing of fMRI data remains to be a complex but important research topic.

2.2 Node definition

The most basic node definition is the voxel in a 3D fMRI image. Each voxel within the brain can be treated as a node and the constructed voxel-based network covers the whole brain. However, since the spatial resolution of fMRI is relatively high (2mm–4mm), the number of voxels is rather large (around the magnitude of 100,000) and the constructed network requires huge computation power for further analysis. Researchers have proposed specialized methods, such as the Parallel Graph-theoretical Analysis (PAGANI) toolkit to accelerate the processing of voxel-based whole-brain networks [15].

On the other hand, the nodes of the brain network can be defined as regions in the brain. The preprocessed data of voxels within a region are averaged spatially as the signal related to this node. The region can be specified manually by drawing regions of interest (ROI). Independent component analysis (ICA) can also reveal the component region but requires specifying parameters, such as the number of components. Both methods require human intervention and depend heavily on expert knowledge.

We proposed a fuzzy node definition method in Ref. [16] for tumor-brain, named “Spatial-Neighborhood and Functional-Correlation (SNFC)” based on fuzzy connectedness. It is a self-adapting method where the network was divided into functional connection and spatial adjacency. In the SNFC method, fuzzy connectedness between two voxels acts as a measurement to decide if they belong to the same node. Each voxel in the brain could be mapped into two feature spaces—structure feature space $S$ and correlation feature space $C$. Let $s_{i,k}$ represent the spatial relationship between voxel $v_i$ and voxel $v_k$, acting as a judgment of the neighboring relationship. $c_{i,k}$ is the correlation coefficient between the BOLD signal of $v_i$ and $v_k$. The features of structural space $S$ guarantee the principle of the spatial neighborhood and the features of correlation space $C$ ensure the principle of consistency. Fuzzy connectedness between two voxels could be defined as the following:

$$F{C}_{i,k}=s_{i,k}\cdot c_{i,k}$$E1

If $F{C}_{i,k}>T$, then $v_i$ and $v_k$ belong to the same node, where $T$ is the correlation threshold determining whether the correlation of two voxels is strong enough to be in the same node.

The nodes can also be defined using regions in the brain atlas to avoid the subjective error caused by human intervention and enable automatic processing for large cohorts of data. The most known brain atlas is the Brodmann atlas, created by the German anatomist Korbinian Brodmann based on cytoarchitecture [17]. Another popular brain parcellation is the Automated Anatomical Labeling (AAL) atlas [18]. The AAL atlas focuses on brain structure and the finer partition of certain cortices was proposed in AAL2 [19] and AAL3 [20]. Apart from structure, the brain atlas derived from diffusion and functional data is getting more attention. The Brainnetome Atlas was proposed based on DTI data with fine-grained parcellation [21]. Researchers also developed functional atlas, such as the Atlas of Intrinsic Connectivity of Homotopic Areas (AICHA) that considered the homolog of regions in both hemispheres [22]. The above-defined network is called a region-based whole-brain network. We can also construct networks within a region. In this scenario, the voxels are defined as nodes, and the network only consists of voxels within a region. The constructed network is called a voxel-based local network, representing the topology within certain regions. We proposed a multilevel brain network joint analysis method on voxel-based whole-brain networks, voxel-based local networks as well as region-based whole-brain networks (Figure 2) [23].

Figure 2.

Construction of multilevel functional brain networks.

Node definition has a fundamental influence on the topology of the brain network. Different atlas parcels the cerebrum and cerebellum based on different information, and it plays a key role in linking physiological regions to abstract brain network nodes. However, similar to the preprocessing of fMRI data, there is no gold standard for the node definition. Several researches have been carried out to investigate the effect of node definition on network analysis [24], resting-state networks [25], and the topology of both functional networks [26] and structural networks [27]. It is still an open question and needs more thorough research.

2.3 Static and dynamic functional connectivity

Edges in brain networks are represented by the connectivity between nodes. One of the most common connectivity measures is functional connectivity (FC). In 1995, Biswal et al firstly reported the correlation of intrinsic low-frequency BOLD signal fluctuation under resting-state and since then, multiple efforts have been devoted to FC analysis [1, 3]. Functional connectivity is commonly defined as the Pearson correlation between the BOLD signal of spatially distant regions. In recent years, researchers realized that FC ignores the dynamics of neural activity and developed dynamic functional connectivity (DFC) or Chronnectome [28, 29, 30]. The research on DFC is becoming popular and has attracted lots of attention.

Technically speaking, FC or static functional connectivity (SFC) is calculated using the whole time series, whereas DFC utilizes a sliding time window and the correlation of signals within the window is calculated. The window then moves from the beginning of the BOLD signal to the end, with a pre-defined step size. As a result, the connectivity shows dynamic fluctuations as the window slides, and each scanning session is associated with a series of brain networks, or a dynamic brain network. In contrast, there is only one static network related to the scanning session. The network is usually represented by a graph adjacency matrix, which is a square symmetric matrix and the $\left(i,j\right)$ value equals the connection of node $i$ to node $j$. For a dynamic network, there is a time axis along with the adjacency matrix.

There are two major parameters regarding DFC calculation—the window length and the sliding step size. With a longer window length, the dynamics of neural activity might be averaged out while a shorter window length can capture transient signal changes. The step size controls the temporal resolution of DFC. Normally it is specified as several TRs. We investigated the optimal window width by using the small-world property as criteria [31]. Node degree distribution has exponential truncated power-law in the small-world network, and the normal human brain network shows a strong small-world property. The reasonable window width range was verified on both SNFC-based and voxel-based whole-brain networks. Results show that the smallest window width is 200 seconds and 260 seconds for normal subjects and brain tumor patients, respectively. Leonardi et al also studied the theory between window length and filter cut-off frequency during preprocessing [32]. Apart from the two window parameters, the shape of the sliding window is another concern. The rectangular window is the simplest solution, but other choices such as tapered window exist. Mokhtari et al also proposed a modulated rectangular (mRect) window to reduce spectral modulations [33].

We also proposed a dynamic network analysis method for enlarging the training samples required by an unsupervised learning classification algorithm [34], such as a classical backpropagation neural network classifier containing a hidden layer. It reached the optimal accuracy of 100% for classifying glioma patients and normal subjects.

Despite controversies, DFC has been used to investigate diseases, such as schizophrenia [35], post-traumatic stress disorder (PTSD) [36], Parkinson’s Disease [37], and autism [38]. It has also been applied to lifespan studies [39] and cognitive research [40]. From either a methodological or application view, the research on DFC is still insufficient.

2.4 Directed connectivity

As the definition implies, both SFC and DFC contain no directional information. Effective connectivity (EC) can measure the directional influence of one region toward another area by calculating the causal relationships between time series. Commonly adopted EC estimation methods are structural equation modeling (SEM) [41], dynamic causal modeling (DCM) [42], and Granger causality analysis (GCA) [43, 44]. The computation cost becomes unacceptable for SEM and DCM as the number of nodes increases [43]. Several amendments have been proposed to reduce the computation requirement of DCM recently [45, 46], but the model complexity is still challenging for clinical applications. We proposed a method based on convergent cross-mapping (CCM) that can reflect the interactions between regions in a dynamic, nonlinear, and deterministic way, which is not covered by GCA [47]. The method overview, together with the extended network-based statistic, is shown in Figure 3.

Figure 3.

CCM-based directed connectivity estimation and extended network-based statistic method.

CCM was originally developed to detect causality in complex ecosystems [48]. It acts as a complement to GCA as CCM assumes the system to be deterministic and dynamical, while GCA works for a stochastic system and requires separability. In GCA, if removing X decreases the predictability of Y, it can be deduced that X causes Y, and in a brain network scenario, there is a directed connection from X to Y. On the other hand, in deterministic dynamic systems where CCM was developed, we can measure how well Y can estimate X to determine the causal relationship from X to Y, which then determines the directed connectivity strength from X to Y. The procedure of estimating X using Y is called cross-mapping. CCM is also applicable under situations where separability is not guaranteed. GCA, on the other hand, may produce erroneous results [49]. As for the brain, it is a dynamic system whose functional organization is poorly understood [50]. Utilizing CCM to estimate directed connectivity between regions could facilitate the investigation of brain activity as well as enable novel clinical applications.

3. Analyzing group differences in brain networks

After brain network construction, for each scanning of each subject, the preprocessed fMRI data were converted to a brain network represented by a graph adjacency matrix. The next question is how to find the difference between groups of brain networks. Here we summarize two popular methods to further analyze brain networks.

3.1 Significance analysis

The most basic method is analyzing functional connectivity directly. Specifically, suppose we are comparing two groups of networks. Each connectivity value is extracted from every network, forming two sets of values. Statistical hypothesis testing can be adopted to decide whether this connection shows a significant difference as well as which group is higher. After performing a comparison on every connection in the network, the group difference network consisting of significant different connections is obtained. All edges with a significant difference were stored in a network for further discussion. We can also select several regions based on prior knowledge, such as the sensorimotor area or visual area, to further filter the set of significant different connections.

Another method is calculating graph theory attributes. Graph theory characterizes the topology of the network by nodal and global attributes. Common node level graph theory attributes are betweenness centrality, clustering coefficient, local efficiency, modularity, and weighted degree, while the network level graph theory attributes include global efficiency and characteristic path length. Small-worldness is also a common index used in brain network analysis. For multilevel brain networks, we define intra-region features as the attributes calculated at voxel-based local networks, and the attributes calculated at region-based whole-brain networks are called inter-region features. We can calculate the global feature of the voxel-based local network (intra-region features), and the nodal feature of the region-based whole-brain network (inter-region features). As a result, for each graph attribute, we obtain a feature vector whose length equals the number of nodes in the network, representing the whole-brain network feature.

After obtaining feature vectors of graph theory attributes, we can perform a statistical comparison on each region similar to FC analysis. The feature at each region is extracted, forming two sets of values; and statistical testing is used to find significant regions or significant different features. Moreover, the clinical relevance of the features can be evaluated by assessing the correlation of features and clinical scores, which produces features with significant correlation. The intersection of significant different and significant correlated features is selected for further discussion and following analysis.

We also investigated methods to analyze dynamic graph theory attributes [51]. For dynamic brain networks, at each sliding window location, the obtained brain network is static, and graph theory attributes can be calculated. As the window slides, graph theory attributes at each window location are estimated, forming the dynamic graph theory attributes of the dynamic network. To combine static and dynamic attributes together with clinical scores, we proposed an analysis framework [51]. The strength and stability of dynamic graph attributes were calculated. We found significant different and correlated features for both static and dynamic networks, as well as their intersection. The resulting features were further analyzed using receiver-operating curves (ROC) to test their ability in classification.

A controversy regarding the above analysis method is the multiple comparison problem. For each single statistical comparison with a 0.05 significance level, there is a 0.05 chance of obtaining a false positive. However, when performing multiple statistical comparisons at the same time, the chance of getting at least one false positive would become higher as the number of comparisons increases. To tackle this problem, correction methods, such as Bonferroni correction and false discovery rate (FDR) correction, were proposed. The basic idea behind these correction methods is to decrease the single comparison significance level according to the number of comparisons. However, since the amount of comparison is related to the number of nodes in the network, and certain features show high within group variance, directly applying correction might result in no significant result. We argue that statistical comparison can be seen as a feature selection procedure. The significant or selected features are then fed into the next module, such as a classifier. During feature selection, we should keep as much useful information as possible. The uncorrected significant features are preliminary scanning results and taking the intersection of significant different and correlated features further select clinically relevant information. Searching for intersected significant features might be an alternative method to multiple comparison correction.

3.2 Network-based statistics

For brain networks, to overcome the multiple comparison issue, network-based statistics (NBS) was proposed, enabling direct comparison of groups of brain networks [52]. NBS assumes that the effect or the group difference forms a certain structure instead of distributed single connections. The edge-wise comparison is performed first and the links are thresholded according to the test statistics or p-values obtained from the edge-wise comparison, producing a binarized difference network. It then searches for structures or connected components in the binarized difference network. The size of the component, defined as the number of edges or nodes, is used to determine if the component is significant by a permutation test, where group labels of samples are randomly shuffled and the same procedure is performed to search for the maximum component size. The permutation is repeated 5000 times and the empirical distribution of the component size is obtained. An empirical p-value can be assigned to the original connected component by calculating the ratio of the number of permutations, where the maximal size is larger than the original size, to the total permutation number.

Compared with edge-wise comparison and direct edge-wise correction, NBS provides higher statistical power at the cost of coarser spatial resolution in detecting differences [52]. In other words, NBS can only declare the connected component as a whole to be significant. It draws no conclusion on the significance of each single connection within the component. However, the original NBS only works for symmetric adjacency matrices, which corresponds to functional connectivity.

Based on directed connectivity, we proposed the extended-NBS (e-NBS) to search for altered connected components in groups of directed networks [47]. The method overview is shown in Figure 3. We search for strongly connected components (SCC) and weakly connected components (WCC) with and without direction information. A classical depth-first search algorithm was adopted when searching for SCCs and WCCs. The edge-wise p-value was utilized to filter for candidate connections and construct a difference network. Since there is no consensus on how to choose the pre-defined p-value threshold, we changed it within a certain range to test method performance. Specifically, an edge is kept if the p-value is less than the pre-define p-value threshold. For edge-wise comparison, we also tried to use two-sample t-test and the non-parametric Mann-Whitney test. The e-NBS method, together with the CCM-based directed connection estimation method, was verified using a dataset of spinal cord injury patients and healthy controls.

Moreover, we note that given the framework of e-NBS, one can define connected components that suit research needs. For example, in a study of motor function alteration following spinal cord injury, researchers are interested in connections related to sensorimotor areas and visual regions. The connected component can be defined as significant different connections related to these regions of interest. Furthermore, we can define two-step connected components that comprise connections directly related to the ROIs in the binarized difference network, and connections related to regions (first level nodes) that connect with ROIs (Figure 4). Either way, the permutation test in e-NBS makes it possible to draw conclusion on the significance of the defined component.

Figure 4.

Two-step connected component. The first level node is directly connected to the ROI in the binarized difference network, while the second level node is connected with the first-level node.

4. Clinical applications

The resting-state fMRI has been applied to clinical research and applications, mainly investigating pathophysiological mechanisms and searching for sensitive biomarkers for early diagnosis [6, 7]. The prognosis predictability of rs-fMRI is intriguing as well [53, 54, 55]. In glioma research, resting-state fMRI has also shown potential in diagnosis and treatment planning. Here we introduce three examples of applications and related works.

4.1 Neurorehabilitation

It has been shown that changes in both brain function and structure occur following central nervous lesions, such as spinal cord injury [56] and cerebral stroke [57]. According to the theory of neuroplasticity, the brain function continues to change during rehabilitation, and it is the theoretical and physiological basis for individualized neurorehabilitation as well as assistive rehabilitation technologies, such as transcranial direct current stimulation (tDCS) [58, 59, 60] and brain-computer interfaces (BCI) [61, 62]. We performed a study on spinal cord injury patients and investigated the alteration of grey matter volume extracted from structural MRI and functional connectivity related to the sensorimotor area, combining clinical assessments [63]. We found that that the alteration of anatomical structure features and the brain network connectivity in the sensorimotor area were non-concomitant following spinal cord injury, and the functional connectivity within the sensorimotor area had a significant correlation with clinical sensory scores, indicating the potential of FC as a prediction biomarker.

Another issue related to neurorehabilitation is the automated objective evaluation of rehabilitation progress. Traditionally, patient recovery is assessed by clinical measurements, which can only reflect behavioral improvements and might include subjective bias. We proposed a distance-based rehabilitation evaluation method that takes resting-state fMRI data of patients and healthy controls as input (Figure 5) [64]. We hypothesize that the sample point distribution of patients and healthy controls in the feature space is dichotomous. A support vector machine (SVM) classifier was first trained using significantly different functional connectivity of healthy controls and the first scanning session of patients. The distance of the patient sample points to the separating hyperplane was calculated and used to evaluate patient recovery. If the patient recovered, the sample point of the patient would move toward healthy controls and the distance would decrease. The method was verified using both group level and individual longitudinal data, and the distance evaluation was consistent with clinical measurements.

Figure 5.

Method overview of the distance-based rehabilitation evaluation framework.

On the other hand, a stroke could lead to certain movement disabilities. Motor-related brain function alteration after stroke and during recovery is of great interest. Brain-Computer Interface (BCI) systems are helpful in motor recovery, possibly by stimulating neuroplasticity following brain activity [65]. The brain network reorganization of stroke patients after BCI training is of great significance. We conducted an experiment to investigate the functional changes after BCI training and their relations to clinical scores [66]. Functional connectivity was calculated using data collected before and after training and we searched for significant increased FC in groups with and without BCI training. The correlation between FC and clinical scores was also calculated. We found increased FC between certain cerebral and subcortical regions and the inter-hemisphere FC was positively correlated with motor scores.

4.2 Multiple system atrophy

Multiple system atrophy (MSA) is a neurodegenerative disease typically characterized by parkinsonism, cerebellar ataxia syndrome, and autonomic nervous dysfunction [67]. It is further divided into two subtypes, MSA with predominant parkinsonism (MSA-P) and MSA with predominant cerebellar ataxia (MSA-C) [67]. Previous studies mainly investigated the structural abnormalities related to MSA patients and compared subtypes of MSA with Parkinson’s Disease (PD) as well as healthy controls [68, 69, 70, 71]. The functional alteration induced by MSA is also studied by calculating regional homogeneity (ReHo) [72], the amplitude of low-frequency fluctuations (ALFF) [73], as well as functional and effective connectivity [74, 75].

The dynamic functional features of MSA-C patients not thoroughly investigated before. We conducted an experiment on MSA-C patients and proposed a method to combine static and dynamic functional connectivity features, as well as clinical scores (Figure 6) [51]. The static and dynamic brain networks were constructed using methods described in section 2.3 and static and dynamic graph theory attributes were calculated. Statistical comparisons and correlation analysis were carried out and significant different and correlated features were found. The significant regions mainly covered the cerebellum and certain cerebral areas, which is consistent with prior knowledge. The dynamic features showed the highest area under the curve (AUC) value during receiver-operating characteristic (ROC) analysis, indicating the potential of dynamic features in disease diagnosis.

Figure 6.

The coalition analysis of rs-fMRI data combining static, dynamic functional connectivity as well as clinical information.

Apart from structural and functional analysis, multimodal research on MSA is getting more attention. We also tried to combine structural, diffusion, tractography, and functional features extracted from T1, DTI, and fMRI to search for sensitive biomarkers for MSA-C patient diagnosis (Figure 7) [76]. The T1 data were processed to produce grey matter and white matter probability maps. We performed tractography on DTI data and counted the number of tracts crossing each brain region. The fraction anisotropy (FA) and mean diffusivity (MD) maps were also obtained. For rs-fMRI, we calculated functional connectivity and constructed brain networks. The extended network-based statistics for the undirected network were adopted to search for significant different connected components between the two groups. By using the AAL atlas, feature maps extracted from different modalities were converted to feature vectors and networks. After that, significant analysis was performed with false discovery rate correction and we identified significant different features, mainly distributed in cerebellar and certain cerebral regions. The correlation of these features with clinical scores was also tested. We also searched for sensitive biomarkers in disease diagnosis by applying a nested leave-one-out cross-validation framework and evaluated classification performance using the significant features of each region with a support vector machine (SVM) classifier, as shown in Figure 7. The identified biomarkers were mainly cerebellar regions. Different modalities contain complementary information. Merging multimodal data and clinical variables together can further reveal the neurological alteration related to the disease as well as increase the accuracy, robustness, and generalization of the disease diagnosis algorithm.

Figure 7.

The multimodal MRI feature fusion framework and the nested leave-one-out cross-validation procedure. GMV: grey matter volume; WMV: white matter volume; FA: fractional anisotropy; MD: mean diffusivity; NBS: network-based statistics; LOOCV: leave-one-out cross-validation.

4.3 Glioma

Glioma stems from the canceration of neurogliocyte and is the most common tumor in the human brain [77]. It has an intensive impact on the structure of the brain and further on the corresponding physiological functions. Different locations of the glioma will result in different functional alterations and prognosis outcomes. For a high-level glioma, it is highly likely to relapse even after being excised in a surgery [78]. As a result, it is necessary to analyze the brain function changes according to the location and volume of glioma for both diagnosis and treatment. We proposed a framework of multilevel functional network analysis to find the functional network characteristics of glioma patients [79]. The multilevel network consists of a hemisphere functional network, glioma voxel local network, and glioma region local network, as illustrated in Figure 8. The hemisphere functional network was constructed based on regions from a single hemisphere in the AAL atlas excluding cerebellar parcellation (Figure 9). The glioma voxel local network is constructed at the voxel level in the region of glioma that is extracted by a tumor segmentation method. And glioma region local network is also constructed at the voxel level, but within each atlas region containing the glioma. A ratio, defined as the number of voxels in an AAL area that belongs to the segmented glioma region over the total voxel number of the area, is used as the threshold for selecting areas containing the glioma in the AAL atlas.

Figure 8.

A framework of multilevel functional network analysis for finding the functional network characteristics of glioma patients.

Figure 9.

The process of the construction of the hemisphere functional networks is based on the AAL atlas of a glioma patient. The green dots stand for the nodes of the functional network. The yellow line segments represent the weighted edges whose thickness reflects the weight. The colored area shows the tumor region and different colors reflect the possibility of whether a voxel belongs to the tumor.

Network features, including connectivity strength, characteristic path length, average nodal betweenness centrality, and average nodal clustering coefficient, were calculated for all networks. The network connectivity strength was defined as the average z-scores of all edges. Network characteristic path length equals the average of shortest paths between each pair of nodes in the network. Nodal attributes, including betweenness centrality and clustering coefficients, are calculated at each node within the network and averaged as network features. For hemisphere functional networks, both static and dynamic functional connectivity were investigated. Since the period of the BOLD signal induced by the hemodynamic response of neuronal activity is about 20s [80], during the reconstruction of dynamic networks, a sliding window with a length of 50s and a step size of 10s was selected. Each glioma patient received functional scanning lasting for 460s. As a result, the sliding time window extracted 46 sub-signals with a length of 50s and constructed dynamic brain networks with 46-time slices.

In this study, 38 patients with tumors in one side of the brain were enrolled. We constructed 38 positive and 38 negative hemisphere functional networks. Among these patients, 15 subjects had glioma area segmentation. Moreover, 15 healthy subjects were collected as the control group. The local network analysis was performed on 15 patients with segmentation and 15 healthy controls. We used the two-sample t-test to evaluate the significant difference of each feature between hemisphere functional networks constructed on the healthy side and the glioma side. The glioma voxel local networks and glioma region local networks were constructed at the same location of glioma segmentation in data collected from healthy controls as well. Statistical comparison was performed to compare network features of glioma voxel local networks and glioma region local networks from patients and healthy controls. There were 41 glioma region local networks constructed from 15 patients, and for comparison, 41 local networks were estimated from healthy controls.

We also investigated the classification performance using hemisphere functional networks. Given that the sample size is small (38 networks with glioma and 38 networks with healthy tissue), linear support vector machine (SVM) was chosen as the classifier. Static and dynamic network features were extracted and aligned into a feature vector of dimensions 4 and 184 ($46\times 4$) as the input to the classifier, and the leave-one-out cross-validation method is employed to evaluate the performance. The results showed that both dynamic and static features can distinguish the normal and abnormal networks. In addition, dynamic features obtained 100% accuracy in our dataset, while static features showed 71.5% accuracy.

Results revealed by the multilevel functional network analysis method showed that the existence of glioma changed certain features of the normal functional networks. Our work finds that glioma weakened the connection strength of the global and local functional networks. Moreover, it decreased the clustering degree of the nodes in both local functional networks, indicating that glioma may destruct the non-randomness and the small-world property of brain networks.

Previous studies have already investigated how glioma alters functional connectivity [80, 81, 82, 83]. We find that glioma attenuates the connectivity of functional networks, which is in accordance with previous studies. Moreover, we also involved network features other than connectivity. Our study emphasized the characteristic features, such as betweenness centrality, clustering coefficient, and characteristic path length, which were not covered by previous research.

5. Future directions

Despite progress in recent years, there are lots of work to be done in developing new methods for constructing and analyzing brain networks, as well as performing group and individualized analysis. In this section, we propose some possible directions in the field of brain network research.

Network science has been used to analyze brain networks and advanced methods need to be developed to characterize the topological features of brain networks. The algebraic topological data analysis (TDA) method provides a new way to analyze the interactions between a set of nodes instead of bilateral connections. TDA could act as a complement to graph theoretical analysis in describing the topology characteristic of brain networks. More advanced network theory concepts, such as algebraic topology, have also been introduced to the analysis of brain networks [5]. Moreover, artificial neural networks and deep learning methods have been shown to be powerful in analyzing graph data. On the one hand, before network construction, models, such as Recurrent Neural Network (RNN) and Transformer, that were originally proposed to process sequential data, such as natural language and voice, can be applied to analyze the BOLD time series, both with and without preprocessing. Since the network perspective mainly models the inter-relationships between signals of spatially distinct regions, applying deep learning models directly to the time series could possibly extract information complementary to statistical dependency, as described by functional connectivity. On the other hand, after constructing brain networks using functional connectivity, directed connectivity, or DTI fiber tracking, Graph Neural Network (GNN) or Graph Convolutional Network (GCN) could be utilized to merge these multimodal networks and combine both edge-wise features (connections) and nodal features, such as graph theory attributes. GNN was proposed to directly analyze graphs that can model relationships between nodes and perform inference on node, edge, or graph level. Applying GNN to brain networks, especially multilevel static and dynamic brain networks, could possibly extract useful features and enable multimodal information fusion.

On the application side, multiple group comparison methods have been developed. However, for clinical application, individualized diagnosis and treatment are crucial. How to transform conclusions derived from group research into individual situations is a challenging question. We define “healthy templates” as a set of methods to delineate characteristics of a healthy population. The healthy templates describe the distribution of features of healthy people and need to be built for each feature extracted from different modalities. In its most basic form, the healthy template can be a value range given a specific feature. Subjects whose feature value falls within this value range would be considered to be normal, similar to the interpretation of a blood test result. Open-source datasets are valuable resources in the construction of healthy templates. However, the site effect of MRI data is a crucial issue and multi-site data harmonization techniques need to be adopted when combining data from different scanning locations. Several methods have been proposed for harmonization but their utility remains to be tested [84, 85]. With low variance healthy templates, individualized precise treatment planning and prognosis prediction would become possible.

6. Conclusion

The human brain is modeled as a functionally inter-connected network. Resting-state functional magnetic resonance imaging enables observing brain spontaneous activity in vivo. In this chapter, we reviewed the process of rs-fMRI data as well as group analysis methods. Different node definitions and edge estimation were discussed during the network construction stage. Nodes can be defined at the voxel level or with the help of a brain atlas. Lesions, such as glioma segmentation result, can also guide node definition. Edges are estimated in static, dynamic as well as directed scenarios. We presented two major methods to compare groups of brain networks data, significance analysis, and network-based statistics. Combined with the brain atlas, whole-brain networks are characterized by graph theory attributes developed in network science. Network-based statistics enables the direct comparison of groups of brain networks. We also discussed the clinical application of rs-fMRI data analysis in neurorehabilitation, multiple system atrophy, and glioma patients. At last, future research directions are discussed, with an emphasis on network science, novel deep learning models, and individualized clinical applications.

Acknowledgments

This work was supported by a grant from the Tsinghua University Initiative Scientific Research Program (No. 20131089382) and the National Natural Science Foundation of China (No. 61171002, 60372023). We would like to thank Zexuan Hao and Ziliang Zhang from the Department of Electronic Engineering, Beijing National Research Center for Information Science and Technology (BNRist), Tsinghua University for useful advice during the experiment and manuscript development.

Conflict of interest

The authors declare no conflict of interest.