{'fraction_non_alphanumeric': 0.0359746835443038, 'fraction_numerical': 0.01139240506329114, 'mean_word_length': 5.007756653992395, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hierarchical Graph Convolutional Network Built by Multiscale Atlases for Brain Disorder Diagnosis Using Functional Connectivity

Functional connectivity network (FCN) data from functional magnetic resonance imaging (fMRI) is increasingly used for the diagnoses of brain disorders. However, state-of-the-art studies used to build the FCN using a single brain parcellation atlas at a certain spatial scale, which largely neglected functional interactions across different spatial scales in hierarchical manners. In this study, we propose a novel framework to perform multiscale FCN analysis for brain disorder diagnosis. We first use a set of well-defined multiscale atlases to compute multiscale FCNs. Then, we utilize biologically meaningful brain hierarchical relationships among the regions in multiscale atlases to perform nodal pooling across multiple spatial scales, namely"Atlas-guided Pooling". Accordingly, we propose a Multiscale-Atlases-based Hierarchical Graph Convolutional Network (MAHGCN), built on the stacked layers of graph convolution and the atlas-guided pooling, for a comprehensive extraction of diagnostic information from multiscale FCNs. Experiments on neuroimaging data from 1792 subjects demonstrate the effectiveness of our proposed method in the diagnoses of Alzheimer's disease (AD), the prodromal stage of AD (i.e., mild cognitive impairment [MCI]), as well as autism spectrum disorder (ASD), with accuracy of 88.9%, 78.6%, and 72.7% respectively. All results show significant advantages of our proposed method over other competing methods. This study not only demonstrates the feasibility of brain disorder diagnosis using resting-state fMRI empowered by deep learning, but also highlights that the functional interactions in the multiscale brain hierarchy are worth being explored and integrated into deep learning network architectures for better understanding the neuropathology of brain disorders.


I. INTRODUCTION
RAIN disorders, such as Alzheimer's disease (AD) and autism spectrum disorder (ASD) propose a severe challenge to the global public health. The AD could induce first mild cognitive impairment (MCI) and eventually dementia in older patients. And, ASD typically affects children's social interaction and communication abilities. As of 2020, 55 million people are living with dementia, and over 75 million people are affected by ASD world-wide according to the statistics of the World Health Organization (https://www.who.int/news-room/fact-sheets). An accurate identification as early as possible is necessary to enable effective interventions, which improve patients' outcomes. However, there is still a lack of clear understanding of the pathology of these brain disorders and well-established methods for the diagnosis.
Recent work focuses on diagnosing the human brain disorders by imaging the brain dynamics with in vivo brain functional imaging techniques. Among many of these techniques, resting-state functional magnetic resonance imaging (rs-fMRI) is a popular neuroimaging method that can probe spontaneous neural activities by measuring fluctuations in blood-oxygen-level-dependent (BOLD) signals. In the rs-fMRI studies, functional connectivity (FC) is the most prevalent measurement to quantitatively characterize spatiotemporal covariation patterns in the brain. Typical FC readings measure the temporal synchrony of neural activity across spatially distant regions by using pairwise Pearson correlations. The FC characterizations among all pairs of regions further lead to a whole-brain FC network (FCN). Many other FCN-constructing methods, such as the 3D maps generated by seed-based correlations or independent component analyses (ICA), can also detect neural synchronizations. FCN was found to be sensitive to capture functional changes induced by brain disorders, and has been used in individualized diagnosis for various brain disorders (see systematic reviews in [1]).
In FCN-based diagnostic studies, a long-lasting debate is about how to efficiently extract disorder-related features from the FCN. Conventional studies typically reshaped FCNs into vectors and then inputted them to classical machine-learning based classifiers, such as a support vector machine (SVM) [2]- [4]. However, such vectorizations may lose the disorder-related information carried on the topological structures of FCNs. Attributing to a rapid development of the deep learning technology, graph convolution neural networks (GCNs) allow implementations of artificial neural network algorithms on graph structures, thus providing suitable representations for the disorder-related information of FCNs. Briefly, a typical GCN deals with individual FCN and integrates nodal features along with the graph-topology neighbors based on the graph Laplacian and generates the representation for the entire FCN, which can be further used in identifications of brain disorders. Notably, another stream to apply GCN is to regard each individual data as a node and inter-individual similarity as a link. They perform nodal classification in the population-based graph to achieve individual diagnosis, which we call as "nodal classification GCN (nc-GCN)". The GCN-based brain disorder diagnosis has earned great attention, in diagnosing ASD [5], [6], major depressive disorder [7], [8], MCI, and AD [9]- [11].
However, so far, most of the existing studies, including those with GCNs, arguably ignored crucial multiscale architectures of the brain functional networks [12], [13]. Actually, in the brain, the nodes at multiple spatial scales are often organized in a hierarchical modular structure, e.g., individual neurons often connect to each other to form a "neural mass", many of which organize as neural circuits at a larger scale; and then neural circuits of many functional subsystems further form large-scale functional networks through inter-subnetwork connectivities [14], [15]. Such multiscale hierarchical brain connectomes could provide a much richer knowledge about the high-level brain functioning and therefore become more meaningful in characterizing brain disorders. First, it has been found that healthy brain connectomes at different scales exhibit different topological attributes, such as degree distribution [16], centrality [13] and small-worldness [17]. Second, brain disorders could modify FCN topology at a specific or multiple scale(s). A recent study specifically investigated modular structures of the lifespan FCNs at different spatial scales [18], indicating that the degree of segregation among network modules exhibits an interaction between age and spatial scales. Third, the relationship or interaction "across" scales could also be affected by atypical brain states, as revealed by a recent study on "cross-scale" functional interactions (namely "associated high-order FC") under the brain disorders [19].
Indeed, the latest studies started to focus on mining multiscale information from FCNs in data-driven manners [20], [21]. The advanced GCN studies have attempted to extract multiscale hierarchical topological features from the given graph, by learning various pooling operations on the nodes. Ying et al. proposed a coarse-graining based pooling "DIFFPOOL" and built a hierarchical GCN [22]. It learns the "cluster assignments" for nodes to aggregate the nodes into a set of clusters and to update the connectivity as well. The coarsened nodes and connectivities are forwarded to the further repeated GCN and DIFFPOOL layers, which can thus estimate multiscale community structures of the graph. Based on the DIFFPOOL, Xing et al. generated multiscale graph representations from FCNs and detected the diagnostic brain region clusters in MCI and obsessive-compulsive disorder patients during automatic diagnoses [20]. Also, Gao et al. proposed a "gPool" operation in a "Graph-U-Net" to generate a hierarchical representation of a given graph [23]. gPOOL learns a hierarchical down-sampling scheme from a single-scale graph based on the "importance" of nodal features, learned via learnable linear projections. It thus gradually focuses on a small set of nodes and their corresponding sub-networks. Following gPOOL, "self-attention graph pooling (SAGPOOL)" is further proposed to refine the hierarchical down-sampling by learning the nodal importance via GCNs on the graph [24]. Li et al adopted the idea of "gPool" and applied it to build hierarchical GCN for ASD detection [21].
However, all the above-mentioned studies can only mine multiscale features from a single-scale FCN, built by a predefined parcellation scheme of the atlas. On the other hand, multiscale brain atlases have been increasingly studied, which naturally provides FCNs at different spatial scales and brain hierarchies with more solid neuroscience knowledge. Recently, Yao et al. proposed a mutual multi-scale triplet GCN (MMTGCN) for brain disorder classifications based on the multiscale atlas information and the inter-subject sample relationship [25]. In addition to the "triplet sampling" strategy to learn inter-subject relationship, they further proposed a "mutual learning strategy" to collaboratively train multiscale encoders to generate aligned and optimal representations of the individual brain networks.
Lastly, the advanced brain atlas study offers refined multiscale parcellation and the biologically-validated prior brain hierarchy as spatial relationship among multiscale atlases [26]. In this work, we propose a novel multiscale-atlases-based hierarchical GCN architecture, to respect and utilize the prior brain network hierarchy under the GCN framework and to more effectively and sensitively capture the disorder-related intraand inter-scale topological changes in FCNs. Specifically, we first implement multiscale atlases to construct multiscale FCNs for avoiding the potential bias in single-scale analysis and to comprehensively characterize the brain network. We also invent an Atlas-guided Pooling (AP) to perform nodal aggregation based on the prior hierarchical relationship defined in the multiscale atlases. Specifically, in each stage of nodal pooling, we update the corresponding graph with the empirical FCNs from multiscale atlases. We thus built a hierarchical architecture by stacking the GCNs and AP layers to extract multiscale disorder-related features for final diagnoses. We utilize three open-access datasets and perform diagnoses of AD, MCI, and ASD to demonstrate the effectiveness of our proposed method.
The main contributions of the presented study are three-fold. First, to the best of our knowledge, our method is the first attempt to embed the guiding priors from multiscale brain atlases into the hierarchical graph convolutional neural networks for brain disorder diagnoses. Our method distinguishes itself from other GCN methods as a novel and biologically-meaningful framework to naturally integrate domain knowledge into a data-driven learning-based method. Second, experiments demonstrate that our method can effectively detect more disorder-related features than other single-scale-based methods (including data-driven graph pooling) and other multiscale-atlases-based methods with classical fusion schemes, and thus promote precise brain disorder diagnoses. Third, our method contributes to the improvement of interpretability for artificial neural networks and also revealing of pathological factors of brain disorder among multiscale brain connectomes.
A part of this work was reported at MICCAI 2021 conference [27]. In this paper, we have added new variants of the proposed method by generalizing its pooling operations. We have also added new analyses to demonstrate feature interpretability. In addition, we have included more detailed descriptions of method implementations, systematic testing of effects of parameter changing, as well as extensive experiments on additional datasets to evaluate our proposed method.

A. Subjects and rs-fMRI data preprocessing
We use the rs-fMRI data from three datasets: i) Alzheimer's Disease Neuroimaging Initiative (ADNI-2 and ADNI-GO; http://adni.loni.usc.edu/) [28]; ii) Open Access Series of Imaging Studies (OASIS-3; https://www.oasis-brains.org/) [29]; and iii) Autism Brain Imaging Data Exchange (ABIDE-I; https://fcon_1000.projects.nitrc.org/indi/abide/) [30]. Table 1 gives the overall information of subjects in the three datasets. Chi-square test and two-sample t-test are used to confirm no significant discrepancies in gender and age across groups (i.e., p-value>0.05), while there are significant differences in MMSE between HC and AD subjects. Specifically, the ADNI dataset is used as the main dataset for illustrations, while the OASIS and ABIDE datasets serve as additional independent datasets to partially evaluate the generalizability of our proposed method.  For the ADNI dataset, a totally of 531 fMRI data from different subjects at baseline scans are selected (with no repeated scan for any single subject, which is same for the OASIS and ABIDE data), which contains 364 normal control (NC), 117 early MCI (eMCI), and 50 AD subjects. In ADNI, eMCI is defined to describe the subjects exhibiting cognitive decline but subtler than those at MCI, measured by the Wechsler Memory Scale neuropsychological test. Accordingly, the eMCI state can be more difficult to detect based on FCN. Each fMRI data is acquired with TR = 3s, TE = 30ms, Flip angle=80°, resolution=3.3×3.3×3.3mm 3 , 48 axial slices, and 420s in duration (140 volumes).
In the OASIS dataset, based on the availability, we include 207 NC and 43 AD subjects in the study. The fMRI acquisition protocols are TR=2.2s, TE=27ms, resolution=4×4×4mm 3 , 36 axial slices, and 372s in duration (169 volumes).
We applied well-accepted toolboxes, AFNI [31] (for ADNI) and DPARSF [32] (for ABIDE and OASIS), to perform a standardized preprocessing procedure for fMRI data. In particular, the first 5 or 10 volumes of each image are discarded due to potential non-equilibrium magnetization. The rigid-body transformation is performed to correct the subject's head motion. There is no large head motion (i.e., larger than 2mm or 2°) found in the used subjects. We do not further perform scrubbing/censoring of data as it may introduce additional artifacts [33]. The signals of white matter, cerebrospinal fluid, and head motion are regarded as nuisance covariates, and are regressed out from individual data. The fMRI images are then normalized to the Montreal Neurological Institute (MNI) space and spatially smoothed with a Gaussian kernel with full width at half maximum (FWHM) of 6×6×6mm 3 . The BOLD signals are further band-pass filtered (0.01 ≤ f ≤ 0.1 Hz) to remove the neural-irrelevant high-frequency noises and low-frequency drift from MRI. For ABIDE dataset, since volumes of scans are different among the collecting sites, we use its minimum common length, i.e., 115 volumes, from the middle of the preprocessed fMRI sequence for further processing.

B. Multiscale atlases for brain parcellation and FCN construction
Schafer et al. provide a set of atlases for multiscale brain parcellation [26], which is used in this paper for generating multiscale FCNs and guiding the node pooling across scales. The atlases are generated by clustering voxels (or vertices) based on FC patterns by considering both global similarity and spatial proximity. Using different resolution parameters during the clustering results in a set of brain functional parcellation at multiple scales, ranging from 100 to 1000 regions of interest (ROIs). These ROIs belong to seven resting-state functional networks (RSNs) defined in [34]. It can be observed that the RSN structures are largely preserved after parcellations at all scales (Fig. 1A). The relationship among the atlases at different scales thus characterizes a biologically-meaningful functional hierarchy.
Given an atlas at a specific spatial scale, the ROI-level signals can be extracted by averaging voxel-level BOLD signals within each ROI. The FCN at the given scale is then computed by Pearson correlation among all pairs of ROI-level signals and denoted as (Fig. 1B). To demonstrate our method and reduce computational and storage costs, in this study, we use the first five scales, i.e., from 100 to 500 ROIs, respectively (note that this method can be easily extended to other scales). Therefore, has a range from 100 to 500. The effect of increasing and decreasing the range of scales to be integrated is explored in later analyses, which can be extended to infer the situations when more scales are included.

C. Multiscale-Atlases-based Hierachical Graph Convolutional Network (MAHGCN) 1) Graph convolutional network
We implement the spectral graph convolution [35] to build GCN, which applies a series of approximation to obtain a simple formula of graph convolution operation and is the most popular version of the GCN. For a given graph with adjacency matrix with nodes, its corresponding symmetric normalized Laplacian is where is the identity matrix and is the degree matrix of , yielded by filling the diagonal elements with the degree (∑ ) of and the off-diagonal elements with 0. The eigenvectors of graph Laplacian defines the Fourier basis of the graph. With , the graph convolution operation on the inputted feature vector ℎ (with dimensions × ℎ , where ℎ refers to the dimension/channel size in the inputted node feature) can be defined in the spectral domain, where denotes the graph convolutional kernel with learnable parameters , and ℎ is the Fourier transform of feature vector ℎ.
Direct estimation of can be computationally expensive. According to [35], we apply a simplified operation where the outputted feature from GCN ℎ ( × ℎ ) is given by (3) ̃= + , and � is the corresponding degree matrix of ̃. is a ReLU activation function. The ℎ depends on the configuration of kernel in GCN. In this paper, ℎ is set to 1 for all built GCNs.

2) Atlas-guided pooling
We further define the "atlas-guided pooling (AP)" to bridge the GCNs at different scales and guide the fine-to-coarse processing of graph features. Multiple small ROIs in fine-scale parcellation could be merged into one ROI in coarse parcellation. The spatial overlapping of ROIs defined at different scales is thus utilized to approximately characterize the inter-scale hierarchy. We construct the mapping matrix between atlases by checking spatial overlapping between ROIs defined at different scales. Formally, from the atlas ℛ (with ROIs/nodes) to the atlas ℘ (with ROIs/nodes, where > ), the atlas mapping matrix ℛ→℘ (a × matrix) is defined as follows. First, we compute how many voxels of ROI in the finer-scale atlas ℛ are spatially overlapping with ROI in the coarser-scale atlas ℘, and further calculate the overlapping ratio between the number of overlapping voxels in ROI and the total number of voxels in ROI . Second, a threshold ℎ is applied to the overlapping ratio to define the inter-scale mapping matrix in the AP: The mapping matrix ℛ→℘ aims to pool the nodal features defined in the finer-scale atlas into the nodal features in the coarser-scale atlas. The GCN-outputted feature vector ℎ defined in the atlas ℛ can be mapped to ℎ for the atlas ℘ by a simple matrix multiplication with the transported mapping matrix, ( ℛ→℘ ) ℎ = ℎ .
Note that the above-defined multiplication-based mapping performs a sum pooling on the nodal feature across scales from the mathematical viewpoint. One can further define various pooling operations, such as average pooling and max pooling, with a slight modification.

3) Stacks of GCNs and APs for hierarchical architecture
After defining the AP, we can then build the multiscaleatlases-based hierarchical graph convolutional network (MAHGCN) with stacks of GCNs and APs in a hierarchical manner. In this paper, the MHAGCN integrating all five scales will be used as a full method for demonstrations. We extract FCNs from five scales (100, 200, …, 500 ROIs) at maximum, which can be used to build a five-layer MAHGCN (as illustrated in Fig. 1C. Different scale ranges are also tested in later experiments). In particular, MAHGCN starts with a graph convolution between the nodal features and the FCN from the finest parcellation (i.e., the 500-ROI atlas in the example of Fig.  1C). As the initial inputted nodal features, a one-hot feature is used here, i.e., a 500 × 500 identity matrix, rather than the features extracted from brain data. This is because the effect of choosing different nodal features is not in the scope of our study. And, such a null configuration (one-hot feature as used in [35]) will make the tested models here focus on topological information in FCNs and thus provide an extendable baseline. Using the AP, the outputted feature vector ℎ 500 (a 500 × 1 vector) is then converted to ℎ 400 (a 400×1 vector), which is further used in the next-layer graph convolution with 400 . The hierarchical graph convolution continues until FCN from the coarsest-scale atlas (i.e., the 100-ROI atlas) is processed. A skip connection design, inspired by [36], is also implemented to allow outputs from all the scales to be involved in the final prediction. The outputs from all the scales are concatenated to form a fused representation for the state of the multiscale brain system and sent to two fully-connected layers (FLs). In Fig. 1C, the fused representation is ℎ 1500 , with its size depending on the range of integrated scale (i.e., 1500×1 at most, and 300×1 at least). The first FL compresses the fused representation from MAHGCN into a 64×1 vector (which is fixed for variable ranges of integrated scale) and the second one converts the outputs from the first FL into the predicted probability for the two classes with a Softmax function (a 2×1 vector). During the implementation, the GCNs are associated with dropout (rate = 0.3), and the FLs are with the batch normalization and ReLU activation functions. These settings for GCN and FLs are kept in other methods for comparison.

D. Model explainability analysis
To examine the explainability of our proposed model, we apply the Gradient-guided Class Activation Map (Grad-CAM) algorithm [37] to MAHGCN. By computing the Grad-CAM from GCNs at different scales, the important brain regions, defined by multiscale atlases during the predictions, can be discovered.
Briefly, the Grad-CAM tracks the gradient of the prediction output for a certain class (i.e., the normal or disorder state) with respect to each element in the feature maps outputted from the intermediate layers in the neural network (i.e., the ℎ ). The gradient can be interpreted as the amount of influence of the intermediate feature elements on the outputted decisionmaking, and thus can be regarded as the importance of each feature element. The class activation map (CAM) is computed by the product between gradient map and feature map at individual data. Since the elements that positively contribute to the final prediction are more focused, the negative part of the individual CAM is further filtered out with a ReLU function.

A. Experimental settings 1) Implementation
The neural network based models are implemented using the open-source framework Pytorch [38] in Python (Version 3.6.2). The training of the network is accelerated by one Nvidia GTX 3080 GPU. To address the imbalance between positive and negative sample sizes, we apply weighted cross-entropy as a loss function, where the weights are set as the inverse of the sample ratio in the training set. Each model is trained with epoch=100, learning rate=0.001, and batch size=30. The time cost for training MAHGCN is around 30 minutes for ADNI and OASIS while 70 minutes for ABIDE. The parameters in each model are initialized with random weights. The Adam [39] is used as an optimizer, with weight decay=0.01 for avoiding overfitting.

2) Validation scheme and evaluation metrics
Five times of holdout cross-validation are carried out to evaluate the models. The dataset is randomly assigned into the training and testing sets with 80-20 split for five times. And then each model is trained and tested in the five different trainingtesting assignments, respectively. The accuracy (ACC), sensitivity (SEN), specificity (SPE), and area under the receiver operating-characteristic curve (Area Under ROC, AUC) are chosen as performance metrics. Since the sample sizes of different classes are imbalanced in ADNI and OASIS datasets, AUC can be a more informative assessment for evaluating performance there. The above-mentioned performance metrics from each method are generally provided as the mean and standard deviation over five times of cross-validation. Where possible, we report the confidence interval (CI) for more comprehensive statistical information.

3) Statistical analysis
To demonstrate the significance of improvement by our proposed method, we conduct statistical comparisons for performance metrics from cross-validations, with both onesided paired t-tests and Wilcoxon signed-rank test to examine the significance. The Wilcoxon signed-rank test is a nonparametric test on the data median, and has no assumption on the normality of the data distribution. Both types of tests are performed based on the build-in functions in Matlab, i.e., ttest and signrank (R2020b, Mathworks Inc, USA).

B. Comparison methods
First, the single-scale GCN method is adopted to provide a baseline to be compared with multiscale atlases based methods.
• GCN: The GCN is based on [35] and the configurations described in the Method parts. The GCN is single-layered, which is followed by two FLs to generate a prediction of the status for each subject. We report the results from GCNs based on FCNs at 500-ROI scales for better a comparison with the following experiments.
Second, several representative data-driven methods are included for comparison including DiffPool [22], gPool [23], and SAGpool [24]. To make fair comparison, we mimic the structure of MHAGCN with five scales and let the DiffPool (DP), gPool (GP) and SAGPool (SAGP) mine the 500-ROI FCN and generate 400-, 300-, 200-, 100-ROI FCN in the same hierarchical manner. After layers of GCN and pooling, the encoded features are concatenated and fed into two FLs, where configurations are the same as the main architecture of MAHGCN as described above.
• DP: The DP is akin to a clustering method, which learns a nodal assignment/clustering matrix based on the features extracted from the nodal and graph features by a GCN. The assignment matrix serves as the mapping matrix ℛ→℘ in our method, and can be used to generate the clustered/coarse-grained graph for further processing with GCNs. Specifically, they further proposed an auxiliary link prediction objective and entropy regularization to ensure the learned assignment matrix to be sparse and modularized.  • GP and SAGP: The GP and SAGP estimated the coarsescale graph in a down-sampling manner. In GP, the nodal features are inputted into FLs and a Softmax operation to estimate the nodal importance, based on which the graph is down-sampled into the subgraph with the top K important nodes. Here, K is the parameter to control the sampling rate. In SAGP, the nodal importance is learned on the nodal and topographical features by a GCN. The down-sampled graph is forwarded to the following GCNs for processing.
Third, we compare multiscale-atlases-based methods with classical fusion schemes to demonstrate the effectiveness of multiscale atlases fusion and the advantage of our proposed method, MAHGCN. In particular, we perform • MV and WV: Decision level fusion is first implemented.
A Majority Voting (MV) among predictions from respectively trained GCN-based classifiers at different scales. The prediction from each classifier is regarded as one vote, and the final predicted label is the most voted label among all the classifiers. In addition, we also perform Weighted Voting (WV), in which the predicted probability from the pre-trained GCN-based classifiers is weighted-summed with the weights learned by training.
• MAGCN: Moreover, feature-level fusion is performed. We train a model with multiple GCNs to independently process the inputs of multiscale FCNs. The outputs from multiple GCNs are concatenated and forwarded to two FLs. We name this model as MAGCN. Note that this model can be also viewed as MAHGCN without AP, and thus can also serve as an ablation study to highlight the importance of including the prior of inter-scale hierarchy among GCNs at different scales.

A. Comparison study
In this experiment, MAHGCN is compared with three datadriven hierarchical graph encoders and three multiscale-atlasesbased methods. The models are trained and tested in ADNI, OASIS and ABIDE datasets, respectively. In Table 2, the performances of all comparison methods on HC-vs-eMCI, HCvs-AD and HC-vs-ASD classification tasks are depicted. In all classification tasks, the following consistent results can be observed. First, both hierarchical graph processing with datadriven pooling methods (DP, GP and SAGP) and the multiscale-fusion based methods outperform the single-scale based GCN, suggesting the advantage of multi-scale analyses. Second, all best performances are achieved by multiscaleatlases-based methods, rather than data-driven methods, which suggests the advantage in biological-prior based methods over data-driven methods in properly capturing multiscale brain features. Third, among all the comparison methods, our proposed method achieves the best ACC and AUC. Most improvements are showing significances. In addition, when compared with MAGCN, i.e., MAHGCN without AP, our MAHGCN provides significantly higher AUCs in all the four tasks. This demonstrates that our proposed MAHGCN based on AP can best integrate information of multiscale-atlases-based FCNs for brain disorder diagnosis.
Besides, the performance on HC-vs-AD classification by MAHGCN on the OASIS dataset is comparable to that yielded on ADNI dataset, which could partially support the generalizability of our proposed model on AD identification across different datasets. In addition, our method outperforms other methods in the HC-vs-ASD classification task using the ABIDE dataset in terms of ACC, SPE and AUC, suggesting that our proposed method is generalizable to different brain disorders other than MCI and AD. We further compare performances of our method with performances from state-of-the-art studies on FCN-based diagnostic analyses. In Table 3, we group the studies based on the classification tasks, and report the number of subjects and the achieved performances. Due to the variation of sample selection, the performances may not be quantitatively comparable. In the HC-vs-eMCI and HC-vs-AD classification tasks, our method obtains better results than the advanced CNN-based methods. And, in the HC-vs-ASD classification task, our methods achieve slightly worse but competitive results when compared to other GCN-based methods. The possible reason could be the usage of the larger dataset that brings more difficulty to the classification task.

B. Comparison with state-of-the-art methods
We further perform a comparison with the MMTGCN method [25] that also used the multiscale atlases. Particularly, in the ADNI dataset, we select 148 MCI subjects (63 males, age=74.2 ± 8.9), match them with 155 HCs (77 males, age=74.9±6.4), and apply our MAHGCN. Note that Ref. [25] also used ADNI dataset and the ratio between HC and MCI in their study is not significantly different from ours (Chi-square  [25] in the HC-vs-MCI classification. However, we emphasize that their idea and our idea are different and could be integrated together. Their idea relies more on the inter-sample relationship, which is different from our design. More specifically, MMTGCN introduces a "triplet sampling" strategy to consider the relationship among subjects at each scale and further aligns the distributions of the encoded sample among multiple scales by a "mutual learning strategy" to benefit the integration of multiscale information. Our method considers the inter-scale correlation in nodal features when hierarchically encoding the FCN from one subject. In addition, the discrepancy in the performance may also attribute to difference in the diagnostic tasks. Their strategy of using intersample relationship could be more beneficial for the MCI detection. Whereas, our work could be more specific to the diagnosis for the eMCI state, which is earlier than MCI and thus is more clinically significant. We have investigated the effect of changing the range of integrated scales, using HC-vs-eMCI and HC-vs-AD classification tasks. In Fig. 2, it can be also observed that the ACCs and AUCs of MAHGCN are generally increasing with more scales included. We also try to integrate non-neighboring scales such as including 100 and 500, as well as 100, 300 and 500 in the HC-vs-eMCI classification, and observe similar increasing trends. . Both observations suggest that MAHGCN can effectively integrate information from multiscale FCNs for diagnosis. However, the performance may not grow continuously with further increase of the range of integrated scales. We expect a ceiling effect could appear, although currently not being observed in our experiments. The artfacts in parcellation at fine scales could affect a meaningful hierarchical encoding. Also, more integrated scales require increased parameters in MAHGCN, leading to the risk of over-fitting on small-sample data. We perform further ablation studies by changing the configurations of the MAHGCN. The HC-vs-eMCI classification, as the most difficult task, is performed to sensitively assess the changes due to modification of the model. We test the performance under different pooling schemes (i.e., using max, average, and sum pooling) in the AP. The corresponding models are called Max AP, Average AP, and Sum AP. Note that the "Sum AP" is a default configuration we use above.

2) Different configurations of atlas-guided pooling
Interestingly, unlike natural image processing where the average pooling and the max pooling are the most efficient pooling schemes, we find that the sum pooling significantly outperforms other types of pooling in terms of ACC, SPE, and AUC (Fig. 3). This may attribute to the essential difference between natural image data and brain network data. Average pooling normalizes differences in a number of nodes under integration at the finer scales, and flattens differences of nodes at the coarser scales; while, the max pooling only relies on the ROI with max features to represent other ROIs at the finer scale.
These pooling schemes could ignore meaningful biological information for the brain system and thus produce decreased performance. In Fig. 4, we systematically test MAHGCNs with different schemes of AP and their corresponding versions without SCs. The SCs significantly improve AUC in MAHGCNs with all schemes of AP. This is reasonable since the final output from MAHGCN only provides an abstract feature shared by FCNs at all scales, while SCs allow the reservation of scale-specific information at all scales. In addition, the SC could overcome the gradient vanishing and benefit to the optimization of MAHGCN.

4) Effects of thresholds for atlas mapping matrix
In this paper, we simply use ℎ=0, under which any sharing between the ROIs in atlases is considered as overlapping ROIs. Other settings of ℎ are also tested, and the respective diagnosis performances are provided in Fig. 5. Results suggest that a tight threshold could diminish the hierarchical relationship among scales and degrade the performance.

5) Effects of pre-processing pipeline
The ABIDE dataset provides pre-processed fMRI data using different pipelines, which are further quality-controlled (see https://preprocessed-connectomes-project.org/abide/). To test the effect of pre-processing on our method, we utilize those preprocessed data under two pipelines, i.e., DPARSF and Connectome Computation System (CCS, based on AFNI and FSL [42]), to investigate the effect of pre-processed pipeline on our MHAGCN. After checking data availability for both pipelines and matching the age and gender, we include 371 HC (311 males and age=16.60±6.70 years) and 325 ASD (284 males and age=17.07±8.32 years) for this analysis.
The Wilcoxon signed-rank test (p=0.526). The reduction of performance compared to results in Table 1 may attribute to the changes of sample size.

D. Explainability analysis
In order to demonstrate the explainability of the prediction from MAHGCN, we apply Grad-CAM to explore important features during the decision making by MAHGCN. In particular, the group-wise CAMs of five models from crossvalidations are first obtained by averaging the CAMs from each individual within brain disorder groups. Then, the resulting group-wise CAMs are weight-averaged according to the AUC of each model, yielding the CAMs as brain patterns (Fig. 6, A-B). In Fig. 6, C-D, we further assign class activation values in ROIs to the corresponding RSNs and compute the RSN-wise averaged class activation values.
For AD diagnosis, a pronounced finding is that the limbic network (LIM) is detected as an important predictor shared by most of the scales (100, 300, 400 and 500). This is consistent with the widely accepted understanding that AD is associated with the damage in the hippocampus related to the limbic network and thus typically exhibits amnestic cognitive impairment [43]. In addition, the default mode network (DMN) (related to executive functions) at scales of 200, 400, and 500 ROIs are found to be the second most contributive RSN, which is also consistent with the findings in the AD literature [43]. These consistencies could partially support the explainability of the features identified by our method.
For eMCI diagnoses, different scales rely on different RSNs as the most important features. Somatomotor (SM), limbic, visual (VIS), attention (ATT) and salience (SAL) networks can be identified as the predictive RSNs at different scales, and there is few RSN being consistently identified as a crucial predictor. This may also be reasonable due to complex symptoms of the eMCI involving impairments in multiple cognitive domains, including executive function, attention, and visuospatial ability or language, besides memory [44], [45]. This could deliver a new insight that eMCI damages the functioning of different cognitive domains at different scales.

V. CONCLUSION
We have presented a novel method for brain disorder diagnosis based on information from multiscale FCNs and predefined multiscale atlases. Using neuroimaging data from ADNI, OASIS and ABIDE, we have demonstrated that our proposed method has a higher effectiveness in identifying brain functional deficits induced by MCI, AD, and ASD. The improvement of diagnosis by our proposed method comes from the usages of prior knowledge on brain parcellation as well as proper designs (based on brain hierarchy) to fuse multiscale spatial information in a deep-layer GCN structure. We have further investigated the features of atypical brain connectomes across different scales, and provide novel clues about the pathology of MCI and AD, paving a new path to explore multiscale-atlases-based FCN in the brain disorder diagnostic studies. Through FCN is focused in this study, we expect our method or its variations could be applied in brain network data based on other modalities [46]- [49] in the future.