Selective peripheral nerve recordings from nerve cuff electrodes using convolutional neural networks

The same in vivo data as in our previous study [32] was used here. Briefly, nine Long-Evans rats (retired breeders) were placed under isoflurane anesthesia. Access to the sciatic nerve was obtained through an oblique incision on the posterior and dorsal aspect of the hip. This provided access through the natural splitting of the fibers of the gluteus maximus after clearing of the skin and deep fascia. A 56-multi-contact spiral polyimide nerve cuff electrode (CorTec GmbH, Freiburg, Germany) was placed on the nerve to record neural activity, just proximal to its branching into the tibial, peroneal and sural nerves. The electrode consisted of seven rings of eight contacts distributed over the length of the electrode. The cuff had a length of 23 mm and a diameter of 1 mm. A needle electrode, placed in the back of the animal, served as a reference.

Data was acquired through a neural data acquisition board (RHD2000, Intan Technologies, USA) using a sampling rate of 30 kHz. Recordings were bandpass filtered on the acquisition board between 256 Hz and 7.5 kHz. Afferent activity consisting of dorsiflexion, plantarflexion and pricking of the heel was selectively evoked in different fascicles of the sciatic nerve using mechanical stimuli. The foot was held by the claws using a clip with plastic tips and dorsiflexion/plantarflexion of the ankle (approximately 60°) was applied manually to evoke proprioceptive activity in the tibial and peroneal branches, respectively [40]. A cutaneous stimulus to the heel, applied using a Von Frey monofilament (300 g), was used to elicit activity in the sural branch [41]. 100 trials of each activity (dorsiflexion, plantarflexion and pricking of the heel), along with an alternating dorsi- and plantar- flexion activity were collected. The application of stimuli was video-recorded using a webcam at 30 frames s⁻¹ and synchronized to the neural activity through a visible light emitting diode (LED) in the videos. The videos were later manually annotated to determine the times at which each type of stimulus was being applied, and this process provided the ground truth nCAP labels for training the classifier.

The experimental procedures were approved by the Animal Care Committee of the University of Toronto and all experiments were performed in accordance with the Animal Care Committee's guidelines.

This framework first used in [32] involves extracting the spatiotemporal signatures of nCAPs associated with neural pathways of interest (figure 1(a)). The defining characteristic of the ESCAPE framework is the use of these spatiotemporal signatures to improve recording selectivity. In this context, we use the term 'neural pathway' to refer to a small group of fibers with a related function and firing approximately synchronously. The resulting nCAP produced by such a group of fibers is conjectured to underlie a 'spike' (i.e. sharp change in amplitude) visible in the nerve cuff recording (as discussed further below).

Zoom In Zoom Out Reset image size

Preprocessing: First, the raw signal is referenced using a referencing montage (e.g. tripolar [42] or common average [43]) and any filtering needed is applied per contact. nCAP Detection: The preprocessed signals are used to obtain an average signal for each ring and delay-and-add operations, as described in the velocity selective recording methods [33–35, 37, 44], are applied obtaining a single average signal of the middle ring with improved SNR. nCAPs are then detected by applying a thresholding method (e.g. median absolute approach [45]) to this average signal. Signature Extraction: nCAPs' spatiotemporal signatures are extracted from the original filtered referenced signal to be used in the classifier training and evaluation. Classification: The extracted spatiotemporal signatures are used as input to a classifier that associates each nCAP with a neural pathway. Specific steps used in this study are explained below, a visual example of the signal throughout the framework is shown in figure 1(b), and spatiotemporal signatures of each class are shown in figure 2.

Zoom In Zoom Out Reset image size

2.2.1. Preprocessing.

2.2.2. CAP detection

Next the average signal in each ring was obtained and these signals were shifted in time against each other by an integer multiple of a time-step dt and summed to form a single signal (delay-and-add). A dt of 1 time sample was used, matching the expected nCAP conduction velocity, such that signals in the different rings sum and increase the amplitude of the delay-and-add output. This cleaned signal was then thresholded using the median absolute deviation estimate approach (equation (1)) [45] to find nCAP locations. Detected nCAPs with peak values above 15 µV in the tripolar referenced recordings (before delay-and-add) were discarded to minimize the inclusion of noise or artefacts in the dataset (previous analysis found a negligible impact of small variations in this threshold value on the end results [32]).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm Threshold}=~4\times \frac{{\rm Median}(\left| X \right|)}{0.6745}.\nonumber \end{align} \tag{ 1 }$$

2.2.3. Spatiotemporal signatures

Once nCAPs were detected, 49 time samples before and 50 time samples after the peak location were used to create the spatiotemporal signatures. These spatiotemporal signatures were extracted from the signal after tripole referencing, but before the delay-and-add operation that was used solely for nCAP detection.

Spatiotemporal signatures for each activity (dorsiflexion, plantarflexion and pricking of the heel) were constructed using detected nCAPs for each stimulus using M contacts at T consecutive time samples. This created an M  ×  T matrix constituting the spatiotemporal signature for that nCAP. By associating each nCAP in a new recording with one of these known signatures, the activities of multiple neural pathways can be discriminated.

Two different spatiotemporal signatures were created by numbering contacts within a ring (radially), creating a representation that had a 'spatial emphasis' (SE), or numbering contacts along rings (longitudinally), creating a representation that had a 'temporal emphasis' (TE) as illustrated in figure 3. In addition, different referencing montages (tripolar, and common average) were applied to create alternative representations of the spatiotemporal signatures. The tripolar reference montage uses the average signal from the two outer rings of contacts (i.e. average signal of the 16 contacts) of the nerve cuff electrode as a reference for all contacts. The common average reference montage takes the average signal across all contacts as the reference. These spatiotemporal signatures were then used as inputs into the CNN.

Zoom In Zoom Out Reset image size

CNNs are comprised of three fundamental building blocks: convolutional layers, pooling layers and a fully connected layer.

A convolutional layer convolves an input image with an N  ×  N filter which slides along the image with a stride of K. The pooling-layer groups neighbouring pixels and represents them as a single pixel. For example, in max pooling, the maximum value from these neighbouring pixels is used. This results in an output image which has a smaller dimensionality than the input image. Lastly, a fully connected layer takes the output of several convolutional and max pooling layers and uses this output as features for training a fully connected neural net.

The CNN architecture used in this study was created in Keras [46]. Two pairs of convolutional and max-pooling layers followed by another convolutional layer was applied before the outputs were flattened and sent into a fully connected neural net (table 1). The convolutional layers used an 8  ×  8, 4  ×  4 and 2  ×  2 filter for the first, second and third convolutional layers, respectively, with a stride length of 1 and zero-padding to maintain the same dimensionality. All max-pooling layers were based on 2  ×  2 groupings. All activation functions used in ESCAPE-NET were the rectified linear unit (ReLu) except for the output layer in which softmax was used. A block diagram of the CNN is shown in figure 4.

Zoom In Zoom Out Reset image size

A separate CNN was trained for each rat, using a 3-fold cross-validation approach. nCAPs that were detected from each activity (the 100 trials each from dorsiflexion only, plantarflexion only, and pricking only) were split into training and test sets. These training sets included data augmentation to increase the number of training examples to 10 000 (a combination of detected and augmented nCAPs) for each class. Testing sets only consisted of detected nCAPs (no augmented nCAPs) for each class.

Data augmentation was performed to equalize the number of nCAPs used for training for each activity per rat. This procedure aimed to enhance the training of the CNN by addressing the class imbalance and provided the CNN with more training examples. New synthetic nCAPs were constructed by generating noisy versions of an averaged nCAP for a particular neural pathway. First, an average nCAP was obtained by averaging all nCAPs detected for a particular neural pathway. Second, the noise at each contact was characterized by subtracting the average nCAP from each individual nCAPs and fitting a Gaussian distribution to the resulting noise-only data. These Gaussian distributions were then sampled to create new noisy nCAPs by adding sampled noise to the average nCAP template of a particular neural pathway. This process was based on the assumption that the neural activity is reproducible and that a large portion of the observed variability is due to noise in the nerve cuff recordings, which have low SNR.

ESCAPE-NET was trained for 1000 epochs with an early stopping criterion of no decrease in the cost function over 15 epochs on the validation set (random 10% subset of the training set). ESCAPE-NET was optimized using stochastic gradient descent, with momentum value of 0.9 and calculated using a categorical cross entropy loss. Each convolution layer used 64 filters and the fully connected layer had a 256 node hidden layer and an output of three nodes (one for each class). The first and second convolutional layers were each trained for 25 epochs to initialize their weights before actual training of the full ESCAPE-NET architecture. We additionally investigated combining information from the SE and TE representations (SE  +  TE) of the spatiotemporal signatures. In this case, the outputs of ESCAPE-NET trained separately on SE and TE were concatenated before the fully connected layer and fed into a 512 (instead of 256) node hidden layer with the rest of the network staying the same.

Classification accuracy was quantified on an individual nCAP basis, assessed through the percentage of correctly identified nCAPs. Additionally, the F₁-score was computed for each class. The F₁-score metric is a better indication of ESCAPE-NET's performance than classification accuracy when classes are imbalanced in the test set. For the multi-class problem, the F₁-score was calculated based on the precision and recall of each class (e.g. dorsiflexion, plantarflexion and pricking). The reported F₁-score is the average of the calculated F₁-scores for all classes.

ESCAPE-NET was compared to algorithms developed previously in [31, 32]. The first algorithm quantifies the similarity of a nCAP with average spatiotemporal templates for each class using a bank of matched filters. The pathway whose corresponding filter has the highest maximum normalized output for a given nCAP is chosen as the active pathway. The second method is a random forest classifier that uses the 1D spatiotemporal signatures as input. This method takes random subsets of the data to produce different decision trees whose outputs can then be combined to obtain a classification. The last method was a feedforward neural network that also uses these 1D spatiotemporal signatures for classifying the activity of different neural pathways.

The contact information metric (CIM) previously defined in [43] was developed to evaluate the selectivity of a multi-contact configuration without being biased by the number of contacts or the configuration of the electrode. The CIM is defined individually for each contact using Shannon entropy to quantify the amount of information that the contact provides for discriminating between different neural pathways.

The original metric calculated the entropy of a single probability mass function derived from a histogram of the peak values of the detected nCAPs. However, this peak value may be drastically affected by noise and thus an enhanced definition of the CIM is described here. It uses multiple points from the nCAP and evaluates the entropy from the probability mass function of each point, allowing for more robustness against noise.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{H}_{kl}}=-\underset{i}{\overset{n}{\mathop \sum}}\,{{P}_{{{x}_{i}}}}{{\log}_{2}}{{P}_{{{x}_{i}}}}.\nonumber \end{align} \tag{ 2 }$$

In equation (2), H_kl is the Shannon entropy calculated for the kth contact at the lth time point, x is the histogram of values from the kth contact at the lth time point observed in noisy recordings from M neural pathways, x_i is the ith bin of x, and n is the number of equally spaced bins in x. Bins were selected using the histogram bin optimization method by Shimazaki and Shinomoto [47]. The entropy value is determined by variability in the measurements, which can be due to both variability between pathways and to noise. To deal with this issue, the H_kl's were then normalized by the Shannon entropy calculated from a single neural pathway, within which the variability is the result of noise only. This was accomplished by taking the mean Shannon entropy obtained from each of the M neural pathways. This resulting normalized value provides the CIM for each contact at each time point.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm CI}{{{\rm M}}_{kl}}=\frac{{{H}_{kl}}}{\frac{1}{M}\sum\nolimits_{i}^{M}{{H}_{i}}}.\nonumber \end{align} \tag{ 3 }$$

Where H_i is Shannon's entropy calculated for a single neural pathway for the kth contact and lth time point. The CIM for the kth contact is then calculated by taking the mean of the CIM at all time points (100 time samples) for the kth contact.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm CI}{{{\rm M}}_{k}}=\frac{1}{l}\underset{1}{\overset{l}{\mathop \sum}}\,{\rm CI}{{{\rm M}}_{kl}}.\nonumber \end{align} \tag{ 4 }$$

This formulation allows for CIM values to range from 0 to  ∞  where a higher CIM would correspond to a more informative contact and a lower CIM would correspond to a less informative contact.

Smaller electrodes with fewer contacts may provide greater flexibility in terms of implantation sites and applications. We therefore explored how the classification accuracy was affected by the use of electrode configurations with fewer contacts. Contacts were removed from the original 7  ×  8 configuration, creating 5  ×  8 and 3  ×  8 configurations. In addition, contacts were also removed using the CIM described above, to examine the effects of removing less informative contacts from the spatiotemporal signature on the nCAP discrimination performance (for example, for a 40-contact configuration, only the 40 contacts with the highest CIMs were used).

The final metric investigated was the ability to reconstruct the angle of the rat's ankle from the reconstructed nCAP trains of each neural pathway. Training sets consisted of all nCAPs from the dorsiflexion, plantarflexion and pricking only trials in combination with augmented nCAPs for each activity to reach 10 000 training examples per class. ESCAPE-NET was then evaluated on the test set data from an alternating dorsi- and plantar- flexion recording never seen in the training set. Detected nCAPs in the test set were classified, and both estimated and ground truth nCAP trains reconstructed for each neural pathway. The ground truth was constructed using manual classification of all detected nCAPs based on the synchronized video data.

Ground truth and estimated nCAP firing rates for each neural pathway were obtained by convolving the ground truth and estimated nCAP trains, respectively, with a Gaussian kernel with a standard deviation of 150 ms [48]. These firing rates were then binned in 25 ms intervals and used to train a recurrent neural network (RNN) in Keras [46]. The dataset was split into 70% training and 30% testing. Inputs to this RNN were ten time sample windows (250 ms windows) of the binned ground truth firing rates and ten shifted versions (e.g.  −500 ms to  −250 ms, −475 ms to  −225 ms, ...  −250 ms to 0 ms). These inputs were used to predict the next time sample's angle.

The RNN architecture had two long short-term memory layers (64 units, 128 units respectively), followed by a 1 unit fully connected layer for regression to output the predicted ankle angle. The RNN was trained for 1000 epochs with an early stopping criterion of no decrease in the cost function over 15 epochs on the validation set (10% random subset of the training set). Weights were initialized using a normal distribution. The network used a mean squared error loss function and was optimized using the adaptive moment estimation (ADAM) approach [49].

This trained network was then used to estimate the angle of the rat's ankle from the firing rates of the testing set by using either the estimated firing rate (to obtain the results with the current algorithm), or the ground truth firing rate (to obtain the best theoretical results). Pearson correlation coefficients were used to describe the relationship between the manually labeled ankle angle obtained from annotating recorded videos (with interpolation performed for points in between) and the estimated ankle angles from the trained RNN. Rat 3's data was removed from this analysis due to the degradation of the plantarflexion signal through the course of the experiment.

A repeated one-way analysis of variance (ANOVA) was performed to compare all six algorithms (table 2), and post hoc pairwise comparisons were done with t-tests, including a Bonferroni correction for multiple comparisons. Subsequent evaluations (effects of reference montage and contact configurations) were applied using a pairwise t-test (table 3) or repeated one-way ANOVAs (tables 4 and 5) were performed on the best performing algorithm only (ESCAPE-NET SE  +  TE). Post hoc comparisons were conducted with t-tests, and Bonferroni correction was applied.

The mean and standard deviations of detected nCAPs obtained from all rats for each activity were $5513\pm 2465,~5615\pm 3638,$ and $7513~\pm 4046$ nCAPs for dorsiflexion, plantarflexion and pricking respectively. The minimum and maximum numbers of nCAPs detected were 1950 and 8808, 293 and 9954, and 1164 and 13 180 for dorsiflexion, plantarflexion and pricking, respectively.

Table 2 shows the mean classification accuracy and corresponding mean F₁-score for the algorithms presented previously in [32] in comparison to our new and improved method from all rats in the 3-class discrimination problem (dorsiflexion, plantarflexion and pricking) using a tripolar reference. A significant main effect was seen in the choice of algorithm on the mean classification accuracy and F₁-score (p   <  0.001 for both). In post hoc pairwise comparisons, ESCAPE-NET showed significant (p   <  0.001, except ESCAPE-NET TE versus Random Forest (p   <  0.01)) improvements in both mean classification accuracy and F₁-score over previously reported methods (table 2), regardless of whether the SE, TE or SE  +  TE representations were used. All post hoc pairwise comparisons can be found in supplementary tables I and II (stacks.iop.org/JNE/17/016042/mmedia).

Table 3 shows the mean classification accuracy and corresponding mean F₁-score for our proposed method from all rats in the 3-class problem for different referencing montages (tripolar referencing, and common average). The reference montage used showed no significant differences in reported mean classification accuracy or F₁-score.

Figure 5 shows an example of the predicted firing pattern of dorsiflexion and plantarflexion using ESCAPE-NET SE  +  TE for a subsection of the alternating dorsi-/plantar- flexion recording in rat 6.

Zoom In Zoom Out Reset image size

Table 4 shows the mean classification accuracy and corresponding mean F₁-score for reduced electrode configurations (5  ×  8, 3  ×  8 configurations) from all rats in the 3-class problem using a tripolar reference. The contact configuration showed a significant main effect (p   <  0.001) on both mean classification accuracy and F₁-score. Additionally, a steady significant decrease in both mean classification accuracy and F₁-score (p   <  0.001) between all pairwise comparisons when using ESCAPE-NET SE  +  TE was found. All post hoc pairwise comparisons can be found in supplementary tables III and IV.

Table 5 shows the mean classification accuracy and corresponding mean F₁-score for reduced electrode configurations from removing the least informative contacts using the CIM (52, 48, 44, 40 contacts) from all rats in the 3-class problem using a tripolar reference. Removing contacts with low CIM scores led to a slight decline in recording selectivity and a main effect of removal of contacts was seen in the F₁-score (p   <  0.05). Significant differences were only found between 44 versus 40 contact configurations for both mean classification accuracy (p   <  0.05) and F₁-score (p   <  0.01). All post hoc pairwise comparisons can be found in supplementary tables V and VI.

Continuing to reduce the number of contacts based on CIM, down to ten contacts, led to gradual decrease in F₁-score, but did not reveal a distinct corner in the curve that could be indicative of a critical number of contacts (supplementary figure 1).

The mean Pearson correlation coefficient was 0.812  ±  0.177 for the ankle angle predicted based on the estimated firing rate from ESCAPE-NET compared to the manually labelled ankle angle. To contextualize this value, the theoretical maximum performance of the RNN that could be achieved, if classification of the nCAPs was perfect, was calculated by predicting the ankle angles using the ground truth nCAPs as inputs to the RNN. This resulted in a mean Pearson correlation coefficient of 0.980  ±  0.009. Two examples of the predicted joint angles are seen in figure 6, representative of the best and worst results achieved.

Zoom In Zoom Out Reset image size

Figure 7 shows the classification performance and the mean Pearson correlation coefficient (i.e. ability to correctly predict the joint angle) for each rat except rat 3, which was excluded as mentioned in section 2.7. The figure shows an expected relationship between the two metrics and suggests that an F₁-score of 0.7–0.8 is required for robust tracking in a 3-class problem.

Zoom In Zoom Out Reset image size