Clustering and control for adaptation uncovers time-warped spike time patterns in cortical networks in vivo

How information in the nervous system is encoded by patterns of action potentials (i.e. spikes) remains an open question. Multi-neuron patterns of single spikes are a prime candidate for spike time encoding but their temporal variability requires further characterisation. Here we show how known sources of spike count variability affect stimulus-evoked spike time patterns between neurons separated over multiple layers and columns of adult rat somatosensory cortex in vivo. On subsets of trials (clusters) and after controlling for stimulus-response adaptation, spike time differences between pairs of neurons are “time-warped” (compressed/stretched) by trial-to-trial changes in shared excitability, explaining why fixed spike time patterns and noise correlations are seldom reported. We show that predicted cortical state is correlated between groups of 4 neurons, introducing the possibility of spike time pattern modulation by population-wide trial-to-trial changes in excitability (i.e. cortical state). Under the assumption of state-dependent coding, we propose an improved potential encoding capacity.

(a) For each single neuron and stimulus condition combination, first spike times are normalised (combination mean subtracted followed by division by combination standard deviation). Means of these normalised first spike times are calculated for neuron groups and 10 trial bins (i.e. trials 0-9, ..., 90-99). These mean values are plotted with a 95\% confidence interval. Values for excitatory and inhibitory neuron groups shown in red and blue respectively. Note that L2/3 is undersampled due to sparse activity and responsiveness of L2/3 neurons (e.g. de Kock and Sakmann, PNAS, 2009). (b) Same as (a) but for normalised spike counts. (c) Table displaying the number of spike sorted neurons recorded over all experimental sessions by layer and E/I type.
de Kock, C. P., & Sakmann, B. (2009). Spiking in primary somatosensory cortex during natural whisking in awake head-restrained rats is cell-type specific. Proceedings of the National Academy of Sciences,106(38), 16446-16450.  Fig. 3 Heterogeneous first spike rasters following cortical onset for five single-neuron single-stimulus combinations. L4 INH (i) and L5A EXC neuron (ii) demonstrate reliable short-latency response between 5 and 10 ms poststimulus with prominent adaptation. L2/3 EXC neurons (iii and iv) also demonstrate adaptation but responses are more variable. In contrast, L5A EXC neuron (v) shows variable responses but not spike-time adaptation.  Fig. 4. Single-neuron single-stimulus first spike trial-lagged autocorrelations and partial autocorrelations. Trial-lagged autocorrelation and partial autocorrelation analysis allows quantification of further autocorrelative structure in the trial-to-trial sequence of first spike times. Here we consider the first spike times s t of single neurons on spiking trials t=1,…,T of stimulus condition where T is the number of spiking trials. Trial-lag-n autocorrelations and partial autocorrelations test the correlations and partial autocorrelations between spike times on trials t with spike times on trial t+n.
Plots show the normalized cumulative sums of trial-lag 1-5 autocorrelation (left) and trial-lag (2-5) partial autocorrelation (center) p-values for 2416 single neuron single stimulus combinations. Trial-lagged autocorrelations (left) were above chance for trial-lags 1-3 supportive of general adaptative trends. This is demonstrated by the fact that the cumulative sums are slightly above y=x in the plot, demonstrating non-uniform distributions of p-values. As this is hard to see, accompanying tables show results (χ 2 statistics and p-values ) of Fisher-method application to the p-values produced for each trial-lag. This confirms that autocorrelations and are above chance for several values. This is not the case for trial-lagged partial autocorrelations, supporting the presence of general adaptative trends without additional structure in the trial-to-trial sequence of single neuron first spike times.  Mean θ 45 angles with 95% empirical confidence intervals for positively correlated Stage 2 clusters detected in angled surrogate control dataset. Blue if significantly greater than 0°(p-value < 0.025) and less than 45°(p-value < 0.025). Red otherwise. Detection of clusters different from 45°appear to be at chance-levels. (c) Histogram of p-values for correlation angles being different from 45°for Stage 2 clusters detected by the main clustering algorithm for the generated angle surrogate control dataset. Correlation angles different from 45°are not found above chance-levels. The KPSS test was used to test whether the cluster first spike times of single neurons were stationary around a deterministic trend. Single neuron cluster first spike time sequences were determined as not stationary around a deterministic trend if the null hypothesis was rejected (p < 0.05). Python`statsmodels' implementation of p-value estimation rounds all p-values > 0.1 to 0.1. Supplementary Fig. 8 (1 of 3). Stationary Stage 2 cluster single-neuron single-stimulus first spike (a) autocorrelations, (b) partial autocorrelations, (c) cross-correlation and (d) partial cross-correlation. Autocorrelation, partial autocorrelation, cross-correlation and partial cross-correlation analysis allows quantification of autocorrelative and cross-correlative structure in the trial-to-trial sequence of first spike times and first spike pairs for the 100stationary Stage 2 clusters. Here we consider the first spike times s j,t of neurons j=0,1 on cluster trials t=1,…,T where T is the number of cluster trials for a single cluster. Trial-lag-n autocorrelations and partial autocorrelations test the correlations and partial autocorrelations between single neuron cluster spike times on trials t with spike times on trial t+n. Plots (a) and (b) show the normalized cumulative sums of trial-lag 1-5 autocorrelation and trial-lag (2-5) partial autocorrelation p-values for 200 (=2 * 100 stationary Stage 2 clusters) single neuron single cluster spike time sequences. Trial-lag-n cross-correlations and partial cross-correlations were similarly tested between neuron pair cluster spike times on trials t with spike times on trial t+n.
Autocorrelations, partial autocorrelations, cross-correlations and partial cross-correlations are at chance-levels for trial-lags 1-5. This is demonstrated by the fact that the cumulative sums lie on y=x, demonstrating uniform distributions of p-values. Accompanying tables show results (χ 2 statistics and p-values ) of Fisher-method application to the p-values produced for each trial-lag, which confirms that autocorrelations, partial autocorrelations, cross-correlations and partial cross-correlations are not above chance.
Subfigures a-d are equivalent to 6b-e but for time series modelled correlated unclustered response distributions.
Supplementary Fig. 10 (1 of 2). Autocorrelations (left) and cross-correlations (right) at chance-levels for criteria fulfilling ARIMA modelled Stage 2 clusters. As with Supplementary Fig. 6 (a) and (c) but for the 41 criteria fulfilling ARIMA modelled clusters. Autocorrelations and cross-correlations are also at chance levels confirming that ARIMA modelling successfully removed autocorrelations introduced by differencing and that there is no additional cross-correlative structure between the neurons in ARIMA modelled clusters.   Fig. 11 (2 of 2). Same as Figure 7 but for stationary and made-stationary correlated unclsutered response distributions (rather than Stage 2 clusters). Spatial separation and E/I-type of all correlated unclustered response distribution neuron pairs Supplementary Fig. 12. Onset of cortical activation following whisker stimulation. Histogram of single neuron spike times to all stimuli with bins of 0.2 ms and between -50 ms before and +100 ms after stimulus. The histogram was smoothed with a gaussian filter (width s.d. 1 ms). The mean and standard deviation of the smoothed bin counts was calculated for the period before cortical activation onset at 0 ms. The cortical activation onset spiking threshold was determined as this mean value + 3 standard deviations. Cortical onset was determined as the time point when the smoothed histogram went above this threshold (red dot). The left and right figures show the same smoothed histogram between -50 and +100 ms (left) and 0 and +25 ms (right).  Fig. 13. Example clustering of sample first spike response distribution for two neurons (L5AI, L4I) for single epsilon value over a 50ms window following cortical onset (1 outlier out of view). (a) Point colours illustrate DBSCAN clustering (epsilon = 1.35). Black points are outliers. Dashed ellipses represent 3σ intersection ellipse for each DBSCAN cluster. Setting the DBSCAN `minimum number of elements in a cluster' parameter to 2 increased the cluster quality control, as intersection ellipses were created for 2 element clusters rather than assigning them as outliers. (b) Dark green flat ellipse illustrates 4σ confidence interval of Stage 1 cluster. Light green angled ellipse illustrates 4σ confidence interval of Stage 2 cluster. Orange ellipse illustrates 4σ ellipse of all spike pairs under the assumption of independence and normality.