Refractoriness Accounts for Variable Spike Burst Responses in Somatosensory Cortex

Refractoriness explains the intraburst intervals

We demonstrated that the STPM with a decaying intensity function and an absolute refractory period can produce a PSTH that agrees qualitatively with the responses of neurons in S1 of macaques. To test whether the STPM can also quantitatively reproduce the fine details of neuronal responses recorded in vivo, we inferred the intensity and recovery functions directly from the data. The two functions were defined on per-bin basis and were treated as the free parameters of the model. These parameters were then fitted to the experimental spike trains using a convex optimisation technique guaranteeing the identification of the most optimal model (see Methods; Fig. 2A).

The fitted intensity function peaks shortly after the stimulus onset (<10 ms) and decays back to baseline when the burst is terminated (Fig. 2A, left). The intensity function still contains three distinct peaks, but they are less prominent compared to the peaks in the PSTH (Fig. 1B). This smoothing can be attributed to the decoupling of synaptic inputs, which are represented by the intensity function, from the refractoriness, which is represented by the recovery function (Fig. 2A, right). Although the maximum of the intensity function is much above the rate at which individual neurons can fire spikes, the refractoriness limits the firing rate of the model neuron. In agreement with the properties of biological neurons, the fitted recovery function is equal to 0 for the first 1 ms after emitting a spike (absolute refractory period), but after a few milliseconds fully recovers from the refractoriness returning to the rest state (w(t) ≈ 1). Interestingly, immediately after the absolute refractory period the recovery function over-shoots for ∼1 ms, largely exceeding the rest value. The fast (≪1 ms) fluctuations following this over-shoot represent statistical noise due to the finite size of the data set. Altogether, the parameters of the STPM disentangle the synaptic inputs from the refractory effects.

The STPM provides a parsimonious description of bursting in S1 cortex

The simulated peristimulus time histogram (Fig. 2B, red line) matches closely the one obtained from the experimental data (Fig. 2B, dark blue line). To demonstrate that this good match is not a result of an over-fitting, we performed cross-validation. First, the data set was divided randomly into two subsets: training data and validation data. The model was fitted only to the first subset, and then the results of the simulation were validated on the second (Fig. 2B, light blue line). We found that the difference of the fitted PSTH from the validation set was of the same magnitude as the variation within the dataset (see Methods; F test, p > 0.01). This test indicated that the model optimally captured the features of both training and validation set without considerable over-fitting.

The parameters of the model were fitted to each of the 17 neurons yielding similar results. Importantly, an application of the cross-validation procedure revealed that in 12 out of the 17 neurons the PSTH simulated with the model was not significantly different from the PSTH calculated from the recorded spike trains (F = 0.65–1.59, p > 0.01). In the remaining five neurons the modeled PSTH deviated significantly from the validation PSTH (F = 2.33–4.88, p < 0.01, F test). This sub-population of neurons may have differing firing properties that would need more sophisticated models (implementing, for example, bursting mechanisms; see Discussion).

To further analyze the cases in which the simulated spike trains differed from the data, we subtracted the model PSTH from the validation data PSTH (Fig. 2B, bottom). The resulting residuals still contained fluctuations aligned to the peaks of the PSTH. This indicated that the model does not fully capture the shape of the PSTH. Indeed, the correlation coefficient between the residuals and the PSTH of the validation set was significantly positive (bootstrapped 95% confidence intervals; Fig. 2D) meaning that the residuals contain some remnants of the averaged neuron response. Altogether, these results show that the STPM model is sufficient for describing trial-averaged responses in the majority of recorded neurons.

Poisson-like variability explains the occurrence of temporal spike patterns in repeated trials

Having shown that the interplay between the intensity and recovery functions of the STPM can account for a large part of the trial-averaged response of a single neuron, we tested whether the model can also explain the trial-to-trial variability of the spiking of cortical neurons.

To quantify the trial-to-trial variability of neuronal responses, we sorted single-trial spike trains according to the occurrence of spikes in predefined temporal windows (Fig. 3A–C). Each spike train was assigned a binary word based on occupancy of preferred firing windows the borders of which were aligned to the troughs of the PSTH (Fig. 3A, windows labeled x, y, and z). As explained above, these troughs reflect the periods of quiescence due to the refractoriness of the neuron. When the single-trial spike trains were re-ordered according to the associated binary word, we could distinguish between several patterns of activity (spike patterns). In most trials the neuron fired in all three windows (triplet, 111) or only the first two (doublet, 110), where the input was the strongest (compare with Fig. 2A), but also doublets with other combinations of spikes and silences were common (Fig. 3C, spike pattern frequency distribution). For example, the doublet 101 corresponds to trials in which the neuron fired at the onset of the stimulation (Fig. 3A, window x), then remained silent during the second window (y) and fired again in the third window (z); the omission of the spike in the window y is the consequence firing late in the window x (see the corresponding line in the raster plot; Fig. 3B), so that the neuron is refractory during the window y. In other neurons the number of discrete firing windows (determined by the number of PSTH peaks) ranged from two to four, and a similar distribution of spike patterns was obtained. The appearance of such spike patterns can be attributed to the chance phenomena (Poisson firing) and their interplay with the structured input and refractoriness.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

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).

We found that the distribution of spike patterns in the experimental data were similar to the distribution obtained from the STPM (Fig. 3D). In contrast, when the recovery function was constrained to 1 for all bins and the intensity function estimated from the data (nonrefractory STPM) some spike patterns appeared at frequencies much different from the experimental data (e.g., spike patterns 010, 110, 001, and 111; Fig. 3D), despite the fact that the overall PSTHs were almost identical (Fig. 3D, inset). The differences of spike pattern frequencies can thus be understood as the effect of refractoriness; without it the probabilities of firing in each window are independent of the occurrence of spikes in the previous windows, in which case the frequency of a spike pattern can be directly predicted from the trial-averaged response (PSTH).

To quantify the similarity between the experimental and modeled spike patterns, we used a cross-validated χ² statistics (see Methods, Eq. 7). In 12 of 17 examined neurons the spike-pattern distribution of the STPM was similar to the experimental distribution, and for five cells they were significantly different (F test, p < 0.01); in two of these five neurons the PSTH was not accurately predicted by the STPM precluding the possibility of predicting the trial-to-trial variations. In the remaining three neurons there were substantial differences in the frequency of selected spike patterns, which might reflect the misestimation of the recovery function. Overall, these results show that in most neurons the STPM with time-dependent inputs and refractoriness can account not only for the trial-averaged but also the trial-to-trial variability of responses to somatosensory stimulation.

Within-burst intervals manifest significant correlations

Next, we investigated whether the correlations between consecutive ISIs (serial correlations) may play a role in the generation of spike patterns. The STPM predicts that the response should be fully determined by the current input and the interval since the last spike. However, the calculation of the serial correlations in the experimental data obtained from S1 showed that two consecutive ISIs are not independent (Fig. 3E). Since significant serial correlations might be induced by the firing-rate variations alone, we compared the experimental serial correlations with the ones obtained with the STPM, which does not assume any correlations between ISIs. In the example shown in Figure 3E the serial correlations are indeed accounted for by the STPM model meaning that the spiking history before the last spike does not affect the response.

In 12 out of 17 neurons the experimental and model serial correlations were not significantly different (two-sided bootstrap test, p > 0.01; Fig. 3F) confirming that for most neurons the spiking memory did not extend over the last spike. In three neurons the coefficient could not be determined because of a low number of triplets identified in responses. In two neurons the correlation coefficients were larger in the data than in the fitted STPM model (bootstrap test, p < 0.01).

We also compared the STPM with a GLM (Fig. 2C), which can account for spiking history extending over the last spike. The GLM showed a similar power in explaining both the average PSTH compared to the STPM (t test, p < 0.01; Fig. 2D, right box). However, it allowed for using larger bins without significant loss of goodness-of-fit (Fig. 2E,F). Finally, the introduction of spike history effects extending over multiple preceding spikes did not explain the significant serial correlations in every neuron. The GLM could account for the measured serial correlations in 13 out of 17 neurons. Overall, these results show that refractoriness is sufficient to explain the statistics in the within burst intervals obtained in most recorded neurons.

Trial-to-trial input variations induce significant serial correlations

The significant serial correlations found in two neurons could result from variability of the inputs that they receive. Although the peripheral stimulation of the median nerve used to evoke the somatosensory responses was well controlled over the duration of the recording, it is possible that the effective input to the cortex was modulated at the early stages of somatosensory pathway (cuneate nucleus, thalamus) and by on-going activity in the cortex. On the other hand, the STPM was fitted under the assumption that the inputs and model parameters do no change in time, i.e., that they are stationary.

To test the effects of trial-to-trial variability on the estimated STPM parameters and the serial correlations, we simulated spike trains from the STPM with a step-like recovery function and an exponentially decaying intensity function (Fig. 4B,C, dashed lines). In addition, in each trial we modulated the amplitude of the intensity function by a multiplicative gain, G, which was randomly drawn from uniform distribution on the interval [0.2, 1.8] (Fig. 4A). Next, we fitted the simulated surrogate data with an STPM assuming that the intensity function was fixed and that the trial-to-trial variability resulted solely from the probabilistic nature of the model. The fitted intensity function (Fig. 4B, red line) reflected the rapid onset and slower decay of the input after the stimulus, but its trace deviated from the “ground-truth” intensity function used in the simulation (compare the solid red and dashed gray curves; Fig. 4B). Importantly, the intensity function contained small ripples akin to the ones visible in the intensity function fitted to experimental data (Fig. 2A, left). Similarly, the fitted recovery function did not capture the step-like transition from refractoriness to baseline, but it manifested a prominent overshoot following the absolute refractory period and slower decay to baseline (Fig. 4C); such time-dependence was reminiscent of the shape of recovery function estimated from the data (Fig. 2A, right).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

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).

We also studied the effects of the gain modulation on the GLM. The intensity function estimated with this model still contained fluctuations absent in the function used for simulation, but their amplitude was reduced compared to the STPM. A greater improvement was observed in the GLM estimate of the recovery function, which approximated well the real function without a visible overshoot. Overall, both STPM and GLM misestimated some model parameters in presence of trial-to-trial variation, but we found that the GLM was more robust (Fig. 4D).

Finally, we estimated the serial correlation between the ISIs in presence of the input modulation. We found that the serial correlation was significantly larger compared to the spike trains simulated with the STPM with no trial-to-trial variations (Fig. 4E). This result shows that positive serial correlations can be obtained when neuronal responses vary from trial to trial reflecting changing inputs or excitability of the neuron. Since in our analysis in Figure 3E,F we compared experimental serial correlations to the ones obtained from the STPM, which does not account for the input variability, our estimate of serial correlations could reflect input modulation.

In summary, we show that the trial-to-trial variations of the input can explain several aspects of the STPM fitted to experimental data, in particular the ripples in the fitted intensity functions, overshoot following the refractoriness in the recovery function, and significant correlations between consecutive ISIs.

Trial-to-trial input variations induce correlations between single-neuron and population responses

Simultaneous recordings of single-neuron spike patterns and macroscopic EEG signals recorded from the surface of dura (high-frequency, >400 Hz, epidural hf-EEG) have shown that the spike patterns are not private to each neuron but that they are coordinated across a population of neurons responding to peripheral stimulation (Telenczuk et al., 2011). Such a coordination could possibly by achieved with a millisecond range-synchronization of the neurons, but the mechanisms of such a synchronization are not clear. Alternatively, it could be produced by the shared modulation of inputs or excitability. To test the latter hypothesis, we applied our probabilistic single-neuron model, the STPM, to a population of neurons receiving common gain modulation (Fig. 5A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

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).

As before, we assumed that the gain varies from trial to trial due to fluctuations in excitability, synaptic strength, or background activity. To investigate the effect of the gain factor on the population response, we simulated 5000 identical, statistically independent model neurons with the parameters estimated from the experimental data. We found that the frequencies of individual spike patterns depended on the value of the gain factor: Some patterns (for example 100) occur more frequently at low gain (G = 0.8), while others (for example 110) tend to occur more often at high gain (G = 1.2; Fig. 5B). Concurrently, the amplitude of the binned spike trains averaged across neurons (population PSTH) increased with the gain (Fig. 5C).

The concurrent dependence of population PSTH and single-neuron spike pattern distribution on the gain factor may also explain the correlation between single-neuron responses and macroscopic population activity found in experimental data. Spike patterns that are more frequent at low gain coincide predominately with a low-amplitude population PSTH whereas spike patterns elicited more frequently at high input gain, on average, coincide more often with a high-amplitude population PSTH. Consequently, the amplitude of the population PSTH might co-vary with single-neuron spike patterns. In particular, we found that the RMS amplitude of the high-pass filtered PSTH (>400 Hz) depends on the spike pattern used for grouping the trials (Fig. 5D).

To test whether gain modulation could explain the experimental results, we simulated the STPM model with trial-varying gain factor (see Methods) and compared the spike-pattern-specific hf-EEG amplitude calculated from experimental data with the simulated population response. We found that already for a modest level of the gain modulation (modulation strength γ = 0.2) the RMS amplitudes of the experimental hf-EEG and high-pass-filtered population PSTH of the model were strongly correlated (an example for one neuron, Pearson’s r = 0.93; Fig. 5E).

We found a positive correlation coefficient in 12 of 16 neurons that produced at least three different patterns. This fraction is significantly above the chance level expected from uncorrelated quantities (two-sided binomial test, p < 0.05). Thus, we conclude that gain modulation can introduce correlations between the single-neuron spike patterns and macroscopic population responses.

Spike patterns emerge as input-driven phenomena in a simplified biophysical model of a cortical neuron

The probabilistic models presented so far are abstract, and their parameters (intensity and recovery functions) cannot be linked directly to biophysical properties of a neuron. To interpret the generation of spike patterns mechanistically, we developed a simplified biophysical model of a cortical neuron based on the LIF model. Although this model does not reproduce faithfully all biological properties of realistic neurons, it captures their integration and spike generation properties, which are essential to the responses analyzed here. We simulated the neuron with two types of synaptic inputs, tonic excitatory and inhibitory inputs, and phasic thalamic excitatory inputs representing the barrage of action potentials triggered by peripheral stimulation.

In absence of thalamic inputs the model neuron elicits only few spikes due to spontaneous threshold crossings. However, in the model the median nerve stimulation is assumed to activate the thalamocortical fibers (28 synapses per cortical neuron), which then fire randomly according to a Poisson distribution with the rate of 700 spikes per second. These massive inputs trigger excitatory postsynaptic currents bringing the membrane quickly to the threshold. This results in a series of spike emissions accompanied by rapid successions of membrane de- and repolarizations (Fig. 6A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

A LIF model produces variable spike patterns. A, Sample traces of the membrane potential V_m(t) of a LIF model (see Methods for details) for three repetitions of the simulation. The ticks mark the threshold crossings, which lead to spike emission (color matched to the color of V_m trace). B, Poststimulus time histogram (firing rate) of spike trains obtained from 500 repetitions of the simulation. Vertical dashed lines delineate the events used for spike pattern analysis in D. C, Spike raster from all repetitions of the simulation. The “stim” arrow denotes the onset of the simulated thalamic inputs. D, Distribution of spike patterns obtained in the simulation of a LIF neuron (inset, PSTH). E, Distribution od spike patterns and PSTH (inset) for a model with modified parameters. Increasing the presynaptic firing rates of intracortical connections leads to higher coincidence of 101 and 011 patterns. In A-D the following parameters were used: w_exc = 0.0072 μs, τ_exc = 0.9 ms, f_exc = 10 Hz, w_inh = 0.02252 μs, w_Th = 0.035 μs, n_Th = 28, f_Th = 700 Hz, U = 0.65. E, Four parameters were modified from this baseline: f_exc = 30 Hz, w_Th = 0.05 μs, f_Th = 300 Hz, U = 0.7. All definitions and values of the remaining parameters are listed in Table 1.

We calculated the PSTH of the model by summing spike responses of n = 500 repetitions of the simulation. In each repetition the intracortical excitatory and inhibitory inputs, as well as the thalamocortical inputs, were drawn randomly from the Poisson distribution. Despite this randomness, the model PSTH is composed of discrete peaks well separated by short valleys showing that the neuron fired precisely at preferred latencies (Fig. 6B). Although the LIF model does not contain an explicit refractoriness, the intervals between the PSTH peaks correspond to the time required to depolarize the membrane from the reset potential (V_reset = –70 mV) to the spiking threshold (V_thr = –40 mV). In Figure 6A, this time is seen as the slow rise time following the rapid downstrokes (reset) of the membrane potential triggered by spikes. Such a hyperpolarized period acts effectively as the refractory period as seen in the STPM and GLM.

The characteristic decay of the response in the somatosensory cortex observed in the experimental data (Fig. 1B) could be driven by the adaptation of the neuron to the intense stimulation either at synaptic (Markram and Tsodyks, 1996) or neuronal level (Benda and Herz, 2003). Here, we model this process by means of short-term synaptic depression, which reflects the depression of thalamocortical synapses due to prolonged activity (Gil et al., 1997). The gradual decrease of synaptic drive makes the subsequent peaks smaller, broader, and separated by longer intervals (Fig. 6B) as observed also in the experimental PSTH (Fig. 1B). After 10 ms of stimulation, the thalamocortical synapses deplete, abolishing further discharges.

In practice, the inputs to somatosensory cortex can also decay after a brief median nerve simulation (0.2 ms) applied to the median nerve invalidating our assumption of sustained synaptic drive. However, it has been found that the thalamocortical projections can respond with prolonged firing to brief presentation of the stimulus (Swadlow and Gusev, 2001). Interestingly, such responses also formed bursts of action potentials. If the axonal delays of multiple thalamocortical neurons are matched at the submillisecond level, such bursts could provide oscillatory inputs cortical level. The effects of such input patterns on the cortical responses should be investigated in the future.

The responses of the LIF model neuron vary from trial to trial (Fig. 6C,D). This variability results from random cortical and thalamocortical inputs, which provide Poisson-distributed input spikes. Increasing the rate of inhibitory and excitatory inputs in a balanced fashion puts the neuron in a so-called noise-driven regime in which spikes are evoked by the random fluctuations over the threshold rather than by the mean depolarization (Destexhe et al., 2003; Zerlaut et al., 2016). In this regime, the responses of the neuron are more variable such that a broader range of different spike patterns is obtained across the trials. In particular, we found that the patterns with long latencies (such as, 011) or spike omissions (101) became more frequent at higher intracortical firing rates (Fig. 6E).

In summary, the LIF model indicates that bursts in the somatosensory cortex can be driven by the input and do not always require intrinsic bursting mechanisms (reviewed by Krahe and Gabbiani, 2004). The number of spikes per burst and the within-burst intervals can be mechanistically explained by the integrating properties of single neurons that are equipped with an intrinsic adaptation process or driven by synapses that show short-term depression. Strong thalamic inputs can produce precise population responses at preferred latencies, which can overcome the variability. At the single-trial level the variability of the thalamic input is expressed in the form of stereotyped spike patterns.