Motor units exhibit more complex action potential waveforms than neurons.

(a) Example recording with a Myomatrix array in a rat triceps brachii (lateral head) during the stance phase of a single stride cycle during treadmill locomotion. The red dashed line indicates the beginning of stance. Channels were spaced in the muscle by 1 mm. (b) Grey region expanded from panel (a), with 4 unique MUAP shapes identified by colored highlights. (c) Extracted mean waveforms from each of the identified MUAPs in this recording. Grey dashed lines illustrate temporal delays across recording channels. (d) Neuropixels 1.0 recording from primary motor cortex layer 5 in a mouse during a decision-making task. The red dashed line indicates movement onset. Channels were spaced in the brain by 40 and 120 μm. Data and resources for this panel are from (International Brain Laboratory et al., 2024; Kristensen & Jörntell, 2023). (e) Grey region expanded from (d), with 3 unique neural waveforms identified by colored highlights. (f) Extracted mean waveforms from this recording. Note the significant conduction delay across channels in motor unit data (c) but not neural data (f) as well as the greater waveform complexity of MUAPs as compared to neural action potential waveforms.

Electrode placement within a muscle determines relative time delays of cross-channel signals for a given MUAP.

By comparing cross-channel correlations across a range of delays, motor unit action potential (MUAP) time alignment can be performed to improve spike sorting. (a) Schematic depicting Myomatrix electrode array and MUAPs propagating along muscle fibers. (b) EMG data from a rat triceps during treadmill locomotion. Four channels from the original recording are shown by red dashed lines. The black lines show the data after the channels have been shifted. (c) Grayscale map of temporal delays that maximize channel correlations for each reference channel. For each row, unshifted reference channels lie along the yellow diagonal. The shading in each square indicates the best channel shift to apply to maximize the correlation with the reference channel, with lighter shades representing higher correlations. The black rectangle shows the best choice of reference channel and shifts that produced the highest sum of correlations. Colored squares show the true position of channels taken from the full recording for use in this figure. (d) Correlation values for each of the 4 selected channels. Peak correlations are marked by colored triangles, indicating the channel time shift that maximized the correlation with the reference channel.

Increasing the number of temporal components provides the capacity to reconstruct a wider range of complex MUAPs across channels.

Temporal components are shown on the left, which are effectively combined during template matching to represent all spike shapes in the dataset. The shape for each motor unit and channel uses the same temporal components but different values for U. Additional temporal components are shown in magenta. (a) Example mismatch between the true multiphasic spike shapes and the reconstruction when using only 6 temporal components with MUAP data. (b) Example of better fitted MUAP spike shapes after increasing the number of temporal components from 6 to 9.

EMUsort provides improved methods for template initialization.

(a) Time series graph showing the default procedure KS4 uses for finding non-overlapping spike examples to use as data for creating “simple templates” (see main text). The non-overlapping crossings (green highlights) are included in formation of the simple templates, but the overlapping crossings are rejected. (b) Depiction showing the use of a single threshold (as in KS4) to extract waveforms. The blue threshold, T1, extracts several examples of MU 1 and 2 (colored yellow and orange), for use with template creation, but none of MU 3 due to the isolation criteria. (c, d) Using additional thresholds (T2; new in EMUsort) captures waveforms of various amplitudes (e.g. larger, more crowded MUAPs like MU 3, colored brown). (e) PCA projections of the 61-point waveforms selected by using all thresholds combined, with the same color coding. (f) All PCA projected crossings from a real rat triceps brachii recording during treadmill locomotion (PCA computed using all isolated waveform shapes from the recording). Using K-means, the cluster centers are identified and used for simple template shapes. The bottom plot shows the 61-point waveforms corresponding to each cluster center. (g) Outlier removal in EMUsort with the same biological data. Low density, outlier spikes are removed with HDBSCAN, and the remaining spikes are clustered as normal with K-means., where outlier removal has resulted in a cleaner set of simple template shapes.

Qualitative sorter performance comparison of biological data and a representative example of simulated data.

Multichannel MUAP spike sorting can be difficult to validate by visual inspection, necessitating quantitative methods for performance assessment. Here we show the sorting of biological data and a sample of the multichannel simulation created from the real dataset. (a) Real EMG data recorded from triceps brachii in a rat during treadmill locomotion. (b) Event plot showing the visually “best” spike sorting result from EMUsort. (c) Example of simulated EMG “data” using simulated spike times (see text).

EMUsort achieves high performance with rat-based simulations using real MUAP recordings from triceps brachii during locomotion.

(a) Sort accuracy for EMUsort and KS4 across a parameter sweep with 25 combinations. Sort accuracies from four MUedit runs are also shown. Each point shows the mean accuracy across ground truth units for each sort run. (b) Same metrics as (a) but the accuracy calculations were computed only data epochs that included highly overlapped spikes, i.e., only MUAPs within at least 2 ms of another MUAP. This assesses how each sorter performs in the most challenging examples in the simulated dataset. In (a, b), the distribution plot is created with kernel density estimation for visualization. (c) Per-unit performance metrics from the best sort for each sorter using all spikes in the rat dataset, arranged from highest to lowest spike count. (d) Per-unit metrics for the same best sorts as in panel (c) but computed only on overlapping MUAPs as in panel (b). (e) Scatter plot showing how well the best sort from each sorter matches the true spike counts, where paired counts for each ground truth unit are shown as a point. The dotted line indicates where the sorted count should match the true count. Note that panels (e, f) reflect total spike counts regardless of sorting accuracy. (f) Same metrics as in (e) but applied only to data with overlapping MUAPs.

Sorter accuracy statistics across multiple spike sorting runs with a rat-based simulation.

Runs were performed across a range of different parameters and the accuracies were computed using all MUAPs. The bracketed values correspond to results using only the overlapping MUAPs. Data points in these distributions are mean accuracies across units in each sort.

Sorter accuracy statistics across multiple spike sorting runs with a monkey-based simulation.

Runs were performed across a range of different parameters, and the accuracies were computed using all MUAPs. The bracketed values correspond to results using only the overlapping MUAPs. Data points in these distributions are mean accuracies across units in each sort.

Parameter sweep values for each sorter with the rat simulated dataset.

With EMUsort and KS4, each pair of Th_universal and Th_learned values were used with each set of Th_single_ch values for a total of 25 parameter combinations. For explanations of Th_universal, Th_learned, and Th_single_ch, see “Kilosort4 spike sorting details”. With MUedit, each Initialization setting was used with each Contrast Function setting listed in the table, totaling 4 parameter combinations.

Parameter sweep values for each sorter with the monkey simulated dataset.

With EMUsort and KS4, each pair of Th_universal and Th_learned values were used with each set of Th_single_ch values for a total of 25 parameter combinations. For explanations of Th_universal, Th_learned, and Th_single_ch, see “Kilosort4 spike sorting details”. With MUedit, each Initialization setting was used with each Contrast Function setting listed in the table, totaling 4 parameter combinations.

EMUsort achieves high performance with monkey-based simulations using real MUAP recordings from biceps brachii during reaching.

(a) Graph of sort accuracy distributions, for EMUsort and KS4 across a parameter sweep with 25 combinations. Sort accuracies from four MUedit runs are also shown. Each point shows the mean accuracy across ground truth units for each sort run. (b) A similar graph as panel (a), but the accuracy calculations were computed only for highly overlapped spikes, i.e., only MUAPs within at least 2 ms of another MUAP were included in the calculations. This assesses how each sorter performs in the most challenging examples in the simulated dataset. In (a, b), the distribution plot is created with kernel density estimation for visualization. (c) Per-unit performance metrics from the best sort for each sorter using all spikes in the rat dataset, arranged from highest to lowest spike count. (d) Per-unit metrics for the same best sorts as in panel (c) but focuses only on overlapping spikes, as in panel (b). (e) Scatter plot showing how well the best sort from each sorter matches the true spike counts, where paired counts for each ground truth unit are shown as a point. The dotted line indicates where the sorted count should match the true count. Please note that panels (e, f) only reflect spike counts, and do not reflect accuracy. (f) A similar scatter plot, where MUAPs are only counted if they overlap with other MUAPs.

EMUsort configuration parameters.

List and explanation of important EMUsort parameters that control its operations. Bolded parameters were introduced in EMUsort, whereas items that are not bolded already existed in KS4, but are made accessible in the EMUsort configuration file for greater usability and parameter sweeping. For shared parameters, default values are for both EMUsort and KS4, unless otherwise noted. Note that the parameters marked “None” or “Results” for the Section field only appear in the copy of emu_config.yaml automatically saved to a sort output folder after a sort run is completed and are simply appended to the bottom for record keeping and usage in future analyses.

Descriptive statistics for distributions of cross-unit mean accuracy scores for each sort run with the rat simulated dataset, using the definition of accuracy as computed in (Pachitariu et al., 2024).

Values in brackets were computed using only MUAPs that overlapped within half the template width of another MUAP.

Descriptive statistics for distributions of cross-unit mean error rates for each sort run with the rat simulated dataset.

This is computed by taking the accuracy scores from Table 1 - table supplement 1 for each sort run and subtracting them into 1. Values in brackets were computed using only MUAPs that overlapped within half the template width of another MUAP.

Per-unit performance statistics from the best EMUsort run, computed with all spikes from the rat simulated dataset.

Per-unit performance statistics from the best sort run with EMUsort, computed only with overlapping spikes from the rat simulated dataset.

Per-unit performance statistics from the best KS4 run, computed with all spikes from the rat simulated dataset.

Per-unit performance statistics from the best sort run with KS4, computed only with overlapping spikes from the rat simulated dataset.

Per-unit performance statistics from the best MUedit run, computed with all spikes from the rat simulated dataset.

Per-unit performance statistics from the best sort run with MUedit, computed only with overlapping spikes from the rat simulated dataset.

Descriptive statistics for distributions of cross-unit mean accuracy scores for each sort run with the monkey simulated dataset, using the definition of accuracy as computed in (Pachitariu et al., 2024).

Values in brackets were computed using only MUAPs that overlapped within half the template width of another MUAP.

Descriptive statistics for distributions of cross-unit mean error rates for each sort run with the monkey simulated dataset.

This is computed by taking the accuracy scores from each sort run and subtracting them into 1. Values in brackets were computed using only MUAPs that overlapped within half the template width of another MUAP.

Per-unit performance statistics from the best EMUsort run, computed with all spikes from the monkey simulated dataset.

Per-unit performance statistics from the best sort run with EMUsort, computed only with overlapping spikes from the monkey simulated dataset.

Per-unit performance statistics from the best KS4 run, computed with all spikes from the monkey simulated dataset.

Per-unit performance statistics from the best sort run with KS4, computed only with overlapping spikes from the monkey simulated dataset.

Per-unit performance statistics from the best MUedit run, computed with all spikes from the monkey simulated dataset.

Per-unit performance statistics from the best sort run with MUedit, computed only with overlapping spikes from the monkey simulated dataset.

EMUsort leverages open-source libraries with a modular organization.

The legend at the top indicates the color coding and symbolic meanings of the various boxes and fonts. The arrows indicate the order and hierarchy of function calls. The main script for the command line interface, “emusort.py”, utilizes a configuration file which stores all settings defining the input dataset properties, sorting parameters, settings defining separate groups of channels to use, and finally KS4 and SpikeInterface parameters. SpikeInterface methods are leveraged for loading, preprocessing, and starting KS4 processes in series or as parallel jobs with different settings. The modifications to KS4 source code, indicated by the legend, were implemented to remove channel delays and improve the initialization of templates.

The composite scores produced by EMUsort for 2 of our simulated rat datasets can serve as a rough estimate of the computed ground truth accuracies (averaged across all units for each sort) with a weak positive correlation for Dataset 1 and a strong positive correlation for Dataset 2.

Dataset 1 was used for the main results generated for Figure 6. Pearson’s R value is shown in the legend of each graph.

Example time stretch in the rat simulated dataset showing the ground truth versus identified spike times.

Each row shows spike times with uniquely colored raster ticks for each motor unit. Darker ticks for each color show the ground truth spike times. Green diamonds mark correct times within +/- 1 ms. Yellow X’s mark false negatives. Red X’s mark false positives. (a) Event plot showing EMUsort evaluations. (b) Event plot showing KS4 evaluations.