EMD-based analysis of complexity with dissociated EEG amplitude and frequency information: a data-driven robust tool -for Autism diagnosis- compared to multi-scale entropy approach

Abbreviations: ASD: Autism spectrum disorder; ML: Machine learning; EEG: Electroencephalogram; AUC: Area under curve; PCA: Principal component analysis; SCG: Scaled conjugate gradient; GD: Conjugate direction; LM: Levenberg–Marquardt; CV: Cross validation; fMRI: Functional magnetic resonance imaging; MSE: Multi-scale entropy; EMD: Empirical Mode Decomposition; IMF: Intrinsic Mode Function; NN: Neural Network; DQ: Direct Quadrature; SE: Sample entropy; SF: Scale factor; CARS: Childhood autism rating scale; ADOS: Autism Diagnostic Observation Schedule; ATEC: Autism Treatment Evaluation Checklist; AM: Amplitude modulated; FM: Frequency modulated; sMSE: Short-time multi-scale entropy; sMRI: Structural magnetic resonance imaging; MEG: Magnetoencephalography; RMSE: Refined multi-scale entropy; RCMSE: Refined composite multi-scale entropy; DFA: Detrended fluctuation analysis; IFAST: Implicit function as squashing time algorithm; RQA: Recurrence quantitative analysis

1.   Introduction

Autism spectrum disorder (ASD) is usually characterised by altered social skills, repetitive behaviours, and difficulties in verbal/nonverbal communication. Early signs of autism are generally noticeable by the age of 2 or 3. Timely intervention has been shown to have positive effects on people with ASD. However, diagnosis can be challenging, as there is no medical test to diagnose ASD. The only available diagnosis method is assessment by an experienced professional through descriptive scaling based on individual's developmental history and behaviour ^[1,2,3,4].

Several approaches have been proposed in the literature with the aim of providing a reliable diagnosis technique for ASD. Among these procedures, machine learning (ML)-based methodologies have demonstrated promising results ^{[5,6,7,8,9,10,11,12]}. The common goal of these ML methods is to identify the best set of quantitative features that can lead to the correct classification of two distinct groups: neurotypical and ASD ^[6]. Hence, all of the explored features must be derived from the top discriminative characteristics (between ASD and other conditions), e.g., the power of specific frequency bands in electroencephalogram (EEG) or functional magnetic resonance imaging (fMRI), connectivity or coherence between brain sites, and brain asymmetry ^{[13,14,15,16,17,18]}.

Any technique implemented to detect ASD via EEG-extracted features should satisfy a number of requirements. First, it should be well adapted to non-stationary processes ^[19], which excludes all "classical" frequency domain methods such as the Fourier transform. Second, an oft-neglected requirement in many approaches is the amplitude information. In ^[20,21], it was shown that the critical information in ASD EEGs is found not only in the frequency data, but also in the electric potential amplitude because both contents comprise the "brain rate", which is the fundamental state characteristic. The amplitude of the EEG reflects the general intensity of cortical neural activities, providing a good indication of the cerebral state ^[22]. Third, it was shown in ^[23] that linear methods such as Fourier, short-time Fourier, or wavelet transform methods exhibit many shortcomings when applied to EEGs complex nature ^{[21,23,24,25,26,27]}, such as nonlinear interactions and complexity. Examples of nonlinear methods implemented in the literature (for the investigation of ASD particularities) include higher-order statistical analysis ^[28], complexity analysis ^[29], and phase synchronisation analysis ^[30]. The bi-spectrum is an example of a higher-order statistical method that can describe the dependence among EEG frequency components. However, the bi-spectrum is usually normalised, which makes it useful for investigating the frequency and phase, but not the amplitude ^[28,31]. The same limitation is found in the concept of phase synchronisation, which considers the neural coupling independently of the amplitude. Therefore, this study focuses on complexity analysis techniques.

It has been noted in the literature that EEGs in ASD are characterised by a special sort of atypical complexity ^[32]. Consequently, EEG randomness and unpredictability are tempting features to consider. The principal techniques for the quantitative assessment of neural network dynamic complexity are the computation of entropy or fractality ^[33,34]. Fractality is a measure of self-similarity ^[35]. In ^[34], Katz's fractal dimension was applied to standard EEG bands to detect ASD. On the other hand, entropy is a physical quantity that refers to the irregularity in the studied signal. The most well-known related technique is multi-scale entropy (MSE), which was developed in ^[36]. MSE measures the entropy at different underlying timescales obtained by the coarse-graining procedure. The timescales are copies of the original time series at different resolutions and lengths, which are obtained through successive averaging of the original signal ^[37]. The concept behind timescales is that neural systems are controlled by multiple complex processes that interact at different temporal scales. In ^[29], MSE was implemented to detect infants at high risk of ASD. Differences in MSE between the high-risk group and control group across all scales were studied. The results provided an encouraging starting point for further investigations highlighting the differences in EEG complexity ^[33]. Hence, a number of studies have considered the application of MSE to ASD data ^{[29,38,39,40]}. For instance, in ^[40], MSE provided evidence of the relationship between ASD and cerebral dysfunction; significant MSE differences were observed between groups in the right hemisphere and central cortex. In addition, in ^[41], differences in entropy values were found in the frontal and temporal regions as well as in the whole-brain entropy. However, in a recent study ^[42], MSE was extensively examined to determine its ability as a measure of signal variance and spectral power. The results indicated that MSE failed to obtain accurate frequency-specific information on timescales, and it reflected dynamics at mismatched temporal scales. In addition, in ^[43], although the obtained results demonstrated differences in EEG complexity between ASD and neurotypical groups over temporo-parietal and occipital regions, the observed changes in the complexity values did not reflect the changes in the EEG power spectra. Furthermore, although MSE is based on the conditional probability sample entropy principle ^[44], which makes it sensitive to both spectral bandwidth and amplitude distribution ^[45], it has several limitations: first, it is dominated by bandwidth dependence; second, its sensitivity to the presence of high-frequency components and to the variance in larger shorter scales is problematic ^[45,46]; and finally, it is not completely suitable for non-stationary signals with amplitude- and fluctuation-based dispersion ^[46,47]. Therefore, MSE suffers from the overlap between frequency and amplitude information, which makes it imprecise for the study of both ^[48,49].

The present study suggests a novel, nonlinear, non-stationary, adaptive, data-driven, and accurate method for the classification of ASD and neurotypical groups based on EEG complexity and entropy without the shortcomings of MSE. Moreover, the new results are compared with the MSE outcomes. (1) The adaptability of the method originates from the fact that the scales are not extracted artificially as in MSE; rather they are extracted based on the inherent physiological components, and thus the frequency response of the method is dictated by the actual brain complexity ^[48]. (2) The adaptive technique overcomes the problems of variance in scales. (3) The proposed approach does not suffer from the short-time issue at larger scales, as all scales have the same length. (4) The proposed method has good sensitivity at all scales. (5) The overlap of the spectral and amplitude information is solved. (6) The dissociation of the frequency and amplitude distributions leads to precise features in both. (7) The classification is robust with accurate outcomes.

2.   Materials and methods

The proposed method is as follows: (a) each ASD and neurotypical EEG is decomposed using empirical mode decomposition (EMD) to obtain the intrinsic components (intrinsic mode functions). (b) The extracted components are normalised through the direct quadrature procedure, and (c) Hilbert transforms of the components are computed. (d) The analytic counterparts of the components (and normalised components) are obtained. (e) The instantaneous frequency function of each analytic normalised component is calculated. (f) The instantaneous amplitude function of each analytic component is calculated. (g) The Shannon entropy values of the instantaneous frequency and amplitude vectors are computed. (h) The entropy values are classified using a neural network (NN). (i) The obtained accuracy is compared to that obtained with the MSE-based classification. (j) The consistency of the entropy 3D mapping results with clinical data is assessed. All of the details of the procedure are described in the Materials and Methods section.

2.1.   Subjects

A total of 122 volunteers (61 ASD and 61 neurotypical subjects) with ages in the range of 4–13 y participated in this study. The participants had the following severity distribution: mild (18), moderate (21), and severe (22). The recording protocol and data collection were approved and registered by the Ethical Committee of Research Involving Human Subjects of the Institutional Review Board (IRB) at Jordan University of Science and Technology (JUST) and King Abdullah University Hospital (KAUH) according to the guidelines and principles embodied in the Declaration of Helsinki.

Each group included 48 males and 13 females, with 41 in the age range of 6–10 y, 12 in the range of 4–6 y, and 8 in the range of 11–13 y. ASD subjects were selected from special education institutions following assessment by neurologists in order to exclude subjects with comorbidities, sensory anomalies, or those taking medication. Neurotypical subjects were chosen from regular schools to match the demographic distribution of the ASD group. Healthy subjects with extreme mental capacities, hyperactivity, or those taking medication were excluded. The level of severity for each ASD subject was assessed by a neurologist and a rehabilitation expert based on medical and behavioural histories (verbal, social, non-verbal, etc.), as well as socio-behavioural scales such as the childhood autism rating scale (CARS), Autism Treatment Evaluation Checklist (ATEC), and Autism Diagnostic Observation Schedule (ADOS), which help evaluate the level of required substantial support related to the severity level.

2.2.   Data collection and pre-processing

The EEG measurement protocol was applied to each volunteer in a resting state for 20 min. The setup included an appropriately sized Waveguard head cap with 64 channels (ANT Neuro®). The electrodes were placed based on the 10–20 system. There was also a connection with an eego^TM sports high-quality amplifier (ANT Neuro®) with a sampling frequency of 500 Hz. The eego^TM sports ASA software was used for the data acquisition. Based on the recommendations in ^[62], the data were not subjected to other pre-processing steps because this could distort the data and influence the analysis. The artefact removal process is described in the discussion section.

2.3.   Proposed approach: EMD-based classification with dissociated amplitude and frequency entropies

2.3.1.   EEG decomposition

The components of a non-stationary nonlinear signal can be obtained via the sifting process of the EMD method ^[63]. These underlying components are called intrinsic mode functions (IMFs). Compared with classical spectral analysis techniques (e.g., Fourier transform), the temporal and spectral resolutions of EMD based on the Heisenberg Box are adaptive to the processes inherent in the signal. In addition, compared to other time–frequency distributions (e.g., wavelet transform), EMD is a data-driven nonlinear method ^[23] that does not require pre-assumptions or pre-models (e.g., a mother wavelet). A resulting component is considered an IMF if it satisfies two essential requirements: [number of extrema – number of zero crossings = 0 or 1], and [mean of the interpolation curves passing by the maxima and minima ≈ zero]. These two conditions imply an amplitude-modulated frequency-modulated mono-component nature.

All 64 EEG channels for each volunteer are decomposed using EMD. In the EMD technique, the extrema are initially interpolated to trace the envelopes. The technique used in the present study is cubic spline interpolation. Then, the local mean of the envelopes is computed. The result of subtracting the mean from the signal can be considered as IMF1 if it satisfies the conditions above; otherwise, the interpolation and local mean computation should be iterated. In the present study, the lowest value of the mean is controlled by an upper limit of iterations for every sifting (max = 2000). IMF2 is hence the first component of the subsequent decomposition of the signal after the removal of IMF1. The remainder (IMF3, etc.) can be obtained by analogy using the same steps. The final IMF represents the trend.

2.3.2.   Normalisation of intrinsic mode functions

Each obtained IMF is normalised using the direct quadrature (DQ) method ^[64]. Because the obtained IMFs are amplitude-modulated frequency-modulated (AM-FM) mono-components, the normalisation implemented in this study represents an important preliminary step toward the subsequent dissociation of the spectral and amplitude contents. The DQ method is performed as follows:

where n is the number of successive normalisations (in the present study, n = 5), the term IMF (k) represents the kth IMF sample, f_n (k) is the frequency-modulated component, and e_n (k) is the envelope passing through the maxima of the absolute values of f_n-1 (k), which is determined using the same approach as that used to calculate the envelope in the EMD.

2.3.3.   Analytic normalised intrinsic mode functions

Each IMF is considered as the real part of its counterpart analytic signal. The imaginary part is found using the Hilbert transform ^[24,65]:

where re(t) and im(t) are the real and imaginary parts of the analytic IMF, respectively; RE(f) and IM(f) are the Fourier transforms of re(t) and im(t), respectively; and ${F}^{-1}$ is the inverse Fourier transform.

2.3.4.   Calculation of the local frequency and amplitude

For each analytic IMF, the local (instantaneous) amplitude is calculated. In addition, the instantaneous pulsation and frequency ^{[21,23,24,25,26,66]} are computed based on the analytic normalised IMFs, as follows:

2.3.5.   Entropy calculation

Shannon entropy values of the resulting instantaneous frequency and amplitude vectors are calculated using MATLAB®. Therefore, two entropy values correspond to each IMF: the entropy of the frequency and the entropy of the amplitude. The entropy is computed as follows:

where pi is the probability of a value appearing in a stream of data.

The entropy values are calculated for each instantaneous frequency and amplitude vector. This is applied to all of the IMFs for all channels and subjects. Thus, an array with a size of [64 × (number of IMFs obtained from one channel) × 2] is obtained for each subject. For instance, if the number of IMFs extracted from each channel is 10, then the entropy array contains 1280 elements for each subject and 156,160 elements for all subjects.

Figure 1 illustrates the procedure used to decompose the channels, normalise the IMFs, obtain the exact instantaneous frequency values in the analytic normalised IMFs, obtain the exact instantaneous amplitude values in the analytic IMFs, and compute the entropy values.

2.4.   Multi-scale entropy (MSE) approach

2.4.1.   MSE calculation

MSE analysis is applied to all of the EEG channel time series. The MSE analysis is realised by computing the sample entropy (SE) at multiple timescales based on the coarse-graining procedure ^[29,38], as follows:

where the SE value is defined as the negative of the logarithmic probability of two similar sequences of m consecutive data points under the condition that the point remains similar at the next point (m + 1).

In the present study, m = 2, SDV is the standard deviation of the signal, r = 0.15 × SDV, N = 40,000 data points, the number of scale factors (SF) is 40, r is the coefficient of tolerance, d is the distance between points ${{\rm{x}}}_{{\rm{i}}}^{{\rm{m}}}$ and ${{\rm{x}}}_{{\rm{j}}}^{{\rm{m}}}$ in m-space, and ${{\rm{C}}}_{{\rm{i}}}^{{\rm{m}}}\left({\rm{r}}\right)$ is the number of all probable pairs (i, j) with [d < r × time series SDV].

In the MSE analysis, given an EEG time series of length N, [x₁, x₂, …, x_N], the signal is constructed into consecutive coarse-grained time series, ${{\rm{y}}}^{{\rm{ \mathsf{ τ} }}}$, where ${\rm{ \mathsf{ τ} }}$ is the scale factor. First, the original time series is divided into non-overlapping windows of length ${\rm{ \mathsf{ τ} }}$. Second, for each window, the average of the data points is calculated. Each element of a coarse-grained time series is computed as follows:

The length of the shortest coarse-grained time-series segment is 1000 data points.

The coarse-graining procedure provides consecutive time series, each with a length equal to that of the original time series (N) divided by the scale factor (${\rm{ \mathsf{ τ} }})$. Thus, for scale one, the time series ${{\rm{y}}}^{1}$ represents the original EEG signal. Finally, the SE is computed for each coarse-grained time series ^[29,36,38].

2.4.2.   MSE features

The present study explores four features of the MSE: the average MSE, high and low scale factor MSEs, and the area under the MSE curve (AUC). The average MSE values are computed by averaging the sample entropy values over the 40 SFs for each electrode. The AUC is commonly used in MSE studies ^[36,39] to compare the relative complexity of normalised EEG time series between groups at different timescales; the AUC is obtained by plotting the sample entropy of each scale factor as a function of the timescale. Therefore, the size of the feature array is 122 × 64 × 4.

2.5.   Dimension reduction

Because the obtained array of extracted features is relatively large, the process requires a dimension reduction process to decrease the processing cost, power, and time ^[67]. Dimension reduction is achieved using principal component analysis (PCA) with the aim of selecting the most significant factors that represent the highest percentage of the variance. The dimension reduction approach is applied to the entire set of features; the first five PCA components representing the highest variances are selected and used for the subsequent analyses.

2.6.   Classification of extracted features

Selected features are classified using an artificial neural network with a scaled conjugate gradient (SCG) backpropagation training algorithm ^[68], which is a type of fast, supervised learning. A confusion matrix is then obtained through cross-validation.

The architecture of the implemented feedforward NN consists of input, hidden, and output layers. The size of the hidden layer is adjusted to optimise the accuracy; the number of neurons is therefore set to 20. Optimisation is conducted on the entire reduced feature dataset as follows: first, the hyper-parameter values (the size of the hidden layer (nodes) as well as the number of layers) in previous successful models in the literature are referenced to define a preliminary set of ranges/combinations of values. Then, systematic experimentation is carried out to examine the performance of the "outlier" combinations, i.e., the values situated away from the combinations' scatter plot. The performance of those points is evaluated using performance analysis parameters such as the prediction accuracy; only the "outlier" values leading to low performance indexes are excluded. Finally, a grid search is performed on the remaining set of values. The values leading to the best cross-validation outcome are chosen.

The hyperbolic sigmoid function is used as the activation function. The performance of the training function is evaluated using performance analysis parameters such as the prediction accuracy. The SCG algorithm is used as the training function. This backpropagation function uses a mix of conjugate directions (GD) and the Levenberg–Marquardt (LM) method. Hence, the SCG algorithm with the Levenberg–Marquardt method does not execute a line search at each iteration, and the step size in the SCG algorithm is determined via a scaling mechanism. In addition, the SCG is fast and does not require user-dependent parameters. The design parameters are independently updated at each iteration.

The 10-fold cross-validation (CV) technique is used to validate and evaluate the outcomes related to the training and testing phases, as well as to prevent overfitting. In 10-fold CV, the original sample is randomly partitioned into 10 equal-sized subsamples; a single subsample is used for testing, and the others are used for training. The procedure is repeated 10 times and then averaged. This study utilises "repeated" k-fold cross validation, in which the procedure is repeated n times. The data are shuffled prior to each repetition, which results in a different split of the sample; finally, the mean over multiple k-fold CVs is considered.

2.7.   Statistical evaluation of outcomes

The classification outcomes are assessed based on the confusion matrix. Table 1 presents the relevant equations.

3.   Results

The main concept of the present study is to provide a robust indicator of EEG complexity that can be applied to ASD and neurotypical EEGs. The first step is data-driven decomposition using EMD. The resulting components are the IMFs that satisfy the required conditions of quasi-mono-components. Figure 2 (top) illustrates three selected IMFs obtained from the EMD of an ASD EEG channel (11 s). The IMFs are ordered from the highest frequency (top IMF) to the lowest frequency (bottom IMF). The IMFs are generated according to inherent physiological processes, not due to an imposed model (e.g., wavelet mother) or artificial smoothing (e.g., MSE). In addition, EMD does not yield components with variable lengths, as in MSE.

The second step is the tracking of the instantaneous frequency and amplitude calculated from the counterpart analytic normalised IMFs and analytic IMFs, respectively (bottom left and right of Figure 2). IMF normalisation by the DQ method ensures accurate tracking of the frequency without the effect of strong amplitude modulation. In contrast, the MSE procedure is sensitive to both frequency and amplitude without dissociation. This may be one of the reasons for the results reported in the literature regarding the unsuccessful compatibility of MSE with non-stationary signals ^[46,47]. EMD and DQ are therefore more "natural" ways of studying the pace of the mechanisms embedded in the brain and manifested in EEG oscillations. Furthermore, the instantaneous amplitude of the EEG clearly offers a very important impression of the profile of oscillations as well as their varying levels without deformation due to the frequency information. It is noteworthy that although the importance of EEG amplitude levels in ASD has been emphasised in the literature ^[69], it has not attracted sufficient attention in research on EEG frequencies, which may obscure an essential indicator.

Table 2 lists the values of the amplitude entropy computed for a certain IMF order across a number of selected EEG channels recorded for two children with neurotypical and ASD pre-diagnoses. Variations in the amplitude entropy across channels and subjects can be observed.

The type of entropy used in the present work is the Shannon entropy. The Shannon entropy measure depends on the values of a signal rather than the bandwidth, in contrast to the sample entropy ^[18,45,75]. It is therefore applied to both the extracted frequency and amplitude values to consider both temporal and spectral information ^[76,77].

The results in Figures 3, 4, and 5 are plotted after calculating the [average ± standard deviation] of the instantaneous frequency and amplitude entropy values for each IMF for all subjects. Hence, the figures illustrate the mega-profile rather than individual profiles.

Figure 3 shows the specific regions in the brains of the ASD volunteers exhibiting lower and higher amplitude entropy values compared to the neurotypical children. FP2, F3, M1, T7, M2, CP5, Pz, P4, F5, P6, and PO3 exhibit lower values. In contrast, Fz, FC1, T8, P7, F1, FCz, C2, PO5, FT7, and PO7 show higher values. The regions with higher entropy are positioned between those of lower values; this may be an indicator of a neural "compensation activity", where zones of higher (or lower) values of entropy suffer from the level of atypical synaptic connectivity, causing the neural dynamic network to be more (or less) stochastic ^[39,70,71]. This conclusion was previously reported by ^[41,72], which showed abnormal sample entropy in the insula and calcarine sulcus as well as in whole brain assessments in children with ASD, indicating increased randomness, less predictability, and less organisation. The indicated zones, together with the inference regarding the connectivity status, are consistent with the literature ^{[29,41,42,43]}, which has underlined the alterations in the frontal, parietal, and temporal regions and emphasised the changes in short-range and long-range connectivity in ASD.

Figure 4 illustrates the differences in frequency entropies of IMFs in the intermediate frequency range between neurotypical and ASD volunteers. ASD volunteers exhibit lower values in FP1, FP2, T7, AF7, F8, and T8, while they have higher values in C5, Cz, C1, FC1, and FCz. It can be observed that there is a good match between the directions and zones of differences in the amplitude and frequency entropy values. However, there are some exceptions, e.g., in temporal regions of one hemisphere. This finding is consistent with ^[40,73], which also noted the main entropy variation in one hemisphere. Figure 5 illustrates the differences in the frequency entropies of IMFs in the high-frequency range. It is clear that the FC zones are more affected by variability. Figure 6 illustrates the percentages of the data variance represented by the first five PCA factors. Figure 7 demonstrates the distribution of the first five PCA factors for 15 neurotypical and 15 ASD subjects based on the features of low-order IMFs; the results show two distinct clusters of points. The data in the figure are from 15 boys with severe ASD and 15 neurotypical boys with ages of 6–10 years.

In the implemented MSE analysis, the parameter values of the MSE algorithm are selected based on the recommendations in ^[18,21,38] to capture the majority of the main scales of the EEG activity. Table 3 lists the p-values for the different features extracted by MSE for selected channels. These values are calculated as follows: the value of the MSE feature is calculated at each channel for every subject. All of the values of that feature in the specific channel for neurotypical subjects are then combined in one column. All of the values of that feature in the specific channel for the ASD subjects are combined in another column. Then, the p-value between the two columns is calculated using the Student's t-test to determine whether the difference is significant. This procedure is repeated for the four MSE features in all channels. The MSE results demonstrate that the discrimination between features does not yield the same conclusions as above regarding the distribution of region-specific entropy, i.e., although the comparison between the counterpart features from a number of EEG channels (in ASD and neurotypical individuals) can provide a significant p-value of less than 0.05, it is not consistent with the outcomes of the suggested approach or those in the literature, which confounds the real mapping of EEG alterations. Indeed, in ^[46,74], it was shown that MSE occasionally led to incorrect low entropy owing to erroneous variance values of the scales. Table 4 lists the number of ASD subjects with entropy values higher than those of neurotypical subjects. It is clear that the low- and high-scale entropies are more discriminative. However, the numbers do not reflect the reported clinical data indicated in the literature, especially those available based on thorough imaging investigations.

Table 5 summarises the classification outcomes based on the MSE of the four extracted features, as well as the outcomes based on the proposed method. The accuracy, sensitivity, and specificity values in the table support the conclusion that the proposed method outperforms MSE, with values of 93.4%, 91.8%, and 95.1% compared to 66.4%, 34.43%, and 98.36%, respectively. Therefore, the present study suggests a novel, nonlinear, accurate classification of ASD and neurotypical groups based on EEG complexity and entropy.

4.   Discussion

It is worth noting that the goal of the inclusion of all severity levels in the selected sample in our study aims to verify the robustness of the proposed method when it is confronted with a pool comprising different severity levels, which has not been sufficiently targeted in the literature to the best of our knowledge. However, according to the analysis of the outcomes of the proposed procedure, the achieved sensitivity (91%) is lower than the specificity (95%) because of the relatively lower capability of the method to detect mild ASD compared to moderate and severe ASD owing to the similar features between mild ASD and neurotypical EEGs. Nevertheless, it is important to highlight that the EMD technique has been previously applied by our team to distinguish mild and severe cases with an accuracy of up to 97% ^[26].

In general, the accuracy values obtained in the literature for ASD EEG analyses using different ML protocols for classification purposes are in the range of 80–96% ^{[14,15,16,17,18]}. Nevertheless, it is noteworthy that the obtained accuracies are based on samples consisting of 20–40 subjects with a narrow range of severity levels, and comparisons of the proposed method with other approaches applied to ASD should thus be conducted carefully. Nonetheless, as the specific scope of this study is the exploration of ASD EEG entropy for classification purposes, the achieved accuracy will only be compared to accuracies from approaches with similar targets. In general, the two main tools used to characterise the entropy variations in ASD are spectral entropy and MSE. Spectral entropy is based on linear techniques that are incompatible with the nonlinear nature of EEGs; this has been concluded in previous studies by the first author ^{[21,23,25,27,66,78]}. In addition, there are shortcomings in MSE approaches; in ^[18], the method yielded an accuracy of 83%–98%. However, the number of subjects with ASD was only 9. In ^[29], the accuracy obtained with modified MSE was 70% for children in the age range of the present study. Several other approaches have focused only on general conclusions based on MSE p-value differences for the purposes of clinical investigation, studies of region-specific/age-specific/scale-specific variations, or therapy ^[79]. In ^[14], MSE was applied to a low number of EEG signals with the addition of the implicit function as squashing time algorithm (IFAST) technique; the obtained accuracy was 84%. In ^[80], the eye movement signal was added to entropy features to improve the accuracy. In ^[38], MSE was explored to distinguish only severity levels. In ^[81], variance values were added as indicators with MSE information to account for both temporal and spectral variations; the obtained accuracy was 63–79%. In ^[22], with the aim of overcoming the drawbacks of MSE, the authors combined the amplitude data with the sample entropy information in EEGs. In addition, in ^[50,51], MSE was applied with the addition of recurrence quantitative analysis (RQA) to create a more rigorous method. Another improvement was suggested in ^[52] by replacing coarse-graining with the moving average. Similarly, in ^[51], the authors added other nonlinear invariant measures derived from detrended fluctuation analysis (DFA) and recurrence quantitative analysis (RQA) to the MSE of ASD EEGs to overcome the limitations of MSE alone. Moreover, a number of modified MSE approaches have been proposed in the literature, such as refined multi-scale entropy (RMSE) ^[48], composite multi-scale entropy ^[53], refined composite multi-scale entropy (RCMSE) ^[54], modified multi-scale entropy for short time series ^[55], short-time multi-scale entropy (sMSE) ^[56], hierarchical entropy ^[57], multi-scale fuzzy sample entropy ^[58], multi-variate multi-scale entropy ^[59], and generalised multi-scale entropy ^[60]. However, these modified versions could only partially solve the issues of sensitivity to short time and noise or the frequency response yielding scales ^[61]. Furthermore, other types of entropy have been applied to ASD in the literature (e.g., Renyi, fuzzy) ^[32]. However, the obtained accuracy values were not higher than those achieved with MSE.

Several methods in previous research involved additional indicators, e.g., functional magnetic resonance imaging (fMRI) or structural magnetic resonance imaging (sMRI). However, while neuroimaging is an attractive and easily obtainable piece of clinical data, the experiments are usually limited by different sources of data, including harmonisation of scans for head motion, different MRI scanners and sequences, high cost, the need for training and expertise, a complex set-up, and fMRI data obtained under different conditions. Such variables limit the utility of this data in building personalised medicine models ^[84]. These disadvantages are generally not found in EEG applications. However, the disadvantage of EEG is its limited topologic usefulness or spatial resolution; EEG recordings do not provide accurate information about the localisation of the sources or whether they are surface or deep activities. Therefore, to overcome this limitation, the present study combines the aspect of computerised 3D mapping to distinguish between normal and abnormal cortex functions. Hence, the proposed solution satisfies the requirements of spatial and temporal resolutions to a certain extent. Furthermore, portable EEG can be more suitable for routine use than MRI or magnetoencephalography (MEG). It is also worth noting that the results achieved in the present study are not contradictory to those obtained by neuroimaging, e.g., in ^[70,72]. In addition, the accuracy is comparable to or higher than that of those approaches ^[85,86].

In the procedure in this study, four principal methodological issues were encountered. The first was the management of the plentiful artefacts detected in children, especially those with ASD, e.g., myogenic or group-specific activities. The recorded EEG segments required careful inspection to correct these artefacts. Second, there was the possibility of hypothesis confirmation and data significance by chance owing to the consideration of many features. Third, the size of the hidden layer needed to be adjusted. Fourth, the type of CV used required consideration. To solve the first problem, (1) artefacts or odd activities were marked in the (eego sport^TM, ANT Neuro®) software, (2) two neurologists involved in the study visually inspected the recorded epochs and confirmed artefact-free segments, and (3) all persisting artefacts were removed via (ASA 4.10.1, ANT Neuro®) software by source separation and statistical regression techniques. To solve the second problem, the PCA analysis technique was used to convert the feature array into significant factors. Then, the most significant factors were selected based on the highest representations of variance, which led to a reduced dimension. However, an important aspect of the application of PCA should not be ignored. The PCA method was used herein as an unsupervised learning technique to reduce the dimension and select the most important basis vectors through eigenvector information. PCA was applied to the whole dataset, which may raise questions regarding the bias and generalizability of the results. The goal of this decision was to establish a first step towards a universal standard number of factors for later use in autism diagnosis approaches based on the present study, which is similar to the concept adopted in several previous pioneering articles about autism ^[30]. Nevertheless, the authors conducted an additional investigation focused on over-fitting and generalizability by applying PCA to the training set in every iteration of the CV and then transforming the fitted model to the test data. The results of that investigation yielded an accuracy of classification of 91.02%, which falls in the confidence interval of the accuracy obtained when PCA was applied to the whole feature dataset, indicating good generalisation performance. In addition, the same conclusion was obtained when PCA was applied to a whole set of training and testing data, and the fitted model was transformed into a validation set. To solve the third problem, the optimisation of the size of the hidden layer (nodes) as well as the number of layers was applied to the whole dataset rather than the training set in the CV iterations because: a) this is similar to the common approach in literature related to autism learning ^[12], and b) the performance indicators were based on the outcomes of a robust "repeated" k-fold CV. Finally, to solve the fourth problem, the study adopted "repeated" k-fold CV. The other types of CV are hold-out and other sub-types of k-fold, such as leave-one-out, leave-p-out, stratified, complete, and nested. The reason the "repeated" type was used is that the criteria of computation cost, generalizability error, performance improvement, suitability for the level of variance between samples, and bias–variance trade-off should all be taken into account. For example, complete and leave-one-out CVs are computationally expensive. The standard k-fold CV does not account for the variance effectively. Nested k-fold CV is performed via k-fold CV in every fold of the CV; it is thus computationally expensive and it is more suitable for cases requiring highly accurate hyper-parameter tuning during model evaluation. "Repeated" k-fold CV is therefore an intermediate solution that is capable of improving the estimate of generalisation error by reducing random uncertainties due to model instability. It also ensured that samples of the same subject were not included in the training and test sets of a CV iteration.

Several aspects led to the increase in accuracy of the proposed approach compared to other entropy-based techniques, particularly MSE. In ^[82,83], it was shown that EMD was more adapted to non-stationary nonlinear processes because its procedure is more data-driven than the linear operations of smoothing used in MSE. It has also been shown that EMD has a better potential to detect high frequencies. In addition to the adaptability and improved performance of EMD, the ^[22] demonstrated that using sample entropy without considering the EEG amplitude intensity level was not useful for the study of EEG in the case of cerebral injury. This is a valid reason for using an accurate amplitude and exact frequency as complementary indicators in the present study.

5.   Conclusions

The proposed study on the classification of ASD and neurotypical subjects outperforms other entropy-based approaches as well as provides high accuracy compared with other related ML-based techniques. Future work will combine entropy and connectivity indicators to cover the two main alterations in ASD and improve the classification outcomes.

Acknowledgments

This study was conducted under the umbrella of the Scientific Research Support Fund (SRSF), project #MPH/1/11/2014, Ministry of Higher Education and Scientific Research, Hashemite Kingdom of Jordan.

Conflict of interest

All authors declare no conflicts of interest in this paper.