An artificial neural network for automated behavioral state classification in rats

Animal care and use

All sample data were obtained from previous experiments conducted at Middlebury College with all methods performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and approved by Middlebury College’s Institution Animal Care and Use Committee (approved research protocol #316-17). Male, Sprague-Dawley rats (3–4 months old, Charles River, Wilmington, MA, USA) were obtained and pair-housed upon arrival in clear, plastic rodent caging (16 in. × 7.5 in. × 8 in.; Teklad TEK-Fresh bedding, Envigo, Indianapolis, IN, USA). Food (Teklad 2020X) and water were provided ad libitum along with in-cage enrichment objects (wooden blocks, pvc tubing, paper towels). After an initial acclimation period (>1 week from initial arrival to the vivarium), each rat underwent stereotactic surgery. Here, rats were anesthetized via isoflurane (3.5% induction, 2–3% maintenance) and given pre-operative analgesic (Meloxicam; 2 mg/kg; MWI, Boise, ID, USA) and antibiotic (Penicillin, 100,000 units/kg) treatments. During surgery, electrodes were implanted for EEG/LFP and nuchal EMG recordings (see sample data set below for additional recording characteristics) and were affixed to the skull using dental acrylic (Lang Dental; Wheeling, IL, USA). A postoperative analgesic (Meloxicam, 2 mg/kg) was administered between 12–24 h after surgery completion. For the duration of the experiment, rats were then single-housed to minimize potential damage to implanted electrodes and the headstage preamplifier (100× amplification, Pinnacle Technologies, Lawrence, KS, USA). Throughout the experiment, each rat was monitored daily for overt signs of pain/distress including immobility, poor grooming, weight loss, porphyrin staining, postural abnormalities, lack of food/water consumption, and signs of infection around the surgical site. Observation of these signs leads to direct consultations with the attending veterinarian and vivarium staff to determine the appropriate course of action (e.g. treatment or euthanasia). In the present study, no rats were euthanized prior to the conclusion of the experiment. In accordance with the 2013 AVMA Guidelines of Euthanasia, all rats were euthanized via CO₂ exposure at the end of the experiment.

Sample data set

Uninterrupted baseline recording days (N = 50 from 11 total rats) were selected from male, Sprague-Dawley rats (3–4 months old). Care was taken to generate a sample of recording days and rats that collectively comprise common recording characteristics across rodent sleep research. Specifically, within each rat, two electrophysiological signals to measure brain activity and one nuchal electromyogram (EMG) were recorded (sampling rate = 250 Hz; 8401 DACS; Pinnacle Technologies, Lawrence, KS, USA). Rats used in our sample data set had a wide range of recording locations (e.g. prefrontal cortex, motor cortex, parietal cortex) and diverse recording modalities (2 EEGs, 1 EEG/1LFP, or 2 LFPs). Likewise, these data differed in overall signal quality (e.g. 4.83 ± 0.56% mean daily artifact prevalence; range: 0.15% to 13.48%). Consequently, this data set represents typical electrophysiological recordings of rodent sleep. Of note, two rats were ultimately excluded from analyses because of very poor EMG signal quality resulting in a final sample data set of 40 baseline days recorded across nine rats. Details for how EMG signal quality was assessed are presented below as part of our description of our artificial neural network for automated sleep classification.

Trained undergraduate sleep researchers visually scored behavioral state offline in 4 s time intervals for each baseline recording day (following conventional sleep research terminology, 4 s time intervals will be referred to as epochs throughout the remainder of the manuscript). Epochs were classified as waking when containing low-voltage, high-frequency EEG/LFP activity and elevated EMG. High-voltage, low-frequency EEG/LFP activity and an absence of EMG activity was characteristic of NREM sleep, while REM sleep epochs were classified when low-voltage, high-frequency EEG/LFP activity and an absence of EMG activity were observed.

Artificial neural network

To develop an effective research tool with minimal entry barriers for use, we sought to design a freely available, open-source implementation to automatically classify behavioral state using electrophysiological features. With this goal in mind, we chose to implement our algorithm using R (R Core Team, 2020), a freely-available statistical computing platform that has been widely adopted across scientific disciplines (e.g. 58% of ecological papers published in 2017 used R as their primary analytic tool Lai et al., 2019). Instructions for full implementation of our algorithm and source code can be found at https://github.com/jellen44/AutomaticSleepScoringTool.

Within R, we implemented an artificial neural network (ANN) with four, fully connected sequential layers (256, 128, 32 and 3 nodes) to automatically classify behavioral state in freely behaving rodents. Our ANN was implemented using the Keras package in R with a Tensorflow backend (Abadi et al., 2016). In each of the hidden layers, a weighted sum of inputs is computed, followed by an activation function rectified linear unit activation (ReLU) that introduces nonlinearity and regulates neuronal activation. These weights are adjusted during the training process to improve classification performance. The output layer of this network takes input from the last hidden layer and outputs a vector containing three values, each representing the probability that the epoch is one of the three sleep states. We accomplished this with the SoftMax function, which is commonly used in multiclass classification problems because it takes a K-dimensional vector (the final layer in our network) as an input and uses it to estimate a range of probabilities over a given number of classes (Duan et al., 2003).

Analyses and statistical approaches

All analyses of model performance were conducted in Mathworks Matlab (Natick, MA, USA) with additional statistical analyses (correlation, repeated-measures ANOVA, t-tests) conducted in SPSS (IBM, Armonk, NY, USA). All data are presented as mean ± standard error.

Performance of the automated sleep scoring algorithm was assessed by comparing its predictions of behavioral state for each epoch with corresponding epochs classified manually. Key metrics calculated include overall accuracy (% of epochs that were scored identically by automated and manual approaches), state-dependent sensitivity (TP/TP + FN), and state-dependent specificity (TN/TN + FP); where TP, true positive; FN, false negative; TN, true negative, and FP, false positive. To determine the extent to which training the ANN using data from one baseline day could generalize to (a) other baseline days from the same rat and (b) other baseline days from different rats, we likewise calculated accuracy, sensitivity, and specificity of the algorithm under these conditions. For within rat generalization, ANNs that were separately trained from each baseline day in the data set were used to predict behavioral state for all other baseline days recorded from that rat. For between rat generalization, ANNs that were separately trained from each baseline day in the data set were used to each predict behavioral state for five baseline days randomly selected from other rats.

Lastly, we calculated the effects of using ANN classification on common sleep parameters. Power spectra were calculated for each 4 s epoch using Welch’s method (Hamming window) and averaged across all epochs of the same behavioral state. To account for potential differences in overall signal strength, spectra were normalized to mean broadband power (i.e. from 0.5–125 Hz). To calculate slow wave activity (SWA), we calculated band limited power (0.5–4 Hz) for each 4 s epoch of NREM sleep and averaged these values across each hour of the light period. Hourly SWA values were then normalized to the mean SWA across the light period.

Electrophysiological characteristics associated with behavioral state are useful features for automated state classification

Behavioral state scoring, be it manual or algorithmic, relies upon the distinct patterns of neuronal activity characteristic of each behavioral state. Waking epochs are typified by low-voltage, high-frequency EEG activity and variable EMG activity. By contrast, as a function of the large slow waves that predominate NREM sleep, EEG activity during NREM epochs is characteristically high-voltage and low-frequency. REM sleep meanwhile, exhibits EEG activity that is similar to that of waking, but with the complete absence of EMG activity due to REM-associated motor atonia. While these electrophysiological differences can typically be discerned visually in the time-domain (see Fig. 1A for example tracings), quantification of state-dependent activity is usually performed through the extraction of key frequency-domain features. As evident in Fig. 1B, EEG band-limited power (BLP) likewise exhibits pronounced state dependency. For example, delta BLP (0.5–4 Hz; i.e. low-frequency) is clearly elevated during NREM sleep while higher frequency BLP (e.g. theta, gamma) is greater during waking and REM sleep than NREM sleep.

State scoring algorithms consistently rely upon a collection of these electrophysiological features to classify behavioral state (e.g. Allocca et al., 2019; Gao, Turek & Vitaterna, 2016; Kreuzer et al., 2015; Yan et al., 2017). For example, by simply using three features (Delta BLP, Gamma BLP and EMG-RMS), we can reliably segregate the majority of epochs from each behavioral state (Fig. 2A). However, as this simple example also shows, many boundary cases remain difficult to classify because a substantial proportion of epochs from each state still overlap within the feature space. To enhance the utility of using electrophysiological features for behavioral state classification, we implemented an artificial neural network (ANN) to predict behavioral state from EEG/EMG features (Fig. 2B). During model training, computational weights are continually altered to produce activation maps that accurately predict behavioral state from electrophysiological feature input (see “Methods” and Fig. 2C). Once trained, the ANN can uniquely and nonlinearly combine electrophysiological features to generate distinct activation maps characteristic of each behavioral state. Figure 2D depicts activation maps across all four layers of our ANN in response to feature input from manually scored waking, NREM, and REM epochs. With each progressive layer of the ANN, activation maps converge for epochs derived from the same behavioral state and diverge for those derived from differing states. In doing so, the ANN efficiently and effectively uses electrophysiological features to classify behavioral state without researchers needing to know how to precisely combine feature input, a priori.

An artificial neural network effectively classifies behavioral state

To assess the performance of the ANN detailed above, we compared its behavioral state predictions with manual scoring of the same electrophysiological data by trained undergraduate researchers. Figures 3A–3B depict behavioral state hypnograms for a single day using manual and algorithmic scoring, respectively. Here we observe that each scoring method produces a qualitatively similar classification of behavioral state across the 24 h day. Quantification of the model’s performance for this day reveals that the ANN predicted behavioral state with high accuracy, sensitivity, and specificity (Fig. 3C). Across 40 recording days from nine rats (Fig. 3D), the ANN predicted behavioral state with an overall accuracy of 90.87 ± 0.31%, high sensitivity (Wake: 89.82 ± 0.50%, NREM sleep: 92.81 ± 0.58% and REM sleep: 86.28 ± 0.99%), and high specificity (Wake: 95.18 ± 0.34%, NREM sleep: 92.71 ± 0.41% and REM sleep: 97.18 ± 0.19%). Thus, our overall model performance is comparable with that of other high-performing sleep scoring algorithms despite (a) the small amount of training data needed (560, 4 s epochs; 2.6% of day’s total) and (b) in contrast to some algorithms, the ANN predicts behavioral state for all test epochs including those manually-scored as artifact (see Table 1). Moreover, the ANN accomplishes this high performance with rapid computational time (average total run time including fast Fourier transform, feature selection, model training, and prediction output of 223.0 ± 1.1 s when run on a 2.4 GHz Quad-Core Intel Core i5 processor with 8GB 2133 MHz LPDDR3 RAM).

Although the epochs used to train the ANN models characterized above were randomly selected as described in the methods, the overall performance of these models could nevertheless be dependent upon the specific training epochs that happened to be selected. To address this possibility, we repeated the analyses described above five times for each recording, each time using a different random selection of training epochs. Random differences in training epochs selected produced minimal alterations to model performance; between models trained on different random epochs, low average standard errors were observed for overall accuracy (0.25 ± 0.02%), sensitivity (Wake: 0.55 ± 0.05%, NREM sleep: 0.53 ± 0.04% and REM sleep: 1.81 ± 0.16%), and specificity (Wake: 0.51 ± 0.05%, NREM sleep: 0.46 ± 0.04% and REM sleep: 0.24 ± 0.02%). Thus, model performance appears largely unaffected by which specific 2.6% of epochs researchers select for model training.

Although our model compares favorably with other automated scoring approaches, its intended purpose is to serve as a functional tool for sleep researchers. Therefore, more direct assessments of the functional consequences of its use may additionally serve as important metrics of performance. Consequently, across our entire data set, we compared common sleep parameters calculated using our ANN scoring with those calculated using manual scoring. Both scoring methods yield characteristic state-dependent power spectra (Fig. 4A); power spectra exhibit clear 1/f power law scaling (Bédard, Kröger & Destexhe, 2006), prominent peaks within theta frequencies during waking and REM sleep, and pronounced slow wave activity (SWA) during NREM sleep. Overall, we observed expected significant differences in power spectra as a function of (1) behavioral state (F(2,16) = 47.26, p =1.9 × 10⁻⁷) and (2) frequency (F(77,616) = 72.18, p = 2.32 × 10⁻²⁶⁰), but did not observe a significant effect of scoring method (F(1,8) = 0.07, p = 0.79). Furthermore, we observed a canonical homeostatic decline in SWA across the light period regardless of scoring method employed (Fig. 4B; effect of light period time: F(11,88) = 12.88, p = 3.51 × 10⁻¹⁴; effect of scoring method: F(1,8) = 3.63, p = 0.09). Thus, key electrophysiological characteristics of behavioral state and homeostatic sleep are not significantly affected by the use of this ANN for behavioral state classification. Lastly, we examined whether ANN classification altered total duration within each behavioral state (Fig. 4C) and average episode duration (Fig. 4D). During the dark period, the amount of time spent within each behavioral state significantly varied as expected (F(2,16) = 618.92, p = 7.03 × 10⁻¹⁶), and this effect was not significantly altered as a function of scoring method (F(2,16) = 1.97, p = 0.17). In the light period, however, observed significant differences in state duration (F(2,16) = 481.88, p = 5.06 × 10⁻¹⁵) were significantly affected by scoring method (F(2,16) = 15.83, p = 1.6 × 10⁻⁴). Post-hoc t-tests reveal that this effect appears largely driven by a significant increase in REM sleep scored by the ANN (24.32 ± 4.60% more light period REM) at the expense of a significant decrease in light period waking (−9.00 ± 1.58%). Mean episode durations across the entire day, however, were not significantly affected by scoring method (Fig. 4D; F(1,8) = 0.42, p = 0.54). Thus, the use of this ANN for behavioral state classification does not appear to significantly affect most common sleep parameters, with the exception of an increased prevalence of light-period REM duration that arises at the expense of waking duration.

Optional approaches enhance scoring accuracy or reduce scoring time

Ideally, automated sleep scoring algorithms afford significant time savings with minimal classification error. To this end, we have shown above that our ANN scoring approach classifies behavioral state with high accuracy, sensitivity, and specificity, while only requiring researchers to manually score 2.6% of the day. Indeed, increasing the size of the initial training set beyond 2.6% of the day did not appreciably increase model performance (Fig. S1). Certain research applications, however, may necessitate enhanced accuracy even at the expense of additional scoring time. Alternatively, analyses of very large data sets may benefit from reduced scoring time as long as classification performance is not severely diminished. Below we characterize two optional approaches available with our ANN to satisfy these contrasting needs.

In addition to producing a categorical behavioral state prediction for each epoch, our ANN generates an estimate of the certainty of each prediction. To explore how model certainty relates to model performance, we calculated sensitivity and specificity metrics for each behavioral state as a function of model certainty (Fig. 5A). Model sensitivities (Wake: r(238) = 0.80, p = 1.10 × 10⁻⁵⁴; NREM: r(238) = 0.90, p = 8.69 × 10⁻⁸⁸; REM: r(238) = 0.81, p = 4.32 × 10⁻⁵⁷) and specificities (Wake: r(238) = 0.75, p = 1.33 × 10⁻⁴⁴; NREM: r(238) = 0.64, p = 4.68 × 10⁻²⁹; REM: r(238) = 0.83, p = 2.75 × 10⁻⁶²) were all significantly correlated with model certainty. Moreover, these analyses reveal that the vast majority of epochs (88.03 ± 1.23%) were predicted with model certainty greater than 90%. These observations raise the possibility that manually rescoring the relatively small proportion of lower certainty epochs could greatly enhance classification performance. As evident in Fig. 5B, manually rescoring epochs associated with lower prediction certainty provides a focused approach that (1) increases sensitivity and specificity across each behavioral state and (2) still affords a significant reduction in researcher time as the ANN automated scores are used for the vast majority of epochs (see Table 1).

Lastly, although the base implementation of our algorithm only requires researchers to manually score 2.6% of the original file, these small manual scoring requirements may nevertheless present significant time constraints when analyzing very large datasets. Since electrophysiological features associated with each behavioral state are likely to share similar characteristics across recording days and subjects, our trained ANN may be generalizable. To test this possibility, we first trained the ANN on one recording day and then used the trained model to predict behavioral state across all other recording days from the same rat (Fig. 6A). Under these conditions, the ANN model maintained high levels of accuracy, sensitivity, and specificity, albeit with a small drop in performance relative to when tested on the same day as trained. Small, but statistically significant (paired t-tests), reductions in overall model accuracy (−1.74 ± 0.41%), NREM sensitivity (−1.95 ± 0.62%) and Waking/NREM specificity (−2.38 ± 0.77% and −1.24 ± 0.57%, respectively) were observed. Waking sensitivity and REM specificity were not significantly affected, while a more pronounced and statistically significant decrease in REM sensitivity was present (−6.42 ± 2.17%). Collectively, these results appear to indicate that the ANN can effectively classify multiple recording days from the same subject when trained using epochs from a single day. We then sought to determine whether an ANN trained on data from one rat could generalize to recordings from other rats (Fig. 6B). While the ANN still maintains moderate performance under these conditions for most parameters of interest (e.g. overall accuracy: 85.33 ± 0.40%), significant decreases in performance were observed relative to ANN models trained on the same day as tested in all performance metrics except REM specificity. REM sensitivity was particularly affected with this approach and significantly decreased by −27.50 ± 2.52%. Thus, the ANN appears limited when generalizing to recordings from different rats.