CNNs are widely utilized in analyzing fMRI imaging data for autism diagnosis. Typically, CNNs employ convolutional filters to detect structural features within magnetic resonance imaging (MRI) scans to classify whether the patient has the disease or not. However, this approach often falls short in capturing the unique spatial differences inherent to individual patients. To overcome this limitation, we use an enhanced residual convolutional neural network. This network is designed with a more extensive and deeper architecture compared with traditional CNNs, enabling it to discern spatial variances more effectively between regions, which vary significantly among patients. We have defined spatial features not in 2-d space but in 3-d voxel space, hence we utilize a 3D ResNet structure.

We have devised a specific residual operation, denoted as $H(\cdot )$, alongside the convolutional operation $F(\cdot )$. These operations process the input at row $i$, column $j$, and layer $l$, represented as $x_{ij}^l$.

(1)$$H(x)=F(x)+x.$$

The key innovation in the proposed method lies in the incorporation of a residual connection that maintains the gradient $F^{\prime }(x)+\mathrm{\hspace{0.25em}1}$ across all layers.

(2)$$F(x)={\mathrm{\Sigma }}_{a=0}^{m-1}{\mathrm{\Sigma }}_{b=0}^{m-1}{w}_{ab}^l{x}_{(i+a)(j+b)}^{l-1}.$$

This feature is particularly critical in addressing the issue of gradient vanishing, a common challenge in deep neural networks. It not only mitigates the gradient vanishing issue but also ensures minimal loss of dynamic FC inherent in convolutional networks. Simultaneously, it maximizes the extraction and analysis of spatial variations.

The spatial features of the brain extracted through the residual network are fed into a BiLSTM-based multi-head attention mechanism for mapping time series attributes. One of the primary objectives of the proposed method is to extract and enhance transient brain features appearing in fMRI images, which is why we employ an attention-based mechanism to map these time series attributes. We use multi-head attention instead of single head attention because the latter is often insufficient to fully extract the dependencies between temporal features in the complex spatial context of the brain. Multi-head attention allows us to model these dependencies more effectively.

When considering the features $x_1,x_2,\mathrm{\dots },x_k$ extracted from the residual CNN blocks, their embedding through a BiLSTM and multi-head attention mechanism can be described by referencing the formula provided. At first, the features $x_1,x_2,\mathrm{\dots },x_k$ of the CNN are input into the BiLSTM network, which processes these features to capture the temporal dependencies, producing a series of hidden states $h_1,h_2,\mathrm{\dots },h_N$. In the context of the multi-head attention mechanism, each pair of hidden states $h_i$ and $h_j$ is evaluated to compute the attention weights $a_{ij}$. This is achieved through the following softmax function:

(3)$$a_{ij}={\text{𝑠𝑜𝑓𝑡𝑚𝑎𝑥}}_j\left(h_{ij}\right)=\frac{\exp \left(\text{𝑠𝑐𝑜𝑟𝑒}(h_i,h_j\cdot W_v{h}_i)\right)}{{\mathrm{\Sigma }}_{n=1}^N\exp (h_i,h_n\cdot W_v{h}_i)}.$$

Here, $score(h_i,h_j)$ is a function that measures the compatibility of the hidden states $h_i$ and $h_j$, which is calculated as:

(4)$$\text{𝑠𝑐𝑜𝑟𝑒}(h_i,h_j)=f(W_q{h}_i,W_k{h}_j)$$

where $W_q$ and $W_k$ are weight matrices for the query and key in the attention mechanism, respectively.

Finally, the output representation $r_i$ for each hidden state $h_i$ is computed as a weighted sum of all hidden states, with the weights being the attention scores:

(5)$$r_i={\mathrm{\Sigma }}_{n=1}^N{a}_{in}\cdot h_n.$$

This process is repeated independently $K$ times, and the results are concatenated to form a single vector. The proposed multi-head attention’s integration of spatial and temporal aspects enables it to focus on specific time frames and effectively model the transient characteristics of diseases. This mechanism allows for the analysis of disease characteristics from multiple perspectives, providing a detailed and multi-faceted understanding of the unique features of the disease.

Park et al. [41] demonstrated that the extraction of spatio-temporal features while preserving the integrity of dynamic connectivity significantly enhances accuracy. Nevertheless, a more direct method to learn the dynamic connectivity across 39 brain regions could offer further improvements. As depicted in Fig. 1e, the method leverages the spatio-temporal model as a foundational, pre-trained model. We transform the output of the intermediate layer, which consists of time-series spatial feature maps, into a connectivity matrix that delineates the relationships between 39 distinct brain regions. Subsequently, this spatio-temporal feature map is merged with a connectivity adjacency matrix, incorporating attention weights between regions. The GNN is then employed to explicitly learn and articulate the connection information among these 39 brain regions.

Upon processing the fMRI data through the LSTM network, enhanced with an attention mechanism, the output consists of a series of feature representations, symbolized as $R=r_1,r_2,\mathrm{\dots },r_T$. Each element $r_t$ in this series represents the features for different brain regions at time $t$, which are then fed into the ‘ConnectivityMeasure’ module to compute the connectivity matrices.

The connectivity matrix, denoted as $C$, is derived from these feature representations. Specifically, for each pair of brain regions, denoted as $i$ and $j$, over the entire time series $T$, the connectivity metric is computed, resulting in a matrix $C$ of $N\times N$ dimensions, where $N$ is the number of brain regions (we use 39 regions). This can be formally represented as:

(6)$$C_{ij}=\mathit{\text{ConnectivityMeasure}}(r_i,r_j).$$

Here, $C_{ij}$ represents the connectivity value between regions $i$ and $j$. The ‘ConnectivityMeasure’ method allows us to apply various metrics such as correlation, partial correlation, tangent, or covariance, depending on analysis needs. This results in different forms of the matrix $C$, each encapsulating unique aspects of inter-regional brain dynamics:

(7)$$C_{metric}=\mathrm{metric}(F),$$ $$\text{where metric}\in [Corr,Pcorr,Tan,Cov].$$

This method leads to the creation of a comprehensive connectivity matrix $C$, which captures the dynamic and complex inter-regional relationships within the brain. Such a matrix not only illustrates the direct interactions between different brain areas but also provides deeper insights into the overall functional network architecture.

Next, we interpret the connectivity matrix $C$ obtained from fMRI data as the adjacency matrix $A$ of a graph, with each element $C_{ij}$ representing the strength of connection between brain regions $i$ and $j$. This matrix is crucial for understanding the complex interconnections within the brain. To facilitate graph-based analysis, we categorize the continuous correlation values in matrix $C$ into discrete edge labels based on quartiles, as follows:

(8)

This labeling method allows us to transform the continuous brain connectivity data into a format suitable for GNN analysis.

In the proposed method, we translate the fMRI-derived connectivity matrix $C$ into an adjacency matrix $A$for graph construction, with discrete edge labels derived from $C$’s quartile-based categorization. The edge-conditioned convolution (ECC) layer plays a crucial role in embedding this graph, capturing intricate connectivity patterns for ASD diagnosis. The ECC layer is defined by the following process:

Given the graph $G=(V,E)$ with nodes $V$ (representing brain regions) and edges $E$ (representing connections between regions), and an edge label function $L:E\to {\text{𝐑}}^s$ assigning labels to each edge, the ECC layer embeds the graph by computing a signal $X^l$ at each node. This signal is a function of the node’s features and its neighborhood structure, mathematically defined as:

(9)$$X^l(v)=\sigma \left({\mathrm{\Sigma }}_{u\in 𝒩(v)}{\mathrm{\Theta }}_{uv}^l\cdot X^{l-1}(u)\right).$$

Here, $X^l(v)$ is the feature vector of node $v$ at layer $l$, $N\left(v\right)$ denotes the set of neighbors of $v$, ${\mathrm{\Theta }}_{uv}^l$ is the edge-conditioned weight matrix for the edge from $u$ to $v$ at layer $l$, and $\sigma$ is a non-linear activation function. The weight matrix ${\mathrm{\Theta }}_{uv}^l$ is computed as a function of the edge label $L(u,v)$, which reflects the strength and characteristics of the connectivity between regions $u$ and $v$:

(10)$${\mathrm{\Theta }}_{uv}^l=f^l(L(u,v)).$$

where $f^l$ is a learnable function specific to layer $l$, typically implemented as a neural network.

Through successive ECC layers, the graph is embedded into a high-dimensional space that captures not only the features of individual brain regions but also the complex connectivity patterns among them. This embedded representation is then used for classification tasks, such as identifying ASD. The ECC’s ability to condition on edge labels enables the model to adaptively learn different patterns of brain connectivity, crucial for accurate disease diagnosis.

The final classification is performed by feeding the embedded graph representation through a series of fully connected layers, culminating in a softmax layer for disease categorization:

(11)$$Y=\text{softmax}(\text{FC}{(X)}^L).$$

where $Y$ is the output classification, $X^L$ is the feature representation from the final ECC layer, and FC denotes the fully connected layers. In the proposed method, we do not use the typical graph classification Readout. Since we have predefined that each graph consists of the same number of nodes (39 in total), we do not employ an integrated graph embedding method that accounts for this.

We employed the preprocessed ABIDE I and II datasets [44], accessible through a web-based platform. We preprocessed the dataset using the Configurable Pipeline for the Analysis of Connectomes (C-PAC) as our preprocessing pipeline. This version of the preprocessing steps includes slice timing and motion correction, intensity normalization (4D global mean = 1000), nuisance signal regression, band-pass filtering (0.01–0.1 Hz), and registration to the 3 mm Montreal Neurological Institute (MNI) standard template. The fMRI volumes were down sampled to 4 mm in the MNI space. During the training phase, the data was augmented through intensity normalization, random Gaussian noise injection, and affine transformations, along with sequence sampling. During inference, a sliding window inference technique was employed, where the model made predictions for a continuous sequence of 10 frames (window size of 10) from user data using a sliding window approach. Soft voting was then conducted using the model’s predictions across the sliding window to aggregate the final prediction results. Table 2 shows the data collection institutions for ABIDE I and II, along with the distributions of ASD/normal and male/female, and the minimum and maximum ages.

Given the variability in the size and length of each fMRI image, we standardized all images to a uniform size of 24 $\times$ 24 $\times$ 24 and applied the sliding window technique to acquire serial images of fixed length. To ensure no overlap in patient information between training and testing, the dataset was partitioned based on individual patients rather than the images themselves. Consequently, the dataset was split into a training set comprising 90% of the cohort and a validation set consisting of the remaining 10%. The validation set was utilized for visualization and various analyses. Detailed information about the dataset used and the details of the proposed model are shown in Tables 3,4.

The results of the 10-fold cross-validation experiments in ABIDE I and II demonstrated the superiority of the proposed method. As shown in Fig. 2 and Table 5 (Ref. [29, 30, 31, 32, 41, 44, 45, 46, 47]) when compared with the graph-based approach Covariance FC with GCN, the method exhibited a performance improvement of approximately 3.7%p in the ABIDE I dataset and about 2.1%p in the ABIDE II dataset. This performance difference is attributed to two main reasons. Firstly, the use of covariance in FC for graph construction in traditional methods may not sufficiently capture the neurological characteristics of the brain. The proposed method, applying multi-head attention, diversifies perspectives in extracting correlations, thereby enhancing the ability to identify and emphasize heterogeneous data features from fMRI data. Secondly, unlike GCN which uses only node features as input, the method benefits from the Edge conditioned layer that models the correlations between nodes, leading to more effective embedding of the constructed FC and subsequently improved performance (detailed in section 4.3). Due to these characteristics, the proposed method outperformed existing 3D CNN neural network-based approaches that achieved state-of-the-art (SOTA) results through additional parameters or enhanced learning algorithms.

The receiver operating characteristic (ROC) curves of the five models are shown in Fig. 3. In contrast to alternative network models, our proposed model demonstrates a more balanced trade-off between accuracy and recall rates. Additionally, the aria under the ROC curve (AUC) index of our model surpasses that of the others, affirming its capacity to capture the intrinsic features of the data. This underscores the superior generalization ability of our model.

Fig. 4 illustrates the brain regions activated as identified through the proposed method’s model internal attention values, based on the presence or absence of various diseases, age, and sex derived from data and additional analyses. Notably, in the activation maps of Fig. 4a–c, where ASD is present, each map shows different regional activations, highlighting that the areas highlighted by the proposed method vary with age and sex. Furthermore, the activation maps of Fig. 4c,d for similar age groups and the same sex show generally similar patterns of highlighted regions, despite the presence or absence of disease. This underlines the necessity for models diagnosing ASD to capture not only specific regions and times of brain activity but to also comprehend the extent of each disease present in the dataset and the unique characteristics of each patient. In essence, this disorder exhibits heterogeneity based on disease progression, sex, and age, indicating that various complex external factors must be considered in diagnosing ASD. Based on our analysis of Fig. 4 and the detailed observations of brain region activations through the proposed method’s attention values, it becomes evident that our approach excels at extracting the nuanced characteristics of ASD and its variability across different ages and sexes. This proficiency is particularly crucial given the heterogeneity of ASD, which manifests distinctly in each individual based on a complex interplay of factors including disease progression, age, and sex.

Fig. 5 illustrates the analysis of the correlation matrix extracted by the proposed method. As observed in this correlation matrix, the time-series features extracted from ASD patients and processed through the attention mechanism closely connect various regions of the brain, reflecting their interdependencies and enabling effective prediction of ASD. In contrast, in a normal brain, the attention mechanism does not establish as strong connections between regions, indicating that the model fails to identify significant characteristics that could highlight connections across different brain areas. This distinction highlights the strength of the proposed method, which leverages attention-based feature extraction, graph transformation, and a GNN. These steps enhance the differentiation of disease-specific features, thereby improving the accuracy of disease diagnosis. The method’s effectiveness is rooted in its ability to accentuate the contrast between the neurobiological patterns of ASD and those of a typical brain, illustrating the nuanced complexities captured by the proposed method. The integration of the attention mechanism with graph-based analyses provides a more comprehensive understanding of brain connectivity patterns, crucial for discerning the subtle yet significant differences in neurological conditions. This advanced methodology demonstrates a significant stride in enhancing diagnostic precision through sophisticated analysis of fMRI data.

Fig. 6 depicts FC overlaid on brain images to analyze differences based on the presence of disease, as well as sex and age. Below each image, details such as the collection site, subject ID, sex, age, and disease status are documented. This allows for a detailed examination of patterns in ASD across different age groups, as observed in the differences between images Fig. 6a,b. This not only suggests that FC characteristics in ASD can vary with age and between individuals but also highlights the importance of deep learning models in capturing these individual variations. Furthermore, a comparison between images Fig. 6a,c demonstrates sex-based variations in brain characteristics. For instance, the posterior region of the left hemisphere is almost inactive in image Fig. 6c compared with Fig. 6a, highlighting the diversity in brain activation patterns in ASD. The differences among images Fig. 6a–d, where no similar patterns are observed, underscore the necessity of the proposed model. Additionally, the similarities observed in images Fig. 6d–f, regardless of the presence or absence of disease, present a unique challenge. Considering that d represents ASD while Fig. 6e,f are from healthy brains, distinguishing these subtle differences becomes a primary function of the model. Therefore, the effectiveness of our proposed method, which employs enhanced FC with multi-head attention and an edge label-based GNN, is validated by its superior performance.

Table 6 shows the statistical significance of various methods in comparison with the proposed method across a range of experimental settings. Each row represents a different method, while the columns delineate the p-values obtained through Wilcoxon tests and their corresponding false discovery rate (FDR)-corrected values. These metrics serve as indicators of the likelihood that the observed differences between the proposed method and its counterparts are due to chance. Upon examining the table, it is evident that the p-values associated with the proposed method consistently exhibit a substantial disparity compared with other techniques. For instance, the p-values for CNN Ensemble, Covariance FC with GCN, ASD-SAEnet, Vision-Transformer, 3D-CNN-LSTM, ASD-DiagNet, Contrastive CNN, and CNN-LSTM-GCN are 3.5 $\times$ 10${}^{-9}$, 6.58 $\times$ 10${}^{-8}$, 2.2 $\times$ 10${}^{-6}$, 1.3 $\times$ 10${}^{-5}$, 3.4 $\times$ 10${}^{-5}$, 4.1 $\times$ 10${}^{-5}$, 1.1 $\times$ 10${}^{-3}$, and 8.4 $\times$ 10${}^{-3}$, respectively, with FDR-corrected values in a similar range. This stark contrast suggests that the differences observed between the proposed method and each alternative are highly unlikely to have occurred by random chance alone. Furthermore, the FDR-corrected p-values reinforce the robustness of these findings, indicating that even after adjusting for multiple comparisons, the statistical significance of the disparities remains pronounced. Thus, based on these analyses, it can be confidently concluded that the proposed method exhibits statistically significant differences compared with the other methodologies evaluated in the study. Such findings underscore the efficacy and potential superiority of the proposed approach in the context of the addressed problem domain.

Our experimental analysis confirms that the proposed method significantly outperforms existing models in detecting ASD, showcasing its potential to improve early diagnosis and contribute to personalized healthcare strategies. The attention-based feature extraction mechanism effectively identifies the most relevant features across diverse patient profiles, illustrating a deep understanding of the disorder’s complexity. This capability is a testament to the model’s design, which integrates advanced machine learning techniques to navigate the intricate patterns of brain activity unique to ASD.

In conclusion, the proposed method introduces a more complex yet significantly more effective approach to ASD diagnosis. Despite the increase in complexity, this method marks a substantial advancement by offering unparalleled precision in understanding and diagnosing the disorder. Its sophisticated use of advanced techniques, including but not limited to multi-head attention, allows for the precise extraction of features from brain activation maps, enhancing diagnostic accuracy and our understanding of ASD. This method’s ability to navigate the disorder’s complexities provides clear advantages, including the potential for more accurate diagnoses, the development of personalized treatment plans, and the facilitation of further research into ASD. By leveraging these advanced techniques, the proposed method sets a new standard for ASD diagnosis, promising improved outcomes for individuals with ASD and paving the way for more effective, personalized care.

Although our method has shown promising results in effectively extracting brain networks and has undergone rigorous performance and qualitative evaluations for disease identification, we acknowledge the existence of multiple challenges that impede our progress towards a truly comprehensive diagnostic tool. In practical scenarios, it is highly unlikely to encounter situations where only a single disease needs detection. This reality encourages us towards the ambitious goal of developing a universal diagnostic model that can seamlessly transition between diagnosing various neurological conditions, such as Alzheimer’s disease, attention deficit/hyperactivity disorder (ADHD), and beyond.

However, a significant hurdle that currently limits our research is the focused nature of existing studies on particular diseases. This specificity necessitates a unique preprocessing routine for each disease dataset, complicating the development of a universal model. The need for a simplified and standardized preprocessing procedure cannot be overstated, as it would significantly facilitate the adaptation of models to handle diverse datasets efficiently.

Moreover, the advent of transformer-based feature extractors, renowned for their efficacy in numerous fields, particularly due to their strengths in transfer learning, presents an opportunity for innovation in fMRI data analysis. These models’ ability to adapt and learn from vast amounts of data could revolutionize how we approach fMRI data, making them an excellent candidate for future research. Nonetheless, developing a transformer-based model that caters specifically to the unique characteristics of fMRI data, such as its high dimensionality and the subtle nature of signal variations indicative of different neurological conditions, poses an additional layer of complexity. Such a model would need to not only extract relevant features from the fMRI data effectively but also be capable of doing so in a way that is computationally feasible and clinically relevant.

To evolve into a universal model that is adaptable to various downstream tasks derived from fMRI data, we must also consider the integration of multimodal data sources, including genetic information, clinical assessments, and behavioral data. This integration would undoubtedly enhance the model’s diagnostic capabilities but also introduces the challenge of managing and analyzing heterogeneous data types.