Voxelwise-based Brain Function Network using Multi-Graph Model

In the research of the fMRI based brain functional network, the pairwise correlation between vertices usually means the similarity between BOLD signals. Our analysis found that the low (0:01–0:06 Hz), intermediate (0:06–0:15 Hz), and high (0:15–0:2 Hz) bands of the BOLD signal are not synchronous. Therefore, this paper presents a voxelwise based multi-frequency band brain functional network model, called Multi-graph brain functional network. First, our analysis found the low-frequency information on the BOLD signal of the brain functional network obscures the other information because of its high intensity. Then, a low-, intermediate-, and high-band brain functional networks were constructed by dividing the BOLD signals. After that, using complex network analysis, we found that different frequency bands have different properties; the modulation in low-frequency is higher than that of the intermediate and high frequency. The power distributions of different frequency bands were also significantly different, and the ‘hub’ vertices under all frequency bands are evenly distributed. Compared to a full-frequency network, the multi-graph model enhances the accuracy of the classification of Alzheimer’s disease.

• The study of brain functional networks using frequency division method found that the BOLD signal similarity matching of voxels in different frequency bands are different. • A multi-graph-based brain functional network model is proposed. Using the multi-graph model to establish a voxel-wise brain functional network with the combination of low, intermediate, and high-band signals, and a multi-band study of brain function signals is implemented. • The analysis of the proposed model shows that the topological characteristics of the voxel-wise based brain functional network in the low, intermediate, and high-band signals of the model have obvious differences, which further demonstrates the effectiveness of the crossover network model.
The rest of the paper is organized as follows. In the section: Statistical analysis of pairwise correlation we introduce the experiment and statistics. The details of the proposed multi-graph model are introduced in the section: Multi-graph brain functional network model. The model analysis section shows the model properties which were used in the experiment. Finally, we conclude this paper in the Section: Conclusion.

Statistical Analysis of Pairwise Correlation
In this section, the data sources and preprocessing is introduced first, and then the similarity of voxel BOLD signals in different frequency bands are discussed.
Functional magnetic resonance imaging data. The experimental data. fMRI data were obtained from the open access NITRC 21 : 1000 Functional Connectomes Project (https://www.nitrc.org/). 662 subjects were used as shown in Table 1. Ethical statements are present in the NITRC, and we confirmed that all experiments were performed in accordance with relevant guidelines and regulations of the access. We also confirmed that all experimental protocols were approved by Sino-Dutch Biomedical & Information Engineering School, Northeastern University, China.
Preprocessing. For each subject, the first 7 volumes of the functional images were discarded for signal equilibrium and to allow the participant to adapt to the experiment, leaving rest volumes for next steps. Standard preprocessing was applied on rs-fMRI dataset of all patients using Data Processing Assistant for Resting State fMRI (DPARSF) toolbox 22 and Statistical Parametric Mapping Software Package package (SPM12) (http://www.fil.ion. ucl.ac.uk/spm) 23 . Slice-timing correction to the last slice was also performed. fMRI time-series realigned using a six-parameter rigid-body spatial transformation to compensate for head movement effects. Then all images were normalized into the Montreal Neurological Institute (MNI) space and re-sampled to 3-mm isotropic voxels. A band-pass filtered (0.01-0.2 Hz) was used to retain the relevant information.
Pairwise correlation between vertices. In the establishment of a voxel-wise brain functional network, the similarity of voxel BOLD signals is usually used as the pairwise correlation between vertices, and the signal is usually regarded as a full frequency band signal (filtered image). First, Fourier transform is used to transform voxel SCIENTIfIC REPORTS | (2018) 8:17754 | DOI: 10.1038/s41598-018-36155-z signals into frequency domain. Second, the frequency domain signal is split into low, intermediate and high ranges. Third, the frequency ranges is restored to the same time domain signal with the same length through inverse Fourier Transform, which is the corresponding frequency band signals, the brain functional networks are established according to the correlation between these separated band signals. Two voxels in the same frequency band are represented as x(x 1 , x 2 , …, x n ) and y(y 1 , y 2 , …, y n ). The similarity of the BOLD signal variation sequence is calculated by Euclidean Distance proposed by Muldoon et al. 1 , and we can get the correlation by setting the threshold in the equation below.
Euclidean distance between vertex x and y is defined as When the d(x, y) is smaller than the threshold d(x, y) <λ, then we can say that there is a functional connection, i.e. the pairwise correlation.
Statistical analysis of similarity of BOLD signal in different frequency bands. In the study of BOLD signal in brain voxel, the low frequency range shows a very high signal intensity, which is difficult to see in other frequency domain, displayed in Fig. 1(a). When the low frequency range is removed, and when the signal is larger than 0.15 Hz, the waveform has an upward trend in the interception signal ranges, shown in Fig. 1(b). in the study of Achard 17 , the small world properties in 0.01-0.06 Hz frequency are similar, and the small world properties in 0.06-0.23 Hz frequency are falling, and the changes in the amplitude around 0.06-0.15 Hz and 0.15-0.2 Hz are different from our experiments, as shown in Fig. 1(c,d).
Since there are obvious differences in different frequency ranges, it is necessary to study the differences of BOLD signals in the frequency bands. We found that there are many situations that are similar in a matching of functional signals on different frequency bands. The voxels with different signal situations were sampled by the  voxel's BOLD signal of brain network. In each situation, the two voxels were represented by red and blue in Fig. 2. As shown in Figure 2(a)I-III, signal matches in the full frequency band signal, when the two voxels of brain are matched at low, intermediate and high frequencies respectively. We can also see that the signal also matches in the full frequency band, that is shown in Fig. 2(a)IV. Meanwhile, we also found that some voxels have another state, as shown in Fig. 2(b), two BOLD signals are not similar in the full frequency band. However, after the signals are separated by low, intermediate and high frequency band, their signals similar in all three bands. We also found that when some signals are similar in the full frequency band, the signals are not similar at low, intermediate and high frequency bands and results are shown in Fig. 2(c). In the study, there are also quite a few signals which are not similar under the low frequency band, but are similar in intermediate and high frequency band, the results show high similarity, and it also shows similarity in full frequency band signal as shown in Fig. 2(d).
Therefore, it is necessary to consider the effect of this BOLD signals similarity difference on the pairwise correlation in functional networks. In order to summarize and analyze the problems found in the experiment, the brain functional networks about 14000 vertices are established at low, intermediate, high frequencies and full frequency bands respectively under the same threshold λ. The low frequency network correlation matrice is represented as H L , the intermediate frequency network is represented as H I , the high frequency network is represented as H H , and the full frequency bands network is represented as H Full .
The results displayed in the Table 2 are the number of pairwise correlations (i.e. edges) between vertices in the brain network mentioned above. First, as shown from rows 1 to 4 in Table 2, the edges contained in H L , H I and H H is different from H Full , which shows the difference of connectivity patterns in different frequency bands. Secondly, from rows 5-15, we can determine that there is a large number of overlaps at the edges of the network in different frequency bands. Especially when compared with H Full , H L covers 96.91% of the edges of H Full . The intersection ratios of H I and H H are 13.18% and 15.63% respectively. From the results above, we can determine that the impact of low frequency information of the signal in H Full is more obvious.
When we assemble different frequency bands together, the coverage ratio is also different. The intersection ratio of H L and H I covers 12.24% of H Full , and intersection ratio of H L and H H is 15.1%. These two situations are higher than those of H I and H H , who has an intersection ratio of 3.54%. It can be seen that the influence of high frequency signal on the full frequency band signal is higher than intermediate frequency signal on full frequency band signal. As shown in Table 2 Table 2, a lot of edges cannot be displayed in the H L , H I , H H network. Therefore, due to the coverage of low frequency information, the use of full band frequency signals to match nodes will lose a large amount of information.
The use of p-values is common in statistical hypothesis testing for significance, the p-values are the probability, reflecting the probability of an event. According to the p-values obtained by the significance test method, p > 0.05 is generally considered to be no statistically significant differences, p < 0.05 is generally considered to be statistically significant differences, p < 0.01 is considered to be extremely statistically significant differences. Therefore, we use p-values to calculate whether there is a significant difference between networks with different frequency bands. In summary, the smaller the p-values, the more significant the result. From Table 3, the differences between H L and H Full are not obvious, as the p-value is 0.0672. And by the p around 0.01, H L differs significantly from H I and H H . There is also a significant difference between H I and H H , but it is lower than H L .

Multi-Graph Brain Functional Network Model
Multi-graph network model. In order to describe the brain network better, a voxel-wise based brain crossover functional network by the multi-graph model is proposed in this paper, in which the low, intermediate and high frequency information is integrated. The matching information between vertices is described by multiple edges, which form the network model by multi-graph. According to the different working modes, the number of edges between different vertices may also be different.
Namely G will be decomposed into k subgraphs, so that H 1 , H 2 , …, H k are the subgraphs of G. H 1 , H 2 , …, H k are pairwise edge disjoint. Using this method we can realize the full coverage of the brain network, and can also decompose the multiple possibilities of the brain network operating mode. H 1 , H 2 , …, H k can be used to represent different brain activities.
There are two patterns of work when studying the connectivity of brain network vertices.

Multi-graph model analysis.
According to the multi-graph model mentioned above, H 1 , H 2 , …, H k can be used to represent different brain activities, whose organization forms represent different brain network models. Therefore, the proposed multi-graph network model based on low, intermediate and high frequency bands H L , H I , H H is discussed with full frequency band signal H Full .
Small-world properties. We can see from the existing research, the brain functional network has small-world properties. Therefore, the small world attributes of the low, intermediate and high frequency network will need to be analyzed 1 . In the analysis of small world networks, there are generally two parameters: Clustering coefficient and the shortest path. The small world property of human brain network has important function meaning.

• Clustering coefficient
According to the graph theory, clustering coefficient C is a coefficient of the vertices aggregation degree in a graph 1 . Clustering coefficient means separation module, that is, the network has the ability of efficient modular processing of information. The network that has a higher average clustering coefficient was found to have a modular structure in the study. The calculation results are shown in Fig. 5.

• The shortest path
The shortest path length L represents the number of edges required by communicating between any two vertices in the network 1 . The length of the shortest path means that there is a fast and efficient information transfers function between the cerebral regions within the network, and the ability to transfer information    is expressed as C L R , C I R , C H R . The shortest path is expressed as L L R , L I R , L H R . Typically, in a small-world network, we expect the ratio γ = C/C R > 1 and the ratio λ = L/L R > 1. A scalar summary of small-worldness is therefore the ratio σ = γ/λ, which is typically >1.
It can be seen that all subgraphs retain better small world properties under voxel-wise network, as shown in Fig. 7.
First, as shown in Figs 5 and 7(a), the low frequency subgraph H L has a higher clustering coefficient than others. Due to the higher clustering coefficient, the node connection appears a higher degree of modular distribution, which means that brain function network mainly forms a modular structure to complete the work under the subgraph H L . The modular connectivity of other subgraphs is less than H L . Second, the subgraph H I and H H have a lower clustering coefficient than the subgraph H L , which means that the modular distribution degree of H I and H H is lower than that of the subgraph H I . Thus it can be seen that the network structures are more global under the subgraph H I and H H . The working subgraph of H I and H H will tend to be the cooperative model of whole brain vertices.
Finally, the H I has a similar clustering coefficient with the H ALL , which is due to the strong low-frequency signal, overlapping other signals.
From the Fig. 7(b) the length of the shortest path of the network varies with different network structures. At the same time of modular communication in the brain, the vertices maintain a high degree of connectivity, proving that the brain network is an efficient and integrated whole. Therefore, as shown in Fig. 7(c), we can see that low, intermediate and high frequency networks and their combinations have small-world properties.
Scale-free analysis. Due to the discussion of scale-free properties of brain networks in existing researches, scale-free networks have serious heterogeneity. In addition, the connection condition (degree) of each vertex in scale-free network has serious uneven distribution: In the network, a few vertices, called Hub vertices, have vast connections, and most nodes only have a small number of connections. Therefore, in this paper, the scale-free properties of low, intermediate, high frequency band and full frequency band are statistically analyzed. As the results of all data sets are similar, only the 'Cambridge Buckner' data-set is used as an example for subsequent experiments and the results are shown in Fig. 8.
The statistics of the power rate characteristics in low, intermediate, high frequency and full frequency band network can be seen from the results. During the analysis of scale-free properties of the proposed network, the scale free curves shown in full frequency band are similar to those shown in low frequency band, but are not strict scale free curves. In the intermediate and high frequency, we can determine that the structure has scale-free characteristics, the low-frequency working state of the brain is more prominent, and that the intermediate and high frequency band also contains a lot of relevant information. Therefore, in the voxel level multi-graph model, the brain network exhibits some scale-free properties.
Similarly, we use the p-value to calculate the significant differences under different frequency bands. As shown in Table 4   The 'hub' vertices. Due to the voxel network structure defined in the study, the vertices cover all voxels in the whole brain. Therefore, it is necessary to analyze the distribution of important vertices (i.e. the 'hubs') 24,25 . In order to calculate the importance of vertices, the degree centrality is adopted. Degree centrality is the most direct measure of vertex centrality in network analysis. The greater the vertex degree of vertices, the higher the degree centrality of the vertex, therefore, the vertex is more important in the network. In actual research, the most commonly used vertices distribution is the prior brain automated anatomical labeling (AAL) distribution structure of 90 cerebral regions used by Salvador 7,26,27 . Therefore, after extracting the vertices with the highest degree centrality as the brain network 'hub' vertices in H L , H I , H H , we compare the distribution data with AAL as shown in Fig. 9.
In order to show the difference of vertices distribution more intuitively, the distance between the vertex i in the obtained network and the central nvertex k in the one of regions of AAL is expressed as d ik . Setting the threshold δ. If d ik < δ, The node i belongs to the cerebral region k. If d ik > δ, vertex i continue to compare with other cerebral regions. If there are i, j vertices, and d ik < δ, d jk < δ, then only one count is retained. After all vertices     Fig. 10. As seen from the Fig. 10, in the H L , H I and H H , the calculated vertices can cover all the anatomical cerebral regions. Among them, the cerebral region of hub vertex is the most dispersed in H L , which is dispersed in different cerebral regions, and the vertices in H I and H H are second. In different frequency bands, the change of the working model of the brain network leads to the difference of the position of the central 'hub' vertex. This shows that nodes tend to work globally in low frequency networks and work more via inter-regional communication in intermediate frequency and high frequency networks.
Vertices connection distribution. In order to investigate the connection patterns among vertices, the connectivity patterns of all vertices in the whole brain are statistically analyzed. First, all vertices are registered into 90 cerebral region of AAL 26 . All the vertices are clustered according to the cerebral region. Next, count the number of connection edges between the vertices in each cerebral region and the vertices in all other cerebral regions. The brightness of each lattice in the matrix reflects the number of relationships between the two cerebral regions, shown in Fig. 11(a). The cerebral region corresponding to the number is shown in Fig. 11(b), label 1-45 represent the left portion of the brain, and label 46-90 represent the right portion of the brain.  The vertices connection distribution matrixes are also test by p-value, as shown in Table 5. According to the table, the vertices connection distribution on full frequency band has no significant difference with the low frequency band by the p = 0.0679. and low-, intermediate-, and high frequency band have the significant differences as the p < 0.05.
As shown in Fig. 11(a)I, in the low frequency brain function network, the connectivity of the vertices in the cerebral region shows a uniform state. Previous studies have shown that cerebral regions show a higher degree of correlation at low frequencies 28 . Therefore, it can be said that in low frequency brain functional networks, the connections between different regions of the brain are mainly low-frequency connections, and this connections are similar to that of Fig. 11(a)IV. We can concluded that in the whole brain, low frequency connection is more significant in the remote cerebral region connection. But in the diagonal dominated region in Fig. 11(a)I, the connectivity is stronger. This situation is more obvious in Fig. 11(a)II,III. By comparing Fig. 11(a)I,IV with Fig. 11(a)III and Fig. 11(a)III, low frequency and full frequency band functional network show high similarity. Therefore, the brain network constructed with full frequency band covers a large amount of information in Fig. 11(a)III and Fig. 11(a)III.
By observing Fig. 11(a)II,III, in the intermediate frequency and high frequency brain functional networks, the connectivity of the brain network shows a higher local connectivity, which cannot be expressed in the full frequency band signal or low frequency signal. In Fig. 11(a)III, we can see that the connectivity structure of the cerebral region on the diagonal presents a state of high brightness. In the network, the connectivity of the brain is dominated by its own internal communications and the proximity of cerebral regions this can be seen from the top right of Fig. 11(a)III. The left cerebral region in the brain corresponds to the right cerebral region in a highly     connected state, which is due to the same or similar function of the left and right corresponding areas in the original AAL partition. This situation, as shown in Fig. 11(a)II, is also reflected in the network. From Fig. 11(a)II,III, we can see the left and right brain connections show a high intensity distribution on the diagonal on the lower left region of the image. It can be concluded that the functional connection between the cerebral region, the adjacent cerebral region and the left and right brain is mainly connected by intermediate and high frequency bands. Therefore, the correlation between voxels and the function of brain regions are taken into account, if only the full frequency band signal is used as the main signal of brain network analysis, then a considerable amount of information in the research of brain function network will be covered by low frequency signals. Therefore, when establishing the brain functional network model, the low frequency, intermediate frequency and high frequency information separation analysis can be used to create a larger, multi-pattern network model, which contributes significantly towards the study of brain function.
Multi-graph network based Alzheimer's disease classification. Brain network alterations in patients with Alzheimer's disease (AD) has been the subject of much research, but the biological mechanisms underlying these alterations remain poorly understood 29,30 . In this section, we verify whether the multi-graph model can be used for classification of these patients from healthy control (HC) subjects and whether we can achieve better classification accuracy than the no frequency division model.
In typical connectivity-networks-based classification approaches, local measures of connectivity networks are first extracted from each region-of-interest as network features, which are then concatenated into a vector for subsequent feature selection and classification. However, some useful structural information of network, especially global topological information, may be lost in this type of approaches. To address this issue, Jie et al. 31 proposed a connectivity-networks-based Alzheimer's classification framework that involves the use of a new recursive feature elimination method based on graph kernel to measure directly the topological similarity between connectivity networks and use kernel Support Vector Machine (k-SVM) to fuse all features from multiple thresholded networks for final classification (namely, RFE-GK).
Jie et al. 31 has already demonstrated the advantage of RFE-GK in the classification of Alzheimer's disease compared with other methods. And all the Alzheimer's disease classification methods (including RFE-GK) in Jie's work are still based on the non-frequency division network model (i.e. Full frequency band network). Therefore we only compared the proposed multi-graph model with the RFE-GK.
For the multi-graph model, we also use the method of graph kernel construction proposed by Jie et al. 31 . However, the difference is that the graph kernel we built is based on the multi-graph model. If there are n subjects (including AD and HC), for the subgraph χ i m of the subject i under the frequency band m, the graph kernel is expressed as The accuracy of AD classification is verified by data obtained from ADNI (Alzheimer's Disease Neuroimaging Initiative, http://adni.loni.usc.edu/) database. Ethical statements are present by ADNI, and we confirmed that all experiments were performed in accordance with relevant guidelines and regulations of the access. We also confirmed that all experimental protocols were approved by Sino-Dutch Biomedical & Information Engineering School, Northeastern University, China.
A total of 184 rs-fMRI data subjects were obtained, including 92 AD patients and 92 HC subjects. The details of the data are shown in Table 6.
In order to verify the superiority of the proposed model and avoid the influence of classifiers, the classification performance was compared with RFE-GK by full frequency band network by the classifiers kernel Extreme Learning Machine (k-ELM) 32 and kernel Support Vector Machine (k-SVM) 33 . The result is shown in the Table 7.
It can be seen from the table that the classification performance of the proposed multi-graph models on the two classifiers is much better than the contrast model. The proposed model has 82.5% classification accuracy on the k-ELM and 77.5% classification accuracy on the k-SVM, while the contrast model has only 70% classification accuracy on the k-ELM and 67.5% classification precision on the k-SVM. It also can be seen that the classification effect of the two methods on k-ELM is better than that of k-SVM classification. The results show that the proposed multi-graph method can make full use of the hierarchical structure of the brain functional network, thus enabling better classifications of AD patients.  Second, the ROC curves are compared on k-ELM with better classification effect, shown in Fig. 12.
The superiority of the multi-graph method can be observed on the ROC curve of the above image, and the comparison between full frequency network model which has 0.68 AUC values and the multi-graph model having 0.81 AUC value through the experiment. It also shows that the brain network can be constructed from multi frequency bands angle using the multi-graph model, and the different frequency bands of brain activity are fully considered. The frequency domain information preserves the local topology of the brain network under different frequency bands, which fully integrates the brain activity information of the multi frequency band, and enhances the accuracy of the computer-aided diagnosis of Alzheimer's disease. Thus, compared with the original full frequency model classification method, the proposed multi-graph model can be helpful to enhance the classification accuracy.

Conclusion
A brain functional network model based on multi-graph is proposed in this paper. With the advantage of multi-graph, multi-band voxelwise brain network is constructed. As the possibility of the relationship between nodes can be described by multi-graph, a network model that can clearly describe the composition of brain function network is established. In this model, the components of the low frequency, intermediate frequency and high frequency bands in brain network are separated and analyzed, and some properties of them in the analysis of the brain network model are discussed. According to the results of network analysis, we can see that multi-graph model has important significance for the analysis and research of brain functional network. In our study, it is found that the distribution patterns of low-, intermediate-and high-frequency brain functional networks are of great significance to the study of the relationship between different cerebral regions.  Table 7. Classification performance of two models.