Introduction

Motor units transduce the neural signals that originate from supraspinal centres, spinal circuits, and sensory systems into force (1). Each action potential propagating along the axon of an alpha motor neuron elicits action potentials in all its innervated muscle fibres. As such, the hundreds of muscles fibres that belong to each motor unit amplify the discharge activity of spinal motor neurons. The unique spatial location of these muscle fibres facilitates the identification of separate motor unit spike trains from complex electromyographic (EMG) recordings. Decoding the discharge characteristics of a population of motor units provides direct information about the neural control of movement (2, 3). Researchers classically identify a few motor units during relatively weak contractions using concentric needle or fine wire recording electrodes (4) and separate the overlapping spike trains with a template-matching algorithm (e.g., (5)). Recent developments of intramuscular (6–9) and surface (3, 10) EMG electrode arrays facilitate both a larger recording zone (i.e., the area from which motor unit action potentials will be recorded), and the recording of the same motor unit action potential across multiple channels. In conjunction with this hardware advance, the development of novel EMG decomposition software/programs, such as spike-sorting (11–14) and blind source separation (15–19) algorithms, enables a relatively large number of individual motor units to be decoded from each recording (7–10).

Blind source separation algorithms invert the generative model of EMG signals and identify motor unit spike trains from the interference EMG (15–19). In short, this class of algorithms retrieves the spike trains by optimizing a set of separation vectors, i.e., motor unit filters based on the action potentials of each motor unit recorded across an array of electrodes, which maximizes the sparseness of the spike trains (15–19). However, most of the current implementations of this approach rely on offline processing, which restricts its ability to be used for neurofeedback and human interfacing technologies.

Recent studies have reported real-time capabilities of motor unit identification by adapting the offline blind source separation algorithm (20–24). These studies used a two-step approach: 1) the separation vector for each motor unit is identified with offline decomposition during the training phase; and 2) the same vectors are applied in real-time to new EMG recordings. The real-time identification of motor units is only used by a few specialized research teams and there are no publicly available algorithms for this task. In addition, the accuracy and boundary capabilities of online decomposition have not been systematically tested. Such information is necessary to better design experimental paradigms demonstrating, for example, the neural constrains on human movement generation and control (21).

Here, we provide an open-source software that can be used to visualise and track motor unit discharge activities in real-time. We document the software capability boundaries on data collected during isometric contractions of varying force levels, using either grids of surface electrodes or intramuscular electrode arrays from five lower limb muscles (gastrocnemius lateralis and medialis, vastus lateralis and medialis, and tibialis anterior). The accuracy of real-time identification of motor neuron discharge activity was determined from the rate of agreement calculated between the motor neuron spike trains identified in real-time and those identified offline after manual editing. Of note, we did not assess the capabilities of the software with simulated signals because the approach we used for identifying motor neurons has previously been validated for offline analysis (18, 20). Data, codes, and a user manual are available at https://github.com/simonavrillon/I-Spin.

Materials and methods

Overview of the approach

An EMG signal represents the sum of trains of action potentials from all the active motor units within the recorded muscle volume (Figure 1A). During stationary conditions, e.g., isometric contractions, the train of motor unit action potentials can be modelled as the convolution of series of discrete delta functions, representing the discharge times, and motor unit action potentials (Figure 1A). The EMG decomposition with blind source separation consists of inversing the generative model of EMG signals by estimating the series of discharge times for each active motor unit (15–17). When EMG signals are recorded with an array of electrodes, the shape of the recorded potential of each motor unit differs across electrodes. This is due to the varying conduction velocity of action potentials among the muscle fibres, and 2) the location/depth of the muscle fibres that belong to each motor unit relatively to the electrodes, which impact the low pass filtering effect of the tissue on the recorded potential. Therefore, the EMG signals recorded with an array of electrodes can be considered as an instantaneous mixture of the original motor unit spike trains and their delayed versions.

Figure 1.

Increasing the number and spatial distribution of recording electrodes increases the likelihood that each motor unit will have a unique motor unit action potential profile (shape), i.e., a temporal and spatial profile that differs from all the other active motor unit within the recorded volume (10, 25). The uniqueness of motor unit action potential profiles is necessary for the blind source separation to accurately estimate the motor unit discharge times. Our software uses a fast independent component analysis (fastICA) to optimize the separation vector (i.e., the motor unit filter) for each motor unit (Figure 1B; (18–20)).

Overview of the software

The software was coded as a MATLAB app (version 2021a, The MathWorks, Inc, USA). It allows researchers to record and process signals from one to four grids of surface electrodes or intramuscular arrays with up to 64 channels. Of note, the current version has been developed to interface with only one commercialised multichannel EMG system (EMG-Quattrocento, 400 channel EMG amplifier; OT Bioelettronica, Italy). However, the blind source separation algorithm in the code can be used with all systems that record EMG signals with grids of surface electrodes or intramuscular arrays of electrodes. As the accuracy of the algorithm relies on the consistency of motor unit filters, it is recommended to record these EMG signals during stationary conditions - e.g., isometric contractions – to limit changes in muscle geometry or position/orientation of the active muscle fibres relative to the electrodes (26, 27). The framework to perform real-time identification of motor neuron activity has four steps. First, the EMG signal is recorded while participants perform a contraction at the requested intensity such that an electrode mask is manually generated to remove channels with artifacts or low signal-to-noise ratio. This electrode mask is then used for the rest of the experimental session. Second, the force offset is measured and removed before performing MVC. The measured MVC is used to standardize all the submaximal isometric contractions. Third, a baseline contraction is performed at a force level close to the intensity of the testing task and the separation vectors (motor unit filters) are identified using offline blind source separation of the EMG signals. Fourth, the separation vectors are applied over incoming segments of EMG signals during a test contraction to identify motor unit discharge activity in real-time. The algorithm used to identify motor units discharge activity is based on that proposed by Negro et al. (18) and Barsakcioglu et al. (20).

Four forms of feedback can be displayed to the participant: a raster plot of the discharge times for each motor unit of a given muscle, a quadrant displaying the cumulative spike trains (CST) of two groups of motor units for a given muscle, a quadrant displaying the firing rate of two motor units, and the smoothed discharge rates of all the identified motor units for a given muscle (Figure 2).

Figure 2.

Electromyographic Recordings

The accuracy and boundary capabilities of the online EMG decomposition algorithm was tested on a series of experimental data collected with either a high-density grid or an intramuscular array of electrodes from five different muscles.

Surface EMG signals were recorded from the gastrocnemius medialis and lateralis (GM and GL; triceps surae protocol) or the vastus medialis and lateralis (VM and VL; quadriceps protocol) with two-dimensional adhesive grids of 64 electrodes [13 × 5 electrodes with one electrode absent on a corner, gold-coated, interelectrode distance: 8 mm; (GR08MM1305, OT Bioelettronica, Italy)]. Surface EMG signals were recorded from the tibialis anterior with a two-dimensional adhesive grid of 64 electrodes [13 × 5 electrodes with one electrode absent on a corner, gold-coated, interelectrode distance: 4 mm; (GR04MM1305, OT Bioelettronica, Italy)]. Before the placement of the grids, the skin was shaved and cleaned with an abrasive gel (Nuprep, Weaver and company, USA). Each adhesive grid was held on the skin using a semidisposable biadhesive foam layer. The cavities within the adhesive layer were filled with conductive paste (SpesMedica, Italy) to facilitate the skin-electrode contact. A 10-cm wide elastic band was placed over the electrodes to ensure good contact between the electrodes and the skin throughout the experiment.

An intramuscular linear array of 16 electrodes on a thin-film (platinum coated, interelectrode distance: 1mm) was inserted into the tibialis anterior in one participant at an approximate angle of 30°. The insertion was guided with a portable ultrasound probe (Butterfly IQ+, Butterfly Network, USA).

A reference electrode was positioned over the tibia of the right limb (triceps surae protocol), over the patella of the right limb (quadriceps protocol), or over the medial malleolus (tibialis anterior protocol). A strap electrode dampened with water was placed around the ankle (ground electrode) for each data collection. The EMG signals were recorded in monopolar mode and digitized together with the torque signal at a sampling rate of 2048 Hz for the grids of surface electrodes and 10240 Hz for the intramuscular array of electrodes (EMG-Quattrocento, 400 channel EMG amplifier; OT Bioelettronica, Italy).

Experimental procedure

Eight physically active male volunteers (mean ± standard deviation; age: 28 ± 5 yr, height: 181 ± 5 cm, body mass: 77 ± 18 kg) participated in the triceps surae protocol, and eight physically active male volunteers (age: 31 ± 6 yr, height; 181 ± 5 cm; body mass: 78 ± 16 kg) participated in the quadriceps protocol; with data recorded using a pair of surface grid electrode. Two individuals participated in both the triceps surae and quadriceps protocols. Five physically active male volunteers (age: 26 ± 4 yr; height: 174 ± 7 cm; body mass: 66 ± 15 kg) participated in the tibialis anterior protocol, with data recorded using a surface grid electrode; and one physically active volunteer (age: 30 yr, height: 176 cm, body mass: 70 kg) participated in the tibialis anterior protocol with an intramuscular array of electrodes.

We specifically recruited males as: i) it is well known within the field that a greater number of motor units can be decomposed from the signals of males compared to females (although the reasons for this discrepancy are yet to be well established (28); and ii) we required a large number of motor units to test the accuracy of real-time identification. None of the participants reported lower limb injury or pain in the 6 months prior to testing. Ethical committees approved the study (triceps surae and quadriceps protocols: CERNI – Nantes Université, n°04022022; tibialis anterior protocol: Imperial College London, no. 18IC4685). All participants provided their informed written consent before the beginning of the experiment.

The right side of the body was tested for all participants and for all protocols. For the triceps surae protocol, participants sat on a dynamometer (Biodex System 3 Pro, Biodex Medical, USA) with their hip flexed at 80°, 0° being the neutral position, and their right leg fully extended. Their ankle angle was set to 10° of plantarflexion, 0° being the foot perpendicular to the shank. For the quadriceps protocol, participants sat on the dynamometer with their hips flexed at 80° and the knee of their right leg flexed at 80°, 0° being the full extension. Inextensible straps were tightened during both tasks to immobilize the torso, pelvis, and thigh on the test side. For the tibialis anterior protocol, participants sat on a chair while their foot was fixed onto the pedal of a dynamometer (OT Bioelettronica, Italy) coupled with a load cell (CCT Transducer, Italy) and positioned at 30° in the plantarflexion direction (0° being the foot perpendicular to the shank). The foot was fixed to the pedal with inextensible straps positioned around the proximal phalanx, metatarsal and cuneiform. Force signals were recorded using the same acquisition system as for the HD-EMG recordings (EMG-Quattrocento; OT Bioelettronica, Italy).

All experiments began with a standardized warm up, which included five 3-s isometric plantar flexion or knee extension contractions at 50%, 60%, 70%, and 80%, and three 3-s contractions at 90% of the participants’ subjective maximal torque. Then, after 2 min of rest, participants performed three maximal voluntary contractions (MVC) for 3–4 s, with 60-s of rest in between. Peak MVC torque was considered as the maximal value obtained from a moving average window of 250 ms.

For the triceps surae and quadriceps protocols, participants performed three trapezoid isometric contractions at 40% of the MVC (referred to as baseline contractions) to identify motor unit filters offline. Each of these contractions involved a 5-s ramp-up, a 20-s plateau, and a 5-s ramp-down phase and was separated by 60 s of rest. The separation vectors (motor unit filters) were identified offline from these contractions and then applied in real-time to estimate the firing activity of each identified motor unit during three additional trapezoid contractions (referred to as the online task). Each of these online tasks involved a 5-s ramp-up, a 30-s plateau and a 5-s ramp-down phase and was separated by 60 s of rest. To test the effect of variations in contraction intensity between the online task and the baseline contraction used to identify the separation vectors (motor unit filters), each plateau consisted of three successive 10-s targets at 35%, 40%, or 45% of the MVC performed in a random order. During these online tasks, feedback of motor unit discharge times and torque output was displayed in real-time on a monitor to the participants.

To test the effect of the baseline contraction intensity on the accuracy of real-time identification of motor unit discharge activity; the procedure was repeated with a baseline intensity of 20% MVC. During the last three trapezoidal contractions, each plateau consisted of three successive 10-s targets at 15%, 20%, or 25% of the MVC performed in a random order.

For the tibialis anterior protocol with grids of surface electrodes, participants performed a trapezoid contraction at 20% of the MVC, involving a 10-s ramp up, a 60-s plateau, and a 10-s ramp down phase (baseline contraction). The separation vectors (motor unit filters) were identified from this contraction and then applied in real-time over a second identical contraction (online task). The same procedure was repeated for the tibialis anterior protocol with an intramuscular array of electrodes, with contractions involving a ramp up phase of 2-s, a plateau of 20-s, and a ramp down phase of 2-s.

Data processing

Identification of separation vectors (motor unit filters) with offline EMG decomposition

The monopolar EMG signals collected during the baseline contractions were extended with an extension factor of , where m is the number of channels free of any noise or artifact. The signals were then demeaned and whitened. A fixed-point algorithm was applied to identify the sources of the EMG signals, i.e., the motor unit spike trains (Figure 1B). Motor unit spike trains can be considered as sparse sources with most samples being 0 (i.e., absence of spikes) and a few samples being 1 (i.e., spikes). In this algorithm, a contrast function was iteratively applied to estimate a separation vector (the motor unit filter) that maximized the level of sparseness of the identified source. Convergence was reached when the level of sparseness did not vary when compared to the previous iteration, with a tolerance fixed at 10^-4. At this stage, the estimated source contained high peaks (i.e., the spikes from the identified motor unit) and low peaks from other motor units and noise. High peaks were separated from low peaks and noise using peak detection and K-mean classification with two classes: ’spikes’ and ’noise’ (Figure 1B). The peaks from the class with the highest centroid were considered as spikes of the identified motor unit. A second algorithm refined the estimation of the discharge times by iteratively recalculating the separation vector and repeating the steps with peak detection and K-mean classification until the coefficient of variation of the inter-spike intervals was minimized. The accuracy of each estimated spike train was assessed by computing the silhouette (SIL) value between the two classes of spikes identified with K-mean classification (18). When the SIL exceeded a predetermined threshold (see below), the motor unit filter was saved for the real-time decomposition, together with the centroids of the ‘spikes’ and ‘noise’ classes (Figure 2A).

Even though a SIL value of 0.9 is generally used for offline decomposition (18), we purposely decreased this threshold to 0.8 for the triceps surae and quadriceps experiments. We made this choice such that we maximised the number of motor units for testing the online decomposition capabilities. It is important to note that the spike trains from all of these retained units were visually checked for false positives and false negatives. The spike trains with no clear separation between the spikes and the noise were removed by the investigator. For the remaining motor units, manual editing was then performed by the investigator. Manual editing requires: i) identifying and removing the spikes of low quality (i.e., those with peaks close in amplitude to the level of noise), ii) re-calculating the motor-unit filter and re-applying it over a portion of the signal, and iii) adding the new spikes recognized as motor unit firings.

For the tibialis anterior experiments, the SIL value was set at 0.9. Visual checking and manual edition of the spike trains were performed using the same process as for the triceps surae and quadriceps experiments.

Real-time identification of motor neuron activities

EMG signals were transmitted by packages of 256 data points for the surface grid recordings (125 ms with a sampling frequency of 2048 Hz) or 1280 data points for the intramuscular array recordings (125 ms with a sampling frequency of 10240 Hz). The electrode mask determined from the baseline contractions was applied to remove the channels with noise or artifacts and the data was extended using the same extension factor as for the baseline contraction (see the section on EMG decomposition above). The matrix of separation vectors (motor unit filters) identified during the baseline contraction was applied over the extended EMG signal. Local peaks were identified for each motor unit using the MATLAB function “islocalmax” with a minimal separation of 25 ms between peaks to limit the number of false positives (Figure 2B). These peaks were considered as spikes when their distance from the centroid of the ‘spike’ class was shorter than the distance from the centroid of the ‘noise’ class (Figure 2B). Both centroids were identified from the offline decomposition made on the baseline contraction.

We calculated the CST of each muscle as the sum of the spikes of all the identified motor units over a moving window of eight consecutive epochs of 125 ms (total duration = 1s; units: spikes per second). Similarly, we calculated the discharge rate (i.e., spikes per second) of individual motor units as the sum of the spikes over a moving window of eight consecutive epochs of 125 ms). We chose this approach to facilitate the smoothness of the visual feedback, in contrast to the instantaneous discharge rate, which is calculated as the instantaneous inversed interspike interval, which oscillates more due to the presence of synaptic noise. While this approach introduces a delay of 500 ms in the estimation of the discharge rate, we identified in pilot testing that this also facilitates the control of motor unit discharge activities by the participant. To further increase the smoothness of the online biofeedback, we added the option in the software to average CSTs or individual discharge rates using a moving window. For the data reported in this paper, we selected a value of four consecutive values, corresponding to a window of 500 ms. Note that researchers who aim to minimize the delay of the visual feedback can disable this option, change its value, or use the real-time raster plots that displays the instantaneous discharge times of each motor unit.

Computational time

The computational time depends on the number of identified motor units during the baseline contraction, the number of peaks sorted during each epoch, and the number of EMG channels retained for the analysis. We considered the computational time for the decomposition as the time between the reception of the EMG signals by the computer and the estimation of the discharge times of all the identified motor units. We considered the computational time for the feedback as the time between the decomposition of the EMG signals and the update of the online feedback. The computational times were calculated on a laptop equipped with an Apple M1 Max chip and 64 GB of RAM.

Accuracy of the real-time identification of motor unit firing activity

To assess the accuracy of the real-time identification of motor unit spike trains, we compared the motor unit spike trains identified in real-time with those obtained after manual edition. The manual edition was performed offline as described above (28, 29).

The accuracy of the real-time decomposition was assessed for each motor unit by computing the sensitivity, the precision, the false negative rate, and the rate of agreement between the manually edited spike train (offline) and the spike train identified in real-time, as follows:

Where TP (true positive) is the number of spikes identified in both the real-time and edited spike trains, FP (false positive) is the number of spikes only identified in the real-time spike train and FN (false negative) is the number of spikes only identified in the edited spike train.

To assess the accuracy of the biofeedback provided by the software, we measured the root mean square error (RMSE) between the path drawn by the smoothed discharge rate of motor units estimated in real-time and the path estimated from the manually edited motor unit spike trains.

Results

Offline EMG decomposition of the baseline contraction and manual editing

After the completion of the isometric baseline contractions, the recorded EMG signals were decomposed over 150 iterations, which took approximately 5 min per grid of electrodes. When considering the baseline contractions performed at 20% MVC, on average 20 ± 9 (VL), 14 ± 5 (VM), 25 ± 11 (GL), and 19 ± 9 (GM) motor units per participant were identified with a SIL>0.8. Their SIL values calculated before any manual editing were 0.89 ± 0.05 for VL, 0.83

± 0.04 for VM, 0.82 ± 0.03 for GL, and 0.87 ± 0.04 for GM. After manual editing, a significant number of motor units was removed as their spike train showed no clear separation between the spikes and the noise. The number of motor units removed was: 6 ± 5 for VL, 10 ± 5 for VM, 21 ± 12 for GL, and 6 ± 5 for GM (Figure 3). The remaining motor units exhibited a SIL value of 0.91 ± 0.04, 0.89 ± 0.06, 0.89 ± 0.03, and 0.90 ± 0.02 for VL, VM, GL, and GM, respectively (Figure 3).

Figure 3.

When considering the baseline contraction performed at 40% of MVC, on average 17 ± 6 (VL), 15 ± 5 (VM), 29 ± 13 (GL), and 13 ± 4 (GM) motor units were identified by automatic decomposition, with 12 ± 8 (VL), 3 ± 3 (VM), 5 ± 6 (GL), and 4 ± 3 (GM) motor units remaining after manual editing (Figure 3). The SIL value of the selected motor units reached 0.91 ± 0.03 for VL, 0.89 ± 0.05 for VM, 0.87 ± 0.02 for GL, and 0.88 ± 0.03 for GM (Figure 3).

The above results highlight that the motor units kept for the online EMG exhibit a SIL value well above 0.8. For the TA experiments, we chose to set a threshold at 0.9 for the minimal SIL for the offline decomposition to significantly decrease the burden of visual checking and manual editing. On average, we identified 15 ± 4 motor units with a SIL value of 0.94 ± 0.01 during the baseline contraction performed with grids of surface electrodes placed over the TA muscle. We identified 10 motor units with an average SIL value of 0.95 ± 0.02 during the baseline contraction performed with the intramuscular array of electrodes inserted into the TA (Figure 3A). Only one motor unit was removed from each dataset after manual editing. It is noteworthy that manual editing did not change the SIL value for any of the motor unit, either recorded with the grid of surface electrodes or the intramuscular array (Figure 3B).

Computational time of online EMG decomposition and visual feedback

On average, the time between the reception of an epoch of EMG signals and the identification of the spikes was 4.9 ± 3.0 ms for data collected with high-density grids of surface electrodes.

To estimate the computational time per motor unit, we performed a linear regression between the number of identified motor units and the computational time per epoch using data from all the experiments made with high-density surface electrodes. The slope of the linear fit, i.e., the computational time per motor unit, was 0.37 ms (Figure 4A). Across all the experiments, the maximal computational time was 13.1 ms for 26 identified motor units. Of note, when considering the experiment made with intramuscular arrays of electrodes, the computational time reached 32.0 ms for 9 identified motor units. This was longer than observed for the same number of identified motor units using high-density surface electrodes (5.2 ms; Figure 4B) and was due to the greater sampling frequency used for the intramuscular recordings.

Figure 4.

We then calculated the additional computational time to display visual feedback in the form of a raster plot or smoothed discharge rates of all the identified motor units. On average, we observed computational times of 8.6 ± 5.7 ms and 11.2 ± 7.6 ms for raster plots and smoothed discharge rates, respectively. The computational time per motor unit was 0.74 ms and 0.99 ms for raster plots and smoothed discharge rates, respectively (Figure 4C and 4D). Of note, the computational time for the quadrant plot made from the activity of two motor units reached 14.8 ± 0.0 ms on average. It is noteworthy that the computational times for experiments with grids of electrodes or intramuscular arrays of electrodes was similar regardless of the type of visual feedback (Figure 4B). Overall, as the total computational time was constantly shorter than the duration of an epoch of EMG signals, the visual feedback was always updated during the recording of the next epoch of EMG signals. Therefore, the only delay was the incompressible recording time per epoch of signals, i.e., 125 ms.

Accuracy of the real-time identification of motor unit firing activity

We assessed the accuracy of the motor unit spike trains identified in real time using their manually edited version as reference (Figure 5). It is noteworthy that some motor units present in the baseline contractions were not identified during the online contractions, either because they were not recruited, or because they were merged with newly recruited motor unit(s), i.e., their separation vector was no longer unique to other vectors. At 20% of MVC, the number of missed units was between 0 (TA intramuscular array) and 4 ± 4 (GM), and at 40% of MVC between 1 ± 1 (GM) and 4 ± 5 (GL).

Figure 5.

For the sake of clarity, we only report in this section the rates of agreement for each intensity and muscle. Values of sensitivity, precision, and rates of false negatives are reported in Figure 5. The highest rate of agreement was observed for the TA (0.94 ± 0.09 with the grid and 0.97 ± 0.05 with the intramuscular array). Those values were lower for the vastii (VL: 0.82 ± 0.20; VM: 0.75 ± 0.18) and gastrocnemii muscles (GL: 0.88 ± 0.08; GM: 0.81 ± 0.18). When considering the contractions performed at 40% of MVC, the rate of agreement was 0.84 ± 0.16 for VL, 0.75 ± 0.21 for VM, 0.82 ± 0.08 for GL, and 0.88 ± 0.09 for GM (Figure 5).

For the triceps surae and quadriceps experiments, we also asked the participants to match force targets set at -5% (i.e., 15% and 35% of MVC for 20% and 40% of MVC conditions) and at +5% (i.e., 25% and 45% of MVC for 20% and 40% of MVC conditions). These contractions aimed to test the robustness of the separation vectors (motor unit filters) to changes in contraction intensity, i.e., in the number of active motor units generating the EMG signal. Results in terms of sensitivity, precision, false negative rate and rate of agreements are displayed on Figure 6.

Figure 6.

Accuracy of the real-time biofeedback

Because the accuracy of the raster plot feedback is directly related to the accuracy of the estimation of the discharge times reported in the previous section, here we only focus on the accuracy of the smoothed discharge rates (Figure 2). It can either be displayed in a quadrant, with a cursor moving within a two-dimensional space according to the discharge rates of two motor units, or with one cursor for each identified motor units moving in the vertical direction and following a scrolling path.

At 20% of MVC, the root mean squared error (RMSE) between the real time smoothed discharge rate and its edited version was consistently below two pulses per second for all the motor units (Figure 7). The lowest RMSE was observed with the TA muscle (0.64 ± 0.77 pps and 0.44 ± 0.43 with grids and intramuscular arrays of electrodes). RMSE was 1.4 ± 1.5 for VL, 1.8 ± 1.5 for VM, 1.7 ± 1.1 for GL, and 1.1 ± 1.1 pps for GM. At 40% of MVC, the RMSE reached 1.1 ± 1.2 for VL, 1.8 ± 1.1 for VM, 0.8 ± 0.5 for GL, and 0.6 ± 0.4 pps for GM. Overall, these RMSE provides strong evidence that the biofeedback accurately reflects the motor unit discharge activity.

Figure 7.

Discussion

We have introduced an open-source software based on blind source separation for decoding the activity of alpha motor neurons in real-time. We first presented the full framework, from the initialisation of the real time decoding using an offline decomposition of EMG signals with fastICA, to the real-time feedback. We demonstrated that: i) the computational time to decompose the EMG signals in real-time is low, i.e., 0.37 ms per motor unit when high-density EMG electrodes are used, ii) the highest accuracy of the real time decoding is observed with intramuscular array comparing to grids of surface electrodes, iii) the highest accuracy of the real time decoding is observed for the tibialis anterior relatively to the other muscles, iv) the smoothed motor unit discharge rates displayed to the participant is accurately estimated from real time EMG decomposition, with a RMSE value lower than 2 pulses per second. Overall, this open-source software provides a set of tools for neuroscientists to design experimental paradigms where participants can receive neurofeedback in real-time on the output of the spinal cord circuits.

There are processing assumptions for the blind-source separation algorithm to accurately identify motor unit discharge times from multichannel EMG signals. Among them, the most important is the uniqueness of the distribution of motor unit action potentials across electrodes (that defines the separation vector) among all the other active motor units within the recording volume (10, 25). One way to satisfy this condition is to increase the selectivity of the electrodes to record the discharge activity of only a few motor units from a small volume of muscle (4). For example, Figure 1A shows a motor unit action potential detected only over three to four electrode sites along the array of intramuscular electrodes, while a motor unit action potential can be observed across many more electrodes with grids of surface electrodes. Therefore, the likelihood of having spatially overlapping motor unit action potentials is lower, which explains why the rate of agreement of motor units identified from intramuscular arrays of electrodes is much higher than with single channel wire electrodes or grids of surface electrodes (7, 8). A second way to increase the percentage of discriminable motor units among all the active motor units in the recording volume is to increase the spatial sampling of their activity using multiple electrodes (3, 10, 17, 25). This has led to the growth of EMG recording systems with high-density grids of surface electrodes (3), which compensate for the lack of specificity of motor unit action potential profiles that are recorded when using a pair of traditional bipolar EMG electrodes (30).

Another necessary condition for EMG decomposition using the blind-source separation algorithm is the consistency of the motor unit filters across time. An obvious reason inducing changes in motor unit filters would be the displacement of the electrodes relatively to the source. Such drifts also exist with intracortical arrays and can be corrected with appropriate methods that track waveforms across electrodes (13). However, the geometry of the muscle tissue is much more variable than that of the cortical tissue, especially during movements. For example, muscle bellies become bulkier while shortening (31), increasing the distance between the surface electrodes and the deep sources. In addition, the pennation angle of muscle fibres can change with contraction intensity (32), modifying the direction of the propagation of motor unit action potentials along the fibres relatively to the position of the electrodes. All these factors impact the recorded motor unit action potential profiles across electrodes, which in turn will reduce the capacity to discriminate the same motor unit from the EMG signal (27, 33, 34). For these reasons, we recommend applying our approach during isometric contractions with a stable level of force, as even large changes in force during isometric contractions can impact the orientation of the muscle fibres relatively to the skin within muscles (35).

The first step of the algorithm consists of identifying motor unit separation vectors (motor unit filters) from EMG signals with fastICA (18, 20, 36). Classically, metrics such as the Pulse to Noise ratio (37) or the silhouette value (18) are used to assess the reliability of the identified motor units by estimating the distance between the spikes and the noise. Here, we purposely lowered the threshold classically recommended for offline decomposition (i.e., 0.9 in (18)) to maximize the number of identified units (Figure 3A). It is noteworthy that such approach must be associated with extensive manual editing to remove the discharge times incorrectly selected and to add missing spikes (28, 29). This necessary step precedes the update of the motor unit filter with all the detected spikes (38), leading to an increase in the silhouette value and a decrease in the coefficient of variation of the interspike intervals (Figure 3B). Alternatively, one could speed up the processing flow by setting a higher threshold for the silhouette value. This would decrease the burden of manual editing at the cost of decreasing the number of identified motor units. However, it is important to note that using more stringent criteria does not preclude the manual editing, even for motor units with a high silhouette value (see motor units 1 and 2 with missed spikes in Figure 2).

The second step of the algorithm consists of identifying motor units discharge times in real time by projecting extended segments of EMG signals into separation vectors (motor unit filters) to estimate motor unit pulse trains (Figure 2B; (20, 22)). This method was effective at identifying motor unit discharge times, with rates of agreement >0.75, regardless of the contraction intensity and muscle (Figure 5). It is noteworthy that the performance was particularly high for the recordings made with an intramuscular array of electrodes (rate of agreement of 0.97 ± 0.05, Figure 5). As mentioned above, the better performance of blind source separation on multichannel intramuscular EMG has already been reported with offline analyses (7, 8, 19). This is explained by the higher spatial selectivity of the electrodes (4), the larger bandwidth of the signal (30), and the higher robustness of the motor unit filter as the signal is less affected by the geometric changes of the volume conductor. In contrast, precision and rates of agreements were lower for motor units identified over the vastii and gastrocnemii muscles when compared to the TA muscle (Figure 5). Even though the reason for this between-muscle difference is unclear, it is possible that the specific activation level of the vastii and gastrocnemii muscles varied more than that of TA between the baseline and test contractions, because of muscle redundancy leading to multiple coordination strategies possible to perform knee flexion or plantarflexion. For example, a decrease in activation would decrease the height of the peaks of the estimated sources, potentially classified in the noise class during the online decomposition (22). Conversely, an increase in activation may activate motor units spatially close to those observed during the baseline contraction, corrupting their pulse train with merged sources (16, 17). This explanation is in line with the lower rate of agreement observed when participants tracked a force target 5% of the MVC higher than the level of the baseline contraction (Figure 6). One way to overcome these challenges would be to dynamically updating the motor unit filters and the centroids of the ‘spikes’ and ‘noise’ classes (22). While appealing, this approach is also computationally demanding, with computational times significantly higher than those observed in this study (22). We propose in our software to update the motor unit filters and the centroids of the ‘spikes’ and ‘noise’ classes during the resting periods. In addition, it is recommended that the operator displays the motor unit pulse trains and identified discharge times between contractions to check for decomposition accuracy across the session.

Overall, the main purpose of our software is to display to the participant a real-time visual feedback based on the activity of individual motor units or populations of motor units. The root-mean square error of the smoothed discharge rates was constantly below two pulses per second, with values as low as 0.44 ± 0.43 pps for the TA muscle recorded with intramuscular electrode arrays. Thus, the movement of the cursors accurately represented the variations in motor unit discharge activity to the participant. In this study we have presented results of control of smoothed CST over a relatively large smoothing window, but the duration of the smoothing filter when computing the CST can be chosen by the user according to the needs and applications. Operators could use the provided software to interact with a virtual environment, such as typing on a keyboard with a cursor moved by modulating motor unit discharge rates (24). In the field of motor control, neuroscientists may train participants to selectively activate individual motor units (21), testing the concept of rigid versus flexible motor control (21, 24, 39), or movement augmentation (40). Generally, we hope that this open-source software will open perspectives for neuroscientists to design experimental paradigms that takes advantage of online EMG decomposition to study the neural control of movements at the motor neuron level.

Supplemental data availability

The entire data set (raw and processed data) codes, and a user manual of the software are available at https://github.com/simonavrillon/I-Spin

Disclosures

DF is inventor of a patent (Neural Interface. UK Patent application no. GB1813762.0. August 23, 2018) and of a patent application (Neural interface. UK Patent application no. GB2014671.8. September 17, 2020), which are partly related to the methods and applications of this work. All other authors have no financial or other relationships that might lead to a conflict of interest.