Noninvasive Biomarkers for Assessing the Excitatory/Inhibitory Imbalance in Children with Epilepsy

Participants

We recruited 49 children and adolescents with a diagnosis of epilepsy (EP group) [age, 7−19 years; mean ± standard deviations (SD), 14.61 ± 2.9; 26 females] and 48 neurotypical healthy controls (control group; age, 10–18 years; mean ± SD, 13.58 ± 2.48; 20 females). Children from the EP group were recruited from the epilepsy clinic at Cook Children's Medical Center (CCMC). Neurotypical healthy controls were recruited from the local community. We assessed age differences between the two groups and found no differences (p > 0.05). In the EP group, 33 patients (67%) were diagnosed with focal epilepsy, while 16 patients (33%) had generalized epilepsy. Among the patients, 27 (55%) had either a normal magnetic resonance imaging (MRI) or findings unrelated to epilepsy, while 10 patients (20%) had an MRI-detected anatomical lesion related to epilepsy. MRI findings were unavailable for the remaining 12 patients (24%). The most common pathology among patients was malformation of cortical development, found in five patients (10%), followed by stroke in two patients (4%), tumor in one patient (2%), and hippocampal sclerosis in one patient (2%). One patient (2%) had both malformation of cortical development and gliosis. Pathology was unknown for 39 patients (80%), including those with a normal MRI and those without MRI data. For patients, epilepsy onset ranged from 6 months - 16 years (mean ± SD, 8.6 ± 4.49 years), whereas epilepsy duration ranged from 6 months to 18 years (mean ± SD, 6.38 ± 4.2 years). Epilepsy onset age was unknown for four patients. Patients received individualized ASM regimens with the number of concurrently prescribed ASMs ranging from one to five. This number reflects the total count of distinct agents per participant. Among the 49 patients, sodium channel blockers were the most common treatment mechanism, prescribed to 32 patients (65%). Additionally, 29 patients (59%) received at least one GABAergic (gamma-aminobutyric acid) ASM, while 14 patients (29%) were treated with at least one calcium channel blockers. The demographics of children with epilepsy are summarized in Table 1 and Extended Data Table S1.

Participants and their legal guardians provided fully informed consent to participate in the study. The experimental procedures were approved by CCMC and North Texas Regional Institutional Review Boards.

Experimental design

Simultaneous MEG and HD-EEG recordings were performed during visual stimuli presentations (Fig. 1A). MEG data were acquired using a 306-sensor (204 planar gradiometers and 102 magnetometers) whole-head system (TRIUX, Elekta Neuromag). Depending on the participant's head size, HD-EEG recordings were obtained using a Geodesic EEG System 400 (Magstim EGI) with either 256 electrodes (n = 77 participants) or 128 electrodes (n = 2 participants). All recordings were conducted at the MEG lab of CCMC. Recordings were performed in a dimly lit, one-layer magnetically shielded room (MSR) with a sampling rate of 1 kHz. Prior to recording, we placed five head position indicator (HPI) coils at known locations on the child’s head (i.e., one on each side of the forehead near the hairline; one on each mastoid bone; and one on the top of the head). The HPI coils were used to accurately determine the child’s head position relative to the MEG sensors. We used a motion tracking device (3SPACE FASTRAK, Polhemus) to digitize the fiducial landmarks, the participants' head shape (capturing over 500 surface points), and the locations of HPI coils. The locations of the positioning coils were recorded before the MEG experiment commenced. More details about the followed experimental procedures are provided previously (Corona et al., 2024).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Data analysis pipeline. A, Simultaneous recordings of high-density EEG (HD-EEG, 256 channels) and MEG (306 channels) were performed in a MSR; visual stimuli consisted of a checkerboard pattern with a diagram character (85 presentations; blurred to comply with copyright rules) at the center of the screen, along with pseudorandom 22 trials of the checkerboard pattern alone per EEG/MEG run with interstimulus interval of 4 ± 1 s. B, Individual MRIs were used to create realistic head models using the boundary element method, which were coregistered with MEG and HD-EEG sensors. Sensor-level and source-level VEFs and VEPs were computed and highlighted components were utilized for further analysis (MEG, M100, M250, M250, M500; HD-EEG, N1, P1, N2, P3). Red and blue traces in the VEFs and VEPs represent the signal envelopes used to identify key components and peak responses. Noise covariance matrices were then computed from VEPs and VEFs separately for MEG and HD-EEG using prestimulus periods. Source maps were reconstructed, and virtual electrodes were implanted at the primary visual cortex (V₁). C, TFR of ESI and MSI at virtual channels was computed using complex Morlet wavelets across 1–100 Hz and 100–250 Hz. TFR maps were generated separately for evoked and induced responses. Evoked TFRs were derived from the average source maps across all trials, whereas induced TFRs were obtained by averaging wavelet transforms of source maps for individual trials. Induced activity was computed in three frequency bands (i.e., 1–30 Hz, 30–100 Hz, 100–250 Hz), with baseline correction applied from −0.2 to 0 s relative to stimulus onset. RPC, relative percentage change.

For the HD-EEG recordings, we placed the Geodesic HD-EEG net on the child’s head and ensured that resistances were <10 kΩ for all electrodes. The coregistration of HD-EEG electrodes relative to the child’s head position was carried out by using the GeoScan Sensor Digitization System (Magstim EGI). During the recordings, participants were seated in a MEG-compatible chair and instructed to remain as still as possible. VEFs and visual-evoked potentials (VEPs) were recorded in response to visual stimuli presented on a screen placed in front of the participants (Fig. 1A). The setup included a MEG-compatible projector (Imedco) connected to a computer. Stimulus presentation was controlled by the Gentask module of the STIM package (NeuroScan). The visual stimuli included 85 age-appropriate captivating diagram characters (i.e., Disney, Minecraft, Airbender, SpongeBob, Mario, Phineas Ferb, Gumball, and Superhero) or faces overlaid on a checkerboard pattern, as well as 22 random trials of checkerboard patterns alone to minimize potential biases and ensure that the responses are not influenced by predictable patterns or sequences (Fig. 1A). To maintain consistency and control for variables known to affect VEFs/VEPs (i.e., contrast, angle, speed, and frequency of visual stimuli), these diagram characters were uniformly presented throughout the experiment for all participants. Although diagram stimuli may introduce some variability in visual responses, their inclusion was essential for maintaining engagement for our participants. To assess whether stimulus differences confounded group-level effects, we conducted separate analyses of both checkerboard-only trials and responses to individual diagram stimuli. For both stimulus conditions, fixation points or diagram characters were centered on the checkerboard to minimize variability from different gaze positions. Stimuli were presented for 1 s, followed by a variable interstimulus interval of 4 ± 1 s, in a pseudorandom order to prevent systematic effects that could arise from the order of presentation, such as anticipatory responses or habituation (Fig. 1A). Participants were instructed to focus on the fixation point or diagram characters to reduce artifacts from eye movements, blinks, and saccades. Foam wedges were placed between each participant's head and the inside of the dewar to provide comfort and to ensure stability and minimize any potential motion. If excessive movements were observed, the session was repeated. A member of the research team was always present inside the MSR during the recordings. The participants were monitored using a MEG-compatible video system to ensure compliance with instructions and to detect any notable artifacts, such as head movements, improper head positioning, or lapses in attention to the screen. The recording lasted ∼7 min per run; the total duration of the actual recording was ∼1 h (including breaks between runs). Electrocardiography (ECG) and electrooculography (EOG) data were simultaneously collected during MEG and HD-EEG recordings at a sampling rate of 1 kHz. ECG data were recorded using two leads placed on the sternum and the fifth intercostal space on the left side of the body, while EOG data were captured with two electrodes positioned above and below the eye on the right side of the face.

MRI data acquisition

High-resolution, T1-weighted structural MRI sequences were acquired using a Siemens Skyra 3 T scanner with a 10-channel head coil. Specifically, a three-dimensional magnetization–prepared rapid acquisition gradient-echo sequence in the transverse plane was acquired with generalized autocalibrating partial parallel acquisition, with the following parameters: repetition time, 2,200 ms; echo time, 2.45 ms; flip angle, 8°; field of view, 250 mm; 176 slices; matrix size, 256 × 256; and acquired resolution, 1.0 × 1.0 × 1.0 mm. MRI data were unavailable for 12 patients and six neurotypical controls (Table 1, Table S1). For these individuals, age-matched MNI152 T1-weighted structural MRI templates were used for further analysis (Grabner et al., 2006).

MEG and HD-EEG data analysis

We analyzed the MEG and HD-EEG data using Brainstorm (Tadel et al., 2011). MEG data were initially processed using temporal signal-space separation to suppress external noise and artifacts from nearby sources enhancing signal-to-noise ratio (Taulu and Simola, 2006). We then combined the MEG and HD-EEG data into a single .fif file by using an in-house Python code that spatially aligns the MEG and HD-EEG sensor locations and synchronize signals and triggers across the two acquisition systems (for details, see Corona et al., 2024). This integration was not performed to merge sensor data but to align and synchronize data streams for parallel analyses. Each modality was then processed independently to preserve its distinct spatial sensitivity and noise characteristics. MEG is primarily sensitive to sources with tangential orientation, while EEG captures sources with both tangential and radial orientation (Chowdhury et al., 2015; Chikara et al., 2023; Bastola et al., 2024). This complementary sensitivity is particularly important in pediatric populations, where higher anatomical variability and immature cortical folding can complicate single-modality source modeling (Ball et al., 2021). Hence, the use of simultaneous HD-EEG and MEG improves spatial resolution and increases the probability of capturing various ranges of source orientation and cortical depths. We further inspected the MEG and EEG data to remove bad channels; any time segments contaminated with movement or muscle tension artifacts were identified through visual inspection and excluded from further analysis. Additionally, biological artifacts, such as heartbeats and eyeblinks, identified through simultaneously recorded ECG and EOG data, were removed using the signal-space projection (SSP) method (Uusitalo and Ilmoniemi, 1997), which is well suited for correcting spatially consistent and time-locked biological artifacts, particularly in evoked paradigms (Atti et al., 2024). Both modalities underwent separate artifact correction procedures (e.g., EOG, ECG removal, channel-wise rejection) to preserve signal integrity before subsequent analysis. We selected SSP over ICA based on its computational efficiency, minimal reliance on assumptions, and its lower risk of removing task-related brain activity in time-locked analyses (Mutanen et al., 2016). Heartbeat and eyeblink artifacts were removed using SSP in a sequential, artifact-specific manner: first, cardiac artifacts were removed using heartbeat-only segments; next, blink artifacts were corrected using blink-isolated data, avoiding cross-contamination. This approach preserved neural signals while minimizing overcorrection. We removed the DC offset and applied a fourth-order bandpass infinite impulse response Butterworth filter with separate frequency ranges of 1–100 Hz and 100–250 Hz, as well as a 60 Hz notch filter and its harmonics (i.e., 120, 180, 240 Hz).

For the visual-evoked and induced cortical oscillations, we segmented each dataset into trials from −200 to 500 ms relative to the stimulus onset. Some MEG/EEG trials were lost due to technical difficulties and signal artifacts, leading to a reduction in the total number of analyzable trials (Table S1). We examined whether the number of artifact-free trials differed between the EP and control groups and found no differences (p > 0.05). These artifact-free trials from MEG and HD-EEG were then used for further analysis. Furthermore, to enhance signal focality and minimize the influence of volume conduction, skull thickness, and other noncortical contributions, HD-EEG trials were processed using common-average referencing (Muthukumaraswamy et al., 2010). In some participants (n = 22), stimulation-related artifacts were observed in the MEG and HD-EEG data. For those participants, we removed the time segment from −10 to 10 ms around stimulation onset and replaced those values with linear interpolations using MATLAB's interp1 function; we excluded this period from further analysis to ensure homogeneity across all participants.

Electric and magnetic source imaging

We performed electric and magnetic source imaging (ESI and MSI) on the averaged data of each participant using the dynamic statistical parameter mapping (dSPM), which is a noise-normalized solution based on the minimum norm estimate (Hämäläinen and Ilmoniemi, 1994) (Fig. 1B). We averaged the source imaging maps across runs for each participant to obtain an individualized average dSPM source model for each participant. We built realistic head models using boundary element model based on each participant's T1-weighted brain MRI scans or age-matched template MRIs to estimate the forward model (Fig. 1B). The models were generated with the OpenMEEG software and consisted of three surface layers (i.e., brain, inner skull, and outer skull). The dSPM solutions were estimated at the cortical surface (∼15,000 sources) reconstructed in Brainstorm using statistical parametric mapping (SPM) functions. The dSPM source solutions were computed separately for EEG and MEG using their respective noise covariance matrices and lead fields. This allowed us to retain the physiological specificity of each modality while using shared anatomical constraints for comparability across participants. The noise covariance was estimated on merged artifact-free prestimulus portions previously selected from the HD-EEG and MEG recordings. Grand-averaged cortical sources for each participant in the EP and control groups were calculated confirming that the maximum cortical sources of VEPs and VPFs were located within the primary visual cortex (V₁). This was further validated by overlaying scouts from the left and right cuneus of the Desikan–Killiany atlas (Desikan et al., 2006) using FreeSurfer (Fischl, 2012). We averaged source maps for each condition and each participant and then projected each source map onto a common 5 mm grid based on the MNI brain (7,743 vertices) for group analysis. For further analysis, we focused on the maximum cortical activation of implanted scout in V₁, calculating mean signals for each cortical vertices in V₁ for each participant (An et al., 2018; Fig. 1B). We then compared the resulting source waveforms between the EP and control groups using an overlapping windowing approach. Specifically, we defined 25 ms windows (after stimulation onset) with 50% overlap, computed the mean amplitude for each window, and performed a “Wilcoxon rank-sum test” to assess group differences, followed by false discovery rate (FDR) correction for each window.

To analyze cortical components (determined from MEG) and peaks (determined from EEG) of cortical activity during visual stimulation in individuals with epilepsy and neurotypical controls, we first identified the upper and lower envelops of VEFs and VEPs at the sensor level from all the channels initially determined with MATLAB function envelop. Then, the maximum of upper envelopes and minimum of lower envelopes across all channels was determined and used to detect peaks and troughs. The MATLAB function findpeaks was applied in the detection of components or peaks in subject-specific grand–averaged VEFs and VEPs, with a minimum peak height constraint set to three standard deviations above the prestimulus amplitude (Song et al., 2024). The identified pairs of peaks and troughs were further evaluated visually. Additionally, we compared the peak amplitude and latency (relative to stimulus onset) of visual cortical components or peaks between individuals with epilepsy and controls in VEFs and VEPs. For VEFs, we analyzed the M100 (40–100 ms), M150 (101–150 ms), M250 (151–250 ms), and M500 (251–500 ms) components (Fig. 1B). For VEPs, in similar time range as VEFs, we examined the N1, P1, N2, and P3 peaks (Shigeto et al., 1998; Meeren et al., 2013). Furthermore, to investigate these components at the source level, we used virtual sensors from the left and right cuneus, defined in the Desikan–Killiany atlas, to reconstruct source VEFs and VEPs (Song et al., 2024). Same parameters were applied to analyze the peak amplitude and latency of the N1, P1, N2, P3 M100, M150, M250, and M500 components at both sensor and source levels (Fig. 1B).

Time–frequency analysis

In both MEG and HD-EEG data, we determined the peak of the initial visual cortical response after stimulus onset for each participant individually. We generated virtual channels at the site of peak cortical activity (Fig. 1B) to enhance source-level accuracy by minimizing head movement artifacts (Velmurugan et al., 2022). These virtual sensors were modality-agnostic in their anatomical placement but derived independently from each participant's source-localized peak activity, as determined from EEG and MEG source models, respectively. Subsequently, we conducted time–frequency analysis of electric and magnetic brain activity at this virtual channel using complex Morlet wavelets to generate time–frequency representation (TFR) maps. To estimate evoked responses, we applied Morlet wavelets with five cycles per wavelet at center frequencies ranging from 1 to 100 Hz (in 1 Hz steps) to the averaged source waveforms (Fig. 1C). For the analysis of the BBG, we used Morlet wavelets with the same five-cycle width at center frequencies between 100 and 250 Hz in 1 Hz steps (Fig. 1C). We then corrected the TFR maps with respect to the baseline activity from −200 to 0 ms to estimate percentage change in power (i.e., relative power change). To estimate visual-induced fields (VIFs) and potential (VIPs) responses, we applied Morlet wavelets (using the same parameters as for evoked responses) to the source waveforms for individual trials from each participant (Fig. 1C). The estimation of TFRs for the induced response was performed using wavelet transforms at three frequency ranges (1–30, 30–100, and 100–250 Hz, separately) to minimize edge effects and 1/f power decay (Kujala et al., 2015). Furthermore, for counteracting 1/f power decay, we applied “spectral flattening” to the TFRs computation in the 1–30 Hz and 30–100 Hz ranges, where power in each frequency bin was multiplied by its corresponding frequency value to balance trade-off between power and frequency precision (Demanuele et al., 2007; Voytek et al., 2015). Without this compensation, low-frequency power could potentially obscure true neuronal responses in 1–100 Hz range (Miller et al., 2009). “Spectral flattening” was applied to the 1–30 Hz and 30–100 Hz frequency ranges and tends to overcompensate for frequencies beyond 100 Hz (Scheffer-Teixeira et al., 2013); thus, we left the 100–250 Hz range unadjusted. We then corrected the TFR maps for baseline activity in each trial across three frequency bands and averaged the baseline-corrected maps separately for each frequency range (Papadelis et al., 2023). Evoked and induced TFRs were averaged for each subject and then grand-averaged for all participants in the EP and control groups (Fig. 1C). Furthermore, to capture true neuronal oscillations in the different frequency bands for each evoked and induced TFRs, we defined oscillations as wavelets with peak amplitudes at least twice that of the background noise. By setting this threshold, we aim to increase the likelihood that the detected oscillations are genuine, minimizing the risk of false positives and improving the accuracy of our analysis (Inoue et al., 2004; Strigaro et al., 2021). Each TFRs for both evoked and induced responses were further evaluated visually by two reviewers (S.R. and F.K.K).

For each participant on subject-specific grand–averaged individual TFRs, we identified the peak frequency, latency (relative to the stimulus onset), and amplitude of the evoked and induced responses, separately for the beta, NBG, and BBG bands (Perry et al., 2014; Kujala et al., 2015). The peak frequency was identified as the frequency showing the largest power change after stimulus onset in each band (i.e., beta, NBG, and BBG); the peak latency was defined as the time–instance after stimulus onset showing the largest power change in each band; the peak amplitude was determined at the time–frequency bin corresponding to the peak frequency and time–instance in each band. Furthermore, we averaged the relative power change within the beta, NBG, and BBG to determine the average power change (in terms of magnitude) for each band; all these measures were used as “TFR features” of evoked and induced responses in subsequent statistical analyses (Fig. 1C). Furthermore, to address intraindividual variability in neuronal responses particularly in average power and peak amplitude, we applied min–max normalization across the participants (controls and EP group) using the formula X_normalized= (X − X_min)  / (X_max − X_min), where X is the original feature value for a participant and X_min and X_max are the minimum and maximum value across participants, respectively (Duma et al., 2024). Additionally, for participants exposed to multiple diagram stimuli, we extracted key response features (i.e., peak frequency, latency, amplitude, and power) and compared them with those from the checkerboard-only stimulus. Trials with features deviating by more than three standard deviations from the checkerboard-only trials were excluded to minimize stimulus-related confounds.

Statistical analysis

We performed statistical analysis of the visual-evoked and induced cortical oscillations data with FieldTrip toolbox in Brainstorm and MATLAB. For the time–frequency analysis, we analyzed the TFRs between subjects (i.e., EP vs control groups). All statistical tests were conducted on the poststimulus window (10–500 ms) using Monte Carlo permutation tests (1,000 permutations), with multiple comparisons corrected via cluster-based permutation testing (Meyer et al., 2021) in FieldTrip. Clusters were formed from neighboring pixels exceeding a p < 0.05 threshold. A reference distribution was generated with 1,000 permutations, and clusters were considered significant if they exceeded the 99th percentile of this distribution. For each participant, we extracted the power of relative responses within the time and frequency window where statistically significant clusters were observed and compared between EP and control groups (“Wilcoxon rank-sum test”). Additionally, we extracted the TFR features (i.e., peak frequency, peak latency, peak amplitude, and average power change of evoked and induced responses) for beta, NBG and BBG frequency bands and compared between the EP and control groups (“Wilcoxon rank-sum test”) followed by FDR correction (p < 0.05). To examine epilepsy subtype-specific effects, we performed subgroup analyses comparing children with focal and generalized epilepsy. TFR features of evoked and induced responses were compared across control, focal, and generalized epilepsy groups using pairwise Wilcoxon rank-sum tests with FDR correction (p < 0.05). Similarly, to address potential effects related to epilepsy pathology, we conducted a subgroup analysis comparing evoked and induced TFRs features between neurotypical controls and patients with nonlesional epilepsy (Table 1). Group differences were assessed using Wilcoxon rank-sum tests followed by FDR correction (p < 0.05). Successively, we performed “Spearman’s correlation” between the EP group TFR features of the beta, NBG, and BBG bands and clinical parameters (i.e., age of onset, epilepsy duration, and the number of ASMs), as well as age of participants during MEG/EEG experiment (correlations coefficient, r; p < 0.05). Furthermore, we performed multiple linear regression to assess the distinct contributions of each clinical parameter and patients' age to neurophysiological features (p < 0.05). Normalized feature values (accounting for interindividual variability) served as dependent variables, with age at epilepsy onset, patient age, the number of ASMs, epilepsy subtype (focal and generalized), and epilepsy duration as predictors (Duma et al., 2024). The model was implemented in MATLAB as follows: model = lm (feature matrix ∼ age of epilepsy onset + patient's age + number of ASMs + epilepsy duration + epilepsy subtype).

Classification model

We implemented a supervised machine learning algorithm to classify participants into specific classes (i.e., epilepsy vs neurotypical controls). More specifically, we trained a support vector machine (SVM) classifier with a radial basis function kernel to distinguish participants based on their neurophysiological response patterns to visual stimulation. The performance of the classifier was evaluated using stratified 10-fold cross-validation, confusion matrix metrics, receiver operating characteristic (ROC) curve analysis, and significance testing by χ² test (p < 0.05). To reduce overfitting, we applied a stratified 80/20 train–test split repeated over 1,000 iterations. In each iteration, models were trained on 80% of the data and tested on the remaining 20%, ensuring independence of the test set and unbiased estimates of model performance. We trained separate SVMs for three participant groups based on their recording modalities, using 80% of the data for training and the remaining 20% for testing: (1) children with both HD-EEG and MEG data available (Group 1), (2) children with only MEG data available (Group 2), and (3) children with only HD-EEG data available (Group 3). Each model used a regularization parameter of 10. For each group, we performed ANOVA on all normalized features that showed group-level differences between the EP group and neurotypical controls, including the peak amplitude and latency of VEP peaks (e.g., N1) and VEF components (e.g., M100, M150, M250). In addition, we analyzed normalized TFR features such as average power change, latency, and amplitude of both evoked and induced responses in the beta, NBG, and BBG bands. Furthermore, we retained only the features that showed statistically significant differences between classes before training the classifiers (Partamian et al., 2025). The ANOVA test provided p values for each feature, adjusted using the FDR correction. Significant features were identified based on these adjusted p values with importance scores calculated as the negative logarithm of the p value −log₁₀(p). To mitigate potential bias from ASM use, we excluded features associated with ASMs in the EP group. The remaining significant features served as predictors, with class labels as the target. SVMs were validated using a stratified 10-fold cross-validation, ensuring balanced class representation in each fold and robust performance estimation across the dataset. To further strengthen reliability, this process was repeated across 1,000 random 80/20 stratified train–test splits, allowing us to assess the stability and generalizability of the model's performance across different data partitions. For each iteration, we defined true positive (TP) as correctly classified patients with epilepsy, true negative (TN) as correctly classified neurotypical controls, false positive (FP) as neurotypical controls incorrectly classified as patients with epilepsy, and false negative (FN) as epilepsy patients incorrectly classified as neurotypical controls. A confusion matrix was created for each group, and the model's performance was assessed using key metrics: (1) sensitivity = TP/(TP + FN); (2) specificity = TN / (TN + FP); (3) positive predictive value (PPV) = TP /  (TP + FP); (4) negative predictive value (NPV) = TN/(TN + FN); and (5) accuracy = (TP + TN) / (TP + TN + FP + FN). ROC curves were then generated to evaluate model's performance by calculating the area under the curve (AUC) for each group, with Youden's index identifying the optimal balance between TPs and FPs. Additionally, Chi-squared test (χ²) was utilized to assess the independence of predicted and actual classifications for all three groups.

Moreover, we used bootstrap-based method to statistically compare the AUCs for all three group models to determine whether the model performance was statistically different across the three groups (Pang et al., 2023). For this, we drew bootstrap sample (n = 10,000 times) from pool of AUC values for all three group models (1,000 iterations) and computed difference in AUC between all three models to build empirical difference distribution which was compared between groups for significance differences (p < 0.05). Moreover, to evaluate whether combining features from multiple modalities improves classification performance, we performed logistic regression on MEG and EEG from participants with simultaneous recordings (Group 1). First, we computed a base model using either EEG-only or MEG-only features from Group 1. Then, we built a full model that included both EEG and MEG features to evaluate whether the addition of the second modality significantly improved classification performance (Taneichi et al., 2011). Finally, a likelihood ratio test was performed to compare model deviances between MEG/EEG-only model and full model (p < 0.05) to evaluate whether integrating features across modalities provides additional predictive value compared with using a single modality alone.