A Fusion-Based Machine Learning Approach for Autism Detection in Young Children Using Magnetoencephalography Signals

In this study, we focused on the age range of 4–7 years (47–86 months) as our primary goal is to identify early pathophysiological biomarkers of ASD in young children. There were two groups of children: (i) ASD and (ii) TD. The ASD group included 30 children (4 females) with autism spectrum disorder with a mean (± SD) age of 64.66 (± 10.12) months, and the TD group included 30 typically developing children (4 females) with a mean (± SD) age of 64.83 (± 10.51) months; the two groups did not significantly differ in age (two-tailed t-test, p > 0.95). The ASD children were diagnosed by an experienced psychologist and clinical psychiatrist using different screening tools; such as the Autism Diagnostic Observational Schedule, Generic (ADOS) (Rutter et al., 2003), the Diagnostic Interview for Social and Communication Disorders (DISCO) (Wing et al., 2002) and the DSM-IV (APA, 1994) criteria before taking the MEG data. Specifically, the ADOS cutoff and parent report on the social communication questionnaire are the popular screening tools with which the ASD was confirmed. The children in the TD group had no behaviour or language problems as reported by their parents. The written informed consents were obtained from the parents prior to the data acquisition. The Ethics Committee of Kanazawa University Hospital, Japan, approved this study protocol, and the experiment was conducted as per the World Declaration of Helsinki.

The MEG signals were acquired using a 151-channel superconducting quantum interference device (SQUID) whole head coaxial gradiometer MEG system (PQ1151R; Yokogawa/KIT, Kanazawa, Japan) that was customized for children in a magnetically shielded room (Daido Steel, Nagoya, Japan). The custom child-sized MEG system ensured the sensors were placed quickly and effectively to minimize head movement (Johnson et al., 2010). An experimenter was also present in the MEG room at the time of recording to ensure that the children remained relaxed to prevent movement throughout the recording. Resting MEG data were recorded from all children while they lay supine on the bed and viewed a silent cartoon of their own choice that was selected before the recording. There was no task involved, and the duration of the resting MEG data was 180 s. The sampling frequency was 1000 Hz and filtered online with a 200 Hz low-pass filter.

In this experiment, the preprocessing and analysis of MEG data were incorporated by MATLAB based toolbox, FieldTrip (Oostenveld et al., 2010), and by custom-made MATLAB scripts. The data were first visually inspected for large artefacts. Afterward, bad sensors that were either flat or noisy were detected and replaced by nearest-neighbour interpolation. Next, the MEG dataset was subsequently filtered between 1 and 100 Hz, and to reduce powerline noise, a notch filter was applied at 50 Hz. Finally, independent component analysis (ICA) was applied to remove eye-blink related artifacts.

Given the role of large scale brain responses characterizing ASD, we divided the broadband MEG signal into six standard frequency bands (fb) (Donner & Siegel, 2011): delta-band (1–4 Hz), theta-band (4–8 Hz), alpha-band (8–13 Hz), beta-band (13–30 Hz), lower gamma-band (30–50 Hz) and higher gamma-band (50–100 Hz). Our primary focus was classifying ASD children from TD children using resting-state MEG signals. Towards this, we adopted a machine learning framework (Barik et al., 2019a, 2019b; Barik et al., 2019a, 2019b); machine learning based approaches are useful in diagnostic and intervention research in clinical neuroimaging (Iniesta et al., 2016; Vieora et al., 2017), including ASD (Hyde et al., 2019). Furthermore, experts have a certain degree of subjectivity (supported by rigorous statistical analysis) in creating the diagnostic instruments (i.e., ADOS, DISCO). Hence, it is reasonable to believe that objective machine learning methods may provide more reliable performance and increased efficiency by reducing redundancy within an instrument. The supervised machine learning framework consists of the following main blocks—feature extraction, feature selection and modelling which leads to the classification of two different classes, TD and ASD. We followed a nonliear artificial neural network (ANN) based modelling scheme (Bishop, 1995; Haykin & Network, 2004) to learn the underlying relationship among the selected fetaures of a particular class to separate it from the other class.

Frequency domain analysis is the most common and familiar analysis for characterizing neural oscillations of a signal. Fourier Transform (FT) is the keystone of this frequency domain analysis. FT decomposes a time-domain signal into its constituent frequencies with different amplitude and phase angles (Oppenheim, 1999). It was found that phase angles are independent of magnitude (Cohen, 2014). We investigated the power spectral density (PSD) and phase-based features to compare classification performance. In this paper, we explained a frequency band-wise preferred phase angle (PPA) as a novel feature in MEG signal processing. The power spectral density of a continuous-time and finite power signal estimates the power of each of the constituent frequency components. For each channel of the MEG, the PSD is estimated using Welch’s method (Welch, 1967), where the whole time-series (180 s) is divided into equal 8 segments with 50% overlap. Hence, the length of each of the 8 segments is 40 s. Then, each segment is windowed with a Hamming window, and its magnitude squared FT is computed. The PSD estimate for each channel is obtained by averaging the periodogram estimates across segments, as shown below:

The frequency domain representation, \(X_{n} \left( f \right)\) is the FTof time signal \(x_{n} \left( t \right)\). After getting all the PSD over 1 to 100 Hz in steps of 1 Hz frequency components, we averaged the power spectrum over the six frequency bands as shown in the following equation.where fb represents each frequency band, and N_fb represents the number of components in each fb. For example, in case of alpha band, fb_min = 8, and fb_max = 13, so N_fb = 6, as our frequency resolution was 1 Hz. The frequency band-wise power spectral density estimates are denoted as PSD_fb.

A frequency band-wise phase angle-based feature named phase angle clustering was initially developed in a previous study (Barik et al., 2020). In this study, we have analyzed this phase-based feature in detail. After the computation of Fast Fourier Transform (FFT), we get phase angles (θ) at each frequency component (f) from 1 to 100 Hz in steps of 1 Hz for each channel (Oppenheim, 1999). The phase angle is calculated using simple Euler's formula \(e^{{i\theta_{\left( f \right)} }}\) (Strang, 1991). Phase angles are commonly represented in polar space, ranging from 0 to 2π. These can be illustrated as vectors on a unit circle. In the polar space, the length of the average phase angle vector is known as phase consistency (Cohen, 2014), and the angle of this resultant phase vector is termed as preferred phase angle (PPA). For each MEG channel, in a given frequency band (\(fb\)), PPA_fb is calculated as:where all the notations are similar to the first equation. Here, Fig. 1 is shown for a detailed illustration of PPA_fb. Figure 1 represents the PPA of theta band (PPA_theta), where each unit length phase angles depict 4 Hz, 5 Hz,…, 8 Hz frequency components. These phase angles of different frequency bins within a frequency band are represented by different colors (line of blue, orange, yellow, violet, green, cyan). The average phase angle represents PPA clustered of theta band frequencies (indicated by red color), and the average phase consistency (indicated by magenta color) is the length of the average phase angle vector in the polar space. Similarly, it is found for the alpha band in the right side from 8 to 13 Hz with the gap of 1 Hz frequency components. All the ASD subjects’ PPA_theta and PPA_alpha are shown in Fig. 1 from a randomly chosen channel (sensor: AG102). Here, in the polar space, PPA_theta value is 77.95°, and PPA_alpha value is 196.85°.

For comparison of classification performance, we computed commonly used power-based along with phase-based spectral domain features. In addition, PSD_fb and PPA_fb features are computed across all MEG channels. Hence, for both features (PSD_fb and PPA_fb), the dimension for each participant is 906 [channel (151) × frequency-band (6)] for the overall brain, whereas, for each hemispheric features, the dimension is 420 [channel (70) × frequency-band (6)]. The features are normalized using z-score normalization (Kochendörffer, 1965) before proceeding to the feature selection step.

To note that as the feature values are normalized, the time-shift-related phase offset does not show up in the normalized PPA features. Here, the normalization confines the values of PPA features to ± 1 rad. Therefore, the normalization reduces the arbitrariness of the range of the computed PPA values which is a standard step to feed extracted fetaures for modelling. Accordingly, it also confines the time-shift of the data (phase shift of the features) to the same range.

The objective of feature selection is to extract a subset of the most relevant features while removing redundant ones (Yu & Liu, 2003). Feature selection is helpful as it reduces feature dimensions. When there are many features, the simple filter feature selection method, such as t-test, performs better than the complex wrapper and embedded techniques (Haury et al., 2011). We used a t-test and ensured the normality assumption (see Fig. 2). Even it is to be noticed that the relevance ranking strategies (e.g., t-test) take moderately less calculation time (Chandrashekar & Sahin, 2014) for feature selection. We have used p-values to rank the features. The p-value is the probability of obtaining the calculated t-statistics. From the ranked features, we selected a subset of the features that were below the chosen p-value. These p-values are used as thresholds to obtain a coarse selection of features to reduce the feature dimension. Different thresholds were employed to investigate the effect of increasing the number of selected features (Mwangi et al., 2014; Wang et al., 2014).

After feature selection step, the classification model finds the underlying associations or relationships among the normalized features to arrive at the decision model. Using t-distributed Stochastic Neighbor Embedding (t-SNE) plots (Hinton & Roweis, 2002), which convert high-dimensional data into its two-dimensional analogues, we present a visualization of the clustering effect for the fetaures extracted from TD and ASD classes in Fig. 3. This t-SNE embeds high-dimensional points in low dimensions solely based on their relative similarities that correspond closely to the true labels. Figure 3, also shows reduction in the number of outliers after z-score normalization which helps in modeling by making the decision boundary simpler, thereby, enhancing the classification model's discriminative ability.

This work aims to study and analyze how the underlying neural mechanisms of autistic children differ from typically developing children based on a machine learning framework. Here, MEG signals are used to classify ASD children from TD children. The detailed classification process through a block diagram is illustrated in Fig. 4. It must be stressed here that the training data set and test data are completely independent in terms of different subjects, i.e., the subjects of the test dataset are different from the subjects of the training dataset. The problem of distinguishing between the TD and the ASD children is posed as a two-class classification problem where the classification accuracy indicates the distinguishable ability of the spectral features. In this framework, a two-layered feedforward back-propagation artificial neural network (ANN) (Bishop, 1995; Haykin & Network, 2004) is used as a classifier. This ANN model finds the underlying relationships among these features for a particular class or category by learning the training data. The ANN is trained through an iterative, backpropagation algorithm which is essentially nonlinear and finds higher order statistical distribution related patterns or relationship embedded in the normalized feature vectors. This underlying relationship comes from the relative positioning of the feature vectors that remains invariant for a particular class and caters to the decision making process. We incorporated fivefold cross-validation (CV) in our classification model to remove the confusion of data biases. In this ANN model, a hidden layer consisted of 10 neurons, and an output layer included two neurons indicating the two classes. The scaled conjugate gradient descent algorithm (Møller, 1993) is used to train the neural network, as it converges fast and does not get stuck on local minima.

For the analysis using ANN, the hyperbolic tangent sigmoid transfer function is set as activation function, and the analysis has been carried out considering 10,000 as an upper limit for cycles and 10e-5 as mean square error. Normalization of feature vector is embedded before classification, and this z-score normalization rescales the feature values such that the mean of all of the values is 0 and the standard deviation is 1 (Kochendörffer, 1965). Early stopping criteria are employed in the validation set to restrict the overfitting of the current model. In this classification process, the three subsets are training, validation, and testing. The available data of all participants are divided into three subsets, with test set subjects differing in each fold of the fivefold CV. The following is the distribution of all 60 patients in each fold: the training set has 36 participants (18 from each class), the validation set contains 12 subjects (6 from each class), and the test set contains the remaining 12 individuals.The first subset is a training set, wherein gradient computation, network weights, and bias updation are performed. The second subset (validation set) validates the errors throughout the training process. Our model tried to minimize the errors through the validation process from the beginning. The continuation of this process is stopped when the validation error showed reversal trends. Here, the stopping criterion is the maximum validation error increasing check up to 6 times. Finally, the model with minimal validation error is used for the testing process. The third subset, i.e., the test set, evaluates the average classification performances over all CV. The fivefold CV technique is used with mutually exclusive test sets to calculate the performance matrices in our classification process. Our model calculated the classification accuracy, sensitivity, specificity, and standard error of the mean (SEM) for all the feature types. Here, sensitivity and specificity referred to the percentages of ASD and TD children correctly identified, respectively. Finally, the percentage to correctly classify the ASD and TD children in their respective classes is indicated as the classification accuracy. In this experiment, we calculated average classification accuracy, which averaged over fivefold classification outcomes.

It is challenging for any pattern classification task to model the perfect one by only one hypothesis. An effective fusion method is necessary to combine information from multiple single modality systems in a multimodal experimental framework. Information from multiple sources of different attributes can be consolidated at several levels, including the feature extraction level, score level, and decision level. This paper examined the feature level's performance and the sum rule-based score level fusion.

The analysis of \(PPA_{all}\), \(PPA_{Left}\) and \(PPA_{Right}\) features were carried out for each participant. To compare the performances with PSD features, we have investigated \(PSD_{all}\), \(PSD_{Left}\) and \(PSD_{Right}\) features. \(PPA_{all}\) and \(PSD_{all}\) features were chosen from all 151 MEG sensors. The \(PPA_{Left}\) and \(PPA_{Right }\) correspond to the left and right hemispheric sensors, respectively. Hemispheric features were similarly considered for \(PSD_{Left}\) and \(PSD_{Right}\) features. In these hemispheric features, 11 midline sensors were excluded. Hence, in each hemisphere, there are 70 sensors only.

The classification outcome of PPA features illustrated in Fig. 6. The classification accuracy was graphically illustrated with the empirical chance level of 60% (Combrisson & Jerbi, 2015) by setting the p-value threshold in between 0.005 to 0.05 with the interval of p-value 0.005. Thus, a suitable threshold was empirically ascertained for selecting the features. With the stricter p-value threshold, selecting discriminating features becomes less, resulting in degradation in the classification accuracy. As the classification accuracy did not show substantial improvement than the chance level, the p-value was gradually increased to find out the optimal threshold beyond which the classification accuracy almost saturated. The optimal p-value was chosen using the classification results of the validation set data only.

Table 1 shows the classification performances of the ANN classifier using these PPA feature types, whereas similarly, Table 2 represents the same for PSD feature types. We only mentioned the decision outcomes at optimal p-values for each feature type. Comparing the results of both the features, it is found that PPA features yield better classification accuracy than PSD features for ASD classification in children from their resting state MEG data. The classification performance of PPA_all is 88.20 ± 3.87%, whereas the same for PSD_all is 82.13 ± 2.11%. For PSD features, right hemispheric features contained more discriminative information than the left hemispheric ones. However, for PPA features based on classification outcome, we have found that the distinguishing information of either left or right hemisphere was comparable. The hemispheric differences were statistically tested for both the PPA and PSD features. The statistical t-test indicated that the hemispheric PPA feature groups (i.e.,\(PPA_{Left}\) and \(PPA_{Right}\)) did not significantly differ in performances (p > 0.95). However, the hemispheric PSD feature groups significantly differed (p < 0.002) in fold-wise test accuracy, according to a paired t-test.

By comparing the results of Tables 1 and 2, this phase based feature (PPA_all) demonstrated superior (88%) classification performance over the one (82%) based on commonly used power spectral density (PSD_all) based feature. As a consequence, this phase-based feature is investigated further. The mean angle of the PPA features of each frequency band is illustrated here in polar plots of Fig. 7 for both ASD and TD classes. It represents frequency band-wise PPA features at the AG102 sensor (located at right parietal area) for 30 participants before normalization. To give an example, we randomly chose one sensor (AG102) for each of the frequency bands. Observations from this Fig. 7 reflect that the PPA value ranges from 0° to 360°. The mean angle of PPA (angle of the magenta line) reflects that the mean angle of PPA of the ASD class is greater than the mean angle of the PPA of the TD class in all frequency bands except the low gamma band.

Using a data-driven approach, we can classify autistic children from normal children. However, this is not the only target of this study. We also explored neural oscillations and spatial pattern analysis to identify the frequency bands and spatial patterns that consistently show discriminative ability between two classes across all the subjects. Most discriminating PPA features are found here from the theta frequency band. Spatial patterns of PPA discriminatory features are located in the central, parietal, and frontal brain regions. However, in this frequency band-wise spatial pattern analysis, discriminating PSD features are found mostly from high gamma-band oscillations, and these discriminating features are mostly selected from central, parietal, and temporal regions. The row-wise representation of Fig. 8 represents the discriminatory features of PPA_all and PSD_all. The frequency band-wise distribution of selected features is illustrated in Fig. 8a,c and the spatial pattern of those features ispresented in topoplots of Fig. 8b,d.

For feature-level fusion, we combined PSD_all and PPA_all features for each participant. Since the feature dimension was large, we performed feature selection as mentioned earlier. Following ANN based classification, we only choose the optimal p-value, i.e., p < 0.05 (the highest threshold) for both the PSD and PPA features. After feature-level fusion, the average classification accuracy of features is 94.46 ± 1.35% for the specific p-value mentioned above. The sensitivity, specificity, and the number of selected features are 94.13 ± 2.90%, 94.80 ± 2.52%, and 508, respectively. Among these 508 discriminative features, 20 discriminative features are common in both PSD and PPA. These discriminative features are mostly from theta-band oscillations, which correlate with autistic symptomatology. The discriminating features, those are common in individual feature types, are shown in Fig. 9. Notably, the feature-level fusion of PSD and PPA features performance is, on average, 6% better than the classifier based on PPA features.

For score level-fusion, the classification performance of sum-rule based score-level fusion is 98.33 ± 0.74% at p-values threshold on < 0.02. Here probability weight is set to 0.2 (w_k) of PSD features, and that for PPA feature is 0.8. Table 3 shows the classification performances of ANN classifier using these fusion-based models, i.e., feature-level fusion and score-level fusion types. We only mentioned the decision outcomes at optimal p-values for each fusion-based model type. Using ANN, the performance of score-level fusion has improved 4% compared to feature-level fusion.

Figure 10 shows the average classification performance of power spectral density (PSD_all) features, preferred phase angle (PPA_all) features, and fused feature vector after feature-level fusion and score-level fusion concerning different p-value thresholds used in the feature selection method. After feature and score-level fusion, using spectral-domain analysis, we achieve a classification accuracy of 94% and 98%, respectively, thereby score-level fusion offering a slightly better classification accuracy than the feature-level fusion. Finally, we show the comparison of classification accuracy, sensitivity, and specificity at individual feature level (PSD and PPA), feature-level fusion, and score-level fusion in Fig. 11 through barplot.

In autism detection using MEG signals, this fusion based model is a novel machine learning approach. To understand how efficient this fusion of PSD and PPA features is, we must check which subjects are misclassified in individual models. It is illustrated using boxplots in Figs. 12 and 13 that the PSD model misclassified 11 subjects and PPA misclassified 7 subjects, respectively. The first 30 boxplots (in ‘blue’color) represent the TD class, and the rest 30 represent the ASD class. The TD subjects are correctly classified if the output scores are ≤ 0.5, and ASD subjects are correctly classified if those scores are ≥ 0.5. Among these sets of subjects, only the 17th subject is commonly misclassified. Except for this (17th) subject, different sets of subjects are misclassified from the ANN model, using each model (with PSD and PPA features). Hence, from the perspective of classification scores or output, it is evident that the PSD feature cannot misclassify a different set of subjects than the PPA characteristics. Instead of criteria of fusion of individual features, here on the other way again, we showed that two individual models are complementary; hence multimodal fusion is more effective to classify ASD and TD in this machine learning framework.

We employed a pattern classification approach using artificial neural network (ANN) based modelling for discriminating ASD and TD children using PPA features. As noted that phase angles are independent of power (Cohen, 2014), we compared our results with power spectrum density (PSD) features that capture the mean amplitude of the neural oscillations in each frequency band. We found that PPA yields better classification accuracy than PSD, suggesting that phase angles aggregated in each frequency reveal significant phase alignment effects that distinguish autistic brain function from typically developing brain. Furthermore, in spectral analysis, the phase-based measures reflect considerable improvement over power-based measures, giving insight into the presence /absence of systematic rhythms across oscillations in each neural frequency band. Hence, this result suggests a new way to quantify temporal synchronization of brain oscillations as an effective measure in discriminating ASD from TD children.

According to the E-I hypothesis, when the neuronal membrane potential is balanced, the relative phase and magnitude of excitation-inhibition (E-I) are also balanced, reflecting systematic synchronous neuronal oscillations in the brain (Sohal & Rubenstein, 2019; Zhou & Yu, 2018). The imbalance of E-I can disturb this signature of relative neural firing time and magnitude, reflecting synchronous neural oscillation’ power and phase. Hence, spectral power and phase are disrupted in functioning. Our proposed feature PPA measures the resultant phase angle of frequency bins over a frequency band. It reflects the phase of excitability of the underlying neuronal assemblies oscillating at different nearby frequencies but belonging to a frequency band. It is not directly linked to the phase of the membrane potential. Instead, PPA may be linked to the cortex’s imbalanced synaptic excitation and inhibition during the resting state and sensory processing. Our experimental outcome shows that the phase of excitability of underlying neurons of the autistic brain is different from the typically developing brain in each frequency band.

The length of the average phase angle vector in polar space within a typical frequency band is known as phase consistency (Cohen, 2014), and the angle of this resultant vector is known as the preferred phase angle. Phase angles are commonly represented in polar space, ranging from 0 to 2π. These can be illustrated as vectors on a unit circle. These PPA ranges are different in each class, particularly for discriminant features. Frequency band-wise observations in phase based feature reflect that the PPA of ASD class is greater than PPA of TD class in all frequency bands except low gamma (Fig. 7). PPA represents how much phase synchronization within a frequency band deviates from the initial phase in a counter-clockwise direction. Using PPA, it is clear that the initial phase angle deviation of phase in autism class is more than in TD class.

Using a data-driven exploratory approach, we showed that it is possible to classify ASD from TD children. However, we also explored the ASD pathophysiology by performing frequency band-wise analysis and spatial pattern analysis of the features over whole cortex MEG sensors. More speculatively, the results of neural oscillations and spatial pattern analysis are presented in Fig. 8. In the case of PPA features, neural oscillations that distinguished the ASD class from the TD class were not specific to any particular frequency band. Instead, the discriminatory PPA features were from all the frequency bands; a marginal dominance (p-value = 0.054) of theta band (4–8 Hz) was observed. In the spatial pattern analysis, the discriminatory features were dominant in the central, parietal, and frontal brain regions. In the case of PSD, discriminatory features were observed in the high-frequency gamma band. Discriminating power-based spectral features are found in central, parietal, and temporal brain regions. Already, NeuroSPECT work concluded that temporal and frontal lobe dysfunction is significantly involved in autism spectrum disorder (Goldberg et al., 1999). For both spectral features, discriminating features are common from central and parietal brain regions. The explanation of the relationship between autistic symptomatology and these brain areas is found in the literature (Courchesne et al., 1993; DeRamus et al., 2014). It was found that compared to TD, ASD children showed significantly enhanced activation in the parietal region (particularly the angular gyrus) while detecting the location of objects on visuospatial processing (DeRamus et al., 2014). Even parietal lobe abnormalities of autistic children were detected in neuroanatomy investigation (Courchesne et al., 1993), where the cerebrum was reported as the origin of this abnormality. The cerebrum region mediates many high cognitive functions; several of those (i.e., social communication, language, reasoning, planning, and organizing) are severely disrupted in the autistic brain.

Previous research has found abnormal high gamma spectral power in young ASD children (Lushchekina et al., 2012; Orekhova et al., 2007; van Diessen et al., 2015). Our study also found that the discriminating PSD features were generally from high gamma oscillations. Likewise, this finding aligns with the studies that identified unusual activation of gamma-band oscillations in ASD children (Lushchekina et al., 2012; Orekhova et al., 2007; van Diessen et al., 2015). In addition, this study demonstrates a rightward-lateralized neural oscillation in power analysis similar to previous literature (Cornew et al., 2012; Kikuchi et al., 2013). In the literature of the last decade, phase-based features were investigated in the connectivity analysis (O'Reilly et al., 2017; Velazquez et al., 2009; Ye et al., 2014), however not yet explored in spectral-domain analysis for the autism classification problem. In spectral analysis, the preferred phase angle is a novel feature, where phase features are implemented even in frequency band-wise clustering. Our proposed feature yields better classification accuracy (88.20 ± 3.87%) than the PSD features (82.13 ± 2.11%). So we suggest that PPA can also be used as a useful measure to identify early ASD in children from their ongoing brain responses.

To explore the frequency band-wise analysis and spatial pattern analysis of the fusion-based system, features over all 151 sensors were also investigated. The average classification performance in the feature-level fusion-based system is 94.46 ± 1.35% for the p-value < 0.05. In this specific p-value threshold (p < 0.05), the number of selected features is 508. Among these 508 discriminative features, 20 discriminative features are common in both PSD and PPA (Fig. 9). These discriminative features are mostly from theta-band oscillations, corresponding to autistic symptomatology. Finally, this fusion-based system estimates are correlated with ASD scores indicating clinical relevance of the feature level fusion of PPA and PSD matrices.

In this study, autism detection in young children is extended using spectral features based on large-scale neural oscillations in a machine learning classifier in a fused system. Overall this work demonstrates the resting-state neural function aberration, a key component of the autism symptom profile (Kessler et al., 2016). Even to distinguish ASD children from TD children, disrupted oscillatory synchronization is found in multiple frequency bands as reported in (Simon & Wallace, 2016). In the current study, theta oscillations are found as discriminative features closely linked to changes in memory state, memory performance, and auditory responses. Hence, this theta-band preferred phase angle can also be an optimal measure to discriminate ASD from TD, as autistic children have specific memory and memory strengths difficulties. Autistic children use fewer strategies spontaneously to retrieve information (e.g., visual and rote memory) (Williams et al., 2006). Multiple studies suggest that memory deficit observed in ASD has a biological explanation linked to abnormalities in the hippocampus and other neural regions responsible for strategic regulation, such as the amygdala and the frontal cortex (Wojcik et al., 2013).

Our approach is not without limitations. First, we restricted our analysis to the MEG sensor space only, so we cannot provide any clear conclusion about the precise involvement of underlying brain sources. Future studies could adopt advanced source localization methods that would provide a more accurate measure of brain activity to examine PPA features in ASD better. Second, both PPA and PSD features were computed at each sensor, so our findings cannot provide any direct information about the nature of information transfer between near and distant sensor regions. Third, we only analyzed a brief 3 min period of brain data, and it would be essential to obtain a longer recording to validate the robustness of our findings. Fourth, our findings suggested some complementarities between PPA and PSD; however, additional studies are needed to substantiate this claim further. For example, future studies with larger sample sizes could match children on PSD values and then determine whether PPA would still discriminate between ASD and TD (a similar approach can be adopted by first matching on PPA values). Finally, we termed our recording as resting-state brain activity while the children watched cartoons of their choice, and such natural viewing conditions might have introduced various heterogeneities in the recording. Though this type of data collection is increasingly being used in studies with young children (Richardson et al., 2018), one cannot distinguish between intrinsic and task-driven contributions to the reported patterns, both in terms of phase and amplitude, of brain oscillations. Future studies could investigate individual contributions of both task-driven and intrinsic by collecting MEG data during the task and actual resting state from the same child.