Synchrony Makes Neurons Fire in Sequence, and Stimulus Properties Determine Who Is Ahead

Surgical procedures.

In five adult cats (three males; two females), anesthesia was induced with a mixture of ketamine (Ketanest; Parke-Davis; 10 mg/kg, i.m.) and xylazine (Rompun; Bayer; 2 mg/kg, i.m.). After tracheotomy, anesthesia was maintained with 70% N₂O and 30% O₂, supplemented with ∼1.0% halothane. After craniotomy, the level of halothane was reduced to ∼0.5%. When anesthesia was stable and sufficiently strong to prevent any vegetative reactions, the animal was paralyzed with pancuronium bromide applied intravenously (pancuronium; Organon; 0.15 mg · kg⁻¹ · h⁻¹). The respiration rate was adjusted manually to keep end-tidal CO₂ at 3–4%. Rectal temperature was kept at 37–38°C. Glucose and electrolytes were supplemented intravenously. All procedures complied with the German law for the protection of animals and were overseen by a veterinarian.

Data acquisition.

Recordings were acquired from area 17 using silicon-based 16-site probes supplied by the Center for Neural Communication Technology at the University of Michigan (illustrated in supplemental Fig. S7, available at www.jneurosci.org as supplemental material), as described previously by Biederlack et al. (2006). The electrode contacts (recording sites) had an impedance of 0.3–0.5 MΩ at 1000 Hz and were organized in a 4 × 4 matrix on four shanks, with a distance of 0.2 mm between the neighboring contacts. Thus, the recording area of one probe spanned ∼0.6 × 0.6 mm. The probes were carefully placed to enter the cortex at a 90° angle (i.e., perpendicularly to the surface of the cortex), and the uppermost row of electrode contacts was lowered 100–200 μm below the cortex surface. Thus, the same probe recorded neurons from different layers and different orientation columns. We always inserted two probes simultaneously into area 17 of the same hemisphere, thus recording from 32 electrode contacts in parallel. In one probe, good signals were available only on nine recording channels. Thus, in total, we recorded from 153 sites.

Signals were amplified 1000×, filtered at 500 Hz to 3.5 kHz, and sampled with a frequency of 32 kHz. Action potentials (spikes) were detected with a two-sided threshold discriminator adjusted manually to yield a signal-to-noise ratio of ∼2:1. For each detected action potential, the time of the event was recorded together with the spike waveform over a duration of 1.2 ms. Single-unit activity was then extracted off-line by spike sorting. The distance between neighboring electrode contacts was sufficient to prevent them from recording activity of the same neuron—in no case did cross-correlation histograms (CCHs) show the high, narrow peaks in the very center of the CCH characteristic of autocorrelations. Thus, spike sorting was applied to the signals from each recording site separately. Spike waveforms were sorted using principal component analysis (PCA), according to the following two criteria: (1) waveform separation: the shape of the spike waveforms should be clearly defined and different from all other waveforms in the multiunit (MU), as indicated by the first three PCA components; (2) refractory period: not more than 0.5% of spikes should occur with an interspike interval of <1 ms. The 126 isolated single units (SUs) that satisfied these criteria (7–20 neurons per recording probe) had an average firing rate of 18 ± 15 spikes/s (mean ± 2 SD) in the most optimal stimulation condition.

Visual stimulation.

Visual stimuli were presented on a 21 inch Hitachi CM813ETPlus monitor with a refresh rate of 100 Hz, positioned at a distance of 57 cm from the eyes. The pupils were dilated with atropine, the nictitating membrane retracted with neosynephrine and the cornea protected with contact lenses containing artificial pupils of 2 mm diameter. After refraction, the eyes were focused on the monitor with correcting lenses. The optical axes of the two eyes were aligned by mapping the receptive fields (RFs) for each eye separately and placing a prism in front of one eye to achieve binocular fusion of the RFs. All RFs were determined by manual mapping with high-contrast white bars on a black background. Because of the spatial proximity between electrode contacts, the RFs recorded from one probe always overlapped, producing clusters that spanned up to ∼10° of visual angle. The stimuli were always positioned such that their centers matched the center of the RF cluster.

Stimuli were always presented binocularly, using the presentation software ActiveSTIM (www.ActiveSTIM.com). In the main set of recordings, for which the largest part of the reported analyses was performed, we presented high-contrast drifting sinusoidal gratings of circular shape, and varied the orientation and drifting direction in 12 30° steps, thereby covering the entire circle of 360°. In three cats, the gratings had a spatial frequency of 2.4°/cycle and covered 12° of visual angle, in two cases moving with a speed of 2.2°/s and in one with 6.6°/s. In the remaining two cats, gratings had a spatial frequency of 1 and 2.5°/cycle, respectively, covering 7 and 18° of visual angle and moving with a speed of 1.5 and 2°/s.

To exclude stimulus-locked spike timing as a cause for the time delays in CCHs, in three animals we additionally recorded responses to high-contrast moving bars. The bars were oriented to match the orientation preference of a majority of the cells and were moved across the receptive fields in the two directions perpendicular to the orientation of the bar. Bars spanned 0.5 × 6° of visual angle, appeared on the screen at an eccentricity of 2° from the center of the RF cluster, and moved across the center of the cluster with a speed of 1°/s, disappearing after 4 s.

All stimuli were presented 40–120 times (trials) each, in randomized order. The recordings were made in two to six blocks, each block consisting of 20–60 trials per stimulus condition. Each block lasted 20–100 min, and often several-hour-long breaks were made between different blocks. This enabled us to investigate the stability of the relative firing times over longer periods of time. The responses obtained from the two simultaneously inserted Michigan probes were analyzed separately (i.e., timing relationships between units recorded from different probes were not investigated here). Thus, for each cat we obtained two datasets (denoted A and B). Thus, the total of 8 experiments (5× drifting gratings + 3× bars) yielded a total of 16 data sets, 10 in which we investigated the stimulus-induced changes in the firing sequences in response to grating stimuli, and 6 with responses to bar stimuli used to analyze the relationship between the preferred spike delays and the external temporal dynamics of the stimulus (i.e., to prove the internal origin of the delays).

Extracting pairwise preferred delays.

The preferred spike delay between pairs of neurons was determined by computing their CCH with 1 ms resolution across the sustained responses to each stimulus (i.e., excluding the responses to stimulus onset and offset to eliminate rate covariation from the CCHs). To ensure that CCHs did not include rate covariation (Brody, 1999), the responses showing the strongest rate variation (i.e., those immediately after the onset and offset of the stimuli) were not included in the CCH analysis. Hence, CCHs were computed only for the sustained responses over a duration of 3.2 or 1.8 s for grating and bar stimuli, respectively. Stimulus-locked rate covariation can be detected by computing CCHs between responses with shuffled trial orders (i.e., shift predictors) (Perkel et al., 1967). We computed shift predictors for a random sample of four datasets and most of these 423 shifted CCHs were flat, indicating that rate covariation had not contributed to the peaks in the CCHs. More importantly, regardless of the shapes of the shift predictors, the estimates of phase shifts were identical with and without subtraction of shift predictors (r = 0.99; n = 423).

We then fitted a Gaussian function to the central peak of the CCH (±15 ms) and located its maximum on the x-axis (see Fig. 1B). When applied on reasonably smooth CCHs, particularly with 1 ms binning, this approach estimates preferred spike delays with submillisecond precision (Havenith et al., 2009). Compared with other typical fitting methods for CCHs, like damped Gaussian (i.e., Gabor) functions (König, 1994; König et al., 1995a; Roelfsema et al., 1997; Nikolić, 2007) or cosine functions (Schneider and Nikolić, 2006; Schneider et al., 2006), the Gaussian function applied here produced the same results as the other two (r > 0.95; 91 CCHs tested) but was least sensitive to initial parameter settings and least likely to enter local minima during the fitting process.

The submillisecond precision of the fitting procedure could be achieved only if the CCH peak consisted of at least three consecutive bins (1 ms width each) with at least N = ∼8 entries (coincident events). With fewer entries, precision decreased rapidly because of fitting errors (Havenith et al., 2009). This requirement had several implications for the selection of the data included in the analyses as follows.

• (1) Not all neurons in all stimulation conditions generated enough spikes to enable a reliable assessment of the preferred time delays from the pairwise CCHs. This was usually a problem with stimuli of nonoptimal orientation, which did not drive the neurons well. If a neuron had overall low firing rates and participated in CCHs with insufficient entries, it was completely excluded from the analysis. However, we did retain neurons that only had “poor” CCHs in some stimulation conditions. The reason is that our analysis required that the networks of pairwise preferred delays were complete (i.e., did not contain missing values). The criterion for the inclusion of a unit was that the Gaussian functions fitted the CCHs with r² ≥ 0.5 in >50% of the cases. Since the overall preferred firing time of a neuron was computed from multiple pairwise delays, this criterion provided sufficiently accurate estimates. This procedure eliminated 48 of the 126 isolated SUs (38%), leaving 78 neurons for additional analysis (5–12 per recording probe; median, 7 neurons). This criterion could not be applied for the polar plots of firing times shown in Figure 5B, and for the time-resolved changes in firing times shown in Figure 5, E and F. These analyses needed to include also responses with very low spike counts (e.g., for nonpreferred grating orientations). Thus, all neurons with nonempty CCHs were included in these analyses, regardless of the goodness of fit of Gaussian functions.

• (2) As mentioned, some grating orientations evoked particularly low firing rates in the recorded neurons, producing very sparse, or sometimes even completely empty, CCHs. Therefore, in some analyses, we had to exclude the responses to certain grating orientations. A stimulation condition was removed from the analysis if it was not possible to compute a firing sequence for five or more neurons after applying the criterion explained in (1). In one data set, neurons responded vigorously enough to compute firing sequences for all 12 stimulation conditions. The smallest number of stimulation conditions that entered the analysis was two, and the median was 3.5. As mentioned above, we did not use these criteria for the computation of the orientation tuning of the relative firing times.

• (3) Finally, it was not feasible to estimate the intertrial variability of firing sequences directly by determining pairwise time delays for single trials and to make at the same time a reasonable comparison to firing rates. This is because, as a result of low firing rates, insufficient entries were available to estimate CCHs from single trials, which lead to a disproportionate increase in the error of time delay estimates compared with estimates of firing rates [for an illustration of the trial-to-trial variability in estimated time delays, see Nikolić (2007) and supplemental Fig. S9, available at www.jneurosci.org as supplemental material]. This lack of distributivity of errors in the estimates of time delays is not an indication of inherent imprecision of time delays but occurs because of the fact that the measure of delay is extracted indirectly, by fitting a Gaussian function. The fitting errors of each parameter do not reduce with equal speed. For example, the baseline of the Gaussian can be estimated with accuracy that is proportional to 1/n, where n would denote the amount of available data (e.g., the number of trials). However, other parameters (e.g., the mean and the width) can be estimated with a similar rate only when the baseline is estimated accurately, and the estimate of the peak delay of the Gaussian is in turn reliant on all three of these other parameters. Consequently, the error with which the mean delay of the Gaussian is measured, is for small numbers of data points, much larger than would be expected from the widths of the center peaks. In other words, with a small number of trials (that is, CCHs containing only a small number of coincident events), the best fit is likely to position the Gaussian at values far outside the region within which the center peaks of CCHs vary (e.g., an estimated delay of 30 ms can be easily encountered). Consequently, an average of k Gaussian fits, each of which is applied to 1/k of the total available amount of data, will be much more inaccurate than a single Gaussian fit calculated for a CCH including all data. This is referred to as the lacking distributivity of errors. Thus, to obtain accurate estimates in single trials, each trial would be required to contain an unrealistically high number of entries (e.g., rich MU instead of SU). Therefore, time delays had to be estimated always from CCHs calculated across multiple trials. Because of these limitations, to estimate the variability of firing times across repeated measurements we randomly split the data sets into two halves 100 times and calculated CCHs from each half of the data. This restriction does not apply to rate responses, as they can be calculated for n individual trials and subsequently averaged to obtain an estimate with the error decreasing with the usual factor of ∼1/n.

Delays in RF activation.

To estimate the delays in activation caused by the RF properties of different units, we used two different techniques. For the data shown in Figure 1D and supplemental Figure S1, B and C (available at www.jneurosci.org as supplemental material), we applied the double sliding window technique for the measurement of onset latencies, based on peristimulus time histograms (PSTHs) (Berényi et al., 2007). We computed PSTHs with a resolution of 5 ms, and the width of the sliding window was 400 ms. The method rarely entered local minima, and the estimated latencies of the response onset corresponded well to those obtained by visual inspection of the PSTHs. For the data shown in supplemental Figure S1, B and C, we computed PSTHs with a bin size of 10 ms, and then generated cross-correlation functions of those PSTHs for delays of ±1000 ms. The delays in PSTH peaks were then estimated by the same procedures based on fitting a Gaussian as used for the regular CCHs.

Firing sequences.

Using the method described by Schneider et al. (2006) and Nikolić (2007), the matrices of cross-correlation delays (also referred to as networks) (see Fig. 2B) between simultaneously recorded neurons were transformed into linear sequences of preferred firing times. This is justified if a data set satisfies the principle of additivity: For any three units A–C, the time delay ϕ between units A and C equals the sum of the time delays between units A and B and units B and C, as follows: If this relationship holds, one can assign each unit k a temporal position indicating the preferred time x_k at which it fires action potentials relative to all other units (see Fig. 2C). x_k is determined as the normalized sum of all time delays ϕ_jk of unit k with respect to all other units j as follows: where n indicates the total number of units. Thus, if a unit tends to fire earlier than others, it is assigned a positive value and vice versa, the sum of all values being always zero. In this way, n(n − 1)/2 pairwise time delays are condensed into n temporal positions on a single time axis (see Fig. 2F). To estimate how accurately this linear representation matched the original pairwise delays, we computed a measure of the error as the normalized difference between the distance of the positions of two neurons on the time axis on the one hand, and the preferred delay measured in their CCH on the other. This error is referred to as additivity error σ_Add because it is directly related to the deviation of the original pairwise delays from additivity (see Eq. 1). The overall error σ_Add takes into account all possible pairs of neurons, whereas the individual error σ_Add(k) of a unit k only takes into account its pairwise connections with all other neurons. For more details, see Schneider et al. (2006).

Error measures of firing sequences.

To estimate the error variance of the preferred firing times across repeated measurements, we used several different error measures. First, we applied a bootstrapping procedure, whereby we computed firing sequences from a randomly chosen subsample of trials. Each subsample contained one-half of the total number of available trials, which were sampled with replacement. For each data set, 100 such samples were made and the SD of the obtained relative firing times was used as an estimate of the measurement error, referred to as “bootstrap error,” and is shown in Figures 3E and 4C, and also used in the computation of mutual information detailed below.

To estimate the stability of the firing sequences over time, we split the data in two halves as they were recorded in time (i.e., early vs late recording blocks) (unlike the random choice of trials used in the bootstrapping procedure). For the number of recordings used in each analysis and other details, see supplemental Table S1 (available at www.jneurosci.org as supplemental material).

Finally, we also computed a measure of the error variance resulting from the imperfect additivity of the measured pairwise delays. Pairwise spiking delays are not guaranteed to be additive. In principle, spike trains can also produce nonadditive time delays (Schneider et al., 2006). Therefore, we needed to ensure that the time delays were sufficiently additive to warrant the application of our analyses. The degree to which the preferred firing sequence correctly represents the measured pairwise time delays is expressed in the average positioning error per unit, σ_Add [defined as σ_Add≡σ̂_x in the study by Schneider et al. (2006)]. To obtain σ_Add, first the differences between the distance δ_ij separating two units i and j on the time axis and their measured time delay ϕ_ij are summed up as follows: Next, Q_Add is normalized by the number of time delays that enter its computation [adapted from the study by Schneider et al. (2006), their Eq. 4] as follows: σ_Add is also referred to as the additivity error because the measured time delay between two units and their distance on the time axis correspond only if a data set is additive (see Fig. 2C–E).

In addition to this average error across all units, we computed an individual positioning error σ_Add(k) for each unit k, which was used to determine the widths of the Gaussians above each unit in the representations of firing sequences shown in Figures 2F, 3A,G, 4A, and 5A. Similar to the average error for all units, σ_Add(k) was also computed as a normalized sum of errors, but only for the time delays ϕ_jk of the unit k with all other units j (j ≠ k). It can be shown that the mean of σ²_Add(k) across all units k equals the overall additivity error σ²_Add.

Statistical testing of firing sequences.

To statistically test the stimulus-dependent changes of firing sequences, we applied three tests based on different types of error variance. First, we tested whether the difference between two firing sequences exceeded the variance resulting from the imperfect additivity of pairwise delays. For this purpose, we used the ANOVA presented in the study by Schneider et al. (2006), termed here ANOVA-Add, which is based on the positioning error σ_Add (defined in the previous section), The resulting F value, computed as the ratio between the stimulus-dependent differences between two mean firing sequences, and the average additivity error σ_Add, is distributed with (n − 1) and (n − 1)(n − 2) degrees of freedom and assumes large values if the positioning error is small and the difference between the firing sequences is large. The ANOVA-Add is designed only for comparisons between two firing sequences. Thus, when comparing changes in firing sequences for more than two stimulus conditions, we conducted multiple pairwise ANOVA-Adds and controlled for the type I error by correcting p values using the Dunn–Sidak correction for multiple comparisons (Dunn, 1961). A stimulus-dependent change of firing sequences was considered significant if across all comparisons the smallest corrected p value did not exceed p = 0.05.

Second, we corroborated this parametric analysis with a nonparametric statistical test also designed to investigate pairwise delays (Nikolić, 2007). This test considers only the direction of time delays, ignoring their magnitudes, and tests whether these directions indicate a consistent rank order of firing across the units (i.e., if unit A fires before B, and B fires before C, then A fires before C), a property known as transitivity. Similar to, for example, Wilcoxon's signed-ranks test, changes between two firing sequences are evaluated by computing first the differences between all pairwise time delays across two stimulation conditions, and then testing the transitivity of this difference network (matrix). To this end, all subnetworks of size three (triples) are considered, the total number of nontransitive triples, C_NTrans, is counted, and its significance for a given network size is taken from a look-up table provided in the study by Nikolić (2007). The test makes no assumptions about normal distribution of errors. To correct for multiple comparisons in the transitivity test, we used an α level of 0.01.

Finally, to compare the stimulus-dependent changes of firing sequences against their error variance across multiple measurements, we used a two-way (stimulus by unit) ANOVA for repeated measurements (ANOVA-RM). We obtained repeated measurements by splitting each dataset into two parts, analyzing either odd or even trials (according to the order of presentation). The data sets were only split in two parts because the estimates of preferred firing delays required fitting CCH with Gaussian functions, and this procedure is particularly susceptible to producing large errors when a small number of data points is used (for details, see above, Extracting pairwise preferred delays). On these split-half measures, we applied the ANOVA-RM. Note that a significant change of firing sequences cannot be detected by testing the main effect of the factor “stimulus.” The reason is that the mean of a firing sequence is, by definition, always zero and, therefore, does not change across stimulus conditions. Instead, significant changes in firing sequences have to be detected through the interaction between the factors “stimulus” and “unit,” because the positions of individual units on the time axis change as a function of the stimulus even though the overall mean of all positions is fixed to zero. The same split-half repeated-measurement design was also used for the analysis of rate responses. A comprehensive table of results for all statistical tests on all data sets is provided in supplemental Table S1 (available at www.jneurosci.org as supplemental material) (also see below, Unitary event analysis, for a statistical test for detection of temporal sequence directly in the spike trains).

Sequence-triggered average of the local field potential.

The sequence-triggered average (SQTA) of the local field potential (LFP) displayed in Figure 3B is computed by first averaging the LFP signals recorded across all electrodes with which we recorded also spiking activity. Next, the classical spike-triggered average (STA) of the LFP is computed for each of the neurons in the sequence. Finally, to obtain the SQTA, the relative positions of the individual STAs are offset according to the relative firing times of the corresponding neurons, and then averaged. For example, if a sequence includes two neurons with the preferred firing times of +2 and −2 ms, the two corresponding STAs would be averaged with an offset of 4 ms.

Oscillation frequency of CCHs.

To determine the dominant frequency of oscillation for each CCH, we first generated an extended CCH over a time window of ±128 ms, and then computed its fast Fourier transform to obtain the power spectrum. The power spectra showed a well defined single peak in the range of 20–60 Hz in almost all cases. The frequency of each CCH was defined as the frequency with the highest power within this frequency range.

Orientation and direction tuning of firing times.

In each of the polar plots in Figure 5B, a constant time, t_c, needed to be subtracted from the firing times of the unit relative to all the other units for all 12 stimulus directions. The reason was that, unlike rate responses, relative firing times do not have an absolute zero. Hence, the choice of t_c was arbitrary and, in turn, affected strongly the degree to which firing times appeared direction selective (i.e., the “sharpness” of the tuning curve). Thus, in Figure 5, B and C, the sharpness of orientation and direction selectivity should not be compared across firing rates and firing times. We chose values of t_c that approximately matched the sharpness of the corresponding rate tuning curves to facilitate visual comparisons of their orientation and direction preferences. Importantly, the subtraction of t_c neither affects the estimates of the preferred stimulus direction nor the Pearson's correlation coefficient, r, calculated between firing rates and firing times across the 12 stimulation conditions.

Time-resolved analysis.

The time-resolved analysis of firing times and firing rates was applied to three data sets of 5–7 SUs (n = 17 units in total) that were taken from the three oscillatory recordings. Firing sequences and mean firing rates were extracted for analysis windows of 450–550 ms duration that were moved in 10 ms steps. All the units in a group had overlapping RFs. To compare the modulation amplitudes of firing rates and firing times, we fitted both types of responses with sine functions of the same temporal frequency as the presented grating stimuli. For firing rates, the modulation amplitude was then determined by comparing the amplitude of the sinusoidal modulation of the rate response by the grating (F1) to the mean change in response amplitude from spontaneous to stimulus-evoked activity (F0), known as F1/F0 ratio (Movshon et al., 1978a,b; Skottun et al., 1991). Rate responses with an F1/F0 ratio >1 were classified as simple, whereas rate responses with an F1/F0 ratio <1 were classified as complex. For firing times, the modulation amplitude of the fitted sine was not normalized, but simply expressed in milliseconds.

Mutual information.

The mutual information (MI) was calculated between rate or time responses (r) of individual neurons on one side, and the visual stimuli (s), on the other side by using the following equation: where p(r) and p(s) are probability distributions of responses and stimuli, respectively, and p(r,s) is the joint probability of responses and stimuli. In the experiment, the p(s) was kept constant for all possible stimuli [i.e., p(s) = 1/n, where n is the number of different stimuli]. p(r) and p(r,s) were calculated based on conditional probability p(r|s) by using the following equations: and The conditional probability p(r | s) for both firing rates and firing times was approximated by assuming that the responses follow Gaussian distribution for which the mean and SD were estimated by a bootstrapping procedure described above (see Error measures of firing sequences).

The maximum possible entropy of stimuli, MI_max, is the upper bound for mutual information with other variables and varied between 1 bit, for two different stimuli, and 3.6 bits, for 12 stimuli. Consequently, the capacity for carrying stimulus-related information, pMI, was expressed as proportion of the maximum possible information that was transferred: The MI and pMI values for rate responses and for relative firing times of a neuron were calculated from identical spike trains.

Oscillation score.

The oscillation strength in a data set was estimated by computing its oscillation score (OS) (Muresan et al., 2008). To this end, an autocorrelation histogram (ACH) was computed with 1 ms resolution for each unit and for each investigated stimulus condition. The OS is defined as the ratio between the power at the strongest frequency within a band of interest, which covers in our case the upper beta/lower gamma band (25–50 Hz), and the average power of all other frequencies in that band. The OS computes a corrected oscillation power by removing the artifacts induced by the presence of the center peak in the ACH and by possible bursting behavior of the cells. An OS of 5 for a preferred frequency of 40 Hz indicates that the oscillatory power at 40 Hz is five times larger than the average power in all the other investigated frequencies.

Within a data set, the units tended to have similar OS, indicating a collective fluctuation in the level of oscillatory activity. Thus, the overall oscillation strength for a data set was estimated simply as the average OS for all individual ACHs. In accordance with the recommendations given by Muresan et al. (2008), data sets with an average OS ≥ 5 were classified as oscillatory, and data sets with an OS < 5 as nonoscillatory. For the responses to gratings drifting in different directions, seven data sets were classified as nonoscillatory (OS, 1.5–4.8; mean, 3.2; median, 3.5) and three as oscillatory (OS, 6.0, 6.3, and 7.2).

Unitary event analysis.

Unitary events were determined by dividing the spike responses to gratings (300 ms after stimulus onset until stimulus offset) into short consecutive time windows (2 and 5 ms; results shown only for 5 ms windows). A unitary event was defined as a group of spikes fired within a time window, and the size of the unitary event was defined as the number of spikes in that window. The frequency of occurrence (events per second) was computed for each size of unitary event.

To investigate the stimulus specificity of unitary event distributions, spike trains were shifted in two ways: In one, the shift of each spike train compensated for the delay in the relative firing time of that neuron in the preferred firing sequence (“correct shift”). These shifts were expected to reduce maximally the average delay between action potentials. Hence, this procedure should maximize the number of large unitary events (and concomitantly reduce the number of small unitary events). In the second type of shift, spike trains were shifted according to the preferred firing sequences obtained in response to the other stimuli (“incorrect shift”). If the temporal structure of the spike trains is stimulus dependent, this type of shift should either not reduce the delays between action potentials (e.g., if stimuli are mostly dissimilar) or reduce them less efficiently (e.g., if similar stimuli are used for correct and incorrect shifts).

To assess the significance of the differences in the sizes of unitary events f, we computed the average frequency of occurrence (number of events per second) in spike trains shifted either correctly, F_c, or incorrectly, F_i, and then computed the normalized frequency based on the average for correct and incorrect shifts: F_nc = F_c/(½ F_c + ½ F_i) and F_ni = F_i/(½ F_c + ½ F_i).

Given that the frequencies with which unitary events of different sizes occurred were mutually dependent—if the number of large unitary events increased, the number of small ones necessarily decreased—we had to test for the significance of a change by examining the interaction between the correctness of the of the shift and the size of the unitary event in a two-way ANOVA. The p value of the interaction term indicated whether shifting the spike trains in the different ways affected the shapes of the distributions of unitary event sizes. The F_n ratios obtained for each stimulus condition were treated as repeated measurements and the null-hypothesis (H₀) was that the distributions of unitary event sizes do not increase as a result of a correct shift of a spike train compared with incorrect shift. The results of this test are listed in supplemental Table S1 (available at www.jneurosci.org as supplemental material).

To investigate whether oscillations had an effect on the sizes of unitary events, this approach was modified such that the two-way ANOVA had a factor strength of oscillation instead of correctness of the shift and was made using the F_nc values only (see Fig. 8I).

The significance of changes for each event size separately (see Fig. 8E–H) was tested nonparametrically: We computed ratios between correctly shifted and original spike trains, F_c/F_o, incorrectly shifted and original spike trains, F_i/F_o, and correctly and incorrectly shifted spike trains, F_c/F_i. A binomial test was made for the number of firing sequences with the ratio >1.0, given the total number of firing sequences with unitary events of that size. The binomial probability was computed for obtaining the observed count by chance (i.e., assuming that each firing sequence has equal probability of 0.5 of producing a ratio >1.0 and <1.0).

As a control, we applied all of the above-described analyses to trial-shuffled spike trains. Trial shuffling eliminates unitary events that are generated because of the internal organization but retains those that emerge from the responses that are time-locked to the external stimulus and those that occur by chance. In this case, the average frequencies of unitary events were estimated from 200 repetitions of the shuffling procedure.