An automatic MEP selection algorithm was written to determine whether a given recording contained an MEP [12]. To remove stimulation artifacts from the train of 5, we excluded the first 400 data points (corresponding to 20 milliseconds). Then, two features were computed: the onset latency and duration of the signal. Onset latency is defined as 1 millisecond (empirically determined) before the trace crosses the mean of the baseline (the last 5 milliseconds of the recording) plus or minus the SD of the entire recording. The end of the signal was calculated similarly, by starting from the end of the recording. We defined the duration as the end of signal latency minus onset latency. In addition, the time interval between the first and last peak was determined using the scipy.signal function find_peaks. A recording was considered to contain an MEP if at least one peak was detected, the duration was less than 40 milliseconds, and the interval from the first to last peak was less than 35 milliseconds (in accordance with clinical experience).

In our analysis pipeline, we used three distinct representations of MEP data (see Figure 1B): time, feature, and time-frequency, each tailored to optimize the performance of our ML classifiers.

A finite-impulse response bandpass filter with 30- and 1000-Hz cutoff frequencies was applied to the 1600-dimensional signal vector. These frequency settings align with MEP visualization practices using the monitoring machine at center T₁ (ie, the software filters). The signal vectors are then normalized with respect to the absolute maximum MEP value in each patient. This filtered and normalized time representation of the data was used to train, validate, and test a random forest (RF) classifier, as well as a 1-dimensional convolutional neural network (1D-CNN).

We used a customized feature extraction algorithm to condense each filtered and normalized MEP signal into characteristic features. The initial choice of predictors is a combination of clinically used predictors (eg, latency, amplitude, minimum, and maximum) and routinely used features in general neurophysiological literature (spectral entropy, frequency, etc). After a correlation analysis (see Multimedia Appendix 1), we chose five features that showed no correlation describing relevant domains of the signal: peak latency, maximum signal value, number of peaks, main frequency, and slope.

Peak values were extracted with the scipy find_peaks function with peak prominence defined as twice the SD of the signal. The main frequency was calculated as the frequency at which the Fourier transform of the filtered data attained its maximum absolute value. Finally, the slope of the signal was computed as the mean of the first derivative of the signal (using the NumPy function gradient). The resulting 5-dimensional feature representation of the data was then used to train, validate, and test a RF classifier.

Finally, employing the Python library PyWavelets (pywt), a continuous wavelet transform with a Mexican hat mother wavelet was applied to the data to transform them into 2D time-frequency representations. Scales ranging from 2 to 30 were logarithmically spaced, while the time dimension was undersampled to yield an array of dimensions 224×224. These 2D time-frequency representations were then used to train, validate, and test a 2D-CNN.

Using custom Python scripts, Student t tests were applied to compare the mean values of two different features. Specifically, for each muscle, the means of the different features were compared across the two centers. Statistical significance was set at P=.05.

The Python library scikit-learn [21] was used for the RF classifier, while tensorflow [22] and keras [23] were used to obtain the 1D- and 2D-CNN models. Hyperparameter tuning was carried out for the RF and 1D-CNN, while we used fixed parameters for the 2D-CNN (see below). The hyperparameters used in the grid search of the RF are the same as described in Multimedia Appendix 1. The architecture of our 1D-CNN is inspired by the model of Ahmed et al [24], and consists of two consecutive 1D convolutional layers, followed by a MaxPooling layer, a dropout layer, and a BatchNormalization layer. The model then gets flattened, before adding a dense layer and another dropout layer and finally ending in a dense output layer. The structure of our 2D-CNN is inspired by the model of Wang et al [25]. It consists of 4 blocks of 2D convolutional layers followed by BatchNormalization layers, with a MaxPooling layer after the second and third blocks and a GlobalMaxPooling layer after the fourth block. The model is capped off by a dense output layer. The specifics of these models are shown in Multimedia Appendix 1.

The dataset from Center T₁ was split into 70% for training and 30% for validation (stratified according to patients), while the whole of the dataset from Center T₂ was used for testing. In all cases, we used class weighting [26] to deal with the class imbalance problem (ie, the number of leg muscle MEPs is lower than the number of arm muscle MEPs, see Wermelinger et al [12] for a discussion of this issue).

All prediction models output probabilities of belonging to the predicted class. The decision thresholds were set at a chance level of 0.25 (likelihood of belonging to 1 of 4 classes of muscles). They were systematically explored and reported (see Figure 2B) for decision thresholds up to 0.9.

Accuracy was used as the primary performance metric, and the confusion matrix was used to evaluate the performance of the classification algorithm. The outcome assessment does not require subjective interpretation, since the muscle identity is objectively assessable, independent of sociodemographic background and clinical outcome.

No model updating or recalibration was performed during the model evaluation. While some variability in model performance was observed across different centers (Figure 2), we opted to retain the original model without adjustments. Future work may explore model updating to enhance performance in these areas.

To elucidate how the RF classifiers made their decisions, we used feature importance and Shapley Additive Explanation (SHAP) values. Feature importance values are provided by the feature_importances_ attribute of the scikit-learn RandomForestClassifier class, while the SHAP values are calculated with the SHAP library [27]. The feature importance values are determined by aggregating (mean and SD) the impurity decrease within each decision tree. SHAP values quantify the impact of each feature on prediction outcomes. Positive values signify a positive influence, while negative values indicate the opposite, with magnitude representing the strength of the effect. SHAP values of the RF on feature representation of the MEP data were computed on a random sample containing 20% (5919 MEPs) of the training data set. In the case of CNNs, we used gradient-weighted class activation mapping (Grad-CAM). This is a type of attention map, a visualization tool highlighting regions within an image considered by the neural network to be pivotal for specific predictions [28]. We adapted preexisting code to generate Grad-CAMs for both 1D- and 2D-CNNs [29,30]. The Grad-CAMs of all signals in the training data set were averaged to obtain the corresponding plots for the 1D-CNN (overall) and 2D-CNN (for each muscle).

Although algorithms using explicit feature representations showed poorer performance (80% test accuracy), they provided key insights into the decision-making process. Feature importance analysis (Figure 4A) revealed that peak latency, a primary factor in clinical decision-making, was the main driver of classification. Main frequency was the next most important, despite typically not being used by clinicians. This was followed in order of importance by maximum signal value, while slope and number of peaks were the least important features.

SHAP values from the RF model using feature representation elucidated how different parameters influenced the model’s decisions for each class (Figure 4B). In all four muscles, peak latencies were crucial. For the upper extremity, short latencies favored correct predictions, whereas long latencies were accurate indicators for the lower extremity classes. Interestingly, the main frequencies did not exhibit this pattern. High main frequency values favored predictions for distal muscles (APB and AH), whereas low frequencies were associated with proximal muscles (EXT and TA).

When presented with complex data, such as time-series and wavelet transforms, ML algorithms use internal feature representations to guide their decision-making. For CNNs, these features can be visualized using attention maps, which highlight the areas of the input data that are most salient and decisive for classification. In the depicted Grad-CAMs, these areas correspond to where the signal occurs (see Figure 5B). The insights from explicit feature representation regarding main frequencies are also evident in the Grad-CAMs, since attention for proximal muscles focuses on lower frequencies compared with distal muscles. Similarly, in the feature importance analysis of the RF using the time representation, and in the average (over all signals and all muscles) of the Grad-CAMs of the 1D-CNN, the highest importance is assigned to the time points when the signals occur (see Figure 5A).

Our study demonstrates the high performance of ML models, such as RF, 1D-CNN, and 2D-CNN, in classifying MEPs recorded during IONM. Notably, the RF classifier achieved 87.9% validation accuracy and 80% test accuracy using time representation data, while the 1D-CNN and 2D-CNN achieved comparable performances with slightly increased variations in accuracy across different datasets.

Furthermore, our analysis revealed that frequency is a critical feature that these algorithms use for decision-making, with different frequency ranges (low vs high frequencies) being decisive depending on the muscle group involved (proximal vs distal, respectively). This is an important finding, which should encourage clinicians to investigate this feature for potential warning criteria. Since there are still disagreements when it comes to warning criteria during intraoperative monitoring of motor evoked potentials, our results may already provide an opportunity to increase patient safety during surgical procedures.

Significant differences in MEP latencies between datasets from centers T₁ and T₂ (see Figure 3A) may highlight the influence of different stimulation techniques. The higher number of MEPs induced by DCS at center T₁ might explain this difference, but other factors might also influence these findings, such as data collection methods, types of surgical procedures, and characteristics of the selected patient population including height, age, disease, etc. For example, the higher proportion of upper extremity MEPs in the data from center T₁ is due to different surgical focuses at center T₁ compared with center T₂. Understanding these differences is crucial for interpreting ML model performance and ensuring generalizability across centers. Including more centers and more extensive data collection can mitigate biases and advance research in the field.

Our findings indicate significant differences in MEP frequencies between distal and proximal muscle groups. Distal muscles exhibit higher MEP frequencies compared with proximal muscles, a trend that was effectively used by various of our tested models in their decision-making processes. The underlying neurophysiological mechanisms contributing to these differences are not fully understood. Although the general pathway of distal and proximal MEPs are similar—upper motor neurons synapsing on lower motor neurons that innervate muscle fibers at the neuromuscular junction—anatomical and physiological differences between distal and proximal muscles may explain the observed frequency variations. The following physiological and anatomical characteristics outline potential factors contributing to these differences.

First, distal muscles, involved in fine motor control, possess a higher density of smaller motor units compared with proximal muscles. The smaller motor units of distal muscles have lower activation thresholds but generate less force than the larger motor units found in proximal muscles [31-33]. Furthermore, distal hand muscles contain a greater proportion of slow motor units, which are more fatigue-resistant [34,35].

Secondly, the temporal dispersion of electrical activity differs between muscle groups. Distal muscles, such as those in the hand and foot, exhibit more synchronous and temporally concentrated MEP responses, whereas proximal muscles display greater temporal dispersion. This increased synchrony in proximal muscle MEPs likely contributes to the higher frequency distribution observed in distal muscle MEPs.

Thirdly, the cortical representation of distal muscles is significantly larger than that of proximal muscles, reflecting the dense corticospinal innervation of these areas. Hand muscles receive among the strongest corticospinal inputs, highlighting their critical role in precise motor control [36-38].

In addition, various motor control pathways interact differently with distal and proximal muscle groups, further influencing the MEP frequency characteristics. These interactions likely involve contributions from both corticospinal and other descending motor pathways, though their exact contributions require further investigation [39].

The functional relevance of high and low-frequency bands within MEPs remains uncertain. While our findings suggest that MEP characteristics are largely determined by muscle-specific neurophysiology, it is essential to consider the potential role of top-down modulation from cortical and subcortical regions. These central mechanisms could influence observed frequency differences and may contribute to the observed MEP variations between distal and proximal muscles.

Given that intraoperative changes in MEPs are considered critical warning signs of upper and lower motor neuron impairment, further investigation is necessary to clarify the relationship between MEP frequency components and both muscle neurophysiology and neuronal modulation mechanisms.

ML models, particularly those designed to handle complex data, consistently achieve high performance, but their lack of transparency can hinder interpretability and therefore acceptance in clinical processes. Specifically, we attained 80% accuracy with our five features per signal, compared with 87% accuracy with 1600 data points per signal. To address interpretability, we applied two complementary explainability techniques: SHAP and Grad-CAM, each offering unique insights into model behavior.

SHAP provides a feature attribution approach, assigning precise numerical contributions to each input feature—such as latency or main frequency—to quantify its role in the prediction process. This method is independent of the choice of ML model and excels at identifying the relative importance of features and offers consistent, interpretable insights, albeit being computationally expensive. In contrast, Grad-CAM generates attention maps that visually highlight which regions of the input signal influenced the model’s decisions. They are model-specific and generally limited to CNNs. These visualizations are particularly useful for validating whether the model focuses on relevant, clinically meaningful areas of the MEP signal.

The combination of SHAP and Grad-CAM allowed us to cross-validate findings, ensuring that the observed importance of main frequency was both quantitatively consistent and visually evident. Compared with other interpretability techniques, such as LIME [40]. SHAP provides more robust and consistent explanations by attributing specific contributions of each feature to individual predictions (local interpretability) while also summarizing feature importance across the entire dataset to reveal overall model behavior (global interpretability).

Our complementary approach highlights the importance of main frequency as a decisive MEP feature and demonstrates how XAI methods can uncover meaningful insights that warrant further testing in basic research and clinical trials. Moreover, the trade-off between explainability and complexity must be carefully considered. In intraoperative settings, where clinical trust is crucial, Grad-CAM’s intuitive attention maps may be favored for transparency. Conversely, SHAP’s precise attributions offer deeper insights for research applications requiring feature-level understanding.

Ultimately, the choice of method depends on the context. Highly accurate yet less interpretable models should not be dismissed if extensively validated on diverse populations. In practice, balancing explainability and accuracy through complementary methods ensures optimal utility, whether for guiding clinical decisions or advancing research.

Accurately identifying muscles and avoiding labeling mistakes is critical in IONM, with previous studies highlighting the consequences of such errors that harm the patients we seek to protect [41,42]. During the IONM setup, mislabeling of muscles has caused false negative alarms, in missing MEP alterations of a presumed unaffected muscle. These incidences have caused potentially avoidable motor deficits during surgery and consequently resulted in legal actions. Nevertheless, we have to acknowledge that the IONM setup in the operating room environment is prone to errors as it is a high-pressure environment [43]. Safety checklists have been implemented; however, an automated ML safety check would increase the avoidance of mislabeling. Those algorithms may be implemented in an existing IONM software.

Further, our results suggest that expanding the search for warning criteria to the frequency domain is essential, as different muscles may require tailored approaches. Once these muscle classification models have been validated on more data and more centers, they could be implemented as a safety mechanism at baseline recordings in surgeries.

When analyzing the confidence of our algorithms (see Figure 2B), it became apparent that signal quality and recording modes might limit high-confidence decisions, even with optimal algorithms. This trade-off between data volume and accuracy necessitates either better-trained models or higher-quality recordings. Investing in improved surveillance methods or stable recording techniques, such as averaging or selecting the best MEPs from multiple recordings, could be an essential step. This might affect how MEP monitoring will be done in the future.

Ensuring the trustworthiness of AI involves addressing ethical and legal implications and incorporating decision confidence metrics could bolster acceptance of AI. The question of responsibility is important in a clinical setting and a transparent decision process for any potential implementation of AI is crucial in this regard. This becomes even more evident when discussing legal aspects and accountability. Integrating ML models with robust explainability attributes has the potential to enhance decision-making accuracy. Disclosing the basis of the algorithmic decisions to the neurophysiologists is key, as it allows them to reason how their understanding differs from the algorithm and oversee the intraoperative decision. In our particular MEP muscle identity scenario, it can provide a safety mechanism against muscle mislabeling and facilitate reliable clinical decisions. By elucidating the prediction bases, XAI supports understanding and trust in AI decisions, which is crucial for the seamless implementation of ML tools in real-time surgical environments. This could significantly increase acceptance of AI and its utility in clinical contexts.