Wavelet-Based Spike Sorting of Muscle Spindle Afferent Nerve Activity Recorded With Thin-Film Intrafascicular Electrodes

– The continuous complex wavelet transform offers a convenient framework for neural spike sorting. Results show that wavelet-based neural spike detection outperforms simple threshold detection, especially with signals with low signal to noise ratio. Classification of action potentials using their signatures in wavelet space performed as well as a classifier based upon principal components analysis, and better than a classifier based upon template matching. Applied on experimental intrafascicular recordings of muscle spindle afferent nerve response to passive muscle stretch, the spike sorting algorithm manages to isolate afferent activity of units having a linear relationship between neural firing rate and muscle length, an important step towards a model-based estimator of muscle length.


Introduction
Closed-loop control of neural prosthetics is seen as a means to compensate for many of the non-linearities related to muscle activation, and the unexpected perturbations one might encounter when using a Functional Electrical Stimulation system.With the advent of advanced implanted neuroprosthetic interfaces, natural sensors, such as the muscle spindle afferent nerve, are being explored as an alternative source for feedback information.In a previous study, a simple feedback controller was implemented to follow a desired joint angle trajectory in the presence of externally applied disturbances [1].Muscle length was estimated using aggregate unspecific muscle afferent activity and a lookup table.In another study a simple linear model was proposed to express the relationship between length and unit neural firing rate [2], [3].In order to improve the estimation of muscle length, separation of units in the multi-unit ENG is necessary.
Action potentials are usually detected by looking for signal activity rising above a preset noise gate [4], [5].The detected spikes are classified using different methods.Principal components analysis (PCA) and template matching are often used [4], [5] because of their relative simplicity which enables fast, real-time implementation.In a recent study it was demonstrated that using a continuous complex wavelet based approach to action potential detection in the auditory nerve outperforms a matched filter approach [6].In this study we expand this detection algorithm so it covers a range of temporal scales of action potentials.Additionally, we introduce a classifier of the detected action potentials.The spike sorting algorithm is evaluated on both synthesized data and experimentally recorded muscle spindle afferent response to imposed muscle stretches.The algorithm performance is compared to methods using PCA and template matching.

Acute animal experiments
Acute experiments were conducted on 10 anesthetized New Zealand white rabbits.One leg of the rabbit was anchored at knee and ankle joints to a fixed mechanical frame The Achilles tendon was attached to the arm of a servocontrolled motor using a ligature.Pulling on the ligature stretched the Medial Gastrocnemius (MG) muscle and resulted in ankle dorsiflexion.Releasing the tension on the ligature resulted in the muscle shortening to its original length and movement of the ankle in the plantar direction.
A thin-film longitudinal intrafascicular electrode (tfLIFE) was implanted in the MG nerve in the rabbit's hind limb.The electrode enabled simultaneous multichannel monitoring of the ENG from the fascicle in which the structure was implanted.In order to have purely muscle afferent activity in the ENG recordings, the sciatic nerve was crushed proximal to the implantation site.Signals from the structure were amplified and filtered before recording them to a multichannel digital tape recorder.
The MG muscle was presented with sinusoidal stretches 4 mm in peak-to-peak amplitude and 250 mHz in frequency.The force, relative position and the 4 neural signals from the tfLIFE were simultaneously recorded.

Generating synthetic signals
In order to be able to evaluate the performance of the spike sorting algorithm, knowledge of the exact timing and class of each action potential in the ENG signal is needed.
Because of uncertainty about this information in experimental data, artificial signals based upon recorded action potentials were synthesized.Around 20 action potentials with different waveforms were visually identified and extracted from the recorded experimental data.From the set of extracted waveforms, 5 action potentials with distinctly different shapes were chosen to represent 5 neural spikes originating from different axons, i.e. different spike classes.The waveforms were normalized and used to synthesize spike trains.Spike train firing rates were randomly chosen from within the range found in literature [7].With the exception of burst firing, muscle spindle afferents can fire with a rate up to about 75 Hz.Spike amplitudes were scaled by integer values ranging from 3 to 6 standard deviations of the background noise level, which corresponds to values found by visually inspecting the data recorded using tfLIFE.These 4 scaling factors represented different SNR levels for which the analysis was performed.Signals were synthesized by adding the spike trains onto experimentally recorded background noise.Signals with different number of units firing simultaneously were synthesized, starting from signals having only 2 units firing, up to signals with 10 classes.Spikes having the same waveform and amplitude were considered to be from the same axon (belonging to the same class).A total of 900 signals were generated: 100 signals with 2 units active, another 100 with 3 units, and so on until having 100 signals with 10 units firing simultaneously.

Neural spike sorting 2.3.1 Spike detection
Neural spike detection was based on computing the continuous wavelet transform (CWT) of the ENG using the complex Gaussian wavelet.This wavelet was chosen because its real and imaginary parts resemble the typical action potential shapes found in experimental data.As the output wavelet coefficients are a measure of the degree of similarity of patterns in the signal to the wavelet waveform, the detection algorithm consisted of searching for wavelet coefficients above a preset threshold.Detection was evaluated on a range of thresholds.
Inspection of extracted action potentials waveforms showed that the durations of neural spikes is in a range between a few hundred microseconds up to 1 ms.A multiscale CWT analysis was performed on extracted action potential waveforms in order to find the scales at which the CWT coefficients were maximal.Only the range of scales for which the magnitude of the coefficients were larger than 95% of the magnitude of the maximal coefficient where used in the detection algorithm.

Spike classification
The multi-scale complex CWT also offers a framework for classifying the detected neural spikes.Action potentials differ in their shape and amplitude and it is necessary to choose a feature set and a distance metric with which they would be distinguishable.The multi-scale CWT offers a redundant representation of the spike's shape in time-scale space.Exploratory data analysis on the extracted neural spike waveforms indicated that the CWT coefficients computed using the same set or scale factors as in the detection phase could be suitable as classification features.Classification was performed by creating feature vectors for the detected spikes and then clustering the data using the Euclidean distance metric.

Evaluation on synthetic signals
Performances of the wavelet-based detector and the detector using simple amplitude thresholding are shown on Figure 1 in the form of receiver operating characteristics (ROC curves).These curves are graphical representations of detector sensitivity vs specificity using a range of detection thresholds.On the whole range of SNR levels the wavelet-based detector outperforms the detector based on amplitude thresholding, i.e. for any given specificity, the corresponding sensitivity is greater for the wavelet-based detector.The performance gap becomes especially prominent with low SNR.
Classification results are shown in the form of classification error rates which are ratios of the number of misclassified spikes to the total number of spikes being classified.The wavelet-based classification is compared to two other methods of classification: principal components analysis and template matching.Results are shown starting from the case when only two different spike classes are present in the signal up until the case where 10 units are simultaneously firing (Figure 2).Classification based on template matching produced the highest classification error rates, while wavelet-and PCA-based approaches showed similar results.

Evaluation on experimental data
The spike sorting technique was applied on experimentally recorded muscle spindle afferent nerve activity.Only periods of plantar flexion were analyzed, since the units are silent during the dorsiflexion period.The detection threshold was chosen to be seven times the standard deviation of the background noise level (in wavelet space).Throughout all the trials this threshold value corresponded to the point on the ROC curves where the specificity starts to rapidly deteriorate and there is little improvement in sensitivity.The detected units were classified into 10 clusters.This approximately matches the number of units being picked up by one recording site of the tfLIFE when the muscle is maximally stretched (preliminary results of recent unpublished work).The analysis was performed on data from 3 rabbits.Preliminary results show that from the 10 classes, up to 2 or 3 spike classes show a linear relationship between their computed neural firing rate and instantaneous muscle length.Results from one rabbit are shown on Figure 3.The left plot shows the relationship between the aggregate firing rate of all detected spikes.The relationship is clearly not linear in the region where the normalized muscle length is close to 1.The right plot shows the same, but this time the firing rate was computed using only the activity of two classes having a linear relationship between unit firing rate and muscle length.A linear regression analysis performed on both shows that the fit on the right plot is obviously better.

Discussion
Results show that the continuous wavelet transform using complex wavelets is the preferred method for neural spike detection.Even though wavelet-based classification does not show improvement in error rates compared to the PCAbased algorithm, the advantage of using the wavelet-based approach is in the fact it is it provides a unique framework for both spike detection and classification, i.e. after computing the complex wavelet coefficients in the detection phase, no additional computation is required in the subsequent classification phase.
Preliminary results show that the spike sorting scheme is capable of isolating the activity of single units having a linear relationship between firing rate and muscle length.The aggregate firing rate from muscle spindle afferents should have some non-linearities where most of the units are being recruited or saturate.There should be a range of length where they behave linearly.From the point where a single unit reaches activation threshold, its activity should have a more linear firing rate to length behavior.This is one possible explanation of the more linear behavior of the firing rate computed using a subset of classes.
On-line estimation of muscle length would require a real-time implementation of the wavelet transform, which has been shown to be feasible with well chosen wavelet families [8].
With each new spike detected, either a new cluster would be formed, either the detected spike would be associated to an existing cluster.In the latter case the cluster centroid would be updated which would also make the algorithm capable of tracking slow changes in spike waveform shape and amplitude (either due to electrode drift, either to electrode fibrous encapsulation leading to a degradation of SNR in the recorded signals).
One of the biggest obstacles in both detecting and classifying neural spikes proves to be the low SNR of the ENG recordings.In parallel to working on the improvement of signal processing methods to minimize the effects of noise, novel recording techniques could also be used to improve signal quality [9].It is expected that, when signals are acquired using such techniques, the overall performance of the spike sorting algorithm should improve.
Results from experimental data presented in this paper are preliminary.Final conclusions can not be drawn before processing the data from all 10 rabbit experiments, which is work currently in progress.Comparison with different muscle stretch profiles is also necessary.

Conclusions
The continuous complex wavelet transform offers a convenient framework for both detection and classification of action potentials using intrafascicular electrodes.The neural spike detection outperforms the simple threshold detection, especially with signals with low SNR.Results of classification of units into classes based on action potential waveform signatures in wavelet space indicate that the classifier manages to isolate afferent activity of units having a linear relationship between neural firing rate and muscle length, which is an important step towards a modelbased estimator of muscle length.
Fig. 3 Afferent neural firing rate vs muscle length.The left plot shows the aggregate activity of all detected spikes.On the right, only activity from 2 clusters, having a good linear fit to the data, are used to compute the firing rate.Linear regression analysis was performed on both (full lines).On both plots muscle length was normalized by 4 mm.

Fig 1 .Fig 2 .
Fig 1. ROC curves for four different SNR levels.Performances of a simple threshold detector (empty-circle line) and the waveletbased detector (full-triangle line) are compared