Cortical Pyramidal and Parvalbumin Cells Exhibit Distinct Spatiotemporal Extracellular Electric Potentials

Experimental animals

Seven PV::ChR2 male mice, generated by crossing homozygous male Ai32 mice (catalog #012569, The Jackson Laboratory) with homozygous female PV-Cre mice (catalog #008069, The Jackson Laboratory), were used for chronic recordings. The animals and data were used for the work by Stark et al. (2013). All animal handling procedures were approved by the Rutgers University and New York University Animal Care and Facilities Committees.

Probes and surgery

Every animal was implanted with a four-shank diode probe as previously described (Stark et al., 2012). Probes were constructed by coupling 470 nm blue LEDs (diameter, 2 mm; model LB P4SG, Osram) to 50 μm multimode optical fibers and attaching every diode–fiber assembly to a single shank of a 32-site/four-shank silicon probe (Buzsaki32, NeuroNexus). Fiber tips were located ∼50 μm above the top recording site. Probes were implanted in the right hemisphere (anteroposterior, −1.6 mm; mediolateral, 1.1 mm) under isoflurane anesthesia. During surgery, the probes were lowered to a depth of 0.4–0.7 mm.

Recordings and photostimulation

After allowing the animals to recover for at least 48 h, recordings were initiated. Recordings were conducted in the home cage during spontaneous behavior. Mice were tethered by one ultralight cable for multichannel neuronal recordings and a second cable for multichannel optical stimulation. Recordings were conducted as the probe was moved gradually from the neocortex to the CA1 pyramidal cell layer. At each location in the brain, neuronal activity was inspected for spontaneous spiking activity, and, if encountered, a baseline period of at least 15 min was recorded followed by photostimulation (peak driving current, 60 mA; mean ± SD peak light power, 35 ± 7 μW; 50–70 ms pulses). Signals were generated by custom code written in MATLAB (MathWorks), converted by a digital signal processor (model RX5 and/or RX6, TDT) to voltage signals, and fed into a linear 16-channel current source. After each session, the probe was either left in place or moved (35–70 μm steps), and the brain was allowed to settle overnight.

Spike sorting and ground truth labels

During recordings, neural activity was filtered (1–5000 Hz), amplified (20× by Plexon headstages; 50× by an RC Electronics system), and digitized (16 bits, 20 kHz) on a 128-channel DataMax recording system (RC Electronics). Applied currents were recorded by the DataMax system. Offline, spike waveforms (32 samples/channel) were extracted from the wide-band records, detrended, and sorted into single units automatically (Harris et al., 2000), followed by manual adjustment. Only well isolated units [(amplitude, >50 μV; L-ratio, <0.05 (Schmitzer-Torbert et al., 2005); inter-spike interval index, <0.2 (Fee et al., 1996)] were considered. A total of 199 neocortical and 781 CA1 units conformed to these criteria.

For connectivity-based tagging, we tagged units that participated as a reference in a cross-correlation histogram (CCH) that exhibited a significant (p < 0.001, Bonferroni-corrected Poisson test) peak in the monosynaptic time range (0–5 ms) as excitatory cells (424 of 980 units). Units that exhibited a significant trough in the monosynaptic time range were tagged as inhibitory (21 of 980 units). For optogenetic-based tagging, units that exhibited a significant (p < 0.001, Poisson test) increase in spiking rate during 50–70 ms DC pulses given on the recording shank were tagged as optically activated cells; 98 of 980 (10%) units conformed to the criterion. Next, we labeled units based on the three tags. Units tagged exclusively as excitatory were labeled as PYR cells (420 units), whereas units that were optically activated or inhibitory (but not excitatory) were labeled as PV cells (102 units). The remaining units (458 of 980, 47%) were not labeled and were discarded from the dataset. Ten of the labeled units (PYR 9 units; PV, 1 unit) were recorded using seven instead of eight channels and were therefore discarded as well, yielding a final dataset that included 512 tagged units (PYR cells, 411 units; PV cells, 101 units; Fig. 1). Of the 101 units referred to as “PV cells,” 93 were optically activated (92%), 13 were both optically activated and inhibitory, and 8 units were only inhibitory. Spike width, firing rate, and bursting behavior were similar for the inhibitory and the optically tagged PV cells. Thus, the 8 inhibitory-tagged PV-like cells were grouped with the 93 optically activated cells, and the entire group was denoted as PV cells. A majority of the units (449 of 512 units) were recorded from CA1. The median [interquartile range (IQR)] number of spikes per unit was 8368 [4494–16,929] for PYR, and 66,850 [13,031–174,802] for PV cells.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

PYR and PV interneurons are tagged in freely moving mice. A, Optical tagging of PV cells. a, Every PV::ChR2 mouse was chronically implanted with a four-fiber/four-shank/32-channel optoelectronic array in the neocortex (nCX). Optical stimuli were applied, in separate sessions, in the nCX and in hippocampal region CA1. Peristimulus time histogram of the PV cell (bottom), triggered by the onset of 50 ms light pulses applied on the optical fiber attached to the recording shank (n = 20; 33 μW). The unit is tagged as PV because of a robust firing rate increase during light (gray) compared with no-light periods. ***p < 0.001, Poisson test. b, Wide-band (1–5000 Hz) recordings from four same-shank channels in CA1. Bottom, Spike trains of a PYR (purple) and the PV cell (green). B, Connectivity-based tagging. Top, Mean (±SD) spike waveforms. For every unit, the channel that exhibits the highest trough-to-peak magnitude is denoted the main channel (boxed). Middle, Auto-correlation histograms (ACHs), showing burst spiking activity of the PYR (purple). Bottom, Cross-correlation histogram (CCH; black) between the spikes of the PYR and the optically tagged PV cell. Gray, Monosynaptic window. The CCH is consistent with monosynaptic excitation of the PV cell by the reference unit, tagging the reference unit as excitatory (PYR). ***p < 0.001, Poisson test. C, Tagged dataset. Of the 512 units in the dataset, 411 (80.3%) are PYR, and 449 (87.7%) are from CA1.

Feature extraction

The shape and timing of the recorded spikes were used to extract a total of 34 features. We derived features of the following three modalities: waveform-based features, derived from a single channel (n = 8); spike-timing features, ignoring the spike waveform (n = 8); and spatial features, derived from the multichannel waveforms (n = 18). All features were based exclusively on spontaneous events that occurred in the lack of any light stimuli. For every spike, the waveform was extracted for 32 samples (1.6 ms) on every channel of the recording shank. The limited duration of the spikes places an upper bound on the classification performance of waveform-based and spatial models. For every channel separately, the waveform was averaged over spikes, and the mean waveform was calculated and upsampled by eightfold using Fourier interpolation to increase the temporal resolution. Since the average and the Fourier transform are linear operators, the order of the two steps does not affect the outcome.

Single-channel feature extraction

For the waveform-based features, the channel with the largest trough-to-peak (TTP) magnitude was denoted as the “main” channel. The waveform in the main channel was scaled by dividing all values by the minimal value (i.e., at the trough). Scaling was performed to remove information about the sampling process (e.g., electrode impedance and neuron–electrode distance). When the absolute value at the trough was smaller than at the peak, waveforms were inverted (multiplied by −1). The outcome is a 256-element vector limited to the −1 to 1 range, with at least one value at −1. To provide a rich description of the waveform, a total of n = 8 waveform-based features were extracted from the main channel (Fig. 2, Table 1): four from the waveform itself, one from the first temporal derivative, and three from the second temporal derivative. For every feature, we compared the distribution of values between all available PYR cells (n = 411) and PV cells (n = 101) and calculated effect sizes. The specific measure of effect size used was the nonparametric estimator for common-language effect size (A_w; Ruscio, 2008) which exhibits a smaller bias compared with alternatives (Li, 2016). A_w estimates the probability that a random sample from one distribution is larger than a random sample from a second distribution. Disregarding the direction of the effect, A_w is thus limited to the 0.5–1 range, taking a value of 0.5 when the two distributions are fully intermixed, and 1 when the two distributions do not overlap at all. All eight features (100%) exhibited a consistent difference (0.64 ≤ A_w ≤ 0.98; p < 0.05, U test; Table 1). Thus, all waveform-based features are potentially useful for classification.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Waveform-based and spike-timing features allow near-perfect classification of PYR and PV cells. A, Features extracted from the mean waveform of the main channel. Voltage values were scaled by the absolute value at the maximal negativity, yielding arbitrary units (AU; the example units are the same as in Fig. 1B). a, Trough-to-peak (TTP) duration feature, defined as the time between the maximal negativity and the ensuing maximal positivity. b, Cumulative distribution function (CDF) of the TTP duration feature for the entire population (PYR, n = 411; PV cells, n = 101; no chunking). Here and in all subsequent CDFs, horizontal lines represent 50%, vertical dashed lines indicate medians, and the *** symbol corresponds to p < 0.001 (U test). The filled circles represent values corresponding to the examples given in a. The difference between PYR and PV cells indicates longer TTP durations for PYR compared with PV cells. B, Waveform-based features allow near-perfect classification. Cross-validated random forest models were trained using the waveform-based features (n = 50 partitions). The chunking method yields improved classification compared with no chunking. The ROC AUCs without chunking (blue) and with 50 spike chunks (orange) are higher than chance level. ***p < 0.001, Wilcoxon test compared with chance level, 0.5. Inset, Confusion matrices (no chunking) based on different decision thresholds (top, 0.1; bottom, 0.9) show variability in prediction, exemplifying the shortcomings of threshold-dependent metrics. C, Features extracted from spike timing. High-frequency features derived from single-sided short-term (0–50 ms) ACHs. a, The Uniform-distance feature is defined as the average absolute difference between the single-sided ACH and a straight line (the example units are the same as in Fig. 1B). b, Cumulative distribution of the Uniform-distance feature for the entire population (no chunking). The larger Uniform-distance values for PYR indicate larger deviations from linear recovery for PYR compared with PV cells. D, Classification based on spike-timing features is not consistently improved by chunking. AUCs were derived from ROC curves based on n = 50 cross-validated random forest models. ROC curves for the test data without chunking (blue) and with 1600 spike chunks (orange) with performance above chance level. All conventions are the same as in B. See also Extended Data Figures 2-1, 2-2.

The addition of a redundant feature would contribute no additional information to the classification process. To estimate relations between features, we first computed rank (Spearman’s) correlation coefficients (CCs) between every possible pair of waveform-based features (Extended Data Fig. 2-1A). To go beyond monotonic relations, we estimated mutual information (MI) between distributions of pairs of features (Timme and Lapish, 2018). If a feature had <10 unique values, the feature was considered naturally discrete; only the spatial dispersion (SPD)-Count feature (see Spatial feature extraction subsection below) was naturally discrete. The use of 10 bins limits the maximal information to log₂10 = 3.3 bits. Deviation from chance level was determined using a permutation test, creating a null distribution based on 5000 iterations of shuffled pairs and evaluating the tail of the null distribution above the observed MI. We found that some features were only weakly correlated with others (e.g., the Break-measure, median [IQR] absolute CC was 0.19 [0.076–0.24]), suggesting that independent information could be gleaned by using the feature. Alternatively, a weakly correlated feature may be dominated by noise. However, the possibility is unlikely since the Break-measure differed for the PYR and PV groups (A_w = 0.64; p = 7.85 × 10⁻⁶, U test; Table 1). Other features were more strongly correlated with the host of other features (e.g., full-width at half-maximum (FWHM), 0.81 [0.45–0.87]; Smile-cry, 0.76 [0.33–0.91]). Since a small number of samples limits the power of standard statistical tests (e.g., the Mann–Whitney U test), when comparing CCs between two groups within a modality we applied a permutation test. We compared a statistic, defined as the difference between the medians of the two groups of CCs, to the 95th percentile of a null distribution created by shuffling the CCs between the groups, and calculating the statistic 1000 times. We did not observe stronger absolute CCs within families: the absolute intrafamily correlation was 0.37 [0.14–0.62], whereas the absolute interfamily CC was 0.47 [0.24–0.79] (p > 0.05, permutation test). Quantifying the interrelations between waveform-based features using MI yielded similar results. The median [IQR] MI between waveform-based feature distributions was 0.658 bits [0.292–0.955] bits, and the rank correlation coefficient between the MI and CCs was 0.816 (p = 0.001, permutation test; Extended Data Fig. 2-1B, inset). The bulk of the variance in the MI (R² = 0.67) could be explained by pairwise rank correlations, suggesting that interrelations between feature pairs are largely monotonic. Thus, based on the feature redundancy analysis, partitioning into families may have only semantic value, and all derived features may contribute to classification.

Spike-timing features

For the features based on spike timing, the autocorrelation histogram (ACH) was calculated for every unit over a range of ±1000 ms using a bin size of 0.5 ms. The ACH depends only on the timing of the spikes and is agnostic to the spike waveforms. The ACH was upsampled eightfold using polyphase filtering to increase the temporal resolution. The ACH is an even function, but edge effects necessarily cause asymmetry in any practical implementation. To obtain an ACH (0–1000 ms) free from edge effects, the value in every positive time bin was averaged with its negative homolog. A total of n = 8 spike-timing features was derived (Fig. 2, Table 2). Three of the features were high-frequency features, derived from the short-term ACH (up to 50 ms). Two were low-frequency features, derived from the long-term ACH (50–1000 ms). The last three were wide-band features: two were derived from the complete ACH (0-1000 ms), and one was derived from the entire recording. For every feature, we compared the values derived for the PYR and PV cells. All eight features (100%) exhibited a consistent difference between PYR and PV cells (p < 0.05, U test; Table 2). Thus, spike-timing features may be useful for classification.

To quantify interrelations, we computed CCs and MI between spike-timing features extracted from the complete spike trains (Extended Data Fig. 2-2). The high-frequency features exhibited high within-family absolute correlations (median [IQR], 0.85 [0.83–0.88]), whereas lower absolute correlation values were observed between the other features (low-frequency and wide-band families together: 0.26 [0.24–0.70]; p = 0.098, permutation test; Extended Data Fig. 2-2A). In contrast to the waveform-based feature families, intrafamily absolute correlations (0.81 [0.49–0.83]) were higher than interfamily absolute correlations (0.50 [0.40–0.61]; p = 0.012, permutation test). MI between pairs of spike-timing features yielded results similar to those for pairwise correlations. The median [IQR] MI between spike-timing features was 0.469 [0.283–0.701] bits, and the rank correlation coefficient between MI and CCs was 0.967 (p = 0.001, permutation test; R² = 0.93; Extended Data Fig. 2-2B, inset). The fact that almost all variance of the MI values is explained by pairwise correlations suggests that the interrelations between the pairs of features are largely monotonic. The correlations between high-frequency features, namely Uniform-distance, Kullback–Leibler distance (D_KL)-Short, and Rise-time suggest that the features provide redundant information. Thus, a small number of spike-timing features may suffice for classification.

Spatial feature extraction

For extracting purely spatial features, an event-based δ-transformation was first applied to the mean upsampled waveform of every channel to remove all waveform-based information (Fig. 3A), as follows. (1) First, positive spikes were inverted as done for the waveform-based process. (2) Next, three events were defined. One event was the time of maximal negativity (NEG). For the additional two events, the median over all the channels was calculated. The second event was the first median crossing (FMC), the first time point before the maximal negativity of the channel in which the global median was crossed. The third event was the second median crossing (SMC), the first time point after the maximal negativity of the channel in which the global median was crossed. Every event was detected on every channel, yielding a total of 24 points. (3) Third, the waveform was replaced by a δ-like function that took the value of the maximal negativity of the channel at the singular event time point and zeros everywhere else. The δ-like functions were scaled by the absolute value of the global minimum over all channels, effectively removing all waveform-based information from every single channel. However, waveform-based information may still be available when using multiple events together (e.g., FMC, NEG, and SMC). (4) To remove all waveform-based information, the δ functions were shifted together to centralize (shift to the 129th sample) the time of the event on the main channel. The process transforms the waveform in the main channel to be nearly identical for all units (Fig. 3B). Residual variability in the time of the trough may remain if the channel with maximal magnitude of the TTP and the channel with the maximal negativity are not the same.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Transforming multichannel spike waveforms to event-based δ-like functions removes all waveform-based information and allows extracting purely spatial features. A, The event-based δ-transformation procedure, illustrated for the FMC event. a, The mean waveforms, with δ-like functions marking the FMCs. The transformation replaces all voltage values with zeros, except for the event points, which are assigned the same value as the trough. In gray are channels for which the magnitude of the TTP is below a predetermined threshold (Materials and Methods). The main channels are boxed. b, Next, waveform-related information that might be recovered by combining multiple δ-transformed events is removed. The δ-like functions are scaled and centralized (arrowheads), placing the event of the main channel at the midpoint (129th sample). B, Left, The scaled waveform of the main channel of all units in the dataset before the transformation, sorted for PYR and PV cells separately by the time of the trough. Right, The same waveforms after event-based δ-transformation. The transformation removes nearly all of the variability between units. C, Cross-validated random forest models (n = 50; no chunking) were trained using waveform-based features extracted from the transformed spikes. The confusion matrix, based on a naive decision threshold of 0.5, yields a constant prediction of one class. n.s.p > 0.05, Wilcoxon test. Numbers in every cell denote the median [IQR]. Performance was quantified by the threshold-independent AUC. The classification yields an AUC of exactly 0.5, corresponding to purely random prediction. D, A time-based feature, FMC-Time-lag-SD, derived from the differences between the times of the FMC event in different channels. The feature quantifies the temporal dispersion of the event, without considering the actual positions of the recording electrodes. a, FMC-Time-lag-SD is defined as the SD of the time differences between the FMC event of the main channel (vertical dotted lines) and the other channels. In gray are ignored channels, for which the magnitude of the TTP was below a predetermined threshold. b, Cumulative distribution of the FMC-Time-lag-SD feature for the entire population (411 PYR, 98 PV cells, no chunking). The smaller FMC-Time-lag-SD values of the PYR indicate higher spatiotemporal synchrony for PYR compared with PV cells. All conventions for the CDFs here and in subsequent panels are the same as in Figure 2A. E, A graph-based feature, FMC-Average-weight, derived from the differences between the FMC event time in different channels and the electrode locations. a, FMC-Average-weight is defined as the average edge weight in the event graph. The event graph is a directed graph with vertices representing the electrodes, and edges representing the transmission speed based on the timing of the events and the location of the electrodes. Only channels that passed the threshold for the magnitude of the TTP were considered. b, Cumulative distribution of the FMC-Average-weight feature (no chunking). The larger values for PYR indicate higher transmission rates for PYR compared with PV cells. F, A value-based feature, SPD-Count, derived from SPD of the maximal negativity on every channel. a, SPD-Count is defined as the number of channels that reached at least 50% of the maximal negativity of the main channel. b, Cumulative distribution of the SPD-Count feature (no chunking). No consistent difference between the PYR and PV cells is observed, suggesting similar spatial distributions of the scaled maximal negativity (p > 0.05, U test). See also Extended Data Figure 3-1.

Overall, n = 18 features were extracted from the transformed waveforms (Table 3) associated with the three events, quantifying the following three dimensions: time based, graph based, and value based. (1) In the time-based dimension, only the timing of the events was considered (e.g., the SD of the FMC; Fig. 3D, FMC-Time-lag-SD). Channels for which the magnitude of the TTP (before δ-transformation) did not pass an arbitrary threshold of 25% of the maximal magnitude of the TTP over all channels were ignored (Fig. 3D, gray). A median [IQR] of 4 [3–5] PYR channels and 3 [2–5] PV cell channels were removed. (2) The graph-based dimension included both event timing and the physical locations of the recording electrodes on the probe. An event graph was generated based on a specific event, with a node for each channel (e.g., Fig. 3E, FMC-Average-weight). Only channels for which the magnitude of the TTP passed the 25% threshold were considered. Directed edges connected every two nonoverlapping events, with a weight representing “transmission speed”: the Euclidean distance between the electrodes, divided by the time difference between the events. (3) The value-based dimension (SPD) ignored timing information and considered the scaled maximal negativity values on every channel, based on the global maximal negativity (e.g., Fig. 3F, SPD-Count). We found that 10 of 18 (56%) of the spatial features exhibited differences between PYR and PV cells (p < 0.05, U test). Although some features do not show consistent differences between the two cell types, classification may benefit from the features because of, for instance, distinct second-order statistics.

To estimate feature redundancy because of high correlations, we computed the CCs between the spatial features extracted from the transformed mean waveforms (Extended Data Fig. 3-1A). The rank correlation matrix of the spatial features showed absolute correlations (median [IQR]: 0.2 [0.1–0.33]), that were weaker than for the waveform-based features (p = 3.5 × 10⁻⁴, U test) and for the spike-timing features (p = 1.1 × 10⁻⁷). Eighty percent of the spatial feature pairs exhibited absolute correlations higher than zero (122 of 153; p < 0.05 permutation test). Intrafamily absolute correlations (0.27 [0.17–0.42]) were higher than interfamily correlations (0.17 [0.08–0.29]; p = 0.006, permutation test). The median [IQR] MI between spatially based feature distributions was 0.19 [0.155–0.266] bits, and the rank correlation coefficient between MI and CCs was 0.84 (p = 0.001, permutation test; R² = 0.706; Extended Data Fig. 3-1B). Since the correlation and the MI analysis agreed, the relatively weak correlations between spatial features may result from large amounts of noise in every feature. Alternatively, the features may provide independent information, useful for classification.

Classification procedure

The classification model was chosen to be random forests (Breiman, 1996, 2001) because of the relative simplicity. Furthermore, several methods are available for understanding the determinants of a specific random forest model prediction (Archer and Kimes, 2008). To achieve good estimation of model performance, a nested cross-validation procedure was applied (Varma and Simon, 2006; Krstajic et al., 2014). For every modality (waveform, spike timing, and spatial), the training procedure was repeated n = 50 times. In every iteration, data were first partitioned in an approximate 80:20 ratio into training and test sets in a stratified fashion. Thus, the training set always included 328 PYR and 80 PV cells, and the test set included 83 PYR and 21 PV cells; only the identity of the units changed between iterations. To handle the imbalance between the number of PYR and PV cells in the dataset, the model weights that control the effect of every class on the impurity score used for training the random forest model were adjusted. Specifically, instead of assigning equal weights, class weights were set to be inversely proportional to the number of samples in every class using the following: total number of training set samples/(number of classes × number of class samples in the training set). Second, a fivefold grid search was conducted on the training set to find the best hyperparameters for the model, optimizing the receiver operating characteristic (ROC) of the area under the curve (AUC). The tested hyperparameters were the number of estimators, the depth of each estimator, the minimal number of samples required to split a node, and the minimal number of samples required to be at a leaf node. Other hyperparameters received default values based on the implementation of the scikit-learn library in Python (Pedregosa et al., 2011). Third, using the optimized hyperparameters, the model was trained on the entire training set. Finally, model performance was evaluated using the test set.

Performance and explainability

Model performance was assessed using a metric that is robust to unbalanced data. Two types of metrics exist: threshold dependent and threshold independent. Threshold-dependent metrics consider only the binary decision: in our case, PYR or PV. Threshold-independent metrics consider the raw prediction, a value between 0 and 1, and not the decision itself. Threshold-dependent metrics require choosing a threshold, and thus the outcome may vary according to the chosen decision threshold. An arbitrarily chosen threshold does not necessarily reflect performance, especially when considering unbalanced datasets (Sheng and Ling, 2006). The choice of a decision threshold is not trivial, and is the subject of active research (Esposito et al., 2021). For these reasons, we used an established threshold-independent metric, the AUC (Fawcett, 2006). The AUC reflects the relation between true positives and false positives for all possible thresholds and is hence threshold independent as well as suitable for unbalanced datasets.

The theoretical chance level for the AUC metric is 0.5. To determine the empirical chance level, the performance of models trained on data with shuffled training-set labels was assessed. For every modality (waveform, spike timing, and spatial), the training procedure was conducted with the labels shuffled only for the training set. The AUC exhibited the expected chance-level behavior, yielding the following median [IQR] values: waveform, 0.46 [0.34–0.57]; spike timing, 0.46 [0.39–0.56]; and spatial, 0.48 [0.44–0.55]. Since chance-level results were obtained on the test sets, further assessment of model performance was conducted with respect to 0.5, the theoretical chance level.

To assess the contribution of every feature to the final classification prediction, we used Shapley additive explanations (SHAPs) values (Lundberg et al., 2020). The SHAP method uses a game-theoretic approach to calculate an importance score for every feature in each sample, considering all interdependencies with other features. Since the predictions of random forest models range from 0 to 1, the contribution of every feature can take values from −1- to 1. For each model, we calculated and averaged the absolute contribution of every feature based on the test set, taking the median over the 50 partitions to avoid idiosyncrasies of a single arbitrary partition. Since the absolute value is taken, SHAP values are necessarily non-negative, creating a skewed distribution that does not have an expected value of zero. To calculate significance, we used a permutation test. The SHAP values for each feature were compared with the SHAP values obtained for the models trained with shuffled training-set labels. Models trained with shuffled labels represent chance-level classifiers, for which the importance of every feature can be considered as the chance-level baseline. To create a null distribution, we partitioned the dataset into a train and test sets 1000 times. For every train-test partition, we shuffled the labels, trained the models as described previously, and calculated the SHAP values. Then, to calculate a p-value for every feature, we compared the original SHAP value to the null distribution.

Chunking method

Models trained on larger nonredundant datasets typically exhibit improved performance. Therefore, data augmentation approaches to synthetically increase the size of the dataset are often applied (Moreno-Barea et al., 2018). Augmentation may be implemented by adding noise or transforming the data (rotation and reflection in image classification; Mikołajczyk and Grochowski, 2018). Here, to augment the size of a given dataset, features were extracted from “chunks” that included subsets of spikes, instead of using all available spikes together. Thus, we increased variability using the natural richness of the data that is otherwise flattened to a single mean waveform. The chunking process increases the number of samples in the dataset, with the possible cost of increased noise.

A total of n = 7 different chunk sizes was used, with 25, 50, 100, 200, 400, 800, or 1600 spikes per chunk. For a given unit with N spikes and a chunk size C, the spikes were randomly split into ⌊N/C⌋ chunks so that every chunk consisted of between C and 2C – 1 distinct spikes. Spikes were randomly assigned to chunks. In a given split, all chunks consisted of the same number of spikes up to a difference of a single spike. Every chunk received the label of the source unit. Because of the different numbers of spikes recorded for PYR and PV cells (5,651,196 PYR spikes and 11,612,978 PV cell spikes), the balance between PYR and PV samples changed compared with the original dataset (411 PYR and 101 PV cells). For instance, using a chunk size of 25 spikes, the total number of PYR samples was 225,850 while the number of PV samples was 464,473.

To determine features for every chunk separately, the waveforms were averaged over all available same-channel spikes within a chunk, from which waveform and spatial features were derived. To derive chunk-specific spike-timing features, the ACH was accumulated by summing over all single-spike ACHs. For every spike, the single-spike ACH was based on all spikes that occurred in the −1000 to 1000 ms time window around the reference spike. To be applicable to the chunking case, the single-spike firing rate was defined as the mean of the inverse interspike intervals before and after the spike.

To provide information about the distribution of the values over chunks, several statistics were extracted from the individual values of every chunk. Extracting statistics based on all the chunks of a unit allows considering intraunit variability as a feature. For every feature, the mean, SD, and the 25%, 50%, and 75% quantiles were extracted (referred to as “chunk statistics”). Note that extracting the mean out of all the individual feature values of the chunks is not the same as not using chunking. Without chunking, first all the waveforms are averaged, transformed, and then features are extracted. In contrast, when using the mean over chunks, averaging happens at the end of the process (averaging the waveforms within a chunk still happens at the beginning). If all steps are linear, the two processes yield the same results. However, since feature extraction is not a linear operator, the mean statistic may contain unique information. The specific statistics were chosen to capture the first and second moments of the across-chunk distribution and for their simplicity. The chunk statistics can be expanded further to capitalize on different properties of the distribution over different chunks. The five new chunk statistics increased the total number of features used by the model sixfold. Notably, the chunk statistic features are the same for all the chunks of a given unit.

When training models using chunks, the data were partitioned based on units. In the training set, every chunk was considered independently (not as part of the unit). In addition, instead of performing a grid search for every chunk size, the hyperparameters for all chunk sizes were chosen by a no-chunking equivalent. Chunk statistics were extracted for the no-chunking dataset as well, so the number of features was equal to the number of features in the chunking method. The models applying chunking used the hyperparameters found using the no-chunking equivalent grid search, based on the same partition of the data to training and test sets. For testing and evaluating, every chunk received an independent prediction. Then, predictions were pooled over all same-unit chunks by casting a majority vote, yielding a final chunk-based prediction for the unit.

When calculating the SHAP values for the chunking-based models, we randomly chose 1000 samples out of the test set. The procedure for computing SHAP values for chunking-based data were otherwise the same as for the no-chunking data. The absolute SHAP value of every feature was summed together with the values of the chunk statistics extracted from the same feature, yielding a single importance value for every original feature.

Generalization analysis

Units in the tagged dataset were recorded from CA1 (449 of 512 units, 88%) and from neocortex. The availability of tagged data from two brain regions allows testing and quantifying interregion generalization. Generalization was determined by the performance of models trained using recordings from one region on test data from the same region (“training region”), and from the other region (“non-trained-on region”). To directly quantify generalization, we partitioned the full dataset into the following three sets for each training region: (1) a training set, containing ∼80% of the training region units (CA1: 301 PYR, 58 PV; neocortex: 27 PYR, 23 PV); (2) a test set, containing the remaining 20% of the training region units (CA1: 76 PYR, 14 PV; neocortex: 7 PYR, 6 PV); and (3) a second test set, containing all units of the non-trained-on region (neocortex: 34 PYR, 29 PV; CA1: 377 PYR, 72 PV). Waveform, spike timing, and spatial models were trained on the reduced CA1 dataset with 50, 1600, and 25 spike chunks, respectively (found to yield the best performance for each modality on the combined dataset). Chunked data were used for the grid search: the training set was further partitioned into a “development set” containing 80% of the units and an “evaluation set” containing 20% of the units. If the development set contained >5000 chunks, the grid search was conducted on a random subset of 5000 chunks, minimizing run time while allowing an efficient search. The difference in the number of units between the CA1 and neocortical training sets leads to an inherent difference in absolute performance between the two training region conditions. However, generalization can be readily compared between the two conditions based on the performance of the test set of the non-trained-on region relative to the performance of the test set from the training region.