Encoding variable cortical states with short-term spike patterns

Neurons in the primary somatosensory cortex (S1) generate synchronised high-frequency (∼600 Hz) bursts in response to peripheral stimulation, and single-cell activity is locked to themacroscopic 600 Hz EEG wavelets. The mechanism of burst generation and synchronisation in S1 is not yet understood. Using a Poisson model with refractoriness that was fitted to to unit recordings from macaque monkeys, we can explain the synchronisation of neurons as the consequence of coincident synaptic inputs, while their high firing precision stems from the large input amplitude combined with a refractory mechanism. This model reproduced also the distribution of temporal spike patterns over repeated presentation of the same stimulus. In addition, the fine temporal details of the spike patterns are representative of the trial-to-trial variations in population excitability and bear upon the mean population activity. The findings are confirmed in a more detailed computational model of a neuron receiving cortical and thalamic inputs through depressing synapses. Our findings show that a simple feedforward processing of peripheral inputs could give rise to neuronal responses with non-trivial temporal and population statistics. We conclude that neural systems could use refractoriness to encode variable cortical states into stereotypical short-term spike patterns amenable to processing at neuronal time scales.


Introduction
Neurons generate highly variable responses to repeated presentation of the same stimulus. This variability might originate from thermal noise in ion channels [Chow andWhite, 1996, Schneidman et al., 1998], recurrent activity in the network [Destexhe et al., 2003, van Vreeswijk andSompolinsky, 1996] or modulations of neuronal excitability [Destexhe et al., 2001, Fontanini and Katz, 2008, Faisal et al., 2008. Over recent years many results have shown that a significant fraction of this variability is shared across large populations of neurons. These shared fluctuations were attributed to the variations of incoming stimuli and modulations of excitability [Brody, 1999, Shadlen and Newsome, 1998, Goris et al., 2014, Ecker et al., 2014. However, most of these studies focused on the spike rate variations over relatively long time scales neglecting millisecond-range spike timing differences. Such short time scales might be especially important for neurons that fire brief bursts of spikes at a frequency reaching several hundred spikes per second separated by longer silences [Evarts, 1964, Llinás and Jahnsen, 1982, Krahe and Gabbiani, 2004. Since the transitions between bursting and tonic firing characterised by long inter-spike intervals is dynamically controlled [Swadlow and Gusev, 2001] neuronal variability at both time scales might be relevant for neuronal processing.
Neurons in somatosensory cortices can encode their sensory inputs in the precise lengths (< 10 ms) of inter-spike intervals [Estebanez et al., 2012, Panzeri et al., 2001, Witham and Baker, 2015, so high firing precision is important for the reliability of stimulus encoding. In the primary somatosensory cortex (S1) of macaques single neurons respond to peripheral stimulation with barrages of spikes elicited at sub-millisecond precision [Baker et al., 2003]. However, when presented repetitively with the same stimulus they produce variable responses in terms of number of elicited spikes and lengths of inter-spike intervals, which might limit the amount of information they can carry. Despite the apparent variability their firing precision might be conserved if the single-trial responses are partitioned into spike patterns that reoccur in multiple trials [Toups et al., 2012]. Such spike patterns have been indeed observed in the S1 [Telenczuk et al., 2011], but neither the mechanism of their generation nor their functional significance has been identified.
We propose a mechanism in which the precise patterns of single-neuron responses stem from the interplay between synaptic inputs and intrinsic refractory properties of the cell [Berry andMeister, 1998, Czanner et al., 2015]. To test this hypothesis, we develop simple models capturing the two processes and fit the parameters of the model to extracellular recordings of single-unit activity in the somatosensory cortex [Baker et al., 2003].

Experimental methods
We reanalysed data from [Baker et al., 2003] in which the experimental protocol is described in more detail. Briefly, neuronal responses were evoked in the hand representation of the primary somatosensory cortex of two awake Maccaca mulatta monkeys by electrical median nerve stimulation at the wrist (pulse width: 0.2 ms; repetition rate: 3 Hz; intensity: 150% motor threshold). Single-unit activity was recorded extracellularly using a 16-channel Eckhorn drive [Thomas Recording GmbH;Giessen, Germany;Eckhorn and Thomas, 1993]. Each of the platinum/glass electrodes (electrode impedance: 1 MΩ) was advanced into cortex (area 3b) until well-isolated neurons were found with one of the electrodes. The receptive fields of the cells were tested by means of manual tapping using a stylus.
Spike waveforms were first band-pass filtered (1 kHz -10 kHz) and then sampled with a frequency of 20 kHz. Action potentials of neurons surrounding the microelectrode were detected in the extracellular recordings by means of amplitude thresholding; the threshold was chosen manually to detect spikes whose amplitude was significantly above noise level. The wave shapes of the detected action potentials were parametrised by their amplitude, width and projection coefficients on two main principal components. The spike timings of single units were determined based on these shape features using a manual cluster cutting method that allowed for identification of clusters of arbitrary shapes [Lewicki, 1998, Hazan et al., 2006. To ensure correct clustering the procedure was performed by two operators using different software packages (GetSpike, S. N. Baker; PySpikeSort, B. Telenczuk) and then checked for consistency.
In order to validate the spike discrimination we checked the extracellular action potentials generated by a putative single cell for the consistency of the wave shape and amplitude. Additionally, we searched for interspike intervals (ISIs) shorter than 1 ms; if such short intervals were found the clustering procedure was repeated. Spike trains with evidence of poor spike sorting (inconsistent wave shapes or ISIs < 1 ms) were excluded from subsequent analysis.

Spike patterns classification
From 46 cells identified in both monkeys we selected 17 cells that responded with bursts of spikes defined by: a response with more than one spike following at least 4% of stimuli, and a mode of the interspike interval histogram shorter than 1.8 ms [Baker et al., 2003].
Among these 17 neurons we identified neurons that also fired spontaneous bursts by counting the number of inter-spike intervals in the post-stimulus period (> 30 ms after the stimulus) that were shorter than 1.8 ms. In this time window the initial response dies out and baseline firing rate is re-established. Neurons that fired at least 10% of bursts in this window were labelled as spontaneous burstsers.
In each cell we summed spikes over all trials and identified prominent peaks in the obtained peri-stimulus time histograms (PSTH; bin width 0.2 ms, Figure 2A). As the within-burst spike composition varied from trial to trial, each trial was described with a binary string whose entries (one or zero) represented the occurrence or non-occurrence of a spike in a sequence of bins bracketing the major peaks of the overall PSTH: the borders between the bins were placed manually in the troughs of the PSTH (Figure 2A, B: vertical lines). Each string corresponded to one spike pattern; the length of the string equalled the total number of peaks in the PSTH.

Spike-train probability model
To reproduce the distribution of emitted spikes in a single neuron we chose a minimal model (spike-train probability model, STPM) replicating high variability in the cortical responses [Softky andKoch, 1993, Destexhe et al., 2001] and manifesting refractoriness (decreased probability of spiking for some time after producing a spike).
We assume that a spike emission is a random point process with the probability: where {t i } denotes the spiking history and λ(t|{t i }) is the conditional intensity.
The conditional intensity λ(t|{t i }) is assumed to have a Markov property i.e. it is conditioned only on the time of occurrence of the last spike t last : λ(t|{t i }) = λ(t, t last ). A further assumption is that the firing rate modulation and refractory effects are multiplicative [Berry and Meister, 1998]: where q(t) is the intensity function, w(∆t) is the recovery function and ∆t = t − t last is the time from the last spike.
The parameters of the model, the intensity function q(t) and the recovery function w(∆t), are defined on a per-bin basis and fitted to experimental data by means of the maximum likelihood approach. To capture fine temporal details of the neuronal responses (response onset and interspike intervals) the intensity and recovery functions were defined with short sampling interval (0.05 ms). The log-likelihood function L(q; w|{t i }) is obtained by log-transforming the probability function of an inhomogeneous Poisson process with the conditional intensity (2) [Dayan andAbbott, 2001, Johnson andSwami, 1983]: where T is the duration of response (T = 30 ms), i is the spike index and t i denotes the times of occurrence of recorded spikes spikes. The likelihood function is maximised by means of an iterative expectation-maximization (EM) algorithm, which guarantees that the global maximum is reached [Miller, 1985]. In addition, we ensure that after 5 ms the model neuron recovers from refractoriness by setting the recovery function to unity for long intervals, such that w(∆t > 5 ms) = 1.
To study the effects of refractoriness on the modelled responses we compared the results to the inhomogeneous Poisson model without refractory period [w(∆t) = 1 for all ∆t > 0; Dayan and Abbott, 2001]. The model is fully characterised by its intensity function q poisson (t), which can be estimated directly from the experimental spike trains with the post-stimulus time histogram (PSTH, bin width: 0.05 ms).

Model validation
To validate the model, the dataset was divided into two non-overlapping subsets of equal size: the training and validation set. The trials for each set were selected randomly from all stimulation repetitions. The parameters of the model were fitted to the training set. Based on these parameters 1000 spike trains were simulated. The goodness-of-fit was evaluated based on the difference between one of two statistics (PSTH, bin size 0.2 ms or spike pattern distribution) calculated separately for the simulated and validation spike trains. The difference was quantified by the sum of squared and normalised residuals ("model error", χ 2 statistics) [cf. Rauch et al., 2003]: where X denotes the statistics vector, either the PSTH or spike pattern distribution; N is the size of the vector and equals the number of bins (N = 150) or the number of identified spike patterns (N ≤ 16); and X model and X validate denote the statistics calculated from the modelled spike trains and validation spike trains, respectively.
The model error Err(X model , X validate ) was compared against the error between the training and validation set Err(X train , X validate ) ("reference error"). The significance of the difference between these two measures was quantified by means of the F-test with (N − 1, N − 1) degrees of freedom [Barlow, 1989], where:

Serial correlations
From the responses of single cells (response window t ∈ [0, 50] ms after stimulus) we identified spike triplets defined as three consecutive spikes separated by intervals shorter than 4 ms. Next, we calculated Pearson's correlation between the interspike intervals (ISIs) of the first and the second spike and the second and the third spike in the triplet (r data ). We compared the estimated r data to the correlation coefficient calculated from surrogate data simulated with the STPM model whose parameters were fitted to the experimental spike trains (r model , same number of trials). The significance of the differences between correlation coefficients found in simulated and experimental ISIs was tested by means of a bootstrap test. To this end, r model was calculated from 1000 independently simulated datasets and the resulting coefficients were compared to r data . The fraction of times for which r model was greater (or smaller) than r data was taken as the p value of the test (two-sided bootstrap test).

Population model
To model the relation between spike patterns of a single neuron and the response of the population we extended the model to a population of uncoupled neurons receiving common inputs.
Population response was calculated from a simulated ensemble of 5000 identical neurons.
The parameters of the spike-train probability model were fitted to the responses of the analyzed cell and shared by all model neurons. In each trial j the intensity function of all neurons was modulated by a multiplicative factor (gain, G j ): q j (t) = G j q(t), where gain G j was drawn from uniform distribution on the interval [0.8, 1.2]. From the obtained single-trial single-neuron responses the total population response was calculated by summing the binned spike responses of all cells (bin size 0.2 ms) and subsequent band-pass filtering (400 -1200 Hz).
Population correlates of single-neuron spike patterns were identified by averaging the population responses sorted according to the concomitant spike patterns of a single cell randomly selected from the ensemble. The spike pattern classification was based on the occurrence/omission of spikes in a discrete sequence of spiking "windows" as described above. The root mean square (RMS) amplitude of the pattern-specific average was compared with the experimentally-obtained high-frequency EEG (>400 Hz, hf-EEG) related to the same spike pattern of a recorded cell [cf. Telenczuk et al., 2011]. The similarity of the values across different spike patterns was quantified by means of Pearson's correlation coefficient.

Biophysical model
In order to understand the mechanisms of burst generation we developed a simplified singleneuron model. The model consists of a linear neuron with spiking threshold (leaky integrate-andfire), which receives conductance-based inputs through depressing synapses (short-term synaptic depression). The details of the model and its parameters are given in Supplementary Information.

Results
Neurons in area 3b of macaque monkeys respond with brief (< 10 ms), highly variable bursts of activity in response to stereotypical electrical stimulation of the median nerve ( Figure 1A-B).
To understand mechanisms underlying such response variability, we propose a phenomenological model of bursting of cortical neurons. The model is based on two experimental observations: 1) upon presentation of strong sensory stimuli, layer IV cortical neurons are bombarded with intense and coincident synaptic inputs from thalamocortical neurons [Bruno and Sakmann, 2006, Hanajima et al., 2004, Swadlow and Gusev, 2001, Gil et al., 1999, Cruikshank et al., 2007; 2) after emitting a spike, neurons are refractory, which limits their maximum firing rate [Berry and Meister, 1998, Gray, 1967, Kara et al., 2000. In order to illustrate the effects of these two phenomena on neuronal responses, we simulated a probabilistic model (the spike-train probability model, STPM, see Methods) with an exponentially decaying intensity function and an absolute refractory period τ ref = 1.4 ms ( Figure 1C). The post-stimulus time histogram (PSTH) of the simulated spike responses reproduces main features of the PSTH obtained from experimental data. Specifically, the absolute refractory period leads to an appearance of multiple peaks in the PSTH that are typical for cortical burst responses triggered by peripheral nerve stimulation (compare with Figure 1B, bottom). Moreover, in agreement with experimental findings the first peak of the PSTH is narrowest, reflecting the precise spike timing that is induced by the high amplitude of the intensity function at the response onset. The subsequent decay of the function leads to a gradual widening of the peaks and decrease of their amplitude.
Refractoriness explains the intra-burst intervals. In order to test whether the STPM can reproduce the fine details of neuronal responses recorded in vivo, we fitted it to the spike trains produced by neurons in the somatosensory cortex. Although the fitted parameters, i.e. the intensity function and recovery function ( Figure 1D) do not have a straightforward biophysical interpretation, they could be related to the synaptic input of the neuron and refractoriness [Berry andMeister, 1998, Johnson andSwami, 1983]. When compared to the observed firing rate (PSTH, Figure 1E, dark blue), the prominent dips of spike emission probability are largely diminished ( Figure 1D, left). Although the maximum of the function is much above the rate at which individual neurons can fire spikes, the maximum firing rate of the model neuron is limited by refractoriness as quantified by the recovery function ( Figure 1D, right). In agreement with the properties of biological neurons, the 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, largely exceeding the rest value.
This phenomenon may be interpreted as a facilitation of a spike initiation, which has been found in neuronal membranes exhibiting subthreshold resonance [Engel et al., 2008, Hutcheon et al., 1996, Verechtchaguina et al., 2006, or artifact introduced by the trial-to-trial response variations (see the population model below and Supplementary Figure 2). The fast fluctuations following the facilitation peak represent statistical noise due to finite size of the data set. Altogether, a simple Poisson model with refracoriness (STPM) can explain the structure of PSTHs of neurons in S1 in response to median nerve stimulation.
The spike-train probability model provides a parsimonious description of bursting in S1 cortex. The simulated peri-stimulus time histogram ( Figure 1E,"Model") matches closely the one obtained from the experimental data ( Figure 1E, "Training"). In order to check whether 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. We found that the difference of the fitted PSTH from the validation set was of the same magnitude that the variation within the dataset (see Methods, F-test, p < 0.05). This 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 available cells yielding similar results.
Importantly, an application of the cross-validation procedure revealed that in 9 out of 17 cells the PSTH simulated with the model was not significantly different from the PSTH calculated from the recorded spike trains (F-test, F = 0.59 − 1.3, p < 0.05). In 4 cells, only small deviations of the modelled PSTH from the real PSTH were found (F = 1.5 − 2) and in the remaining 4 cells the deviations were large (F = 2.96 − 9.1).
To further analyse the cases in which the simulated spike trains diverged from the data, we subtracted the model PSTH from the validation data PSTH ( Figure 1E, bottom). The resulting residuals still contained some minor fluctuations aligned to the peaks of the PSTH. This might be a sign that the model does not capture completely the shape of the PSTH. Indeed, the correlation coefficient between the residuals and the PSTH of validation set was significantly positive (bootstrapped 95% confidence intervals, Figure 1F) meaning that the residuals contain the remnants of the PSTH. However, the correlation coefficient was substantially decreased compared to the model with no refractoriness ( Figure 1F, middle box).
Poisson-like variability explains occurrence of temporal spike patterns in repeated trials. Having shown that Poisson-like firing statistics can account for large part of the trialaveraged response of single neuron, we tested whether it can explain the trial-to-trial variability.
In order to quantify the trial-to-trial variability of neuronal responses, we sorted single-trial spike trains according to the occurrence of spikes in pre-defined temporal windows (Figure 2A-C).
Each spike train was assigned a binary word based on occupancy of "preferred firing windows" aligned to the peaks of the PSTH. We call this string of spikes and silences a spike pattern. The distribution of spike patterns in the experimental data was similar to the distribution from the spike-train probability model (STPM, Figure 2D). In contrast, spike patterns produced without a refractory period appeared at frequencies much different from the experimental data, despite the fact that the overall PSTHs were almost identical ( Figure 2D, inset).
To quantify the similarity between the experimental and modelled spike patterns, we used a cross-validated chi-square statistics (see equation (4) in the Methods). In 11 of 17 examined cells the spike-pattern distribution of the model was not significantly different from experimental distribution (F-test, p < 0.05).
Within-burst intervals manifest significant correlations. Next, we investigated whether the correlations between consecutive inter-spike intervals (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 from the last spikes. However, the calculation of the serial correlations in the experimental data obtained from S1 showed that two consecutive inter-spike intervals are not independent ( Figure 2E). Since significant serial correlations might be induced also 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 interspike intervals ( Figure 2E). In 8 out of 17 cells the experimental and model serial correlations were not significantly different (two-sided bootstrap test, Figure 2F) confirming that the neuron spiking memory did not extend over the last spike. In 3 cells the coefficient could not be determined because of a low number of triplets identified in responses. In the remaining 6 cells the correlation coefficients were larger in the data than in the fitted STPM model, which suggest that an adaptation processes or synaptic plasticity could shape the recorded responses, which was not included in our model so far.
To check whether such processes operating at longer time scales might affect the spike responses, we also compared the STPM with the generalised linear model (GLM, for details see Supplementary Information). Although the GLM allowed for better control of the number of parameters, both models showed similar power in explaining both the average PSTH, serial correlations and distribution of spike patterns (Supplementary Figure 1). In particular, we found that the model could not account for the measured serial correlations in 4 out of 17 neurons.
Overall, these results show that the refractoriness can explain the correlations observed in the experimental spike patterns.
Trial-to-trial input variations induce correlations between single-neuron and population responses. Simultaneous recordings of single-neuron spike patterns and macroscopic population signals (high-frequency, > 400 Hz, epidural EEG) have shown that the spike patterns are not private to each neuron but they are coordinated across population of neurons responding to peripheral stimulation [Telenczuk et al., 2011]. Such coordination would require millisecond range synchronisation of the neurons, but the mechanisms with which it could be achieved are not clear. Alternatively, it could be produced by the shared modulations of inputs or excitability. In order to test the latter hypothesis, we extended our probabilistic single-neuron model, STPM, to a population of neurons receiving common inputs ( Figure 3A). These inputs might reflect, for example, feed-forward inputs from lower stages of the processing pathway [Shadlen and Newsome, 1998], mean activity in a recurrent network [Shu et al., 2003] or neuromodulatory signals [Gil et al., 1997].
We assume that the effective input arriving to the neurons varies from trial to trial due to fluctuations in excitability, synaptic strength or total synaptic activity. In the model, these variations translate to a multiplicative gain factor, G, that acts on the intensity function (see Methods). In order 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, while others (for example "110") tend to occur more often at high gain ( Figure 3B).
Similarly, the amplitude of the binned spike trains averaged across neurons (population PSTH) increased with the gain ( Figure 3C).
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 will coincide predominately with a low-amplitude population PSTH, whereas those elicited more frequently at high input gain will, 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 ( Figure 3D).
To test whether gain modulations could explain the experimental results, we compared the spike-pattern-specific high-frequency EEG (hf-EEG) amplitude calculated from experimental data with the model predictions. We found that the root-mean-square amplitudes of the experimental hf-EEG and high-pass-filtered population PSTH of the model were strongly correlated (Pearson's r = 0.99, Figure 3E).
Overall, we found a positive correlation coefficient in 12 of 16 cells that produced at least 3 different patterns. This fraction is significantly above the chance level expected from uncorrelated quantities (two-tailed binomial test, p<0.05). Thus, we conclude that the gain modulations can introduce the 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) can not be linked directly to biophysical properties of a neuron. To interpret the generation of spike patterns mechanistically we developed a sim-plified biophysical model of a cortical neuron based on the leaky integrate-and-fire (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 analysed 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, massive inputs from the thalamus trigger excitatory post-synaptic currents which bring the membrane quickly to the threshold. This results in a series of spike emissions accompanied by rapid successions of membrane de-and re-polarisations ( Figure 4A).
We calculated the PSTH of the model by summing spike responses of n = 500 repetitions of the simulation ( Figure 4B). The model PSTH is composed of discrete peaks separated by short intervals. In spite of the random inputs to the neuron, these peaks are well separated showing that the neuron fired precisely at preferred latencies.
The characteristic decay of the response in the somatosensory cortex observed in the experimental data ( Figure 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 make the subsequent peaks smaller, broader and separated by longer intervals as observed in the experimental PSTH ( Figure 1B). After 10 ms of stimulation the thalamocortical synapses deplete completely abolishing further discharges.
The responses of the model neuron vary from trial to trial ( Figure 4C). This variability results from random cortical and thalamocortical inputs which provide Poisson-distributed input spikes.
The distribution of spike patterns over trials is well captured by the STPM, which means that the spike patterns can be well represented by a time-varying intensity and refractoriness ( Figure   4D).
This simple model shows that bursts in the somatosensory cortex can be input driven and do not require intrinsic bursting mechanisms (reviewed by Krahe and Gabbiani [2004]). The number of spike per burst and the within-burst intervals can be mechanistically explained by the integrating properties of single neurons coupled with adaptation process such as short-term synaptic depression. The strong thalamic inputs can overcome the variable discharges due to ongoing cortical activity producing precise populations response at preferred latencies. However, at single-trial level the variability is expressed in the form of stereotyped spike patterns.
By means of simplified phenomenological models and a point-neuron model we showed that within-burst variability of cortical S1 neurons can be decomposed into the private variability of each neuron and multiplicative input modulation that is shared by the entire population.
The private variability explains most of the differences between responses elicited in single trials and underlies the re-appearance of the same spike patterns over multiple repetitions of the stimulus. The shared gain modulation coordinates the responses of many responding neurons and explains the puzzling co-variability between single-neuron and macroscopic population responses demonstrated in experimental recordings [Telenczuk et al., 2011]. The models shed also light on the mechanism of S1 burst generation, their synchronisation across neurons and suggest that spike patterns may encode time-varying cortical state at fast temporal scales.

Mechanism of bursting
By means of a simple phenomenological model, the STPM, we showed that bursting in the primary somatosensory cortex results from the combination of intense synaptic bombardment and refractory period. Such fast bursting triggered and sustained by an intense synaptic input has been termed "forced bursting" [Izhikevich, 2006].
The model parameters reproducing S1 responses are consistent with the physiological properties of chattering neurons -a class of bursting neurons found in cortical layer III [Gray andMcCormick, 1996, Steriade et al., 1998]. The chattering neuron responses are generated by the interaction of the after-depolarising sodium current (ADP) and after-hyperpolarising pottassium current (AHP) [Brumberg et al., 2000]. The after-depolarising current is activated after each spike and produces a depolarisation on top of which another spike can be triggered. Similarly, the after-hyperpolarising current is triggered by a spike, but it produces a valley in the membrane potential suppressing further spikes.
The shape of fitted recovery function in both models agrees well with the contribution of the two currents: The initial dip, which we interpret as refractoriness, might reflect the AHP and the inactivation of sodium channels, whereas the subsequent over-shoot might correspond to the ADP. We note, however, that the over-shoot might be an artifact due to trial-to-trial variability (Supplementary Figure 2) and it was not necessary to produce bursting responses -in fact we demonstrated in a toy model that the absolute refractory period combined with intense but transient inputs is sufficient to produce bursts with similar (but not identical) statistics ( Figure   1C).
The STPM could also account for most of the correlations between inter-spike intervals (serial correlations). However, in a few cells we found serial correlations differing from the ones it predicted. Since in these neurons processes occurring at long time scales could shape the spike patterns, we fitted them with the GLM (Supplementary Materials), which considers spike history effects extending over multiple inter-spike intervals. We found that the GLM with a horizon of 8 ms provided an optimal fit to these data in agreement with the time scales of short-term synaptic plasticity [Tsodyks and Markram, 1997] and firing-rate adaptation [Benda and Herz, 2003]. The latter is often mediated by the slow AHP currents providing another link between a biophysical process and the recovery function of our phenomenological model.
We were able to reproduce qualitatively both the average and single-trial features of the burst responses in a more realistic leaky integrate-and-fire neuron. Although such models are gross simplification of the real neurons both in terms of spiking mechanism and morphological features, it has been suggested that LIF may faithfully reproduce some features of spike generation [Brette, 2015]. In the model, the within-burst interval was controlled by the time required to reach the threshold from the hyperpolarised state (membrane time constant) and the gradual decay of the amplitude of PSTH peaks was due to the short-term synaptic depression. The latter mechanism can be related to the depletion of the available vesicles in the pre-synaptic terminal [Markram and Tsodyks, 1996]. However, it would be possible to replace it with some other form of adaptation [Brette and Gerstner, 2005]. Both mechanism lead to extinction of the initial synaptic drive, which explains the burst-like transient response to the step-like thalamic inputs.
We note, however, that without the recordings from the thalamocortical projection neurons we can not infer the inputs of the cortical neurons. Our models are still compatible with temporally structured inputs.
The trial-to-trial variability of the model was due to variable arrival times of the thalamic inputs, but also due to the intra-cortical inputs. The latter were configured such that the neuron was in the "high-conductance state" reproducing the property of cortical neurons receiving constant bombardment of inhibitory and excitatory inputs [Destexhe et al., 2003]. Apart from decreasing the membrane time constant thus allowing for rapid repeated discharges, these intracortical inputs introduced substantial trial-to-trial variability that could explain the observed spike pattern distribution.
Previous studies have shown that most of the bursting neurons in the S1 macaque cortex are characterised with broad spikes, which suggests that they are pyramidal neurons or spiny stellate cells [Baker et al., 2003]. This is confirmed by intracellular recordings in barrel cortex showing that regular spiking pyramidal neurons but not intrinsic bursting neurons followed the phase of high-frequency oscillations in surface recordings [Jones et al., 2000]. Our results are consistent with these findings and strengthen the evidence that a subclass of S1 neurons activated by median nerve stimulation belongs to the chattering or regular spiking pyramidal neurons. However, a subset of neurons analysed here (5 of 17 neurons) did also fire spontaneous burst showing that at least some of them may belong to intrinsic bursting class.

Burst synchronisation
A striking feature of the S1 bursting is that the signature of the burst appears also in macroscopic EEG signal. The visibility of the burst in the surface recordings was interpreted as a sign of strong synchronisation between the neurons [Jones et al., 2000], which could be mediated, for example, by fast synaptic potentials or gap junctions [Draguhn et al., 1998]. By extending our model to a population of uncoupled neurons, we demonstrated that the sub-millisecond synchronisation between multiple neurons does not require a fast coupling mechanisms, but results from shared synaptic inputs arriving through thalamocortical fibers. Provided that the biophysical properties of the receiving population and axonal conduction times vary in a narrow range, these inputs will elicit synchronous bursts of spikes. The required precision in the arrivals of afferent spikes could be achieved by means of a plasticity rule that selects inputs arriving synchronously at the cortical synapses [Gerstner et al., 1996].

Role of spike patterns
Trial-to-trial variations in S1 responses can be classified into a set of spike patterns defined by the occurrences of spikes within 10-ms-long bursts [Telenczuk et al., 2011]. Such temporal patterns of neuronal responses were first identified in cat striate cortex and crayfish claw [Dayhoff and Gerstein, 1983], and later in the temporal cortex of monkeys, cat lateral geniculate nucleus [Fellous et al., 2004] and in the rat hippocampus Buzsaki, 2007, Schmidt et al., 2009].
We proposed a model in which the occurrence of spike patterns is regulated by the input intensity, that is the rate of incoming spikes and not their precise timing. The temporal information stored in the spike patterns is complementary to the output rate (spike count) in the sense that the spike patterns with identical number of spikes (and therefore the same output rate) could still provide extra information concerning its inputs. For example, the early (110) doublet is more common for high input intensity; the opposite is true for the late (011) doublet ( Figure 3B). This mechanism could be especially useful for encoding inputs that would normally exceed the maximum firing rate set by the refractory period.
In one study the stimulus intensity was related to the within-burst intervals of spike responses recorded in the dorsal lateral geniculate nucleus (dLGN) [Funke and Kerscher, 2000]. Our results are consistent with this hypothesis. In the STPM, the within-burst intervals are constrained by the refractory period, but their length can also vary as a function of the synaptic drive (intensity function). In addition, the length of refractory period may not be fixed but it might be modulated by the firing rate. It has been shown that models allowing for this modulation may better describe the spike times in response to the time-varying stimulation [Koyama and Kass, 2008].
Short trains of spikes are also well suited to evoke specific synaptic response or trigger synaptic plasticity [Lisman, 1997, Song et al., 2000, Swadlow and Gusev, 2001, Tsodyks and Markram, 1997, Maass and Zador, 1999, they are optimally placed to represent neuronal variables in a form that is easily processed, stored and transmitted [Leibold et al., 2008, Tiesinga et al., 2008. In this spike-timing-based view neural systems take advantage of the temporal information encoded into spike patterns to represent slowly-changing cortical states (such as attention or waking).
Alternatively, spike patterns could also allow for more reliable representation of neuronal inputs [Toups et al., 2012]. These rate-based and spike-timing-based interpretations of spike patterns are not contradictory and could even act as independent communication channels [Tiesinga et al., 2008].
We showed that the distribution of spike patterns over neurons and the amplitude of the averaged population signal are regulated by input magnitude, which could reflect gating of neuronal signals through attention, expectation, sleep and waking [Fontanini and Katz, 2008, Shu et al., 2003, Steriade et al., 2001. Similar gain control mechanisms were implemented in realistic neural models through, for example, concurrent modulation of excitation and inhibition [Chance et al., 2002, Vogels andAbbott, 2009] or short-term synaptic depression [Rothman et al., 2009]. More generally, multiplicative noise can account for the variability and co-variability of neuronal responses in the thalamus and many cortical areas, including the lateral geniculate nucleus, V1, V2 and MT [Goris et al., 2014].
Macroscopic signatures of the bursts were shown to match the somatosensory-evoked potentials in monkey epidural EEG and human scalp EEG, so the high-frequency EEG burst might link the non-invasive macroscopic recordings and microscopic neuronal activity [Curio, 2000].
Our modelling shows how the characteristic features of the underlying spike burst, i.e., its frequency and amplitude, can be related to the biophysical properties, such as refractory period, whereas the small spike timing differences could be under control of the dynamical cortical state. Single-cell responses averaged over all trials (PSTH, same as in Figure 1B) reveal that spikes occur preferentially at discrete latencies (delimited by vertical lines). (B) In single trials multiple spikes are elicited in diverse combinations of preferred latencies resulting in significant trial-totrial response variability. Spike combinations are classified into spike patterns: The time axis was first divided into three windows aligned to the peaks of the overall PSTH. Each trial was then assigned a binary string (spike pattern 'xyz', from '000' to '111'), where 1 represents he occurrence and 0 the absence of a spike in consecutive windows. 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 patterns distribution obtained from the data (white bars, same as C), spike-train probability model (STPM, red bars) and the STPM without refractory period (blue bars). The firing rate of the Poisson model was estimated by a PSTH with bin size 0.05 ms. (E) Scatter plot of two consecutive interspike intervals (ISIs) within spike triplets calculated from the experimental data (filled circles) and simulated responses (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 model (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 of the STPM fitted to experimental data provide the distribution of serial correlations consistent with the model (empty bars); experimental correlation coefficient (vertical arrow, r data ) is likely to be drawn from the estimated distribution (two-sided bootstrap test, p = 0.81).  Figure 1D) with inputs multiplied by a random gain factor G (uniform distribution between 0.8 and 1.2). Only a single neuron from the population is recorded (gray dot) and its spike pattern is identified. From the simulated spike trains the population PSTH was calculated and then high-pass filtered to obtain an estimate of the high-frequency EEG population response. (B) Distributions of spike patterns of a single cell in 1000 repetitions of the simulation with a low (0.8, blue) and a high gain (1.2, red). (C) The amplitude of the population PSTH before (top panel) and after high-pass filtering (bottom panel) varies with the gain. (D) The concomitant changes in the spike pattern distributions and high-pass-filtered population PSTH result in the correlation between single-neuron spike pattern and root-mean-square (RMS) amplitude of the high-pass filtered population PSTH. (E) The simulated population RMS amplitudes correlate with experimental hf-EEG RMS related to the same pattern (hf-EEG RMS) [Telenczuk et al., 2011].  Table 1.