The EEG data and diagnosis details were obtained from Zhongshan Bo’ai Hospital’s internal medical record system, located in Zhongshan City, China. The inclusion criteria for the ASD patient group were: (1) diagnosed as ASD according to Diagnostic and Statistical Manual of Mental Disorders-V (DSM-V) criteria; and (2) aged 1–6 years. Exclusion criteria for the ASD patients were: (1) a history of epilepsy seizures; and (2) patients who had suffered from other neurological, developmental diseases, or extra-cranial injury. Controls were selected from children who were hospitalized for non-neurological disease, e.g., diarrhea, in the same age range. Altogether, 20 ASD patients and 25 health controls were recruited. Informed consent was obtained from the parents of all subjects involved in the study.

Nicolet EEG machines (Host: PN100973M ; Amplifier: V32-09472364, Natus Medical Incorporated, Pleasanton, CA, USA) with a 125-Hz sample rate were used to record the resting EEG signals of the subjects when they were sleeping. According to the international 10–20 montage system [28], a total of 19 channels (Fp1, Fp2, F7, F3, Fz, F4, F8, T3, C3, Cz, C4, T4, T5, P3, Pz, P4, T6, O1, and O2) were fastened onto the patient’s scalp [29]. The following were the EEG pre-processing steps: (1) the average re-reference and band-pass filtering between 0.1 Hz and 45 Hz were carried out for EEG signals, using the “pop_eegfiltnew” function in EEGLab toolbox v.14.1.2.0 (https://sccn.ucsd.edu/eeglab/index.php) [30] to remove artifacts; (2) independent component analysis (ICA) was carried out to remove artifacts originating from eye and muscle movement. The IClabel plug-in unit (https://labeling.ucsd.edu/auth/login) [31] carried out the identification of these components; (3) Manual elimination of the problematic segments from the eye, muscle, or body’s gross movement, that ICA was unable to identify; (4) recheck, by two experts, of the pre-processed EEG data, and extraction of the 100-s signal from the beginning of the pre-processed EEG signal; and (5) EEGs were further filtered into 5 frequency sub-bands, Delta (0.5–4 Hz), Theta (4–8 Hz), Alpha-1 (8–10 Hz), Alpha-2 (10–12 Hz), Beta (12–30 Hz) by the “pop_eegfiltnew” function in EEGLab. The pre-processing described above was carried out in Matlab v. 2019a using the EEGLab toolbox v.14.1.2.0 (https://sccn.ucsd.edu/eeglab/index.php).

In order to investigate the linear relationship between two signals in the fixed frequency, we computed magnitude-squared coherence, or coherence. By assuming that Y (t) and X (t) represent EEG signals in two brain regions or regions of interest (ROIs), first Y (t) and X (t) were converted from the time domain to the frequency domain using a fast Fourier transform:

(1)$$F(f)=\sum _{t=0}^{T-1}x[t]e^{-i2\pi tf}$$

Then, for each frequency point f, power spectral densities S${}_{x\mathit{}x}$ (f) and S${}_{y\mathit{}y}$ (f) and their cross power spectral density S${}_{x\mathit{}y}$ (f) were estimated; finally, the coherence function K${}_{x\mathit{}y}$ (f) at that frequency point f was further calculated as follows:

(2)$$K_{XY}(f)=\frac{S_{xy}(f)}{\sqrt{S_{xx}(f)S_{yy}(f)}}$$

A connection matrix could be formed once all of the functional connectivity had been computed for each pair of ROIs. It is important to apply a threshold, in order to eliminate low-connectivity edges that are driven by noise. A high threshold would produce an over-sparse connection network, whereas a low threshold would construct a regular lattice network. Determining the appropriate threshold is an arbitrary and contentious task. In light of this, we selected an intermediate criterion of 0.8, which preserved the top 20% of the graph’s edges, in accordance with other research [32], to preserve this delicate balance.

As in previous research [33], node strength was used in the present study as a static graph-theory metric for each ROI. The node strength, which represents a node’s overall capacity to convey information throughout the network, characterizes a node’s connectivity strength with all other nodes. A brain area that exhibits a high node strength is considered significant in the brain network and may be regarded as a hub of connectivity. This brain network would sustain more harm if that hub were destroyed. The examination of anomalies in the connections and interactions among different brain regions were facilitated by this indicator [34]. The following is the formula used to obtain the node strength for node I in the graph G:

(3)$$S_{\mathit{\text{nodal}}}(i)=\sum _{i\ne j\in G}{w}_{ij}$$

Where the w${}_{i\mathit{}j}$ represents the weight of the edge linking the node j to node i.

Each subject could obtain 5 (frequency sub-band number) $\times$ 19 (EEG channel number) = 95 features in node strength, for modeling.

All the node strengths were calculated using the “igraph” package v.1.3.5 [35] in R (v.2.3.2, Foundation for Statistical Computing, Vienna, Austria).

A sequence of functional-connectivity matrices can be obtained over time by using a sliding window, which can then be used to compute further dynamic graph-theory indicators. The movement step was set to be 50% or 100% of the window width, and the window width was set as 1 s, 3 s, or 5 s. We described the dynamic features along the extracted windows using six types of dynamic indicators: mean, standard deviation, median, kurtosis, inter-quartile range, and skewness of the node strength [29]. Similar to the static node-strength indicator, each subject’s 5 (frequency sub-band number) $\times$ 19 (EEG channel number) = 95 characteristics could be generated for each of the six dynamic graph-theory indicators. Fig. 2 depicts the process of feature extraction for both static and dynamic indicators.

In the present study, support vector machine-recursive feature elimination (SVM-RFE) was used to conduct simultaneous feature reduction, ranking, and classification, because the number of features was greater than the number of subjects. The objective was to rank the importance of features and extract the most important ones to form the optimal feature subset.

Because of its great efficiency, SVM-RFE is a popular feature-selection method that uses the support vector machine’s (SVM) classification performance as the standard [36]. SVM-RFE is a quick heuristic search approach that requires a straightforward process and relatively little computation. The feature-ranking and selection process happens nearly concurrently with model creation.

The foundation of SVM-RFE is the SVM classification algorithm, which looks for a hyperplane that fully divides samples into positive and negative classes. The sample set can be expressed as follows: The m = sample size, xi = a d-dimension vector of features, and y = the label of that sample.

(4)$$D=\{(x_1,y_1),(x_2,y_2),\mathrm{\dots },(x_m,y_m)\},$$ $$x_i\in R^d,y_i\in \{+1,-1\}$$

The hyperplane w can be formulated as:

(5)$$w^T+b=0$$

(6)$$\{\begin{array}{c}{w}^T{x}_i+b\ge +1,y_i=+1\hfill \\ w^T{x}_i+b\le -1,y_i=-1\hfill \end{array}$$

Finally, for the jth feature, the importance of that feature can be calculated as the square of that coefficient in w, namely:

(7)$$c_j=w_j^2,j=1,2,\mathrm{\dots },d$$

After acquisition of each feature’s importance and rating, the SVM-RFE was used to determine which feature subset had the best classification accuracy. The following is the process for the SVM-RFE: (1) features were arranged according to the coefficient size of each feature, as determined by formulation (7); (2) the SVM model was established, preserving all of the features, and the average accuracy of the test set was determined using leave-one-out cross-validation (LOOCV); (3) the least important feature was eliminated, and the step (2) was repeated to obtain the average accuracy on the test set, until the program reached the end condition; (4) the feature subset that had the highest average accuracy in test set was extracted. The subset with the fewest features would be chosen if the accuracy of the multiple feature subsets was the same. (5) To further characterize the classification performance for the selected feature subset, further classification measures (steps 8–12, below) would be computed. Fig. 3 depicts the SVM-RFE workflow.

(8)$$\mathit{\text{Accuracy}}=\frac{TP+TN}{TP+TN+FP+FN}$$

(9)$$\mathit{\text{Sensitivity}}=\frac{TP}{TP+FN}$$

(10)$$\mathit{\text{Specificity}}=\frac{TN}{TN+FP}$$

(11)$$\mathit{\text{Presicion}}=\frac{TP}{TP+FP}$$

(12)$$\mathit{\text{Recall}}=\frac{TN}{TN+FN}$$

TP: True positive; TN: True negative; FP: False positive; FN: False negative.

We only included the top 30 significant features in the recurrence so as to minimize the amount of computation required. The least significant characteristic would be eliminated in each repetition until the subset contained just the top five significant features. SVM-RFE was run using self-edited code in Matlab v. 2019a (MathWorks, Natick, MA, USA). The SVM model was established using the “fitcsvm” function, and the parameter was left at its default value of C = 1.0.

The median age of the ASD group and the control group were 37.5 months and 32.0 months, respectively, and showed no difference ($\alpha$ = p $\le$ 0.05) according to the Wilcoxon rank-sum test (p = 0.192). The proportions of males were 20% and 28% for the ASD group and control group, respectively, which showed no difference by Fisher’s exact probability test (p = 0.730). Of 20 ASD children, 19 reported an Aberrant Behavior Checklist (ABC) total score; the median was 88 (83.50, 118.5).

Table 1 displays the performance metrics for the static node strength and its six derivative dynamic local graph-theory indicators, under window-width conditions of 1 s, 3 s, and 5 s, with a 50% moving step. Table 2 presents the results of a significance test conducted using the Wilcoxon signed rank test to compare the accuracy and number of selected features of average dynamic graph theory with its static counterpart. Fig. 4 depicts the changing process of accuracy and area under the curve (AUC) along with changing of the number of features, and the Receiver Operating Characteristic (ROC) curve for the optimal feature subset.

According to Tables 1,2 data, the average classification accuracy of dynamic graph-theory indicators, under window widths of 1 and 5 s, was 0.859 and 0.870, respectively, and was lower than that of the static equivalent, 0.933. When the width was 3 s, however, the average accuracy was 0.952, which was greater than the static counterpart (0.933). Among the indicators, the mean, median, and inter-quartile range of node strength even reached the accuracy and AUC equal to 1.000. The significance test for the average accuracy of these 6 dynamic indicators in 3-s windows indicated that there were no significant differences with their static equivalents (p = 0.525). This may have been because the kurtosis and skewness indicators showed comparably lower classification accuracy (0.889 and 0.844, respectively). But as we could see, in the selected optimal feature subset, 6 dynamic indicators required an average of 15.5 features in the 3-s window-width condition. This was the lowest number of features required, compared to the 1-s (22.5) and 5-s (17.3) window-width conditions, and it was also significantly lower than the number of features required in the static condition (22, p = 0.036).

Fig. 4 shows the classification performance for six different types of dynamic indicators, represented by broken lines of different colors. The black line represents the static node strength indicator, and the vertical line indicates the location where this indicator reached its maximum accuracy. Regardless of the number of characteristics (X-axis), we can observe that the accuracy of the red line (Mean), green line (STD), dark blue line (Median), and purple line (IQR) was typically higher than that of the black line under the 3-s window-width condition. However, the vertical lines of the six dynamic indicators were all positioned to the left of the black vertical line, indicating that these indicators could achieve the best accuracy with the fewest features. Fig. 4g–i depicts the ROC based on the optimal-feature subset selected, and we can see that the ROC of 6 dynamic graph-theory indicators were higher than the black ROC line in the 3-s window-width condition. In the 1-s and 5-s window-width conditions, we could not see the advantages for 6 dynamic graph-theory indicators, because the colored lines crossed the black line, and were not always higher than the black line.

The above analytical process was repeated under the condition of 100% moving step, which meant that there was no overlap between two adjacent windows. The results are displayed in Tables 3,4, and Fig. 5.

Table 3 shows that although all conditions were below the static condition accuracy (0.933), the average accuracy of dynamic graph-theory indicators under the 3-s window width remained the greatest (0.930) when compared to the 1-s (0.863) and 5-s (0.851) conditions. The average accuracy in 1-s window width was significantly lower than its static counterpart (p = 0.036), whereas the number of selected features was higher than that in static condition (p = 0.036) (Table 4).

Fig. 5 shows that, in most cases, the accuracy and AUC of the mean (red line) and kurtosis (orange line) were slightly higher than the static node strength (black line) under the 3-s window width condition. However, under the 1-s and 5-s window-width conditions, we were unable to discern any clear advantages of the dynamic indicators over their static counterparts.

Table 5 and Fig. 6 generalized the overall classification ability measures across all scenarios in order to compare the classification abilities of the six dynamic graph-theory indicators. Tables 5,6 show that the IQR indicator outperformed the others, with an average accuracy of 0.922. The median indicator was next, with an average accuracy of 0.907.

The greatest results were obtained with a window width of 3 s and a movement step of 50%. Therefore, we investigated further by examining feature rankings and the distribution of particular features on scalp regions and frequency sub-bands. Table 6 presents the feature rankings and weights for 6 indicators of dynamic graph theory, in this condition. Features in the Beta and Alpha2 frequency sub-bands held the bulk of the top places, indicating that features in the high-frequency sub-bands had a greater contribution to the classification of individuals with ASD. In terms of mean, median, and IQR (all three indicators achieving 100% categorization), Beta-F3, Beta-F3, and Beta-Fp1 were the most significant features.

The proportion distribution of chosen characteristics across 6 dynamic graph-theory indicators in 5 frequency sub-bands is provided in Table 7. The results indicated that, with an unweighted proportion of 28%, and a weighted proportion of 40%, the features from the Beta sub-band had the dominant position. Similar to this, we investigated the distribution of particular traits on 5 regions of the scalp: parietal (P3, Pz, P4, P7, P8), temporal (T4, T5), occipital (O1, and O2), central (C3, Cz, C4), frontal (Fp1, Fp2, F7, F3, Fz, F4, F8), and parietal (P3, Pz, P4, P7, P8). It was evident that the frontal-area features that were chosen had the highest weighted proportion (46.153%) and unweighted proportion (36.364%). We also determined the relative relevance for the scalp region, taking into account that each region had a different number of channels. This was done by dividing the channel proportion (k/19, where k is the number of channels in a region). For instance, the frontal region’s unweighted relative relevance was 0.987, when 0.364 was divided by 7/19. Table 7 shows that, even after accounting for channel percentage, the weighted relative relevance of characteristics from the frontal region remained the highest (1.253).

The results showed that dynamic graph-theory indicators performed best when the window width and moving step were adjusted to 3 s and 50%, respectively. Some of them actually approached 100% accuracy.

The window-width setting in dynamic functional-connectivity analysis is still debatable and mostly dependent on researchers’ personal preferences. Previous research on dynamic functional connectivity demonstrated that although an excessively wide window width would obscure the information pertaining to dynamic transient shifting, an excessively small window width would raise the possibility of adding artifacts and lowering the frequency resolution [24]. Accordingly, it appears that a satisfactory compromise between frequency resolution and the capacity to capture variability was struck in the present study, since the 3-s window width yielded the maximum accuracy. The performance of the 3-s window-width dynamic graph-theory indicators marginally declined when the moving step was set to 100%, but it was still better than the 1-s and 5-s settings, which may provide additional support for choosing 3-s as the window width.

Another benefit in terms of preventing over-fitting was that we found that the number of selected features was lowest in the 3-s window-width and 50% moving-step conditions. One suggestion has been that when there are more features than subjects, over-fitting is more likely to happen [37]. The experimental results of Guyon et al. [36] demonstrated that the SVM algorithm benefits from feature reduction, even if SVM uses regularization approaches to partially avoid the over-fitting problem. Reducing the number of characteristics in the model enhances its interpretability while also reducing the possibility of overfitting. A model with fewer features is likewise simpler, requires less storage space, takes less time to compute, and possibly improves the model’s accuracy, because it contains fewer false positives. Features with higher interpretability and lower computation costs can offer a better understanding of the relationship between the input and output features, which is advantageous in situations in which resource efficiency is crucial [38, 39]. In summary, the 3-s window-width and 50% moving-step conditions allowed for the best classification performance with the fewest features needed. This finding should aid future research that uses dynamic functional-connectivity analysis regarding window-width and moving-step selection.

The IQR, an indicator that characterizes the dispersion and variance of the data distribution, had the highest accuracy among the 6 dynamic graph-theory indicators. This suggests that the value representing the variability of node strength had a significant ability to differentiate ASD patients from the control group. Only dFCA showed a significant difference in the study [40], suggesting that the functional-connectivity variability described by dFCA may offer distinct temporal information. The variability of functional connectivity reflects the instability of the information-transfer process within and between regions, and the flexibility of reconfiguration of one region with remaining regions [41].

Consistent with earlier research, the distribution analysis of a subset of features revealed a clear advantage of the frontal area. The frontal region was commonly identified as the aberrant portion using EEG or fMRI data in earlier research. The medial prefrontal cortex (mPFC) was identified as a hub of dysconnectivity in a study [42] that used lag-phase synchronization to construct functional connectivity on the cortical level, using a standardized low-resolution brain electromagnetic-tomography method. The mPFC plays an important role in the social cognition process. Similarly, several previous investigations identified aberrant functional connectivity in the default-mode network (DMN) in ASD patients [43, 44, 45]. Higher resting-state fronto-parietal-network (FPN)-DMN dynamics were linked to lower cognitive flexibility, according to one article [46]. According to Chen et al. [47], frontal-temporal connections in ASD connectivity demonstrated larger signal fluctuations. That suggested that increased variability of connection between these two regions could hinder the processing of social and cognitive information. The above results suggested that the there is an increased variability of functional connectivity between frontal region and other regions in ASD patients. However, reports of the opposite outcome also surfaced. Findings of Chen et al. [43] indicated that mPFC-insula connection variability was lower in ASD patients. Chen et al. [43] noted that the insula served as a core region of the salience network (SN), which suppresses executive function in the absence of external input, by strengthening its functional connectivity to the mPFC to preserve an internally focused state. However, in order to be alert for any modifications from the outside, the SN would sometimes alter the network configurations. Accordingly, a decreased variability in mPFC-insula connections could explain why ASD patients exhibit more internally directed cognition and respond less to the external world [43]. In summary, our findings indicated that the dynamic graph-theory indicators of node strength of the frontal region had the greatest significance in the diagnosis of ASD, suggesting that there are abnormal interactions between the frontal and other brain regions. These results may also reflect decreased and increased variability of the functional connectivity of the frontal region with other regions simultaneously, as well as an overall abnormality in the variability of topological characteristics in the frontal region.

We found that the Beta sub-band was the most significant of 5 frequency sub-bands (unweighted proportion 28.283%; weighted proportion 40.739%), indicating that high-frequency sub-bands may contain more information about the abnormality in ASD. The most significant features for mean, median, and IQR (these 3 indicators achieved 100% classification) were Beta-F3, Beta-F3, and Beta-Fp1. Beta-frequency coherence is believed to be related to cognitive and attentional functions [48, 49], and typically occurs in frontal and central region. Boersma et al. [50] observed that children with ASD have connection decrease in 51 ROI pairs, mainly in the Beta frequency. Additionally, they found that patients with ASD had significantly lower clustering coefficients and whole-brain average connection strengths in the Beta band than did controls. Research has indicated that people with ASD are less likely to be able to recruit Beta-band synchronization in large-scale networks, which could be a factor in the cognitive impairments that are common in ASD [51]. When the aforementioned information was combined with the findings of our feature-selection study, it became apparent that the aberrant variability of functional connectivity in the Beta sub-band may have the greatest impact on ASD diagnosis.

In a previous study [52], scientists classified ASD using the SVM algorithm, based on variables derived from brain functional connectivity across various frequency sub-bands obtained from fMRI. The best accuracy, according to the results, was 0.792. The Slow-4 (0.027–0.073 Hz) sub-band was found to include the majority of the discriminative features, and they found abnormal connections between the default mode network, the fronto-parietal network, and the cingulo-opercular network. Other research [53] classified children with ASD using the SVM model and used EEG to create a functional-connectivity network by partial correlation. They obtained the maximum accuracy at 0.800. When it came to dFCA, Price and his colleagues [27] extracted features from static and dynamic functional connectivity alternately, to match the multi-kernel support vector machine (MK-SVM) algorithm, and discovered that the dynamic feature could obtain a greater accuracy (0.900). Nevertheless, these studies only used the original functional connectivity as features to fit the model, and the ranking of feature importance was not involved. In order to differentiate ASD patients from normals, one study [54] used an ensemble classification model in conjunction with the SVM-RFE, and discovered that the posterior cingulate gyrus and precuneus had the highest ranking in terms of connection. Other researchers [55] used 5 global graph-theory indicators to attain an accuracy of 0.958 in their research using graph theory for the diagnosis of ASD. Another study [21] used a feature priority ranking by computing the node degree to fit the SVM model. The default mode network demonstrated a comparatively high network degree and discriminative capacity, and they were able to get 0.958 accuracy. In contrast to the previous research, ours has the advantage of combining indicators from local graph theory, dynamic functional connectivity, and ranking of future importance, all at once. This allowed us to make use of the dynamics information from dFCA and rank the selected features to identify the most significant ones, offering a fresh viewpoint on ASD diagnosis.

There are several limitations to our investigation. First, because a 19-channel EEG was used, source reconstruction could not be performed to find the EEG signals at the brain level. In order to improve this, we advise that future research employs EEG recording devices with greater spatial resolution in order to gather proof of convergence between fMRI and EEG. Second, the cross-sectional-data type affected the possibility of clarifying the mechanisms underlying ASD. We anticipate that future research will incorporate longitudinal follow-up EEG to compile data from the development track, in order to improve classification performance and to offer an alternative viewpoint for evaluating the course of the disorder. Third, there is a need for external validation and an increase in sample size. We intend to perform our study on a bigger external-validation dataset in the future.