The dataset was obtained from the barrel cortex of four head-fixed water-restricted mice while performing a go/no-go whisker stimulus detection task. Mice were trained to respond to whisker stimulation by licking at a water spout (Figure 1A). During experimental sessions, a silicon probe with two shanks, 32 electrodes each, was used to obtain electrophysiological recordings from cortical layers (L) 2/3, 4, 5, and 6. Based on a current-source-density (CSD) analysis, we identified different cortical layers and assigned all electrodes to their corresponding layer, selecting the electrode that was most central to derive a layer-specific signal from the MEA recordings (see Methods). Using the signals from different layers allowed us to investigate which layer provides the most important information for stimulus detection and whether integrating the signals from all layers leads to more accurate detection.

In total, 19 experimental sessions were performed with 300 trials each. Each trial of the behavioral task consisted of different epochs (Figure 1B). Before whisker stimulation, the mice had to refrain from licking for at least 400 ms (silence period). Since stimuli were only presented when animals refrained from licking for at least 0.4–3 s (varying randomly from trial to trial), we were able to use the 400 ms before stimulus onset, called pre-stimulus, to characterize spontaneous activity in the barrel cortex. The stimulus itself consisted of a single-whisker 60-Hz sinusoidal vibration of 100 ms duration, delivered in the rostrocaudal direction to the C1 whisker. We used this peri-stimulus time window to analyze stimulus-evoked activity. If the mouse responded within 500 ms after stimulus onset, it was rewarded with a water droplet. Each recording session included 50 catch-trials without stimulation (0% intensity to assess spontaneous licking) and stimulus trials with different intensities ranging from 20 to 100% of maximum amplitude (corresponding to amplitudes between 93 and 463 μm and angular velocities ranging from 324 to 1,600 °/s at 2 mm from the whisker base, respectively; see Vandevelde et al., 2023, for more details). Stimuli in this range have been found to induce weak to moderately strong activity in barrel cortex and therefore span the range of barely to easily detectable stimuli (Yang et al., 2016; Vandevelde et al., 2023). Each stimulus was presented 50 times, and the stimulus sequence was randomized. The use of different intensities allowed us to investigate the response rate of mice during the behavioral experiment, as well as investigate the shape and strength of the evoked response in the barrel cortex. Using a randomized sequence was important to prevent mice from learning any patterns in whisker stimulation, which could potentially bias the response. The response rate increased monotonically with increasing stimulus intensity (Figure 1C), which is expected considering similar experiments (Stüttgen et al., 2006; Stüttgen and Schwarz, 2008; Stüttgen and Schwarz, 2010).

For more details about the behavioral experiment, we refer to our previous publication (Vandevelde et al., 2023).

As we integrated data from different recording sessions, we could not use the spiking activity as each session provides a different set of single- and multi-units, which precludes the assembly of a common feature vector to describe trials across sessions, a prerequisite for the ML algorithms. Besides this, it has been demonstrated on a single-session dataset that (i) LFP features are as informative as features derived from multi-unit activity and (ii) it is not necessary to derive more sophisticated features from the LFP, such as the response peak amplitude or response onset latency, because using the raw signal provides sufficient information for the ML methods used (Petschenig et al., 2022). In general, the LFP recording is more stable in comparison to action potentials, and therefore, it is recommended to use the LFP in the context of neuroprosthetics (Andersen et al., 2004; Markowitz et al., 2011). Following these previous studies, we use the raw LFP traces at 1,000 Hz resolution, low-pass filtered at 150 Hz but without further feature engineering, referring to them as RAW features, hereafter. Additionally, we investigate FFT features, consisting of the frequency components we obtain by applying the Fast Fourier Transform to the filtered signal (Brigham, 1988). Including the FFT features was inspired by Sederberg et al. (2019), who demonstrated that field potentials in the barrel cortex can be informative with respect to how whisker stimuli are perceived in the barrel cortex of awake mice.

We chose to use the LFP for two purposes: first, stimulus detection (SD), i.e., to decode whether a stimulus was presented to the animal, and second, response prediction (RP), i.e., to predict whether the animal was responding to the whisker stimulus. To that end, we benchmarked six ML classification algorithms, also called classifiers, known for good performance in various applications (see Table 1 for an overview of the selected classifiers). The Decision Tree (DT) achieves short runtimes due to its simplicity (Loh, 2011). Compared to a single decision tree, Random Forest (RF), which makes use of an ensemble of decision trees, and XGBoost (XGB), which uses stacked decision trees, are known to provide more accurate models while maintaining low runtime requirements (Breiman, 2001; Chen and Guestrin, 2016). Similar characteristics are associated with generalized linear models (GLM) that are, once trained, very fast during inference as they rely on a linear combination of weighted variables, so the logistic regression, a GLM with underlying logit function, has been included in the benchmark (Kleinbaum and Klein, 2010). Interestingly, comparable studies investigating brain activity classification did not explore algorithms developed for time series classification. We, therefore, included ROCKET (RCKT), a classification algorithm dedicated to time series input (Dempster et al., 2020; Pantiskas et al., 2022). As mentioned above, spiking neural networks are promising due to their ability to run efficiently on low-powered neuromorphic hardware, which motivated the inclusion of the Liquid State Machine (LSM). More complex artificial neural network architectures, such as long short-term memory (LSTM), have not been explored due to their high computational requirements and because they did not provide better models in comparable studies (Petschenig et al., 2022).

All algorithms were compared based on different classification scenarios. The term classification scenario is here used to describe a specific combination of which classification task (SD or RP) is considered and which features are used, further specifying the scenario by type (RAW, FFT), cortical layer, and time window from which the recording was derived (Figure 2A). For each algorithm, hyperparameter tuning was applied by evaluating a grid of parameter settings, also called grid search, using 3-fold cross-validation on the training data (see Supplementary Table T1 for details). Furthermore, the model evaluation was carefully done, splitting up experimental trials based on individual mice to ensure cross-individual validation in a way that trials from one individual are separated by the training set, used for hyperparameter tuning and training the classification model and testing set, used to evaluate the model performance (Figure 2B). Training–testing splits have been repeated four times, using four individuals and for each split, 10 bootstraps were drawn to increase the reliability of the aggregated accuracy measures. For each bootstrapping, the size of the training and testing set differs because the number of experimental sessions available differs between mice (see Supplementary Table T2 for details about the set sizes).

In all classification scenarios, we used perfectly balanced datasets, meaning each subset contained the same number of positive trials (stimulus trial or response trial) and negative trials (catch trial or no-response). This was achieved by down-sampling trials from the overrepresented, and due to the expected 50% positives in the data, we can assume that the baseline accuracy achieved by random predictions converges to 50%. Note that including some scenarios was mainly motivated by creating a baseline scenario in which a low accuracy (close to 50%) is expected. For instance, we would not expect pre-stimulus LFPs to encode information useful for stimulus detection. Hence, an accuracy of 50% would be expected for this scenario from any classifier. For response prediction, the expectation can be slightly different. Even though the barrel cortex is not involved in decision-making, its spontaneous activity might encode different states related to attention and excitement, which is relevant to the behavioral experiment (Cowley et al., 2020). Hence, it might be possible for classifiers to leverage weak patterns predictive for the response, even from pre-stimulus (spontaneous) activity.

The first analysis addresses the question of the extent to which ML classification algorithms, also called classifiers, can leverage information from LFP features to train accurate classification models for response prediction (RP) or stimulus detection (SD) and which classifier trains the most accurate models.

Figure 3 provides an overview of the accuracy obtained for all classifiers (innermost y-axis) applied to the classification scenarios, specified by the two classification tasks (SD, RP), the time window for the recording (PRE, PREI, FULL), feature type (RAW, FFT), cortical layers (innermost x-axis), and the stimulus intensities used for trial selection (see subtitles). Brighter colors represent higher accuracy observed for SD in general. Using the PERI and FULL time window results in higher accuracy for both tasks, SD and PR. Regarding the feature types, the RAW features are more informative than the FFT, shown by higher accuracy for both SD and RP. The lower half of the heatmap shows accuracy measures obtained for SD, for which a higher level of accuracy can be observed in general. The accuracy for SD reached up to 90% when trials with 0% (catch-trials) and maximum stimulus intensity of 100% were contrasted.

Features derived from L4 are more informative than those derived from other layers. Considering a feature vector combining all layers (ALL on the x-axis) slightly improves the accuracy over the performance of L4 alone. As expected, the peri-stimulus activity (PERI), describing the stimulus-evoked response in the LFP, is more informative for both classification tasks in comparison to the pre-stimulus activity. Classification models trained on pre-stimulus activity features are unable to perform accurately for SD, as expected. For RP, however, the accuracy slightly increased over the 50% baseline when RF is used on pre-stimulus FFT features (Figure 3).

Regarding the classifiers, RF performs well in scenarios exploring FFT features. For the RAW features that are more informative for the SD task, we observed that the GLM and LSM achieved the highest performance, with RCKT slightly lower. We do not see an improvement in model performance for SD when integrating the FULL time window over using peri-stimulus alone (Figure 3).

However, from the median accuracies discussed in this section, it is not clear if it is worth integrating all layers, expected to increase model training time, over using solely L4. Moreover, a runtime analysis is required for the best algorithms, considering the practical application of these models. Therefore, a more detailed comparison of the most accurate algorithms follows, focusing on SD based on the peri-stimulus RAW features.

For this comparison of the best-performing classifiers, the SD classification task has been selected with 0 and 100% intensity trials only, because the highest classification accuracy was achieved in this scenario. Considering GLM, RCKT, and LSM as the top-performing classifiers in terms of accuracy, the LSM is the strongest competitor based on the runtime comparison (Figure 4A). RCKT was faster than the GLM for model training but much slower during inference. In general, the runtime for model training and inference was higher when all layers were used, which is expected as algorithms are processing four times more ML features than using only a single layer. Focusing on LSM and GLM as the fastest algorithms, we investigated the difference in accuracy for the SD task using all layers and only L4 (Figure 4B). In this comparison, the LSM achieved the lowest performance when all layers were integrated, which suggests that the GLM is more efficient in dealing with high-dimensional data as it achieved the highest median accuracy overall when it integrated all layers. However, as the accuracy achieved by the GLM with all layers is not significantly higher than that of the LSM with only L4, we conclude that the LSM is the best choice for stimulus detection, especially due to its superior performance in the runtime comparison with the lowest model training and inference runtimes using only L4. Furthermore, it achieved slightly more robust accuracy distributions across all individuals than the GLM, for which more outliers and lower median accuracy measures were observed overall (Figure 4C). To assess whether the accuracy scores of the GLM and LSM (Figure 4B) differed significantly, the Wilcoxon signed-rank test was applied on 40 pairs (four mice times ten bootstraps). A pair of values consists of the accuracy values from the two classifiers, obtained on the same mouse-specific testing set and the same bootstrap sample, using random seeds to ensure bootstrap samples are equal.

For practical applications, it is important to know how accurate the detection is for stimuli of different intensities (Petschenig et al., 2022; Wang et al., 2018b). Furthermore, it is crucial to investigate how long the LFP recordings need to be to maintain accuracy at the highest level possible. The shorter this recording time is, the earlier the model can be applied after stimulus onset. To address these two aspects, an additional SD experiment was performed using subsets of the data with catch-trials and stimulus trials of intensities from 20 to 100%. Note that this includes the results from the benchmark obtained for 100% stimulus intensity and the complete 100 ms of peri-stimulus activity. These SD experiments were further specified by different lengths of the RAW LFP signal as input for the classification model, using L4 as the most informative layer. All additional classification experiments were done for the LSM using the validation strategy applied during the initial ML benchmark analysis (Figure 2B).

We observed that LSM’s detection accuracy is lower for stimuli of low intensities (Figure 5). In comparison to the accuracy achieved on the complete (100 ms) peri-stimulus LFP, approximately 35–40 ms is required to achieve a high level of accuracy for the higher intensities greater than or equal to 60%. For 40% intensity, even though the overall accuracy was much lower than for higher intensities, the LFP can be as short as 15 ms to achieve high accuracy for this case. Note that the RAW features were used for these tests, representing the LFP after low-pass and band-stop filtering, which further delays this process by a total of 0.2 ms (measured on a common CPU, expected to be slightly higher on neuromorphic hardware).

To mitigate potential impacts of class imbalance (over- or underrepresentation of the positive class), datasets have been downsampled to train models on data with 50% positive (stimulus present) and 50% negative (stimulus absent) samples. For deployment, however, it is important to investigate how algorithms perform on imbalanced data. Using L4 as the most informative layer and RAW as the most informative feature type for stimulus detection, we benchmarked the GLM and LSM again and added comparisons on data without downsampling. We used the dataset with all stimulus intensities, resulting in 83% stimulus trials, and the dataset containing only 0% stimulus intensity (catch) trials and 100% stimulus intensity trials, for which downsampling had almost no impact, as the fraction of positives changed from 50 to 49.3%. Due to the class imbalance in the data, we added ROC and Precision-Recall curves to evaluate model performance (Saito and Rehmsmeier, 2015). Additionally, we used the Brier loss score to evaluate model calibration and studied the impact of Platt Calibration on the model probabilities from GLM and LSM (Platt, 1999; Niculescu-Mizil and Caruana, 2005).

The results from this analysis (see Supplementary Figure S1) show that (i) the GLM is well-calibrated and does not benefit from Platt Calibration, (ii) the LSM needs Platt Calibration to improve its calibration performance, especially for imbalanced data, and (iii) the overall performance of GLM and LSM is very similar, both achieving area under ROC and Precision-Recall curve clearly better than a random guess or naïve model.

Note that testing data is not used for Platt Calibration. A default Logistic Regression (scikit-learn) is trained as a calibration model on the classification probabilities and class labels of a 33% hold-out dataset from the training set. Consequently, the model evaluated on the mouse-specific testing set has been fitted on only 67% of the training set when using Platt Calibration. The fitted calibration model is applied to transform the classification probabilities from the testing set before computing the evaluation metrics.

Besides investigating neural activity with respect to how whisker stimuli are processed in the barrel cortex, the dataset allowed us to search for patterns predictive of the behavioral response. This is a critical aspect as we record from a somatosensory brain region believed to be involved in the process of evaluating an incoming stimulus and provide the basis for sensory-driven decision-making (Staiger and Petersen, 2021). Considering the response prediction analysis as part of the ML benchmark, we observe a moderate accuracy for response prediction models trained on the RAW features from the peri-stimulus or full signal (Figure 3). However, an additional analysis based on the predictions of the corresponding models strongly suggests that the response prediction models actually perform stimulus detection, which results in quite accurate models due to the strong correlation between the presence and absence of stimulus and the occurrence of behavioral response (see Supplementary Figure S3).

Interestingly, our benchmark revealed that patterns exist in the harmonic oscillations derived from spontaneous (pre-stimulus) activity that are slightly predictive toward the response of mice (see RP with pre-stimulus FFT features, all intensities, Figure 3). For the FFT features derived from L2/3, L4 or L5 we observed accuracy slightly, but significantly better than the random baseline when all stimulus intensities are considered in the data (p-value from Fischer’s test < 1e⁻²⁰ in all three cases). Even though this difference in accuracy is small, these results indicate that neuronal activity in the barrel cortex is also related to attention or anticipation of the expected whisker stimulus (Lee and Dan, 2012; Stüttgen and Schwarz, 2018). Such weak patterns in the spontaneous activity could reflect brain states in which the mouse is more attentive to perceive the stimulus, which impacts the response, especially when not only the full-intensity trials but also those with a small or intermediate intensity are analyzed. While these are indications relevant to the role of the barrel cortex in the behavioral response of mice, future work is required to further investigate the FFT features, which is out of scope for this study and might require further behavioral experiments.

Previous in vivo recordings in the rat barrel cortex have identified inhibitory interneurons in L4 as the units carrying the highest amount of sensory stimulus-related information (Reyes-Puerta et al., 2015; Reyes-Puerta et al., 2015). Thus, neuronal activity in L4 represents a good target to obtain information on the spatio-temporal properties of the sensory input. Measuring the activity from L4 is not only possible but also very reliable with invasive technologies as used in our experimental model. However, in a practical application, it will be more beneficial to use non-invasive technologies to avoid an intervention in the patient’s brain. Resolving activity in a specific portion of tissue without having direct access could be realized by laminar inference based on MEG recordings, although MEG is not applicable in neuroprosthetics (Bonaiuto et al., 2018).

Our study confirms that ML models can be highly accurate on whisker stimulus classification, even for LFP signals recorded from awake mice. However, the datasets used for such explorations, including the one we used, are based on precisely controlled single-whisker stimuli, certainly different from how whiskers bend and move when animals are actively palpating during unrestrained exploration (Petschenig et al., 2022; Wang et al., 2018a; Wang et al., 2018b). Therefore, for future work, it would be relevant to apply and test the suggested algorithms on recordings from microchips implanted in freely moving animals. The hardware required to do this could use wireless and battery-free devices (Martinez et al., 2018; Won et al., 2023). Recordings from such experiments might result in LFP traces that are less distinguishable, as we would expect that whisker deflections from freely moving animals are less distinct.

Additionally, our benchmark for imbalanced data covers only the case of overrepresentation of the positive class, which is the opposite of a real-life scenario in which stimulus events are less frequent compared to non-stimulus recordings. Our dataset provided well-defined trials for events with an induced stimulus, with five different intensities over catch trials without a stimulus. Thus, using all trials of the dataset without downsampling results in overrepresentation of the positive class. Nevertheless, this analysis demonstrates that both algorithms perform well on imbalanced data in general and reveals that additional calibration of classification probabilities is required when using the LSM, particularly on imbalanced data. Interestingly, on the imbalanced data, the GLM and LSM achieved accuracy values close to those expected from a naïve model that always returns the positive class (stimulus present). However, both algorithms yielded high area-under-the-curve scores when using the ROC and Precision-Recall curves as threshold-independent evaluation metrics. We therefore conclude that, in imbalanced deployment settings, both algorithms require decision threshold optimization, which is left for future work.

To realize the closed-loop system, the output of the ML classifier could be used to control an electrode located in another brain region to induce intracortical microstimulation that can convey artificial signals that may be as complex as the somatic perception of tactile properties from different objects (Jackson et al., 2006; O’Doherty et al., 2011). However, a complete experimental setup would be required to verify whether the whole process of reading in the MEA recording, applying the ML classifier, and creating a feedback signal is fast enough to restore brain functionality, establishing timely signal processing with other cortical areas.

Another important aspect is the robustness of the findings concluded from our ML-based analyses. The sample size is a key factor in statistical analyses, and machine learning models tend to be more meaningful and generalizable when trained on large datasets (Carlini and Wagner, 2017; Rice et al., 2020). Even though electrophysiological recordings are challenging from different perspectives, we could explore a comparatively large dataset using >3,400 trials for six stimulus intensities. The final subsets used for the ML-based analysis were restricted by a selection of trials for the different intensities, resulting in a sample size of at least 600 trials (see Supplementary Table T2), which is approximately three times higher than in comparable studies using only one recording session, e.g., (Petschenig et al., 2022; Wang et al., 2018a; Wang et al., 2018b).

In summary, based on a large electrophysiological dataset enabling a robust model evaluation, we demonstrated that SNNs are accurate and efficient in interpreting neural activity toward neuroprosthetics applications related to the barrel cortex. Using the LSM classifier, model training could be realized on a small microchip to be implemented on the prosthesis, which can be beneficial to retrain a model for adjusting to small drifts of the implanted electrode, for instance. While our tests were done on an infrastructure with common CPUs with the main perspective to compare the model training time between the different classifiers, it has been shown recently that LSMs can efficiently be implemented and executed on microchips, as demonstrated on a DYNAP-SE neuromorphic processor (Petschenig et al., 2022). Another opportunity is to train a model on an external machine, transferring the model and only applying it on the prosthesis (Wang et al., 2018b). If future technologies enabled a more stable implantation of the electrodes, retraining on the prosthesis might not be necessary. Our generalizability analysis confirms that using a pre-trained model is possible even if the underlying data was derived from an individual different from the patient who needs the prosthesis.

The data investigated for this study were obtained from experiments with head-fixed mice performing a go/no-go whisker stimulus detection task (Vandevelde et al., 2023). Mice were water-restricted and learned that they were rewarded by a drop of water if licking within 500 ms after the onset of a 100-ms sinusoidal vibration of a single whisker. Before mice were trained to perform this task, they were habituated to the head-fixed setting, which was important to enable the recording of neural activity in the cortical column in the barrel cortex associated with the stimulated whisker. This was done with multi-electrode-array (MEA) silicon probes through all cortical layers using 2-shank-64-channel probes with a distance between the shanks of either 200 or 250 μm. (Neuro Nexus, Ann Arbor, United States, or Cambridge Neurotech, Cambridge, United Kingdom, respectively). Each shank has 32 electrodes, and the distance between electrodes is 25 μm, which enables the observation of neural activity over an extent of ~0.8 mm.

For more details about the mouse line, behavioral task training, habituation, surgery, MEA recordings, and the layer-electrode association using current source density (CSD) plots, we refer to Vandevelde et al. (2023). After inspecting the CSD plots, we used one electrode per layer to derive layer-specific LFPs from the MEA recordings. The location and extent of layer 4 could be easily identified by a prominent short-latency current sink in the CSD, computed from the average LFP response to repeated high-intensity C1 whisker deflections, and therefore could be separated from layers 2/3 above and layer 5 below. The border between layers 5 and 6 was estimated based on recording depth (for more details, see Vandevelde et al., 2023, as well as Reyes-Puerta et al., 2015; Yang et al., 2017, 2018). Psychometric curves (Figure 1C) differ slightly from those shown in our previous study (Vandevelde et al., 2023) because we excluded trials with a lick response within the first 100 ms after stimulus onset as the muscle movement related to licking caused artifacts in the recordings.

To exclude LFP traces contaminated by artifacts, we applied trial filters for the full range from −420 ms to +120 ms according to the following criteria. Trials in which animals were presumably moving and/or licking randomly had large voltage fluctuations and sometimes burst-like oscillations in activity. Trials in which the signal reached saturation were identified by applying a cutoff of ±2 mV at any time point within the time range specified above. Bursts were identified by moving time windows of 20 ms in which all electrodes strongly correlated, computed by Pearson’s correlation (Schober et al., 2018). Electric shorts were partly captured by the MEA electrodes, resulting in a continuous 50 Hz sinusoidal signal, and trials with a strong 50 Hz frequency power were removed. Trials for which a lick was detected within the above-mentioned time window were also excluded.

The RAW features describe the LFP trace at 1,000 Hz resolution, resulting in 400 features describing, for instance, the pre-stimulus activity of 400 ms from −400 ms to 0 ms (stimulus onset). Hence, the peri-stimulus activity from 0 to +100 ms is described by a feature vector of length 100, while 500 features represent the full time window from −400 ms to +100 ms. Using causal Butterworth filtering, the LFP traces were low-pass filtered at 150 Hz, and the 47–53 Hz band was removed, using a first-order low-pass and a sixth-order band-stop filter (Selesnick and Burrus, 2002). Using the same time-windows, the Fast-Fourier-Transform (FFT) was applied to derive the strength of the frequency components (Brigham, 1988). According to the length of the signal derived from the investigated time windows, the resolution is 2.5 Hz and 10 Hz going up to 147.5 Hz and 140 Hz for the pre- and peri-stimulus activity, respectively. The first FFT feature (0 Hz) represents the constant, which is simply the average over the signal. Eventually, 60 pre-stimulus and 15 peri-stimulus features describe the FFT feature vectors, resulting in 75 features for the full range. Due to the strong differences between ongoing and evoked activity, the FFT was applied to those time windows separately, and the resulting FFT features were merged afterward in order to create the full (pre- and peri-stimulus) FFT feature vector.

We used MOGON II, the high-performance computing system from the University of Mainz. The whole analysis was implemented in Python. We ran DT, RF, and GLM from the scikit-learn package (Pedregosa et al., 2011) and XGB using a separate library called xgboost (Chen and Guestrin, 2016). The algorithm ROCKET (RCKT) was integrated via the sktime package (Löning et al., 2019), implementing its univariate and multivariate variants (Dempster et al., 2020; Pantiskas et al., 2022). The LSM algorithm was used from the GitHub repository https://github.com/IGITUGraz/LSM, published by Maass et al. (2002). The CPU specifications on the SMP compute nodes on MOGON are Intel® Xeon® E5 v4, 2.2GHz.