Refractoriness Accounts for Variable Spike Burst Responses in Somatosensory Cortex

Abstract Neurons in the primary somatosensory cortex (S1) respond to peripheral stimulation with synchronized bursts of spikes, which lock to the macroscopic 600-Hz EEG waves. The mechanism of burst generation and synchronization in S1 is not yet understood. Using models of single-neuron responses fitted to unit recordings from macaque monkeys, we show that these synchronized bursts are the consequence of correlated synaptic inputs combined with a refractory mechanism. In the presence of noise these models reproduce also the observed trial-to-trial response variability, where individual bursts represent one of many stereotypical temporal spike patterns. When additional slower and global excitability fluctuations are introduced the single-neuron spike patterns are correlated with the population activity, as demonstrated in experimental data. The underlying biophysical mechanism of S1 responses involves thalamic inputs arriving through depressing synapses to cortical neurons in a high-conductance state. Our findings show that a simple feedforward processing of peripheral inputs could give rise to neuronal responses with nontrivial 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 (tens of milliseconds).


Introduction
Neurons usually generate highly variable responses to repeated presentations of the same stimulus. This vari-ability might originate from thermal noise in ion channels (Chow and White, 1996;Schneidman et al., 1998), recurrent activity in the network (van Vreeswijk and Sompolin-

Significance Statement
Neurons in the hand area of the primary somatosensory cortex respond to repeated presentation of the same stimulus with variable sequences of spikes, which can be grouped into distinct temporal spike patterns. In a simplified model, we show that such spike patterns are product of synaptic inputs and intrinsic neural properties. This model can reproduce both single-neuron and population responses only when a private variability in each neuron is combined with a multiplicative gain shared over whole population, which fluctuates over trials and might represent the dynamical state of the early stages of sensory processing. This phenomenon exemplifies a general mechanism of transforming the ensemble cortical states into precise temporal spike patterns at the level of single neurons. sky, 1996;Destexhe et al., 2003) or modulation of neuronal excitability (Destexhe et al., 2001;Faisal et al., 2008;Fontanini and Katz, 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 modulation of excitability (Shadlen and Newsome, 1998;Brody, 1999;Ecker et al., 2014;Goris et al., 2014). However, most of these studies focused on spikerate variations over long time scales, neglecting millisecondrange 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 much longer intervals of silences (Evarts, 1964;Llinás and Jahnsen, 1982;Krahe and Gabbiani, 2004). Since the transitions between bursting and tonic firing characterized by longer interspike intervals are dynamically controlled (Swadlow and Gusev, 2001) 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 interspike intervals (Panzeri et al., 2001;Estebanez et al., 2012;Witham and Baker, 2015), which suggests that high firing precision is important for the reliability of stimulus encoding. In the primary somatosensory cortex (S1) of macaque, single neurons respond to peripheral stimulation with barrages of spikes elicited at sub-millisecond precision (Baker et al., 2003;Fig. 2). However, when presented repetitively, the same stimulus produces variable responses in terms of the number of elicited spikes and the lengths of interspike intervals, which might limit the amount of information they can carry. It is, however, possible that such trial-to-trial variability represents an alternation between several classes of reliable responses, called spike patterns (Toups et al., 2012). Such spike patterns have been indeed observed in S1 (Telenczuk et al., 2011), but neither the mechanism of their generation nor their functional significance has been identified.
Here, we propose a mechanism that explains the precise patterns of single-neuron responses as an interplay between synaptic inputs and intrinsic refractory properties of the neuron (Berry and Meister, 1998;Czanner et al., 2015). To test this hypothesis, we develop simple models capturing the two processes, and we are able to fit the parameters of the models to extracellular recordings of single-unit activity in the somatosensory cortex.

Experimental methods
Neuronal responses were evoked in the hand representation of the S1 of two awake Maccaca mulatta female monkeys by electrical median nerve stimulation at the wrist (pulse width: 0.2 ms; repetition rate: 3 Hz; intensity: 150% motor threshold; see also Fig. 1A). Single-unit activity was recorded extracellularly using a 16-channel Eckhorn drive (Thomas Recording GmbH;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 neurons were tested by means of manual tapping using a stylus.
In addition, we recorded EEG signals from the surface of the dura (epidural EEG) with two electrodes placed in the vicinity of the micro-electrode array. The signals were then high-pass filtered (Ͼ400 Hz) to obtain the highfrequency EEG (hf-EEG).
All experimental procedures were performed according to Home Office UK (Scientific Procedures) Act 1986 regulations and institutional ethical guidelines.

Spike sorting
From the extracellular recording we obtained spike waveforms that were first bandpass filtered (1-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 parametrized 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 and then checked for consistency.
To validate the spike discrimination, we checked the extracellular action potentials generated by a putative single neuron 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 analyses.

Spike pattern classification
From 46 neurons identified in the two monkeys, we selected 17 neurons that responded with bursts of spikes. Bursting neurons were defined by responses with more than one spike for at least 4% of stimuli and a mode of the ISI 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 ISIs in the poststimulus 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 reestablished. Neurons that fired at least 10% of bursts in this window were labeled as spontaneous bursters.
In each neuron we summed spikes over all trials, and we identified prominent peaks in the obtained poststimulus time histograms (PSTH; bin width 0.2 ms; Fig. 3A). As the within-burst spike composition varied from trial to trial, each trial was described with a binary string whose entries remarks on earlier version of this paper.
(one or zero) represented the occurrence or nonoccurrence 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 (Fig. 3A,B, vertical lines). Each binary string corresponded to one spike pattern; the length of the string equalled the total number of peaks in the PSTH.
In addition, we averaged the concomitant hf-EEG responses over trials concurring to each of the identified spike patterns of a single neuron.

Spike-train probability model (STPM)
To reproduce the distribution of emitted spikes in a single neuron, we chose a minimal model (STPM) that could replicate the observed high variability in the cortical responses (Softky and Koch, 1993;Destexhe et al., 2001) and manifest refractoriness (decreased probability of spiking for some time after producing a spike).
We assumed that a spike emission is a random point process with the probability (1) where ͕t i ͖ denotes the spiking history earlier than time t, 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 t last of occurrence of the last spike at time: ͑t Խ ͕t i ͖͒ ϭ ͑t, t last ͒. A further assumption is that the firing-rate modulation and refractory effects are multiplicative, thus reflecting the reduction of spike probability due to, for example, inactivation of sodium channels or hyperpolarization caused by opening of potassium channels (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 interval since 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 they are fitted to experimental data by means of a maximum likelihood approach. To capture fine temporal details of the neuronal responses (for example, response onset and ISIs) the intensity and recovery functions were defined with a 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 above conditional intensity (2; Johnson and Swami, 1983;Dayan and Abbott, 2001): ( 3) 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. The likelihood L of obtaining the experimental spike train t i is maximized with respect to the parameters q(t) and w(⌬t) 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, i.e., we require w(⌬t Ͼ 5 ms) ϭ 1 ( Fig. 2A).
To study the effects of refractoriness on the modeled responses, we compared the results to the STPM without refractory period (nonrefractory STPM, w(⌬t) ϭ 1 for all ⌬t Ͼ 0). The model is fully characterized by its intensity function q(t), which can be estimated directly from the experimental PSTH (bin width set to 0.05 ms to allow for sufficient temporal precision).

Generalized linear model (GLM)
One limitation of the STPM is that the history effects are restricted to the last spike only. To evaluate effects evoked beyond the last spike, we considered the GLM (Truccolo et al., 2005;Czanner et al., 2015) with conditional intensity GLM ͑t Խ ͕t i ͖͒ of the form where s(t) is the driving force and h() is the spike history kernel.
Note that the intensity function q(t) of the STPM can be identified with exp ͑s͑t͒͒, and the recovery function w(⌬t) corresponds to exp͑ ͚ i h͑t Ϫ t i ͒͒. In contrast to the STPM, in the GLM the effects of the previous spikes can extend infinitely back in time. In practice, we reduce the number of free parameters of the GLM by restricting the history horizon above which the spikes cannot contribute to the responses anymore; we thus set h͑t Ͼ t max ͒ ϭ 1. The horizon t max ϭ 8 ms was selected to maximize the Akaike Information Criterion (AIC), which balances the goodness of fit with the number of free parameters of the model (Fig.  2C).
The likelihood of the GLM is defined analogously to the STPM: where the sums go over all spikes.
Since the log-likelihood function is a convex function of the parameters, they can be found using standard optimization techniques. In the results presented here we used the conjugate gradient optimization.
We compared the goodness-of-fit of the STPM and the GLM using the time-wrapping method (Brown et al., 2002): The ISIs in the experimental data were rescaled to account for temporal variations in firing probability. If the model perfectly reproduced the data the distribution of the rescaled ISIs would be uniform (Fig. 2E, diagonal).

Model validation
To validate the model, the dataset was divided into two nonoverlapping subsets of equal size: a training and a 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 separately for two statistics X, that is, the PSTH (with bin size 0.2 ms) and the spike pattern distribution. For each of the two statistics, the model error was evaluated as the normalised differences between the simulated X model and validations spike trains X validate (cf. Rauch et al., 2003): where X ϭ X iϭl N is either the PSTH or the spike pattern distribution of model (X model ) or validation (X validate ) set; N is the size of the vector and equals the number of bins (N ϭ 70 for T ϭ 14 ms and 0.2-ms bins) or the number of identified spike patterns (N Յ 16 for binary words of length less or equal to 4).
The model error Err͑X model , X validate ͒ was compared against the error between the training and validation sets Err͑X train , X validate ͒ (reference error). The significance of the difference between the model and reference errors 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 neurons we identified spike triplets defined as three consecutive spikes sepa-rated by intervals shorter than 4 ms. In this analysis, to increase the number of intervals, we broadened the analysis window to 50 ms after the stimulus. Next, we calcu- New Research lated Pearson's correlation between the 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 (r model ) for the same number of trials, which were generated by the STPM model with parameters fitted to the experimental spike trains. The significance of the differences between correlation coefficients found in simulated and experimental ISIs was evaluated by means of a bootstrap test. To this end, 1000 estimates of r model were obtained from independently simulated datasets, and the resulting coefficients were compared to r data. The p value was taken as the smaller of two values multiplied by 2: (1) the fraction of bootstrap trials for which r model was greater than r data ; or (2) the fraction of bootstrap trails for which r model was smaller than r data (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 STPM to a population of uncoupled neurons receiving common inputs.
The population response was calculated from a simulated ensemble of 5000 identical neurons. The parameters of the STPM were fitted to the responses of the analyzed neuron, and these parameters were shared by all 5000 model neurons. In each trial j the intensity function of all neurons was modulated by a multiplicative gain factor G j that was drawn from a uniform distribution on the interval ͓1 Ϫ ␥, 1 ϩ ␥͔, where 0 Յ ␥ Յ 1 is the strength of modulation. The intensity function in trial j was then q j (t) ϭ G j q(t). From the obtained single-trial single-neuron responses the total population response was calculated by summing the binned spike responses of all neurons (bin size 0.2 ms) and subsequent bandpass filtering (400-1200 Hz) corresponding to the analysis of EEG data.
Next, we randomly selected a single neuron from the population and used its spikes for further analysis. We classified the spike patterns of this neuron in single trials based on the occurrence/omission of spikes in a discrete sequence of spiking ЉwindowsЉ. The bandpass filtered population response was then averaged over trials with respect to the type of concomitant spike pattern. This procedure, when applied to the model, reproduced the analysis that was applied to the experimental data and described above (see below, Spike pattern classification). The root mean square (RMS) amplitude of the patternspecific average was compared with the experimentallyobtained hf-EEG related to the same spike pattern (Telenczuk et al., 2011). The similarity of the values across different spike patterns was quantified by means of Pearson's correlation coefficient.

Biophysical model
To understand the mechanisms of burst generation, we developed a simplified single-neuron model. The model consists of a linear neuron with a spiking threshold (leaky integrate-and-fire, LIF), which receives conductancebased inputs through depressing synapses (short-term synaptic depression). The membrane potential in the model follows the standard membrane equation: where C m is the membrane capacitance, g L is the leak conductance, V rest is the resting potential, and I syn are the synaptic currents. When the membrane potential reaches the threshold V thr a spike is generated and the potential is reset to V reset putting the cell into a hyperpolarized state. The total synaptic current in the leaky integrate-and-fire (LIF) neuron is a sum of intracortical and thalamocortical currents: I syn ͑t͒ ϭ I Cortex ͑t͒ ϩ I Th ͑t͒.
The cortical synaptic currents are conductance-based inputs from n inh inhibitory and n exc excitatory neurons. The excitatory g exc and inhibitory g inh synaptic conductance are a sum of contributions mediated by each spike, such The times of the excitatory t exc and inhibitory t inh synaptic inputs are drawn from a homogeneous Poisson process with equal rates for excitatory and inhibitory inputs f exc ϭ f inh . Each spike results in a transient increase of the synaptic conductance with an exponential time course: for t Ն t sp and 0 otherwise. Here t sp is the time of the spike, w is the synaptic weight and is the synaptic time constant. The reversal potentials for excitation and inhibition are E exc and I inh, respectively. With these definitions the total current of cortical synapses is: In addition to the intracortical inputs, the neuron receives excitation from n Th thalamocortical excitatory neurons. The thalamocortical neurons are silent in absence of peripheral stimulation and generate Poisson-distributed spikes 7.68 ms after the onset of the median nerve stimulus (the delay takes account of the propagation delays from periphery to the cortex). The strength of thalamocortical excitatory synapses providing the feedforward inputs to the model decays with the presynaptic activity following the short-term synaptic depression mechanism (Tsodyks and Markram, 1997): where y i and z i are fractions of synpatic resources in active and inactive states, 1 is the decay constant of synaptic conductance, rec is the recovery time from synaptic depression and U describes the fraction of available resources used by each presynpatic spike. The postsynaptic current due to the thalamocortical inputs is then: where the total conductance due to thalamic inputs is given by: g Th,i stands for conductance of a single synapse and y i its efficiency.
Eight parameters of the model were adjusted to reproduce the experimental PSTH: weights of excitatory (w exc ), inhibitory (w inh ) and thalamocortical (w Th ) synaptic inputs, excitatory synaptic time constant ( exc ), firing rates of cortical (f exc ) and thalamocortical (f Th ) presynaptic neurons, number of thalamocortical synapses (n Th ), and use of synaptic resources by thalamocortical synapse release (U). Other parameters were fixed to values found in the literature. The values of other parameters are given in Table 1.

Results
Neurons in area 3b of macaque monkeys (S1) show brief (Ͻ10 ms) bursts of activity in response to stereotypical electrical stimulation of the median nerve (0.2-ms pulse, 1.5 time motor threshold applied transcutaneously to the median nerve; see also Fig. 1A). In a dataset of 46 neurons recorded extracellularly using movable platinumglass electrodes (Eckhorn drive, Thomas Recordings) we found 17 neurons that responded with burst of spikes (defined as trains of two or more spikes with ISIs shorter than 1.8 ms). When averaged over several repetitions of the stimulation the responses gave rise to a poststimulus time histogram (PSTH) with prominent peaks coincident with the within-burst spikes (Fig. 1B). The appearance of such PSTH peaks points to the precision of the burst timing with respect to the onset of the stimulus.
Some of these bursting neurons also elicited spikes in absence of median nerve stimulations (five neurons fired at least 10% of bursts in the window [30, 300] ms after the stimulus). The evoked and spontaneous bursts differed slightly with respect to mode of the within-burst interval distribution [evoked: 1.82 (1.71) ms; spontaneous: 1.32 (0.41) ms, mean (SD) across neurons] and burst length [evoked: 2.76 (1.26) spikes per burst; spontaneous: 2.18 (0.39) spikes per burst], but these differences were not statistically significant (t test, p Ͼ 0.01).
To understand the mechanisms underlying bursting of neurons in the S1, we propose a phenomenological model of the single-neuron response to the median nerve stimulation. The model is based on two experimental observations. (1) On presentation of strong sensory stimuli, layer IV cortical neurons are bombarded with intense and coincident synaptic inputs from thalamocortical neurons (Gil et al., 1999;Swadlow and Gusev, 2001;Hanajima et al., 2004;Bruno and Sakmann, 2006;Cruikshank et al., 2007). (2) After emitting a spike, neurons are refractory, which limits their maximum firing rate (Gray, 1967; Berry  Gil et al. (1997) Value column indicates typical parameter values or ranges found in the literature (where available); ‫ء‬ denotes the parameters which were adjusted to fit the experimental data. and Meister, 1998; Kara et al., 2000). To illustrate the effects of these two phenomena on neuronal responses, we simulated a probabilistic model (the STPM; see Methods) with an exponentially decaying intensity function and an absolute refractory period ref ϭ 1.2 ms (Fig. 1C). The PSTH of the simulated spike responses qualitatively 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 (three peaks visible in Figure 1C: at 6, 7.8, and 9.5 ms after the stimulus) separated by deep troughs corresponding to periods of quiescence during which the neuron is refractory. The presence of such peaks and troughs in the trial-averaged PSTH is possible because the responses of the neuron are reliable across trials. The first peak of the PSTH reflects the initial spike triggered by the sharp transient of the intensity function (Fig. 1C, red line). This initial response is highly reliable, giving rise to the narrowest and tallest PSTH peak (half-amplitude width: 0.75 ms; peak-to-trough amplitude: 1773 spikes/s in Fig.  1C, black line). The refractory state following the first spike leads to a pronounced decrease of firing probability and gives rise to the deep PSTH trough following the initial PSTH peak. The subsequent PSTH peaks become wider and are of smaller amplitude due to the gradual decay of the intensity function (second peak: 1.25 ms, 1002 spikes/s; third peak: 1.25 ms, 290 spikes/s). The PSTH obtained from this simulation is qualitatively similar to cortical burst responses triggered by peripheral nerve stimulation (compare Figure 1C, left, with B, bottom).

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 A B C D

E F
x y z 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).
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. 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 2 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 timedependent 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 (twosided 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 steplike 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).
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-totrial 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).
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-passfiltered 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).
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. 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. 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 Poissondistributed 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 shortterm 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.

Discussion
By means of simplified phenomenological models and a biophysical point-neuron model, we showed that withinburst 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 synchronization 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 S1 results from the combination of intense synaptic bombardment and a refractory period. Such fast bursting triggered and sustained by an intense synaptic input has been termed "forced bursting" (Izhikevich, 2006).
The shape of a fitted recovery function in both models agrees well with the contribution of an after-hyperpolarization (AHP)-mediated potassium current and an after-depolarization (ADP) due to either persistent (Bal and McCormick, 1996;Brumberg et al., 2000) sodium or low-threshold calcium current (Jahnsen and Llinás, 1984): 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 (Fig. 4C). We also demonstrated in a toy model that the over-shoot is not critical for bursting responses, the absolute refractory period combined with intense but transient inputs is sufficient to produce bursts with similar (but not identical) statistics (Fig. 1C).
The STPM could also account for the correlations between ISIs (serial correlations). However, in a few neurons 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, which considers spike-history effects extending to multiple ISIs. 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. The significant serial correlations could be also explained by a model which includes trial-to-trial variations of the input intensity (gain modulation). We found that introducing such variations in the model results in the over-estimation of the serial correlations estimated from the simulated data. In addition, these variations may lead to the over-estimation of the recovery function in form of the overshoot appearing briefly after the absolute refractory period. Although such an overshoot is also present in the recovery function estimated from the data, we believe that it is not an artifact of the estimation method. First, the modulation must be strong (␥ ϭ 0.8) to produce a visible overshoot, whereas we found that modest modulation (␥ ϭ 0.2) is consistent with the serial correlation and EEG correlation estimated in the data. Secondly, we found that GLM is robust with respect to such modulation introduced in the simulated model, but still it uncovers an overshoot in the experimental data. Nevertheless, in the future it would be instructive to extend the spike-train models (STPM and GLM) with the fluctuating gain factor and fit it directly to the data.
We were able to reproduce qualitatively both the average and single-trial features of the burst responses in a more realistic LIF neuron. Although such models are a gross simplification of the real neurons both in terms of spiking mechanism and morphologic features, it has been suggested that the 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 hyperpolarized 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 presynaptic terminal (Markram and Tsodyks, 1996). However, it would be possible to replace it with some other form of adaptation (Brette and Gerstner, 2005). Both mechanisms 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 recordings from thalamocortical projection neurons we cannot 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 intracortical 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 characterized by broad spikes, which suggests that they are pyramidal neurons or spiny stellate neurons (Baker et al., 2003). This is confirmed by intracellular recordings in barrel cortex showing that regular spiking 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 regular spiking neurons. However, a subset of neurons analyzed here (5 of 17 neurons) did also fire bursts that were not locked to the median nerve stimulation showing that at least some of them may belong to the intrinsic bursting class.

Burst synchronization
A striking feature of the S1 bursting is that the signature of the burst also appears in macroscopic signals such as the EEG. The visibility of the burst in the surface recordings was interpreted as a sign of strong synchronization 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 synchronization 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 (Diba and Buzsáki, 2007;Schmidt et al., 2009).
Here, we proposed a model in which the occurrence of spike patterns is regulated by the input intensity, that is the rate of incoming spikes; in contrast precise timing of the input was not necessary. 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 (Fig. 5B). 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 (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;Tsodyks and Markram, 1997;Maass and Zador, 1999;Song et al., 2000;Swadlow and Gusev, 2001), they are optimally placed to represent neuronal variables in a form that is easily processed, stored and transmitted (Leibold et al., 2008;. In this spike-timing-based view, neural systems take advantage of the temporal information encoded into spike patterns to represent slowlychanging 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 .
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 (Steriade et al., 2001;Shu et al., 2003;Fontanini and Katz, 2008). Similar gain control mechanisms were implemented in realistic neural models through, for example, concurrent modulation of excitation and inhibition (Hô and Destexhe, 2000;Chance et al., 2002;Vogels and Abbott, 2009) or short-term synaptic depression (Rothman et al., 2009). More generally, multiplicative noise can account for the variability and covariability of neuronal responses in the thalamus and many cortical areas, including the lateral geniculate nucleus, V1, V2, and MT (Goris et al., 2014).