Refractoriness Accounts for Variable Spike Burst Responses in Somatosensory Cortex

Models with refractoriness can reproduce the experimental spike trains. A, Intensity function (left) and recovery (right) functions of the STPM fitted to experimental data (an example for a single neuron). B, Comparison of PSTHs of the training data (top, dark blue line), validation data (light blue line), and model data (red line). Note the overlap between the lines, which is a sign of the match between the model and both the training and validation sets. The difference between the model PSTH from the validation PSTH (model residuals, bottom) is equivalent to the intrinsic variation between the training and validation set (F = 1.02, p > 0.01, see Methods for definition). C, Fitted intensity (left) and recovery functions (right) of the GLM (bin size, 0.25 ms). D, Correlation coefficients between the residuals (for the STPM shown in B, bottom panel) and the validation PSTH (for the STPM shown light blue in B) calculated for three different models: the STPM, the STPM without refractoriness (non-refr. STPM), and the GLM. Box plots represent the distribution of bootstrapped correlation coefficients: boxes - quartiles, notch - median with confidence intervals, whiskers - range. The horizontal dashed line denotes the correlation coefficient between the difference of PSTH of validation and training dataset with the training dataset PSTH. E, The empirical cumulative distribution of the ISIs of the experimental spike trains rescaled according to the conditional intensity function of all three fitted models (time-wrapping test). If the model perfectly reproduced the experimental ISIs the cumulative distribution should line up with the diagonal. This procedure was repeated for two different bin sizes (0.05 ms, left; and 0.25 ms, right). F, The Kolmogorov-Smirnov (K-S) distance of the model (maximum divergence of the model’s cumulative distribution from diagonal in E) decreased with an increasing bin size in both STPM and GLM. This dependence on bin size was less pronounced for the GLM.

The STPM explains trial-to-trial variability of the data. A, Single-neuron responses averaged over all trials (PSTH, same as in Fig. 1B) reveal that spikes occur preferentially at discrete latencies (delimited by vertical lines and indexed by x for the first peak, y for the second peak, and z for the third peak). B, In single trials, multiple spikes are elicited in diverse combinations of preferred latencies resulting in significant trial-to-trial response variability. Spike combinations are classified into spike patterns. The time axis was first divided into three windows aligned to the peaks of the PSTH. Each trial was then assigned a binary string (spike pattern xyz, from 000 to 111), where 1 represents the occurrence and 0 the absence of a spike in a window. Spike timings of eight representative sample responses assigned to each pattern are shown as raster plots. C, Frequency at which the spike patterns occurred over repeated trials for the neuron in A. D, Firing pattern distribution obtained from the data (white bars, same as C), the STPM (red bars) and the nonrefractory STPM (blue bars). The firing rate of the Poisson model was estimated by a PSTH with bin size 0.05 ms. Inset compares the PSTHs obtained from each model (color-coded like the bars in main panel). E, Scatter plot of two consecutive ISIs within spike triplets calculated from the experimental data (filled circles) and responses simulated with STPM (empty circles). Serial correlations (Pearson’s correlation coefficient) found in the experimental intervals (r_data) differ only slightly from the respective correlations predicted by the STPM (r_model, see values in the legend, solid and dashed lines represent the best linear fit to the experimental and model data, respectively). F, Repeated Monte-Carlo simulations (n = 1000) of the STPM fitted to experimental data provide the distribution of serial correlations consistent with the STPM (empty bars); the serial correlations estimated directly from experimental data (vertical arrow, r_data) are likely to be drawn from the same distribution (two-sided bootstrap test, p = 0.81).

Input modulation may explain deviations of spike train statistics from the STPM and GLM. A, The STPM was extended by including a multiplicative gain factor, which acts on the input function. The gain factor was randomly selected from a uniform distribution [0.2, 1.8] in each trial. The model was simulated with an exponentially decaying intensity function (B, dashed line); maximum amplitude, 4000 spikes/s, time constant, 3 ms; and a step recovery function (C, dashed line); refractory period, 1.4 ms. B, Intensity function of the STPM (red) and GLM (blue) fitted to the simulated spike trains. The intensity function manifests deviations from the real intensity function used in the simulation (dashed line; gray-shaded area corresponds to the amplitude range of intensity function taking into account the gain factor). C, Recovery function of the STPM (red) and GLM (blue) fitted to the simulated spike trains. The STPM-estimated recovery function displays a characteristic overshoot soon after end of the absolute refractory period (1.4 ms; dashed line, real recovery function underlying the spike trains). D, RMS error of the intensity and recovery functions estimated with the STPM and GLM. E, The serial correlation of the model with gain modulation (arrow) is significantly larger (p < 0.01) than predicted in absence of modulation (bar plot, histogram of 1000 serial correlation coefficient obtained from Monte Carlo simulations of the STPM with the intensity and recovery functions shown in B, C, red line).

Coordination of spike patterns in the population. A, Simulation of 5000 identical units described by the STPM (Fig. 2C) with gain modulation of the strength γ = 0.2. From the simulated spike trains of all neurons (short ticks represent sample spike times) the population PSTH was calculated (bottom, black line) and then high-pass filtered to obtain an estimate of the hf-EEG population response (bottom, red line). B, Distributions of spike patterns of a single neuron in 1000 repetitions of the simulation with a low (0.8, blue) and a high (1.2, red) gain. C, The population PSTH before (top panel) and after high-pass filtering (bottom panel) varies with the gain (blue: 0.8, red: 1.2). D, Single-neuron spike pattern and root-mean square (RMS) amplitude of the high-pass filtered population PSTH are correlated because both the spike pattern and the PSTH depend on the gain. Box plots represent distribution of single-trial RMS amplitudes for each spike pattern of a single neuron (boxes - quartiles, horizontal lines - medians, whiskers - range excluding outliers, crosses - outliers). E, The simulated population RMS amplitudes correlate with experimental hf-EEG RMS related to the same pattern (hf-EEG RMS; Telenczuk et al., 2011).