Skip to main content

Estimating the neural spike train from an unfused tetanic signal of low-threshold motor units using convolutive blind source separation

Abstract

Background

Individual motor units have been imaged using ultrafast ultrasound based on separating ultrasound images into motor unit twitches (unfused tetanus) evoked by the motoneuronal spike train. Currently, the spike train is estimated from the unfused tetanic signal using a Haar wavelet method (HWM). Although this ultrasound technique has great potential to provide comprehensive access to the neural drive to muscles for a large population of motor units simultaneously, the method has a limited identification rate of the active motor units. The estimation of spikes partly explains the limitation. Since the HWM may be sensitive to noise and unfused tetanic signals often are noisy, we must consider alternative methods with at least similar performance and robust against noise, among other factors.

Results

This study aimed to estimate spike trains from simulated and experimental unfused tetani using a convolutive blind source separation (CBSS) algorithm and compare it against HWM. We evaluated the parameters of CBSS using simulations and compared the performance of CBSS against the HWM using simulated and experimental unfused tetanic signals from voluntary contractions of humans and evoked contraction of rats. We found that CBSS had a higher performance than HWM with respect to the simulated firings than HWM (97.5 ± 2.7 vs 96.9 ± 3.3, p < 0.001). In addition, we found that the estimated spike trains from CBSS and HWM highly agreed with the experimental spike trains (98.0% and 96.4%).

Conclusions

This result implies that CBSS can be used to estimate the spike train of an unfused tetanic signal and can be used directly within the current ultrasound-based motor unit identification pipeline. Extending this approach to decomposing ultrasound images into spike trains directly is promising. However, it remains to be investigated in future studies where spatial information is inevitable as a discriminating factor.

Introduction

The smallest functional element that can be voluntarily activated in the neuromuscular systems comprises a motoneuron and a group of muscle fibres innervated by the motoneuron. This functional element is called a motor unit (MU). The gold standard for measuring the characteristics of a large population of MUs is based on high-density surface electromyography (sEMG) [1]. High-density sEMG consists of a grid of electrodes that records mixed MU activities superficially from the skin (< 20 mm). Then, the activity is decomposed into single MU spike trains by blind source separation (BSS) techniques, which use many electrodes [2, 3]. This approach provides access to the neural drive of the spinal cord via the motoneurons to specific muscles [4]. Although sEMG has been used successfully for MU analysis, it is well-known that it has limited spatial selectivity and field of view.

Ultrafast ultrasound has been shown to image and analyse voluntarily activated MUs for a large field of view in the muscle (40 × 40 mm) [5,6,7,8,9,10,11], providing spatiotemporal mechanics at a high resolution (< 1 mm and > 1 kHz). This technique is based on recording radio frequency signals (B-mode images) (Fig. 1A, B), calculating displacement velocity images (Fig. 1C) and separating the images into spatiotemporal components, i.e. each component is associated with a (1) spatial map and (2) a time signal (Fig. 1D). A subset of the components is putative estimates of the (1) MU territory and (2) a sequence of MU twitches (unfused tetanic signal) evoked by the spike trains [5, 6, 9, 10]. The spike trains are estimated from the unfused tetanic signal (Fig. 1E). The spike train estimation is currently based on a Haar wavelet method (HWM), which detects the twitch onsets due to the large gradient changes, and it has good performance [12]. Although this ultrasound-based pipeline (Fig. 1) has great potential to provide comprehensive access to the neural drive of the spinal cord to muscles for a large population of MUs simultaneously, the method currently has a limited identification rate of the active MUs [6].

Fig. 1
figure 1

The methodological pipeline to identify voluntarily activated motor units (MUs) using ultrafast ultrasound. A The ultrasound probe is placed on the skin to record data from the muscle’s cross-section (transverse view). B The recorded B-mode images. C Calculating displacement velocity images. D Separating the velocity images into spatiotemporal components using instantaneous blind source separation (BSS) with a focus on spatial sparsity, where each component is associated with a time signal and a spatial image. A subset of the components is putative estimates of the (1) MU territory in the spatial maps and (2) a sequence of MU twitches (unfused tetanic signal) evoked by the spike trains. E The spike trains are estimated based on unfused tetani. Previously using a Haar wavelet method (HWM). In this study, we estimate the spikes trains in E using the unfused tetanic signals from D using convolutive blind source separation (CBSS) and compare the performance against another method, i.e. the Haar wavelet method (HWM)

The estimation of spikes partly explains the current pipeline’s limited identification rate of active MUs from each estimated unfused tetanic signal that often is noisy [6, 13]. Since HWM is based on large gradients, the noise in the unfused tetanic signals may induce false positives. Therefore, we must consider alternative methods with at least similar performance and robust against noise, among other factors. An alternative method with great potential is based on convolutive blind source separation (CBSS), which has been used with great success on multichannel sEMG signals [2, 3, 14]. Since CBSS uses multiple samples over time to estimate each spike, it should be able to estimate spike trains from single-channel unfused tetanic signals robustly with high performance.

The main challenge of using CBSS to estimate the spike train from an unfused tetanic signal is that the unfused tetanus of a MU will consist of variable successive twitches [15,16,17], where the twitch duration is longer than the time between two succeeding spikes (i.e. an ISI). In contrast, the MU action potential is more consistent, and its duration is shorter than the ISI. Previous studies of spike train estimation of an unfused tetanic signal, i.e. the HWM, have exploited the similarity of each twitch rise (onset gradient) and observed that the twitch rise duration is shorter than the ISI [12]. Assuming the twitch rise is the action potential counterpart, and the remaining twitch activity is another additive noise component, we hypothesise that estimating the spike train using CBSS should have a high agreement with the ground truth. If this is the case, this method can be used directly in the current ultrasound-based pipeline to assess the neural drive to the muscle (Fig. 1). Indeed, it also implies the feasibility of extending the method to spatiotemporal data such that spikes can be estimated directly from the ultrasound-based image sequences.

This study aimed to estimate spike trains from simulated and experimental unfused tetani using a CBSS algorithm and compare it against the previously optimised and evaluated HWM [12]. We implemented the CBSS algorithm requiring three parameters based on fixed-point iterations [2, 3, 14] and peak detection [12]. We evaluated the parameters of CBSS using the estimated spike trains’ rate of agreement with the ground truth for the simulations. We also explored the algorithm’s parameters’ relation with the spike delta (time difference between truth or reference spikes and the estimated spikes) and its variation. Finally, we compared CBSS against the HWM based on their rate of agreement performance using simulated and experimental unfused tetanic signals from voluntary contractions of humans and evoked contraction of rats.

Methods

The overview of the methods includes (1) simulations describing the model and the parameters, (2) the experimental data, (3) the CBSS algorithm to estimate the spike train from an unfused tetanic signal, (4) the parameter evaluation of CBSS using simulated data, and finally, (4) the comparison between CBSS and HWM using simulated and experimental data.

Simulations

Simulation model

We used a simulation model that generated a motor unit (MU) unfused tetanus \(Y\left(t\right)\) at time \(t\ge 0\):

$$Y\left(t\right)=\sum_{i=1}^{n}{F}_{i}\left(t-{\Delta }_{i}\right)+\varepsilon \left(t\right),$$
(1)

where \(t=\{\mathrm{0,1},\dots ,T\}\) is the time, and \(n\) is the number of times that the MU fires. \({F}_{i}\) is the \(i\)th twitch generated by the MU at its \(i\)th spike at time instant \({\Delta }_{i}\). The duration of \({F}_{i}\) is longer than the inter-spike interval (ISI) \({\Delta }_{i+1}-{\Delta }_{i}\). \(\varepsilon (t)\) is additive noise. Note that the twitch shape commonly varies within a contraction, which has been explained in other studies for voluntary contractions [15,16,17]. We express the corresponding spike train \(s\left(t\right)\ge 0\) as:

$$s\left(t\right)=\sum_{i=1}^{n}\delta \left(t-{\Delta }_{i}\right),$$
(2)

where \(\delta \) is the Dirac Delta function.

Simulation parameters

Similar to a previous study [12], each simulated unfused tetanus \(Y\left(t\right)\) was based on simulating \(n\) twitches and \(n\) spikes (Fig. 2). The twitches \({F}_{i},i\in \{\mathrm{1,2},\dots ,n\}\) were simulated using a twitch model [17] with five parameters (Fig. 2A) that were randomly sampled to generate varying twitch-shapes (Fig. 2B) suitable for low-threshold MUs at low force levels [13] (see Additional file 1: Table S1). The \(n\) spike times \(\Delta =({\Delta }_{1},{\Delta }_{2},\dots ,{\Delta }_{n})\) were simulated using a Gaussian renewal process such that \({\Delta }_{i+1}-{\Delta }_{i} \sim N({\mu }_{ISI},\mathrm{CV})\), where \({\mu }_{ISI}\) and \(\mathrm{CV}=\frac{{\sigma }_{ISI}}{{\mu }_{ISI}}\) is the mean ISI and coefficient of variation [18] (Fig. 2C). We simulated combinations of different average firing rates (8, 12, and 16 Hz) [6, 19] and ISI CV values (5, 20, and 40%) [20] feasible for low-threshold MUs.

Fig. 2
figure 2

The different steps of the convolutive blind source separation algorithm. A The input signal is unfused tetanus. B The unfused tetanus is extended using the extension factor \(R\), where the extended signal has a dense covariance matrix. C Performing whitening such that the covariance matrix equals the identity matrix. D The fixed-point algorithms find the separation (projection) vector (red line), and by projecting it on the whitened extended signal in C, the resulting output is the source estimate (blue line). E The time instants of the local maxima of the source estimate were detected using peak detection with two parameters: the minimal peak distance (\(\mathrm{MPD}\)) and the number of mean absolute deviations (\(n\mathrm{MAD}\)). F The time instants of the local maxima are the estimated firing times of the spike train \({\varvec{s}}\) (red line)

Each twitch was summed to the respective spike resulting in an unfused tetanic signal (Fig. 2C). The unfused tetanus was differentiated with respect to time (from force to yank or displacement to velocity) [21] because of the stability of its mean value and velocity is used in MU identification using ultrasound. The last step was adding normally distributed noise so that the signal-to-noise ratio (SNR) was at three levels, i.e. 10, 20, and 30 dB. For more information on the simulations, see Additional file 1.

Experimental signals

We retrospectively included, from previous studies, two datasets, including n = 21 unfused tetani [12]. The first dataset included three 2-s-long unfused tetani from the biceps brachii of three healthy human subjects (28.3 ± 0.6 years; one male and two females) at very low isometric force levels using the ultrasound-based MU identification pipeline with a (concentric) needle EMG electrode as [6] (see Additional file 1). The ultrasound and EMG systems were synchronised (sampling rates 2 kHz and 64 kHz).

The second dataset included 18 unfused tetani from the medial gastrocnemius of five adult female Wistar rats [22]. After a surgical procedure resulting in functionally isolated MUs, the Achilles tendon was connected to a force transducer while stretched to achieve isometric conditions where the force was measured simultaneously as a wire electrode (sampling rates 1 kHz and 10 kHz). Individual axons were stimulated based on average ISI times: 60, 70, 80, and 100 ms. For two rats, the last three stimuli intervals were used. The intervals between the individual stimuli were randomly set (mean ISI ± 50%). Before further analysis, the force-based signals were (1) filtered using a sixth-order zero-phase Butterworth bandpass filter with a high- and low-pass cut-off equal to 3 and 100 Hz and (2) differentiated. For detailed information about the functional isolation of MUs, see Additional file 1.

The convolutive blind source separation algorithm to estimate the spike train

Previously a shift-invariant model has been used to describe a multichannel sEMG signal assuming (1) constant action potential waveforms for the same MU in one channel and (2) that the waveform duration is shorter than every ISI [2, 3, 14]. This study considers a single-channel unfused tetanus \(Y\left(t\right)\) \(,\) which has the following two characteristics: (1) the twitch waveforms for the same MU vary [15,16,17], and (2) the twitch waveform duration is longer than the ISI [16, 23]. Given this, we express the time derivative of the unfused tetanic signal as:

$$\widetilde{Y}\left(t\right)=\sum_{i=1}^{n}\sum_{l=1}^{L}{h}_{i}\left(l\right)s\left(t-l\right)+\widetilde{\varepsilon }\left(t\right),$$
(3)

where \(i=\{1,\dots ,n\}\) denotes the \(i\mathrm{th}\) spike, \({h}_{i}\left(l\right)\) is a twitch rise waveform of length \(L\) as a response to the \(i\)th firing, \(s\left(t\right)\) is the spike train at time \(t\), and \(\widetilde{\varepsilon }(t)\) is additive noise including the remaining signal of the twitch waveforms after subtracting the twitch rise waveform. The convolutive mixture with finite impulse response filters in Eq. (3) can be represented as a linear and an instantaneous mixture by extending the observation signal, source signal, and noise using their \(R\) delayed versions [2, 3, 14]:

$$\widetilde{{\varvec{Y}}}\left(t\right)={\left[\widetilde{Y}\left(t\right),\widetilde{Y}\left(t-1\right),\dots ,\widetilde{Y}\left(t-R+1\right)\right]}^{\mathrm{T}},$$
$$\widetilde{{\varvec{s}}}\left(t\right)={\left[s\left(t\right),s\left(t-1\right),\dots ,s\left(t-L-R+2\right)\right]}^{\mathrm{T}},$$
$$\widetilde{{\varvec{\varepsilon}}}\left(t\right)={\left[\widetilde{\varepsilon }\left(t\right),\widetilde{\varepsilon }\left(t-1\right),\dots ,\widetilde{\varepsilon }\left(t-R+1\right)\right]}^{\mathrm{T}},$$
(4)

where \(L\) denotes the twitch rise waveform length and \(R\) denotes the extension factor. Then, the linear instantaneous mixture model is defined as:

$$\widetilde{{\varvec{Y}}}=\widetilde{{\varvec{A}}}\widetilde{{\varvec{s}}}+\widetilde{{\varvec{\varepsilon}},}$$
(5)

where \(\widetilde{{\varvec{A}}}\) contains the twitch rise waveforms \({h}_{i}\). To estimate the spike train \({\varvec{s}}=s\left(t\right),t=\{\mathrm{0,1},\dots ,T\}\) using Eq. (5), we used a linear instantaneous BSS model referred to as independent component analysis (ICA) [24].

Algorithm

Given an unfused tetanic signal (Fig. 3A), the spike train \({\varvec{s}}\) was estimated by solving the separation problem in Eq. (5) using a fixed-point algorithm [25], which has been used for decomposing multichannel EMG signals [2, 3]. After extending the signal (Fig. 3B), the next step, common to most BSS algorithms [24], was to spatially whiten the extended observation matrix using eigenvalue decomposition (Fig. 3C). The whitened extended observation matrix had a covariance matrix equal to the identity (Fig. 3C). After whitening, the fixed-point algorithm estimates the separation (projection) vector and, thereby, the source estimate (Fig. 3D). Here, the separation vector was initialised using a normally distributed random vector. The fixed-point algorithm is based on maximising the non-Gaussian distribution of the source vector using a cost function \(G\left(\bullet \right)\) [25]. We used a contrast function associated with sparseness since the spike train is sparse (mostly zeros with a few ones). Here, we choose \(G\left(x\right)=\mathrm{log}\left(\mathrm{cosh}\left(x\right)\right)\) due to its kurtotic distribution in line with the characteristics of a spike train that is kurtotic (and skewed) in contrast to a normal distribution.

Fig. 3
figure 3

Simulating an unfused tetanic signal using a twitch model and a spike train. A A twitch (blue line) is generated by its five parameters (black dashed lines) of the model are the half-contraction time Thc (ms), the contraction time Tc (ms), the half-relaxation time Thr (ms), the twitch duration Ttw (ms), and the maximal amplitude Amax (arbitrary unit, a.u.). B Generated twitches using the twitch model in A by randomly sampling their parameters from a uniform distribution. C Time instants (spikes, black vertical lines) for each twitch were simulated using a Gaussian renewal process. The unfused tetanus was obtained by summating the twitches in B at their respective time instants with added normally distributed noise to achieve a certain signal-to-noise (SNR) level. Note that the unfused tetanic signal has been differentiated with respect to time, providing unfused tetanus concerning the rate of change (e.g., yank or velocity instead of force or displacement)

After obtaining the source estimate, the spike train was estimated by peak detection to identify the time instants of the local maxima. The peak detection was based on (1) the height of the peaks and (2) the distance between consecutive peaks in milliseconds (Fig. 3E). The height of the peaks was based on the number of mean absolute distances (\(nMAD\)) of the source estimate. The distance between the consecutive peaks was based on the minimum peak distance (\(MPD\)) where the highest peak within that distance was selected. The resulting time instants of the local maxima were considered the estimated firing times of the spike train \({\varvec{s}}\) (Fig. 3F). The overall CBSS algorithm is summarised in pseudo-code below:

figure a

Parameter evaluation for CBSS

We evaluated the algorithms’ three parameters jointly (based on simulations): the extension factor \(R\) and peak detection parameters \(nMAD\) and \(MPD\). After initial tests to locate the global maxima, we evaluated the following parameters: \(R\in \left[\mathrm{10,15},\dots ,40\right]\), \(nMAD\in \left[\mathrm{1.0,1.5},\dots ,4.0\right]\), and \(MPD\in \left[\mathrm{10,15},\dots ,40\right]\) ms. For each parameter combination, 2700 unfused tetani were generated based on the combinations of SNR values (30, 20, and 10 dB), ISI CV values (5, 20, and 40%), and firing rates (8, 12, and 16 Hz) as explained in a previous study [12].

The rate of agreement (RoA) between the simulated (or reference) and estimated spikes was based on the following metric: \(\mathrm{RoA}=\frac{100\times \mathrm{CS}}{\mathrm{CS}+\mathrm{FS}+\mathrm{MS}}\), where \(CS\) was the number of correctly identified spikes, \(\mathrm{FS}\) was the number of false spikes, and \(\mathrm{MS}\) was the number of missed spikes. The tolerance window for correctly identified spikes was set between 0 and \(R+10\) (extension factor delay) ms. This tolerance was selected based on maximal RoA (see Additional file 1: Fig. S1). We also analysed the time between a correctly estimated spike (\(\mathrm{CS}\)) and the simulated spike (referred to as spike delta) and its variability (spike delta variability) for different extension factors \(R\). The difference in spike delta variation between different extension factors (in parameter evaluation) was tested using Levene’s test with a significance level of 5% [26].

Comparison between CBSS and HWM—simulated and experimental data

The parameters that maximised the mean RoA of CBSS were extracted based on averaging the RoA of all simulated parameter combinations. For HWM, we used the optimised parameters from a previous study with a similar parameter evaluation procedure using the same RoA threshold (between − 5 and 25) [12]. Note that HWM consists of three steps: (1) computing the continuous Haar wavelet transform of the signal at a given pseudo-frequency, (2) standardising the wavelet transform coefficients based on z-score, and (3) finding the time instants of the local maxima of the wavelet coefficients using peak detection [12].

We reported the mean and standard deviation of the RoA values for each method, and we tested the pairwise differences between the spike trains of CBSS and HWM in three ways. First, we compared them using the same simulation procedure in the parameter evaluation for various SNR, ISI CV, and firing rate values and tested the pairwise differences using Wilcoxon signed rank test. Second, we compared them using the experimental unfused tetani. Third, we calculated the computational time to run each algorithm on the simulated data (2700 times) and tested the pairwise differences between the algorithms. The significance level was set to 5%, and the p-values were corrected using Bonferroni correction.

Results

Parameter evaluation for CBSS—simulations

The highest mean RoA for CBSS was achieved when \(R=20\) with 97.5 ± 1.6% (Fig. 4A). Given extension factor \(R=20\), we found that the peak parameters associated with the highest RoA values and smallest standard errors were \(n\mathrm{MAD}=2\) and \(\mathrm{MPD}=20\) (Fig. 4B, C).

Fig. 4
figure 4

Parameter evaluation of the convolutive blind source separation (CBSS) algorithm using many simulations. A The mean rate of agreement (RoA) and its 95% standard error for different extension factors \(R\). Highest mean RoA was for \(R=20\). B Given the extension factor \(R=20\), the peak detection parameters (\(n\mathrm{MAD}\) and \(\mathrm{MPD}\)) that had the highest mean RoA were \(n\mathrm{MAD}=2\) and \(\mathrm{MPD}=20\). C The standard errors associated with the mean RoA values in B

We found that the spike delta differed for different extension factors (Fig. 5A), but also for the same extension factor at different iterations (Fig. 5B). The larger the extension factor, the greater variability in the mean spike delta (Fig. 5C). However, there was no difference (p > 0.05) in spike delta variability for the different extension factors (Fig. 5A, B, D). The spike delta variation was 2.2 ± 0.5 ms.

Fig. 5
figure 5

Spike delta and spike delta variability for different extension factors (\(R\)) regarding CBSS. A The spike delta differed for different \(R\), but there was no difference in spike delta variability. B The spike delta differed for the same \(R\) for five reiterations of the algorithm, but there was no difference in spike delta variability. C The spike delta differed for different \(R\), and there was a difference in spike delta variability. D Subtracting the spike delta for each simulated signal: the spike delta did not differ for different \(R\), and there was no difference in spike delta variability

Comparison between CBSS and HWM—simulated data

Considering all parameter combinations, CBSS had a higher RoA than HWM (97.5 ± 2.7 vs 96.9 ± 3.3, p < 0.001). By considering each simulation parameter separately, CBSS had a higher RoA than HWM in 12 of the 27 simulation parameter combinations (Table 1). In contrast, HWM had a higher RoA than CBSS in 6 parameter combinations. For the other nine combinations, there was no difference in RoA. Finally, CBSS was more robust to high ISI CV values than HWM (Table 1). However, there was no clear pattern with respect to noise levels or firing rates.

Table 1 Comparing the performance of CBSS and HWM to estimate spikes based on simulated signals with varying parameters

The computational time was much lower for HWM than for CBSS (p < 0.001). However, there is no practical difference since the computational time of HWM is approximately 2 ms and CBSS 11 ms for a signal that is about 10 s (Fig. 6).

Fig. 6
figure 6

The computational time for the convolutive blind source separation (CBSS) and the Haar wavelet method (HWM). Each boxplot is based on 2700 computational time values, i.e. 100 signals containing 100 spikes for 27 different simulation parameters. ***p < 0.001

Comparison between CBSS and HWM—experimental data

The three unfused tetani from human biceps brachii had an average firing rate of 11.4 Hz and an average ISI CV of 13.4%. The estimated spike trains from CBSS and HWM highly agreed with the EMG reference spikes (98.0 ± 3.4% and 96.3 ± 3.2%, respectively). See Table 2.

Table 2 Comparing the performance of CBSS and HWM to estimate spikes based on experimental signals (voluntary and evoked contractions)

The 18 unfused tetani from rat gastrocnemius had an average firing rate of 14.4 Hz and an average ISI CV of 30.6%. The estimated spike trains from CBSS and HWM highly agreed with the EMG reference spikes (98.0 ± 2.8% and 96.5 ± 5.6%, respectively). Given this, there was no difference between the two methods (p = 0.24). See Table 2.

Discussion

This work aimed to estimate simulated and experimental unfused tetani spike trains using CBSS and compare it against HWM [12]. For this purpose, we evaluated the parameters of CBSS using simulations and explored the algorithm’s parameters’ relation with the spike delta and its variation. Then, we compared CBSS against the HWM based on their RoA using simulated and experimental unfused tetanic signals from voluntary contractions of humans and evoked contraction of rats. The main finding was that CBSS had (on average) a higher performance than HWM with respect to the simulated firings than HWM (97.5 ± 2.7 vs 96.9 ± 3.3, p < 0.001). In addition, we found that the estimated spike trains from CBSS and HWM highly agreed with the experimental spike trains (98.0% and 96.4%).

We have shown that CBSS can be used to estimate the spike train of an unfused tetanic signal since the firings highly agree with the simulated (see Fig. 4A and Table 1) and EMG reference spike trains from two different experimental datasets (see Table 2). We found that the spike delta variability did not depend on the extension factor, although the spike delta differed for different extension factors and the same extension factor (Fig. 5A, B). This observation suggests that there are multiple local maxima, and the convergence to different local maxima may depend on the initialisation of the separation vector [2]. Also, a larger extension factor leads to a larger spike delta variance (between trials, Fig. 5C), suggesting that the convergence issue increase with the extension factor value. This problem could be overcome using a more standardised initiation of the separation factor or finding a standardised way to shift the estimated spikes without knowing the ground truth. Although there is potential for improvement, both CBSS and HWM used optimal RoA thresholds of the same length (0 to 30 ms and − 5 to 25 ms, respectively), where CBSS had, on average, better performance than HWM (Table 1).

The computational time was about five times slower for CBSS than for HWM. This finding is explained by CBSS extending the signal from a vector to a matrix to make the convolutive approach an instantaneous linear problem suitable for independent component analysis (ICA) [2, 3]. Then, fixed-point iterations are used to find a projection vector. However, this computational time difference is of no practical difference since the computational time of HWM is approximately 2 ms and CBSS 11 ms for a signal that is about 10 s (Fig. 6). For example, a time difference of 200 ms between visual feedback and movement leads to the perception of a delay [27]. The computational times are well within that time difference without optimisation.

This study considered observations of unfused tetani. Although CBSS can be used directly to estimate firings based on the estimated unfused tetanus from the ultrasound-based pipeline [6, 12], this study indicates the potential of either extending or including the temporal CBSS approach to the current spatially focused BSS method [6, 28] to improve the separation of displacement velocity images from ultrasound to increase the identification rate [8]. Another solution would bypass the estimation of an unfused tetanic signal (Fig. 1D) and go directly to spikes (from Fig. 1C to Fig. 1E), thereby reducing the pipeline’s exposure to error propagation, i.e. estimating spikes based on occasionally poor estimates of an unfused tetanic signal. Given this, the next challenge emerges when expanding to ultrasound images, i.e. successive twitches within the same MU may differ. There is a possibility that twitches from a MU may be highly similar to the ones of other MUs. However, one could overcome this challenge by including spatial information, which has a high resolution (< 1 mm) as the MU relates to a physical component (muscle unit) in the spatial domain. A potential solution may be expanding current sEMG decomposition algorithms [2, 3] to include spatial dependence or sparsity in addition to the temporal deconvolution and validating it using an authentic simulation model [29]. Yet, all these approaches need to deal with the motion of non-MU-related structures that hides a large part of the movement caused by a MU in ultrasound images [11]. One potential solution could be to use the spatiotemporal clutter filtering approach to improve the sensitivity to detect microvascular networks or blood flows corrupted by significant tissue or probe motion artefacts [30]. Nevertheless, the implementation and validation of these approaches will be investigated in the future.

In conclusion, this study estimated simulated and experimental unfused tetani spike trains using a CBSS algorithm and compared its performance against a previously optimised method, i.e. HWM. We found that the estimated spike trains from CBSS and HWM highly agreed with the simulated and EMG reference spike trains, and CBSS had, on average, higher performance. This result implies that the CBSS of an unfused tetanic signal can be used to estimate its spike train, and it can be used directly within the current ultrasound-based MU identification pipeline. Extending this approach to decomposing ultrasound images into spike trains directly is promising. However, it remains to be investigated in future studies where spatial information is inevitable as a discriminating factor.

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article (and its additional file). The code and data are available on request from the corresponding author RR. For requests about the experimental dataset 2 (on rats), the corresponding author RR will pass the request to Professor Jan Celichowski.

References

  1. Farina D, Merletti R, Enoka RM. The extraction of neural strategies from the surface EMG: an update. J Appl Physiol. 2014;117:1215–30.

    Article  Google Scholar 

  2. Chen M, Zhou P. A Novel Framework Based on FastICA for High Density Surface EMG Decomposition. IEEE Trans Neural Syst Rehabil Eng. 2016;24:117–27.

    Article  Google Scholar 

  3. Negro F, Muceli S, Castronovo AM, Holobar A, Farina D. Multi-channel intramuscular and surface EMG decomposition by convolutive blind source separation. J Neural Eng. 2016;13:26027.

    Article  Google Scholar 

  4. Farina D, Holobar A, Merletti R, Enoka RM. Decoding the neural drive to muscles from the surface electromyogram. Clin Neurophysiol. 2010;121:1616–23.

    Article  Google Scholar 

  5. Rohlén R, Stålberg E, Stöverud KH, Yu J, Grönlund C. A method for identification of mechanical response of motor units in skeletal muscle voluntary contractions using ultrafast ultrasound imaging - simulations and experimental tests. IEEE Access. 2020;8:50299–311.

    Article  Google Scholar 

  6. Rohlén R, Stålberg E, Grönlund C. Identification of single motor units in skeletal muscle under low force isometric voluntary contractions using ultrafast ultrasound. Sci Rep. 2020;10:1–11.

    Article  Google Scholar 

  7. Ali H, Umander J, Rohlén R, Grönlund C. A deep learning pipeline for identification of motor units in musculoskeletal ultrasound. IEEE Access. 2020;8:170595–608.

    Article  Google Scholar 

  8. Rohlén R, Yu J, Grönlund C. Comparison of decomposition algorithms for identification of single motor units in ultrafast ultrasound image sequences of low force voluntary skeletal muscle contractions. BMC Res Notes. 2022;15:207.

    Article  Google Scholar 

  9. Carbonaro M, Meiburger KM, Seoni S, Hodson-Tole EF, Vieira T, Botter A. Physical and electrophysiological motor unit characteristics are revealed with simultaneous high-density electromyography and ultrafast ultrasound imaging. Sci Rep. 2022;12:1–14.

    Article  Google Scholar 

  10. Carbonaro M, Zaccardi S, Seoni S, Meiburger KM, Botter A, Detecting anatomical characteristics of single motor units by combining high density electromyography and ultrafast ultrasound: a simulation study. 44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). IEEE. 2022;2022:748–51.

    Google Scholar 

  11. Lubel E, Grandi-Sgambato B, Barsakcioglu DY, Ibanez J, Tang M-X, Farina D. Kinematics of individual muscle units in natural contractions measured in vivo using ultrafast ultrasound. J Neural Eng. 2022;34:78.

    Google Scholar 

  12. Rohlén R, Antfolk C, Grönlund C. Optimization and comparison of two methods for spike train estimation in an unfused tetanic contraction of low threshold motor units. J Electromyogr Kinesiol. 2022;67:78.

    Article  Google Scholar 

  13. Rohlén R, Raikova R, Stålberg E, Grönlund C. Estimation of contractile parameters of successive twitches in unfused tetanic contractions of single motor units – A proof-of-concept study using ultrafast ultrasound imaging in vivo. J Electromyogr Kinesiol. 2022;45:102705.

    Article  Google Scholar 

  14. Holobar A, Zazula D. Multichannel blind source separation using convolution kernel compensation. IEEE Trans Signal Process. 2007;55:4487–96.

    Article  MathSciNet  MATH  Google Scholar 

  15. Burke RE, Rudomin P, Zajac Iii FE. The effect of activation history on tension production by individual muscle units. Brain research Elsevier. 1976;109:515–29.

    Article  Google Scholar 

  16. Raikova R, Celichowski J, Pogrzebna M, Aladjov H, Krutki P. Modeling of summation of individual twitches into unfused tetanus for various types of rat motor units. J Electromyogr Kinesiol. 2007;17:121–30.

    Article  Google Scholar 

  17. Raikova R, Pogrzebna M, Drzymała H, Celichowski J, Aladjov H. Variability of successive contractions subtracted from unfused tetanus of fast and slow motor units. J Electromyogr Kinesiol. 2008;18:741–51.

    Article  Google Scholar 

  18. Fuglevand AJ, Winter DA, Patla AE. Models of recruitment and rate coding organization in motor-unit pools. J Neurophysiol. 1993;70:2470–88.

    Article  Google Scholar 

  19. Milner-Brown HS, Stein RB, Yemm R. Changes in firing rate of human motor units during linearly changing voluntary contractions. J Physiol. 1973;230:371.

    Article  Google Scholar 

  20. Tracy BL, Maluf KS, Stephenson JL, Hunter SK, Enoka RM. Variability of motor unit discharge and force fluctuations across a range of muscle forces in older adults. Muscle Nerve. 2005;32:533–40.

    Article  Google Scholar 

  21. Lin DC, McGowan CP, Blum KP, Ting LH. Yank: the time derivative of force is an important biomechanical variable in sensorimotor systems. J Exp Biol. 2019;222:e180414.

    Article  Google Scholar 

  22. Drzymała-Celichowska H, Celichowski J. Functional isolation of single motor units of rat medial gastrocnemius muscle. J Visual Exp. 2020;34:e61614.

    Google Scholar 

  23. Krutki P, Pogrzebna M, Drzymała H, Raikova R, Celichowski J. Force generated by fast motor units of the rat medial gastrocnemius muscle during stimulation with pulses at variable intervals. J Physiol Pharmacol. 2008;59:85–100.

    Google Scholar 

  24. Hyvärinen A, Karhunen J, Oja E. Independent Component Analysis. New York: Wiley; 2001.

    Book  Google Scholar 

  25. Hyvärinen A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans Neural Netw. 1999;10:626–34.

    Article  Google Scholar 

  26. Levene H. Robust tests for equality of variances. Contributions to probability and statistics Essays in honor of Harold Hotelling. New York: Stanford University Press; 1960. p. 279–92.

    Google Scholar 

  27. Zbinden J, Lendaro E, Ortiz-Catalan M. A multi-dimensional framework for prosthetic embodiment: a perspective for translational research. J Neuroeng Rehabil. 2022;19:122.

    Article  Google Scholar 

  28. Rohlén R, Lundsberg J, Malesevic N, Antfolk C. A fast blind source separation algorithm for decomposing ultrafast ultrasound images into spatiotemporal muscle unit kinematics. Bioengineering. 2022. https://doi.org/10.1101/2022.11.22.517488.

    Article  Google Scholar 

  29. Ali H, Umander J, Rohlén R, Röhrle O, Grönlund C. Modelling intra-muscular contraction dynamics using in silico to in vivo domain translation. Biomed Eng Online. 2022;21:1–19.

    Article  Google Scholar 

  30. Demené C, Deffieux T, Pernot M, Osmanski B-F, Biran V, Gennisson J-L, et al. Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity. IEEE Trans Med Imaging. 2015;34:2271–85.

    Article  Google Scholar 

Download references

Acknowledgements

We thank Professor Jan Celichowski and Rositsa Raikova for providing experimental data on rats’ irregularly stimulated single motor units.

Funding

Open access funding provided by Lund University. This work was supported by the Swedish Research Council for Sport Science (D2023-0003), the Promobilia Foundation, Stiftelsen för bistånd åt rörelsehindrade i Skåne and the Swedish Research Council (DNR 2019–05601).

Author information

Authors and Affiliations

Authors

Contributions

RR designed and implemented the study; RR, JL, and CA analysed and interpreted the data and implementation; RR wrote the first draft; RR, JL, and CA revised the paper. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Robin Rohlén or Christian Antfolk.

Ethics declarations

Ethics approval and consent to participate

For experimental dataset 1: The subjects gave informed consent verbally and in written form before the experimental procedure. The project conformed to the Declaration of Helsinki and was approved by the Swedish Ethical Review Authority (Reference 2019–01843).

For experimental dataset 2: The experimental procedures followed the EU animal care guidelines, and the principles of Polish Law on the Protection of Animals were approved by the Local Bioethics Committee. This dataset was obtained (and consent for using it was given) by Professor Jan Celichowski, Department of Neurobiology, Poznan University of Physical Education, Poland.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no financial and non-financial competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1. Table S

1. The simulation parameters. Figure S1. Sensitivity analysis of the tolerance interval for rate of agreement (RoA) calculation using different values of extension factor \(R\) (10, 20, and 30). We considered a simulation of 5,000 spikes and subsequent twitches using a firing rate of 12 Hz, an ISI CV of 20%, and an SNR of 20 dB. The red crosses denote the selected tolerance intervals.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rohlén, R., Lundsberg, J. & Antfolk, C. Estimating the neural spike train from an unfused tetanic signal of low-threshold motor units using convolutive blind source separation. BioMed Eng OnLine 22, 10 (2023). https://doi.org/10.1186/s12938-023-01076-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12938-023-01076-0

Keywords