We used 8-week-old male Wistar rats (283–464 g; n = 6; Australian Phenomics Facility, Canberra, ACT, Australia). There were 1–3 animals housed per cage, with a 12/12-h light/dark cycle. Food and water were available for animals to access ad libitum. All procedures were approved by the Australian National University Animal Experimentation Ethics Committee (A2014/52) and adhered to the Australian code of practice for the care and use of animals for scientific purposes.

Rats were anesthetized with urethane (1.4g kg⁻¹ i.p.). A tracheotomy was performed, and a breathing tube inserted to aid natural respiration. Rats were placed in a stereotaxic frame with head flexion at −20°; (Stoelting Instruments). The dorsal skin and muscles of the neck, and the dura and arachnoid mater were excised between the foramen magnum and the C1 vertebra to expose the dorsal surface of the brainstem. In most cases, the C1 vertebra was cut away with rongeurs to aid placement of an adapted surface multi-electrode array (sMEA; Nucleus 22 Auditory Brainstem Implant, Cochlear Limited). Details of the electrode dimensions are shown in the inset in Figure 1, and further details of the electrode can be found in Chelvanayagam et al. (2008). A flexible plastic rod held by a micromanipulator on the stereotaxic frame was lightly pressed onto the sMEA to hold it symmetrically over the brainstem midline.

We applied two types each of tactile- and proprioception-dominated stimuli to the rat limbs amounting to 16 (4 stimuli × 4 limbs) possible mechanical stimulus conditions. The tactile-dominated stimuli were applied by pressing either a wooden dowel rod (diameter: 2.0 mm) or brush into the palmar/plantar pads of the paws. The brush was 25 mm in length, with diameter: 5.0 mm that tapered to a point over the last 10 mm. The rod and brush were fixed at the end of a flexible tube to deliver 20 g of force. The proprioception-dominated stimuli were applied to the rat by flexing or extending the rat’s limbs. The rats’ limbs were fixed to a dowel rod with ethyl cyanoacrylate instant adhesive (Loctite^®, Henkel Adhesive Technologies) which enabled the experimenter to directly manipulate the limbs. Both flexion and extension involved the movement of the hip, knee, and ankle joints for the lower limbs, and the shoulder, elbow, and wrist joints in the forelimb. Therefore, the proprioceptive-dominated stimuli activated afferents across the entire limb, spanning three joints, and the associated muscles, in addition to the cutaneous afferents activated throughout the limbs by the movement. We do not describe the stimuli as exclusively tactile or proprioceptive, because the application of the tactile-dominated stimuli might also activate receptors in the joints and intrinsic muscles of the paw, and the proprioception-dominated stimuli undoubtedly activated tactile receptors associated with hair and skin throughout the limbs, in addition to the proprioceptive receptors in the skin, muscles, tendons, and joints. All stimuli were applied by the experimenter who was triggered by the timing of a metronome.

Stimuli were applied for 2.4 s with a rest period of 2.4 s. Tactile-dominated stimuli had a stimulus on/off pattern (4.8 s per trial), but proprioception-dominated stimuli also had a 2.4 s period which was used to return the limb to the resting state, where it remained in resting state for 2.4 s (7.2 s per trial). Stimuli were applied in 10 sets of 10 trials (100 stimuli per type), with 30 s rest between sets.

DCN electrical signals were acquired through seven electrodes of the sMEA and filtered (50 Hz notch filter; 10 kHz low-pass filter) through custom-built amplifiers. The signals were then digitally recorded (40 kHz sample rate) through a PowerLab 16/35 acquisition system and viewed in LabChart Pro software (Version 8.1.1, AD Instruments, Bella Vista, NSW, Australia). The seven electrodes of the sMEA are referred to from rostral to caudal as follows: left side electrodes: e1, e2; midline electrodes: e3, e4, e5; right side electrodes: e6, e7 (Figure 1A, right insert), and all seven electrodes combined are referred to as e1–7.

The large electrodes in this surface array were able to cover most, if not all, of the gracile nuclei, approximately three-quarters of the cuneate nuclei, and parts of the external cuneate nuclei. These electrodes allow acquisition of neural activity from large populations of neurons simultaneously and thereby permit finding the most salient and information-rich signal features that dominate somatosensory-evoked activity across most of the DCN, and how these features vary over time.

Signal processing, feature extraction, and analysis were performed offline (MATLAB version R2018a, MathWorks). We used 22 DCN signal features as artificial neural network (ANN) inputs. The names, descriptions, and examples of feature extraction for all 22 features are shown in Figure 2. Features were categorized into LF (<200 Hz) and HF (>200 Hz) bands for analysis. Others have used a similar frequency (usually between 100 and 300 Hz) to separate LF activity that typically represents local field potentials, calcium spikes, and intrinsic neural membrane potentials, from HF activity that typically represents a single unit or multiunit spiking activity (Buzsáki et al., 2012). Four features were extracted from HF filtered (bandpass 0.55–3.3 kHz; 5-order Butterworth filter) signals (HF features; Figures 2A–D), and five from LF signals (LF features; Figures 2E–I), which were quantified from the HF signals after they had been rectified and low-pass filtered (<80 Hz; 5-order Butterworth filter). For all time-domain signals in which peak amplitudes were measured, we used the prominence output variable from the findpeaks MATLAB function. Thirteen features were extracted from frequency spectrograms (spectrogram MATLAB function) of DCN signals bandpass filtered between 4–5,000 Hz (8-order Butterworth). Eight of these features were quantified from the peak power spectral density of HF bands (HF PSD features): 200–600 Hz, 600–1,000 Hz, 1,000–1,500 Hz, 1,500–2,000 Hz, 2,000–2,500 Hz, 2,500–3,000 Hz, 3,000–3,500 Hz, 3,500–4,000 Hz (Figure 2J). The other five spectral features were quantified from the peak power spectral density of LF bands (LF PSD features): 4–8 Hz, 8–13 Hz, 13–40 Hz, 40–80 Hz, 80–200 Hz (Figure 2K). For all filtering, we used a zero-phase response filter: filtfilt function (MATLAB).

HF and LF features contain multiunit activity information, which in some cases outperforms spiking activity or local field potentials in offline decoding of motor tasks (Stark and Abeles, 2007) and has been used to control motor neural prostheses in tetraplegic patients (Flint et al., 2013; Bouton et al., 2016). The HF PSD features have been minimally investigated, but similar frequency band ranges have been used to decode motor signals in a neural prosthesis capable of effecting limb movement through neuromuscular stimulation (Bouton et al., 2016). The LF PSD features have been widely investigated in the motor cortex during reaching and grasping movements or used as decoding features for electrocorticographic or electroencephalographic neural prosthetic control (Wolpaw et al., 1991; Kostov and Polak, 2000; Leuthardt et al., 2004; Rickert et al., 2005; Flint et al., 2012; Chen et al., 2013; Marathe and Taylor, 2013; Bundy et al., 2016).

Individual signal features extracted from each of the seven electrodes were paired to the stimuli that generated them, to create input/output pairs for machine-learning. All machine-learning experiments used a standardized ANN with a supervised learning classification algorithm (patternnet MATLAB function). The standardized ANN comprised 42 hidden neurons and 16 output neurons corresponding to the 16 possible stimuli (four stimulus types, presented to four different limbs). We selected the number of hidden neurons by increasing their number between 16 (the number of outputs) and 154 (the highest possible number of inputs from our feature set). For each number of hidden neurons, we determined the average feature-learnability with the minimum (7) and maximum (154) number of features included in the input set and found that 42 hidden neurons produced the highest average feature-learnability. Thus, only the input neurons were altered, depending on the number of input features required for each experiment. The ANN used gradient descent with momentum and the adaptive learning rate backpropagation training function traingdx. The hyperbolic tangent sigmoid function, tansig, was used in the hidden layer units, and a softmax transfer function (softmax) was used in the output layer. Both inputs and output targets were normalized such that all values fit between −1 and 1, as per the patternnet default setting. For training, cross-validation, and testing, the inputs were separated into training, validation, and test subsets, using the patternnet default settings.

Feature-learnability provides a stable measure of relevant information content provided by input features including a measure of the information content variability (Loutit et al., 2019). We used the Within Individual Animal (WIA) approach (Loutit et al., 2019) for all feature-learnability testing unless otherwise specified. In this approach, input/output pairs from 1,600 stimuli were generated for each animal. For each animal, ten repeated training (70%), validation (15%), and testing (15%) machine-learning cycles were performed under random initializing conditions. These ten confusion matrices were averaged to produce a single representative confusion matrix for each animal, then the representative confusion matrices for all animals were averaged. Feature-learnability was determined by the mean ± SEM derived from the diagonal of this matrix. To determine the feature-learnability of an individual feature, the inputs to the standardized ANN were restricted to the single input feature in question, resulting in a total of seven inputs (the signal feature extracted from each electrode).

To obtain a measure that represents the maximum possible information content contained in signal features, a benchmark feature set was determined that produces the highest feature-learnability outcomes. This feature-learnability benchmark facilitates comparisons of information content contained by individuals or subsets of input features. Individual features were extracted from all seven electrodes over the first 1,000 ms from stimulus onset, and feature-learnability determined. Individual features were then ranked from highest to lowest feature-learnability, prioritized by the largest means and smallest SEM (Loutit et al., 2019). Input feature selection was performed using an adapted sequential forward floating search algorithm (Whitney, 1971; Pudil et al., 1994). Our approach differed from a typical sequential forward floating search algorithm as we added four features (one from each feature category) to the input set per round of feature-learnability testing, rather than one feature. We also used a sequential backward search in each round to test whether removing features improved feature-learnability. The search was continued until the possible feature combinations presented by this method were exhausted and the maximum feature-learnability was determined.

A second benchmark was determined for features extracted over an optimal short-time window length (described below). The identical approach was used as for the 1,000 ms window benchmark, except that features were extracted from the shorter time course.

To determine the generalisability of the features across animals we used the 1,000 ms window benchmark feature set to compare feature-learnability using a leave-one-out (LOO) approach. We randomly assigned input/output pairs from five animals into training (70%) and validation (30%) sets, and testing was performed on the remaining animal (100%). The LOO approach was applied such that each of the six animals was examined as the test data set once. For more details on these methods see Loutit et al. (2019).

We sought to determine feature-learnability based on a shorter time-window for two purposes. Firstly, the features we investigated may be useful for decoding neural signals unrelated to DCN-specific signal features, such as neural signals from motor systems. Typical motor decoding algorithms use time windows ranging from 50–100 ms (Lebedev and Nicolelis, 2017), which trades classification accuracy for a feasible reduction in time required for real-time applications. Secondly, some peripheral somatosensory neural prostheses use stimulus patterns that mimic natural peripheral afferent firing properties, like variations in firing rate shown by fast and slowly-adapting afferents at the onset and offset of a stimulus. These biomimetic approaches have succeeded in eliciting more natural somatosensory percepts than simpler stimulus patterns that use only amplitude or frequency modulation to encode stimulus intensity (Valle et al., 2018; George et al., 2019). Therefore, determining feature-learnability over shorter time windows may help determine which features could inform the construction of artificial stimulus patterns during different stimulus events, such as contact onset, static holding, and contact offset. Extraction of features over longer time windows might mix those encoding onset, hold, and offset, and therefore be less informative to a decoding algorithm.

To determine an optimal short-time window that provides the highest feature-learnability in the shortest time, we tested windows between 20 ms to 150 ms with 10 ms increments, 250 ms, and 500 ms, in addition to the 1,000 ms windows described above. All windows started at the time of stimulus onset. We plotted feature-learnability against time windows for each feature and used the findchangepts MATLAB function to determine the most abrupt change in feature-learnability. We then averaged the changepoint values across all features and rounded to the nearest time window. A benchmark set of signal features was determined for the optimal short-time window (described above).

Once we established the optimal short-time window, we sought to use the short rolling window to extract features throughout stimulus presentation and determine how feature-learnability varied over time. To do this we extracted the 22 features from the optimal short-time window, sampled at 1,000 ms before the stimulus, then every 10 ms from 200 ms pre-stimulus to 4,000 ms post-stimulus. Finally, we tested feature-learnability for each of the features extracted from each of the 421 windows.

We used R (version 3.4.4; R Core Team, 2018) for all statistical analyses with the RStudio integrated development environment (version 1.1.442). For comparison of the feature-learnability benchmark within and across animals we used a linear model (LM; R lm function, stats package). For all other comparisons, we used a linear mixed-effects model (LMER; lmer function from the lmerTest package; Kuznetzova et al., 2019). Multiple comparisons were performed using estimated marginal means comparisons with Tukey p-value adjustments [R emmeans function, emmeans package (Lenth et al., 2019)]. We verified the model fitting by examining model residuals for normality. Where required, log or logit transformations (car package; Fox and Weisberg, 2019) were applied to the data before statistical modeling. All data are expressed as means ± SEM unless otherwise stated. Probabilities of p < 0.05 were deemed significant.

Preliminary observations of the data revealed neural activity was greatest on midline electrodes and those ipsilateral to the site of stimulus. Tactile-dominated stimuli evoked activity that was greatest at stimulus onset and/or offset. Dowel stimuli generally evoked a short and sharp burst of neural activity that peaked within 10 ms of stimulus contact/removal, whereas the brush stimuli evoked a longer, ramped burst of neural activity that was comparatively delayed at onset/offset and to reach maximum (20–30 ms). In some cases, brush stimuli evoked two initial bursts at stimulus onset. Proprioception-dominated stimuli evoked more neural activity than tactile-evoked stimuli. In general, flexion resulted in greater neural activity compared to extension, however, the time to reach maximum neural activity was similar for both proprioception-dominated stimuli, which was at approximately the midpoint of the movement. Examples of filtered signals (0.55–3.3 kHz) acquired from seven electrodes in response to the four types of stimuli applied to the left forelimb are shown in Figure 1B. Features representing neural signals are shown in Figure 2.

We sought to rank individual features for their capacity to be informative of stimulus type and location. Figure 3A shows the rank order of feature-learnability for all 22 features extracted from the first 1,000 ms following the stimulus onset across all seven electrodes. All 22 features performed significantly greater than chance levels of 6.25%. Feature-learnability ranking clustered into three groups: (1) highest-performing features, comprising three of the four HF and one of the five LF features; (2) middle performing features, comprising the remaining LF and all but one of the HF PSD features; and (3) the lowest-performing features, comprising mainly LF PSD feature and the remaining HF and HF PSD features, all of which still performed three times greater than chance levels. These findings indicate that time-domain HF features are the most informative for determining a combination of stimulus location and type.

We previously demonstrated that combinations of LF and HF signal features improve machine-learning outcomes (Loutit et al., 2019). We, therefore, sought to determine a combination of input features that produced the highest feature-learnability and to establish a benchmark for subsequent comparisons. Combinations derived from the best performing pair of features from within the HF, LF, and HF PSD categories resulted in significant feature-learnability improvements, compared to single features (p ≤ 0.01, LMER, Tukey), but tended to plateau after adding the second feature (Figure 3B). The best performing pair of frequency-domain features in the LF PSD category was not better than the best performing single LF PSD feature (p = 1.0, LMER, Tukey). The highest-ranked learnability from within-category feature combinations was achieved by combining all four HF features (86.5%). However, this HF 4-feature combination was not significantly greater than the 2 or 3-feature HF combinations (p = 1.0, LMER, Tukey), all of which also significantly out-performed the best individual HF feature and all other within-category combinations (p ≤ 0.0003, LMER, Tukey; Figure 3B).

To determine if combinations of features from different categories improved feature-learnability, we combined the highest-ranked feature from each category (i.e., a combination with 4 features). This combination yielded significantly greater feature-learnability than all LF, LF PSD, and HF PSD highest-ranked feature combinations (p ≤ 0.0063, LMER, Tukey). Compared to the HF feature category, the 4-feature across-category combination was only significantly greater than the highest-ranked single HF feature (p < 0.0001, LMER, Tukey), but was not significantly different when additional HF features were added (p = 1.0, LMER, Tukey; Figure 3B).

To find the combination with the highest-ranked feature-learnability, we sequentially added the next best features in groups of four (i.e., the next best individual feature from each of the four categories). The highest-ranked feature-learnability achieved was 87.2 ± 1.3% with 13 features which we defined as our 1,000 ms window feature-learnability benchmark for subsequent comparisons. Despite improved feature-learnability with 13 features compared to 4, 8, or 11 features, there were no significant differences between feature-learnability outcomes of any of the across-category combined feature sets (p ≥ 0.66, LMER, Tukey; Figure 3B). The individual features that contribute to the 1,000 ms benchmark combination are indicated in Figure 3A, and their confusion matrices are shown in Figures 4A–D.

In summary, time-domain HF features resulted in the highest feature-learnability rankings, and although the benchmark feature set included 13 features from all four categories, benchmark feature-learnability was not significantly higher than the combination of two HF features (HF spike amplitude and HF spike count).

To determine what information the 1,000 ms benchmark input features contribute to feature-learnability, we plotted confusion matrices for all 13 features (Figure 4). Correct predictions (i.e., the diagonal) of each matrix are replotted with their SEM in the right panels of Figures 4E–H to facilitate performance comparisons of each feature and to provide a measure of variability among animals.

The HF feature category (Figure 4A) significantly outperformed all other categories (p ≤ 0.001, LMER, Tukey). Proprioception-dominated stimuli were significantly better classified than tactile-dominated stimuli by all categories (p < 0.0001, LMER, Tukey), except the LF PSD category (Figure 4H, p = 0.74, LMER, Tukey). The forelimbs were significantly better classified than hindlimbs across all categories (p ≤ 0.004, LMER, Tukey), and for both tactile- and proprioceptive-dominated stimuli (p < 0.0001, LMER, Tukey). HF spike amplitude was the best predictor of tactile-dominated stimuli and significantly outperformed all other features (p ≤ 0.0026, LMER, Tukey; Figures 4E–H). HF integral was the best predictor of proprioceptive-dominated stimuli and significantly outperformed most other features (p ≤ 0.0019, LMER, Tukey), except HF spike count and LF sum burst amplitudes (p 1.0, LMER, Tukey; Figures 4E–H). Interestingly, the LF burst count predicted proprioceptive-dominated stimuli significantly better when evoked from right limbs compared to left limbs (p = 0.035, LMER, Tukey; Figure 4F).

To investigate how HF spike amplitude predicted tactile stimuli significantly better than all other features, including the next highest performing feature HF spike count, we quantified HF spike amplitude and compared this to HF spike count acquired from each electrode and animal. We previously demonstrated that a feature’s learnability is correlated to the number of instances different stimuli evoke significantly different magnitudes of that feature (Loutit et al., 2019). We, therefore, quantified the two features from anatomically relevant electrodes, i.e., features acquired from stimuli applied to left forelimbs were quantified from left electrodes (e1 and e2), right forelimbs from right electrodes (e6 and e7), and hindlimbs from midline electrodes (e3, e4, e5).

Although the effect size was small, HF spike amplitude evoked by dowel stimuli (30.2 ± 1.2 μV) were significantly higher than brush stimuli (29.7 ± 1.0 μV; p = 1.3e-4, paired t-test), as were HF spike counts (dowel, 41.6 ± 2.5 events; brush, 40.0 ± 1.9 events; p = 0.041, paired t-test). HF spike amplitudes evoked by flexion (32.4 ± 1.2 μV) were significantly higher than when evoked by extension (31.1 ± 1.2 μV; p = 2.1e-12, paired t-test), as were HF spike counts (flexion, 129.1 ± 8.2 events; extension 87.9 ± 5.9 events; p = 3.3e-14). Moreover, HF spike amplitude of proprioception-dominated stimuli (31.7 ± 0.5 μV) was significantly higher than when evoked from tactile-dominated stimuli (29.8 ± 0.5 μV; p = 0.013, Student’s t-test), and proprioception-evoked HF spike counts (108.5 ± 5.7 events) were more than 2.5 times larger than when tactile-evoked (40.8 ± 2.0 events; p = 8.4e-24, Student’s t-test).

In summary, despite small effect sizes, the difference in the HF spike amplitude evoked from dowel vs. brush stimuli resulted in a much smaller probability value compared to HF spike counts, which may account for the improved tactile classification outcome from HF spike amplitude. Dowel, flexion, and proprioception-dominated stimuli evoke higher HF spike amplitudes and HF spike counts than a brush, extension, and tactile-dominated stimuli, respectively.

How generalizable are signal features of the 1,000 ms benchmark across different animals? To answer this question, we used the LOO approach to measure how features extracted from an individual animal perform when presented to a neural network trained by the same features derived from the remaining cohort of animals. This approach measures how features from one animal generalize to all other animals (Loutit et al., 2019). Compared to the WIA approach derived from neural networks optimized for individual animals using 13 features over 1,000 ms (Figure 5A), feature-learnability determined by the LOO approach (Figure 5B) was significantly reduced, by almost half to that of the WIA approach (Figure 5C; p < 2.2e-16, LMER). More than 50% of the reduction in feature-learnability derived by the LOO approach was accounted for by greater confusion errors associated with tactile-dominated stimuli, specifically, that dowel and brush stimuli were generally poorly discriminated, which were exacerbated by left/right errors for the hindlimb (Figures 5B,C). Compared to the WIA approach, proprioception-dominated stimuli from the LOO approach demonstrated an insignificant 11% reduction in classification accuracy at the forelimb (p = 0.25, LMER, Tukey), but a 43% reduction associated with the hindlimb (p < 0.0001, LMER, Tukey; Figure 5C). These results indicated that the 13 features contain information that is highly generalizable among animals for proprioceptive-dominated stimuli of the forelimb, but much less so for hindlimb proprioception and forelimb and hindlimb tactile-dominated stimuli.

Combining the two highest-ranked HF features (Figure 3) resulted in WIA feature-learnability (85.6 ± 1.2; Figure 5D) not significantly different from the combination of 13 features (87.2 ± 1.3; Figure 5A, p = 0.13, LMER). To determine how these two features alone generalize across animals, the LOO approach was applied (Figures 5E,F). The LOO approach restricted to the two highest-ranked HF features revealed an almost identical pattern as the 13 features data set, but with the addition of forelimb/hindlimb confusion for tactile-dominated stimuli (Figures 5E,F). Compared to the 13-feature combination, LOO feature-learnability of the 2-feature combination was overall significantly reduced (p = 0.002, LMER; Figure 5F), which resulted because of reduced tactile stimulus classification accuracy that was generalized across both hind- and forelimbs (p = 0.025, Tukey; Figures 5E,F). Interestingly, there was no significant reduction in proprioception performance in the LOO approach between the 13- and 2-feature sets (p = 0.37, LMER, Tukey).

These results indicate that the two HF features extracted over 1,000 ms can provide almost all of the information provided by the larger 13-feature set for forelimb proprioceptive-dominated stimulation, but the other 11 features contribute some additional information, common across all animals, which is important to discriminate hindlimb proprioceptive-stimuli as well as fore-and hindlimb tactile-dominated stimuli.

Features extracted over a 1,000 ms period have little practical value for neural prosthetic feedback control. A DCN-targeted neural prosthetic device would need to provide stimulus features over shorter time windows to enable rapid updating of limb sensory status. To determine how learnable DCN signal features are over shorter periods, feature-learnability was calculated for each signal feature extracted over time windows ranging from 20–500 ms starting from stimulus onset (Figures 6A–D). A window of 60 ms was determined as the time of most abrupt change in feature-learnability across all features (see arrows, Figures 6A–D) before reaching a plateau, where feature-learnability improved slightly, or not at all until much larger time windows (≥250 ms). Feature-learnability rank order derived from features extracted from 60 ms windows (Figure 6E) was altered compared to 1,000 ms windows (Figure 3A). Notably, HF integral out-ranked HF spike amplitude (the best feature from 1,000 ms windows), although these features were not significantly different (p = 0.052, LMER, Tukey; Figure 6E), and the second-ranked feature was LF sum burst amplitudes (ranked 4^th from 1,000 ms windows). HF’s mean width remained the lowest-ranked feature.

Combining the best two features (HF integral + LF sum burst amplitude) did not significantly improve feature-learnability (41.8 ± 2.2) compared to HF integral alone (41.0 ± 2.5, p = 0.09, paired t-test). A combination of four signal features that included the highest-ranked feature from each category extracted from 60 ms, significantly improved feature-learnability compared to most within-category combinations, except for the combination of two or more HF features (p ≥ 0.76) and five or more HF PSD features (p ≥ 0.35, LMER, Tukey; Figure 6F). Also noteworthy was that feature-learnability from combining HF integral and HF spike amplitude was significantly greater than all other within-category feature combinations, except for three or more combined HF PSD features, all of which were not significantly different (p ≥ 0.6, LMER, Tukey). The 60 ms benchmark feature set, determined by the continual addition of the next best features from each category (where possible), resulted in a feature-learnability outcome of 59.4% from 17 features (gold arrows, Figures 6E,F).

To measure how informative signal features, extracted over a short period, are throughout stimulus presentation, we determined a feature-learnability time-series by extracting DCN signal features from 60 ms rolling windows (Figure 7). We started by examining the feature-learnability temporal profile for the 60 ms benchmark 17-features set (indicated by gold arrows, Figure 6E). This produced a complex waveform (gold trace, Figure 7A) with three distinct peaks: the 1st is an abrupt, relatively sharp peak coinciding with the stimulus onset; the 2nd and 3rd peaks were broadened and peaked at the approximate midpoint of the stimuli and beginning of the stimulus-off/return phase, respectively. To establish which stimuli contributed to these peaks, neural networks were restricted to input/output data sets of tactile- or proprioception-dominated stimuli (Figure 7A, red and blue traces, respectively). This revealed that the 1st peak arose from tactile-dominated stimuli, the 2nd peak arose from proprioceptive-dominated stimuli, whereas the 3rd peak arose from both tactile- and proprioceptive-dominated stimuli. Another interesting observation was that during the stimulus-off/return phase, feature-learnability was significantly elevated from chance levels, indicating that the input features during this period were informing the neural network of some information about the stimulation being performed.

To determine how individual features contribute to feature-learnability over this time course, we examined feature-learnability time profiles for each of the 60 ms benchmark features (Figures 7B–D). For all features, the 3rd peak coincided with the commencement of the off-stimulus period and their amplitudes never exceeded their respective 1st or 2nd peaks. Of the four time-domain HF features, the maximum peak for three of these coincided with the 1st (tactile) peak, whereas the HF integral maximum fell on the 2nd (proprioceptive) peak (Figure 7B). Feature-learnability returned to near-chance levels during the stimulus-off/return period for all HF features except HF spike amplitude, which remained significantly elevated. Of the five time-domain LF features, the maximum peak for two of these coincided with the 1st peak (LF max burst amplitude and LF mean burst amplitude), whereas the other three features had similar 1st and 2nd peak magnitudes (Figure 7C). Feature-learnability returned to chance levels during the stimulus-off/return period for all LF features.

Of the eight HF PSD features, the maximum peak coincided with the 1st peak for two features derived from frequencies between 0.2–1.0 kHz, whereas maximum feature-learnability coincided with the 2nd peak for features derived from frequencies between 1.5–4.0 kHz (Figure 7D). All HF PSD features remained elevated above chance levels during the stimulus-off period, with frequencies closer to the 2 kHz range demonstrating greater feature-learnability during the stimulus-off/return period. For all five LF PSD features, the maximum peak coincided with the 1st peak, and all remained slightly elevated above chance levels during the stimulus-off/return period (Figure 7E).

To gain insight into why feature-learnability remained above chance levels during the stimulus-off/return period, we inspected confusion matrices from all animals for features associated with above-chance performance at 1-s pre-stimulus. Figure 8 shows examples of two features (HF spike amplitude and HF PSD 2.0–2.5 kHz), which demonstrates that increased feature-learnability 1 s before stimulus arose for different reasons across different animals, but similar reasons within animals. One common theme was that the correct limb was identified for the two stimulus categories. In all but one animal, one or more outputs were correctly classified with >35% accuracy. To investigate if there were learning differences at the beginning vs the end of trials, we then divided data sets into thirds for: (i) each of the 10 trials; and (ii) for the 100 sequential stimulus presentations, to see if repeated stimuli within and across the trials contributed to improved outcomes. We found no difference in feature-learnability between the first-third and last-third of stimulus presentation within-trial sets, or across all 100 trials, indicating that there were no changes to feature-learnability as a result of repeating trials.

In summary, for 60 ms rolling windows some features performed better over periods rich in tactile-stimulus information and others over the proprioception-rich periods. PSD features derived from <1kHz, contributed more during tactile-rich periods, whereas that >1.5kHz contributed more toward proprioception-rich periods, while feature-learnability during the stimulus-off/return period was greatest for PSD features around 2 kHz and reduced for lower and higher frequencies. The reasons for above-chance performance during the stimulus-off/return period appears non-generalizable because it was different for each animal.

Across all individual features and our benchmark feature set, proprioception-dominated stimuli were generally predicted better than tactile-dominated stimuli, and forelimbs were predicted better than hindlimbs. More accurate proprioception classification compared to tactile is likely to have resulted because there was a clear significant separation in the neural activity evoked by flexion and extension (indicated by the HF spike count feature alone), and both these proprioceptive stimuli evoked activity that was significantly different to both tactile stimuli, whereas, the significant difference and effect size was much less between the two tactile stimuli. Greater proprioception evoked activity may have resulted from activating more receptors and/or evoking more action potentials from each activated receptor, as described above. Furthermore, tactile stimuli evoked most activity at stimulus onset, then appeared to quickly adapt (Figures 1B, 7A), whereas moving the limb evoked high activity levels throughout the entire duration of both proprioceptive stimuli, thus the potentially higher number of proprioceptors were also active for a longer duration.

Forelimbs may have been predicted better than hindlimbs because most hindlimb activities were acquired by the same midline electrodes which spanned across the two gracile nuclei on both sides, whereas each cuneate nucleus on either side had its electrode, thereby facilitating better spatial discrimination. Some of the hindlimb errors resulted from confusing tactile- and proprioceptive-dominated stimuli of the same limb (Figures 4A–D, 5B,E). Most hindlimb proprioceptive afferents either project onto DCN neurons in the ventral gracile nuclei or to nuclei X and Z, which are rostral nuclei that form part of the DCN-complex (for a comprehensive review see Loutit et al., 2020). Although both hindlimbs tactile and proprioceptive afferents project to the DCN, the proprioceptive DCN regions may be too deep to acquire HF features, while nuclei X and Z were not covered by the placement of our electrode array. Thus, prediction of hindlimb stimuli may have relied mostly on tactile information to discriminate between the tactile and proprioceptive-dominated stimuli and were confined to midline electrodes, thereby reducing overall lower hindlimb prediction accuracy. Meanwhile, the external cuneate nuclei, part of the DCN-complex that exclusively receives forelimb proprioceptive afferents (Loutit et al., 2020), were located partially under the lateral electrodes. Thus, tactile evoked forelimb activity was mainly acquired from both lateral and midline electrodes, whereas proprioceptive evoked forelimb activity was acquired mainly from lateral electrodes. The segregation of proprioceptive information across separate electrodes for forelimbs, but not hindlimb-evoked activity, is likely to explain the reduced feature-learnability derived from hindlimb stimuli.

To determine if benchmark features are generally useful for robust decoding of mechanical stimuli in different animals and not unique to individual animals, we quantified feature-learnability using the LOO approach. This approach trains the machine-learning algorithm on features extracted from all other animals and tests on the remaining animal. We previously established that comparing feature-learnability under LOO and WIA conditions provide insight into how well, or poorly, DCN-signal extracted features generalize across animals (Loutit et al., 2019).

Forelimb flexion and extension were minimally perturbed under LOO conditions, indicating that 1,000 ms benchmark features were highly relevant across all animals for forelimb proprioception stimuli. However, LOO conditions reduced feature-learnability by 37%, which was mostly due to a more than 50% reduction in the ability to predict tactile stimuli and a 40% reduction in the ability to predict hindlimb proprioception-dominated stimuli. This indicates that a significant portion of the information contained within benchmark features are no longer relevant, or as generalizable, for tactile- and hindlimb-presented stimuli across the different animals.

Hindlimb stimuli classification errors (both tactile- and proprioception-dominated) mainly resulted from confusing the side of the body that a stimulus was presented. Signal asymmetry across the DCN surface, as we have previously shown for hindlimb stimulated nerves (Loutit et al., 2017, 2019), could lead to larger signal variations in left and right limb-derived activity, across the different animals we observed in the present study.

HF spike amplitude was the best performing feature for the 1,000 ms windows, and fourth in the 60 ms window, yet, to our knowledge, this signal feature has not been previously used for neural decoding. The HF spike amplitude feature resulted in fewer errors when classifying brush and dowel stimuli of the same limb, compared to all other features (Figure 4), which implies that the machine-learning algorithm detected greater contrast in this feature’s magnitude under the two tactile stimulus conditions compared to other features. We observed that at stimulus onset, dowel stimuli produced more precisely timed bursts of activity than brush stimuli. This may have facilitated the improved dowel/brush classification from the HF spike amplitude feature because spikes from multiple neurons or afferent fibers arriving at an electrode simultaneously (i.e., responses evoked by dowel stimuli) will summate, and therefore show higher average peak amplitudes than single spikes that have less temporal overlap (i.e., brush stimulus-evoked responses). This explanation is supported by the observed significantly greater HF spike amplitude evoked by a dowel, compared to brush stimuli.

The combination of HF spike amplitude and HF spike count features, extracted over 1,000 ms, was sufficient to achieve feature-learnability performance equivalent to the 87% feature-learnability benchmark. Furthermore, this 2-feature input configuration resulted in no confusion errors under WIA, and very few under LOO conditions, between the tactile and proprioception stimulus categories (i.e., no and few confusion errors in the lower left and upper right quadrants of Figures 5D,E respectively). This result can be explained by the significant difference for HF spike amplitude and the >2.5-fold difference in HF spike count when tactile-, compared to proprioception-dominated stimuli were used.

For a fixed time-window, HF spike count provides information about spike frequency and the total number of action potentials generated by peripheral afferents. The large difference in this feature between tactile- and proprioception-dominated stimuli may have resulted from the relatively small vs. larger receptive field sizes respectively that these stimuli engaged. HF spike amplitude captures spiking temporal alignment (as discussed above) and spike magnitudes. Spike magnitudes are influenced by neuron and or axon size, as well as the event’s distance from the electrode (Gasser and Grundfest, 1939; Hunt, 1951; Nelson, 1966; Buchwald and Grover, 1970; Grover and Buchwald, 1970). Group I proprioceptive afferents generally have larger diameters than Aβ tactile afferents (Gasser, 1941). In addition to the laterally placed external cuneate nucleus, proprioceptive afferents preferentially terminate deep on ventral DCN neurons (Campbell et al., 1974), which have a high proportion of large somas (Cheema et al., 1983) that are located 500–700 μm below the brainstem surface in rats (Li et al., 2012). Tactile afferents terminate on smaller neurons in a cluster zone that is located at approximately half the depth as the ventral DCN neurons (Li et al., 2012). How these two neuronal populations (larger somas located deeper vs smaller somas located more superficially) affect spike magnitudes measured from the surface is difficult to assess without dedicated experiments, however, these anatomically segregated neuronal populations likely evoke different shaped spike waveforms. The combination of stimulus-evoked spike-timing, the signal location, and neuronal types responsible for generating the signal may, therefore, result in unique amplitude profiles for each of the stimuli, and thereby contribute to improved stimulus classification.

LF features are measures of the rectified signal envelope and therefore capture various aspects of spike bursts. For example, LF sum burst amplitude, LF max burst amplitude, and LF burst count provide information about burst frequency and magnitudes combined, the largest burst, and the frequency of bursts respectively. The LF feature category demonstrated varied feature-learnability, with its best performer, LF sum burst amplitudes, ranking in the top four features in the 1,000 ms and 60 ms windows. This feature likely represents very similar activity to HF integral, as both are extracted from rectified HF signals, and summing envelope peaks captures related information to the integral of the signal. Congruently, LF sum burst amplitudes and HF integral showed very similar feature-learnability characteristics (compare Figures 2A, 6E, 7B,C) and confusion errors (Figures 4A,B), and combining these two features did not significantly improve feature-learnability from HF integral alone. These features may capture a combination of HF events, slower synaptic activity, or subthreshold events, representing a measure of the signal energy in a neural population.

The HF PSD feature set was similar to that used by Bouton et al. (2016) who extracted features over 100 ms windows for motor signal decoding in a brain-machine interface capable of effecting limb movement through neuromuscular stimulation. This feature set represents frequency power in the multiunit activity range, which is typically considered to be about 300–6,000 Hz (Stark and Abeles, 2007; our range: 200–4,000 Hz). We found these features generally to be good predictors of somatosensory DCN signals. In our 60 ms window data, the frequency bands of 1.5–2.5 kHz were of interest because: (1) the highest decoding for both tactile- and proprioception-dominated stimuli were found for this feature category over this range; and (2) frequency bands above 1.5 kHz showed a peak in the proprioception-dominated phase, while bands below 1.5 kHz, including the LF PSD features, had peaked in the tactile-dominated phase. Each of the eight HF PSD bands appeared to add some unique information not captured by the other bands, as successive additions of these features continued to improve feature-learnability until all the HF PSD features were exhausted (see the yellow line, Figure 6F). However, feature-learnability derived from all eight HF bands was not significantly greater than combining HF integral and HF spike amplitude, indicating that these features could offer superior neural decoding compared to all HF PSD features.

Figure 7 provides insight into how neural activity captured by individual features is altered within 60 ms time windows across the dynamic stimuli. The best performing features over a specific period of the stimuli indicates that there are statistically different neural spiking behaviors (Loutit et al., 2019) across the different stimuli for the 60 ms window over which the feature was quantified. Tactile-dominated stimuli evoked small bursts of activity at stimulus onset, while proprioception-dominated stimuli evoked little activity at stimulus onset, but larger amounts during the middle of the stimulus and return periods on anatomically relevant electrodes (Figure 1B). These phases of activity were captured by the peaks of the feature-learnability temporal profile (gold curve, Figure 7A). The 1st peak of this temporal profile reflects the period where considerable tactile information is presented to the ANN because this peak is abolished when tactile data were omitted from the learning algorithm (blue curve, Figure 7A), but remains when tactile data is present and proprioceptive data were omitted (red curve, Figure 7A). The same logic dictates that the 2nd peak is attributed to the duration when maximum proprioceptive information is present. Interestingly, the tactile-dominated feature-learnability temporal profile coincides with the stereotyped characteristics of fast and slowly adapting afferent activity, whereas the proprioception feature-learnability temporal profile is consistent with the dynamic movement of the proprioceptive stimuli, with its peak coinciding when the limbs were maximally flexed or extended (i.e., at the midpoint of the Move phase, Figure 7).

The feature-learnability temporal profiles of individual features (Figures 7B–E) provide insight into how significantly different a feature is for the different stimuli (and locations) throughout the progression of the mechanical stimuli (Loutit et al., 2019). HF features capture information about single or multiunit spiking activity, for example, HF integral, HF spike count and HF spike amplitude provide information about signal energy, spiking frequency and amplitude respectively; the latter may relate to the anatomical location relative to the recording site. As LF features capture measures of spike bursts, greater feature-learnability of LF features over the time course of the stimuli indicates the phases of the stimulus where the bursting activity feature is significantly different across the stimuli. For example, LF PSD signal features below 200 Hz do not differ greatly across the proprioceptive-evoked regions of the stimuli, whereas their differences are greater and therefore contribute more toward decoding tactile events. The data suggests that most of the features were able to extract some level of information from single-/multi-unit spikes as well as burst behavior of neural activity during the initial contact of the tactile stimuli, which might reflect the afferent volley of sensory information arriving at the DCN.

Feature-learnability remained above chance levels 1 s before stimuli were presented (Figure 7) for the HF spike amplitude and all HF PSD features. This indicates that during the rest periods, some information specific to the limb and/or the stimulus being presented was encoded in these DCN signal features. That the more successful predictions were not consistent across all animals (Figure 8), suggests that the prediction is not based on some residual activity from particular stimuli, but rather, multiple features in the same animal contributed information that was unique to that animal for a combination of stimulus and location. Furthermore, it was unlikely that above-chance level classification arose from some learned effect during repetitive stimulus presentations because feature-learnability was not different when comparing the first and last third of trials, although, it should be noted that training/testing data sets were also reduced to one third, making accurate classification more challenging.

What then could account for this predictive capacity before the stimuli? One possibility is that the limb undergoing repetitive stimulation experienced a different afferent activation status during the stimulus-off period compared to the remaining limbs. This could have arisen from continued reactivation of slowly adapting afferent activity during the off-stimulus period because of the continued interruption of the repetitive stimulation cycles, whereas the slowly adapting activity from the non-stimulated limbs would have greater opportunity to cease firing. An alternative explanation could be from activity resulting from long-lasting membrane potential depolarizations seen in DCN cells in response to sensory stimulation (Canedo et al., 1998), or rhythmic activity that outlasts stimulation periods (Nuñez and Buño, 1999). We have not found evidence from other studies finding similarly high decoding accuracy in stimulus off periods. Nevertheless, this phenomenon demonstrates remarkable sensitivity of the HF PSD and HF spike amplitude features for capturing status differences in afferent populations.